/ChIP-seq-analysis

ChIP-seq analysis notes from Tommy Tang

Primary LanguagePythonMIT LicenseMIT

ChIP-seq-analysis

A pipeline

I developed a Snakemake based ChIP-seq pipeline: pyflow-ChIPseq.

Resources for ChIP-seq

  1. ENCODE: Encyclopedia of DNA Elements ENCODExplorer: A compilation of metadata from ENCODE. A bioc package to access the meta data of ENCODE and download the raw files.
  2. ENCODE Factorbook
  3. ChromNet ChIP-seq interactions
    paper: Learning the human chromatin network using all ENCODE ChIP-seq datasets
  4. The International Human Epigenome Consortium (IHEC) epigenome data portal
  5. GEO. Sequences are in .sra format, need to use sratools to dump into fastq.
  6. European Nucleotide Archive. Sequences are available in fastq format.
  7. Data bases and software from Sheirly Liu's lab at Harvard
  8. Blueprint epigenome
  9. A collection of tools and papers for nucelosome positioning and TF ChIP-seq
  10. review paper:Deciphering ENCODE
  11. EpiFactors is a database for epigenetic factors, corresponding genes and products.
  12. biostar handbook. My ChIP-seq chapter is out April 2017!
  13. ReMap 2018 An integrative ChIP-seq analysis of regulatory regions. The ReMap atlas consits of 80 million peaks from 485 transcription factors (TFs), transcription coactivators (TCAs) and chromatin-remodeling factors (CRFs) from public data sets. The atlas is available to browse or download either for a given TF or cell line, or for the entire dataset.

Papers on ChIP-seq

  1. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia

  2. Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data

  3. Systematic evaluation of factors influencing ChIP-seq fidelity

  4. ChIP–seq: advantages and challenges of a maturing technology

  5. ChIP–seq and beyond: new and improved methodologies to detect and characterize protein–DNA interactions

  6. Beyond library size: a field guide to NGS normalization

  7. ENCODE paper portol

  8. Enhancer discovery and characterization

  9. 2016 review Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation

  10. bioinformatics paper:Features that define the best ChIP-seq peak calling algorithms compares different peak callers for TFs and histones.

  11. Systematic comparison of monoclonal versus polyclonal antibodies for mapping histone modifications by ChIP-seq The binding patterns for H3K27ac differed substantially between polyclonal and monoclonal antibodies. However, this was most likely due to the distinct immunogen used rather than the clonality of the antibody. Altogether, we found that monoclonal antibodies as a class perform as well as polyclonal antibodies. Accordingly, we recommend the use of monoclonal antibodies in ChIP-seq experiments.

  12. A nice small review: Unraveling the 3D genome: genomics tools for multiscale exploration

  13. Three very interesting papers, Developmental biology: Panoramic views of the early epigenome

    Protocols

  14. A computational pipeline for comparative ChIP-seq analyses

  15. Identifying ChIP-seq enrichment using MACS

  16. Spatial clustering for identification of ChIP-enriched regions (SICER) to map regions of histone methylation patterns in embryonic stem cells

  17. ENCODE tutorials

  18. A User's Guide to the Encyclopedia of DNA Elements (ENCODE)

  19. A toolbox of immunoprecipitation-grade monoclonal antibodies to human transcription factors The data portal https://proteincapture.org/

Quality Control

Data downloaded from GEO usually are raw fastq files. One needs to do quality control (QC) on them.

  • fastqc
  • multiqc Aggregate results from bioinformatics analyses across many samples into a single report. Could be very useful to summarize the QC report.

Peak calling

Be careful with the peaks you get:
Active promoters give rise to false positive ‘Phantom Peaks’ in ChIP-seq experiments

It is good to have controls for your ChIP-seq experiments. A DNA input control (no antibody is applied) is prefered. The IgG control is also fine, but because so little DNA is there, you might get many duplicated reads due to PCR artifact.

For cancer cells, an input control can be used to correct for copy-number bias.

A quote from Tao Liu: who develped MACS1/2

