\documentclass[10pt]{article}
\usepackage[hmargin=1.5cm,top=2cm,bottom=2cm]{geometry}
\usepackage{multicol}
\setlength\columnsep{15pt}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{array}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage[auth-sc]{authblk}
\usepackage{longtable,array,ragged2e}
\usepackage{multirow}
\usepackage{hyperref}
\usepackage{enumerate}
\usepackage[labelfont=bf]{caption}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage[framemethod=TikZ]{mdframed}
\usepackage{graphics}
\usepackage{yfonts,color}
\usepackage{tikz}
\usepackage{pifont}
\usepackage{rotating}
\usepackage{array}
\usepackage{lscape}
\usepackage{caption}
\usepackage{cprotect}
\usepackage{breakurl}
\usepackage{todonotes}
\usepackage{hanging}
\usepackage[final]{pdfpages}
\usepackage[leftFloats,CaptionAfterwards]{fltpage}
\usepackage{abstract}
\usepackage{enumitem}
\usepackage{soul}
\usepackage{titlesec}
\usepackage[numbers,super,sort&compress]{natbib}
\titleformat{\section}[block]{\large\bfseries\filcenter}{\thesection.}{0.4em}{}
\titleformat{\subsection}[block]{\normalsize\sc\bfseries\filcenter}{\thesubsection.}{0.4em}{}
\setcounter{secnumdepth}{5}
% \usepackage{fancyhdr}
% \pagestyle{fancy}
% \renewcommand{\sectionmark}[1]{
% \markright{\thesection\ #1}}
% \fancyhf{} % delete current header and footer
% \fancyhead[CO]{Third-generation sequencing and functional genomics}
% \fancyhead[CE]{Georgi K. Marinov}
% \fancyhead[RO,LE]{\thepage}
% \fancyhead[LO]{\bfseries\rightmark}
% \fancyhead[RE]{\bfseries\leftmark}
% \renewcommand{\headrulewidth}{0pt}
% \renewcommand{\footrulewidth}{0pt}
% \addtolength{\headheight}{0.5pt} % space for the rule
% \fancypagestyle{plain}{%
% }

\newcolumntype{P}[1]{>{\RaggedRight\arraybackslash\hspace{0pt}}p{#1}}

\setlength{\bibsep}{0pt plus 0.3ex}

\makeatletter
\def\@biblabel#1{\@ifnotempty{#1}{#1.}}
\makeatother

\newcommand{\filllastline}[1]{
\setlength\leftskip{0pt}
\setlength\rightskip{0pt}
\setlength\parfillskip{0pt}
#1}

\newenvironment{Figure}
{\par\medskip\noindent\minipage{\linewidth}}
{\endminipage\par\medskip}

% \title{\bf Genome-wide screening reveals modest importance of individual loop anchor CTCF sites for cellular fitness and gene expression in differentiated mammalian cells}
% \title{\bf Massively parallel characterization of perturbations of TAD loop anchor CTCF sites}
\title{\bf Automated quality control and reproducible peak calling for transcription factor ChIP-seq data}
\renewcommand\Authfont{\scshape\normalsize}
\author[1,2,*,$\#$]{Anshul Kundaje}
\author[3,*]{Youngsook L. Jung}
\author[2,*]{Jin Wook Lee}
\author[4,*]{Qunhua Li}
\author[2,*]{Nathan Boley}
\author[2,*]{Georgi K. Marinov}
\author[2]{J. Seth Strattan}
\author[3,5]{Peter V. Kharchenko}
\author[11]{Yuchun Guo}
\author[12]{Arif Harmanci}
\author[13]{Dongjun Chung}
\author[19]{Tao Liu}
\author[21]{Robert E. Thurman}
\author[20]{Noam Shoresh}
\author[14]{Russell Bonneville}
\author[15]{Flo Pauli}
\author[2]{Trupti Kawli}
\author[6]{Joel Rozowsky}
\author[7]{James B. Brown}
\author[2]{Cricket A. Sloan}
\author[2]{Benjamin C. Hitz}
\author[2]{Esther T. Chan}
\author[2]{Jason A. Hilton}
\author[2,10]{Arend Sidow}
\author[1]{Serafim Batzoglou}
\author[2]{Michael P. Snyder}
\author[15]{Richard Myers}
\author[16]{Peggy J. Farnham}
\author[17]{Victor Jin}
\author[12]{Mark Gerstein}
\author[11]{David Gifford}
\author[18]{Sunduz Keles}
\author[8]{Barbara J. Wold}
\author[2]{J. Michael Cherry}
\author[10]{Peter J. Bickel}
\author[3,$\#$]{Peter J. Park}
\renewcommand\Affilfont{\itshape\small}
\affil[1]{Department of Computer Science, Stanford University, Stanford, California, 94305, USA}
\affil[2]{Department of Genetics, Stanford University, Stanford, California, 94305, USA}
\affil[3]{Department of Biomedical Informatics, Harvard Medical School, Boston, Massachusetts, 02115, USA}
\affil[4]{Department of Statistics, Penn State University, University Park, Pennsylvania, 16802, USA}
\affil[5]{Division of Hematology Oncology, Boston Children’s Hospital, Massachusetts, 02115, USA}
\affil[6]{Department of Molecular Biophysics and Biochemistry, Yale University, Yale, Connecticut, 06520, USA}
\affil[7]{Lawrence Berkeley Laboratory, Berkeley, California, 94710, USA}
\affil[8]{Division of Biology and Beckman Institute, California Institute of Technology, California, 91125, USA}
\affil[9]{Department of Pathology, Stanford University, Stanford, California, 94305, USA}
\affil[10]{Department of Statistics, University of California at Berkeley, Berkeley, California, 94710, USA}
\affil[11]{Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA}
\affil[12]{Program in Computational Biology and Bioinformatics, Yale University, New Haven, Connecticut, 06520, USA}
\affil[13]{Division of Biostatistics and Bioinformatics, Medical University of South Carolina, Charleston, South Carolina, 29425, USA}
\affil[14]{Comprehensive Cancer Center, The Ohio State University College of Medicine, Columbus, Ohio, 43210, USA}
\affil[15]{HudsonAlpha Institute for Biotechnology, Huntsville, Alabama, 35806, USA}
\affil[16]{Department of Biochemistry and Molecular Medicine, University of Southern California, Los Angeles, California, 90089, USA}
\affil[17]{Department of Molecular Medicine, University of Texas Health Science Center, San Antonio, Texas, 78229, USA}
\affil[18]{Department of Biostatistics and Medical Informatics, University of Wisconsin, Madison, Wisconsin, 53705, USA}
\affil[19]{Department of Biochemistry, The State University of New York, Buffalo, New York, 14203, USA}
\affil[20]{Broad Institute, Cambridge, Massachusetts, 02142, USA}
\affil[21]{Department of Genome Sciences, University of Washington, Seattle, Washington, 98195, USA}
\affil[*]{These authors contributed equally to this work}
\affil[$\#$]{Correspondence should be addressed to A.K. (\href{mailto:anshul@kundaje.net}{anshul@kundaje.net}) and P.J.P. (\href{mailto:peter\_park@harvard.edu}{peter\_park@harvard.edu}).}
\date{}

\begin{document}
\maketitle

% \centerline{}
% \centerline{}
\begin{abstract}

\noindent {\normalsize \textbf{Chromatin immunoprecipitation with high-throughput sequencing (ChIP-seq) is the predominant assay for generating genome-wide maps of protein-DNA interactions. We present AQUAS (\underline{A}utomated \underline{QU}ality control and reproducibility \underline{A}nalysis of ChIP-\underline{S}eq data), an end-to-end protocol for robust data processing and quality control of transcription factor ChIP-seq experiments. AQUAS supports several complementary data quality indicators including the strand cross-correlation enrichment measure and uses the Irreproducibility Discovery Rate (IDR) framework for identifying reproducible binding events from replicate experiments. We provide standalone and cloud-based software implementations that support resumable distributed workflows, user-friendly analysis reports and visualization via genome-browser sessions. The protocol has been used by the ENCODE project to analyze over 6,500 ChIP-seq datasets. AQUAS enables investigators to critically assess and process their TF ChIP-seq data in concordance with ENCODE standards. The workflow takes $\sim$1 day to process a typical pair of single-end or paired-end TF ChIP-seq replicate experiments. 
}}
\centerline{}
\centerline{}
\end{abstract}

\begin{multicols}{2}

\section*{Introduction}

Chromatin immunoprecipitation with high-throughput sequencing (ChIP-seq) is an assay for profiling protein-DNA interactions and chromatin modifications genome-wide. It was one of the first applications of next-generation sequencing technologies for functional genomics and remains one of the most widely used assays for regulatory and epigenomic profiling of genomes\cite{Ren2000,Barsk2007}. However, ChIP-seq experiments are sensitive to a myriad experimental parameters as well as the properties of the target protein/marker. Hence, ChIP-seq experiments can exhibit significant variation in data quality and default analysis workflows can often result in unstable or poorly calibrated results \cite{Marinov2014}. As the number of ChIP-seq datasets continues to grow rapidly in individual laboratories\cite{Dreos2017} and large coordinated consortia \cite{ENCODE2012,modENCODE2010,Gerstein2010,Roadmap2015}, there is a need for robust protocols for automated large-scale data processing and objective, calibrated measures of data quality that enable systematic critical evaluation of datasets. Here, we present AQUAS (\underline{A}utomated \underline{QU}ality control and reproducibility \underline{A}nalysis of ChIP-\underline{S}eq data), an end-to-end protocol with stand-alone and cloud-based software implementations for large-scale, automated data processing and quality control of ChIP-seq datasets targeting proteins with punctate binding landscapes.

\begin{figure*}
\begin{center}
\includegraphics[width=16cm]{Fig1.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\cprotect\caption{{\bf The workflow of the AQUAS pipeline}. 
Cross-correlation measures and IDR analysis in AQUAS are integrated with other quality control measures. \verb|Rep0| is the merged file from replicates 1 (\verb|Rep1|) and 2 (\verb|Rep2|); \verb|Ctl0| is the merged file from control replicates 1 (\verb|Ctl1|) and 2 (\verb|Ctl2|); \verb|pr1| or \verb|pr2| stands for pseudo-replicate 1 or 2. $N_t$ and $N_p$ are the replicate and the pooled-pseudoreplicate thresholds, respectively. $N_1$ and $N_2$ are self-consistency thresholds for Rep1 and Rep2, respectively (see Experimental design and Procedure for details).}
\label{Fig1}
\end{figure*}

A typical ChIP-seq experiment consists of a series of steps including antibody selection and validation, sample preparation, co-immunoprecipitation, DNA library preparation, and sequencing\cite{Park2009,Pepke2009}. Critical factors affecting data quality are antibody specificity and efficiency, quality of sonication, size selection protocols and the resulting distribution of immunoprecipitated DNA fragment sizes, library complexity, presence of PCR artifacts, sequencing read quality and sequencing depth\cite{Landt2012,Jung2014}. Downstream computational analysis involves mapping of sequenced reads to the genome, creating genome-wide signal coverage maps and identifying regions of enrichment (i.e., peaks) often in comparison to a corresponding control experiment (i.e., input DNA). An ideal ChIP-seq experiment should have high specificity and sensitivity, i.e., significant clustering of immunoprecipitated (IP) fragments centered at the genuine binding (or target) sites of the protein (or chromatin mark) of interest with a low background signal throughout the rest of the genome. A truly comprehensive ChIP-seq experiment must also be sufficiently ``saturated'' as to allow discovery of the a majority of true target sites. Replicate experiments should also be performed to ensure the identification of reproducible target sites. Biological replicates are strongly preferred, since most of the variation comes from variability in biological conditions and lack of specificity for the antibody.

\hl{Users have frequently} resorted to visualizing the processed data in a genome browser, such as the UCSC Genome Browser\cite{Tyner2017}, WashU Epigenome Browser\cite{Zhou2011} and IGV\cite{Thorvaldsdottir2013}, to make a subjective judgment on data quality and reproducibility between replicates. While this practice is often useful, a systematic and automated analysis is necessary to ensure that the assessment is unbiased, especially when less is known about the protein or modification of interest. Objective measures also allow researchers to compare and evaluate \hl{samples generated using different reagents and by different laboratories}. The need for such analyses is more acute in projects that generate a large number of ChIP-seq experiments, since it is not feasible to visually inspect each profile. 

A large number of tools commonly referred to as ``peak callers'' have been developed for identifying in vivo binding events of transcription factors (which typically have punctate, sharp peak of enrichment) and other chromatin-associated proteins as well as genomic landscapes of histone variants and modifications (which typically exhibit broad regions of enrichment)\cite{Zhang2008,Kharchenko2008}. However, the choice of peak caller and its parameters can influence analysis results significantly. Default analysis parameters can produce vastly differing results from different peak callers applied to the same data % (See Section \hl{XX}). 
Moreover, default parameters of peak callers often cannot generalize to the diversity of binding landscapes of DNA binding proteins in terms of number of binding sites and signal-to-noise ratios.  Users are often left with no choice but to manually tune parameters to satisfy subjective definitions of optimality -- an approach that is not scalable. Hence, there is a need for\hl{ }automated methods for robust and reliable peak calling that can adapt to the properties to the data.

The protocol we describe here has been developed as part of the Encyclopedia of DNA Elements (ENCODE) project. ENCODE has been generating a massive reference compendium of ChIP-seq data for a wide range of target proteins and histone modifications across multiple cellular contexts and organisms\cite{ENCODE2012,modENCODE2010,Gerstein2012,Gerstein2010}, necessitating the development of a robust and automated quality control and statistical analysis pipeline. Many of the guidelines and practices of the consortia with respect to ChIP-seq data have been summarized previously\cite{Landt2012}. Here we describe the quality control and reproducibility steps in greater detail as well as providing \hl{an overview of the results of applying these steps to the Consortia data to create a reference set of values}. 

The overview of \hl{the AQUAS protocol is shown in Figure}  \ref{Fig1}. It has two key components:

(1)	 Quality measures based on cross-correlation analysis of strand-specific read density profiles (hereafter referred to as ``cross-correlation measures'')

(2)	Irreproducible discovery rate (IDR) framework for measuring reproducibility across experiments\cite{Li2011}

\filllastline{The cross-correlation \hl{metrics} capture the degree of\hl{ }fragment clustering \hl{in a dataset} by accounting for the asymmetric, stranded structure of read distributions expected near true binding sites. These \hl{metrics} can be efficiently computed after read alignment but before peak calling, thereby providing early feedback on data quality independent of peak calling algorithms and parameters. The IDR framework provides an approach to measure the} 

\end{multicols}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 1 $|$ Irreproducible Discovery Rate (IDR)}}
\centerline{}
\centerline{}
IDR is a statistical method to assess the reproducibility of entries on rank lists\cite{Li2011}. Given a set of peak calls, the peaks can be ranked based on a criterion of significance, such as the $p$-value, the $q$-value, or the fold enrichment. Significant peaks generally are ranked more consistently across the replicates than the peaks with low significance. This change of consistency provides an indicator of the transition from real signal to noise. IDR quantifies this transition by modeling the peaks as a mixture of a reproducible and an irreproducible group using a copula mixture model, where the peaks in the reproducible group are ranked higher and more consistently across replicates than the ones in the irreproducible group. It assigns each peak a reproducibility index, called the local IDR, to estimate its probability to be reproducible, and also reports the expected rate of irreproducible discoveries in the selected peaks (called IDR) in a fashion analogous to that of false discovery rate (FDR).  Similar to FDR, peaks with lower IDR are more reliable. Hence, an IDR threshold can be used as a reproducibility-based criterion to determine the rank threshold of the peaks at which signal transits to noise. The IDR thresholds exhibit significantly more stable behavior than FDR thresholds provided by peak callers. A single IDR threshold can be used comparably across datasets as it reflects the inherent sampling reproducibility of a dataset and is not significantly affected by antibodies, labs, platforms and analysis protocols\cite{Landt2012}.  When using IDR, a relatively relaxed peak-calling threshold is advised because the IDR algorithm requires sampling of both signal and noise distributions to assess the reproducibility of peaks.\par}
\end{mdframed}
\centerline{}

\begin{multicols}{2}

\noindent reproducibility of the peaks from replicate experiments and generates highly stable thresholds based on reproducibility\cite{Li2011} (more details below). This protocol has been applied to over 2000 diverse TF ChIP-seq datasets from the ENCODE and modENCODE\cite{ENCODE2012,modENCODE2010,Gerstein2012,Gerstein2010} \hl{consortia}, as well as other independent data sets\cite{Kasowski2013,Marinov2014,Motamedi2014}.

\section*{Applications of the method}

This protocol was developed and tested primarily for ChIP-seq data. \hl{However, parts of it can be used for other data types or other purposes. The cross-correlation function and metrics can be applied to assess the fragment length distribution and the extent of fragment clustering of a wide variety of data types, while the peak calling procedures centered around reproducibility analysis can be effectively used to identify consensus sets of peaks in chromatin accessibility profiling assays such as DNase-seq}\cite{Crawford2006}, FAIRE-seq\cite{Giresi2007} and ATAC-seq\cite{Buenrostro201320}.

% \hl{SE vs PE}

\section*{Comparison with other methods}

There are several methods for evaluating the \hl{quality} of high-throughput sequencing experiments\cite{Patel2012,Daley2013,Yang2013,}. However many measures, while useful, are not specific to ChIP-seq and provide only limited information about the overall quality of the experiment. For example, the distribution of quality scores along the read positions or the fraction of reads that map confidently to the target genome is often used to determine sequencing quality\cite{Patel2012}, \hl{but provides no information about the success of the ChIP reaction. A} number of ChIP-seq specific quality metrics have been developed\cite{Diaz2012,Planet2012,Mendoza2013}. Diaz and colleagues estimated the strength of ChIP signals by decomposing \hl{them} into real IP signal and background signal, and calculating the percentage of the IP signals\cite{Diaz2012}. Planet et al. \hl{assessed} IP enrichment using the measure of the dispersion of read density by estimating the standard deviation and the Gini index\cite{Planet2012}. Mendoza-Parra and colleagues estimated the stability of peaks by resampling\cite{Mendoza2013}. \hl{However, cross-correlation measures have numerous advantages over these approaches: they directly evaluate the extent of fragment clustering, they provide information about the fragment length distribution of a dataset, and they} can be easily implemented as an intermediate step in a typical peak calling pipeline as we demonstrate here. In addition, researchers can compare the cross-correlation measures with those \hl{for the large collection of TFs that have been profiled by the ENCODE and modENCODE consortia} (\burl{http://genome.ucsc.edu/ENCODE/qualityMetrics.html}).

To ensure that experimental results are reproducible, one may examine the reproducibility either at the level of reads or identified peaks, or both\cite{Landt2012}. The reproducibility of the reads is often measured by computing the Pearson correlation coefficient of the (mapped) read counts at each genomic position\cite{Bardet2012}. However, this quantity can be dominated by a small number of very highly enriched regions. In our approach, we assess the reproducibility of identified peaks between replicates using IDR analysis\cite{Li2011} (Box 1). \hl{The IDR pipeline} reports the number of peaks that pass a user-specified reproducibility threshold (e.g., IDR = 0.05), allowing one to select a reproducible set of peaks and to identify experiments with low reproducibility. A major advantage of the IDR method is that it is rank-based; thus, it is independent of \hl{peak-calling algorithms} and can be applied to a range of significance criteria across labs and platforms. It has been shown to produce a stable threshold that is more consistent across laboratories, antibodies, and analysis protocols (e.g., peak callers) than FDR measures\cite{Li2011}.

\end{multicols}

\begin{FPfigure}
\begin{center}
\includegraphics[width=16.5cm]{Fig2.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf The anatomy of a cross-correlation profile (CCP) for a ChIP-seq dataset}. 
\textbf{(a)} Computation of the cross-correlation function and the biochemical rationale behind its application to ChIP-seq. Transcription factors generally bind to specific sequence motifs precisely positioned within DNA. Crosslinking and sonication result in DNA fragments centered around the site of transcription factor binding (as well as a background population of fragments from the rest of the genome). Sequencing library preparation adapters are randomly ligated to the two ends of these fragments, which results in the typical asymmetric read distribution on the forward and reverse DNA strands around the transcription factor binding site after immunoprecipitation, with the peaks on the two strands being separated by the average DNA fragment size $f$. The cross-correlation function is computed by shifting the read distributions on the two strand relative to each other and computing the cross-correlation function at each shift value. The bottom panel shows that two peaks exhibit maximum overlap when reads on the reverse strands are shifted towards the positive strand (5$^{\prime}$ to 3$^{\prime}$ direction) by $f$.
\textbf{(b)} Cross-correlation profiles (CCP) computed as the Pearson correlation between the positive strand profiles and negative strand profiles at different strand shift distances $d$ (for detailed description, see the Experimental design section and Box 2). Typical CCPs are shown for high-quality, medium-quality, low-quality, and input datasets, respectively, from top to bottom. Blue dashed vertical lines: $d$ = read-length ($r$); blue solid horizontal lines: the magnitude of the read-length peak; red vertical lines: $d = f$; orange solid horizontal lines: the magnitude of the estimated fragment length peaks; gray solid horizontal lines: minimum values of CCP in the given range. Note that the estimated fragment length peak becomes stronger relative to the read length peak, for datasets with greater enrichment. Also note that in the actual implementation of the pipeline cross-correlation is computed over a range of shift values $d \in \left[-500,1500\right]$, but the narrower $d \in \left[0,400\right]$ range is shown here for clarity. The high-quality dataset is BAM file ID ENCFF901DWV from ENCODE experiment ID ENCSR000BMY, the medium-quality one is ENCFF757ZWI from ENCSR000EBQ, the low-quality one is ENCFF000VSZ from ENCSR000EYZ, the control is ENCFF487SIK from ENCSR000EAF.
\textbf{(c)} NSC (Normalized Cross-correlation Coefficient), RSC (Relative Cross-correlation Coefficient) are measures of data quality based on parameters derived from CCPs. FRiP (Fraction of Reads in Peaks) is another measure of ChIP enrichment. The values of these metrics for the datasets shown in (b)) are displayed here.}
\label{Fig2}
\end{FPfigure}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 2 $|$ The definitions of cross-correlation measures (NSC, RSC and Qtag)}}
\centerline{}
\centerline{}
Let $CC(d)$ be the strand cross-correlation function at a shift $d$ (see Figure  \ref{Fig2}). \hl{It is calculated as follows:}

\begin{equation*}
 CC(d) = \sum_{c \in C} \cfrac{|R_{c}|}{\displaystyle\sum_{c \in C}|R_{c}|} CC(d)_{c}
\end{equation*}

where $C$ is the set of all individual chromosomes $c$, and $R_{c}$ is the set of all reads uniquely mapping to a chromosome $c$. 

The normalized strand cross-correlation (NSC) coefficient is defined as:

\begin{equation*}
 NSC = \cfrac{max_{d \in (d_{min},r - \Delta)\cup(r + \Delta,d_{max})}(CC_d)}{min_{d \in (d_{min},d_{max})}(CC_d)}
\end{equation*}

where $r$ is the read length, $\Delta$ is a flanking region near the read length, $d_{min}$ and $d_{max}$ are the minimum and maximum value of strand shift examined, respectively. This coefficient represents the peak height of the CCP ignoring any peak at a strand shift equal to the read length, normalized by the minimum CCP value over a large range of strand shifts. NSC represents the enrichment of localized clustering of relatively fixed-sized IP fragments and is an effective continuous measure of ChIP-seq data quality.
   
The relative strand cross-correlation coefficient (RSC) is defined as:

\begin{equation*}
 RSC = \cfrac{max_{d \in (d_{min},r - \Delta)\cup(r + \Delta,d_{max})}(CC_d) - min_{d \in (d_{min},d_{max})}(CC_d)}{CC(r) - min_{d \in (d_{min},d_{max})}(CC_d)}
\end{equation*}

where $r$ is the read length, $\Delta$ is a franking region near the read length, $d_{min}$ and $d_{max}$ are the minimum and maximum value of strand shift examined, respectively. RSC represents the relative value of the cross-correlation maximum at the estimated fragment length (excluding the read-length peak) to the value at the read length with respect to the minimum cross-correlation across all strand shifts. For convenience, we also define a discretized version we call a quality tag (Qtag) as follows:  

If $RSC \geq 1.5 \rightarrow Qtag = 2$ (Very High)

If $RSC \in [1,1.5) \rightarrow Qtag = 1$ (High)

If $RSC \in [0.5,1) \rightarrow Qtag = 0$ (Medium)

If $RSC \in [0.25,0.5) \rightarrow Qtag = -1$ (Low)

If $RSC < 0.25 \rightarrow Qtag = -2$ (Very Low). \par}
\end{mdframed}
\centerline{}

\begin{multicols}{2}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{FigJSD.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf Summary of and correlation between quality control metrics used in AQUAS.} A set of $\sim$500 ENCODE transcription factor ChIP-seq datasets were processed using the AQUAS pipeline and the quality control metrics part of it were calculated for each
\textbf{(a, b, and c)} Library complexity quality control metrics ($NRF$, $PBC1$, and $PBC2$).
% \textbf{(b)} Comparison of FRiP scores for all ENCODE TF ChIP-seq data sets using two different peak callers and different methods of thresholding peak calls. Left: the SPP and MACS peak callers were used with default FDR thresholds of 1\% to call peaks on all ENCODE TF ChIP-seq data sets. FRiP scores were computed on these sets of peaks. Several data sets (marked in green ellipsoids) exhibited many peaks based on one peak caller but very few to no peaks based on the other peak caller. Middle: similar to the left panel but using the default MACS $p$-value threshold to call peaks. There are several data sets with no peaks based on SPP thresholds but with large number of peaks based on the MACS threshold. Right: SPP and MACS peak calls for all data sets were thresholded based on self-consistency analysis using the IDR framework. FRiP scores were then computed based on the self-consistent peak calls. The instability in the FDR and $p$-value thresholded peak sets is no longer observed. In addition, the correlation between FRiP scores is markedly higher.
\textbf{(d, e, and f)} Read clustering and ChIP quality control metrics based on cross correlation profiles ($NSC$, $RSC$) as well as the synthetic Jensen-Shannon Divergence $JSD$.
}
\label{FigJSD}
\end{figure*}


\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig3.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf The irreproducible discovery rate (IDR) framework for assessing reproducibility of ChIP-seq data sets.} 
% \textbf{(a)} The left figure shows a scatterplot of the significance scores of peaks identified based on two replicates. The IDR method classifies peaks into a reproducible (black) and an irreproducible (red) group, and computes for each peak the probability that the peak belongs to the irreproducible group. Peaks are ranked and selected according to this probability, and IDR, i.e. the amount of expected irreproducible discoveries in the selected peaks, is computed. THe figure on the right shows estimated IDR values at different rank thresholds. These panels were adapted from Bailey et al. \cite{Bailey2013}
\textbf{(a)} Example (ENCODE experiment ID ENCSR000DYI) of the output of the IDR (version 2.0).
% \textbf{(b)} Comparison of FRiP scores for all ENCODE TF ChIP-seq data sets using two different peak callers and different methods of thresholding peak calls. Left: the SPP and MACS peak callers were used with default FDR thresholds of 1\% to call peaks on all ENCODE TF ChIP-seq data sets. FRiP scores were computed on these sets of peaks. Several data sets (marked in green ellipsoids) exhibited many peaks based on one peak caller but very few to no peaks based on the other peak caller. Middle: similar to the left panel but using the default MACS $p$-value threshold to call peaks. There are several data sets with no peaks based on SPP thresholds but with large number of peaks based on the MACS threshold. Right: SPP and MACS peak calls for all data sets were thresholded based on self-consistency analysis using the IDR framework. FRiP scores were then computed based on the self-consistent peak calls. The instability in the FDR and $p$-value thresholded peak sets is no longer observed. In addition, the correlation between FRiP scores is markedly higher.
\textbf{(b and c)} The IDR framework improves consistency between peak calls, minimizes variation and eliminates systematic differences 	caused by the specific parameter and algorithmic choices implemented by different peak calling packages. A set of $\sim$500 ENCODE transcription factor ChIP-seq datasets were processed using the MACS2 and SPP peak callers individually and also using the IDR framework; the FRiP concordance between the two peak callers is shown in \textbf{(b)}, and the FRiP concordance between the two peak callers after applying IDR is shown in \textbf{(c)}. 
}
\label{Fig3}
\end{figure*}

