Peak Calling and BED format
Peak calling is the method used to detect regions of the genome enriched in mapped reads. Those regions are called peaks.
In brief, the genome is scanned looking for regions where the number of mapped reads is higher than expected by random picking.
For avoiding false positives, a number of measures are taken such as the removal of PCR artifacts and the detection of “paired peaks” that are a diagnostic of good peaks.
There are a number of tools available for doing peak calling, one of the most popular being MACS.
For the detailed description about the MACS algorithm see the original paper and this course.
The ChIP-Seq reads is normally less than 100 bases but it covers only the ends of the ChIP fragments that is normally longer (between 200 and 400 bps dependings on the sonication size).
Since ChIP-DNA fragments are equally likely to be sequenced from both 5’ and 3’ ends, the read density around a true binding site should show two peaks, one at the 5’ and one at 3’. MACS takes advantage of this bimodal pattern to empirically model the fragment size d and shifts all reads d/2 toward the 3’ ends.
| MACS |
|---|
 |
| Zhang Y et al. Model-based analysis of ChIPSeq (MACS). Genome Biol. 2008;9(9):R137. |
MACS can be easily used either for the ChIP sample alone, or along with a control sample which increases specificity of the peak calls.
In the latter case, a sample (or, treatment) is immunoprecipitated, while a control, is not.
MACS then slides across the genome using a window size of 2d to find candidate peaks, estimates for each peak the probability to be found randomly in the local background, and calculates the adjusted p-values corrected for multiple comparison using the Benjamini-Hochberg correction.
Let’s look at one of the output files of MACS, the file H3K4me_peaks.xls. First let’s create the folder to store data on peaks.
pwd
cd ..
mkdir peaks
cd peaks
wget https://biocorecrg.github.io/PhD_course_genomics_format_2023/data/H3K4me_peaks.xls
more H3K4me_peaks.xls
# This file is generated by MACS version 2.2.4
# Command line: callpeak -t H3K4me1.bam -c input.bam -g 46709983 -n H3K4me
# ARGUMENTS LIST:
# name = H3K4me
# format = AUTO
# ChIP-seq file = ['H3K4me1.bam']
# control file = ['input.bam']
# effective genome size = 4.67e+07
# band width = 300
# model fold = [5, 50]
# qvalue cutoff = 5.00e-02
# The maximum gap between significant sites is assigned as the read length/tag size.
# The minimum length of peaks is assigned as the predicted fragment length "d".
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 1000 bps and 10000 bps
# Broad region calling is off
# Paired-End mode is off
# tag size is determined as 50 bps
# total tags in treatment: 514197
# tags after filtering in treatment: 201517
# maximum duplicate tags at the same position in treatment = 1
# Redundant rate in treatment: 0.61
# total tags in control: 612920
# tags after filtering in control: 421614
# maximum duplicate tags at the same position in control = 1
# Redundant rate in control: 0.31
# d = 276
# alternative fragment length(s) may be 276,300 bps
chr start end length abs_summit pileup -log10(pvalue) fold_enrichment -log10(qvalue) name
21 6565481 6565885 405 6565634 22.00 7.98323 3.66432 3.26622 H3K4me_peak_1
21 8426561 8427181 621 8426758 35.00 5.64064 2.27604 2.09468 H3K4me_peak_2
21 8433166 8433530 365 8433248 30.00 6.48973 2.65288 2.49481 H3K4me_peak_3
21 8438082 8438492 411 8438329 31.00 7.12466 2.77893 2.84380 H3K4me_peak_4
--More--(42%)
Another output file is H3K4me_peaks.narrowPeak, which, as you can see contains the information only about peaks. This is a BED file.
wget https://biocorecrg.github.io/PhD_course_genomics_format_2023/data/H3K4me_peaks.narrowPeak
more H3K4me_peaks.narrowPeak
21 6565480 6565885 H3K4me_peak_1 32 . 3.66432 7.98323 3.26622 153
21 8426560 8427181 H3K4me_peak_2 20 . 2.27604 5.64064 2.09468 197
21 8433165 8433530 H3K4me_peak_3 24 . 2.65288 6.48973 2.49481 82
21 8438081 8438492 H3K4me_peak_4 28 . 2.77893 7.12466 2.84380 247
BED format
This tab-separated text format is called Browser Extensible Data (BED) and is normally contains 3, 6 or 12 columns.
Here is the description of each column:
| Column number | Column name | Details |
|---|---|---|
| 1 | chrom | name of the chromosome or scaffold |
| 2 | chromStart | The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0 |
| 3 | chromEnd | The ending position of the feature in the chromosome or scaffold |
| 4 | name | Defines the name of the BED line |
| 5 | score | A score between 0 and 1000 |
| 6 | strand | Defines the strand. Either “.” (=no strand) or “+” or “-“ |
| 7 | thickStart | The starting position at which the feature is drawn thickly |
| 8 | thickEnd | The ending position at which the feature is drawn thickly |
| 9 | itemRgb | An RGB value of the form R,G,B (a color code) |
| 10 | blockCount | The number of blocks (exons) in the BED line |
| 11 | blockSizes | A comma-separated list of the block sizes |
| 12 | blockStarts | A comma-separated list of block starts |
In the case of our narrowPeak file, the first six fields correspond indeed to the first six fields of a BED file, while the remaining ones are MACS-specific columns:
| Column number | Column name | Details |
|---|---|---|
| 7 | fold_enrichment | fold enrichment for the region |
| 8 | -LOG10(pvalue) | pvalue (-log10) |
| 9 | -LOG10(qvalue) | qvalue (i.e. false discovery rate) (-log10) |
| 10 | dist summit | The distance between the summit of the peak and the start |
More on MACS output files can be found here.
Let’s count how many peaks were found:
wc -l H3K4me_peaks.narrowPeak
We can check how many peaks have q-value (in column #9) larger than 0.01 (i.e., -log10(0.01) = 2).
awk -F"\t" '{if ($9>=2) print}' H3K4me_peaks.narrowPeak | wc -l
28
EXERCISE
- How many peaks in the file H3K4me_peaks.narrowPeak have fold enrichment larger than 5 ?