The following introduces several utilities useful for ChIP-Seq data. They are not part of the actual workflow.

Rle object stores coverage information

library(rtracklayer); library(GenomicRanges); library(Rsamtools); library(GenomicAlignments)
aligns <- readGAlignments(outpaths(args)[1])
cov <- coverage(aligns)
cov

Resizing aligned reads

trim(resize(as(aligns, "GRanges"), width = 200))

Naive peak calling

islands <- slice(cov, lower = 15)
islands[[1]]

Plot coverage for defined region

library(ggbio)
myloc <- c("Chr1", 1, 100000)
ga <- readGAlignments(outpaths(args)[1], use.names=TRUE, param=ScanBamParam(which=GRanges(myloc[1], IRanges(as.numeric(myloc[2]), as.numeric(myloc[3])))))
autoplot(ga, aes(color = strand, fill = strand), facets = strand ~ seqnames, stat = "coverage")
Jump to: next_page