\section*{Experimental Design}

\textbf{Aligning reads.} \hl{Raw sequencing data in FASTQ format} can be obtained \hl{from a sequencing facility or from} public repositories such as GEO (\burl{http://www.ncbi.nlm.nih.gov/geo/}) and \hl{the} ENCODE portal (\burl{encodeproject.org}). \hl{Including a large number of reads with poor quality scores in an analysis} can \hl{potentially} result in mapping artifacts \hl{thus it is advisable to perform a QC step at the level of the raw reads}. Read quality histograms are a simple way to summarize overall read quality. One can compute the histogram of PHRED scores per position in all reads that map to some genomic location as well as the mean score per position. High quality sequenced reads from the Illumina sequencing platform tend to have average PHRED scores $>$30. Average PHRED scores $<$20 reflect problems with read quality. \hl{Reads can then be mapped against a reference genome }using a variety of aligners\hl{, such as BWA}\cite{Li2009}, Bowtie\cite{Langmead2009} and Bowtie2\cite{Langmead2012}. 

\textbf{Mapping and sequencing statistics.} For each dataset, we compute (i) the total number of reads, (ii) the number and fraction of reads that map to the genome, and (iii) the number ($N_u$) and fraction ($f_u$) of uniquely-aligned reads. Human ENCODE datasets with $N_u < 8$ million are flagged as having suboptimal sequencing depth. Datasets with abnormally low $f_u$ values ($<$60\%) tend to have poor read quality or other data quality issues \hl{(although exceptions are possible, for example, a highly enriched ChIP-seq experiment targeting repressive histone modificatons associated with repetitive elements can exhibit a high fraction of multimapping reads, and, accordingly, a low fraction of uniquely mappable reads)}.

\textbf{Post-alignment filtering.} All analyses are subsequently performed on alignment files containing only uniquely mapped reads. We remove \hl{alignments} that are either unmapped or non-primary alignment. The multi-mapped reads are removed using the samtools. 

\textbf{Library complexity measures}. If a sequencing library has an inadequate representation of unique DNA molecules, a single molecule may give rise to a very large number of clones during the PCR amplification process. This PCR ``bottlenecking'' results in many more reads mapping to the same location than expected from a library with greater complexity\hl{, which can potentially lead to peak calling and other analytical artifacts. We use two metrics to assess library complexity}:

\begin{itemize}
\item PCR Bottleneck Coefficient ($PBC$). Let $M_0$, $M_1$, and $M_2$ be the numbers of distinct genomic locations to which at least one, exactly one, and exactly two reads map uniquely, respectively. We define two simple heuristic PCR Bottleneck Coefficients ($PBC1$ and $PBC2$) as follows:

\begin{equation}
PBC1 = M_1/M_0
\end{equation}

\begin{equation}
PBC2 = M_1/M_2
\end{equation}

High $PBC1$ and $PBC2$ scores ($PBC1 \in [0,1]$, $PBC2 \geq 1$) indicate that the dataset has low bottlenecking (e.g., high-molecular complexity), while low scores imply that it has high bottlenecking (i.e. low complexity). 

\item Non-Redundant Fraction ($NRF$). Let $U_P$ be the set of genomic positions to which 5$^{\prime}$ ends of reads map uniquely. Let $U_R$ be the total number of uniquely mapped reads. We define the Non-Redundant Fraction ($NRF$) as follows:

\begin{equation}
NRF = U_P/U_R 
\end{equation}

A high $NRF$ score ($NRF \in \left(0,1\right)]$) indicates that the dataset has a high molecular complexity.
\end{itemize}

\textbf{Converting file formats}. Although the cross-correlation analysis can take BAM files, we recommend to convert them to a tagAlign file format since tagAlign files are more convenient in the downstream process in the IDR analysis.

\textbf{Strand cross-correlation measures}. \hl{To calculate cross-correlation metrics}, stranded coverage vectors are \hl{first} computed for each chromosome by calculating the number of uniquely-mapped high quality reads at each genomic location on the positive and negative strands (Figure \ref{Fig2}). We mask genomic locations that \hl{show} artifactual stacking of reads, i.e., singular genomic locations \hl{with} very high read counts as compared to surrounding locations\cite{Kharchenko2008}. The strand cross correlation profile (CCP) for each chromosome is computed as the Pearson correlation coefficient between coverage vectors of the positive and negative strands shifted by specific distances towards (positive shifts) or away (negative shifts) from each other. We typically use strand shifts ranging from $-500$ bp to 1500 bp at steps of 5 bp. We then compute a weighted average of the CCPs over all chromosomes, using the fraction of reads mapped to each chromosome as weights. All computations are performed using \verb|R| scripts in AQUAS based on the functions in the \verb|SPP| peak caller package\cite{Kharchenko2008}.

We define a Normalized Strand Cross-correlation coefficient (NSC) as the peak height of the CCP divided by a background cross-correlation value which is defined as the minimum cross-correlation over a sufficiently wide range of strand shifts (typically $-500$ bp to $\sim$3 times the expected size of DNA fragments based on sonication and size selection protocols; see also Box 2). NSC represents the enrichment of localized clustering of relatively fixed-size IP fragments at potential target sites, and the normalization with respect to the background makes the measure comparable across datasets. The smallest possible NSC value is 1 and larger values represent greater enrichment. Based on calibration using a comprehensive analysis of control datasets (input DNA), for the human genome, NSC values lower than 1.1 typically represent datasets with low enrichment. Datasets with tight fragment size distributions and a large number of narrow peaks with strong signal compared to background show higher NSC values.
 
Poorly-enriched or failed datasets typically exhibit a low degree of localized fragment clustering. In the extreme case, we observe a predominant peak in the CCP at a strand shift equal to the read length and a significantly weaker peak at a strand shift equal to the predominant fragment length (Figure  \ref{Fig2}b, lower \hl{two panels}). This ``read-length peak'' has been previously noted and attributed to variable local mappability of the genome\cite{Sasaki2009}. To elaborate, if reads map uniquely to a genomic position $x$ on the positive strand, then it implies that reads can also map uniquely to the position ($x$ + read length) on the negative strand\hl{. }Hence, in the absence of any significant fragment clustering and variable mappability across the genome, at random, on average, there is a higher probability of observing reads mapping to positions on the positive and negative strand separated by a distance equal to the read length than any other strand shift \hl{when restricting analysis to uniquely mappable reads only}.
 
This read-length cross-correlation peak can serve as an internal reference point that can be used to normalize the fragment-length peak. Thus, we define the Relative Strand Cross-correlation coefficient (RSC) as the height of the fragment length peak divided by the height of the read-length peak, after adjusting for the background from each (See Box 2). Highly enriched ChIP-seq datasets show RSC values greater than 1 (a dominant fragment-length peak);\hl{ }RSC values markedly smaller than 1 (a dominant read-length peak) \hl{are a sign of potential problems }with ChIP enrichment\hl{ }or sub-optimal sequencing depths. \hl{RSC scores can be discretized into easy to interpret integer scores, which we refer to as ``Qtag'' (see See Box 2).} High quality datasets tend to display higher values of NSC, RSC or Qtag as well as a larger number of detected peaks and a higher fraction of reads residing within the peaks compared to low quality data\hl{sets} (Figure  \ref{Fig2}c).

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18cm]{Fig4.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf Study case 1: Comparison of cross-correlation coefficients between two pairs of replicates from ENCODE datasets.} 
\textbf{(a)} Two pairs of replicates (Rep1 and Rep2) of H3K4me3 ChIP-seq experiments from two ENCODE production groups (Lab1 and Lab2) in the GM12878 cell-line. Genome-browser shots (upper) show normalized ChIP-seq signals (in RPM, Read Per Million). The highlighted regions (in gray) were zoomed in to view of signal at a specific peak location (lower left) and background signal at non-peak regions (lower right). A red box highlights that Rep1 from Lab2 displays higher background noise compared to other replicates.
\textbf{(b)} RSC and NSC scores for the 4 datasets. Replicates from Lab1 show higher NSC and RSC scores than replicates from Lab2. 
\textbf{(c)} Cross-correlation metrics (CCPs) for the 4 datasets.}
\label{Fig4}
\end{figure*}

