ChIP-seq Analysis with Bioconductor

October 22, 2012

(This article was first published on mintgene » R, and kindly contributed to R-bloggers)

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
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 input ## base calling quality of the input and chip dataset input_qa <- quality(input) chip_qa <- quality(chip) input_qa ## 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)

To leave a comment for the author, please follow the link and comment on their blog: mintgene » R. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)