ChIP-Seq Bioinformatics Analysis

Quality Control

For any NGS analysis method, our first step in the workflow is to explore the quality of our reads prior to aligning them to the reference genome and proceeding with downstream analyses. A quality control will be used for this purpose for instance FastQC is a tool that provides a simple way to do quality control checks on raw sequence data coming from high throughput sequencing pipelines.

Mapping

The sequenced reads (in FASTQ format) are mapped using tools such as Bowtie2 or BWA etc. Bowtie2 and BWA can consider indels (insertions and deletions) by gapped alignments which is appropriate for long or paired-end reads. There are several output formats for map files, such as SAM, BAM, CRAM and tagAlign while the BAM format is the most widely used so far. After alignment, reads mapped to the same genomic positions are filtered as redundant reads, and the remaining non redundant reads are used for analysis.

Filter BAM

We filter the BAM file to keep only uniquely mapping reads. ChIP-seq datasets typically show mirror peaks of mapped reads on the ‘+’ and ‘−’ strand around binding sites of the target protein separated by a characteristic distance equal to the average lengths of the immunoprecipitated DNA fragments. Therefore, in order to represent the target of interest and generate signal tracks that show peak signal at the target binding site, it is important to account for this strand-specific read-shift. For example bamCoverage is a tool that takes an alignment of reads or fragments as input (BAM file) and generates a coverage track (bigWig or bedGraph) as output.

Peak Calling

Peak calling is a computational method used to identify areas in the genome that have been enriched with aligned reads as a consequence of performing a ChIP-sequencing experiment. Peak-calling results are generally returned in BED format. Although ChIP-seq peaks do not have strand information, it can be estimated from the gene information when focusing on the histone marks that are enriched around TSS (transcription start size). Peak Calling tools include BayesPeak, cisGenome, MACS2 etc. MACS2 is the most commonly used peak calling tool.

Quality Assessment of ChIP

Prior to performing any downstream analyses with the results from a peak caller, it is best practice to assess the quality of your ChIP-seq data. What we are looking for is good quality ChIP enrichment over background. To ensure the reproducibility of the experimental results, at least two biological replicates of each ChIP-seq experiment are recommended to be performed. The reproducibility of both reads and identified peaks should be examined. Any quality assessment tool can be used for instance ChIPQC (a bioconductor tool).

Peak Annotation

In order to identify sites that are differentially enriched between two or more sample groups we will be using a tool for instance; DiffBind that is an R Bioconductor package. The core functionality of DiffBind is the differential binding affinity analysis, which enables binding sites to be identified that are statistically significantly differentially bound between sample groups. The core analysis routines are executed, by default using DESeq2 with an option to also use edgeR. BED files will be created for each set of significant regions identified by DESeq2, separating them based on the gain or loss of enrichment. We will write these regions to file and use as input for downstream visualization.

In order to compare peak calls within groups, we use HOMER software that can detect the peaks that are constantly called in a single state or a single type replicate.

Visualization

In visualization first we create bigWig (standard file format commonly used for ChIP-seq data visualization) files and then we visualize enrichment patterns at particular locations in the genome. Individual software programs allow detailed analysis of peaks, biological interpretation, and visualization of ChIP-seq results. There are various strategies for visualizing enrichment patterns. Interactive visualization tools such as Integrated Genome Viewer (IGV) or SeqMonk are available. Several web servers (e.g UCSC genome browser and WashU Epigenome Browser) can integrate the obtained ChIP-seq results with other annotation data, such as evolutionary conservation and gene expression in various tissues.

Motif Discovery

Once you have your set of enriched regions which represent binding sites from your protein of interest, a next logical question is within those regions are there any similar patterns of sequences that are over-represented. These sequences are often referred to as motifs. Tools for motif analysis often require sequence information for each of your binding regions, which we will obtain using bedtools (a powerful toolset for genome arithmetic). Motif analysis tools include MEME-chip, HOMER motif analysis etc.

Annotation & Functional Enrichment

After identifying the regions of the genome that are enriched in the number of aligned reads for each of our transcription factors of interest, these enriched regions represent the likely locations of where these proteins bind to the genome. After we obtain a list of peak coordinates, it is important to study the biological implications of the protein–DNA bindings.

Many cis-regulatory elements are close to TSSs of their targets, it is common to associate each peak to its nearest gene, either upstream or downstream. ChIPseeker is an R Bioconductor package for annotating peaks. Additionally, it has various visualization functions to assess peak coverage over chromosomes and profiles of peaks binding to TSS regions.

Once we have obtained gene annotations for our peak calls, we can perform functional enrichment analysis to identify predominant biological themes among these genes by incorporating knowledge from biological ontologies such as Gene Ontology, KEGG and Reactome. The gene lists we have obtained through annotation can be interpreted using freely available web- and R-based tools, for instance ClusterProfiler(R-based tool)and GREAT(web-based tool).