\textbf{Additional enrichment metrics}. In addition to cross-correlation, we also compute the Jensen-Shannon distance (JSD) for a sample as previously defined\cite{Diaz2012b} and implemented by the \verb|deepTools|\cite{Ramirez2016} \hl{package. The JSD metric is defined as follows:}

\begin{equation*}
JS = \sqrt{JSD(P||Q)}
\end{equation*}

Where $JSD$ is the Jensen-Shannon divergence:

\begin{equation*}
JSD(P||Q) = \cfrac{1}{2}D_{KL}(P||M) + \cfrac{1}{2}D_{KL}(Q||M)
\end{equation*}

where 

\begin{equation*}
M = \cfrac{1}{2}(P+Q)
\end{equation*}

and

\begin{equation*}
D_{KL}(X||Y) = \displaystyle\sum_{i}X_i log\left(\cfrac{X_i}{Y_i}\right)
\end{equation*}

Where $D_{KL}$ is the Kullback-Leibler divergence, and $P$ and $Q$ are the two statistical distributions; in this case the JS metric evaluates the difference between the evaluated sample and a ``perfect'' (i.e. a Poisson probability distribution with $\lambda$ equal to the mean coverage of the sample) input sample sequenced to the same depth. The distributions compared are the following: for both the actual and the ideal input sample, the genome is split into bins $b_1,....,n$, and for each bin $b_i$, the ratio of sequencing reads for this bin relative to bin with the highest coverage, i.e. $y_i = b_i/max_{1,...,n}(b_i)$ is calculated; bins are then ranked and the resulting rank versus fractional coverage $y_i$ curves are used to evaluate the JS divergence and calculate the synthetic JS distance. % (Figure XXXXX). 

\textbf{Determining reproducible peaks using IDR analysis.} The IDR framework (Box 1; Figure  \ref{Fig3}) is used to obtain more stable, objective estimates of the number of peaks or regions of enrichment based on the reproducibility across biological replicates\cite{Landt2012,Li2011,Bailey2013}. To evaluate\hl{ }reproducibility across biological replicates, a peak caller is used to call peaks with relaxed thresholds on each biological replicate. If a peak in one replicate overlaps with a peak in another replicate by at least 1 bp, then the pair is considered as calling the same region. The significance scores of such pairs (peak signal score for SPP and peak $p$-value for MACS) \hl{are then} used as \hl{input to the} IDR framework to evaluate the reproducibility.  The number of peaks that pass a use\hl{r}-specified IDR threshold is recorded and is used to determine the number of reproducible peaks to keep. 
 
For ENCODE datasets, we \hl{use} SPP\cite{Kharchenko2008} and MACS\cite{Zhang2008} with FDR = 0.9 for SPP and $p$-value $\leq$ 0.01 for MACS, to obtain 150,000 to 300,000 peaks. We use\hl{ }the default peak width for SPP and a fixed peak width of +/- 200 bp around the peak summit for MACS since MACS tends to call wider peaks at relaxed peak calling thresholds.  We determine\hl{ }the number of reproducible peaks ($N$) at IDR=0.05 and \hl{keep} the most significant $N$ peaks on each replicate sample. 

\end{multicols}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 3 $|$ Cross-correlation analysis in AQUAS}}
\centerline{}
\centerline{}
1. Input file format

 (1) BAM format
 
This is a binary alignment format specified in \burl{https://genome.ucsc.edu/goldenPath/help/bam.html}

 (2) TagAlign files
 
 This a text-based BED3+3 alignment format that is easier to manipulate. It contains 6 tab delimited columns.
\centerline{}
\centerline{}
 \indent COL1: \verb|chrom       string  | Name of the chromosome
 
 COL2: \verb|chromStart  int     | The starting position of the feature in the chromosome. The first base in a chromosome is numbered 0.
 
 COL3: \verb|chromEnd    int     | The ending position of the feature in the chromosome or scaffold. The \verb|chromEnd| base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as \verb|chromStart| = 0, \verb|chromEnd| = 100, and span the bases numbered 0--99.
 
 COL4: \verb|sequence    string  | Sequence of this read
 
 COL5: \verb|score       int     | Indicates uniqueness or quality (preferably 1000/\verb|alignmentCount|).
 
 COL6: \verb|strand      char     | Orientation of this read (+ or -)  
\centerline{}
\centerline{}
2. Output file 

The out file is a tab delimited text file, containing the following columns.
\centerline{}
\centerline{}
\indent COL1: \verb| Filename        |: tagAlign/BAM filename

COL2: \verb| numReads        |: effective sequencing depth i.e. total number of mapped reads in input file

COL3: \verb| estFragLen      |: comma separated strand cross-correlation peak(s) in decreasing order of correlation.

\verb|      | The top 3 local maxima locations that are within 90% of the maximum cross-correlation value are output.

\verb|      | In almost all cases, the top (first) value in the list represents the predominant fragment length.

\verb|      | If you want to keep only the top value simply run

\verb|      |\verb|sed -r 's/,[^\t]+//g' <outFile> > <newOutFile>|

COL4: \verb| corr_estFragLen |: comma separated strand cross-correlation value(s) in decreasing order (col2 follows the same order)

COL5: \verb| phantomPeak     |: Read length/phantom peak strand shift

COL6: \verb| corr_phantomPeak|: Correlation value at phantom peak

COL7: \verb| argmin_corr     |: strand shift at which cross-correlation is lowest

COL8: \verb| min_corr        |: minimum value of cross-correlation

COL9: \verb| | Normalized strand cross-correlation coefficient (NSC) = COL4 / COL8

COL10: \verb| |Relative strand cross-correlation coefficient (RSC) = (COL4 - COL8) / (COL6 - COL8)

COL11: \verb|QualityTag      |: Quality tag based on thresholded RSC (codes: -2:veryLow, -1:Low, 0:Medium, 1:High, 2:veryHigh)
\centerline{}
\centerline{}
3. Parameter options

Optional Arguments 
\centerline{}
\centerline{}
\indent \verb|-s=<min>:<step>:<max>|, strand shifts at which cross-correlation is evaluated, default=\verb|-500:5:1500|

\verb|-speak=<strPeak>|, user-defined cross-correlation peak strandshift

\verb|-x=<min>:<max>|, strand shifts to exclude (This is mainly to avoid region around phantom peak) default=\verb|10:(readlen+10)|

\verb|-p=<nodes>|, number of parallel processing nodes, default=0

\verb|-fdr=<falseDisoveryRate>|, false discovery rate threshold for peak calling

\verb|-npeak=<numPeaks>|, threshold on number of peaks to call

\verb|-tmpdir=<tempdir>|, Temporary directory (if not specified R function \verb|tempdir()| is used)

\verb|-filtchr=<chrnamePattern>|, Pattern to use to remove tags that map to specific chromosomes e.g. \verb|_\verb| will remove all tags that map to chromosomes with \verb|_| in their name
\centerline{}
\centerline{}
\indent Output arguments 
\centerline{}
\centerline{}
\indent \verb|-odir=<outputDirectory>| name of output directory (If not set same as ChIP file directory is used)

\verb|-savn=<narrowpeakfilename>| OR \verb|-savn| NarrowPeak file name (fixed width peaks)

\verb|-savr=<regionpeakfilename>| OR \verb|-savr| RegionPeak file name (variable width peaks with regions of enrichment around peak summits)

\verb|-savd=<rdatafile>| OR \verb|-savd|, save Rdata file

\verb|-savp=<plotdatafile>| OR \verb|-savp|, save cross-correlation plot

\verb|-out=<resultfile>|, append peakshift/phantomPeak results to a file

\verb|-rf|, if plot or rdata or narrowPeak file exists replace it. If not used then the run is aborted if the plot or Rdata or narrowPeak file exists

\verb|-clean|, if used it will remove the original chip and control files after reading them in. CAUTION: Use only if the script calling. \par}
\end{mdframed}
\centerline{}

\begin{multicols}{2}


\textbf{Flagging \hl{poorly reproducible datasets }using IDR analysis.} In the event that one replicate is significantly worse than the other replicate, the reproducibility between the replicates is low and the number of reproducible peaks is lower bounded by the worse replicate. To differentiate \hl{such cases from the scenario in which the targeted protein genuinely has a small number of occupancy} sites, we \hl{pool} the uniquely aligned reads from the biological replicates and divide the pooled reads into two randomly sampled, non-overlapping sets, referred to as \hl{``pseudo-replicates''. We then assess reproducibility between pseudo-replicates using IDR as described above with some modifications (See Box 6)}.  If the number \hl{of} reproducible peaks identified on the pooled pseudo replicates ($N_p$) is significantly more than the number of reproducible peaks identified on the biological replicates ($N_t$), the dataset is flagged \hl{for poor reproducibility}. For ENCODE data, a dataset is flagged if $N_p/N_t > 2$. 

\textbf{Fraction of Reads in Peaks (FRiP).} Peak callers are typically used to identify statistically significant peaks or regions of enrichment, with a false discovery rate (FDR)\cite{Kharchenko2008,Rozowsky2009} or $p$-value threshold\cite{Zhang2008} for selecting high confidence sites. One measure of IP enrichment is the fraction of reads that lie within peaks or regions of enrichment, \hl{referred to as} ``fraction of reads in peaks'' (FRiP)\cite{Landt2012}. These values \hl{can} range from 0 to 1 \hl{(although in practice values significantly above 0.5 are very rarely observed)}, with higher values indicating greater IP enrichment. The NSC and RSC thresholds were calibrated using input DNA datasets as a gold standard negative set. For FRiP scores, calibration is difficult because peak calling is ill-defined for the input control. Hence, we used NSC and RSC thresholds to define a corresponding FRiP score threshold. We \hl{have} found that a self-consistency FRiP threshold of 1\% works well for human ENCODE datasets and showed consistency with NSC and RSC thresholds.

\section*{Study cases}

\subsection*{Study case 1:}

\hl{We illustrate how cross-correlation analysis can be used to }

\subsection*{Study case 2:}

XXXX \hl{We} illustrate how\hl{ }cross-correlation analysis can be used to \hl{provide} feedback \hl{to experimentalists with ChIP-seq experiments targeting the transcriptional coactivator p300} in the HepG2 cell\hl{ }line (Figure  \ref{Fig5}). The cross-correlation measures flagged one of the replicates (Figure  \ref{Fig5}b lower left). This \hl{was} confirmed by visualization of the signal tracks. The production lab promptly replaced these datasets with new ones. The new datasets showed high NSC/RSC scores comparable to the high quality older dataset (Figures  \ref{Fig5}b and  \ref{Fig5}c). Hence, cross-correlation analysis can be used to provide data producers \hl{with} instant feedback about data quality and effectively flag potentially low-quality datasets for replacement resulting in higher-quality data\hl{sets}. XXXX 

\subsection*{Study case 3:}

