################################################### ### chunk number 1: setup ################################################### library(chipseq) library(lattice) ################################################### ### chunk number 2: ################################################### load("../data/alignedLocs.rda") alignedLocs ################################################### ### chunk number 3: ################################################### alignedLocs$sample1 ################################################### ### chunk number 4: ################################################### str(head(as(alignedLocs$sample1, "list"), 4)) ################################################### ### chunk number 5: eval=FALSE ################################################### library(BSgenome.Mmusculus.UCSC.mm9) mouse.chromlens <- seqlengths(Mmusculus) ################################################### ### chunk number 6: ################################################### head(mouse.chromlens) ################################################### ### chunk number 7: ################################################### ext <- extendReads(alignedLocs$sample1$chr5, seqLen = 200) head(ext) ################################################### ### chunk number 8: ################################################### cov <- coverage(ext, start = 1, end = mouse.chromlens["chr5"]) cov ################################################### ### chunk number 9: ################################################### islands <- slice(cov, lower = 1) islands ################################################### ### chunk number 10: ################################################### viewSums(head(islands)) viewMaxs(head(islands)) nread.tab <- table(viewSums(islands) / 200) depth.tab <- table(viewMaxs(islands)) head(nread.tab, 10) head(depth.tab, 10) ################################################### ### chunk number 11: ################################################### islandReadSummary <- function(x) { g <- extendReads(x, seqLen = 200) s <- slice(coverage(g, 1, max(end(g))), lower = 1) tab <- table(viewSums(s) / 200) ans <- data.frame(nread = as.numeric(names(tab)), count = as.numeric(tab)) ans } ################################################### ### chunk number 12: ################################################### head(islandReadSummary(alignedLocs$sample1$chr5)) ################################################### ### chunk number 13: ################################################### nread.islands <- gdApply(alignedLocs, islandReadSummary) nread.islands <- as(nread.islands, "data.frame") head(nread.islands) ################################################### ### chunk number 15: ################################################### xyplot(log(count) ~ nread | sample, nread.islands, subset = (chromosome == "chr5" & nread <= 40), layout = c(1, 3), pch = 16, type = c("p", "g")) ################################################### ### chunk number 17: ################################################### xyplot(log(count) ~ nread | sample, nread.islands, subset = (chromosome == "chr5" & nread <= 40), layout = c(1, 3), pch = 16, type = c("p", "g"), panel = function(x, y, ...) { panel.lmline(x[1:2], y[1:2], col = "black") panel.xyplot(x, y, ...) }) ################################################### ### chunk number 19: ################################################### islandDepthSummary <- function(x) { g <- extendReads(x, seqLen = 200) s <- slice(coverage(g, 1, max(end(g))), lower = 1) tab <- table(viewMaxs(s)) ans <- data.frame(depth = as.numeric(names(tab)), count = as.numeric(tab)) ans } depth.islands <- gdApply(alignedLocs, islandDepthSummary) depth.islands <- as(depth.islands, "data.frame") xyplot(log(count) ~ depth | sample, depth.islands, subset = (chromosome == "chr5" & depth <= 20), layout = c(1, 3), 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, ...) }) ################################################### ### chunk number 21: ################################################### peaks <- slice(cov, lower = 8) peaks ################################################### ### chunk number 22: ################################################### peak.depths <- viewMaxs(peaks) cov.pos <- coverage(extendReads(alignedLocs$sample1$chr5, strand = "+", seqLen = 200), start = 1, end = mouse.chromlens["chr5"]) cov.neg <- coverage(extendReads(alignedLocs$sample1$chr5, strand = "-", seqLen = 200), start = 1, end = mouse.chromlens["chr5"]) peaks.pos <- copyIRanges(peaks, cov.pos) peaks.neg <- copyIRanges(peaks, cov.neg) which(peak.depths >= 40) ################################################### ### chunk number 23: ################################################### plotPeak(peaks.pos[6], peaks.neg[6]) ################################################### ### chunk number 25: ################################################### plotPeak(peaks.pos[96], peaks.neg[96]) ################################################### ### chunk number 27: ################################################### plotPeak(peaks.pos[126], peaks.neg[126]) ################################################### ### chunk number 29: ################################################### plotPeak(peaks.pos[173], peaks.neg[173]) ################################################### ### chunk number 31: ################################################### extRanges <- gdApply(alignedLocs, extendReads, seqLen = 200) peakSummary <- diffPeakSummary(extRanges$sample1, extRanges$sample2, chrom.lens = mouse.chromlens, lower = 10) xyplot(asinh(sums2) ~ asinh(sums1) | chromosome, data = peakSummary, subset = (chromosome %in% c("chr1", "chr2")), panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(median(y - x), 1) }, type = c("p", "g"), alpha = 0.5, aspect = "iso") ################################################### ### chunk number 33: ################################################### peakSummary <- within(peakSummary, { diffs <- asinh(sums2) - asinh(sums1) resids <- (diffs - median(diffs)) / mad(diffs) up <- resids > 2 down <- resids < -2 }) head(peakSummary) ################################################### ### chunk number 34: ################################################### data(geneMouse) gregions <- genomic_regions(genes = geneMouse, proximal = 2000) gregions$gene <- as.character(gregions$gene) str(gregions) ################################################### ### chunk number 35: ################################################### ans <- contextDistribution(peakSummary, gregions) head(ans) ################################################### ### chunk number 36: ################################################### sumtab <- as.data.frame(rbind(total = xtabs(total ~ type, ans), promoter = xtabs(promoter ~ type, ans), threeprime = xtabs(threeprime ~ type, ans), upstream = xtabs(upstream ~ type, ans), downstream = xtabs(downstream ~ type, ans), gene = xtabs(gene ~ type, ans))) cbind(sumtab, ratio = round(sumtab[, "down"] / sumtab[, "up"], 3)) ################################################### ### chunk number 37: ################################################### sessionInfo()