Using ChIP-Seq to identify and/or quantify bound regions (peaks)
Reviews
Evaluation of ChIP-Seq performance Wilbanks EG, and Facciotti MT (2010; PLoS ONE)
A manually curated ChIP-seq benchmark demonstrates room for improvement in current peak-finder programs (NAR Nov. 2010)
Step 1: Map reads
- See the mapping SOP for more details.
Step 2: Strand cross correlation analysis
- The goal of this step is to asses the quality of the IP and to estimate the fragment size of the DNA immunoprecipitated.
- For a detailed explanation on strand cross-correlation analysis see box 2 of this paper (Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data ).
/nfs/BaRC_Public/phantompeakqualtools/run_spp.R -c=TreatmentIP.bam -savp -out=TreatmentIP.run_spp.out
- After this analysis a good ChIP-seq experiment will have a second peak (reflecting the fragment size) at least as tall as the first peak (reflecting read length). This is how the graph should look: (Fig4E ). If the second peak is smaller that the first, (like the example shown in Fig4G Marginal ), macs will not estimate fragment size correctly. In that case we recommend running macs with parameters "--nomodel" and "--shiftsize=half_of_the_fragment_size". The fragment size is detected on the strand cross correlation analysis.
Step 3: Call peaks (bound regions)
Some of the parameters to consider when comparing programs are:
- Adjustment of sequence tags to better represent the original DNA fragment (by shifting tags in the 3′ direction or by extending tags to the estimated length of the original fragments)
- Background model used
- Use of strand dependent bimodality
BaRC performed a ChIP-Seq bake off in the summer of 2011, with these ChIP-Seq bake off results.
Based on our ChIP-Seq bake off and on a published review (Evaluation of ChIP-Seq performance), MACs and SISSRs are good programs to try.
MACS14
- For MACS to work the header of the sequences has to have no spaces.
- macs points to macs14 on WIBR local machines
- MACS may have trouble with a SAM file from bowtie if it contains unmapped reads (which it generally does). As a result, you may need to filter out unmapped reads with a command like
samtools view -hS -F 4 all_reads.sam > mapped_reads.sam
MACS (README) (Zhang et al., 2008)
- MACS is appropriate for both proteins like transcription factors that may have narrow peaks, as well as histone modifications that may affect broader regions. However, for broader peaks changing values of the parameters may be needed: eg. using --nomodel, --nolambda (if there's no control), --call-subpeaks. Running MACS with different parameters and viewing in IGV the results can help in choosing the optimal values.
Sample commands to run MACS (current version as of March 5 2012): 1.4.2 using mapped reads in map or sam format:
macs -t IP_mapped.map -c Control_mapped.map -g 1.87e9 --name=outputName --format=BOWTIE --tsize=36 --wig --space=25 --mfold=10,30 macs -t IP_mapped.sam -c Control_mapped.sam -g 1.87e9 --name=outputName --format=SAM --tsize=36 --wig --space=25 --mfold=10,30 macs -t IP_mapped.sam -c Control_mapped.sam -g 1.87e9 --name=outputName --format=SAM --tsize=36 --wig --space=25 --nomodel --shiftsize=100
The parameters used on the command are:
- -t TFILE Treatment file
- -c CFILE Control file
- --name=NAME Experiment name, which will be used to generate output file names. DEFAULT: "NA"
- --format=FORMAT Format of tag file, "BED" or "ELAND" or "ELANDMULTI" or "ELANDMULTIPET" or "SAM" or "BAM" or "BOWTIE". DEFAULT: "BED"
- --tsize=TSIZE Tag size. DEFAULT: auto detected tag size
- --wig Whether or not to save shifted raw tag count at every bp into a wiggle file
- --mfold=MFOLD Select the regions within MFOLD range of high-confidence enrichment ratio against background to build model. The regions must be lower than upper limit, and higher than the lower limit. DEFAULT:10,30
- -g GSIZE Effective genome size. It can be 1.0e+9 or 1000000000, or shortcuts:'hs' for human (2.7e9), 'mm' for mouse (1.87e9), 'ce' for C. elegans (9e7) and 'dm' for fruitfly (1.2e8), Default:hs
- --keep-dup=1 Controls the MACS behavior towards duplicate tags at the exact same location. DEFAULT: 1 in MACS 1.4; auto in MACS2.
- --nomodel whether or not to build the shifting model. If True, MACS will not build model. by default it means shifting size = 100.
- --shiftsize The arbitrary shift size in bp. When nomodel is true, MACS will use this value as 1/2 of fragment size. DEFAULT: 100.
bsub "macs14 -t IP_mapped.map -c Control_mapped.map --name=outputName --format=BOWTIE --tsize=36 --wig --space=25 --mfold=10,30" bsub "macs14 -t IP_mapped.sam -c Control_mapped.sam --name=outputName --format=SAM --tsize=36 --wig --space=25 --mfold=10,30"
Note: The wig files that macs14 generates are not normalized.
MACS2
- MACS2 is appropriate for both proteins like transcription factors that may have narrow peaks, as well as histone modifications that may affect broader regions. For broader peaks we recommend using --nomodel, --nolambda (if there's no control), and using the fragment size calculated on the strand cross correlation analysis. We recommend using macs2 rather than macs14 for broad peaks.
macs2 callpeak -t IP_reads.mapped_only.bam -c Control_reads.mapped_only.bam -f BAM -g mm -n Name --nomodel -B
- --nomodel whether or not to build the shifting model. If True, MACS will not build model. by default it means shifting size = 100.
- -f Input format
- -B create bedgraph output files
SISSRs
SiSSRs (Reference) (Manual) ChIP-Seq bake off
SISSRs uses strand bimodality to try to find the summit of the peak. The summit should be very close to the DNA bound by the TF. It is more appropriate for TFs because they tend to bind in specific narrow regions.
Map with Bowtie, use --sam parameter to get a SAM output file
bsub "bowtie -t -m 3 -n 3 -l 36 --strata --best --solexa1.3-quals --sam inputSeq bowtieOutput.sam"
SISRs input is a bed file. Convert mapped reads from SAM to BAM and from BAM to bed format
bsub "samtools view -S -b -o bowtieOutput.bam bowtieOutput.sam" -S input is SAM -b output BAM
bsub "bamToBed -i bowtieOutput.bam > bowtieOutput.bed"
Run SISSRs with a sample command like
sissrs.pl -i bowtieOutput.bed -o outputFile -s 2716965481 -b Background.bed -L 200
The parameters used in the sample command:
- -s is the size of the genome
- -L is the maximum length of the fragment
- -m is the percentage of mappable bps. Default is .8 for Eland in human.
For more detailed description of parameters see ChIP-Seq bake off
Step 4: Linking bound regions to genes
Both MACS and SISSRs provide bed files with the set of peaks, presumably indicating bound regions.
To link this regions to genes check out this SOP: Linking genome regions to genome annotations
Below it is a detailed example of how to do the following annotation:
- Find the genes that overlap with your peaks
- Find the genes that have a peak within -D distance upstream of the TSS/gene
- Find the genes that have a peak -D distance downstream of the gene
- Find the genes that have a peak -D distance downstream of the TSS
Note: We assume that the TSS is the 5' end of the transcript (gene), but you can also go to the UCSC Bioinformatics table browser and download the coordinates of the transcription start sites (TSS).
Comparing binding across different samples
- Method 1: Intersect bound regions (peaks)
- This binary comparison asks simply, "Is each peak in sample A present in sample B (and vice versa)?"
- This ignores the enrichment/score of each peak and gives each non-peak a score of 0.
- Method 2: Compare affinity (number of mapped reads) across bound regions
- Filter out redundant reads from each BAM file with a command like 'samtools rmdup'.
- Merge peaks of all samples (and make non-redundant):
- mergeBed -i Peaks_all_redundant.bed > Peaks.bed
- Count reads in each sample that overlaps each peak
- intersectBed -bed -wb -abam Reads_1.noDups.bam -b Peaks.bed | cut -f13-15 | mergeBed -i stdin -n > Peak_counts.reads_1.bed
- Modify counts per peak to account for differing numbers of mapped reads.
- Method 3: Use ngsplot to make stacked heatmaps of peaks for each sample
Step 1
Get the coordinates of the annotation features that you are going to use. For analysis A, B and C those coordinates can be downloaded from the UCSC tables.
i.e., on the table browser select: Mammal->Mouse group: Genes and Gene Prediction Tracks track: RefSeq Genes table: refGene region: genome output format: bed (name your file) get output
for A: Select Whole Gene for B: Select Upstream by D bases for C: Downstream by D bases
You will have now 3 annotation files named for example:
mm9_REFseqgenes.bed (gene bodies) mm9_REFseqgenes1kbUp.bed (1Kb upstream of the gene) mm9_REFseqgenes1kbDownofGene.bed (1Kb downstream of the gene)
The gene bodies file has more columns that you need, so keep only the first 6 columns.
cut -f1-6 mm9_REFseqgenes.bed > mm9_REFseqgenes_TRIM.bed
For annotation D, to get the coordinates of the regions from the TSS to a certain distance downstream of the TSS (beginning of the gene):
- Take the genes files and keep the coordinate that represents the start of the gene (based on the strand),
awk -F "\t" '{if ($6 == "+") {print $1"\t"$2"\t"$2+1"\t"$4"\t"$5"\t"$6} else {print $1"\t"$3-1"\t"$3"\t"$4"\t"$5"\t"$6} }' mm9_REFseqgenes.bed >| mm9_REFseqgenesTSS.bed
- Then use the tool "slopBed" to increase the range of genomic coordinates to include downstream of the TSS
The tool "slopBed" needs a file with the chromosome sizes that can be created with this command
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from mm9.chromInfo" > mm9.genome
Slop the file with the TSS coordinates to get the desired distance downstream (3') of the TSS
slopBed -i mm9_REFseqgenesTSS.bed -g mm9.genome -l 0 -r 1000 -s >| mm9_REFseqgenes_1kbdownofTSS.bed
Step 2: Intersect your peaks and your annotation files
intersectBed -a genes.bed -b peaks.bed -wa -wb -wa returns the list from file A -wb returns the list from file B
intersectBed -a mm9_1KbUpDownofTSS.bed -b peaks.bed -wa -wb >| 1KbUpDownofTSS_peaks.bed