XXXX  \hl{We} illustrate how IDR can be used to identify suboptimal \hl{datasets with an ENCODE ChIP-seq dataset targeting the YY1 transcription factor, generated in the GM12878} cell line. This dataset consists of two pairs of biological replicates, each of which was generated by \hl{a different} lab (Lab 1 and Lab 2, respectively). \hl{Peaks were called} using SPP and the reproducibility of identified peaks on the replicates within a lab \hl{was} assessed using IDR. At any given IDR \hl{threshold}, a much smaller number of reproducible peaks \hl{were} reported on samples from Lab 1 than on the samples from Lab 2, as shown in Figure  \ref{Fig6} \hl{(e.g. 1,939 in Lab 1 vs 22,286 for Lab 2 at IDR = 0.01)}. The small number of reproducible peaks indicates that the samples from Lab 1 \hl{are of} relatively low quality, which is also reflected by other evidence, such as a high proportion of peaks that are identified only on one replicate \hl{(86.7\% in Lab 1 vs 63.6\% in Lab 2).} XXXX 

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig5.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf Study case 2: Identifying and replacing low quality p300 ENCODE ChIP-seq datasets resulting in significant improvement in data quality and reproducibility.} 
\textbf{(a)} Genome-browser shots of normalized ChIP-seq signal (RPM) targeting the p300 enhancer binding protein. The lower-panel shows a zoomed view of signals at a specific peak location (highlighted in gray). The bottom two tracks represent two replicates where one replicate (old Rep2) was found to show significantly lower NSC and RSC scores due to a strong read-length peak in the CCP and significantly weaker ChIP-seq peak signal compared to the other replicate (old Rep1). These datasets were reported to the production group and replaced with new replicates (new Rep1 and Rep2).
\textbf{(b)} CCPs generated from datasets in (a). Replaced datasets show substantially higher fragment length peaks compared to the old dataset (old Rep2). 
\textbf{(c)} Corresponding RSC, NSC and FRiP scores. A significant increase of RSC, NSC and FRiP scores of replaced datasets indicates the improvement of data quality of the datasets.}
\label{Fig5}
\end{figure*}

\section*{Limitations}

The quality control metrics and IDR analysis \hl{described here} have been largely tested on TFs \hl{and on or} active histone marks which have a \hl{localized} binding pattern, such as H3K4me3 and H3K27ac. Due to the nature of broad histone modifications, determining enriched regions is challenging compared to \hl{peak calling for TFs and is not yet a fully solved problem}. Thus, this protocol is not optimized for broad marks. We are \hl{in the process of} developing the analysis pipeline for broad histone modifications. 

\hl{We also note that i}f the read-length for a dataset is within the range of the predominant fragment length, then the read-length cross-correlation peak can be hard to distinguish from the fragment-length cross-correlation peak. \hl{The evaluation of read clustering by cross-correlation also relies on the assumption that reads on opposite strands are independent of each other. This assumption is violated if paired-end data is to be analyzed by cross-correlation, as in such cases there is always another read at roughly a mean fragment length distance away from each read; as a result cross-correlation are always very high on paired-end data and provide little meaningful information regarding ChIP quality. Of note, this issue cannot be resolved by only considering the first end in mapped read pairs if reads have been mapped in a paired-end format. This is because, as discussed above, the phantom peak arises due to local mappability variations in the genome; when reads are mapped as pairs, the additional information provided by the second read in a pair makes many of the repetitive regions in the genome that give rise to the phantom peak in the cross-correlation function uniquely mappable, thus suppressing the phantom peak relative to the fragment-length peak. Another consideration to keep in mind is that if reads are trimmed to a variable length by removing adapters, and this results in a significant proportion of alignments being of less than the maximum length, the phantom peak will again be lowered somewhat compared to the fragment length peak, for the same reasons having to do with unique mappability in the genome discussed above. Finally, the assumptions behind cross-correlation analysis are also violated if multimapping reads are included in the input alignment files, and it should therefore only be carried out on unique alignments alone for the purposes of ChIP-seq quality control.}

% \filllastline{ }

\end{multicols}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 4 $|$ Code description for IDR analysis}}
\centerline{}
\centerline{}
1. Input file formats

(1) Peak files (specified with \verb|--samples [PEAK_FILE_1] [PEAK_FILE_2]|)

Peak files must be in narrowPeak format. NarrowPeak files are in BED6+4 format. It consists of 10 tab-delimited columns
\centerline{}
\centerline{}
\indent COL1:\verb|  chrom       string| Name of the chromosome

COL2:\verb|	 chromStart  int   | The starting position of the feature in the chromosome. The first base in a chromosome is numbered 0.

COL3:\verb|	 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 \verb|chromStart=0|, \verb|chromEnd=100|, and span the bases numbered 0-99.

COL4:\verb|	 name        string| Name given to a region (preferably unique). Use '.' if no name is assigned

COL5:\verb|	 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 \verb|signalValue| per base spread between 100-1000.

COL6:\verb|	 strand	     char  | +/- to denote strand or orientation (whenever applicable). Use '.' if no orientation is assigned.

COL7:\verb|	 signalValue float | Measurement of overall (usually, average) enrichment for the region.

COL8:\verb|	 pValue      float | Measurement of statistical signficance ($-log_{10}$). Use -1 if no \verb|pValue| is assigned.

COL9:\verb|	 qValue      float | Measurement of statistical significance using false discovery rate ($-log_{10}$). Use -1 if no \verb|qValue| is assigned.

COL10:\verb| peak        int   | Point-source called for this peak; 0-based offset from \verb|chromStart|. Use -1 if no point-source called.
\centerline{}
\centerline{}
\indent The $p$-value and $q$-value columns should be at $-log_{10}$ scale. 

The narrowPeak format has 3 columns that can be used to rank peaks (\verb|signal.value|, \verb|p.value| ($-log_{10}$) and \verb|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 which measure you are using rank peaks is relatively continuous without too many ties. You can choose a measure with \verb|--rank [RANK_MEASURE]| Recommend to use \verb|signal.value|, \verb|p.value| and \verb|q.value| for the SPP, MACS2 and PeakSeq peak callers, respectively.

(2) Oracle peak list (specified with \verb|--peak-list [ORACLE_PEAK_LIST]|)

Oracle peak list is usually peaks called from merged replicates. For each oracle peak a single peak from each replicate is chosen that overlaps the oracle peak. If there are multiple peaks that overlap the oracle, then ties are broken by applying the following criteria in order: 1) choose the replicate peak with a summit closest to the oracle peak's summit 2) choose the replicate peak that has the largest overlap with the oracle peak 3) choose the replicate peak with the highest score. If an oracle peak list is not given, Peaks are grouped by overlap and then merged. The merged peak aggregate value is determined by \verb|--peak-merge-method|. Peaks that don't overlap another peak in every other replicate are not included unless \verb|--use-nonoverlapping-peaks| is set.

2. Output file formats

(1) IDR peak

The output format mimics the input file type, with some additional fields. The first 6 columns are a standard BED6, the first 10 columns are a standard narrowPeak. Also, for columns 7-10, only the score that the IDR code used for ranking will be set -- the remaining two columns will be set to -1. \verb|--idr-threshold [IDR_THRESHOLD]| only returns peaks with a global IDR threshold below this value while \verb|--soft-idr-threshold [IDR_THRESHOLD]| reports statistics (including plots) for peaks with a global IDR below this value but return all peaks.
\centerline{}
\centerline{}
\indent COL11:\verb|  localIDR          float| $-log_{10}$(Local IDR value)

COL12:\verb|  globalIDR         float| $-log_{10}$(Global IDR value)

COL13:\verb|  rep1_chromStart   int  | The starting position of the feature in the chromosome or scaffold for common replicate 1 peaks, shifted based on offset. The first base in a chromosome is numbered 0.

COL14:\verb|  rep1_chromEnd     int  | The ending position of the feature in the chromosome or scaffold for common replicate 1 peaks. The chromEnd base is not included in the display of the feature.

COLl15:\verb|  rep1_signalValue float| Signal measure from replicate 1. Note that this is determined by the \verb|--rank| option. e.g. if \verb|--rank| is set to \verb|signal.value|, this corresponds to the 7th column of the narrowPeak, whereas if it is set to \verb|p.value| it corresponds to the 8th column.

COL16:\verb|  rep1_summit       int  | The summit of this peak in replicate 1.
\centerline{}
\centerline{}
(2) Plots

Four plots are generated in a single PNG format suffixed by \verb|.png|. These plots are informative about the quality and similarity of the replicates.

Upper Left: Replicate 1 peak ranks versus replicate 2 peak ranks -- peaks that do not pass the specified IDR threshold are colored red.

Upper Right: Replicate 1 $log_{10}$ peak scores versus replicate 2 $log_{10}$ peak scores -- peaks that do not pass the specified IDR threshold are colored red.

Bottom Row: Peaks rank versus IDR scores are plotted in black. The overlaid boxplots display the distribution of IDR values in each 5\% quantile. The IDR values are thresholded at the optimization precision 1e-6 by default. \par}
\end{mdframed}
\centerline{}

\begin{multicols}{2}

\hl{For these reasons, when working with single-end reads longer than 50 bp, and when working with paired-end reads, we carry out peak calling, IDR analysis, and signal track generation on the full-length reads, but we map reads separately in a single-end 50 bp format for the purposes of cross-correlation analysis.}

\hl{Another important consideration is that the thresholds presented in this protocol are mainly focused on mammalian ENCODE datasets. Other genomes can exhibit very different properties (due to smaller or larger size, a very different repetitive element and unique mappability structure, and other factors), and the parameters used for mammalian genomes may not be appropriate in such cases. Careful calibration for each new genome needs to be carried out in such situations.}

\hl{Finally, the quality control metrics described here are designed with conventional ChIP-seq datasets in mind. In recent years, novel methods for ChIP library preparation, such as ChIPmentation}\cite{Schmidl2015}, TCL\cite{Zarnegar2017}, and others, novel variations of the ChIP assay, such as ChIP-exo\cite{Rhee2011} and ChIP-nexus\cite{He2015}, and alternatives to the ChIP assay, such as ChEC-seq\cite{Zentner2015} and CUT\&RUN\cite{Skene2017}\hl{ have been described. While the general IDR peak calling and reproducibility framework is applicable to such datasets, the correspondence between fragment ends and protein occupancy sites in the libraries generated by these assays may not conform to the assumptions of the cross-correlation quality control analysis as designed for conventional ChIP libraries, therefore cross-correlation results should be interpreted with caution with respect to the degree of target enrichment in such datasets.}

\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig6.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\caption{{\bf Study case 3: Comparison of IDR results and identified peaks for four YY ChIP-seq replicates from two ENCODE production groups.} 
\textbf{(a)} Concordance of scores and the ranks of scores on the replicates from each lab. Only peaks that are identified on both replicates are shown. Black: peaks selected at IDR$\leq$0.01; red: peaks with IDR$>$0.01. For Lab2 samples, only the top 30,000 peaks are plotted in the rank scatter plot to keep the plot on the same scale as Lab1 samples. 
\textbf{(b)} IDR curve.}
\label{Fig6}
\end{figure*}


\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=18.5cm]{Fig7.png}
\end{center}
\captionsetup{singlelinecheck=off,justification=justified}
\cprotect\caption{{\bf Distribution of quality control metrics over ENCODE datasets.} Summarized here is the set of $\sim$500 transcription factor datasets (corresponding to $\sim$1000 individual ChIP-seq experiments) from one of the main ENCODE production labs also shown in previous figures.
\textbf{(a)} Library complexity metrics
\textbf{(b)} ChIP enrichment metrics evaluating the extent of read clustering in datasets.}
\label{Fig7}
\end{figure*}

\end{multicols}

\section*{Materials}

\subsection*{Equipment}

\subsection*{Operating system}

This protocol assumes it runs on a Unix-like operating system (e.g., Linux). Although each step could be run on a Windows system, the provided scripts and protocols are based on a Unix-like operating system with bash shell commands.

\subsection*{Hardware requirement}

A system with high speed internet connection, Intel/AMD 64-bit CPU cores $\geq 4$, memory $\geq 16$ GB and disk space $\geq$ 300 GB.

% \begin{figure*}
% \begin{center}
% \includegraphics[width=18cm]{Fig8.png}
% \end{center}
% \captionsetup{singlelinecheck=off,justification=justified}
% \caption{{\bf Graphical outputs generated from the IDR procedure for illustrating the reproducibility of the samples.} 
% \textbf{(a)} Correspondence curves of the three pairs of replicates for all the identified peaks.  
% \textbf{(b)} Similar in (a) but for the peaks that are identified on both replicates. Among all curves, the one for Rep1-Rep2 show the largest departure from the line of slope 1 in (a) and the strongest parabola shape in (b), which indicates that Rep1-Rep2 is the least reproducible pair among all the three pairs.
% \textbf{(c)} IDR curve: Irreproducible discovery rate (IDR) at different number of selected peaks, plotted at various IDR cutoffs for each pair of replicates. Only peaks that are commonly identified are included. At a given IDR threshold, Rep1-Rep2 has the least number of reproducible peaks, which confirms that Rep1-Rep2 pair is the least reproducible pair.
% \textbf{(d)} The change of correspondence curve for the three pairs of replicates for all the identified peaks. 
% \textbf{(e)} Similar in (d) but for the peaks that are commonly identified on both replicates. Although the breakdown point is difficult to determine, Rep1-Rep2 shows the largest departure to the pattern at perfect consistency, which agrees with (a-c). Prototypical shapes and interpretation of the curves can be found in Box 5.}
% \label{Fig8}
% \end{figure*}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 5 $|$ Interpretation of the output of the IDR statistical analysis}}
\centerline{}
\centerline{}
The output of the IDR statistical analysis consists of two main elements, a graphical tool and the irreproducible discovery rate for each entry on the rank lists. 

The graphical tool provides a quick illustration for visualizing and inspecting the reproducibility of the replicates and the transition from reproducible to irreproducible signals, without making any model assumptions. This tool is convenient for diagnosis and quality control. This tool provides two curves, namely, the correspondence curve and the change of correspondence curve. For the correspondence curve, a perfect rank consistency between peaks shows as a straight line of slope 1, and lack of consistency shows as departures from the diagonal line, such as curvature bending toward the $x$-axis. If almost no association is present, the curve shows a parabolic shape. For the change of correspondence curve, a perfect consistency shows as a straight line of slope 0 with intercept 1, and lack of consistency shows as a line with a positive slope. The transition of the shape of the curves, if present, indicates the breakdown of consistency, which can be used for determining the threshold. Prototypical plots and detailed interpretation of the two curves can be found in Li et al.\cite{{Li2011}}.  The patterns in real datasets usually are more complicated than the prototypical ones. Yet the prototypical plots are still helpful for interpreting the actual patterns.  

Both the local IDR and global IDR are reported in the output for peaks that are commonly identified on replicates. To illustrate the relationship between the IDR thresholds and the number of top peaks, an IDR curve is produced to plot the numbers of selected peaks at various IDR cutoffs. Reproducible peaks can be selected according to a user-specified IDR threshold. \par}
\end{mdframed}
\centerline{}

\subsection*{Software}

\begin{itemize}
\item \textbf{Conda} (Anaconda Python $\geq$ 4.0 or Miniconda $\geq$ 4.0.5): Conda is a cross-platform package/environment manager. All dependencies for this protocol are installed and managed by Conda. Conda is included in Anaconda Python or Miniconda. Softwares installed by Conda are listed in Box X with their version information and description for manual installation.

\item \textbf{BigDataScript} (BDS $\geq$0.999999l): BDS is a cross-platform workflow language that manages subtasks of this protocol. Such subtasks include mapping of raw FASTQs, filtering of mapped BAMs, peak-calling on filtered BAM/\hl{tagAlign} files, and post peak-calling procedures such as IDR. This protocol can also work without BDS by manually running each subtask in a correct order. However, BDS is strongly recommended to automate the entire procedure as \hl{a single-command line} solution. BDS allows a pipeline to be run on stand-alone machines or compute clusters with different submission engines (Sun Grid Engine and SLURM). BDS automatically creates a dependency tree for its subtasks such that a failed pipeline can be resumed at the point of failure and all subtasks can be run in parallel efficiently. \hl{A d}etailed HTML report including task graph and quality metrics is also provided at the end of the pipeline run.

\item \textbf{Java} (JRE $\geq$ 1.8 or JDK $\geq$ 1.8): BDS runs on Java Runtime $\geq$ 1.8.
\end{itemize}

\subsection*{Equipment Setup}

\textbf{Java:} You need to have a super-user privilege to install Java on your system or locally install Java and make sure that \hl{the} Java binary path is appended to your \verb|$PATH| for all login/interactive bash shells.