I remember in a PloS One paper last year by Elizabeth G. Wilbanks et al., authors pointed out the best way to sort results in MACS is by -10*log10(pvalue) then fold enrichment. I agree with them. You don't have to worry about FDR too much if your input data are far more than ChIP data. MACS1.4 calculates FDR by swapping samples, so if your input signal has some strong bias somewhere in the genome, your FDR result would be bad. Bad FDR may mean something but it's just secondary.

  1. The most popular peak caller by Tao Liu: MACS2. Now --broad flag supports broad peaks calling as well.

  2. TF ChIP-seq peak calling using the Irreproducibility Discovery Rate (IDR) framework and many Software Tools Used to Create the ENCODE Resource

  3. SICER for broad histone modification ChIP-seq

  4. HOMER can also used to call Transcription factor ChIP-seq peaks and histone modification ChIP-seq peaks.

  5. MUSIC

  6. permseq R package for mapping protein-DNA interactions in highly repetitive regions of the genomes with prior-enhanced read mapping. Paper on PLos Comp.

  7. Ritornello: High fidelity control-free chip-seq peak calling. No input is required!

  8. Tumor samples are heterogeneous containing different cell types. MixChIP: a probabilistic method for cell type specific protein-DNA binding analysis

  9. Detecting broad domains and narrow peaks in ChIP-seq data with hiddenDomains tool

  10. BroadPeak: a novel algorithm for identifying broad peaks in diffuse ChIP-seq datasets

  11. epic: diffuse domain ChIP-Seq caller based on SICER. It is a re-writen of SICER for faster processing using more CPUs. (Will try it for broad peak for sure).

  12. Cistrome: The best place for wet lab scientist to check the binding sites. Developed by Shierly Liu lab in Harvard.

  13. Accounting for GC-content bias reduces systematic errors and batch effects in ChIP-Seq peak callers tool in github

  14. [SUPERmerge]((https://github.com/Bohdan-Khomtchouk/SUPERmerge):ChIP-seq coverage island analysis algorithm for broad histone marks

Different parameters using the same program can produce drastic different sets of peaks especially for histone modifications with variable enrichment length and gaps between peaks. One needs to make a valid argument for parameters he uses

An example of different parameters for homer findPeaks:

Binding does not infer functionality

The Functional Consequences of Variation in Transcription Factor Binding

" On average, 14.7% of genes bound by a factor were differentially expressed following the knockdown of that factor, suggesting that most interactions between TF and chromatin do not result in measurable changes in gene expression levels of putative target genes. "

We analyzed the dependence of the ChIP signal on the duration of formaldehyde cross-linking time for two proteins: DNA topoisomerase 1 (Top1) that is functionally associated with the double helix in vivo, especially with active chromatin, and green fluorescent protein (GFP) that has no known bona fide interactions with DNA. With short time of formaldehyde fixation, only Top1 immunoprecipation efficiently recovered DNA from active promoters, whereas prolonged fixation augmented non-specific recovery of GFP dramatizing the need to optimize ChIP protocols to minimize the time of cross-linking, especially for abundant nuclear proteins. Thus, ChIP is a powerful approach to study the localization of protein on the genome when care is taken to manage potential artifacts.

Gene set enrichment analysis for ChIP-seq peaks

The Gene Ontology Handbook Read it for basics for GO.

  1. Broad Enrich
  2. ChIP Enrich
  3. GREAT predicts functions of cis-regulatory regions.
  4. ENCODE ChIP-seq significance tool. Given a list of genes, co-regulating TFs will be identified.
  5. cscan similar to the ENCODE significance tool.
  6. CompGO: an R package for comparing and visualizing Gene Ontology enrichment differences between DNA binding experiments
  7. interactive and collaborative HTML5 gene list enrichment analysis tool
  8. GeNets from Broad. Looks very promising.
  9. Bioconductor EnrichmentBrowser
  10. clusterProfiler by Guangchuan Yu, the author of ChIPseeker.
  11. fgsea bioconductor package Fast Gene Set Entrichment Analysis.
  12. paper: A Comparison of Gene Set Analysis Methods in Terms of Sensitivity, Prioritization and Specificity

Chromatin state Segmentation

  1. ChromHMM from Manolis Kellis in MIT.

In ChromHMM the raw reads are assigned to non-overlapping bins of 200 bps and a sample-specific threshold is used to transform the count data to binary values

  1. Segway from Hoffman lab. Base pair resolution. Takes longer time to run.
  2. epicseg published 2015 in genome biology. Similiar speed with ChromHMM.
  3. Spectacle: fast chromatin state annotation using spectral learning. Also published 2015 in genome biology.
  4. chromstaR: Tracking combinatorial chromatin state dynamics in space and time
  5. epilogos visualization and analysis of chromatin state model data.
  6. Accurate promoter and enhancer identification in 127 ENCODE and Roadmap Epigenomics cell types and tissues by GenoSTAN
  7. StatePaintR StateHub-StatePaintR: rules-based chromatin state annotations.
  8. [IDEAS(https://github.com/yuzhang123/IDEAS/): an integrative and discriminative epigenome annotation system http://sites.stat.psu.edu/~yzz2/IDEAS/

deep learning in ChIP-seq

  • Coda uses convolutional neural networks to learn a mapping from noisy to high-quality ChIP-seq data. These trained networks can then be used to remove noise and improve the quality of new ChIP-seq data. From Ashul lab.
  • DeepChrome is a unified CNN framework that automatically learns combinatorial interactions among histone modification marks to predict the gene expression. (Is it really better than a simple linear model?)
  • deep learning in biology

Peak annotation

  1. Homer annotatePeak
  2. Bioconductor package ChIPseeker by Guangchuan Yu
    See an important post by him on 0 or 1 based coordinates.

Most of the software for ChIP annotation doesn't considered this issue when annotating peak (0-based) to transcript (1-based). To my knowledge, only HOMER consider this issue. After I figure this out, I have updated ChIPseeker (version >= 1.4.3) to fix the issue.

  1. Bioconductor package ChIPpeakAnno. There is a bug with this package, not sure if it is solved or not. Still a post from Guangchuan Yu: Bug of R package ChIPpeakAnno.

I used R package ChIPpeakAnno for annotating peaks, and found that it handle the DNA strand in the wrong way. Maybe the developers were from the computer science but not biology background.

  1. annotatr Annotation of Genomic Regions to Genomic Annotations.

  2. geneXtendeR computes optimal gene extensions tailored to the broadness of the specific epigenetic mark (e.g., H3K9me1, H3K27me3), as determined by a user-supplied ChIP-seq peak input file. As such, geneXtender maximizes the signal-to-noise ratio of locating genes closest to and directly under peaks

Differential peak detection

Look at a post and here describing different tools. A review paper A comprehensive comparison of tools for differential ChIP-seq analysis

  1. MultiGPS

  2. PePr. It can also call peaks.

  3. histoneHMM

  4. diffreps for histone. developed by Shen Li's lab in Mount Sinai who also develped ngs.plot.

  5. diffbind bioconductor package. Internally uses RNA-seq tools: EdgR or DESeq. Most likely, I will use this tool.

  6. ChIPComp. Very little tutorial. Now it is on bioconductor.

  7. csaw bioconductor package. Tutorial here

  8. chromDiff. Also from from Manolis Kellis in MIT. Similar with ChromHMM, documentation is not that detailed. Will have a try on this.

  9. MACS2 can detect differential peaks as well

  10. paper Identifying differential transcription factor binding in ChIP-seq

Motif enrichment

  1. HOMER. It has really detailed documentation. It can also be used to call peaks.

For TF ChIP-seq, one can usually find the summit of the peak (macs14 will report the summit), and extend the summit to both sides to 100bp-500bp. One can then use those 100bp-500 bp small regions to do motif analysis. Usually, oen should find the motif for the ChIPed TF in the ChIP-seq experiment if it is a DNA binding protein.

It is trickier to do motif analysis using histone modification ChIP-seq. For example, the average peak size of H3K27ac is 2~3 kb. If one wants to find TF binding motifs from H3K27ac ChIP-seq data, it is good to narrow down the region a bit. MEME and many other motif finding tools require that the DNA sequence length to be small (~500bp). One way is to use findPeaks in homer turning on -nfr(nucleosome free region) flag, and then do motif analysis in those regions.

suggestions for finding motifs from histone modification ChIP-seq data from HOMER page:

Since you are looking at a region, you do not necessarily want to center the peak on the specific position with the highest tag density, which may be at the edge of the region. Besides, in the case of histone modifications at enhancers, the highest signal will usually be found on nucleosomes surrounding the center of the enhancer, which is where the functional sequences and transcription factor binding sites reside. Consider H3K4me marks surrounding distal PU.1 transcription factor peaks. Typically, adding the -center >option moves peaks further away from the functional sequence in these scenarios.

Other strategy similar to -nfr was developed in this paper: Dissecting neural differentiation regulatory networks through epigenetic footprinting. In the method part of the paper, the authors computed a depletion score within the peaks, and use the footprinted regions to do motif analysis. (Thanks kadir for pointing out the paper)

http://homer.ucsd.edu/homer/ngs/peakMotifs.html

Region Size ("-size <#>", "-size <#>,<#>", "-size given", default: 200) The size of the region used for motif finding is important. If analyzing ChIP-Seq peaks from a transcription factor, Chuck would recommend 50 bp for establishing the primary motif bound by a given transcription factor and 200 bp for finding both primary and "co-enriched" motifs for a transcription factor. When looking at histone marked regions, 500-1000 bp is probably a good idea (i.e. H3K4me or H3/H4 acetylated regions). In theory, HOMER can work with very large regions (i.e. 10kb), but with the larger the regions comes more sequence and longer execution time. These regions will be based off the center of the peaks. If you prefer an offset, you can specify "-size -300,100" to search a region of size 400 that is centered 100 bp upstream of the peak center (useful if doing motif finding on putative TSS regions). If you have variable length regions, use the option "-size given" and HOMER will use the exact regions that were used as input.

I just found PARE. PARE is a computational method to Predict Active Regulatory Elements, specifically enhancers and promoters. H3K27ac and H3K4me can be used to define active enhancers.

  1. MEME suite. It is probably the most popular motif finding tool in the papers. protocol:Motif-based analysis of large nucleotide data sets using MEME-ChIP
  2. JASPAR database
  3. pScan-ChIP
  4. MotifMap
  5. RAST Regulatory Sequence Analysis Tools.
  6. ENCODE TF motif database
  7. oPOSSUM is a web-based system for the detection of over-represented conserved transcription factor binding sites and binding site combinations in sets of genes or sequences.
  8. my post how to get a genome-wide motif bed file
  9. Many other tools here
  10. A review of ensemble methods for de novo motif discovery in ChIP-Seq data
  11. melina2. If you only have one sequence and want to know what TFs might bind there, this is a very useful tool.
  12. STEME. A python library for motif analysis. STEME started life as an approximation to the Expectation-Maximisation algorithm for the type of model used in motif finders such as MEME. STEME’s EM approximation runs an order of magnitude more quickly than the MEME implementation for typical parameter settings. STEME has now developed into a fully-fledged motif finder in its own right.
  13. CENTIPEDE: Transcription factor footprinting and binding site prediction. Tutorial
  14. msCentipede: Modeling Heterogeneity across Genomic Sites and Replicates Improves Accuracy in the Inference of Transcription Factor Binding
  15. DiffLogo: A comparative visualisation of sequence motifs
  16. Weeder (version: 2.0)
  17. MCAST: scanning for cis-regulatory motif clusters Part of MEME suite.
  18. Sequence-based Discovery of Regulons iRegulon detects the TF, the targets and the motifs/tracks from a set of genes.
  19. Regulatory genomic toolbox
  20. Parse TF motifs from public databases, read into R, and scan using 'rtfbs'
  21. Romulus: Robust multi-state identification of transcription factor binding sites from DNase-seq data: Romulus is a computational method to accurately identify individual transcription factor binding sites from genome sequence information and cell-type--specific experimental data, such as DNase-seq. It combines the strengths of its predecessors, CENTIPEDE and Wellington, while keeping the number of free parameters in the model robustly low. The method is unique in allowing for multiple binding states for a single transcription factor, differing in their cut profile and overall number of DNase I cuts.
  22. moca: Tool for motif conservation analysis.
  23. gimmemotifs Suite of motif tools, including a motif prediction pipeline for ChIP-seq experiments. looks very useful, will take a look!

Super-enhancer identification

The fancy "supper-enhancer" term was first introduced by Richard Young in Whitehead Institute. Basically, super-enhancers are enhancers that span large genomic regions(~12.5kb). The concept of super-enhancer is not new. One of the most famous example is the Locus Control Region (LCR) that controls the globin gene expression, and this has been known for decades.

A review in Nature Genetics What are super-enhancers?

paper: Genetic dissection of the α-globin super-enhancer in vivo

By generating a series of mouse models, deleting each of the five regulatory elements of the α-globin super-enhancer individually and in informative combinations, we demonstrate that each constituent enhancer seems to act independently and in an additive fashion with respect to hematological phenotype, gene expression, chromatin structure and chromosome conformation, without clear evidence of synergistic or higher-order effects.

paper: Hierarchy within the mammary STAT5-driven Wap super-enhancer
paper: Enhancers and super-enhancers have an equivalent regulatory role in embryonic stem cells through regulation of single or multiple genes

From the HOMER page How finding super enhancers works:

Super enhancer discovery in HOMER emulates the original strategy used by the Young lab. First, peaks are found just like any other ChIP-Seq data set. Then, peaks found within a given distance are 'stitched' together into larger regions (by default this is set at 12.5 kb). The super enhancer signal of each of these regions is then determined by the total normalized number reads minus the number of normalized reads in the input. These regions are then sorted by their score, normalized to the highest score and the number of putative enhancer regions, and then super enhancers are identified as regions past the point where the slope is greater than 1.

Example of a super enhancer plot:

In the plot above, all of the peaks past 0.95 or so would be considered "super enhancers", while the one's below would be "typical" enhancers. If the slope threshold of 1 seems arbitrary to you, well... it is! This part is probably the 'weakest link' in the super enhancer definition. However, the concept is still very useful. Please keep in mind that most enhancers probably fall on a continuum between typical and super enhancer status, so don't bother fighting over the precise number of super enhancers in a given sample and instead look for useful trends in the data.

Using ROSE from Young lab
ROSE: RANK ORDERING OF SUPER-ENHANCERS

**imPROSE - Integrated Methods for Prediction of Super-Enhancers

Bedgraph, bigwig manipulation tools

WiggleTools
bigwig tool
bigwig-python
samtools
bedtools my all-time favorite tool from Araon Quinlan' lab. Great documentation!
Hosting bigWig for UCSC visualization
My first play with GRO-seq data, from sam to bedgraph for visualization
convert bam file to bigwig file and visualize in UCSC genome browser in a Box (GBiB)

Peaks overlapping significance test

The genomic association tester (GAT)
poverlap from Brent Pedersen. Now he is working with Aaron Quinlan at university of Utah.
Genometric Correlation (GenometriCorr): an R package for spatial correlation of genome-wide interval datasets
Location overlap analysis for enrichment of genomic ranges bioconductor package.
regioneR Association analysis of genomic regions based on permutation tests similaRpeak: Metrics to estimate a level of similarity between two ChIP-Seq profiles

RNA-seq data integration

Beta from Shirley Liu's lab in Harvard. Tao Liu's previous lab.

Heatmap, mata-plot

Many papers draw meta-plot and heatmap on certain genomic regions (2kb around TSS, genebody etc) using ChIP-seq data.

See an example from the ngs.plot:

Tools

  1. deeptools.It can do many others and have good documentation. It can also generate the heatmaps, but I personally use ngs.plot which is esy to use. (developed in Mount Sinai).

  2. you can also draw heatmaps using R. just count (using either Homer or bedtools) the ChIP-seq reads in each bin and draw with heatmap.2 function. here and here. Those are my pretty old blog posts, I now have a much better idea on how to make those graphs from scratch.

  3. You can also use bioconductor Genomation. It is very versatile.

  4. ChAsE

  5. Metaseq

  6. EnrichedHeatmaps from Zuguang Gu based on his own package ComplexHeatmaps. This is now my default go-to because of the flexiability of the package and the great user support. Thx!

  7. A biostar post discussing the tools: Visualizations of ChIP-Seq data using Heatmaps

  8. A bioconductor package to produce metagene plots

  9. Fluff is a Python package that contains several scripts to produce pretty, publication-quality figures for next-generation sequencing experiments I just found it 09/01/2016. looks promising especially for identifying the dynamic change.

One cavet is that the meta-plot (on the left) is an average view of ChIP-seq tag enrichment and may not reflect the real biological meaning for individual cases.

See a post from Lior Patcher How to average genome-wide data

I replied the post:

for ChIP-seq, in addition to the average plot, a heatmap that with each region in each row should make it more clear to compare (although not quantitatively). a box-plot (or a histogram) is better in this case . I am really uncomfortable averaging the signal, as a single value (mean) is not a good description of the distribution.

By Meromit Singer:

thanks for the paper ref! Indeed, an additional important issue with averaging is that one could be looking at the aggregation of several (possibly very distinct) clusters. Another thing we should all keep in mind if we choose to make such plots..

A paper from Genome Research Ubiquitous heterogeneity and asymmetry of the chromatin environment at regulatory elements

Enhancer databases

Interesting Enhancer papers

Enhancer target prediction

Allele-specific analysis

SNPs affect on TF binding

co-occurring TFs

  • In-silico Search for co-occuring transcription factors: INSECT
  • INSECT 2
  • CO-factors associated with Uniquely-bound GEnomic Regions:COUGER

Conservation of the peak underlying DNA sequences

Integration of different data sets

methylPipe and compEpiTools: a suite of R packages for the integrative analysis of epigenomics data

Copy number information from targeted sequencing using off-target reads bioconductor CopywriteR package.

3CPET: Finding Co-factor Complexes in Chia-PET experiment using a Hierarchical Dirichlet Process

ATAC-seq

Some may notice that the peaks produced look both like peaks produced from the TF ChIP-seq pipeline as well as the histone ChIP-seq pipeline. This is intentional, as ATAC-seq data looks both like TF data (narrow peaks of signal) as well as histone data (broader regions of openness).

DNase-seq

Chromatin Interaction data (ChIA-PET, Hi-C)