################################################### ### chunk number 1: Loading library ################################################### library(Starr) ################################################### ### chunk number 2: Reading bpmap file ################################################### dataPath <- system.file("extdata", package="Starr") bpmapChr1 <- readBpmap(file.path(dataPath, "Scerevisiae_tlg_chr1.bpmap")) ################################################### ### chunk number 3: Data read-in ################################################### 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 ################################################### rpb3Chr1 ################################################### ### chunk number 5: assayData ################################################### head(exprs(rpb3Chr1)) ################################################### ### chunk number 6: phenoData ################################################### pData(rpb3Chr1) ################################################### ### chunk number 7: featureData ################################################### featureData(rpb3Chr1) head(featureData(rpb3Chr1)$chr) head(featureData(rpb3Chr1)$seq) head(featureData(rpb3Chr1)$pos) ################################################### ### chunk number 8: Reconstruction of the array image eval=FALSE ################################################### ## plotImage(file.path(dataPath,"Rpb3_IP_chr1.cel")) ################################################### ### chunk number 9: Reconstruction of the array image ################################################### 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 ################################################### ## par(mfcol=c(1,2)) ## plotDensity(rpb3Chr1, oneDevice=T, main="") ## plotBoxes(rpb3Chr1) ################################################### ### chunk number 11: boxplots and density plots ################################################### 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 ################################################### ## plotScatter(rpb3Chr1, density=T, cex=0.5) ################################################### ### chunk number 13: Scatterplot matrix ################################################### 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 ################################################### ## ips <- rpb3Chr1$type == "IP" ## controls <- rpb3Chr1$type == "CONTROL" ## plotMA(rpb3Chr1, ip=ips, control=controls) ################################################### ### chunk number 15: MA plot of raw data ################################################### 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 ################################################### ## 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 ################################################### 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: ################################################### rpb3_loess <- normalize.Probes(rpb3Chr1, method="loess") ################################################### ### chunk number 19: MA-plot of the normalized data eval=FALSE ################################################### ## plotMA(rpb3_loess, ip=ips, control=controls) ################################################### ### chunk number 20: MA-plot of the normalized data ################################################### png("maNorm.png", height=400, width=720) plotMA(rpb3_loess, ip=ips, control=controls) dev.off() ################################################### ### chunk number 21: Calculating ratio ################################################### 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 ################################################### ## 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) ################################################### 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: ################################################### probeAnnoChr1 <- bpmapToProbeAnno(bpmapChr1) ################################################### ### chunk number 25: Remapping probes ################################################### newbpmap <- remap(bpmapChr1, path=dataPath, reverse_complementary=TRUE, return_bpmap=TRUE) ################################################### ### chunk number 26: Summary of bpmap ################################################### str(newbpmap) ################################################### ### chunk number 27: Summary of bpmap eval=FALSE ################################################### ## writeTpmap("newbpmap.tpmap", newbpmap) ## tpmap2bpmap("newbpmap.tpmap", "newbpmap.bpmap") ## ## pA <- bpmapToProbeAnno(newbpmap) ################################################### ### chunk number 28: ################################################### 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 ################################################### 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 ################################################### means_rpb3 <- getMeans(rpb3_loess_ratio, probeAnnoChr1, transcriptAnno[which(transcriptAnno$name %in% filteredIDs),], info) ################################################### ### chunk number 31: correlationPlot eval=FALSE ################################################### ## 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 ################################################### 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 ################################################### ## 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 ################################################### ## labs = paste("cluster",kmeans(clus,2)$cluster) ################################################### ### chunk number 35: profileplotExampleData eval=FALSE ################################################### ## 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 ################################################### 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: ################################################### 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: ################################################### profile <- getProfiles(rpb3_loess_ratio, probeAnnoChr1, tssAnno, 500, 500, feature="TSS", borderNames="TSS", method="basewise") ################################################### ### chunk number 39: plotProfiles eval=FALSE ################################################### ## clust <- rep(1, dim(tssAnno)[1]) ## names(clust) <- tssAnno$name ## plotProfiles(profile, cluster=clust) ################################################### ### chunk number 40: plotProfiles ################################################### 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: ################################################### peaks <- cmarrt.ma(rpb3_loess_ratio, probeAnnoChr1, chr=NULL, M=NULL, frag.length=300) ################################################### ### chunk number 42: diagnostic plots cmarrt eval=FALSE ################################################### ## plotcmarrt(peaks) ################################################### ### chunk number 43: diagnostic plots cmarrt ################################################### png("cmarrt.png", height=800, width=720) plotcmarrt(peaks) dev.off() ################################################### ### chunk number 44: ################################################### peaklist <- cmarrt.peak(peaks, alpha = 0.05, method = "BH", minrun = 4) ################################################### ### chunk number 45: ################################################### str(peaklist) ################################################### ### chunk number 46: smoothing ################################################### 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 ################################################### 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 ################################################### chers <- relateChers(chers, transcriptAnno, upstream=500) ################################################### ### chunk number 49: ChIP-enriched regions ################################################### chersD <- as.data.frame.cherList(chers) chersD <- chersD[which(chersD$feature != ""),] chersD[order(chersD$maxLevel, decreasing=TRUE)[1:5],] ################################################### ### chunk number 50: plotCher ################################################### plot(chers[[11]], rpb3_ratio_smooth, probeAnno=probeAnnoChr1, gff=transcriptAnno, paletteName="Spectral") ################################################### ### chunk number 51: sessionInfo ################################################### toLatex(sessionInfo())