\begin{itemize}
\item For Debian/Ubuntu ($>$14.10) based Linux:
\begin{verbatim}$ sudo apt-get install git openjdk-8-jre
\end{verbatim}
\item For Fedora/Red-Hat based Linux:
\begin{verbatim}
$ sudo yum install git java-1.8.0-openjdk
\end{verbatim}
\item For Ubuntu 14.04 (Trusty):
\begin{verbatim}
$ sudo add-apt-repository ppa:webupd8team/java -y
$ sudo apt-get update
$ sudo apt-get install oracle-java8-installer
\end{verbatim}
\end{itemize}

\textbf{Conda}: Install Miniconda3 (= 4.0.5 or $\geq$ 4.3.12) on your system. Answer \verb|yes| for the final question. If you choose \verb|no|, you need to manually add Miniconda3 to your shell startup script (\verb|$HOME/.bashrc| or \verb|$HOME/.bash_profile|).

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box 6 $|$ Recommendation for the IDR thresholds}}
\centerline{}
\centerline{}
For self-consistency and comparison of true replicates:

\begin{itemize}
\item If you started with $\sim$150 to 300K relaxed pre-IDR peaks for large genomes (human/mouse), then threshold of 0.01 or 0.02 generally works well. 
\item If you started with $<$ 100K pre-IDR peaks for large genomes (human/mouse), then threshold of 0.05 is more appropriate. This is because the IDR sees a smaller noise component and the IDR scores get weaker. This is typically for use with peak callers that are unable to be adjusted to call large number of peaks (e.g., PeakSeq).
\item For smaller genomes such as worm, if you start with $\sim$15K to 40K peaks then once again IDR thresholds of 0.01 or 0.02 work well.
\item For self-consistency analysis of datasets with shallow sequencing depths, you can use an IDR threshold as relaxed as 0.1 if you start with $<$ 100K pre-IDR peaks.
For pooled-consistency analysis
\item If you started with $\sim$150 to 300K relaxed pre-IDR peaks for large genomes (human/mouse), then threshold of 0.0025 or 0.005 generally works well. We use a tighter threshold for pooled-consistency since pooling and subsampling equalizes the pseudo-replicates in terms of data quality. So we err on the side of caution and use more stringent thresholds. The equivalence between a pooled-consistency threshold of 0.0025 and original replicate consistency threshold of 0.01 was calibrated based on a gold-standard pair of high quality replicate datasets for the CTCF transcription factor in human.
\item You can also use 0.01 for less stringent thresholds.
\item If you started with $<$ 100K, pre-IDR peaks for large genomes (human/mouse), then a threshold of 0.01 is more appropriate. Since the IDR sees a smaller noise component and the IDR scores become weaker, you should relax the thresholds. This is typically for use with peak callers that are unable to be adjusted to call large number of peaks (e.g. PeakSeq)
\item For smaller genomes such as worm, if you start with $\sim$15K to 40K, the IDR thresholds of 0.01 work well.
\end{itemize}
}
\end{mdframed}
\centerline{}

\begin{verbatim}
$ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ bash Miniconda3-latest-Linux-x86_64.sh
\end{verbatim}

Remove any other Anaconda Python from your \verb|$PATH|. Check your loaded modules with a command \verb|$ module list| and unload any Anaconda Python modules in your shell initialization script (\verb|$HOME/.bashrc| or \verb|$HOME/.bash_profile|). Open a new terminal after installation.

\textbf{BigDataScript:} Install BigDataScript v0.99999e (forked) on your system.

\begin{verbatim}
$ wget https://github.com/leepc12/BigDataScript/blob/master/distro/bds_Linux.tgz?raw=true 
  -O bds_Linux.tgz
$ mv bds_Linux.tgz $HOME && cd $HOME && tar zxvf bds_Linux.tgz
\end{verbatim}

Add \verb|export PATH=$PATH:$HOME/.bds| to your initialization script (\verb|$HOME/.bashrc| or \verb|$HOME/.bash_profile|). If Java memory occurs, add \verb|export| \verb|_JAVA_OPTIONS="-Xms256M| \verb|-Xmx728M| \verb|-XX:ParallelGCThreads=1| too.

\textbf{AQUAS pipeline (this protocol):} Get the latest version of this protocol.

\begin{verbatim}
$ git clone https://github.com/kundajelab/chipseq_pipeline
\end{verbatim}

\textbf{Dependencies:} Run the following command and then it creates two Conda environments (\verb|aquas_chipseq| and \verb|aquas_chipseq_py3|).

\begin{verbatim}
$ bash install_dependencies.sh
\end{verbatim}

\textbf{Genome data:} This protocol provides a genome database installer for various species (\verb|mm9|, \verb|hg19|, \verb|hg38|, \verb|mm10|, \verb|hg38| and \verb|mm10_ENCODE|). Use \verb|hg38| and \verb|mm10_ENCODE| for genome data associated with the ENCODE project. The genome database includes reference sequences, \verb|bwa| indexes, chromosome sizes file, unique mappability tracks, and blacklisted regions for each species. Run the following command and then it automatically updates a configuration file (\verb|default.env|) on the pipeline \verb|git| directory. It is recommended to use the same \verb|[DESTINATION_DIR]| for multiple species.

\begin{verbatim}
$ bash install_genome_data.sh [SPECIES_NAME] [DESTINATION_DIR]
\end{verbatim}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18.5cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box V $|$ Configuration JSON files for two test cases (ENCSR000DYI for SE and ENCSR936XTK for PE).}}
\centerline{}
\begin{verbatim}
ENCSR000DYI_SE.json
{
        "out_dir" : "ENCSR000DYI",
        "title" : "ENCSR000DYI_SE_CEBPB_A549",
        "se" : true,
        "fastq1" : "rep1.fastq.gz",
        "fastq2" : "rep2.fastq.gz",
        "ctl_fastq1" : "ctl1.fastq.gz",
        "ctl_fastq2" : "ctl2.fastq.gz",
        "species" : "hg38",
        "nth" : 8,
        "use_pooled_ctl" : true
}
ENCSR936XTK_PE.json
{
        "out_dir" : "ENCSR936XTK”,
        "title" : "ENCSR936XTK_PE_ZNF143_GM12878",
        "pe" : true,
        "fastq1_1" : "rep1-R1.fastq.gz",
        "fastq1_2" : "rep1-R2.fastq.gz",
        "fastq2_1" : "rep2-R1.fastq.gz",
        "fastq2_2" : "rep2-R2.fastq.gz",
        "ctl_fastq1_1" : "ctl1-R1.fastq.gz",
        "ctl_fastq1_2" : "ctl1-R2.fastq.gz",
        "ctl_fastq2_1" : "ctl2-R1.fastq.gz",
        "ctl_fastq2_2" : "ctl2-R2.fastq.gz",
        "species" : "hg38",
        "nth" : 16,
        "use_pooled_ctl" : true
}
\end{verbatim}
}
\end{mdframed}
\centerline{}

