### R code from vignette source 'vignettes/chipseq/inst/doc/Workflow.Rnw' ################################################### ### code chunk number 1: setup ################################################### library(chipseq) library(GenomicFeatures) library(lattice) ################################################### ### code chunk number 2: preprocess (eval = FALSE) ################################################### ## qa_list <- lapply(sampleFiles, qa) ## report(do.call(rbind, qa_list)) ## ## spend some time evaluating the QA report, then procede ## filter <- compose(chipseqFilter(), alignQualityFilter(15)) ## cstest <- seqapply(sampleFiles, function(file) { ## as(readAligned(file, filter), "GRanges") ## }) ## cstest <- cstest[seqnames(cstest) %in% c("chr10", "chr11", "chr12")] ################################################### ### code chunk number 3: Workflow.Rnw:85-87 ################################################### data(cstest) cstest ################################################### ### code chunk number 4: convert-cstest (eval = FALSE) ################################################### ## ## code used to convert the GenomeDataList to a GRangesList ## cstest <- seqapply(cstest, function(gd) { ## gr <- do.call(c, lapply(names(gd), function(chr) { ## pos <- gd[[chr]] ## starts <- c(pos[["-"]] - 23L, pos[["+"]]) ## GRanges(chr, IRanges(starts, width = 24), ## rep(c("-", "+"), elementLengths(pos))) ## })) ## }) ################################################### ### code chunk number 5: Workflow.Rnw:105-106 ################################################### cstest$ctcf ################################################### ### code chunk number 6: Workflow.Rnw:115-116 ################################################### library(BSgenome.Mmusculus.UCSC.mm9) ################################################### ### code chunk number 7: estimate.mean.fraglen ################################################### fraglen <- estimate.mean.fraglen(cstest$ctcf) fraglen[!is.na(fraglen)] ################################################### ### code chunk number 8: Workflow.Rnw:151-153 ################################################### ctcf.ext <- resize(cstest$ctcf, width = 200) ctcf.ext ################################################### ### code chunk number 9: Workflow.Rnw:164-166 ################################################### cov.ctcf <- coverage(ctcf.ext) cov.ctcf ################################################### ### code chunk number 10: Workflow.Rnw:176-178 ################################################### islands <- slice(cov.ctcf, lower = 1) islands ################################################### ### code chunk number 11: Workflow.Rnw:183-191 ################################################### viewSums(islands) viewMaxs(islands) nread.tab <- table(viewSums(islands) / 200) depth.tab <- table(viewMaxs(islands)) nread.tab[,1:10] depth.tab[,1:10] ################################################### ### code chunk number 12: Workflow.Rnw:207-217 ################################################### islandReadSummary <- function(x) { g <- resize(x, 200) s <- slice(coverage(g), lower = 1) tab <- table(viewSums(s) / 200) df <- DataFrame(tab) colnames(df) <- c("chromosome", "nread", "count") df$nread <- as.integer(df$nread) df } ################################################### ### code chunk number 13: Workflow.Rnw:221-222 ################################################### head(islandReadSummary(cstest$ctcf)) ################################################### ### code chunk number 14: Workflow.Rnw:227-230 ################################################### nread.islands <- seqapply(cstest, islandReadSummary) nread.islands <- stack(nread.islands, "sample") nread.islands ################################################### ### code chunk number 15: Workflow.Rnw:236-239 ################################################### xyplot(log(count) ~ nread | sample, as.data.frame(nread.islands), subset = (chromosome == "chr10" & nread <= 40), layout = c(1, 2), pch = 16, type = c("p", "g")) ################################################### ### code chunk number 16: Workflow.Rnw:242-243 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 17: Workflow.Rnw:259-266 ################################################### xyplot(log(count) ~ nread | sample, as.data.frame(nread.islands), subset = (chromosome == "chr10" & nread <= 40), layout = c(1, 2), pch = 16, type = c("p", "g"), panel = function(x, y, ...) { panel.lmline(x[1:2], y[1:2], col = "black") panel.xyplot(x, y, ...) }) ################################################### ### code chunk number 18: Workflow.Rnw:269-270 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 19: Workflow.Rnw:277-306 ################################################### islandDepthSummary <- function(x) { g <- resize(x, 200) s <- slice(coverage(g), lower = 1) tab <- table(viewMaxs(s) / 200) df <- DataFrame(tab) colnames(df) <- c("chromosome", "depth", "count") df$depth <- as.integer(df$depth) df } depth.islands <- seqapply(cstest, islandDepthSummary) depth.islands <- stack(depth.islands, "sample") xyplot(log(count) ~ depth | sample, as.data.frame(depth.islands), subset = (chromosome == "chr10" & depth <= 20), layout = c(1, 2), pch = 16, type = c("p", "g"), panel = function(x, y, ...) { lambda <- 2 * exp(y[2]) / exp(y[1]) null.est <- function(xx) { xx * log(lambda) - lambda - lgamma(xx + 1) } log.N.hat <- null.est(1) - y[1] panel.lines(1:10, -log.N.hat + null.est(1:10), col = "black") panel.xyplot(x, y, ...) }) ## depth.islands <- summarizeReads(cstest, summary.fun = islandDepthSummary) ################################################### ### code chunk number 20: Workflow.Rnw:309-310 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 21: islandDepthPlot ################################################### islandDepthPlot(cov.ctcf) ################################################### ### code chunk number 22: peakCutoff ################################################### peakCutoff(cov.ctcf, fdr = 0.0001) ################################################### ### code chunk number 23: Workflow.Rnw:337-339 ################################################### peaks.ctcf <- slice(cov.ctcf, lower = 8) peaks.ctcf ################################################### ### code chunk number 24: peakSummary ################################################### peaks <- peakSummary(peaks.ctcf) ################################################### ### code chunk number 25: Workflow.Rnw:357-368 ################################################### peak.depths <- viewMaxs(peaks.ctcf) cov.pos <- coverage(ctcf.ext[strand(ctcf.ext) == "+"]) cov.neg <- coverage(ctcf.ext[strand(ctcf.ext) == "-"]) peaks.pos <- Views(cov.pos, ranges(peaks.ctcf)) peaks.neg <- Views(cov.neg, ranges(peaks.ctcf)) wpeaks <- tail(order(peak.depths$chr10), 4) wpeaks ################################################### ### code chunk number 26: Workflow.Rnw:375-376 ################################################### coverageplot(peaks.pos$chr10[wpeaks[1]], peaks.neg$chr10[wpeaks[1]]) ################################################### ### code chunk number 27: Workflow.Rnw:378-379 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 28: Workflow.Rnw:382-383 ################################################### coverageplot(peaks.pos$chr10[wpeaks[2]], peaks.neg$chr10[wpeaks[2]]) ################################################### ### code chunk number 29: Workflow.Rnw:385-386 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 30: Workflow.Rnw:391-392 ################################################### coverageplot(peaks.pos$chr10[wpeaks[3]], peaks.neg$chr10[wpeaks[3]]) ################################################### ### code chunk number 31: Workflow.Rnw:394-395 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 32: Workflow.Rnw:398-399 ################################################### coverageplot(peaks.pos$chr10[wpeaks[4]], peaks.neg$chr10[wpeaks[4]]) ################################################### ### code chunk number 33: Workflow.Rnw:401-402 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 34: Workflow.Rnw:414-429 ################################################### ## find peaks for GFP control cov.gfp <- coverage(resize(cstest$gfp, 200)) peaks.gfp <- slice(cov.gfp, lower = 8) peakSummary <- diffPeakSummary(peaks.gfp, peaks.ctcf) xyplot(asinh(sums2) ~ asinh(sums1) | space, data = as.data.frame(peakSummary), panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(median(y - x), 1) }, type = c("p", "g"), alpha = 0.5, aspect = "iso") ################################################### ### code chunk number 35: Workflow.Rnw:431-432 ################################################### plot(trellis.last.object()) ################################################### ### code chunk number 36: Workflow.Rnw:437-446 ################################################### peakSummary <- within(peakSummary, { diffs <- asinh(sums2) - asinh(sums1) resids <- (diffs - median(diffs)) / mad(diffs) up <- resids > 2 down <- resids < -2 change <- ifelse(up, "up", ifelse(down, "down", "flat")) }) ################################################### ### code chunk number 37: Workflow.Rnw:462-465 ################################################### library(TxDb.Mmusculus.UCSC.mm9.knownGene) gregions <- transcripts(TxDb.Mmusculus.UCSC.mm9.knownGene) gregions ################################################### ### code chunk number 38: Workflow.Rnw:469-470 ################################################### promoters <- flank(gregions, 1000, both = TRUE) ################################################### ### code chunk number 39: Workflow.Rnw:474-476 ################################################### peakSummary$inPromoter <- peakSummary %over% promoters xtabs(~ inPromoter + change, peakSummary) ################################################### ### code chunk number 40: Workflow.Rnw:479-481 ################################################### peakSummary$inUpstream <- peakSummary %over% flank(gregions, 20000) peakSummary$inGene <- peakSummary %over% gregions ################################################### ### code chunk number 41: Workflow.Rnw:484-492 ################################################### sumtab <- as.data.frame(rbind(total = xtabs(~ change, peakSummary), promoter = xtabs(~ change, subset(peakSummary, inPromoter)), upstream = xtabs(~ change, subset(peakSummary, inUpstream)), gene = xtabs(~ change, subset(peakSummary, inGene)))) ##cbind(sumtab, ratio = round(sumtab[, "down"] / sumtab[, "up"], 3)) ################################################### ### code chunk number 42: rtracklayer-upload (eval = FALSE) ################################################### ## library(rtracklayer) ## session <- browserSession() ## genome(session) <- "mm9" ## session$gfpCov <- cov.gfp ## session$gfpPeaks <- peaks.gfp ## session$ctcfCov <- cov.ctcf ## session$ctcfPeaks <- peaks.ctcf ################################################### ### code chunk number 43: rtracklayer-view (eval = FALSE) ################################################### ## peak.ord <- order(unlist(peak.depths), decreasing=TRUE) ## peak.sort <- as(peaks.ctcf, "GRanges")[peak.ord] ## view <- browserView(session, peak.sort[1], full = c("gfpCov", "ctcfCov")) ################################################### ### code chunk number 44: rtracklayer-view-5 (eval = FALSE) ################################################### ## views <- browserView(session, head(peak.sort, 5), full = c("gfpCov", "ctcfCov")) ################################################### ### code chunk number 45: Workflow.Rnw:534-535 ################################################### sessionInfo()