################################################### ### chunk number 1: setup ################################################### library(ChIPseqTutorial) ################################################### ### chunk number 2: ################################################### data(ctcf, package="ChIPseqTutorial") data(gfp, package="ChIPseqTutorial") ################################################### ### chunk number 3: ################################################### ctcf gfp ################################################### ### chunk number 4: ################################################### table(strand(ctcf)) table(chromosome(gfp)) ################################################### ### chunk number 5: chromlens ################################################### library(BSgenome.Mmusculus.UCSC.mm9) mouse.chromlens <- seqlengths(Mmusculus) head(mouse.chromlens) ################################################### ### chunk number 6: coverage ################################################### cov.ctcf <- coverage(ctcf, width = mouse.chromlens, extend = 126L) cov.ctcf cov.ctcf$chr10 ################################################### ### chunk number 7: ################################################### islands <- slice(cov.ctcf$chr10, lower = 1) islands ################################################### ### chunk number 8: viewSums ################################################### viewSums(head(islands)) viewMaxs(head(islands)) nread.tab <- table(viewSums(islands) / 150L) depth.tab <- table(viewMaxs(islands)) head(nread.tab, 10) head(depth.tab, 10) ################################################### ### chunk number 9: ################################################### islandReadSummary <- function(cov) { s <- slice(cov, lower = 1) tab <- table(viewSums(s) / 150) ans <- data.frame(nread = as.numeric(names(tab)), count = as.numeric(tab)) ans } ################################################### ### chunk number 10: ################################################### head(islandReadSummary(cov.ctcf$chr10)) ################################################### ### chunk number 11: make.groups ################################################### nread.islands <- lapply(cov.ctcf, islandReadSummary) nread.islands <- do.call(make.groups, nread.islands) head(nread.islands) ################################################### ### chunk number 12: xyplot ################################################### xyplot(log(count) ~ nread | which, data = nread.islands, subset = (nread <= 20), pch = 16, type = c("p", "g")) ################################################### ### chunk number 13: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 14: ################################################### xyplot(log(count) ~ nread | which, data = nread.islands, subset = (nread <= 20), pch = 16, panel = function(x, y, ...) { panel.grid(h = -1, v = -1) panel.lmline(x[1:2], y[1:2], col = "black") panel.xyplot(x, y, ...) }) ################################################### ### chunk number 15: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 16: ################################################### islandDepthSummary <- function(cov) { s <- slice(cov, lower = 1) tab <- table(viewMaxs(s)) ans <- data.frame(depth = as.numeric(names(tab)), count = as.numeric(tab)) ans } depth.islands <- lapply(cov.ctcf, islandDepthSummary) depth.islands <- do.call(make.groups, depth.islands) xyplot(log(count) ~ depth | which, depth.islands, subset = (depth <= 20), pch = 16, panel = function(x, y, ...) { panel.grid(h = -1, v = -1) 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, ...) }) ################################################### ### chunk number 17: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 18: peaks ################################################### peaks <- slice(cov.ctcf$chr10, lower = 8) peaks ################################################### ### chunk number 19: ################################################### peak.depths <- viewMaxs(peaks) peak.areas <- viewSums(peaks) xyplot(peak.areas ~ peak.depths) ################################################### ### chunk number 20: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 21: ################################################### wpeaks <- tail(order(peak.depths), 4) peaks[wpeaks] ################################################### ### chunk number 22: ################################################### coverageplot <- function (peaks, xlab = "Position", ylab = "Coverage", ...) { pos1 <- seq(start(peaks[1]), end(peaks[1])) cov1 <- as.integer(peaks[[1]]) pos1 <- c(head(pos1, 1), pos1, tail(pos1, 1)) cov1 <- c(0, cov1, 0) xyplot(cov1 ~ pos1, ..., panel = panel.polygon, col = "lightgrey", xlab = xlab, ylab = ylab) } ################################################### ### chunk number 23: ################################################### coverageplot(peaks[wpeaks[1]]) ################################################### ### chunk number 24: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 25: all-peaks ################################################### allPeaks <- slice(cov.ctcf, lower = 8) ################################################### ### chunk number 26: ctcf-counts ################################################### ctcf.counts <- aggregate(cov.ctcf, allPeaks, sum) ################################################### ### chunk number 27: gfp-counts ################################################### cov.gfp <- coverage(gfp, width = mouse.chromlens, extend = 126L) gfp.counts <- aggregate(cov.gfp, allPeaks, sum) ################################################### ### chunk number 28: countsdf ################################################### counts <- RangedData(allPeaks, CTCF=as.integer(ctcf.counts), GFP=as.integer(gfp.counts)) ################################################### ### chunk number 29: ################################################### xyplot(asinh(CTCF)~asinh(GFP)|space, as.data.frame(counts)) ################################################### ### chunk number 30: ################################################### plot(trellis.last.object()) ################################################### ### chunk number 31: ctcf-tops ################################################### rd0 <- counts[counts[["GFP"]] == 0,]["chr10"] head(rd0[order(rd0[["CTCF"]], decreasing=TRUE),]) ################################################### ### chunk number 32: ################################################### toLatex(sessionInfo())