\textbf{Installation for internet-free computers:} AQUAS pipeline in this protocol does not need internet connection but installers (\verb|install_dependencies.sh| and \verb|install_genome_data.sh|) do need it. Prepare a computer that has the same OS version as in the internet-free computer and follow the instruction above. Move your miniconda3, genome database and BDS binary (\verb|$HOME/.bds|) directories to the home directory on the internet-free computer. On your internet-free computer, install Java $\geq$ 1.8 and modify the [default] section in \verb|default.env| (\verb|conda_bin_dir=$HOME/miniconda3/bin|; \verb|species_file=$HOME/genome_data/aquas_chipseq_species.conf|) and make sure that all paths in the species configuration file (\verb|aquas_chipseq_species.conf|) in the genome database directory correctly point to the right files. See more details on \burl{https://github.com/kundajelab/chipseq_pipeline#installation-for-internet-free-computers}

\section*{PROCEDURE}

Move to the root directory of outputs for test cases and then download the latest pipeline for this protocol.

\begin{verbatim}
$ git clone https://github.com/kundajelab/chipseq_pipeline
\end{verbatim}

Please make sure that you have installed Java, Conda, BigDataScript and other dependencies described in the previous section. For the demonstration, download raw reads for two replicates of ChIP-seq experiments (\burl{https://www.encodeproject.org/experiments/ENCSR936XTK/}).

\begin{verbatim}
$ wget https://www.encodeproject.org/files/ENCFF960TNP/@@download/ENCFF960TNP.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF640CBP/@@download/ENCFF640CBP.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF246DIP/@@download/ENCFF246DIP.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF616WSS/@@download/ENCFF616WSS.fastq.gz
\end{verbatim}

Download replicates for controls (\burl{https://www.encodeproject.org/experiments/ENCSR398JTO/}).

\begin{verbatim}
$ wget https://www.encodeproject.org/files/ENCFF002EFT/@@download/ENCFF002EFT.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF002EFQ/@@download/ENCFF002EFQ.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF002EFU/@@download/ENCFF002EFU.fastq.gz
$ wget https://www.encodeproject.org/files/ENCFF002EFS/@@download/ENCFF002EFS.fastq.gz
\end{verbatim}

There are two replicates (indices 2 and 3) for IP and two replicates for controls. For convenience, let's rename sample files and change indices for IP replicates (2 and 3 to 1 and 2).

\begin{verbatim}
$ mv ENCFF960TNP.fastq.gz rep1-R1.fastq.gz; mv ENCFF640CBP.fastq.gz rep1-R2.fastq.gz
$ mv ENCFF246DIP.fastq.gz rep2-R1.fastq.gz; mv ENCFF616WSS.fastq.gz rep2-R2.fastq.gz
$ mv ENCFF002EFT.fastq.gz ctl1-R1.fastq.gz; mv ENCFF002EFQ.fastq.gz ctl1-R2.fastq.gz  
$ mv ENCFF002EFU.fastq.gz ctl2-R1.fastq.gz; mv ENCFF002EFS.fastq.gz ctl2-R2.fastq.gz  
\end{verbatim}

Make soft links for the single-ended (SE) test case.

\begin{verbatim}
$ ln -s rep1-R1.fastq.gz rep1.fastq.gz; ln -s rep2-R1.fastq.gz rep2.fastq.gz
$ ln -s ctl1-R1.fastq.gz ctl1.fastq.gz; ln -s ctl2-R1.fastq.gz ctl2.fastq.gz
\end{verbatim}

Install genome database for a species \verb|hg38| on the genome database directory.

\begin{verbatim}
$ mkdir genome_data
$ bash chipseq_pipeline/install_genome_data.sh hg38 genome_data
\end{verbatim}

\subsection*{Procedure using BigDataScript (recommended)}

All subtasks described in the next subsection are run in a correct order with one command line. Specify your own output directory with \verb|‘-out_dir’|. The flag \verb|‘-use_pooled_ctl’| forces the pipeline to use pooled control replicates for relaxed peaks called with a peak caller SPP. There are more parameters to customize the pipeline. Details of all parameters are described in Box Y.

\begin{itemize}
\item For single-ended replicates,
\begin{verbatim}
$ NTH=8 # use 8 threads for the test case SE
$ FASTQ1=rep1.fastq.gz; FASTQ2=rep2.fastq.gz
$ CTL1=ctl1.fastq.gz; CTL2=ctl2.fastq.gz
$ python chipseq_pipeline/chipseq.py -se -species hg38 
  -fastq1 $FASTQ1 -fastq2 $FASTQ2 
  -ctl_fastq1 $CTL1 -ctl_fastq2 $CTL2 
  -out_dir ENCSR000DYI/SE -nth $NTH -use_pooled_ctl
\end{verbatim}

Or you can also run a pipeline with a configuration JSON file that has all the parameters (key names are the same as those in command line arguments) specified within. See Box V for the contents of such JSON file (\verb|ENCSR000DYI_SE.json|).

\begin{verbatim}
$ python chipseq_pipeline/chipseq.py ENCSR000DYI_SE.json
\end{verbatim}

\item For paired-end replicates,
\begin{verbatim}
$ NTH=16 # use 16 threads for the test case PE
$ FASTQ1_1=rep1-R1.fastq.gz; FASTQ1_2=rep1-R2.fastq.gz
$ FASTQ2_1=rep2-R1.fastq.gz; FASTQ2_2=rep2-R2.fastq.gz
$ CTL1_1=ctl1-R1.fastq.gz; CTL1_2=ctl1-R2.fastq.gz
$ CTL2_1=ctl2-R1.fastq.gz; CTL2_2=ctl2-R2.fastq.gz
$ python chipseq_pipeline/chipseq.py -pe -species hg38 
  -fastq1_1 $FASTQ1_1 -fastq1_2 $FASTQ1_2 
  -fastq2_1 $FASTQ2_1 -fastq2_2 $FASTQ2_2 
  -ctl_fastq1_1 $CTL1_1 -ctl_fastq1_2 $CTL1_2 
  -ctl_fastq2_1 $CTL2_1 -ctl_fastq2_2 $CTL2_2 
  -out_dir ENCSR936XTK/PE -nth $NTH -use_pooled_ctl
\end{verbatim}

Or you can also run a pipeline with a configuration JSON file that has all the parameters specified within. See Box V for the contents of such JSON file (\verb|ENCSR936XTK_PE.json|).

\begin{verbatim}
$ python chipseq_pipeline/chipseq.py ENCSR936XTK_PE.json
\end{verbatim}
\end{itemize}

\subsection*{Procedure using bash shell commands}

\begin{enumerate}
\item \textit{Create a directory structure for outputs}. Create directories to store mapped/filtered BAMs, TAGALIGNs and quality metrics. For convenience, the directory structure for this demonstration is simpler than that of the pipeline (\burl{https://github.com/kundajelab/chipseq_pipeline#output-directory-structure-and-file-naming}).

\begin{verbatim}
$ mkdir align; mkdir qc
\end{verbatim}

Create directories to store called peaks.

\begin{verbatim}
$ mkdir peak; mkdir peak/reps; mkdir peak/pooled; 
$ mkdir peak/self_pseudo_reps; mkdir peak/pooled_pseudo_reps
\end{verbatim}

Create directories to store signal tracks

\begin{verbatim}
$ mkdir signal; mkdir signal/reps; mkdir signal/pooled;
\end{verbatim}

Create directories to store IDR peaks.

\begin{verbatim}
$ mkdir idr; mkdir idr/true_reps; mkdir idr/pooled; 
$ mkdir idr/self_pseudo_reps; mkdir idr/pooled_pseudo_reps
\end{verbatim}

\item \textit{Load Conda environment for the pipeline.}

Please make sure that \verb|‘(aquas_chipseq)’| is prepended to the bash prompt after loading Conda environments.

\begin{verbatim}
$ source activate aquas_chipseq
\end{verbatim}

\textbf{Map raw reads for AQUAS}: TIMING $\sim$6 h (when fully parallelized)

\item \textit{Map raw reads in FASTQs to BAMs}. Use two threads per FASTQ. Set an environment variable for bwa index.

\begin{verbatim}
$ NTH=2; BWA_IDX=genome_data/hg38/bwa_index/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta
\end{verbatim}

Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
$ bwa aln -q 5 -l 32 -k 2 -t $NTH $BWA_IDX $REP.fastq.gz > align/$REP.sai
$ bwa samse $BWA_IDX align/$REP.sai $REP.fastq.gz 
  | samtools view -Su - | samtools sort - align/$REP
$ samtools index align/$REP.bam
$ samtools flagstat align/$REP.bam > qc/$REP.flagstat
\end{verbatim}

For paired-end replicates:

\begin{verbatim}
$ bwa aln -q 5 -l 32 -k 2 -t $NTH $BWA_IDX $REP-R1.fastq.gz > align/$REP-R1.sai
$ bwa aln -q 5 -l 32 -k 2 -t $NTH $BWA_IDX $REP-R2.fastq.gz > align/$REP-R2.sai
$ bwa sampe $BWA_IDX align/$REP-R1.sai align/$REP-R2.sai $REP-R1.fastq.gz $REP-R2.fastq.gz 
  | gzip -nc > align/$REP.sam.gz
\end{verbatim}

Remove read pairs with bad CIGAR strings and sort by position for the paired-end case.

\begin{verbatim}
$ zcat align/$REP.sam.gz | awk 'BEGIN {FS="\t" ; OFS="\t"} ! /^@/ && $6!="*" 
  { cigar=$6; gsub("[0-9]+D","",cigar); n = split(cigar,vals,"[A-Z]"); s = 0;
   for (i=1;i<=n;i++) s=s+vals[i]; seqlen=length($10) ; if (s!=seqlen) 
   print $1"\t"; }' | sort | uniq > align/$REP.badCigar
$ if [ $(cat align/$REP.badCigar | wc -l) -gt 0 ]; then zcat align/$REP.sam.gz | 
  grep -v -F -f align/$REP.badCigar | samtools view -Su - | samtools sort - 
  align/$REP; else samtools view -Su align/$REP.sam.gz | samtools sort - align/$REP; fi
$ samtools flagstat align/$REP.bam > qc/$REP.flagstat
$ samtools index align/$REP.bam
\end{verbatim}

Repeat for other IP/control replicates (\verb|REP=rep2|, \verb|REP=ctl1| and \verb|REP=ctl2|).

Remove temporary files.

\begin{verbatim}
$ rm -f align/*.badCigar align/*.sai align/$REP.badCigar align/*.sam.gz
\end{verbatim}

\textbf{Prepare sequences for AQUAS}: TIMING $\sim$2 hours

\item \textit{Filter multiply mapped reads}. If the BAM file contain multiply aligned reads, you can filter them (i.e., those with MAPQ$<$30) using \verb|samtools -q|. We use \verb|sambamba| instead of \verb|samtools| for some part of the pipeline to speed up. \verb|sambamba| is a multi-threaded version of \verb|samtools|, which is faster and more stable. Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
$ samtools view -F 1804 -q 30 -u align/$REP.bam | 
  sambamba sort /dev/stdin -o align/$REP.filt.bam
$ export _JAVA_OPTIONS="-Xms256M -Xmx$12G -XX:ParallelGCThreads=1"
\end{verbatim}

Mark duplicate reads for the single-ended case.

\begin{verbatim}
$ java -jar ${PICARDROOT}/MarkDuplicates.jar INPUT="align/$REP.filt.bam" 
  OUTPUT="align/$REP.dupmark.bam" METRICS_FILE="qc/$REP.dup" 
  VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true REMOVE_DUPLICATES=false
\end{verbatim}

Remove duplicate reads for the single-ended case.

\begin{verbatim}
$ samtools view -F 1804 -b align/$REP.dupmark.bam > align/$REP.nodup.bam
$ sambamba index align/$REP.nodup.bam
$ sambamba flagstat align/$REP.nodup.bam > qc/$REP.nodup.flagstat
\end{verbatim}

Compute library complexity for the single-ended case.

\begin{verbatim}
$ bedtools bamtobed -i align/$REP.dupmark.bam | awk 'BEGIN{OFS="\t"}{print $1,$2,$3,$6}' 
  | grep -v 'chrM' | sort | uniq -c | awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} 
  ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} END{m1_m2=-1.0; if(m2>0) m1_m2=m1/m2; 
  printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1_m2}' > qc/$REP.pbc
\end{verbatim}

For paired-end replicates,

\begin{verbatim}
$ samtools view -F 1804 -f 2 -q 30 -u align/$REP.bam 
  | sambamba sort -n /dev/stdin -o align/$REP.tmp.bam
$ samtools fixmate -r align/$REP.tmp.bam align/$REP.fixmate.bam
$ samtools view -F 1804 -f 2 -u align/$REP.fixmate.bam 
  | sambamba sort /dev/stdin -o align/$REP.filt.bam
\end{verbatim}

Mark duplicate reads for the paired-end case.

\begin{verbatim}
$ export _JAVA_OPTIONS="-Xms256M -Xmx$12G -XX:ParallelGCThreads=1"
$ java -jar ${PICARDROOT}/MarkDuplicates.jar INPUT="align/$REP.filt.bam" 
  OUTPUT="align/$REP.dupmark.bam" METRICS_FILE="qc/$REP.dup" 
  VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true REMOVE_DUPLICATES=false
\end{verbatim}

Remove duplicate reads for the paired-end case.

\begin{verbatim}
$ samtools view -F 1804 -f 2 -b align/$REP.dupmark.bam > align/$REP.nodup.bam
$ sambamba index align/$REP.nodup.bam
$ sambamba flagstat align/$REP.nodup.bam > qc/$REP.nodup.flagstat
$ sambamba sort -n align/$REP.dupmark.bam -o align/$REP.dupmark.tmp.bam
\end{verbatim}

Compute library complexity for the paired-end case.

\begin{verbatim}
$ bedtools bamtobed -bedpe -i align/$REP.dupmark.tmp.bam 
  | awk 'BEGIN{OFS="\t"}{print $1,$2,$4,$6,$9,$10}' | grep -v 'chrM' | sort 
  | uniq -c | awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} ($1==2){m2=m2+1} 
  {m0=m0+1} {mt=mt+$1} END{m1_m2=-1.0; if(m2>0) m1_m2=m1/m2; 
  printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1_m2}' > qc/$REP.pbc
\end{verbatim}

Repeat for other IP/control replicates (\verb|REP=rep2|, \verb|REP=ctl1| and \verb|REP=ctl2|).

Remove temporary files.

\begin{verbatim}
$ rm -f align/*.tmp.bam* align/*.fixmate.bam* align/*.dupmark.bam* align/*.filt.bam*
\end{verbatim}

\item \textit{Transfer file format}. For the IDR rescue strategy later, you need to shuffle and subsample version of the alignment files. This is easier to implement on tagAlign text files using the UNIX \verb|'shuf'| command. Thus, it is recommended to transfer a BAM format to a tagAlign format. Start with IP replicate 1. 

\begin{verbatim}
$ REP=rep1
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
$ bedtools bamtobed -i align/$REP.nodup.bam 
  | awk 'BEGIN{OFS="\t"}{$4="N";$5="1000";print $0}' 
  | gzip -nc > align/$REP.tagAlign.gz
\end{verbatim}

For paired-end replicates:

\begin{verbatim}
$ sambamba sort -n align/$REP.nodup.bam -o align/$REP.nodup.nmsrt.bam
$ bedtools bamtobed -bedpe -mate1 -i align/$REP.nodup.nmsrt.bam 
  | gzip -nc > align/$REP.bedpe
$ zcat align/$REP.bedpe | awk 'BEGIN{OFS="\t"}
  {printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}' 
  | gzip -nc > align/$REP.tagAlign.gz
\end{verbatim}

Repeat for other IP/control replicates (\verb|REP=rep2|, \verb|REP=ctl1| and \verb|REP=ctl2|).

Remove temporary files.

\begin{verbatim}
$ rm -f align/*.nodup.nmsrt.bam
\end{verbatim}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box X $|$ Software requirements (installed by Conda package manager) in AQUAS}}
\centerline{}
\centerline{}
\textbf{Aligner}

\begin{itemize}
\item \verb|bwa| $\geq$ 0.7.13
\textbf{Samtools and bedtools}
\end{itemize}

\begin{itemize}
\item \verb|samtools| ($\geq$ 1.2), \verb|htslib|, \verb|sambamba| ($\geq$ 0.6.5) and \verb|bedtools| ($\geq$ 2.26.0)
\end{itemize}

\textbf{UCSC tools}

\begin{itemize}
\item \verb|ucsc-fetchchromsizes|, \verb|ucsc-wigtobigwig|, \verb|ucsc-bedgraphtobigwig|, \verb|ucsc-bigwiginfo|, \verb|ucsc-bedclip|, \verb|ucsc-bedtobigbed| and \verb|ucsc-twobittofa|
\end{itemize}

\textbf{PICARD tools}

\begin{itemize}
\item \verb|picard| (= 1.126)
\end{itemize}

\textbf{Python and its packages}

\begin{itemize}
\item Python2 and packages: \verb|python2| ($\geq$ 2.7), \verb|numpy| (= 1.10.2), \verb|matplotlib| (= 1.5.1), \verb|macs2| (= 2.1.0) and \verb|pyfaidx| (= 0.4.7.1).
\item Python3 and packages: \verb|python3| ($\geq$ 3.5) and \verb|idr| ($\geq$ 2.0.3).
\end{itemize}

\textbf{R and its packages}
\begin{itemize}
\item \verb|R| ($\geq$ 3.2.2), \verb|r-snow|, \verb|r-snowfall|, \verb|r-bitops|, \verb|r-catools|, \verb|r-spp| (= 1.14) and \verb|bioconductor-rsamtools|.
\end{itemize}

\textbf{Others}

\begin{itemize}
\item \verb|boost| (= 1.57.0), \verb|openssl| (= 1.0.2g-0), \verb|graphviz| ($\geq$ 2.38.0), \verb|libtool|, \verb|ghostscript|, \verb|pigz| and \verb|zlib|
\end{itemize}
\par}
\end{mdframed}
\centerline{}

\textbf{Cross-correlation measures}: TIMING $\sim$30 min

\item Subsample tagAlign and BEDPE files to have 15M reads and compute cross-correlation measures (NSC and RSC). Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
$ zcat align/$REP.tagAlign.gz | grep -v "chrM" 
  | shuf -n 15000000 --random-source=align/$REP.tagAlign.gz 
  | gzip -nc > align/$REP.15M.tagAlign.gz
\end{verbatim}

For paired-end replicates:

\begin{verbatim}
$ zcat align/$REP.bedpe | grep -v "chrM" 
  | shuf -n 15000000 --random-source=align/$REP.bedpe 
  | awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"N","1000",$9}' 
  | gzip -nc > align/$REP.15M.tagAlign.gz    
\end{verbatim}

Run R script for the cross-correlation analysis. For the details of the options and output format, please see Box 3.

\begin{verbatim}
$ Rscript $(which run_spp.R) -rf -c=align/$REP.15M.tagAlign.gz 
  -filtchr=chrM -savp=qc/$REP.cc.pdf -out=qc/$REP.cc
$ sed -r 's/,[^\t]+//g' qc/$REP.cc > qc/$REP.cc.tmp
$ mv qc/$REP.cc.tmp qc/$REP.cc
\end{verbatim}

\setcounter{table}{0}
\renewcommand{\tablename}{}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{ 
% \small \textbf{\underline{Box Y $|$ Parameters for AQUAS}}
% \centerline{}
% \centerline{}
\begin{small}
\begin{center}
\begin{longtable}{|P{3.7cm}|P{4cm}|P{8cm}|}
\caption*{\textbf{\underline{Box Y $|$ Parameters for AQUAS}}}\\
\hline
category & parameter & description \\
\hline
\endfirsthead
\multicolumn{3}{c}%
{\tablename\ \thetable\ -- \textit{Continued from previous page}} \\
\hline
category & parameter & description \\
\hline
\endhead
\hline \multicolumn{3}{r}{\textit{Continued on next page}} \\
\endfoot
% \hline
\endlastfoot
\multirow{6}{*}{general} & screen \verb|<str>| & Create a screen with name \verb|[SCREEN].BDS| and attach BDS pipeline to it. \verb|STDOUT/STDERR| will be redirected to \verb|[SCREEN].log| on the working directory \\
 & type \verb|<str>| & Type of ChIP-Seq pipeline. {\verb|TF|, \verb|histone|} \\
 & \verb|final_stage <str>\verb| & Final stage for pipeline. {\verb|bam|, \verb|filt_bam|, \verb|tag|, \verb|xcor|, \verb|peak| or \verb|idr|} \\
 & \verb|out_dir <str>\verb| & Output directory. \\
 & \verb|title <str>| & Prefix for HTML report and outputs without given prefix. \\
 & \verb|dry_run| & Dry-run.\\
\hline
\multirow{4}{*}{input-endedness} & \verb|-se| & Set endedness of all inputs as single-ended.\\
 & \verb|-pe| & Set endedness of all inputs as paired-end.\\
 & \verb|-se[REP_ID]|, \verb|-ctl_se[CTL_ID]| & Set endedness of each IP replicate/control as single-ended.\\
 & \verb|-pe[REP_ID]|, \verb|-ctl_pe[CTL_ID]| & Set edndedness of each IP replicate/control as paired-end.\\
\hline
\multirow{10}{*}{input-files} & \verb|-fastq[REP_ID]| \verb|<path>| & Single-ended FASTQs for each IP replicate.\\
 & \verb|-ctl_fastq[CTL_ID]| \verb|<path>| & Single-ended FASTQs for each contol.\\
 & \verb|-fastq[REP_ID]_[PAIR]| \verb|<path>| & Paired end FASTQs for each IP replicate.\\
 & \verb|-ctl_fastq[CTL_ID]_[PAIR]| \verb|<path>| & Paired end FASTQs for each contol.\\
 & \verb|-bam[REP_ID]| \verb|<path>| & Raw BAM for each IP replicate.\\
 & \verb|-ctl_bam[CTL_ID]| \verb|<path>| & Raw BAM for each contol.\\
 & \verb|-bam[REP_ID]| \verb|<path>| & Filtered BAM for each IP replicate.\\
 & \verb|-ctl_bam[CTL_ID]| \verb|<path>| & Filtered BAM for each contol.\\
 & \verb|-tag[REP_ID]| \verb|<path>| & \verb|TAGALIGN| for each IP replicate.\\
 & \verb|-ctl_tag[CTL_ID]| \verb|<path>| & \verb|TAGALIGN| for each contol.\\
\hline
\multirow{7}{*}{species} & species \verb|<str>| & Species.\\
 & \verb|species_file| \verb|<path>| & Species file path.\\
 & \verb|species_browser| \verb|<str>| & Species name in WashU genome browser.\\
 & \verb|ref_fa| \verb|<path>| & Reference genome sequence fasta.\\
 & \verb|chrsz| \verb|<path>| & Chromosome sizes file path (use \verb|fetchChromSizes| from UCSC tools).\\
 & \verb|blacklist| \verb|<path>| & Blacklist bed.\\
 & \verb|gensz| \verb|<path>| & Genome size; \verb|hs| for human, \verb|mm| for mouse. Or sum of numbers in the 2nd column in \verb|chrom.sizes| file.\\
\hline
\multirow{4}{*}{cluster} & system \verb|<str>| & Force to use a system. {\verb|local|, \verb|sge|, \verb|slurm|}\\
 & \verb|nice\verb| \verb|<int>| & Set process priority for all tasks; \verb|-20| (highest) \verb|~19| (lowest).\\
 & \verb|retrial\verb| \verb|<int>| & Number of retrial for failed tasks.\\
 & \verb|q| & Submit tasks to a specified cluster queue.\\
\hline
\multirow{14}{*}{resource} & \verb|nth\verb| \verb|<int>| & Maximum number of threads for a pipeline.\\
 & \verb|no_par| & Serialize all tasks (individual tasks can still use multiple threads up to \verb|'-nth'|).\\
 & \verb|wt\verb| \verb|<str>| & Walltime for all single-threaded tasks (example: \verb|8:10:00|, \verb|3h|, \verb|3600|).\\
 & \verb|memory| \verb|<str>| & Maximum memory for all single-threaded tasks (equivalent to \verb|'-mem'|, example: \verb|4.5G|, \verb|1024M|).\\
 & \verb|unlimited_mem_wt| & Use unlimited max. memory and walltime.\\
 & \verb|wt_dedup| \verb|<str>| & Walltime for post-alignment filtering.\\
 & \verb|mem_dedup| \verb|<str>| & Max. memory for post-alignment filtering.\\
 & \verb|mem_shuf| \verb|<str>| & Max. memory for UNIX shuf.\\
 & \verb|wt_bwa| \verb|<str>| & Walltime for \verb|bwa|.\\
 & \verb|mem_bwa| \verb|<str>| & Max. memory for \verb|bwa|.\\
 & \verb|wt_macs2| \verb|<str>| & Walltime for \verb|MACS2|.\\
 & \verb|mem_macs2| \verb|<str>| & Max. memory for \verb|MACS2|.\\
 & \verb|wt_spp| \verb|<str>| & Walltime for \verb|spp|. \\
 & \verb|mem_spp\verb| \verb|<str>| & Max. memory for \verb|spp|. \\
alignment & aligner \verb|<str>| & Aligner to map raw reads in FASTQs. \verb|{bwa}| \\
\hline
\multirow{2}{*}{alignment-bwa} & \verb|param_bwa_aln| \verb|<str>| & Parameters for \verb|bwa aln|. \\
 & \verb|bwa_idx| \verb|<path>| & BWA index (full path prefix of \verb|*.bwt| file). \\
\hline
\multirow{5}{*}{alignment-filtering} & \verb|dup_marker| \verb|<str>| & Dup marker for filtering mapped reads in BAMs. {\verb|picard|, \verb|sambamba|} \\
 & \verb|anon_filt_bam| & Generate an annomymized filtered BAM. This will not affect tasks downsteam. \\
 & \verb|mapq_thresh| \verb|<int>| & Threshold for low MAPQ reads removal. \\
 & \verb|rm_chr_from_tag| \verb|<str>| & Perl style reg-ex to exclude reads from tag-aligns. (example: \verb|'other|$|$\verb|ribo|$|$\verb|mito|$|$\verb|_'|). \\
 & \verb|no_dup_removal| & No dupe removal when filtering raw BAM. \\
\hline
\multirow{2}{*}{alignment-subsampling} & \verb|subsample_chip| \verb|<str>| & Number of reads to subsample IP replicate. Subsampled tagAlign will be used for analysis (peak and IDR) downstream (0: no subsampling). \\
 & \verb|subsample_ctl| \verb|<str>| & Number of reads to subsample control if non-zero (recommended: 40M or lower). \\
\hline
\multirow{4}{*}{cross correlation analysis} & \verb|no_xcor| & Disable cross-correlation analaysis. \\
 & \verb|subsample_xcor| \verb|<str>| & Number of reads to be subsampled for cross correlation analysis only. \\
 & \verb|speak_xcor| \verb|<str>| & Set user-defined cross-correlation peak strandshift (\verb|-speak=| in \verb|run_spp.R|). Use \verb|-1| to disable. \\
 & \verb|extra_param_xcor| \verb|<str>| & Set extra parameters for \verb|run_spp.R|. \\
\hline
\multirow{5}{*}{callpeak} & \verb|peak_caller| \verb|<str>| & Peak caller for IDR analysis. \verb|spp| for TF ChIP-seq and \verb|macs2| for Histone ChIP-seq. {\verb|spp|, \verb|macs2|} \\
 & \verb|ctl_depth_ratio| \verb|<float>| & Cut-off ratio of two control tagAligns for pooling. \\
 & \verb|use_pooled_ctl| & Force to use pooled control (ignoring criteria to choose control for each IP replicate). \\
 & \verb|true_rep| & Call peaks on true replicates only. \\
 & \verb|no_pseudo_rep| & Do not call peaks on self pseudoreplicates. \\
\hline
\multirow{4}{*}{callpeak-spp} & \verb|npeak_spp| \verb|<int>| & Threshold on maximum number of peaks (\verb|-npeak| in \verb|run_spp.R|). \\
 & \verb|max_ppsize_spp| \verb|<str>| & Set \verb|R| stack size (\verb|R| parameter \verb|--max-ppsize=|; between \verb|5000| and \verb|5000000|) for SPP. \\
 & \verb|speak_spp| \verb|<int>| & Set user-defined cross-correlation peak strand shift (\verb|-speak=| in \verb|run_spp.R|). Use \verb|-1| to get from upstream cross-correlation analysis. \\
 & \verb|extra_param_spp| \verb|<str>| & Set extra parameters for \verb|run_spp.R|. \\ 
\hline
\multirow{5}{*}{callpeak-macs2} & \verb|pval_thresh_macs2| & \verb|--pvalue| for \verb|macs2| \verb|callpeak| (\burl{https://github.com/taoliu/MACS#-p--pvalue}). \\
 & \verb|keep_dup_macs2| \verb|<str>| & \verb|--keep-dup| for \verb|macs2| \verb|callpeak| (\burl{https://github.com/taoliu/MACS#--keep-dup}) {all, any positive interger} \\
 & \verb|extsize_macs2\verb| \verb|<int>| & \verb|--extsize| for \verb|macs2| \verb|callpeak| (\burl{https://github.com/taoliu/MACS#--extsize}). Use \verb|-1| to get from upstream cross-correlation analysis. \\
 & \verb|shift_macs2| \verb|<int>| & \verb|--shift| for \verb|macs2| \verb|callpeak| (\burl{https://github.com/taoliu/MACS#--shift}). \\
 & \verb|extra_param_macs2| \verb|<str>| & Set extra parameters for \verb|macs2| \verb|callpeak|. \\
\hline
\multirow{3}{*}{callpeak-idr} & \verb|idr_suffix| & Append IDR threshold to IDR output directory. \\
 & \verb|idr_rank| \verb|<str>| & Scoring column in narrow peak files for IDR. If not specified, \verb|signal.value| for SPP peaks (TF) and \verb|p.value| for MACS2 peaks (histone) are used. {\verb|signal.value|, \verb|p.value|, \verb|q.value|} \\
 & \verb|idr_thresh| \verb|<float>| & IDR threshold : $-log_{10}$(score). \\
\verb|callpeak-naive| \verb|overlap| & \verb|nonamecheck| & \verb|bedtools| \verb|intersect| \verb|-nonamecheck| (\verb|bedtools|$\geq$2.24.0, use this if you get \verb|bedtools| \verb|intersect| naming convenction warnings/errors). \\
\hline
signal track & \verb|sig_trk_for_pooled_rep_only| & Generate MACS2 signal tracks for pooled replicate only. \\ 
\hline
\multirow{2}{*}{visualization} & \verb|url_base| \verb|<url>| & Equivalent URL address for your output directory (\verb|-out_dir|). \\
 & \verb|viz_genome_coord| \verb|<str>| & WashU genome browser genome coordinate (e.g. \verb|chr7:27117661-27153380|). \\
\hline
\multirow{3}{*}{shell environment-conda} & \verb|conda_env| \verb|<str>| & Conda environment name for all softwares including Python2 (for all others). \\
 & \verb|conda_env_py3| \verb|<str>| & Conda environment name for Python3 (for IDR only). \\
 & \verb|conda_bin_dir| \verb|<path>| & Conda bin directory. Useful when you have multiple Conda's in your \verb|$PATH| and want to specify one. \\
\hline
\multirow{3}{*}{shell environment-modules} & \verb|mod| \verb|<str>| & Modules separated by \verb|;| (example: \verb|'bowtie/2.2.4;| \verb|bwa/0.7.7;| \verb|picard-tools/1.92'|). \\
 & \verb|shcmd| \verb|<str>| & shell commands separated by \verb|;\verb|. Shell variables must be written as \verb|${VAR}| not as \verb|$VAR| (example: \verb|'export| \verb|PATH=${PATH}:/usr/test;| \verb|VAR=test'|). \\
 & \verb|addpath| \verb|<str>| & Path separated by \verb|;| or \verb|:| to be PREPENDED to \verb|$PATH| (example: \verb|'/bin/test:${HOME}/utils'|). \\
\hline
\label{TableS}
\end{longtable}
\end{center}
\end{small}
}
\end{mdframed}
\centerline{}

\centerline{}
\begin{mdframed}[roundcorner=2pt,userdefinedwidth=18cm,align=center,linewidth=0pt,backgroundcolor=gray!18]
{\small \textbf{\underline{Box Z $|$ Criteria to choose control in AQUAS}}
\centerline{}
It is important to choose an appropriate control replicate for each IP replicate when calling peaks (using peak-caller SPP) and generating signal tracks (using peak-caller MACS2). There are five criteria.

\begin{itemize}
\item If the flag \verb|‘-use_pooled_ctl’| is activated in the pipeline, choose pooled control replicate for all IP replicates. Such flag is activated for the demonstration test cases.
\item If there is only one control replicate provided, choose it for all IP replicates.
\item If number of reads in any control replicate differ by a factor of the control depth ratio (1.2 by default), choose pooled control replicate for all IP replicates.
\item If there are fewer reads in control $i$ than IP replicate $i$, choose pooled control replicate for IP replicate $i$.
\item Otherwise, choose control replicate $i$ for IP replicate $i$. 
\end{itemize}
\par}
\end{mdframed}
\centerline{}

Repeat for IP replicate 2 (\verb|REP=rep2|).

You can see the results below and in \hl{Figure XXXXXXX}. Please note that QualityTag values for controls are -1, while 1 for IPs.

\begin{verbatim}
$ cat qc/*.cc
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
rep1.15M.tagAlign.gz 15000000 220 0.220570692618442 105 0.2058052 1500 0.1689261 1.305723 1.400376 1
rep2.15M.tagAlign.gz 15000000 190 0.197236457317112 105 0.1905799 1500 0.1726345 1.142509 1.370933 1
\end{verbatim}

For paired-end replicates:

\begin{verbatim}
rep1.15M.tagAlign.gz 15000000 225 0.255919454649913 105 0.2114665 1500 0.1503787 1.701833 1.727691 2
rep2.15M.tagAlign.gz 15000000 220 0.234245076242424 105 0.2077840 1500 0.1676446 1.397272 1.659227 2
\end{verbatim}

\textbf{CRITICAL STEP} Datasets with low cross-correlation coefficients do not necessarily represent experimental failures. The biology of the target factor under study must be considered since ChIP-seq datasets targeting factors that bind sparsely due to the inherent biology of the factor or cell-type specific activity can result in low cross-correlation coefficients (see Figure 7b).

\item \hl{XXXXX Compute JD metric XXXXX}

\textbf{Prepare input files for the IDR framework}: TIMING $<$ 20 min

For self-consistency analysis, IDR analyses will be performed over each replicate, merged replicates (pooled data) and pseudo-replicates for them (for the details, see Experimental design).

\item Subsample to make pseudo-replicates. For self-consistency analysis, randomly split each IP replicate and pooled IP data. Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1
\end{verbatim}

For single-ended replicates:

\begin{verbatim}
$ nlines=$( zcat align/$REP.tagAlign.gz | wc -l )
$ nlines=$(( (nlines + 1) / 2 ))
$ zcat align/$REP.tagAlign.gz 
  | shuf --random-source=align/$REP.tagAlign.gz 
  | split -d -l $((nlines)) - align/$REP.bam_pr1.
$ gzip -nc align/$REP.bam_pr1.00 > align/$REP.pr1.tagAlign.gz
$ gzip -nc align/$REP.bam_pr1.01 > align/$REP.pr2.tagAlign.gz
$ rm -f align/$REP.bam_pr1.00 align/$REP.bam_pr1.01
\end{verbatim}

For paired-end replicates:

\begin{verbatim}
$ nlines=$( zcat align/$REP.bedpe | wc -l )
$ nlines=$(( (nlines + 1) / 2 ))
$ zcat align/$REP.bedpe | shuf --random-source=align/$REP.bedpe 
  | split -d -l $((nlines)) - align/$REP.bam_pr1.
$ awk 'BEGIN{OFS="\t"}
  {printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}' 
  "align/$REP.bam_pr1.00" | gzip -nc > align/$REP.pr1.tagAlign.gz
$ awk 'BEGIN{OFS="\t"}
  {printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}' 
  "align/$REP.bam_pr1.01" | gzip -nc > align/$REP.pr2.tagAlign.gz
$ rm -f align/$REP.bam_pr1.00 align/$REP.bam_pr1.01
\end{verbatim}

Repeat for IP replicate 2 (\verb|REP=rep2|).

\item Merge replicates. pool the data for IP samples and for control samples.

For IP samples:

\begin{verbatim}
$ zcat align/rep1.tagAlign.gz align/rep2.tagAlign.gz | gzip -nc > align/rep0.tagAlign.gz
$ zcat align/rep1.pr1.tagAlign.gz align/rep2.pr1.tagAlign.gz 
                                                 | gzip -nc > align/rep0.pr1.tagAlign.gz
$ zcat align/rep1.pr2.tagAlign.gz align/rep2.pr2.tagAlign.gz 
                                                 | gzip -nc > align/rep0.pr2.tagAlign.gz
\end{verbatim}

For control samples:

\begin{verbatim}
$ zcat align/ctl1.tagAlign.gz align/ctl2.tagAlign.gz | gzip -nc > align/ctl0.tagAlign.gz
\end{verbatim}

\textbf{IDR analysis}: TIMING $<$ 1 day

\item Relaxed peak calls. each IP replicate, pooled IP and its pseudo-replicates, perform peak calling with relaxed thresholds to call 300,000 peaks using SPP peak caller. ($<$ 1 day). Run R script for the IDR analysis. For the details of the options and output format, please see Box 3. Read fragment length for each of the original replicates from the 3rd column of cross-correlation quality analysis log. Use mean fragment length for call peaks for pooled replicates and pooled pseudo replicates.

\begin{verbatim}
$ fraglen_rep1=$(cat qc/rep1.cc | awk '{print $3}')
$ fraglen_rep2=$(cat qc/rep2.cc | awk '{print $3}')
$ fraglen_mean=$(( (fraglen_rep1+fraglen_rep2)/2 ))
\end{verbatim}

For convenience, we call peaks with a pooled control in this demonstration. However in the pipeline, an appropriate control is chosen according to the criteria described in Box Z. Use a flag \verb|‘-use_pooled_ctl’| to force the pipeline to use pooled control for calling peaks. Use 2 threads for each case. Use pooled control replicate for all cases.

\begin{verbatim}
$ NTH=2; CTL=ctl0
\end{verbatim}

Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1; FRAGLEN=$fraglen_rep1; OUTDIR=peak/reps
$ Rscript $(which run_spp.R) -c=align/$REP.tagAlign.gz -p=$NTH 
  -i=align/$CTL.tagAlign.gz -npeak=300000 -odir=$OUTDIR -speak=$FRAGLEN -savr -savp -rf
$ rpeakfile_raw=$OUTDIR/$REP.tagAlign_VS_$CTL.tagAlign.regionPeak.gz
$ rpeakfile=$OUTDIR/$REP.regionPeak.gz
$ rpeakfile_filt=$OUTDIR/$REP.filt.regionPeak.gz
\end{verbatim}

Remove peaks with illegal genome coordinate.

\begin{verbatim}
$ zcat $rpeakfile_raw | awk 'BEGIN{OFS="\t"}{ if ($2<0) $2=0; 
  print $1,int($2),int($3),$4,$5,$6,$7,$8,$9,$10;}' | gzip -f -nc > $rpeakfile
$ mv $OUTDIR/$REP.tagAlign.pdf qc/$REP.relaxed_peak.pdf
\end{verbatim}

Filter out raw relaxed peaks in blacklisted region.

\begin{verbatim}
$ bedtools intersect -v -a <(zcat -f $rpeakfile) -b 
  <(zcat -f genome_data/hg38/hg38.blacklist.bed.gz) | awk 'BEGIN{OFS="\t"} 
  {if ($5>1000) $5=1000; print $0}' | grep -P 'chr[\dXY]+[ \t]' | gzip -nc > $rpeakfile_filt
\end{verbatim}

Remove temporary files.

\begin{verbatim}
$ rm -f $rpeakfile_raw
\end{verbatim}

Repeat for:

IP replicate 2 (\verb|REP=rep2|; \verb|FRAGLEN=$fraglen_rep2|; \verb|OUTDIR=peak/reps|), \\
pooled IP replicate (\verb|REP=rep0|; \verb|FRAGLEN=$fraglen_mean|; \verb|OUTDIR=peak/pooled|), \\
the first pseudo replicate from IP replicate 1 (\verb|REP=rep1.pr1|; \verb|FRAGLEN=$fraglen_rep1|; \\ \verb|OUTDIR=peak/self_pseudo_reps|), \\
the second pseudo replicate from IP replicate 1 (\verb|REP=rep1.pr2|; \verb|FRAGLEN=$fraglen_rep1|; \\ \verb|OUTDIR=peak/self_pseudo_reps|), \\
the first pseudo replicate from IP replicate 2 (\verb|REP=rep2.pr1|; \verb|FRAGLEN=$fraglen_rep2|; \\ \verb|OUTDIR=peak/self_pseudo_reps|), \\
the second pseudo replicate from IP replicate 2 (\verb|REP=rep2.pr2|; \verb|FRAGLEN=$fraglen_rep2|; \\ \verb|OUTDIR=peak/self_pseudo_reps|), \\
the first pooled pseudo replicate (\verb|REP=rep0.pr1|; \verb|FRAGLEN=$fraglen_mean|; \\ \verb|OUTDIR=peak/pooled_pseudo_reps|), \\
and the second pooled pseudo replicate (\verb|REP=rep0.pr2|; \verb|FRAGLEN=$fraglen_mean|; \\ \verb|OUTDIR=peak/pooled_pseudo_reps|).

\item Signal track generation for each IP replicate and pooled IP. Perform peak calling with p-value threshold of 0.01 using MACS2 peak caller. Use pooled control for all IP replicate.

\begin{verbatim}
$ CTL=ctl0
\end{verbatim}

Start with IP replicate 1.

\begin{verbatim}
$ REP=rep1; FRAGLEN=$fraglen_rep1; OUTDIR=signal/reps
\end{verbatim}

Generate preliminary signal tracks.

\begin{verbatim}
$ export LC_COLLATE=C
$ macs2 callpeak -t align/$REP.tagAlign.gz align/$CTL.tagAlign.gz -f BED -n $REP -g hs 
  -p 0.01 --nomodel --shift 0 --extsize $FRAGLEN --keep-dup all -B --SPMR
\end{verbatim}

Generate fold enrichment signal tracks in bedgraph format.

\begin{verbatim}
$ macs2 bdgcmp -t "$REP"_treat_pileup.bdg -c "$REP"_control_lambda.bdg 
                                                      --outdir . -o "$REP"_FE.bdg -m FE
\end{verbatim}

Remove coordinates outside chromosome sizes.

\begin{verbatim}
$ slopBed -i "$REP"_FE.bdg -g genome_data/hg38/hg38.chrom.sizes -b 0 
  | awk '{if ($3 != -1) print $0}' 
  | bedClip stdin genome_data/hg38/hg38.chrom.sizes $REP.fc.bedgraph
\end{verbatim}

Convert \verb|bedgraph| to \verb|bigwig| for fold enrichment signal tracks.

\begin{verbatim}
$ sort -k1,1 -k2,2n $REP.fc.bedgraph > signal/$REP.fc.srt.bedgraph
$ bedGraphToBigWig signal/$REP.fc.srt.bedgraph genome_data/hg38/hg38.chrom.sizes  
  signal/$REP.fc.signal.bigwig
\end{verbatim}

Generate $-log_{10}$($p$-value) signal tracks in \verb|bedgraph| format.

\begin{verbatim}
$ chipReads=$(zcat align/$REP.tagAlign.gz | wc -l | awk '{printf "%f", $1/1000000}')
\end{verbatim}

Compute $sval = min$(no. of reads in ChIP, no. of reads in control) / 1,000,000

\begin{verbatim}
$ controlReads=$(zcat align/$CTL.tagAlign.gz_tag | wc -l 
  | awk '{printf "%f", $1/1000000}')
$ sval=$(echo "${chipReads} ${controlReads}" 
  | awk '$1>$2{printf "%f",$2} $1<=$2{printf "%f",$1}')"
$ macs2 bdgcmp -t "$REP"_treat_pileup.bdg -c "$REP"_control_lambda.bdg 
  --outdir . -o "$REP"_ppois.bdg -m ppois -S "${sval}"
\end{verbatim}

Remove coordinates outside chromosome sizes.

\begin{verbatim}
$ slopBed -i "$REP"_ppois.bdg -g genome_data/hg38/hg38.chrom.sizes -b 0 | awk '{if ($3 != -1) 
  print $0}' | bedClip stdin genome_data/hg38/hg38.chrom.sizes $REP.pval.bedgraph
\end{verbatim}

Convert \verb|bedgraph| to \verb|bigwig| for $p$-val signal tracks.

\begin{verbatim}
$ sort -k1,1 -k2,2n $REP.pval.bedgraph > signal/$REP.pval.srt.bedgraph
$ bedGraphToBigWig signal/$REP.pval.srt.bedgraph genome_data/hg38/hg38.chrom.sizes 
  signal/$REP.pval.signal.bigwig
\end{verbatim}

Repeat for IP replicate 2 (\verb|REP=rep2|; \verb|FRAGLEN=$fraglen_rep2|; \verb|OUTDIR=signal/reps|) and pooled replicate (\verb|REP=rep0|; \verb|FRAGLEN=$fraglen_mean|; \verb|OUTDIR=signal/pooled_reps|).
Remove temporary files.

\begin{verbatim}
$ rm -f *.bdg *.bedgraph signal/*.bedgraph *.xls *.narrowPeak *.bed
\end{verbatim}

\item Perform IDR analysis.  

Please load Conda environment for IDR analysis. Please make sure that (\verb|aquas_chipseq_py|) is prepended to the bash prompt after loading Conda environments.

\begin{verbatim}
$ source activate aquas_chipseq_py3
\end{verbatim}

For each set of peak calls of IP replicate, its pseudo-replicates and pooled pseudo-replicates, perform IDR analysis.  For the details of the input and output format, please see Box 4. Start with IP replicate 1.

\begin{verbatim}
$ PEAK1=peak/reps/rep1.regionPeak.gz; PEAK2=peak/reps/rep2.regionPeak.gz; 
  PEAK_POOLED=peak/pooled/rep0.regionPeak.gz; OUTDIR=idr/true_reps; PREFIX=rep1_vs_rep2
$ idr --samples $PEAK1 $PEAK2 --peak-list $PEAK_POOLED --input-file-type narrowPeak 
  --output-file $PREFIX.relaxed_peak.18-col.bed --rank signal.value 
  --soft-idr-threshold 0.05 --plot --use-best-multisummit-IDR
$ idr_thresh_transformed=$(awk -v p=0.05 'BEGIN{print -log(p)/log(10)}')
$ awk 'BEGIN{OFS="\t"} $12>='"${idr_thresh_transformed}"' {if ($2<0) $2=0; 
  print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10}' $PREFIX.relaxed_peak.18-col.bed 
  | sort | uniq | sort -k7n,7n | gzip -nc > $OUTDIR/$PREFIX.narrowPeak.gz
\end{verbatim}

Filter out peaks in blacklisted regions.

\begin{verbatim}
$ bedtools intersect -v -a <(zcat -f $OUTDIR/$PREFIX.narrowPeak.gz) -b <(zcat -f 
  genome_data/hg38/hg38.blacklist.bed.gz) | grep -P 'chr[\dXY]+[ \t]' 
  | awk 'BEGIN{OFS="\t"} {if ($5>1000) $5=1000; print $0}' 
  | gzip -nc > $OUTDIR/$PREFIX.filt.narrowPeak.gz
\end{verbatim}

Rename plot file names. IDR plots include number of significant peaks and IDR values (Figure 8; see Box 4 for the interpretation of each panel).

\begin{verbatim}
$ mv $PREFIX.relaxed_peak.18-col.bed.noalternatesummitpeaks.png qc/IDR.$PREFIX.idr_plot.png
\end{verbatim}

Repeat for pseudo-replicates of IP replicate 1 \\ 
(\verb|PEAK1=peak/self_pseudo_reps/rep1.pr1.regionPeak.gz;| \\ 
 \verb|PEAK2=peak/self_pseudo_reps/rep1.pr2.regionPeak.gz;| \\ 
 \verb|PEAK_POOLED=peak/reps/rep1.regionPeak.gz;| \\ 
 \verb|OUTDIR=idr/self_pseudo_reps;| \verb|PREFIX=rep1_pr|), \\ 
 pseudo-replicates of IP replicate 2 \\ 
 (\verb|PEAK1=peak/self_pseudo_reps/rep2.pr1.regionPeak.gz;|\\ 
  \verb|PEAK2=peak/self_pseudo_reps/rep2.pr2.regionPeak.gz;| \\ 
  \verb|PEAK_POOLED=peak/reps/rep2.regionPeak.gz;| 
  \verb|OUTDIR=idr/self_pseudo_reps;| \verb|PREFIX=rep2_pr|),\\ 
 and pooled pseudo-replicates \\
 (\verb|PEAK1=peak/pooled_pseudo_reps/rep0.pr1.regionPeak.gz;| \\ 
 \verb|PEAK2=peak/pooled_pseudo_reps/rep0.pr2.regionPeak.gz;| \\ 
 \verb|PEAK_POOLED=peak/pooled/rep0.regionPeak.gz;| \\ 
 \verb|OUTDIR=idr/pooled_pseudo_reps;| \verb|PREFIX=ppr|).

Remove temporary files.

\begin{verbatim}
$ rm -f *.bed
\end{verbatim}

\item Calculate thresholds to truncate peak lists. We will use the IDR thresholds of 0.05 for the original replicates threshold, self-consistency thresholds and pooled-pseudoreplicate threshold. For the recommended IDR thresholds, please see Box 6.

Reload Conda environment for postprocess IDR outputs.

\begin{verbatim}
$ source activate aquas_chipseq
\end{verbatim}

For the original (true) replicates:

\begin{verbatim}
$ Nt=$(zcat idr/true_reps/rep1_vs_rep2.filt.narrowPeak.gz | wc -l)
\end{verbatim}

Let’s call the maximum value of these three values $N_t$ (see Figure 1). In this example, $N_t$ is .....

For self pseudo replicates:

\begin{verbatim}
$ N1=$(zcat idr/self_pseudo_reps/rep1_pr.filt.narrowPeak.gz | wc -l)
$ N2=$(zcat idr/self_pseudo_reps/rep2_pr.filt.narrowPeak.gz | wc -l)
\end{verbatim}

These two values are the self-consistency thresholds, $N_1$ and $N_2$, respectively (see Figure 1) and should be in a similar range. If any of these is substantially different from the others, it points that replicate being different in some way. We use it to flag replicates with low consistency, as we will explain in the next step. In this example, $N_1$ and $N_2$ are ... and ..., respectively.

For the pooled-pseudoreplicate:

\begin{verbatim}
$ Np=$(zcat idr/pooled_pseudo_reps/ppr.filt.narrowPeak.gz | wc -l)
\end{verbatim}

This value is the pooled pseudo replicate, $N_p$ in Figure 1. In this example, $N_p$ is ....

\item Flag replicates for low consistency. Typically the self-consistency thresholds for all the original replicates should be reasonably similar (within a factor of 2). $N_t$ (comparing original replicates) and $N_p$ (comparing pooled-pseudoreplicates) should be similar (within a factor of 2). We typically flag datasets that fail both the above mentioned criteria above. We also flag datasets that fail either one of these criteria by a significant margin (i.e the above these ratios are significantly $\gg$ 2). For example, in this example, the rescue ratio $N_p/N_t$ does not exceed 2.

\textbf{CRITICAL STEP:} You should also check other measures such as the strand cross-correlation based quality metrics, sequencing depth differences and library complexity metrics for the replicates to see if they can further support data quality problems.

\item Obtain the final set of peak calls. You can now generate a conservative and an optimal final set of peak calls for the dataset.

For the optimal set, if $N_t > N_p$ choose IDR peaks (\verb|idr/true_reps/rep1_vs_rep2.filt.narrowPeak.gz|) from true replicate. Otherwise choose IDR peaks (\verb|idr/pooled_pseudo_reps/ppr.filt.narrowPeak.gz|) from pooled pseudo replicates.

For the conservative set, always choose IDR peaks (\verb|idr/true_reps/rep1_vs_rep2.filt.narrowPeak.gz|) from true replicate.

\item Check fraction of reads in peaks (FRiP). You can compute the fraction of reads from each replicate tagAlign and pooled tagAlign that fall within the $N_p$ and $N_t$ peak sets. Start with self pseudo replicates of IP replicate 1.

\begin{verbatim}
$ REP=rep1; FRAGLEN=$fraglen_rep1; PEAK=idr/self_pseudo_reps/rep1_pr.filt.narrowPeak.gz; 
  FRiP=qc/IDR.rep1_pr.FRiP.txt
$ HALF_FRAGLEN=$(( (FRAGLEN+1)/2 ))
$ val1=$(bedtools slop -i align/$REP.tagAlign.gz -g genome_data/hg38/hg38.chrom.sizes -s 
  -l -$HALF_FRAGLEN -r $HALF_FRAGLEN | awk '{if ($2>=0 && $3>=0 && $2<=$3) print $0}' 
  | bedtools intersect -a stdin -b <(zcat -f $PEAK) -wa -u | wc -l)
$ val2=$(zcat align/$REP.tagAlign.gz | wc -l)
$ awk 'BEGIN {print '${val1}'/'${val2}'}' > $FRiP
\end{verbatim}

Repeat for self pseudo replicates of IP replicate 2 \\ 
(\verb|REP=rep2;| \verb|FRAGLEN=$fraglen_rep2;| \\ 
\verb|PEAK=idr/self_pseudo_reps/rep2_pr.filt.narrowPeak.gz;|\\ 
 \verb|FRiP=qc/IDR.rep2_pr.FRiP.txt|), 
 pooled pseudo replicate \\ 
 (\verb|REP=rep0|; \verb|FRAGLEN=$fraglen_mean;| \\ 
 \verb|PEAK=idr/pooled_pseudo_reps/ppr.filt.narrowPeak.gz;| \\ 
 \verb|FRiP=qc/IDR.ppr.FRiP.txt|),\\ 
  and true IP replicates \\ 
  (\verb|REP=rep0|; \verb|FRAGLEN=$fraglen_mean;| \\ 
  \verb|PEAK=idr/true_reps/rep1_vs_rep2.filt.narrowPeak.gz|; \\ 
  \verb|FRiP=qc/IDR.rep1_vs_rep2.FRiP.txt|).

\end{enumerate}

\section*{Acknowledgments}

\hl{We thank Sol Katzman, Frank Jacobs and David Haussler for providing the paired-end ChIP-seq datasets. We thank Ewan Birney, Ian Dunham, Michael Snyder, Rick Myers, Terry Furey, Greg Crawford, Jason Lieb, Georgi Marinov, Florencia Pauli, Preti Jain, Ross Hardison, Kate Rosenbloom, Jim Kent and all members of the ENCODE consortium for the invaluable data, comments and feedback. We also thank the European Bioinformatics Institute for providing extensive data storage space and computational resources. We would like to acknowledge fly modENCODE chromatin group and regulatory elements group. We also thank Psalm Mizuki for helping with data preparation and Nils Gehlenborg for the discussion on figures. This work was supported by the National Institutes of Health Grant [RC2 HG005639 to P.J.P. and , R01 GM109453 to Q.L.], and the National Science Foundation [grant number 0640211 to S.B.]. Funding for open access charge: National Institutes of Health.}

\section*{Author contributions} 

\hl{A.K., P.V.K., Y.L.J. and P.J.P. developed the cross-correlation measures. A.K. and Q.L. developed the IDR framework for a ChIP-seq pipeline. J.R., S.W., and J.B.B. performed some tests for the IDR ChIP-seq pipeline. A.K., Y.L.J. and Q.L. performed the analysis for the study cases. A.K., Y.L.J., Q.L. and P.J.P. wrote the manuscript. N.B., B.J.W, A.S. and S.B. contributed to the discussion and protocol development.}

\section*{COMPETING FINANCIAL INTERESTS} 

The authors declare no competing financial interests.

\begin{multicols}{2}

\begin{thebibliography}{100}

% \section*{References}

\input{references}

\end{thebibliography}

\end{multicols}

\end{document}

\clearpage

\setcounter{table}{0}
\renewcommand{\tablename}{Supplementary Table}
\setcounter{figure}{0}
\renewcommand{\figurename}{Supplementary Figure}

\setcounter{page}{1}
\renewcommand\thepage{{SM }\arabic{page}}

\begin{center}
% {\LARGE \textbf{\begin{spacing}{1.1}The Cost of a Gene. \\ Supplementary Materials\end{spacing} }}
{\LARGE \textbf{Supplementary Materials}}
\end{center}

\section*{Supplementary Tables}

\section*{Supplementary Figures}

