ChIP-seq Analysis with Bioconductor

Often scientists are interested in finding genome-wide binding site of their protein of interest. R offers easy way to load and process the sequence files coming from ChIP-seq experiment. During the next weeks I’m going to present a pipeline that uses several Bioconductor packages which contain necessary functions for working on NGS datasets. To keep blog organised, posts are going to be divided into several parts. As content is uploaded, links will be made in this summary post for easier access.

1. Loading and quality control for sequence alignment maps
2. Filtering and coverage calculation
3. Peak finding using local lambda fit for Poisson distribution
4. Motif discovery
5. Motif scanning

So let’s start with the first topic – loading and QA of sequence reads.

## packages for loading BAM files</pre>
library(ShortRead); library(chipseq)

## loading alignment map
input <- readAligned(dirPath = "/Directory/", pattern = "Input.bam", type = "BAM")
chip <- readAligned(dirPath = "/Directory/", pattern = "Foxa2.bam", type = "BAM")

## content of AlignedRead class object

## base calling quality of the input and chip dataset
input_qa <- quality(input)
chip_qa <- quality(chip)


## convert quality encoding to integer matrix and plot the scores
input_matrix_qa <- as(input_qa[sample(1:length(input_qa), 10000)], "matrix")
chip_matrix_qa <- as(chip_qa[sample(1:length(chip_qa), 10000)], "matrix")

plot(apply(input_matrix_qa, 2, mean), pch = 16, col = "red", ylim = c(0, 40), ylab = "Phred Quality Score", xlab = "Base Position")
lines(apply(chip_matrix_qa, 2, mean), pch = 16, col = "black", type = "p")
legend("topright", c("chip", "input"), pch = 16, col = c("black", "red"), cex = 2, bty = "n")

## nucleotide frequency
sread_chip <- as.character(sread(chip)[sample(1:length(chip), 10000)])
sread_chip <-"rbind", lapply(sread_chip, function(i) strsplit(i, "")[[1]]))
sread_chip <- apply(sread_chip, 2, function(i) table(i)[c("A", "C", "G", "T")]) / 10000

barplot(sread_chip, xlim=c(0,50), ylab = "Nucleotide Distribution")
legend("right", fill=gray(1:4/4), c("A", "C", "G", "T"), bty = "n", cex = 1.5)

About these ads

One thought on “ChIP-seq Analysis with Bioconductor

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s