Period 4 - Analyses of ChIP Seq data¶
Abstract¶
ChIP-seq is a protocol for inferring the locations of proteins bound or associated with DNA. The raw data looks quite different than DNA- or RNA-seq, in that the NGS reads form tall "peaks" at the locations where the proteins were tightly bound to DNA in the cells which were used to create the sample. More specifically, ChIP-seq results in two peaks of reads of different strands (plus/minus also referred to as Watson/Crick), as shown in Figure 1 of the MACS manuscript: Zhang 2008.
Most materials come from the PH525x lab.
Approaches to peak calling and merging¶
Peak calling¶
In the first lab, we use the MACS software to call peaks. The code for this is in the MACS.txt file.
There are many different algorithms for calling peaks, which have varying performance on different kinds of experiments. As mentioned in the lecture, for ChIP of proteins with broad peaks (such as modified histones), algorithms other than those for detecting sharp peaks might perform better.
After peak calling¶
A number of analyses might be of interest following peak calling. In this lab we will focus on differential binding across samples, by focusing on the peak regions and counting the number of ChIP-seq reads which fall into the peaks for each sample.
Motif-finding is common ChIP-seq analysis which is not explored in this course, as we do not cover the basics of analysis of sequences. Motif-finding refers to the task of looking for common strings of DNA letters contained within peaks. These are biologically meaningful, as a number of proteins which are bound to DNA have conformations which make certain strings of DNA letters more preferable for binding. For more references, see the Footnotes.
Calling differential binding across samples¶
The following lab will go over the functionality of the DiffBind package, mostly using code from the vignette. This package is useful for manipulating ChIP-seq signal in R, for comparing signal across files and for performing tests of differential binding.
Reading peak files into R¶
The DiffBind package provides some called peaks that we can analyze. From the DiffBind vignette:
The dataset for this example consists of ChIPs against the transcription factor ERa using five breast cancer cell lines (C.S. Ross-Innes et al, Nature 481(7381):389–393, 2012). Three of these cell lines are responsive to tamoxifen treatment, while two others are resistant to tamoxifen. There are at least two replicates for each of the cell lines, with one cell line having three replicates, for a total of eleven sequenced libraries.
Of the five cell lines, two are based on MCF7 cells: the regular tamoxifen responsive line, as well as MCF7 cells specially treated with tamoxifen until a tamoxifen resistant cell line is obtained. For each sample, we have one peakset originally derived using the MACS peak caller (Y. Zhang et al, Genome Biol, 9(9):R137, 2008), for a total of eleven peaksets. Note that to save space in the package, only data for chromosome 18 is used for the vignette.
We're going to take this opportunity to demonstrate the use of BiocFileCache to organize and keep track of files:
suppressPackageStartupMessages({
library(DiffBind)
library(BiocFileCache)
})
bfc <- BiocFileCache(cache="~/chipseq-data")
DBfiles <- list.files(system.file("extra", package="DiffBind"),
recursive = TRUE, full.names = TRUE)
DBfiles <- DBfiles[!DBfiles %in% bfcinfo(bfc)$rname]
for (i in seq_along(DBfiles))
bfcadd(bfc, rname=DBfiles[i], rtype="local", action="copy")
Peaks are represented as a .bed file per sample, in this example summarized in tamoxifen.csv:
tamfile <- bfcquery(bfc, "tamoxifen.csv")$fpath
read.csv(tamfile)
Just as a note, we now have the paths of all these files stored in our BiocFileCache:
bfcquery(bfc, "peaks")$rpath
Peak merging and differential binding setup¶
The dba function creates the basic object for an analysis of Differential Binding Affinity. The sample sheet specifies a data frame of file with certain required columns. Note that columns have restricted names, including Tissue, Factor, Condition, etc., which will be referred to later in analysis.
This function will automatically create a correlation plot showing the overlap of the peaks for all the samples.
setwd(system.file("extra", package="DiffBind")) #necessary because `tamfile` contains relative paths
ta <- dba(sampleSheet=tamfile)
ta
From the DiffBind vignette, we have:
This shows how many peaks are in each peakset, as well as (in the first line) the total number of unique peaks after merging overlapping ones (3795), and the dimensions of dba.plotPCA (default binding matrix of 11 samples by the 2845 sites that overlap in at least two of the samples).
We can access the peaks for each file:
names(ta)
class(ta$peaks)
head(ta$peaks[[1]])
Differential binding¶
The following code chunk will count the reads from the BAM files specified in the samples slot:
ta$samples
options(repr.plot.width=5, repr.plot.height=5)
# the next line does not actually work, because the BAM files are not included in the package
# ta <- dba.count(ta, minOverlap=3)
# instead we load the counts:
data(tamoxifen_counts)
plot(tamoxifen)
We can perform a test by specifying to contrast over the levels of condition. This will call DESeq software in order to normalize samples for sequencing depth and perform essentially the same analysis as a differential expression analysis for RNA-Seq counts:
tamoxifen$config$AnalysisMethod
The plot produced then looks at correlation only for those peaks which showed evidence of differential binding.
ta2 <- dba.contrast(tamoxifen, categories=DBA_CONDITION)
ta2 <- dba.analyze(ta2)
ta2
Note: We could have included the tissue as a blocking factor, by providing DBA_TISSUE to the block argument of dba.contrast.
From the DiffBind vignette, we have:
By default, dba.analyze plots a correlation heatmap if it finds any significantly differentially bound sites, shown in Figure 3. Using only the differentially bound sites, we now see that the four tamoxifen resistant samples (representing two cell lines) cluster together, although the tamoxifen-responsive MCF7 replicates cluster closer to them than to the other tamoxifen responsive samples."
Finally, we can generate the results table, which is attached as metadata columns to the peaks as genomic ranges. By specifying bCounts = TRUE, we also obtain the normalized counts for each sample.
tadb <- dba.report(ta2)
tadb
counts <- dba.report(ta2, bCounts=TRUE)
Annotating and plotting peaks¶
First, do any differentially bound peaks overlap with genes? Yes, 299 of the 629 peaks do:
library(Homo.sapiens)
gn <- genes(Homo.sapiens, columns="SYMBOL")
summary(counts %over% gn)
How many genes does each peak overlap with?
table(countOverlaps(counts, gn))
Let's look at those peaks that overlap with two genes:
count2 <- counts[countOverlaps(counts, gn) == 2]
gn2 <- gn[gn %over% count2]
width(gn2) / 1e3 #width in kb
gn2 <- gn2[order(ranges(gn2))]
gn2
Plot peaks in the UCSC genome browser¶
Although the DiffBind package doesn't specify the genome, let's assume it is hg19, add a track to the hg19 genome in the UCSC genome browser, then start the browser centered on the peak with greatest fold-change:
genome(counts) <- "hg19"
library(rtracklayer)
session <- browserSession("UCSC")
genome(session) <- "hg19"
track(session, "counts") <- counts
rangemaxFC <- counts[which.max(abs(counts$Fold))]
browserView(session, range=rangemaxFC * 0.75) #0.75 zoom-factor
Finally, let's plot the peak that overlapped with AQP4-AS1 on the sense strand and CHST9 on the antisense strand:
gn2[3:4]
(plotrange <- reduce(gn2[3:4], ignore.strand=TRUE))
browserView(session, range=plotrange * 0.75)
Annotating peaks¶
The ChIPseeker library makes it easy to annotate peaks by overlapping genes, UTRs, promoters, ...
library(ChIPseeker)
peakAnno <- annotatePeak(counts,
TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene)
peakAnno
options(repr.plot.width=6, repr.plot.height=5)
plotAnnoPie(peakAnno)
Exercise: Look at the ChIPseeker vignette and perform additional visualizations of the peak annotations
Blacklisted regions¶
Certain genomic regions are known to produce artefactual ChIP-seq peaks, and are most safely removed or ignored if a decent blacklist is available. For human, a good blacklist is available from ENCODE as a .bed file.
Exercise: Import the ENCODE blacklist file using rtracklayer. Do any of the peaks in our counts object overlap with blacklisted regions? The file import code is provided below.
suppressPackageStartupMessages(library(rtracklayer))
bl <- import("https://www.encodeproject.org/files/ENCFF419RSJ/@@download/ENCFF419RSJ.bed.gz", genome="hg19")
Footnotes ¶
Model-based Analysis for ChIP-Seq (MACS)¶
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. "Model-based Analysis of ChIP-Seq (MACS)". Genome Biol. 2008. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2592715/
Software:
http://liulab.dfci.harvard.edu/MACS/
Motif finding¶
Wikipedia's article on DNA sequence motifs: http://en.wikipedia.org/wiki/Sequence_motif
A non-comprehensive list of software for motif finding:
A survey of motif finding algorithms: http://www.biomedcentral.com/1471-2105/8/S7/S21/
GenomicRanges cheat-sheet¶
