Last Updated: Jul 7, 2013 MAILING LISTPlease join the IDR mailing list https://groups.google.com/group/idr-discuss for FAQs, discussions and updates on software. SUMMARYReproducibility is essential to reliable scientific discovery in high-throughput experiments. The IDR (Irreproducible Discovery Rate) framework is a unified approach to measure the reproducibility of findings identified from replicate experiments and provide highly stable thresholds based on reproducibility. Unlike the usual scalar measures of reproducibility, the IDR approach creates a curve, which quantitatively assesses when the findings are no longer consistent across replicates. In layman's terms, the IDR method compares a pair of ranked lists of identifications (such as ChIP-seq peaks). These ranked lists should not be pre-thresholded i.e. they should provide identifications across the entire spectrum of high confidence/enrichment (signal) and low confidence/enrichment (noise). The IDR method then fits the bivariate rank distributions over the replicates in order to separate signal from noise based on a defined confidence of rank consistency and reproducibility of identifications i.e the IDR threshold.The method was developed by Qunhua Li and Peter Bickel's group and is extensively used by the ENCODE and modENCODE projects and is part of their ChIP-seq guidelines and standards. A publication describing the statistical details of the IDR framework is http://www.stat.washington.edu/qli/IDR-AOAS.pdf or http://projecteuclid.org/euclid.aoas/1318514284 Below we provide an automated pipeline around the IDR framework with various tweaks to systematically analyze ChIP-seq data. One important addition is the use of the IDR framework for assessing the self-consistency of a single dataset and for "rescuing" datasets (pooled-consistency) when replicates are not very comparable in terms of data quality (i.e. getting the maximum juice out of your datasets). This happens more often that you might imagine. It is not trivial to obtain "truly equivalent and comparable replicates". In the event, that one replicate is significantly better than the other, any evaluation or thresholding based on reproducibility will be lower bounded by the worst replicate. Ofcourse, the ideal scenario is to repeat your experiments until you get high quality comparable replicates. However, due to limitations of resources (financial or experimental material) this is often not accomplished. Hence from a practical stand point, the rescue strategy allows one to be less conservative and "smooth out" artificial differences between replicates. A publication on this pipeline is in the works ... NOTE: The current pipeline listed below is heavily tested and verified for human transcription factor ChIP-seq data and punctate chromatin marks (e.g. H3k9ac, H3k27ac, H3k4me3, H3k4me1 etc.) using the SPP peak caller. We do NOT recommend using it as is for broad chromatin marks ChIP-seq. We are still developing a robust pipeline for broad chromatin marks (e.g. H3k9me3, H3k27me3, H3k36me3 etc.). The main issue is that most peak callers do not work well with broad marks and so we are switching to a generic genome-wide segmentation or windowing approach for broad marks. We are also testing a few callers such as RSEG and MACS2 which claim to be decent for broad marks. INTUITIVE EXPLANATION OF IDR AND IDR PLOTSFor an intuitive description of the IDR method and a description of the various plots generated by the IDR package read this IDR-101 document compiled by Qunhua Li. A detailed explanation of the model and the underlying statistics is in the AoAS paper CODE FOR IDR ANALYSISIDR CODE: Download this Code archiveNOTE: This is different from the CRAN version of the IDR package. The CRAN package implements a slightly different underlying algorithm and is more general purpose for pairs of ranked lists. It has not yet been fully tested in the ChIP-seq setting and does not have the necessary wrappers to work directly with peak files. So I recommend using the above linked code even though it is a bit clunky. We will be releasing a unified version soon. SPP PEAK CALLER CODE: http://code.google.com/p/phantompeakqualtools/ Modified from the official version to add some additional functionality and an easy to use unified script. I recommend using this version.
=========================== Qunhua Li and Anshul Kundaje (Oct,2010) =========================== This set of programs are used for consistency analysis on peak calling results on multiple replicates of a dataset ================ DEPENDENCIES ================ unix, R version 2.9 or higher ================ FILES: ================
================ INPUT FILE FORMATS ================ (1) genome_table.txt It contains two space delimited fields Col1: chromosome name (These MUST match the chromosome names in the peak files) Col2: chromosome size (in bp) (1) Peak Files Peak files MUST be in narrowPeak format (and unzipped ... the code currently doesnt handle gzipped peak files directly) NarrowPeak files are in BED6+4 format. It consists of 10 tab-delimited columns 1.chrom string Name of the chromosome 2.chromStart int The starting position of the feature in the chromosome. The first base in a chromosome is numbered 0. 3.chromEnd int The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as chromStart=0, chromEnd=100, and span the bases numbered 0-99. 4.name string Name given to a region (preferably unique). Use '.' if no name is assigned 5.score int Indicates how dark the peak will be displayed in the browser (1-1000). If '0', the DCC will assign this based on signal value. Ideally average signalValue per base spread between 100-1000. 6.strand char +/- to denote strand or orientation (whenever applicable). Use '.' if no orientation is assigned. 7.signalValue float Measurement of overall (usually, average) enrichment for the region. 8.pValue float Measurement of statistical signficance (-log10). Use -1 if no pValue is assigned. 9.qValue float Measurement of statistical significance using false discovery rate (-log10). Use -1 if no qValue is assigned. 10.peak int Point-source called for this peak; 0-based offset from chromStart. Use -1 if no point-source called.
*NOTE*: the p-value and q-value columns MUST be in -log10() scale
The narrowPeak format has 3 columns that can be used to rank peaks
(signal.value, p.value (-log_10) and q.value (-log_10)).
The peak summit column must have values relative to the start coordinate of the peaks.
You can use any of these columns but make sure that whichever measure you are using rank peaks is relatively continuous without too many ties.
e.g. For the SPP peak caller it is recommended to use signal.value column
e.g. For MACS2 it is recommended to use the p.value column
e.g. PeakRanger/PeakSeq peak caller has relatively continuous q.values without too many ties. So for PeakSeq it is better to use q.value================ RUNNING INSTRUCTIONS ================ (1) batch-consistency-analysis.r This is used to run pairwise consistency analysis on a pair of replicate peak files ---------------- GENERAL USAGE: ---------------- [outfile.prefix] is a prefix that will be used to name the output data for this pair of replicates. The prefix must also include the PATH to the directory where you want to store the output data. e.g. /consistency/reps/chipSampleRep1_VS_chipSampleRep2 [min.overlap.ratio]: fractional bp overlap (ranges from 0 to 1) between peaks in replicates to be considered as overlapping peaks. Set to 0 if you want to allow overlap to be defined as >= 1 bp overlap. If set to say 0.5 this would mean that atleast 50% of the peak in one replicate should be covered by a peak in the other replicate to count as an overlap. IMPORTANT: This parameter has not been tested fully. It is recommended to set this to 0. [is.broadpeak]: Is the peak file format narrowPeak or broadPeak. Set to F if it is narrowPeak/regionPeak or T if it is broadPeak. BroadPeak files do not contain Column 10.
[ranking.measure] is the ranking measure to use. It can take only one of the following values
signal.value , p.value or q.value
OUTPUT:
The results will be written to the directory contained in [outfile.prefix]
a. The output from EM fitting: suffixed by -em.sav
b. The output for plotting empirical curves: suffixed by -uri.sav
Note: 1 and 2 are objects that can be loaded back to R for plotting or other purposes (e.g. retrieve data)
c. The parameters estimated from EM and the log of consistency analysis, suffixed by -Rout.txt
d. The number of peaks that pass specific IDR thresholds for the pairwise analysis: suffixed by npeaks-aboveIDR.txt
e. The full set of peaks that overlap between the replicates with local and global IDR: suffixed by overlapped-peaks.txt
(2) batch-consistency-plot.r
This is used to plot the IDR plots and diagnostic plots for a single or multiple pairs of replicates.
----------------
GENERAL USAGE:
----------------
Rscript batch-consistency-plot.r [npairs] [output.prefix] [input.file.prefix1] [input.file.prefix2] [input.file.prefix3] ....
[n.pairs] is the number of pairs of replicates that you want to plot on the same plot
e.g. 1 or 3 or ...
[output.prefix] is a prefix that will be used to name output data from this analysis.
NOT TO BE CONFUSED with [outfile.prefix] in batch-consistency-analysis.r
The prefix must also include the PATH to the directory where you want to store the output data.
e.g. /consistency/plots/chipSampleAllReps
[input.file.prefix 1, 2, 3 ...] are the [outfile.prefix] values used to name the output from pairwise analysis on all replicates
e.g. /consistency/reps/chipSampleRep1_VS_chipSampleRep2
/consistency/reps/chipSampleRep1_VS_chipSampleRep3
/consistency/reps/chipSampleRep2_VS_chipSampleRep3
OUTPUT:
1. summary consistency plots in .ps format: suffixed by -plot.ps
These plots are very informative about the quality and similarity of the replicates.
===================================================
GETTING NUMBER OF PEAKS THAT PASS AN IDR THRESHOLD
===================================================
For each pairwise analysis, we have a *overlapped-peaks.txt file
The last column (Column 11) of the overlapped-peaks.txt file has the global IDR score for each pair of overlapping peaks
To get the number of peaks that pass an IDR threshold of T (e.g. 0.01) you simply find the number of lines that have a global IDR score <= T
IDR PIPELINEMost of these instructions are for use with the SPP peak caller. I have created a slightly modified version of the official SPP package (with some additional functionality and a pretty easy to use automated script for peak calling) that you can download here. I have added code snippets for the MACSv2 peak caller as well where necessary in the README below. I highly recommend NOT using MACS1.4 with IDR due to 3 fundamental problems with MACS1.4.
We recommend using the SPP peak caller with IDR since it has been thoroughly tested. If you do plan to use MACS, use MACSv2 which has fixed each of these 3 problems listed above. We would highly recommend NOT using MACS1.4 with IDR. Ok so now lets start with the analysis pipeline: CALL PEAKS ON INDIVIDUAL REPLICATES- Lets say we have a particular TF ChIP-seq dataset. I will refer to it as chipSample from here on. - Lets say we have 3 replicates (read alignment files) for chipSample named - chipSampleRep1.tagAlign.gz - chipSampleRep2.tagAlign.gz - chipSampleRep3.tagAlign.gz NOTE: The tagAlign format is a simple BED format for If your alignment files are in BAM format, convert them to tagAlign using SAMtools and the bamToBed function from bedTools - Each of the ChIP replicates might be paired with its own input/control tagAlign file or a single tagAlign file. - controlSampleRep1.tagAlign.gz - controlSampleRep2.tagAlign.gz - controlSampleRep3.tagAlign.gz - The first thing we need to do is pool all control datasets that correspond to all replicates of chipSample i.e. simply concatenate the control tagAlign files. NOTE: If you have removed duplicates from your sample use run_spp_nodups.R instead of run_spp.R otherwise you will get errors - This will generate 3 peak files in directory /peaks/reps/ chipSampleRep1_VS_controlSampleRep0.regionPeak.gz chipSampleRep2_VS_controlSampleRep0.regionPeak.gz chipSampleRep3_VS_controlSampleRep0.regionPeak.gz These file are in the standard UCSC/ENCODE narrowPeak format which is an extended BED format. - Strand cross-correlation based data quality measures for the 3 runs will be appended to file /stats/phantomPeakStats.tab (See README.txt in the SPP peak calling package for details about the strand-cross-correlation based data quality measures files) IMPORTANT: It is important to check the cross-correlation plot that is produced by SPP to make sure the estimated fragment length is appropriate. Often for weaker datasets or datasets with few binding sites, the estimated fragment length may be incorrectly estimated as read-length or a value close to read-length (due to a phantom peak in the cross-correlation profile). The code automatically tries to account for this. However, it can make mistakes. For most ChIP-seq datasets, it is good to use the exclusion parameter (-x) with run_spp.R which will help the code skip any cross-correlation peak in the exclusion zone (values that are impossible as fragment lengths). I recommend that you use -x=-500:85 it is relatively unlikely that your fragment lengths will be in that range. #MACS2: If you are using MACS2 then use the following commands macs2 callpeak -t
IMPORTANT: It is important to check the cross-correlation plot that is produced by SPP to make sure the estimated fragment length is appropriate. Often for weaker datasets or datasets with few binding sites, the estimated fragment length may be incorrectly estimated as read-length or a value close to read-length (due to a phantom peak in the cross-correlation profile). The code automatically tries to account for this. However, it can make mistakes. For most ChIP-seq datasets, it is good to use the exclusion parameter (-x) with run_spp.R which will help the code skip any cross-correlation peak in the exclusion zone (values that are impossible as fragment lengths). I recommend that you use -x=-500:85 it is relatively unlikely that your fragment lengths will be in that range. - This will generate a peak file in directory /peaks/pooled chipSampleRep0_VS_controlSampleRep0.regionPeak.gz #MACS2: Use analogous commands are for the individual original replicates (above). IMPORTANT: It is important to check the cross-correlation plot that is produced by SPP to make sure the estimated fragment length is appropriate. Often for weaker datasets or datasets with few binding sites, the estimated fragment length may be incorrectly estimated as read-length or a value close to read-length (due to a phantom peak in the cross-correlation profile). The code automatically tries to account for this. However, it can make mistakes. For most ChIP-seq datasets, it is good to use the exclusion parameter (-x) with run_spp.R which will help the code skip any cross-correlation peak in the exclusion zone (values that are impossible as fragment lengths). I recommend that you use -x=-500:85 it is relatively unlikely that your fragment lengths will be in that range.- This will create two peak files in the directory /peaks/selfPseudoReps chipSampleRep1.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep1.pr2_VS_controlSampleRep0.regionPeak.gz #MACS2: Use analogous commands are for the individual original replicates (above). IMPORTANT: It is important to check the cross-correlation plot that is produced by SPP to make sure the estimated fragment length is appropriate. Often for weaker datasets or datasets with few binding sites, the estimated fragment length may be incorrectly estimated as read-length or a value close to read-length (due to a phantom peak in the cross-correlation profile). The code automatically tries to account for this. However, it can make mistakes. For most ChIP-seq datasets, it is good to use the exclusion parameter (-x) with run_spp.R which will help the code skip any cross-correlation peak in the exclusion zone (values that are impossible as fragment lengths). I recommend that you use -x=-500:85 it is relatively unlikely that your fragment lengths will be in that range.- This will create two peak files in the directory /peaks/pooledPseudoReps chipSampleRep0.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep0.pr2_VS_controlSampleRep0.regionPeak.gz #MACS2: Use analogous commands are for the individual original replicates (above). For IDR analysis, we always compare a pair of peak files For dataset chipSample, we will have ended up with the following files PEAKS ON ORIGINAL REPLICATES chipSampleRep1_VS_controlSampleRep0.regionPeak.gz chipSampleRep2_VS_controlSampleRep0.regionPeak.gz chipSampleRep3_VS_controlSampleRep0.regionPeak.gz PEAKS ON SELF-PSEUDOREPLICATES chipSampleRep1.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep1.pr2_VS_controlSampleRep0.regionPeak.gz chipSampleRep2.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep2.pr2_VS_controlSampleRep0.regionPeak.gz chipSampleRep3.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep3.pr2_VS_controlSampleRep0.regionPeak.gz PEAKS ON POOLED PSEUDOREPLICATES chipSampleRep0.pr1_VS_controlSampleRep0.regionPeak.gz chipSampleRep0.pr2_VS_controlSampleRep0.regionPeak.gz #MACS2: If you are using peaks based on the MACS2 peak caller, then use p.value as the ranking measure.
#MACS2: If you are using peaks based on the MACS2 peak caller, then use p.value as the ranking measure.You can plot the IDR plots using
The following IDR thresholds are recommended for the different types of IDR analyses. For self-consistency and comparison of true replicates
For pooled-consistency analysis
Each of these comparisons will give us a certain number of peaks that pass their respective IDR thresholds. We will refer to them as followsORIGINAL REPLICATE THRESHOLD
#MACS2: Sort peaks by p.value before selecting the top 'N' peaks #MACS2: Sort peaks by p.value before selecting the top 'N' peaks
|
|
Projects >
ENCODE: TF ChIP-seq peak calling using the Irreproducibility Discovery Rate (IDR) framework |
