################################################### ### chunk number 1: Loading library ################################################### #line 38 "vignettes/Starr/inst/doc/Starr.Rnw" library(Starr) ################################################### ### chunk number 2: Reading bpmap file ################################################### #line 50 "vignettes/Starr/inst/doc/Starr.Rnw" dataPath <- system.file("extdata", package="Starr") bpmapChr1 <- readBpmap(file.path(dataPath, "Scerevisiae_tlg_chr1.bpmap")) ################################################### ### chunk number 3: Data read-in ################################################### #line 59 "vignettes/Starr/inst/doc/Starr.Rnw" cels <- c(file.path(dataPath,"Rpb3_IP_chr1.cel"), file.path(dataPath,"wt_IP_chr1.cel"), file.path(dataPath,"Rpb3_IP2_chr1.cel")) names <- c("rpb3_1", "wt_1","rpb3_2") type <- c("IP", "CONTROL", "IP") rpb3Chr1 <- readCelFile(bpmapChr1, cels, names, type, featureData=T, log.it=T) ################################################### ### chunk number 4: ExpressionSet ################################################### #line 69 "vignettes/Starr/inst/doc/Starr.Rnw" rpb3Chr1 ################################################### ### chunk number 5: assayData ################################################### #line 78 "vignettes/Starr/inst/doc/Starr.Rnw" head(exprs(rpb3Chr1)) ################################################### ### chunk number 6: phenoData ################################################### #line 85 "vignettes/Starr/inst/doc/Starr.Rnw" pData(rpb3Chr1) ################################################### ### chunk number 7: featureData ################################################### #line 94 "vignettes/Starr/inst/doc/Starr.Rnw" featureData(rpb3Chr1) head(featureData(rpb3Chr1)$chr) head(featureData(rpb3Chr1)$seq) head(featureData(rpb3Chr1)$pos) ################################################### ### chunk number 8: Reconstruction of the array image eval=FALSE ################################################### ## #line 108 "vignettes/Starr/inst/doc/Starr.Rnw" ## plotImage(file.path(dataPath,"Rpb3_IP_chr1.cel")) ################################################### ### chunk number 9: Reconstruction of the array image ################################################### #line 112 "vignettes/Starr/inst/doc/Starr.Rnw" jpeg(file="image.jpeg", quality=100) plotImage(file.path(dataPath,"Rpb3_IP_chr1.cel")) dev.off() ################################################### ### chunk number 10: boxplots and density plots eval=FALSE ################################################### ## #line 129 "vignettes/Starr/inst/doc/Starr.Rnw" ## par(mfcol=c(1,2)) ## plotDensity(rpb3Chr1, oneDevice=T, main="") ## plotBoxes(rpb3Chr1) ################################################### ### chunk number 11: boxplots and density plots ################################################### #line 135 "vignettes/Starr/inst/doc/Starr.Rnw" png("boxdens.png", height=400, width=720) par(mfcol=c(1,2)) plotDensity(rpb3Chr1, oneDevice=T, main="") plotBoxes(rpb3Chr1) dev.off() ################################################### ### chunk number 12: Scatterplot matrix eval=FALSE ################################################### ## #line 153 "vignettes/Starr/inst/doc/Starr.Rnw" ## plotScatter(rpb3Chr1, density=T, cex=0.5) ################################################### ### chunk number 13: Scatterplot matrix ################################################### #line 157 "vignettes/Starr/inst/doc/Starr.Rnw" png("densscatter.png", height=400, width=360) plotScatter(rpb3Chr1, density=T, cex=0.5) dev.off() ################################################### ### chunk number 14: MA plot of raw data eval=FALSE ################################################### ## #line 175 "vignettes/Starr/inst/doc/Starr.Rnw" ## ips <- rpb3Chr1$type == "IP" ## controls <- rpb3Chr1$type == "CONTROL" ## plotMA(rpb3Chr1, ip=ips, control=controls) ################################################### ### chunk number 15: MA plot of raw data ################################################### #line 181 "vignettes/Starr/inst/doc/Starr.Rnw" png("maRaw.png", height=400, width=720) ips <- rpb3Chr1$type == "IP" controls <- rpb3Chr1$type == "CONTROL" plotMA(rpb3Chr1, ip=ips, control=controls) dev.off() ################################################### ### chunk number 16: Sequence-specific hybridization bias eval=FALSE ################################################### ## #line 201 "vignettes/Starr/inst/doc/Starr.Rnw" ## par(mfcol=c(1,2)) ## plotGCbias(exprs(rpb3Chr1)[,1], featureData(rpb3Chr1)$seq, main="") ## plotPosBias(exprs(rpb3Chr1)[,1], featureData(rpb3Chr1)$seq) ################################################### ### chunk number 17: Sequence-specific hybridization bias ################################################### #line 207 "vignettes/Starr/inst/doc/Starr.Rnw" png("posGC1.png", height=400, width=720) par(mfcol=c(1,2)) plotGCbias(exprs(rpb3Chr1)[,1], featureData(rpb3Chr1)$seq, main="") plotPosBias(exprs(rpb3Chr1)[,1], featureData(rpb3Chr1)$seq) dev.off() ################################################### ### chunk number 18: ################################################### #line 227 "vignettes/Starr/inst/doc/Starr.Rnw" rpb3_loess <- normalize.Probes(rpb3Chr1, method="loess") ################################################### ### chunk number 19: MA-plot of the normalized data eval=FALSE ################################################### ## #line 236 "vignettes/Starr/inst/doc/Starr.Rnw" ## plotMA(rpb3_loess, ip=ips, control=controls) ################################################### ### chunk number 20: MA-plot of the normalized data ################################################### #line 240 "vignettes/Starr/inst/doc/Starr.Rnw" png("maNorm.png", height=400, width=720) plotMA(rpb3_loess, ip=ips, control=controls) dev.off() ################################################### ### chunk number 21: Calculating ratio ################################################### #line 255 "vignettes/Starr/inst/doc/Starr.Rnw" description <- c("Rpb3vsWT") rpb3_loess_ratio <- getRatio(rpb3_loess, ips, controls, description, fkt=median, featureData=F) ################################################### ### chunk number 22: Sequence-specific hybridization bias (normalized data) eval=FALSE ################################################### ## #line 264 "vignettes/Starr/inst/doc/Starr.Rnw" ## par(mfcol=c(1,2)) ## plotGCbias(exprs(rpb3_loess_ratio)[,1], featureData(rpb3_loess)$seq, main="") ## plotPosBias(exprs(rpb3_loess_ratio)[,1], featureData(rpb3_loess)$seq, ylim=c(-0.5,0.5)) ################################################### ### chunk number 23: Sequence-specific hybridization bias (normalized data) ################################################### #line 270 "vignettes/Starr/inst/doc/Starr.Rnw" png("posGC2.png", height=400, width=720) par(mfcol=c(1,2)) plotGCbias(exprs(rpb3_loess_ratio)[,1], featureData(rpb3_loess)$seq, main="") plotPosBias(exprs(rpb3_loess_ratio)[,1], featureData(rpb3_loess)$seq, ylim=c(-1,1)) dev.off() ################################################### ### chunk number 24: ################################################### #line 299 "vignettes/Starr/inst/doc/Starr.Rnw" probeAnnoChr1 <- bpmapToProbeAnno(bpmapChr1) ################################################### ### chunk number 25: Remapping probes ################################################### #line 313 "vignettes/Starr/inst/doc/Starr.Rnw" newbpmap <- remap(bpmapChr1, path=dataPath, reverse_complementary=TRUE, return_bpmap=TRUE) ################################################### ### chunk number 26: Summary of bpmap ################################################### #line 319 "vignettes/Starr/inst/doc/Starr.Rnw" str(newbpmap) ################################################### ### chunk number 27: Summary of bpmap eval=FALSE ################################################### ## #line 327 "vignettes/Starr/inst/doc/Starr.Rnw" ## writeTpmap("newbpmap.tpmap", newbpmap) ## tpmap2bpmap("newbpmap.tpmap", "newbpmap.bpmap") ## ## pA <- bpmapToProbeAnno(newbpmap) ################################################### ### chunk number 28: ################################################### #line 353 "vignettes/Starr/inst/doc/Starr.Rnw" transcriptAnno <- read.gffAnno(file.path(dataPath, "transcriptAnno.gff"), feature="transcript") filteredIDs <- filterGenes(transcriptAnno, distance_us = 0, distance_ds = 0, minLength = 1000) ################################################### ### chunk number 29: means ################################################### #line 363 "vignettes/Starr/inst/doc/Starr.Rnw" pos <- c("start", "start", "start", "region", "region","region","region", "stop","stop","stop") upstream <- c(500, 0, 250, 0, 0, 500, 500, 500, 0, 250) downstream <- c(0, 500, 250, 0, 500, 0, 500, 0, 500, 250) info <- data.frame(pos=pos, upstream=upstream, downstream=downstream, stringsAsFactors=F) ################################################### ### chunk number 30: means ################################################### #line 376 "vignettes/Starr/inst/doc/Starr.Rnw" means_rpb3 <- getMeans(rpb3_loess_ratio, probeAnnoChr1, transcriptAnno[which(transcriptAnno$name %in% filteredIDs),], info) ################################################### ### chunk number 31: correlationPlot eval=FALSE ################################################### ## #line 386 "vignettes/Starr/inst/doc/Starr.Rnw" ## info$cor <- sapply(means_rpb3, mean, na.rm=T) ## level <- c(1, 1, 2, 3, 4, 5, 6, 1, 1, 2) ## info$level <- level ## correlationPlot(info, labels=c("TSS", "TTS")) ################################################### ### chunk number 32: correlationPlot ################################################### #line 393 "vignettes/Starr/inst/doc/Starr.Rnw" png("corPlot.png", height=400, width=360) info$cor <- sapply(means_rpb3, mean, na.rm=T) level <- c(1, 1, 2, 3, 4, 5, 6, 1, 1, 2) info$level <- level correlationPlot(info, labels=c("TSS", "TTS")) dev.off() ################################################### ### chunk number 33: profileplotExampleData eval=FALSE ################################################### ## #line 420 "vignettes/Starr/inst/doc/Starr.Rnw" ## sampls = 100 ## probes = 63 ## at = (-31:31)*14 ## clus = matrix(rnorm(probes*sampls,sd=1),ncol=probes) ## clus= rbind( t(t(clus)+sin(1:probes/10))+1:nrow(clus)/sampls , t(t(clus)+sin(pi/2+1:probes/10))+1:nrow(clus)/sampls ) ################################################### ### chunk number 34: profileplotExampleData eval=FALSE ################################################### ## #line 430 "vignettes/Starr/inst/doc/Starr.Rnw" ## labs = paste("cluster",kmeans(clus,2)$cluster) ################################################### ### chunk number 35: profileplotExampleData eval=FALSE ################################################### ## #line 437 "vignettes/Starr/inst/doc/Starr.Rnw" ## par(mfrow=c(1,2)) ## profileplot(clus,label=labs,main="Clustered data",colpal=c("heat","blue"),add.quartiles=T,fromto=c(0.05,0.95)) ################################################### ### chunk number 36: profileplot ################################################### #line 442 "vignettes/Starr/inst/doc/Starr.Rnw" png("profileplot.png", height=400, width=720) sampls = 100 probes = 63 at = (-31:31)*14 clus = matrix(rnorm(probes*sampls,sd=1),ncol=probes) clus= rbind( t(t(clus)+sin(1:probes/10))+1:nrow(clus)/sampls , t(t(clus)+sin(pi/2+1:probes/10))+1:nrow(clus)/sampls ) labs = paste("cluster",kmeans(clus,2)$cluster) par(mfrow=c(1,2)) profileplot(clus,label=labs,main="Clustered data",colpal=c("heat","blue","red","topo"),add.quartiles=T,fromto=c(0.05,0.95)) dev.off() ################################################### ### chunk number 37: ################################################### #line 467 "vignettes/Starr/inst/doc/Starr.Rnw" tssAnno <- transcriptAnno watson <- which(tssAnno$strand == 1) tssAnno[watson,]$end <- tssAnno[watson,]$start crick <- which(tssAnno$strand == -1) tssAnno[crick,]$start <- tssAnno[crick,]$end ################################################### ### chunk number 38: ################################################### #line 477 "vignettes/Starr/inst/doc/Starr.Rnw" profile <- getProfiles(rpb3_loess_ratio, probeAnnoChr1, tssAnno, 500, 500, feature="TSS", borderNames="TSS", method="basewise") ################################################### ### chunk number 39: plotProfiles eval=FALSE ################################################### ## #line 483 "vignettes/Starr/inst/doc/Starr.Rnw" ## clust <- rep(1, dim(tssAnno)[1]) ## names(clust) <- tssAnno$name ## plotProfiles(profile, cluster=clust) ################################################### ### chunk number 40: plotProfiles ################################################### #line 489 "vignettes/Starr/inst/doc/Starr.Rnw" png("sumPlot.png", height=400, width=720) clust <- rep(1, dim(tssAnno)[1]) names(clust) <- tssAnno$name plotProfiles(profile, cluster=clust, type="l", lwd=2) dev.off() ################################################### ### chunk number 41: ################################################### #line 516 "vignettes/Starr/inst/doc/Starr.Rnw" peaks <- cmarrt.ma(rpb3_loess_ratio, probeAnnoChr1, chr=NULL, M=NULL, frag.length=300) ################################################### ### chunk number 42: diagnostic plots cmarrt eval=FALSE ################################################### ## #line 528 "vignettes/Starr/inst/doc/Starr.Rnw" ## plotcmarrt(peaks) ################################################### ### chunk number 43: diagnostic plots cmarrt ################################################### #line 532 "vignettes/Starr/inst/doc/Starr.Rnw" png("cmarrt.png", height=800, width=720) plotcmarrt(peaks) dev.off() ################################################### ### chunk number 44: ################################################### #line 553 "vignettes/Starr/inst/doc/Starr.Rnw" peaklist <- cmarrt.peak(peaks, alpha = 0.05, method = "BH", minrun = 4) ################################################### ### chunk number 45: ################################################### #line 557 "vignettes/Starr/inst/doc/Starr.Rnw" str(peaklist) ################################################### ### chunk number 46: smoothing ################################################### #line 569 "vignettes/Starr/inst/doc/Starr.Rnw" rpb3_ratio_smooth <- computeRunningMedians(rpb3_loess_ratio, probeAnno=probeAnnoChr1, allChr = "chr1", winHalfSize = 80, modColumn="type") sampleNames(rpb3_ratio_smooth) <- paste(sampleNames(rpb3_loess_ratio),"smoothed") y0 <- apply(exprs(rpb3_ratio_smooth), 2, upperBoundNull) ################################################### ### chunk number 47: ChIP-enriched regions ################################################### #line 577 "vignettes/Starr/inst/doc/Starr.Rnw" distCutOff <- max(transcriptAnno$end - transcriptAnno$start) chers <- findChersOnSmoothed(rpb3_ratio_smooth, probeAnno=probeAnnoChr1, thresholds=y0, allChr="chr1", distCutOff=distCutOff, cellType="yeast", minProbesInRow = 10) ################################################### ### chunk number 48: ChIP-enriched regions ################################################### #line 584 "vignettes/Starr/inst/doc/Starr.Rnw" chers <- relateChers(chers, transcriptAnno, upstream=500) ################################################### ### chunk number 49: ChIP-enriched regions ################################################### #line 587 "vignettes/Starr/inst/doc/Starr.Rnw" chersD <- as.data.frame.cherList(chers) chersD <- chersD[which(chersD$feature != ""),] chersD[order(chersD$maxLevel, decreasing=TRUE)[1:5],] ################################################### ### chunk number 50: plotCher ################################################### #line 595 "vignettes/Starr/inst/doc/Starr.Rnw" plot(chers[[11]], rpb3_ratio_smooth, probeAnno=probeAnnoChr1, gff=transcriptAnno, paletteName="Spectral") ################################################### ### chunk number 51: sessionInfo ################################################### #line 615 "vignettes/Starr/inst/doc/Starr.Rnw" toLatex(sessionInfo())