diffHicUG.Rnw

\documentclass{report}

<<style-knitr, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex(relative.path=TRUE)
@

\bioctitle[\Biocpkg{diffHic} User's Guide]{\Biocpkg{diffHic}: Differential analysis of Hi-C data \\ User's Guide}
%% also: \bioctitle{Title used for both header and title page}
%% or... \title{Title used for both header and title page}
\author[1]{Aaron T. L. Lun}
\affil[1]{The Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia}

\usepackage{framed,color}
\newenvironment{combox}
{ \definecolor{shadecolor}{RGB}{255, 240, 240} \begin{shaded}\begin{center}\begin{minipage}[t]{0.95\textwidth} }
{ \end{minipage}\end{center}\end{shaded} \definecolor{shadecolor}{RGB}{240,240,240} }

\usepackage{bibentry}
\nobibliography*

\begin{document}
\maketitle
            
\begin{abstract}
This document contains instructions on how to use \Biocpkg{diffHic} for differential interaction analyses of Hi-C data.
It covers counting of read pairs into bin pairs, filtering out low-abundance bin pairs, normalization of sample-specific trended biases, 
variance modelling and hypothesis testing, and summarization and visualization of bin pairs.

\vspace{0.05in}
\textit{First edition:} 12 December 2012 \\
\textit{Last revised:} 10 October 2017 \\
\textit{Last compiled:} \Sexpr{format(Sys.Date(), "%d %B %Y")}
\end{abstract}

\packageVersion{\Sexpr{BiocStyle::pkg_ver("diffHic")}}

\tableofcontents
\newpage

<<results="hide",echo=FALSE>>=
picdir <- "plots-ug/"
if (file.exists(picdir)) { unlink(picdir, recursive=TRUE) }
dir.create(picdir)
knitr::opts_chunk$set(error=FALSE, warning=FALSE, message=FALSE)
knitr::opts_chunk$set(fig.path=picdir, dev='png', dpi=150)
@

\chapter{Introduction}
\section{Scope}
This document describes the differential analysis of Hi-C data with the \Biocpkg{diffHic} package.
Differential interactions (DIs) are defined as interactions with significant changes in intensity between experimental conditions.
DIs are identified in a statistically rigorous manner using methods in the \Biocpkg{edgeR} package \cite{edgeR}.
Knowledge of \Biocpkg{edgeR} is useful but is not necessary for this guide.

\section{How to get help}
Most questions about individual functions should be answered by the documentation.
For example, if you want to know more about \Rfunction{preparePairs}, you can bring up the documentation by typing \Rcode{?preparePairs} or \Rcode{help(preparePairs)} at the \R{} prompt.
Further detail on the methods or the theory can be found in the references at the bottom of each help page.
The entire \Biocpkg{diffHic} pipeline is also described in its own publication \cite{lun2015diffhic}.

The authors of the package always appreciate receiving reports of bugs in the package functions or in the documentation. 
The same goes for well-considered suggestions for improvements. 
Other questions about how to use \Biocpkg{diffHic} are best sent to the Bioconductor support site at \url{https://support.bioconductor.org}.
Please send requests for general assistance and advice to the support site, rather than to the individual authors. 
Users posting to the support site for the first time may find it helpful to read the posting guide at \url{http://www.bioconductor.org/help/support/posting-guide}.

\section{A brief description of Hi-C}
The Hi-C protocol \cite{lieberman2009comprehensive} is used to study chromatin organization by identifying pairwise interactions between two distinct genomic loci.
Briefly, chromatin is cross-linked and digested with a restriction enzyme.
This releases chromatin complexes into solution, where each complex contains multiple restriction fragments corresponding to interacting loci.
Overhangs are filled in with biotin-labelled nucleotides to form blunt ends.
Proximity ligation is then performed, which favors ligation between blunt ends in the same complex.
The ligated DNA is sonicated and the sheared fragments containing ligation junctions are purified by a streptavidin pulldown.
Paired-end sequencing is performed on the purified ligation products, and pairs of interacting loci are identified by mapping the paired reads.
Of course, some caution is required due to the presence of non-specific ligation between blunt ends in different complexes.

\section{Quick start}
A typical differential analysis of Hi-C data is described below.
For simplicity, assume that the BAM files have already been processed into ``index'' files in \Robject{input}.
Let \Robject{design} contain the design matrix for this experiment.
Also assume that the boundaries of the relevant restriction fragments are present in \Robject{fragments}.
The code itself is split across several steps: 

<<echo=FALSE>>=
input <- c("merged_flox_1.h5", "merged_flox_2.h5", "merged_ko_1.h5", "merged_ko_2.h5")
fragments <- readRDS("mm10-hindIII.rds")
design <- model.matrix(~factor(c("flox", "flox", "ko", "ko")))
@

% saveRDS(fragments, file="mm10-hindIII.rds")

\begin{enumerate}
\item Converting BAM files to index files.
\item Counting read pairs into pairs of genomic bins.
<<>>=
library(diffHic)
param <- pairParam(fragments=fragments)
data <- squareCounts(input, width=1e6, param=param)
@
\item Filtering out uninteresting bin pairs.
<<>>=
library(edgeR)
keep <- aveLogCPM(asDGEList(data)) > 0
data <- data[keep,]
@
\item Normalizing for sample-specific biases.
<<>>=
data <- normOffsets(data, type="loess", se.out=TRUE)
y <- asDGEList(data)
@
\item Modelling biological variability between replicates.
<<>>=
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design, robust=TRUE)
@
\item Testing for significant differences between groups.
<<>>=
result <- glmQLFTest(fit)
clusters <- diClusters(data, result$table, target=0.05, 
                       cluster.args=list(tol=1))
@
\end{enumerate}
In the various examples for this guide, data will be used from three studies. 
The first dataset examines the chromatin structure in K562 and GM06990 cell lines \cite{lieberman2009comprehensive}.
The second compares interaction intensities between wild-type and cohesin-deficient murine neural stem cells \cite{sofueva2013cohesin}. 
The final study compares ERG-overexpressing RWPE1 cells with a GFP-expressing control \cite{rickman2012oncogene}.
Obviously, readers will have to modify the code for their own analyses.

\chapter{Preparing index files from BAM files}
\label{chap:prep} 

\begin{combox}
This box will appear at the start of each chapter, describing the objects required from previous chapters.
As we're starting out here, we don't really need anything.
\end{combox}

\section{A comment on aligning Hi-C libraries}
\label{sec:align}

In a typical Hi-C library, sequencing will occasionally be performed over the ligation junction between two restriction fragments.
This forms a chimeric read that contains sequences from two distinct genomic loci.
Here, correct alignment of the 5$'$ end of the chimeric read is more important.
This is because the location of the 3$'$ end is already provided by mapping the 5$'$ end of the mate read.
Direct application of alignment software will not be optimal as only one mapping location will be usually reported for each read.
This means that the 5$'$ location will be ignored if the 3$'$ alignment is superior, e.g., longer or fewer mismatches.

Instead, chimeric alignment can be performed with approaches like iterative mapping \cite{imakaev2012iterative} or read splitting \cite{seitan2013cohesin}.
In the latter, we define the ligation signature as the sequence that is obtained after ligation between blunt ends derived from cleaved restriction sites.
For example, the \textit{Hin}dIII enzyme cleaves at \Rcode{AAGCTT} with a 4 bp overhang.
This yields a signature sequence of \Rcode{AAGCTAGCTT} upon ligation of blunt ends.
The ligation signature in each chimeric read is identified with \software{cutadapt} \cite{martin2011cutadapt}, and the read is split into two segments at the center of the signature.
Each segment is then aligned separately to the reference genome using \software{Bowtie2} \cite{langmead2012bowtie}.
Mapping by read splitting is performed in \Biocpkg{diffHic} using a custom Python script:

<<results='hide'>>=
system.file("python", "presplit_map.py", package="diffHic", mustWork=TRUE)
@

Users are strongly recommended to synchronize mate pair information and to mark duplicate read pairs in the resulting BAM file.
This can be achieved using the various tools in the \software{Picard} software suite (\url{http://broadinstitute.github.io/picard}).

\section{Matching mapped reads to restriction fragments}
The Hi-C protocol is based on ligation between restriction fragments.
Sequencing of the ligation product is performed to identify the interacting loci --  or, more precisely, the two restriction fragments containing the interacting loci.
The resolution of Hi-C data is inherently limited by the frequency of restriction sites and the size of the restriction fragments.
Thus, it makes sense to report the read alignment location in terms of the restriction fragment to which that read was mapped.
The boundaries of each restriction fragment can be obtained with the \Rfunction{cutGenome} function, as shown below for the human genome after digestion with the \textit{Hin}dIII restriction enzyme (recognition site of \Rcode{AAGCTT}, 5$'$ overhang of 4 bp).

<<>>=
library(BSgenome.Hsapiens.UCSC.hg19)
hs.frag <- cutGenome(BSgenome.Hsapiens.UCSC.hg19, "AAGCTT", 4)
hs.frag
@

These fragments should be stored in a \Rclass{pairParam} object.
The constructor below checks that the fragment ranges are properly ordered.
Later, this object will also hold other parameters for counting.
This simplifies coordination of the various steps in the \Biocpkg{diffHic} pipeline, as the same \Rclass{pairParam} object can be easily passed between different functions. 

<<>>=
hs.param <- pairParam(hs.frag)
hs.param
@

The \Rfunction{preparePairs} function matches the mapping location of each read to a restriction fragment in the reference genome.
Mapping results for paired reads should be provided in a name-sorted BAM file.
The function converts the read position into an index, pointing to the matching restriction fragment in \Robject{hs.frag}.
The resulting pairs of indices are stored in an index file using the HDF5 format.
The fragments to which the reads mapped are referred to as ``anchors'' -- the larger index is defined as the first anchor fragment whereas the smaller is the second.
This process is demonstrated using Hi-C data generated from GM06990 cells.

<<>>=
diagnostics <- preparePairs("SRR027957.bam", hs.param, file="SRR027957.h5", 
                            dedup=TRUE, minq=10)
names(diagnostics)
@

The function itself returns a list of diagnostics showing the number of read pairs that are lost for various reasons.
Of particular note is the removal of reads that are potential PCR duplicates with \Rcode{dedup=TRUE}.
This requires marking of the reads beforehand using an appropriate program such as \software{MarkDuplicates} from the \software{Picard} suite.
Filtering on the minimum mapping quality score with \Rcode{minq} is also recommended to remove spurious alignments.

<<>>=
diagnostics$pairs
@

Read pairs mapping to the same restriction fragment provide little information on interactions between fragments \cite{belton2012hic}.
Dangling ends are inward-facing read pairs that are mapped to the same fragment.
These are uninformative as they are usually formed from sequencing of the restriction fragment prior to ligation.
Self-circles are outward-facing read pairs that are formed when two ends of the same restriction fragment ligate to one another.
Interactions within a fragment cannot be easily distinguished from these self-circularization events.
Both structures are removed to avoid confusion in downstream steps.

<<>>=
diagnostics$same.id
@

\section{Processing of chimeric reads}
In the BAM file, chimeric reads are represented by two separate alignments for each read.
Hard clipping is used to denote the length trimmed from each sequence in each alignment, and to determine which alignment corresponds to the 5$'$ or 3$'$ end of the read.
Only the 5$'$ end(s) is used to determine the restriction fragment index for that read pair.
The total number of chimeric read pairs will be reported by \Rfunction{preparePairs}, along with the number of pairs where both 5$'$ ends are mapped (\Rcode{mapped}) and the nuber where both 5$'$ ends and at least one 3$'$ end are mapped (\Rcode{multi}).
Of course, the function will also work if the mapping location is only given for the 5$'$ end, though chimeric statistics will not be properly computed.

<<>>=
diagnostics$chimeras
@

The proportion of invalid chimeric pairs is also calculated by \Rfunction{preparePairs}.
Invalid pairs are those where the 3$'$ location of a chimeric read disagrees with the 5$'$ location of the mate.
The invalid proportion can be used as an empirical measure of the mapping error rate - or, at least, the upper bound thereof, given that error rates are likely to be lower for longer, non-chimeric alignments.
High error rates may be indicative of a fault in the mapping pipeline.

<<>>=
diagnostics$chimeras[["invalid"]]/diagnostics$chimeras[["multi"]]
@

Invalid chimeric pairs can be discarded by setting \Rcode{ichim=FALSE} in \Rfunction{preparePairs}.
However, this is not recommended for routine analyses.
Mapping errors for short 3$'$ ends may result in apparent invalidity and loss of the (otherwise correct) 5$'$ alignments.

\section{Filtering artifactual read pairs}

\subsection{Reprocessing index files for quality control}
The \Rfunction{prunePairs} function removes read pairs that correspond to artifacts in the Hi-C procedure.
The returned vector contains the number of read pairs removed for each type of artifact.
Values of \Rcode{length}, \Rcode{inward} and \Rcode{outward} correspond to removal by \Rcode{max.frag}, \Rcode{min.inward} and \Rcode{min.outward}, respectively.
Retained read pairs are stored in another index file for later use.

<<>>=
min.inward <- 1000
min.outward <- 25000
prunePairs("SRR027957.h5", hs.param, file.out="SRR027957_trimmed.h5", 
           max.frag=600, min.inward=min.inward, min.outward=min.outward)
@

The \Rcode{max.frag} argument removes read pairs where the inferred length of the sequencing fragment (i.e., the ligation product) is greater than a specified value.
The length of the sequencing fragment is computed by summing, for both reads, the distance between the mapping location of the 5$'$ end and the nearest restriction site on the 3$'$ side.
Excessively large lengths are indicative of off-site cleavage, i.e., where the restriction enzyme or some other agent cuts the DNA at a location other than the restriction site.
While not completely uninformative, these are discarded as they are not expected from the Hi-C protocol.
The threshold value can be chosen based on the size selection interval in library preparation, or by examining the distribution of inferred fragment lengths from \Rfunction{getPairData}.

<<fraglenhist, fig.small=TRUE>>=
diags <- getPairData("SRR027957.h5", hs.param)
hist(diags$length[diags$length < 1000], ylab="Frequency", 
     xlab="Spacing (bp)", main="", col="grey80")
@

The insert size is defined as the linear distance between two paired reads on the same chromosome.
The \Rcode{min.inward} paramater removes inward-facing intra-chromosomal read pairs where the insert size is less than the specified value.
The \Rcode{min.outward} parameter does the same for outward-facing intra-chromosomal read pairs.
This is designed to remove dangling ends or self-circles involving DNA fragments that have not been completely digested \cite{jin2013highres}.
Such read pairs are technical artifacts that are (incorrectly) retained by \Rfunction{preparePairs}, as the two reads involved are mapped to different restriction fragments.
Appropriate thresholds for both parameters can be determined using strand orientation plots.

\subsection{Setting size thresholds with strand orientation plots}
\label{sec:strorient}
The strand orientation for a read pair refers to the combination of strands for the two alignments.
These are stored as flags where setting \Rcode{0x1} or \Rcode{0x2} means that the read on the first or second anchor fragment, respectively, is mapped on the reverse strand.
If different pieces of DNA were randomly ligated together, one would expect to observe equal proportions of all strand orientations. 
This can be tested by examining the distribution of strand orientations for inter-chromosomal read pairs.
Each orientation is equally represented across these read pairs, which is expected as different chromosomes are distinct pieces of DNA.

<<>>=
intra <- !is.na(diags$insert)
table(diags$orientation[!intra])
@

This can be repeated for intra-chromosomal read pairs, by plotting the distribution of insert sizes for each strand orientation \cite{jin2013highres}.
The two same-strand distributions are averaged for convenience.
At high insert sizes, the distributions will converge for all strand orientations.
This is consistent with random ligation between two separate restriction fragments.
At lower insert sizes, spikes are observed in the ouward- and inward-facing distributions due to self-circularization and dangling ends, respectively.
Thresholds should be chosen in \Rfunction{prunePairs} to remove these spikes, as represented by the grey lines.

<<strorient, fig.small=TRUE>>=
# Constructing the histograms for each orientation.
llinsert <- log2(diags$insert + 1L)
intra <- !is.na(llinsert)
breaks <- seq(min(llinsert[intra]), max(llinsert[intra]), length.out=30)
inward <- hist(llinsert[diags$orientation==1L], plot=FALSE, breaks=breaks)
outward <- hist(llinsert[diags$orientation==2L] ,plot=FALSE, breaks=breaks)
samestr <- hist(llinsert[diags$orientation==0L | diags$orientation==3L], 
                plot=FALSE, breaks=breaks)
samestr$counts <- samestr$counts/2

# Setting up the axis limits.
ymax <- max(inward$counts, outward$counts, samestr$counts)/1e6
xmax <- max(inward$mids, outward$mids, samestr$mids)
xmin <- min(inward$mids, outward$mids, samestr$mids)

# Making a plot with all orientations overlaid.
plot(0,0,type="n", xlim=c(xmin, xmax), ylim=c(0, ymax),
     xlab=expression(log[2]~"[insert size (bp)]"), ylab="Frequency (millions)")
lines(inward$mids, inward$counts/1e6, col="darkgreen", lwd=2)
abline(v=log2(min.inward), col="darkgrey")
lines(outward$mids, outward$counts/1e6, col="red", lwd=2)
abline(v=log2(min.outward), col="darkgrey", lty=2)
lines(samestr$mids, samestr$counts/1e6, col="blue", lwd=2)
legend("topright", c("inward", "outward", "same"), 
       col=c("darkgreen", "red", "blue"), lwd=2)
@

As an aside, the position of the spikes in the above plots can be used to estimate some fragment lengths.
The $x$-coordinate of the outward-facing spike represents the average length of the DNA fragments after restriction digestion.
This is useful as it provides a lower bound on the spatial resolution of any given Hi-C experiment.
The position of the inward-facing spike represents the average length of the fragments after sonication.
This should be lower than the size selection thresholds used in library preparation.

\section{Merging technical replicates}
Hi-C experiments often involve deep sequencing as read pairs are sparsely distributed across the set of possible interactions.
As a result, multiple index files may be generated from multiple technical replicates of a single Hi-C library.
These can be merged together using the \Rfunction{mergePairs} function prior to downstream processing.
This is equivalent to summing the counts for each pair of restriction fragment indices, and is valid if one assumes Poisson sampling for each sequencing run \cite{marioni2008rnaseq}.
An example is provided below that merges several technical replicates for a Hi-C library generated from GM06990 cells \cite{lieberman2009comprehensive}. 

<<>>=
prepped <- preparePairs("SRR027958.bam", hs.param, file="SRR027958.h5", 
                        dedup=TRUE, minq=10)
counted <- prunePairs("SRR027958.h5", hs.param, file.out="SRR027958_trimmed.h5", 
                      max.frag=600, min.inward=min.inward, 
                      min.outward=min.outward)
mergePairs(files=c("SRR027957_trimmed.h5", "SRR027958_trimmed.h5"), "merged.h5")
@

In addition, any Hi-C dataset that is processed manually by the user can be stored in an index file using the \Rfunction{savePairs} function.
This takes a dataframe with the first and second anchor indices, as well as any additional information that might be useful.
The idea is to provide an entry point into the \Biocpkg{diffHic} analysis from other pipelines.
If the dataset is too large, one can save chunks at a time before merging them all together with \Rfunction{mergePairs}.

<<>>=
anchor1.id <- as.integer(runif(100, 1, length(hs.param$fragments)))
anchor2.id <- as.integer(runif(100, 1, length(hs.param$fragments)))
dummy <- data.frame(anchor1.id, anchor2.id, other.data=runif(100))
savePairs(dummy, "example.h5", hs.param)
@

For full compatibility, users should include the alignment positions and lengths for each read pair as \Rcode{anchorX.pos} and \Rcode{anchorX.len},
    where \Rcode{X} is 1 or 2 for each of the paired reads.
The alignment position refers to the 1-based coordinate of the left-most base of the alignment.
The alignment length refers to the span of the alignment relative to the reference, and should be negative for alignments on the reverse strand.
This information will be used in downstream \Biocpkg{diffHic} functions, such as read counting around blacklisted regions.

\section{Handling DNase Hi-C experiments}
DNase Hi-C is a variant of the standard protocol whereby the DNase~I enzyme is used to fragment the genome instead of restriction enzymes \cite{ma2015dnase}.
Random fragmentation provides resolution beyond that of individual restriction fragments.
However, cleavage sites for DNase~I cannot be easily predicted to construct \Rcode{param\$fragments}.
Fragment indices have no meaning here because there are no restriction fragments for reads to be assigned to.

Instead, \Biocpkg{diffHic} handles this type of data by operating directly on the alignment position of each read.
This reflects the theoretical base-pair resolution of the data during quantification.
To indicate that the reads come from a DNase Hi-C experiment, an empty \Rclass{GRanges} object should be supplied as the \Rcode{fragments} in the \Rclass{pairParam} object.
Most \Biocpkg{diffHic} functions will detect this and behave appropriately to exploit the improved resolution.
An example of this approach is shown below, where the SRR027957 library is treated as a DNase Hi-C sample.

<<>>=
seg.frags <- emptyGenome(BSgenome.Hsapiens.UCSC.hg19)
preparePairs("SRR027957.bam", pairParam(seg.frags), file="SRR027957.h5", 
             dedup=TRUE, minq=10)
@

Unlike \Rfunction{preparePairs}, no diagnostic information regarding self-circles or dangling ends is reported here.
Such metrics are irrelevant when restriction fragments are not involved.
Similarly, the \Rcode{length} field in the output of \Rfunction{getPairData} is set to \Rcode{NA} and should be ignored.
This is because the length of the sequencing fragment cannot be generally computed without knowledge of the fragmentation sites.
As a result, the \Rcode{max.frag} argument should also be set to \Rcode{NA} in \Rfunction{prunePairs}.
Metrics for inward- and outward-facing read pairs are unaffected, as these are computed independently of the fragments.
See the documentation for individual functions to determine if/how their behaviour changes for DNase Hi-C data.

Users should also note that the alignment script described in Section~\ref{sec:align} is not appropriate for DNase Hi-C experiments.
This approach is based on splitting chimeric reads at the ligation signature, which is constructed from the recognition site of a restriction enzyme.
The sequence around the ligation junction is not well-defined when DNase~I is used for cleavage.
Instead, alignment programs should be used that can directly handle chimeric reads with arbitrary breakpoints in the genome, e.g., \software{BWA} \cite{li2010fast}.
A Python script for iterative mapping \cite{imakaev2012iterative} is also provided for this purpose.

<<results='hide'>>=
system.file("python", "iter_map.py", package="diffHic", mustWork=TRUE)
@


\chapter{Counting read pairs into interactions}
\label{chap:counting}

\begin{combox}
A different dataset will be used here, so we don't need anything from the last chapter.
Horses for courses; this dataset's a lot nicer for detecting differential interactions.
\end{combox}

\section{Overview}
Prior to any statistical analysis, the read pairs in a Hi-C library must be summarized into a count for each interaction.
This count is used as an experimental measure of the interaction intensity.
Specifically, each pairwise interaction is parameterized by two genomic intervals representing the interacting loci.
The count for that interaction is defined as the number of read pairs with one read mapping to each of the intervals.
Counting is performed for each sample in the dataset, such that each interaction is associated with a set of counts.

The interaction space is defined as the genome-by-genome space over which the read pairs are distributed.
Recall that each paired read is assigned to a restriction fragment index.
The interaction space contains all index pairs $(x, y)$ for $x, y \in [1 .. N]$, where $x \ge y$ and $N$ is the number of restriction fragments in the genome.
This can be visualized as the triangular space between $y=x$, $y=0$ and $x=N$ on the Cartesian plane.
A rectangular area in the interaction space represents a pairwise interaction between the genomic intervals spanned by the two adjacent sides of that rectangle.
The number of read pairs in this area is used as the count for the corresponding interaction.
Non-rectangular areas can also represent interactions, but these are more difficult to interpret and will not be considered here.

The examples shown here will use the neural stem cell dataset \cite{sofueva2013cohesin} in which wild-type cells are compared to cohesin-deficient cells.
Read processing has already been performed to construct an index file for each library.
Some additional work is required to obtain the restriction fragment coordinates for the \textit{Hin}dIII-digested mouse genome.

<<>>=
library(BSgenome.Mmusculus.UCSC.mm10)
mm.frag <- cutGenome(BSgenome.Mmusculus.UCSC.mm10, "AAGCTT", 4)
input <- c("merged_flox_1.h5", "merged_flox_2.h5", 
           "merged_ko_1.h5", "merged_ko_2.h5")
@

\section{Counting into bin pairs}

\subsection{Overview}
Here, the genome is partitioned into contiguous non-overlapping bins of constant size.
Each interaction is defined as a pair of these bins.
This approach avoids the need for prior knowledge of the loci of interest when summarizing Hi-C counts.
Counting of read pairs between paired bins is performed for multiple libraries using the \Rfunction{squareCounts} function.

<<>>=
bin.size <- 1e6
mm.param <- pairParam(mm.frag)
data <- squareCounts(input, mm.param, width=bin.size, filter=1)
data
@

This generates an \Rclass{InteractionSet} object containing information for multiple genomic interactions \cite{lun2016infrastructure}.
Each row of the object corresponds to an interaction, i.e., bin pair, while each column represents a library.
For each interaction, the coordinates of its two ``anchor'' bins are returned.
The first/second anchor notation is used for these bins, whereby the first anchor bin is that with the ``higher'' genomic coordinate.

<<>>=
head(anchors(data, type="first"))
head(anchors(data, type="second"))
@

The returned \Rclass{InteractionSet} object contains a count matrix with the number of read pairs for each interaction in each library.
It also contains a \Rcode{totals} vector in the \Rcode{colData}, which stores the total number of read pairs for each library.

<<>>=
head(assay(data))
data$totals
@

Bin pairs can also be filtered to remove those with to a count sum below \Rcode{filter}.
This removes uninformative bin pairs with very few read pairs, and reduces the memory footprint of the function.
A higher value of \Rcode{filter} may be necessary for analyses of large datasets in limited memory.
More sophisticated filtering strategies are discussed in Chapter~\ref{chap:filter}.

\subsection{Choosing a bin width}
The \Rcode{width} of the bin is specified in base pairs and determines the spatial resolution of the analysis.
Smaller bins will have greater spatial resolution as adjacent features can be distinguished in the interaction space.
Larger bins will have greater counts as a larger area is used to collect read pairs.
Optimal summarization will not be achieved if bins are too small or too large to capture the (changes in) intensity of the underlying interactions.

For this analysis, 1 Mbp bins are used to capture broad structural features.
This is also useful for demonstration purposes, as the counts are large enough for clear manifestation of biases in later chapters.
Of course, \Biocpkg{diffHic} is more than capable of handling smaller sizes (e.g., below 20 kbp) for higher-resolution analyses of looping interactions between specific elements.
In such cases, \Rcode{filter} should be increased to avoid excessive memory usage.

<<>>=
head(regions(data))
@

The boundary of each bin is rounded to the closest restriction site in \Rfunction{squareCounts}.
This is due to the inherent limits on spatial resolution in a Hi-C experiment.
The number of restriction fragments in each bin is recorded in the \Rcode{nfrags} field of the metadata.

Determination of the ideal bin size is not trivial as the features of interest are not usually known in advance.
Instead, repeated analyses with multiple bin sizes are recommended.
This provides some robustness to the choice of bin size.
Sharp interactions can be detected by pairs of smaller bins while diffuse interactions can be detected by larger bin pairs.
See Section~\ref{sec:mergebins} for more information on consolidating results from multiple bin sizes.

\section{Counting with pre-defined regions}

\subsection{Connecting counts between pairs of regions}
For some studies, prior knowledge about the regions of interest may be available.
For example, a researcher may be interested in examining interactions between genes.
The coordinates can be obtained from existing annotation, as shown below for the mouse genome.
Other pre-specified regions can also be used, e.g., known enhancers or protein binding sites.

<<>>=
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
gene.body <- genes(TxDb.Mmusculus.UCSC.mm10.knownGene)
strand(gene.body) <- "*" # Removing strand info, for simplicity.
@

Counting is directly performed for these defined regions using the \Rfunction{connectCounts} function.
Interactions are defined between each pair of regions in the pre-specified set.
This may be easier to interpret than pairs of bins if the interacting regions have some biological significance.
The count matrix and the vector of totals are defined as previously described.

<<>>=
redata <- connectCounts(input, mm.param, regions=gene.body)
redata
@

Again, first/second anchor notation applies whereby the interval with the larger start coordinate in the genome is defined as the first anchor.
Note that the anchor may not have a larger end coordinate if the supplied \Rcode{regions} are nested.
In addition, each region is rounded to the nearest restriction site.
Resorting is also performed, though the indices of the original regions can be found in the metadata as \Rcode{original} if back-referencing is necessary.

<<>>=
head(regions(redata))
@

One obvious limitation of this approach is that interactions involving unspecified regions will be ignored.
This is obviously problematic when searching for novel interacting loci.
Another issue is that the width of the regions cannot be easily changed.
This means that the compromise between spatial resolution and count size cannot be tuned.
For example, interactions will not be detected around smaller genes as the counts will be too small.
Conversely, interactions between distinct loci within a large gene body will not be resolved.

\subsection{Connecting counts between two sets of regions}
The \Rfunction{connectCounts} function can also be used to identify interactions between two sets of regions, by specifying a value for the \Rcode{second.regions} argument.
This only considers interactions between one entry in \Rcode{regions} and another in \Rcode{second.regions}.
This differs from the standard application of the function, which would consider an interaction between any pair of entries in \Rcode{regions}.
If an integer scalar is supplied as \Rcode{second.regions}, the second set is automatically defined as contiguous bins of that size across the genome.

The use of \Rcode{second.regions} is particularly useful in cases where there are defined ``viewpoints'' of interest, e.g., 4C-seq, Capture-C.
These viewpoints can be specified in \Rcode{regions}, as shown below for a set of mock probe locations for a hypothetical Capture-C experiment \cite{hughes2014analysis}.
Specifically, the viewpoint is defined as a 100 kbp bin centred at each capture probe.
The interaction profile across the rest of the genome can then be extracted by setting \Rcode{second.regions} to some bin size.
In this case, 100 kbp bins are used.

<<>>=
probe.loc <- GRanges(c("chr1", "chr2", "chr3"),
                     IRanges(c(1e7, 2e7, 3e7), c(1e7, 2e7, 3e7)))
viewpoints <- resize(probe.loc, fix="center", width=1e5)
viewdata <- connectCounts(input, mm.param, regions=viewpoints, 
                          second.regions=1e5)
head(anchors(viewdata, type="first"))
head(anchors(viewdata, type="second"))
@

As these results demonstrate, interactions are only considered if exactly one interacting locus is from the specified \Rcode{regions}.
The identity of the other locus (i.e., from \Rcode{second.regions}) can be determined based on the \Rcode{is.second} field in the \Rclass{GRanges} object.
This approach avoids loading irrelevant interactions when only specific viewpoints are of interest.

\section{Counting into single bins}
\label{sec:marginal}
The \Rfunction{marginalCounts} function counts the number of reads mapped inside each bin.
This effectively treats each Hi-C library as single-end data to quantify the genomic coverage of each bin.
One can use these ``marginal'' counts to determine whether there are systematic differences in coverage between libraries for a given bin.
This implies that copy number variations are present, which may confound detection of differential interactions.

<<>>=
margin.data <- marginCounts(input, mm.param, width=bin.size)
margin.data
@

Each row of the \Rclass{RangedSummarizedExperiment} contains the marginal read counts for a genomic bin.
For this dataset, there are no major changes in coverage for the vast majority of bins.
The most extreme events occur at low abundances and are unlikely to be precise. 
This suggests that a direct comparison of interaction intensities will be valid.
Remedial action in the presence of copy number changes is not trivial and will be discussed in Section~\ref{sec:copy}.

<<mamargin, fig.small=TRUE>>=
adjc <- cpm(asDGEList(margin.data), log=TRUE, prior.count=5)
smoothScatter(0.5*(adjc[,1]+adjc[,3]), adjc[,1]-adjc[,3],
              xlab="A", ylab="M", main="Flox (1) vs. Ko (1)")
@

\section{Additional parameter options}

\subsection{Restricting the input chromosomes}
\label{sec:restrictchr}

Users can restrict counting to particular chromosomes by setting the \Rcode{restrict} slot in the \Rclass{pairParam} object.
This is useful to ensure that only interactions between relevant chromosomes are loaded.
Sequences such as the mitochondrial genome, unassigned contigs or random chromosome segments can be ignored in routine analyses.

<<>>=
new.param <- reform(mm.param, restrict=c("chr1", "chr2"))
new.param
@

In addition, if \Rcode{restrict} is a $n$-by-2 matrix, count loading will be limited to the read pairs that are mapped between the $n$ specified pairs of chromosomes.
The example below considers all read pairs mapped between chromosomes 2 and 19.
This feature is useful when memory is limited, as each pair of chromosomes can be loaded and analyzed separately.

<<>>=
new.param <- reform(mm.param, restrict=cbind("chr2", "chr19"))
new.param
@

\subsection{Specifying regions to ignore}
Users can also discard alignments that lie within blacklisted regions by setting the \Rcode{discard} slot.
The aim is to eliminate reads within known repeat regions.
Such regions are problematic, as reads from several repeat units in the real genome may be collated into a single representative unit in the genome build.
This results in a sharp, spurious spike in interaction intensity.
The problem is exacerbated by different repeat copy numbers between biological conditions, resulting in spurious differential interactions due to changes in coverage.
Removal of reads in these repeats may be necessary to obtain reliable results.

<<>>=
dummy.repeat <- GRanges("chr1", IRanges(10000, 1000000))
new.param <- reform(mm.param, discard=dummy.repeat)
new.param
@

Coordinates of annotated repeats can be obtained from several different sources.
A curated blacklist of problematic regions is available from the ENCODE project \cite{dunham2012encode}, and can be obtained at \url{https://sites.google.com/site/anshulkundaje/projects/blacklists}.
This list is constructed empirically from the ENCODE datasets and includes obvious offenders like telomeres, microsatellites and some rDNA genes.
Alternatively, repeats can be predicted from the genome sequence using software like RepeatMasker.
These calls are available from the UCSC website (e.g., \url{hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/chromOut.tar.gz} for mouse) or they can be extracted from an appropriate masked \Rclass{BSgenome} object. 
Experience suggests that the ENCODE blacklist is generally preferable.
Repeat predictions tend to be aggressive such that too much of the genome (and interactions therein) will be discarded.

\subsection{Capping the read pairs per restriction fragment pair}
Incomplete removal of PCR duplicates or read pairs in repeat regions may result in spikes of read pairs within the interaction space.
The effect of these artifacts can be mitigated by capping the number of read pairs associated with each pair of restriction fragments.
This is done by specifying a value for the \Rcode{cap} slot.
Diffuse interactions should not be affected, as the associated read pairs will be distributed sparsely across many fragment pairs.
More caution is required if sharp interactions are present, i.e., interactions between 5 - 10 kbp regions.

<<>>=
new.param <- reform(mm.param, cap=5)
new.param
@

\section{Loading counts from existing count matrices}
It is possible to use counts from existing count matrices, by using coercion methods available to \Rclass{ContactMatrix} objects from the \Biocpkg{InteractionSet} package \cite{lun2016infrastructure}.
Assume that there are two \Rclass{ContactMatrix} objects (\Robject{cm1} and \Robject{cm2}) spanning the same region of the binned interaction space.
The \Rfunction{mergeCMs} function converts these two objects into a single \Rclass{InteractionSet} object for use in the \Biocpkg{diffHic} analysis pipeline.
Only bin pairs with a count sum greater than 10 are retained, analogous to the output obtained with a \Rfunction{squareCounts} call with \Rcode{filter=10}.

<<echo=FALSE,results="hide">>=
# Just for testing.
N <- 30
all.starts <- round(runif(N, 1, 100))
all.ends <- all.starts + round(runif(N, 5, 20))
all.regions <- GRanges(rep(c("chrA", "chrB"), c(N-10, 10)),
    IRanges(all.starts, all.ends))

Nr <- 10
Nc <- 20
all.anchor1 <- sample(N, Nr)
all.anchor2 <- sample(N, Nc)
counts <- matrix(rpois(Nr*Nc, lambda=10), Nr, Nc)
cm1 <- ContactMatrix(counts, all.anchor1, all.anchor2, all.regions)
counts <- matrix(rpois(Nr*Nc, lambda=10), Nr, Nc)
cm2 <- ContactMatrix(counts, all.anchor1, all.anchor2, all.regions)
@

<<>>=
data.com <- mergeCMs(cm1, cm2, filter=10)
@

This procedure facilitates input from other Hi-C data processing pipelines that return count matrices rather than BAM files.
It is easily extended to datasets involving more than two samples by simply including additional \Rclass{ContactMatrix} objects in the \Rfunction{mergeCMs} call.

\section{Summary}
Counting into bin pairs is the most general method for quantifying interaction intensity.
It does not require any prior knowledge regarding the regions of interest.
The bin size can be easily adjusted to obtain the desired spatial resolution.
It is also easier/safer to compare between bin pairs (e.g., during filtering) when each bin is roughly of the same size.
Thus, bin-based counting will be the method of choice for the rest of this guide.

For simplicity, all counting steps will be performed here with the default settings, i.e., no values for \Rcode{restrict}, \Rcode{discard} or \Rcode{cap}.
However, users are encouraged to use non-default values if it can improve the outcome of their analyses.
Non-default values should be recorded in a single \Rclass{pairParam} object for consistent use across all functions in the \Biocpkg{diffHic} pipeline.
This ensures that the same read pairs will always be extracted from each index file.

\chapter{Filtering out uninteresting interactions}
\label{chap:filter}

\begin{combox}
Here, we'll need the \Robject{data} object that was loaded in the previous chapter.
We'll also need \Robject{bin.size} and \Robject{mm.param}, as well as the file names in \Robject{input}.
\end{combox}

\section{Overview}

\subsection{Computing the average abundance in a NB model}
Filtering is often performed to remove uninteresting features in analyses of high-throughput experiments. 
This reduces the severity of the multiplicity correction and increases power among the remaining tests. 
The filter statistic should be independent of the $p$-value under the null hypothesis, but correlated to the $p$-value under the alternative \cite{bourgon2010independent}. 
The aim is to enrich for false nulls without affecting type I error for the true nulls.

Assume that the counts for each bin pair are sampled from the negative binomial (NB) distribution.
In this model, the overall NB mean across all libraries is (probably) an independent filter statistic.
The log-mean-per-million is known as the average abundance and can be computed with the \Rfunction{aveLogCPM} function in \Biocpkg{edgeR} \cite{mccarthy2012glm}.

<<avehistplot, fig.small=TRUE>>=
library(edgeR)
ave.ab <- aveLogCPM(asDGEList(data))
hist(ave.ab, xlab="Average abundance", col="grey80", main="")
@

Any bin pair with an average abundance less than a specified threshold value will be discarded.
At the very least, the threshold should be chosen to filter out bin pairs with very low absolute counts.
This is because these bin pairs will never have sufficient evidence to reject the null hypothesis.
The discreteness of low counts will also reduce the accuracy of approximations that are used in normalization and statistical modelling.
The example below identifies those bin pairs where the overall NB mean across all samples is not less than 5.

<<>>=
count.keep <- ave.ab >= aveLogCPM(5, lib.size=mean(data$totals))
summary(count.keep)
@

The \Robject{count.keep} vector is then used to subset the \Rclass{InteractionSet} object.
The resulting \Robject{dummy} object will contain only the interesting bin pairs for downstream analysis.
The same procedure can be used for any logical or integer vector that specifies the bin pairs to be retained.

<<>>=
dummy <- data[count.keep,]
@

This count-based approach is fairly objective yet is still effective, i.e., it removes a large number of bin pairs that are likely to be uninteresting.
However, it will be less useful with greater sequencing depth where all bin pairs will have higher counts.
More sophisticated strategies can be implemented where the choice of threshold is motivated by some understanding of the Hi-C protocol.
These strategies will be described in the rest of this chapter.

\section{Directly removing low-abundance interactions}
\label{sec:direct}

\subsection{Computing the threshold from inter-chromosomal counts}
The simplest definition of an ``uninteresting'' interaction is that resulting from non-specific ligation. 
These are represented by low-abundance bin pairs where no underlying contact is present to drive ligation between the corresponding bins. 
Any changes in the counts for these bin pairs are uninteresting and are ignored.
In particular, the filter threshold can be defined by mandating some minimum fold change above the level of non-specific ligation.

The magnitude of non-specific ligation is empirically estimated by assuming that most inter-chromosomal contacts are not genuine. 
This is reasonable given that most chromosomes are arranged in self-interacting territories \cite{bickmore2013spatial}.
The median abundance across inter-chromosomal bin pairs is used as an estimate of the non-specific ligation rate. 
Here, filtering is performed to retain only those bin pairs with abundances that are at least 5-fold higher than this estimate.
This aims to remove the majority of uninteresting bin pairs.

<<>>=
direct <- filterDirect(data)
direct$threshold
direct.keep <- direct$abundances > log2(5) + direct$threshold
summary(direct.keep)
@

The \Rcode{direct.keep} vector can then be used to filter \Rcode{data}, as shown previously.
This approach is named here as ``direct'' filtering, as the average abundance for each bin pair is directly compared against a fixed threshold value.
In practice, the direct filter can be combined with \Rcode{count.keep} to ensure that the retained bin pairs have large absolute counts.

<<>>=
direct.keep2 <- direct.keep & count.keep
@

\subsection{Computing the threshold from larger bin pairs}
The procedure above assumes that no additional filtering has been performed during count loading with \Rfunction{squareCounts}, i.e., \Rcode{filter} is set to unity.
Any pre-filtering that removes low-abundance bin pairs will inflate the median and lead to overestimation of the filter threshold. 
However, it may not be practical to load counts without pre-filtering.
For example, at small bin sizes, too many non-empty bin pairs may be present to fit into memory.
Counts that are too small will also yield imprecise estimates of the background ligation rate.

To overcome this, counts are loaded for a larger bin size without pre-filtering.
This reduces memory usage as the interaction space is partitioned into fewer bin pairs.
The filter threshold is estimated from the inter-chromosomal interactions as previously described, exploiting the presence of large counts for large bin sizes to achieve precise estimates.
The estimated threshold for the larger bin pairs is then converted into a corresponding threshold for the original (smaller) bin pairs, after adjusting for the differences in bin sizes.

To illustrate, imagine that a bin size of 100 kbp is of interest.
We load the counts for the smaller bin pairs, using a reasonably large \Rcode{filter} to avoid excessive memory consumption.

<<>>=
new.bin.size <- 1e5
smaller.data <- squareCounts(input, mm.param, width=new.bin.size, filter=20)
@

Direct filtering of these smaller bin pairs is performed using \Rfunction{filterDirect}, using the larger (1 Mbp) bin pairs to estimate the filter threshold.
This is done by passing the unfiltered counts for the larger bin pairs in \Robject{data} as the \Rcode{reference} parameter.
The returned statistics can then be used to directly filter the smaller bin pairs.
Here, a minimum 5-fold change over the threshold is required for the retention of each 100 Mbp bin pair.

<<>>=
direct <- filterDirect(smaller.data, reference=data)
direct$threshold
small.keep <- direct$abundances > direct$threshold + log2(5)
summary(small.keep)
@

\section{Filtering as a function of interaction distance}

\subsection{Computing the threshold directly from the bin pair counts}
A more complex filter involves adjusting the threshold according to the distance between the bins in each bin pair. 
In typical Hi-C datasets, larger counts are observed at lower interaction distances.
This is probably driven by non-specific compaction of chromatin into a 3D ``globule'' conformation.
If such compaction is uninteresting, a concomitantly higher threshold is necessary to offset the increase in counts for these local interactions \cite{lin2012global}.

In the trended strategy, a trend is fitted to the abundances for all intra-chromosomal bin pairs using the log-distance as the covariate. 
The bin size is added to the distance as a prior, to avoid undefined values upon log-transformation when distances are zero.
The fitted value is then used as the threshold for each bin pair.
For inter-chromosomal bin pairs, the threshold is set to that from the direct filtering approach, for completeness.

<<>>=
trended <- filterTrended(data)
@

The effect of this strategy can be visualized by plotting the interaction distance against the normalized abundance.
A power-law relationship between distance and abundance is usually observed in Hi-C data \cite{lieberman2009comprehensive}. 
The average abundance (and thus, the filter threshold) decreases as the distance between the interacting loci increases.

<<avedistplot, fig.small=TRUE>>=
smoothScatter(trended$log.distance, trended$abundances, 
              xlab="Log-Distance", ylab="Normalized abundance")
o <- order(trended$log.distance)
lines(trended$log.distance[o], trended$threshold[o], col="red", lwd=2)
@

The assumption here is that the majority of interactions are generated by non-specific packaging of the linear genome.
Each bin pair is only retained if its abundance is greater than the corresponding fitted value at that distance, i.e., above that expected from compaction. 
This favors selection of longer-range interactions, compared to the direct filter.

<<>>=
trend.keep <- trended$abundances > trended$threshold
summary(trend.keep)
@

Of course, the threshold can also be defined at some minimum fold change above the fitted value.
This effectively increases the stringency of the filter.
The example below retains bin pairs with abundances that are two-fold higher than the expected compaction intensity.

<<>>=
trend.keep2 <- trended$abundances > trended$threshold + log2(2)
summary(trend.keep2)
@

The trended filter can also be combined with \Rcode{count.keep} to ensure that the absolute counts are large.
This is particularly important at large distances, where the drop in the threshold may lead to the inappropriate retention of very low abundance bins.

<<>>=
trend.keep3 <- trend.keep & count.keep
@

The distance between bins can also be obtained directly with the \Rfunction{pairdist} function.
Distances for inter-chromosomal bin pairs are marked as \Rcode{NA}.
These tend to dominate the output as they constitute most of the interaction space.
Of course, the majority of these bin pairs will have low counts due to the sparseness of the data.

\subsection{Using larger bin pairs to save memory}
The procedure above assumes that no pre-filtering was performed on the counts for small bin pairs.
In situations with limited memory, the counts for larger bin pairs can again be used to define the thresholds.
The trend is fitted to the abundances and distances of the larger bin pairs, and interpolation (or extrapolation) is performed to obtain fitted values at the distances of the smaller bin pairs.
The interpolated trend can then be applied to filter the smaller bin pairs based on their abundances.
This entire process is done automatically within \Rfunction{filterTrended} when a \Rcode{reference} argument is supplied, as shown below.

<<>>=
trended <- filterTrended(smaller.data, reference=data)
summary(trended$abundances > trended$threshold)
@

Another advantage of this approach is that threshold estimation is more precise with larger counts.
Thus, even if there is enough memory for loading smaller bin pairs with \Rcode{filter=1}, larger bin sizes may still be preferred for computing the filter threshold.
Otherwise, median calculation or trend fitting will be confounded by discreteness at low counts.

\section{Peak-calling in the interaction space}

\subsection{Defining enrichment values for each bin pair}
Peaks in the interaction space are bin pairs that have substantially more reads than their neighbors.
The \Rfunction{enrichedPairs} function counts the number of read pairs in the ``neighborhood'' of each bin pair.
For a bin pair $x$, the function considers several neighborhood regions including bin pairs above or below $x$; in pairs to the left and right of $x$; bin pairs in a box centered at $x$; and for intra-chromosomal bin pairs, the quadrant closest to the diagonal \cite{rao2014kilobase}.
The size of each neighborhood is determined by \Rcode{flank}, defined in terms of bins.
In this case, each bin is around 1 Mbp so the flank width has an actual size of \textasciitilde{}5 Mbp.

<<>>=
flank.width <- 5
en.data <- enrichedPairs(data, flank=flank.width)
en.data
@

The \Rfunction{filterPeaks} function uses the neighborhood counts to compute an enrichment value for each bin pair.
The average abundance across libraries for each neighborhood region is computed, and the region with the largest abundance is chosen.
The enrichment value is defined as the difference between the abundance of the bin pair and that of its chosen neighborhood region (adjusted for the number of bin pairs in the latter).
Bin pairs with high enrichment values are putative peaks that should be retained.
This extends the concept of peak calling on the linear genome (e.g., in ChIP-seq) to the 2-dimensional interaction space \cite{rao2014kilobase}.
Neighborhood regions are defined to capture high-intensity structural features like looping domains, TADs or banding patterns.
This ensures that the structural features themselves do not drive high enrichment values.
Rather, to be called a peak, a bin pair in a structural feature must be enriched relative to other bin pairs in the same feature.

<<>>=
enrichments <- filterPeaks(en.data, get.enrich=TRUE)
summary(enrichments)   
@

In this example, an enrichment value above 0.5 is required for retention.
This is a modest threshold that errs on the side of caution, as weak peaks may still be DIs and should be retained for testing.
In practice, filtering of enrichment values should be combined with filtering on absolute abundances.
This eliminates low-abundance bin pairs with high enrichment values, e.g., when the neighborhood is empty.
Near-diagonal elements should also be removed, as these tend to have high enrichment scores without actually being peaks.
All of these steps can be conveniently applied through the \Rfunction{filterPeaks} wrapper function.

<<>>=
peak.keep <- filterPeaks(data, enrichments, min.enrich=0.5, 
                         min.count=5, min.diag=2L)
sum(peak.keep)
@

\subsection{Examining some peak-calling diagnostics}
This filtering strategy can be evaluated by examining the interaction space around each putative peak (see Chapter~\ref{chap:plaid}).
Briefly, the color intensity of each ``pixel'' is proportional to the number of read pairs between the corresponding loci on each axis.
Red boxes mark the bin pairs retained by filtering, where each side represents a bin.
In both examples, the bin pairs correspond nicely to punctate high-intensity contacts in the interaction space.

<<peakdiagplot,echo=FALSE,fig.wide=TRUE,fig.asp=0.6>>=
#keep2 <- peak.keep & !is.na(dist)
#which(keep2)[order(enrichments[keep2], decreasing=TRUE)][1:100]
#plotPlaid(input[3], param=mm.param, width=2e5, cap=20, 
#	first=resize(anchors(data[chosen], type="first"), fix="center", width=10e6),
#	second=resize(anchors(data[chosen], type="second"), fix="center", width=10e6))
a1.list <- GRanges(c("chr11", "chr4"), IRanges(
			c(80999224, 153002112),
			c(81998782, 154006051)))
a2.list <- GRanges(c("chr11", "chr4"), IRanges(
			c(66001477, 130998845),
			c(67011476, 131997373)))
capped <- c(50, 50)
par(mfrow=c(1,2))
for (chosen in seq_along(a1.list)) { 
	if (!suppressWarnings(any(overlapsAny(anchors(data[peak.keep], type="first"), a1.list[chosen], type="equal") & 
			overlapsAny(anchors(data[peak.keep], type="second"), a2.list[chosen], type="equal")))) { stop("can't find example for peak-calling") }
	plotPlaid(input[3], param=mm.param, width=2e5, max.count=capped[chosen], 
			first.region=resize(a1.list[chosen], fix="center", width=10e6),
			second.region=resize(a2.list[chosen], fix="center", width=10e6),
			main=paste("Example", chosen))
	box()
	rect(start(a1.list[chosen]), start(a2.list[chosen]),
		end(a1.list[chosen]), end(a2.list[chosen]), border="red")
}
@

Another diagnostic is based on the sparsity of putative peaks.
Bin pairs nested within larger ``parent'' bin pairs are identified using the \Rfunction{boxPairs} function.
Most of these parents should contain low numbers of nested bin pairs.
This is because each peak, by definition, should be isolated in its neighborhood.
The example below considers parent bin pairs of size equal to the neighborhood used to compute the enrichment values,
and calculates the percentage of parents that contain a given number (from 1 to 10, or higher) of nested bin pairs.

<<>>=
neighborhood <- (2*flank.width + 1) * metadata(data)$width
boxed <- boxPairs(data[peak.keep], reference=neighborhood)
out <- tabulate(tabulate(boxed$indices[[1]]))
out <- c(out[1:10], sum(out[11:length(out)])) # sum all >10
setNames(round(out/sum(out)*100, 1), c(1:10, ">10"))
@

Most parents should contain few bin pairs ($\le$ 10, as a rule of thumb), as shown above.
In the second example below, a larger proportion of parents with many nested bin pairs is observed.
This suggests that more aggressive filtering on the enrichment score is required.

<<>>=
peak.keep2 <- filterPeaks(data, enrichments, min.enrich=0, 
                          min.count=5, min.diag=2L)
boxed <- boxPairs(data[peak.keep2], reference=neighborhood)
out <- tabulate(tabulate(boxed$indices[[1]]))
out <- c(out[1:10], sum(out[11:length(out)])) # sum all >10
setNames(round(out/sum(out)*100, 1), c(1:10, ">10"))
@

Obviously, peak calling assumes that changes in intensity of broad interactions are not interesting.
This may not be appropriate for every study.
Indeed, the definition of ``broad'' depends on the bin size, whereby a peak called at large bin sizes may be discarded at smaller sizes.
Users should try out the other, more general filters before proceeding to peak calling, to check that potentially interesting differences are not discarded by the latter.

\subsection{Efficient peak calling at high resolution}
The \Rfunction{enrichedPairs} function requires that all bin pairs are loaded into memory, i.e., with \Rcode{filter=1} during \Rfunction{squareCounts}.
This may not be feasible for high resolution analyses with small bin sizes, where there is a large number of bin pairs with low counts.
In such cases, it may be more practical to use the \Rfunction{neighborCounts} function.
This re-counts the read pairs for each bin pair, storing only a limited portion of the interaction space to compute the neighbourhood counts.
The need to load all non-empty bin pairs is avoided.
Only bin pairs with count sums greater than or equal to \Rcode{filter} are returned, along with neighborhood counts that can be used in \Rfunction{filterPeaks}. 
This reduces memory usage at the cost of some speed.

<<>>=
en.data <- neighborCounts(input, param, width=1e5, filter=20, flank=5)
en.data
enrichments <- filterPeaks(en.data, get.enrich=TRUE)
summary(enrichments)
@

\section{Filtering for pre-specified regions}
Filtering for pairs of arbitrary regions is complicated by the potential irregularity of the regions.
In particular, there is no guarantee that the supplied regions will cover the entirety of the interaction space.
Filter thresholds may not be estimated accurately if the covered area is not representative of the rest of the space.
At worst, genuine interactions between all specified regions would preclude direct or trended filtering, which assume that most interactions are uninteresting.
That said, threshold estimation from a subset of the interaction space may actually be desirable in some cases, e.g., using only the captured regions to compute a relevant estimate of the non-specific ligation rate in a Capture-C experiment.

Another consideration is that different regions will have different widths.
The resulting areas used for read pair counting will probably be different between region pairs, such that some will have systematically higher counts than others.
Whether or not this is a problem depends on the user's perspective.
The position of this guide is that it would be inappropriate to penalize larger areas for having larger counts. 
As long as there are sufficient counts for an analysis, the size of the area involved should be irrelevant.
Thus, use of \Rfunction{aveLogCPM} alone is recommended for calculation of the average abundance when filtering region pairs.

<<>>=
summary(aveLogCPM(asDGEList(redata)))
@

\section{Filtering out diagonal elements}
Incomplete removal of artifacts (e.g., dangling ends, self-circles) will generally manifest as counts for diagonal bin pairs, i.e., pairs of the same bin.
If artifact generation and/or removal is not consistent between libraries, the behavior of the diagonal elements will differ markedly from the other bin pairs.
This can be diagnosed using MA plots between libraries, where the diagonal elements show up as a distinct mass that is unconnected to the other points.
As such, they cannot be smoothly normalized and should be removed prior to further analysis or analyzed separately.
This is done by generating \Rcode{ndiag.keep} for filtering, as shown below.

<<>>=
ndiag.keep <- filterDiag(data)
summary(ndiag.keep)
@

Removal of diagonal elements is also motivated by the difficulty of intepreting interactions within a bin.
For small bin sizes, read pairs corresponding to internal interactions are indistinguishable from those driven by non-specific packaging of a locus.
The latter is mostly uninteresting and can be avoided by discarding the bin pairs on the diagonal.
For large bin sizes, internal interactions may be more interesting -- these can be recovered by repeating the analysis with smaller bins.
Short-range interactions will then be represented by the off-diagonal pairs, while artifacts will still be restricted to diagonal bin pairs, so long as the bins are larger than \textasciitilde{}25 kbp (see Section~\ref{sec:strorient}).
Count sizes will also decrease but this should not be a major problem as counts should be inherently larger at low interaction distances. 

\section{Summary of the filtering strategies}
Each filtering strategy can be tuned by increasing or decreasing the minimum fold change required for retention.
This can be driven by the biological knowledge available to the user, based on features that are biologically interesting.
Even when such knowledge is unavailable, the filters are still useful as they can guide the user towards a sensible interpretation of the filter threshold.
This would not be possible if the threshold value was chosen arbitrarily.

The choice of filtering method depends on the features that are most likely to be of interest in each analysis.
In general, less aggressive filtering should be performed if these features are not well defined.
Here, a tentative recommendation is provided for direct filtering as it is reasonable to discard non-specific ligation events.
On a practical level, the simplicity of the direct approach is attractive and will be used throughout the rest of the guide.

<<>>=
original <- data
data <- data[direct.keep2,] 
@

The filtered results are assigned back into \Robject{data}.
This is mostly for consistency, so that all operations are performed on \Robject{data} in the rest of this guide.
The original unfiltered data is retained for later use in Section~\ref{sec:itercor}, but can otherwise be ignored.
Direct filtering is also performed for the smaller bin pairs, and the results stored in \Robject{smaller.data}.

<<>>=
smaller.data <- smaller.data[small.keep,]
@

\chapter{Normalization strategies for Hi-C data}

\begin{combox}
Here, we're using the \Robject{data} object that was filtered in the previous chapter.
We'll also need the \Robject{original} object, as well as \Robject{margin.data} from Section~\ref{sec:marginal}.
Another human dataset will be loaded here, so \Robject{hs.frag} from Chapter~\ref{chap:prep} will be required.
\end{combox}

\section{Normalizing with scaling methods}
\label{sec:simplenorm}

The simplest approach is to use scaling methods such as library size normalization.
This accounts for differences in library preparation efficiency and sequencing depth between samples.
In \Biocpkg{edgeR}, scaling is performed using the effective library size, defined as the product of the library size and the normalization factor for each sample (see the \Biocpkg{edgeR} user's guide).
Library size normalization is thus achieved by setting all normalization factors to unity.

<<>>=
data.lib <- data # Making a copy to avoid side-effects.
data.lib$norm.factors <- 1
@

Alternatively, TMM normalization \cite{oshlack2010tmm} can be applied on the counts for the bin pairs.
This accounts for composition biases by assuming that most interactions do not change in intensity between conditions.
Here, the normalization factors differ from unity as the library size alone is not sufficient to describe the bias in each sample.

<<>>=
data.tmm <- normOffsets(data, se.out=TRUE)
data.tmm$norm.factors 
@

In practice, scaling methods are usually too simplistic.
More sophisticated approaches are necessary to handle the complex biases observed in real Hi-C data.

\section{Removing trended biases between libraries}

\subsection{Using non-linear normalization}
Trended biases may be observed in Hi-C data, caused by uncontrolled differences in sample preparation in a complex protocol.
Changes in cross-linking efficiency or ligation specificity can lead to a systematic redistribution of read pairs throughout the interaction space. 
For example, reduced specificity may result in more counts for weak non-specific interactions, and fewer counts for strong genuine interactions.
Such biases may manifest as an abundance-dependent trend in a MA plot between libraries.
This is especially true for a trend observed between replicates, which must be technical in origin.
The code below generates a MA plot comparing one library from each group, which can be repeated for each pair of libraries.

<<echo=FALSE,eval=FALSE,label=maunnorm>>=
ab <- aveLogCPM(asDGEList(data))
o <- order(ab)
adj.counts <- cpm(asDGEList(data), log=TRUE)
@

<<echo=FALSE,eval=FALSE,label=makema>>=
mval <- adj.counts[,3]-adj.counts[,2]
smoothScatter(ab, mval, xlab="A", ylab="M", main="KO (1) vs. Flox (2)")
fit <- loessFit(x=ab, y=mval)
lines(ab[o], fit$fitted[o], col="red")
@

<<unnormmaplot, fig.small=TRUE>>=
<<maunnorm>>
<<makema>>
@

Trended biases are problematic as they can inflate the variance estimates or fold-changes for some bin pairs. 
They must be eliminated with non-linear normalization prior to further analysis.
We use LOESS normalization with the \Rfunction{normOffsets} function from the \Biocpkg{csaw} package \cite{lun2016csaw}, adapted from existing non-linear methods to handle discrete count data. 

<<>>=
data <- normOffsets(data, type="loess", se.out=TRUE)
@

This function computes an offset term for each bin pair in each library.
When fitting a generalized linear model (GLM) to the counts for a particular bin pair, a large offset for a library is equivalent to downscaling the corresponding count relative to the counts of other libraries. 
The matrix of offsets has the same dimensions as the count matrix and is stored as an element of the \Rcode{assays} slot of the \Rclass{InteractionSet} object.
(Specifying \Rcode{se.out=FALSE} returns the offset matrix directly, instead of a modified \Rclass{InteractionSet} containing the offsets.)

<<>>=
nb.off <- assay(data, "offset")
head(nb.off)
@

In most cases, the above code is sufficient for normalization.
However, as previously mentioned, bin pairs near the diagonal may exhibit irregular behavior relative to the rest of the interaction space.
This manifests in the previous MA plot as a distinct mass of points at high abundances, preventing smooth normalization of the trended bias.
To overcome this, the near-diagonal bin pairs can be normalized separately from the others.

<<>>=
neardiag <- filterDiag(data, by.dist=1.5e6)
nb.off <- matrix(0, nrow=nrow(data), ncol=ncol(data))
nb.off[neardiag] <- normOffsets(data[neardiag,], type="loess", se.out=FALSE)
nb.off[!neardiag] <- normOffsets(data[!neardiag,], type="loess", se.out=FALSE)
assay(data, "offset") <- nb.off # updating the offset matrix
@

We examine the MA plot after adjusting the log-counts with the offsets.
Most of the trend is removed, indicating that non-linear normalization was successful.

<<echo=FALSE,eval=FALSE,label=manorm>>=
adj.counts <- log2(assay(data) + 0.5) - nb.off/log(2)
@

<<normmaplot, fig.small=TRUE>>=
<<manorm>>
<<makema>>
@

\subsection{Requirements of non-linear normalization}
Filtering on average abundance prior to normalization is strongly recommended.
This removes low counts and avoids problems with discreteness during trend fitting.
Filtering also improves the sensitivity of span-based fitting algorithms like LOESS at higher abundances.
Otherwise, the fitted trend will be dominated by the majority of low-abundance bin pairs.

Non-linear normalization also assumes that the majority of bin pairs at each abundance do not represent differential interactions.
Any systematic differences between libraries are assumed to be technical in origin and are removed.
If this assumption does not hold, genuine differences may be lost upon normalization.
For example, a global increase in the compaction of the genome would manifest as an upward trend in the MA plot, as low-abundance distal interactions are weakened in favor of high-abundance short-range interactions.
In such cases, use of library size normalization in Section~\ref{sec:simplenorm} may be more appropriate.

\section{Iterative correction of interaction intensities}
\label{sec:itercor}
While this is not the intended purpose of the \Biocpkg{diffHic} package, a method is also provided for the removal of biases between genomic regions.
The \Rfunction{correctedContact} function performs iterative correction of the interaction space \cite{imakaev2012iterative} with some modifications.
Namely, if multiple libraries are used to present, correction is performed using the overall NB mean for each bin pair, rather than on the counts themselves.
Winsorizing through \Rcode{winsor.high} is also performed to mitigate the impact of high-abundance bin pairs.

<<>>=
corrected <- correctedContact(original, winsor.high=0.02, ignore.low=0.02)
head(corrected$truth)
@

The returned \Rcode{truth} contains the ``true'' contact probability for each bin pair in \Robject{data}.
This is designed to account for differences in sequencibility, mappability, restriction site frequency, etc. between bins.
Comparisons can then be directly performed between the contact probabilities of different bin pairs.
Some \Rcode{NA} values will be present due to the removal of low-abundance bins that do not exhibit stable behavior during correction.
The convergence of the correction procedure can then be checked by examining the maximum fold change to the truth at each iteration.
This should approach unity, i.e., no further change.

<<>>=
corrected$max
@

Note that \Rcode{original} is used as the input to \Rfunction{correctedContact}.
No filtering should be performed prior to iterative correction.
All non-empty bin pairs are needed as information is collected across the entire interaction space.
The contribution of many bin pairs with low counts might end up being as substantial as that of a few bin pairs with large counts.

In addition, normalization from iterative correction can be fed naturally into the DI analysis via the GLM machinery in \Biocpkg{edgeR}.
The log-transformed product of the estimated genomic biases for both bins in each bin pair can be used as the offset for that bin pair.
This is computed separately for each library by setting \Rcode{average=FALSE}, to correct for sample-specific genomic biases.
\Rcode{NA} offsets will be obtained for some bin pairs, but these should not be too problematic given that the affected bin pairs will generally have low counts (as at least one interacting partner will be of low abundance) and be filtered out during the analysis proper.

<<>>=
corrected <- correctedContact(original, average=FALSE)
anchor1.bias <- corrected$bias[anchors(original, type="first", id=TRUE),]
anchor2.bias <- corrected$bias[anchors(original, type="second", id=TRUE),]
iter.off <- log(anchor1.bias * anchor2.bias)
@

Of course, iterative correction only removes biases between different bins.
It is not guaranteed to remove (trended) biases between libraries.
For example, two replicates could have the same genomic biases but a different distribution of read pairs in the interaction space, e.g., due to differences in ligation specificity.
The latter would result in a trended bias, but iterative correction would have no effect due to the identical genomic biases.
In short, normalization within a library is a different problem from normalization between libraries.

\section{Accounting for copy number variations}
\label{sec:copy}

\subsection{Eliminating CNVs with multi-dimensional smoothing}
Copy number variations (CNVs) in the interacting regions will also affect the interaction intensity. 
These CNV-driven differences in intensities are generally uninteresting and must be removed to avoid spurious detection of DIs.
Normalization of CNV-based biases is achieved using the \Rfunction{normalizeCNV} function, which performs  multi-dimensional smoothing across several covariates with the \CRANpkg{locfit} package \cite{loader1999local}.
It requires both bin pair and marginal counts, obtained with the same \Rcode{width} and \Rcode{param} in calls to \Rfunction{squareCounts} and \Rfunction{marginCounts}.

<<>>=
cnv.offs <- normalizeCNV(data, margin.data)
head(cnv.offs)
@

Three covariates are defined for each bin pair.
For a pair of libraries, the ratio of marginal counts for each bin can be used as a proxy for the relative CNV between libraries in that bin.
Each bin pair will be associated with two of these marginal log-ratios to use as covariates.
The third covariate for each bin pair is that of the average abundance across all libraries.
This will account for any abundance-dependent trends in the biases.

The response for each bin pair is defined as the log-ratio of interaction counts between a pair of  libraries.
A locally weighted surface is fitted to the response against all three covariates for all bin pairs.
At any combination of covariate values, most bin pairs are assumed to represent non-differential interactions.
Any systematic differences between libraries are attributed to CNV-driven (or trended) biases and are removed.
Specifically, GLM offsets are returned and can be supplied to the statistical analysis to eliminate the bias.

As with trended biases, filtering by average abundance is strongly recommended prior to running \Rfunction{normalizeCNV}.
This reduces the computational work required for multi-dimensional smoothing.
Discreteness and domination of the fit by low-abundance bin pairs is also avoided.

\subsection{Visualizing the effect of CNV removal}
To demonstrate, the RWPE1 dataset \cite{rickman2012oncogene} is used here as it contains more CNVs.
Interaction counts are loaded for each bin pair, and marginal counts are loaded for each bin.
Some filtering is performed to eliminate low-abundance bin pairs, as previously described.

<<>>=
count.files <- c("merged_erg.h5", "merged_gfp.h5")
rick.data <- squareCounts(count.files, hs.param, width=1e6)
rick.marg <- marginCounts(count.files, hs.param, width=1e6)
rick.data <- rick.data[aveLogCPM(asDGEList(rick.data)) > 0,]
@

Our aim is to plot the log-ratio of the counts for each bin pair between the first two libraries, as a function of the log-ratio of the marginal counts of the corresponding bins.
The code below computes the marginal log-ratios after matching the bins in \Robject{rick.marg} to the bin pairs in \Robject{rick.data} with \Rfunction{matchMargins}.
As two marginal log-ratios are present for each bin pair, these are added together to simplify visualization.

<<>>=
m.adjc <- cpm(asDGEList(rick.marg), log=TRUE)
margin.lr <- m.adjc[,1] - m.adjc[,2]
matched <- matchMargins(rick.data, rick.marg)
margin.lr <- margin.lr[matched$anchor1] + margin.lr[matched$anchor2]
@

We generate plots before and after subtracting the offsets computed by \Rfunction{normalizeCNV}.
Before normalization, we observe that a decrease in copy number results in a decrease in interaction intensity.
This is expected given that fewer copies should result in fewer chances for interaction.
After normalization, the trend in the interaction log-ratios is removed.
This indicates that the CNV effect has been eliminated and can be ignored in downstream analyses.
Note that increasing \Rcode{maxk} may be necessary to obtain accurate results in the internal call to \CRANpkg{locfit}.

<<cnvplotter, fig.wide=TRUE, fig.asp=0.5>>=
before <- cpm(asDGEList(rick.data), log=TRUE)
cnv.offs <- normalizeCNV(rick.data, rick.marg, maxk=1000)
after <- log2(assay(rick.data)+0.5) - cnv.offs/log(2)
par(mfrow=c(1,2), cex.axis=1.2, cex.lab=1.4)
smoothScatter(margin.lr, before[,1]-before[,2], ylim=c(-4, 4), main="Before",
    xlab="Sum of marginal log-ratios", ylab="Interaction log-ratio")
smoothScatter(margin.lr, after[,1]-after[,2], ylim=c(-4, 4), main="After",
    xlab="Sum of marginal log-ratios", ylab="Interaction log-ratio")
@

\subsection{Using alternative copy number definitions}
The method above uses the marginal counts as a proxy for genomic coverage \cite{imakaev2012iterative}.
Changes in the marginal counts can be used as a proxy for change in the copy number between samples.     
In most cases, this is an effective and sufficient strategy as exact quantification of the copy number difference is not required.
However, external information can also be used, based on copy number arrays or genomic sequencing data for each sample.
The (relative) copy number for each bin/region in each sample can be stored in a \Rclass{RangedSummarizedExperiment} and passed to \Rfunction{normalizeCNV} as shown above.
This may be more accurate as it avoids potential problems from interaction-specific effects on the marginal counts.

\chapter{Modelling biological variability}

\begin{combox}
In this chapter, the \Robject{data} object is again required.
The computed offsets in \Robject{nb.off} will also be used from the last chapter.
Methods from the \Biocpkg{edgeR} package should already be loaded into the workspace, but if they aren't, then \Rcode{library(edgeR)} will do the job.
\end{combox}

\section{Overview}
The differential analysis in \Biocpkg{diffHic} is based on the statistical framework in the \Biocpkg{edgeR} package \cite{edgeR}.
This models the counts for each bin pair with NB distributions.
The NB model is useful as it can naturally accommodate low, discrete counts.
It can also consider extra-Poisson variability between biological replicates of the same condition.
Here, biological replicates refer to Hi-C libraries prepared from independent biological samples.

The magnitude of biological variability is empirically determined from these biological replicates.
In \Biocpkg{edgeR}, variability is modelled by estimating the dispersion parameter of the NB distribution.
This is used during testing to reduce the significance of any detected differences when the counts are highly variable.
Similarly, estimation of the quasi-likelihood (QL) dispersion can be performed to model variability of the dispersions \cite{lund2012ql}. 

Dispersion estimation requires the fitting of a GLM to the counts for each bin pair \cite{mccarthy2012glm}. 
To do so, a design matrix must be specified to describe the experimental setup.
For the neural stem cell dataset, a simple one-way layout is sufficient.
The code below specifies two groups of two replicates, where each group corresponds to a genotype.
The aim is to compute the dispersion from the variability in counts within each group.

<<>>=
design <- model.matrix(~factor(c("flox", "flox", "ko", "ko")))
colnames(design) <- c("Intercept", "KO")
design
@

Obviously, more complex designs can be used if necessary.
This includes designs with blocking factors for batch effects, pairing between samples or multiple groups.

The \Rclass{InteractionSet} object also needs to be converted into a \Rclass{DGEList} object for analysis with \Biocpkg{edgeR}.
If the normalization factors or offsets are already stored in \Robject{data}, as previously described, then they will be automatically extracted by the \Rfunction{asDGEList} function and stored in the output \Rclass{DGEList}.
Alternatively, they can be passed explicitly to \Rfunction{asDGEList}.
Note that the offset matrix will take precedence over scaling factors in \Biocpkg{edgeR} analyses.

<<>>=
y <- asDGEList(data)
@

\section{Estimating the NB dispersion}
Estimation of the NB dispersion is performed by maximizing the Cox-Reid adjusted profile likelihood (APL) \cite{mccarthy2012glm} for each bin pair.
Of course, when replication is limited, there is not enough information per bin pair to estimate the dispersion.
This is overcome by computing and sharing APLs across many bin pairs to stablize the estimates.

<<>>=
y <- estimateDisp(y, design)
y$common.dispersion
@

A more sophisticated strategy is also used whereby an abundance-dependent trend is fitted to the APLs.
This should manifest as a smooth trend in the NB dispersion estimates with respect to the average abundances of all bin pairs.
The aim is to improve modelling accuracy by empirically modelling any non-NB mean-variance relationships.

<<bcvplot, fig.small=TRUE>>=
plotBCV(y)
@

In most cases, the relationship should be monotonic decreasing as the counts become more precise with increasing size.
Minor deviations are probably due to the imperfect nature of non-linear normalization.
Major increases are indicative of batch effects.
For example, a cluster of outliers indicates that there may be copy number changes between replicates.

\section{Estimating the QL dispersion}
The QL dispersion for each bin pair is estimated from the deviance of the fitted GLM.
This may seem superfluous given that the NB dispersion already accounts for biological variability.
However, it is mathematically easier to model variability in the QL dispersions across bin pairs, compared to variability in the NB dispersions.
GLM fitting and estimation of the QL dispersion for each bin pair are performed with the \Rfunction{glmQLFit} function.

<<>>=
fit <- glmQLFit(y, design, robust=TRUE)
@

Again, there is not enough information for each bin pair for precise estimation.
Instead, information is shared between bin pairs using an empirical Bayes (EB) approach.
Per-bin-pair QL estimates are shrunk towards a common trend across all bin pairs.
This stabilizes the QL dispersion estimates and improves precision for downstream applications.

<<qlplot, fig.small=TRUE>>=
plotQLDisp(fit)
@

The extent of the EB shrinkage is determined by the variability of the dispersions.
If the true dispersions are highly variable, shrinkage to a common value would be inappropriate.
On the other hand, more shrinkage can be performed to increase precision (and thus detection power) if the true dispersions are not variable.
This is quantified as the prior degrees of freedom, for which smaller values correspond to more variability and less shrinkage.

<<>>=
summary(fit$df.prior)
@

It is important to use the \Rcode{robust=TRUE} argument in \Rfunction{glmQLFit} \cite{phipson2016robust}.
This protects the EB shrinkage against large positive outliers corresponding to highly variable counts.
It also protects against large negative outliers.
These are formed from near-zero deviances when counts are identical, and are not uncommon when counts are low.
Both types of outliers inflate the apparent variability and decrease the estimated prior degrees of freedom.

\section{Further information}
More details on the statistical methods in \Biocpkg{edgeR} can be found, unsurprisingly, in the \Biocpkg{edgeR} user's guide.
Of course, \Biocpkg{diffHic} is compatible with any statistical framework that accepts a count matrix and a matrix of log-link GLM offsets.
Advanced users may wish to use methods from other packages.
This author prefers \Biocpkg{edgeR} as it works quite well for routine analyses of Hi-C data.
He also has a chance of being cited when \Biocpkg{edgeR} is involved \cite{chen2014differential}.

\chapter{Testing for significant interactions}

\begin{combox}
This chapter brings everything together.
We need the \Robject{fit} object from the last chapter, along with the \Robject{data} object (as usual).
We also require the \Robject{smaller.data} object from Chapter~\ref{chap:filter}.
The \Robject{bin.size} and \Robject{mm.param} objects are required from Chapter~\ref{chap:counting}.
\end{combox}

\section{Using the quasi-likelihood F-test}
The \Rfunction{glmQLFTest} function performs a QL F-test for each bin pair to identify significant differences.
Users should ensure that the correct contrast is performed by explicitly specifying the \Rcode{coef} or \Rcode{contrast} arguments.
In this case, the coefficient of interest refers to the change in the KO counts over the WT counts.
The null hypothesis for each bin pair is that the coefficient is equal to zero, i.e., there is no change between the WT and KO groups.

<<>>=
result <- glmQLFTest(fit, coef=2)
topTags(result)
@

More savvy users might wonder why the likelihood ratio test (LRT) was not used here.
Indeed, the LRT is the more obvious test for any inferences involving GLMs. 
However, the QL F-test is preferred as it accounts for the variability and uncertainty of the QL dispersion estimates \cite{lund2012ql}. 
This means that it can maintain more accurate control of the type I error rate compared to the LRT in the presence of limited replication.

It is also convenient to store the significance statistics in the \Rcode{rowData} of the \Rclass{InteractionSet} object.
This ensures that any operations like subsetting are applied on both the genomic coordinates and statistics simultaneously, which simplifies data management.

<<>>=
rowData(data) <- cbind(rowData(data), result$table)
@

\section{Multiplicity correction and the FDR}

\subsection{Overview}
In a \Biocpkg{diffHic} analysis, many bin pairs are tested for significant differences across the interaction space.
Correction for multiple testing is necessary to avoid excessive detection of spurious differences.
For genome-wide analyses, this correction is typically performed by controlling the false discovery rate (FDR) with the Benjamini-Hochberg (BH) method \cite{benjamini1995fdr}. 
This provides a suitable comprimise between specificity and sensitivity.
In contrast, traditional methods of correction (e.g., Bonferroni) are often too conservative.

Interpretation of the differential interactions (DIs) is complicated by spatial dependencies between adjacent bin pairs.
If many adjacent bin pairs were separately reported as DIs, this would complicate the interpretation of the results by introducing unnecessary redundancy.
To mitigate this effect, \Biocpkg{diffHic} provides several methods for clustering bin pairs and controlling the relevant FDR.
Each of these methods is described below, along with its relative strengths and weaknesses.
For routine analyses, the approach based on clustering significiant bin pairs is recommended -- see Section~\ref{sec:sigclust} for details.

\subsection{Direct application of the BH method}
The simplest approach is to directly apply the BH method to the set of $p$-values for all tested bin pairs.
The FDR here refers to the proportion of detected bin pairs that are false positives.
DIs are defined as those bin pairs detected at a given FDR threshold, e.g., 5\%.
Their coordinates are saved to file along with the bin pair-based test statistics.
Resorting on $p$-value is performed to prioritize strong DIs, thus simplifying inspection of the results.

<<>>=
adj.p <- p.adjust(result$table$PValue, method="BH")
sum(adj.p <= 0.05)
useful.cols <- as.vector(outer(c("seqnames", "start", "end"), 1:2, paste0))
inter.frame <- as.data.frame(interactions(data))[,useful.cols]
results.r <- data.frame(inter.frame, result$table, FDR=adj.p)
o.r <- order(results.r$PValue)
write.table(results.r[o.r,], file="binpairs.tsv", sep="\t", 
            quote=FALSE, row.names=FALSE)
@

This approach is best suited for analyses using large bins, i.e., greater than or equal to 1 Mbp.
Here, adjacent bin pairs need not be considered redundant as they cover distinct areas of the interaction space.
Thus, no clustering is required to simplify interpretation.

\subsection{Clustering to reduce redundancy in the results}
Adjacent bin pairs can be aggregated into larger clusters to reduce redundancy in the results.
This is especially useful at smaller bin sizes where multiple bin pairs may overlap a single underlying interaction.
Each cluster of bin pairs will then represent the underlying interaction.
This is demonstrated below with the 100 kbp bin pairs.
The \Rfunction{clusterPairs} function puts two bin pairs in the same cluster if they are no more than \Rcode{tol} bp apart in any dimension.

<<>>=
clustered.small <- clusterPairs(smaller.data, tol=1, upper=1e6)
@

In practice, this approach requires aggressive filtering to avoid chaining effects.
Otherwise, clustering will be confounded in high-density areas of the interaction space, e.g., at short distances or in TADs.
This leads to the formation of very large clusters that are difficult to interpret.
Some protection is provided by specifying the maximum dimensions of each cluster in \Rcode{upper}.
This will break up very large clusters, albeit in a somewhat arbitrary manner.

There is also no guarantee that the cluster will form a regular shape in the interaction space.
In \Rfunction{clusterPairs}, an approximate solution is used whereby the minimum bounding box for each cluster is reported.
This refers to the smallest rectangle in the interaction space that contains all bin pairs in the cluster.
The coordinates of this rectangle can be easily recorded, whereas it is more difficult to store the detailed shape of the cluster.
Identification of the top-ranking bin pair within each cluster may also be desirable (see Section~\ref{sec:nesting}).

<<>>=
length(clustered.small$interactions)
head(clustered.small$interactions)
@

The FDR across clusters is not the same as that across bin pairs.
The former is more useful as the results are usually interpreted in terms of clusters, where each cluster corresponds to one interaction event.
However, the FDR across bin pairs is easier to control by applying the BH method directly to the bin pair $p$-values.
Treating the FDR across bin pairs as that across clusters is usually inappropriate and results in loss of FDR control \cite{lun2014denovo}.
The following text will describe some strategies for controlling the cluster-level FDR.

\subsection{Independent clustering of adjacent bin pairs}
\label{sec:independentcluster}

The basic idea of independent clustering is that it is done independently of the differential test statistics.
This means that other metrics must be used to identify relevant bin pairs.
If there are pre-defined areas of interest in the interaction space
    (e.g., based on existing interactions from earlier experiments or when studying interactions between annotated features),
    identification of these bin pairs can be performed with methods like \Rfunction{findOverlaps} from the \Biocpkg{InteractionSet} package.
All bin pairs overlapping a pre-defined area can then be defined as a single cluster.
Otherwise, filtering on the average abundance will implicitly restrict the clustering to high-abundance bin pairs.
This may be satisfactory when the interaction space is sparse such that the clusters will be small enough to be interpretable.
To demonstrate, a quick-and-dirty differential analysis is performed using counts for the 100 kbp bin pairs.

<<>>=
smaller.data <- normOffsets(smaller.data, type="loess", se.out=TRUE)
y.small <- asDGEList(smaller.data)
y.small <- estimateDisp(y.small, design)
fit.small <- glmQLFit(y.small, design, robust=TRUE)
result.small <- glmQLFTest(fit.small)
@

A combined $p$-value is computed for each cluster from the $p$-values of its bin pairs, using Simes' method \cite{simes1986} as implemented in the \Rfunction{combineTests} function from the \Biocpkg{csaw} package.
Each combined $p$-value represents the evidence against the global null hypothesis, i.e., that none of the constituent bin pairs are significant in the cluster.
Cluster-level FDR control is maintained by applying the BH method on the combined $p$-values for all clusters.
Significant clusters are then defined at a certain FDR threshold, e.g., 5\%.

<<>>=
library(csaw)
tabcluster <- combineTests(clustered.small$indices[[1]], result.small$table)
head(tabcluster)
sum(tabcluster$FDR <= 0.05)
@

The \Rfunction{combineTests} function also reports details about each cluster, including the total number of bin pairs and the number changing in each direction.
These results are saved to file below, along with the coordinates of the minimum bounding box for each cluster.

<<>>=
inter.frame <- as.data.frame(clustered.small$interactions)[,useful.cols]
results.i <- data.frame(inter.frame, tabcluster)
o.i <- order(results.i$PValue)
write.table(results.i[o.i,], file="independent.tsv", sep="\t", 
            quote=FALSE, row.names=FALSE)
@

The statistics can also stored in the metadata of the \Rclass{GInteractions} object returned by \Rfunction{clusterPairs}. 
This ensures that the genomic coordinates and statistics are synchronised.

<<>>=
mcols(clustered.small$interactions) <- cbind(
    mcols(clustered.small$interactions), tabcluster)
@

\subsection{Clustering based on significant bin pairs}
\label{sec:sigclust}
The interaction space is often dominated by high-abundance structural features like domains and compartments. 
In such high-density areas, too many bin pairs may be retained after independent filtering on abundance.
Clustering would then result in excessive chaining, yielding very large clusters for which the coordinates provide no information on DIs.

To avoid this, it may be necessary to cluster using only those bin pairs that are significantly different.
These bin pairs are more sparsely distributed throughout the interaction space, such that the coordinates of the resulting clusters can be easily interpreted.
This is done using the \Rfunction{diClusters} function, as shown below for the results with the 1 Mbp bin pairs.
Each cluster represents a DI containing only significant bin pairs.
This is more useful than reporting the coordinates for an entire domain if only a portion of the domain is changing.

<<>>=
clustered.sig <- diClusters(data, result$table, target=0.05, 
                            cluster.args=list(tol=1))
length(clustered.sig$interactions)
head(clustered.sig$interactions)
@

The \Rfunction{clusterPairs} function is used internally to cluster significant bin pairs. 
No limits are placed on the maximum dimensions as the sparsity of the selected bin pairs should avoid chaining effects.
Additional sophistication can be obtained by setting \Rcode{fc.col="logFC"} above.
This will cluster bin pairs separately if they are changing in opposite directions.

The \Rfunction{diClusters} function will attempt to control the cluster-level FDR below \Rcode{target}, i.e., 5\%.
Specifically, the cluster-level FDR is estimated from the FDR threshold used to define significant bin pairs (see the \Rfunction{clusterFDR} function in the \Biocpkg{csaw} package).
The bin pair-level threshold is then adjusted until the estimate of the cluster-level FDR lies close to \Rcode{target}.
The estimated cluster-level FDR is stored in the \Rcode{FDR} field of the output.
Formally speaking, this procedure is not entirely rigorous and will yield biased estimates for small numbers of clusters.
However, this may be a necessary price to pay for interpretable results.

<<>>=
clustered.sig$FDR
@

The identities of the bin pairs in each cluster are also returned in the \Rcode{indices} field of the output list.
These can be used to compute additional statistics for each cluster using the \Rfunction{combineTests} and \Rfunction{getBestTest} functions in \Biocpkg{csaw}.
The latter will identify the bin pair with the lowest $p$-value, which is useful for finding the strongest changes in large clusters.
Users are advised to ignore the $p$-value and FDR fields as these assume independent clustering.

<<>>=
tabcom <- combineTests(clustered.sig$indices[[1]], result$table)
head(tabcom)
tabbest <- getBestTest(clustered.sig$indices[[1]], result$table)
head(tabbest)
@

Test statistics are saved to file for later examination, along with the coordinates of each cluster's bounding box.
The combined $p$-value from \Rfunction{combineTests} is stored but is only used for sorting and should not be interpreted as a significance measure.

<<>>=
tabstats <- data.frame(tabcom[,1:4], logFC=tabbest$logFC, FDR=clustered.sig$FDR)
results.d <- as.data.frame(clustered.sig$interactions)[,useful.cols]
results.d <- cbind(results.d, tabstats)
o.d <- order(results.d$PValue)
write.table(results.d[o.d,], file="DIclusters.tsv", sep="\t", 
            quote=FALSE, row.names=FALSE)
@

Again, the statistics can stored in the metadata of \Rclass{GInteractions} object.

<<>>=
mcols(clustered.sig$interactions) <- cbind(
    mcols(clustered.sig$interactions), tabstats)
@

\section{Merging results from different bin widths}
\label{sec:mergebins}

\subsection{Clustering with different bin widths}
The optimal choice of bin size is not clear when there are both sharp and diffuse changes in the interaction space.
Smaller bins provide greater spatial resolution and can identify sharp DIs that would be lost within larger bin pairs.
Conversely, larger bin pairs have larger counts and greater power to detect diffuse DIs.
Comprehensive detection of DIs can be achieved by combining analyses from several bin sizes.
For example, \Rfunction{clusterPairs} accepts multiple \Rclass{InteractionSet} objects to cluster results of analyses with different sizes.

<<>>=
clustered.mult <- clusterPairs(larger=data, smaller=smaller.data, 
                               tol=1, upper=1e6)
head(clustered.mult$indices$larger)
head(clustered.mult$indices$smaller)
@

Alternatively, the \Rfunction{boxPairs} function can be used to identify all smaller bin pairs that are nested within each of the larger ``parent'' bin pairs.
This is a form of independent clustering where all nested bin pairs are defined as a single cluster.
Each set of nested bin pairs will only be reported once, reducing redundancy and limiting the potential for misinterpretation of the FDR.
Note that the larger bin size \textit{must} be an integer multiple of the smaller bin size(s).
This is necessary to simplify the interpretation of the nesting procedure.

<<>>=
boxed <- boxPairs(larger=data, smaller=smaller.data)
head(boxed$indices$larger)
head(boxed$indices$smaller)
@

\subsection{Computing consolidated cluster-level statistics}
Regardless of whether \Rfunction{clusterPairs} or \Rfunction{boxPairs} is used, the statistics for each cluster can be computed with the \Rfunction{consolidatePairs} function.
This uses a weighted version of Simes' method to ensure each analysis contributes equally to the significance of each cluster \cite{benjamini1997multiple}.
For each bin pair, the weight of its $p$-value is inversely proportional to the number of bin pairs of the same size in the same cluster. 
This ensures that the combined $p$-value calculation for each cluster is not dominated by smaller, more numerous bin pairs.
The BH method is then applied to the combined $p$-values to control the cluster-level FDR.

<<>>=
merged.results <- list(result$table, result.small$table)
cons <- consolidatePairs(boxed$indices, merged.results)
head(cons)
sum(cons$FDR <= 0.05)
@

Here, the number of detections is greater than that found with large bins alone.
This suggests that the different bin sizes complement each other by detecting features at different resolutions.
Statistics for each of the larger bin pairs can then be stored to file. 
Reordering is performed using the combined $p$-value to promote the strongest changes.

<<>>=
results.b <- data.frame(as.data.frame(boxed$interactions)[,useful.cols], cons)
o.b <- order(results.b$PValue)
write.table(results.b[o.b,], file="boxed.tsv", sep="\t", 
            quote=FALSE, row.names=FALSE)
@

\subsection{Clustering significant bin pairs of different widths}
The \Rfunction{diClusters} function will also accept multiple \Rclass{InteractionSet} objects in the form of a list.
It will simply pass the arguments to \Rfunction{clusterPairs} and control the cluster-level FDR as described in Section~\ref{sec:sigclust}.
It will also use a frequency-weighted version of the FDR to ensure an equal contribution from each bin size when defining significant bin pairs.

<<>>=
merged.data <- list(data, smaller.data)
clustered.mult <- diClusters(merged.data, merged.results, target=0.05, 
                             cluster.args=list(tol=1))
length(clustered.mult$interactions)
@

The output is also directly compatible with \Rfunction{consolidatePairs}.
Again, users should ignore the $p$-values when clustering based on significant bin pairs. 

<<>>=
cons.sig <- consolidatePairs(clustered.mult$indices, merged.results)
@


\section{Reporting nested bin pairs}
\label{sec:nesting}
It is often convenient to identify the top-ranked bin pair nested within each of the larger features, i.e., parent bin pairs or clusters.
Here, the top-ranked bin pair is identified by \Rfunction{getBestTest} as the one with the smallest individual $p$-value.
This means that any high-resolution changes nested within a large feature can be easily identified.
However, keep in mind that the FDR is computed with respect to the features, not the nested bin pairs.

<<>>=
inside <- getBestTest(boxed$indices$smaller, result.small$table)
best.interactions <- interactions(smaller.data)[inside$best,]
inter.frame <- as.data.frame(best.interactions)[,useful.cols[-c(1,4)]]
nested <- data.frame(inter.frame, inside[,c("logFC", "F")])
head(nested)
expanded <- rep(NA, nrow(results.b)) # For parents with no nested elements.
expanded[as.integer(rownames(inside))] <- seq_len(nrow(inside))
results.b <- data.frame(results.b, best=nested[expanded,])
write.table(results.b[o.b,], file="boxed_best.tsv", sep="\t", 
            quote=FALSE, row.names=FALSE)
@

The above code only reports the top-ranked nested bin pair within each large feature.
This may not be sufficient when many internal changes are occurring.
An alternative approach is to store the entirety of the \Robject{smaller.data} in a \R{} save file, along with \Rcode{cons} and \Rcode{data}.
Any interesting nested changes can then be interactively identified for a given feature.

\section{Visualization with plaid plots}
\label{chap:plaid}

\subsection{Using conventional plaid plots}
Plaid plots are widely used to visualize the distribution of read pairs in the interaction space \cite{lieberman2009comprehensive}. 
In these plots, each axis is a chromosome segment. 
Each ``pixel'' represents an interaction between the corresponding intervals on each axis. 
The color of the pixel is proportional to the number of read pairs mapped between the interacting loci.
We demonstrate below using large bin pairs detected in the consolidated analysis with small nested bin pairs.

<<plaidplot1,fig.wide=TRUE,fig.asp=0.5>>=
# Setting up the interaction space to be plotted.
chosen <- o.b[1]
chosen.a1 <- anchors(boxed$interactions[chosen], type="first")
chosen.a2 <- anchors(boxed$interactions[chosen], type="second")
expanded1 <- resize(chosen.a1, fix="center", width=bin.size*5)
expanded2 <- resize(chosen.a2, fix="center", width=bin.size*5)
cap.wt <- 100
cap.ko <- cap.wt*data$totals[3]/data$totals[1]

# Plotting the WT library.
par(mfrow=c(1,2))
plotPlaid(input[1], first=expanded1, second=expanded2, max.count=cap.wt, 
          width=5e4, param=mm.param, main="Flox")
rect(start(chosen.a1), start(chosen.a2), end(chosen.a1), end(chosen.a2))

# Plotting the KO library.
plotPlaid(input[3], first=expanded1, second=expanded2, max.count=cap.ko, 
          width=5e4, param=mm.param, main="KO")
rect(start(chosen.a1), start(chosen.a2), end(chosen.a1), end(chosen.a2))
@

<<results='hide',echo=FALSE>>=
if (as.character(seqnames(chosen.a1))!="chr11" || start(chosen.a1)!=30001298 || end(chosen.a1)!=31000010 ||
    as.character(seqnames(chosen.a2))!="chr11" || start(chosen.a2)!=29000942 || end(chosen.a2)!=30001301) {
	warning("first plaid plot hits the wrong spot")
}
@

Expansion of the plot boundaries ensures that the context of the interaction can be determined by examining the features in the surrounding space. 
It is also possible to tune the size of the pixels through a parameter that is, rather unsurprisingly, named \Rcode{width}.
In this case, the side of each pixel represents a 50 kbp bin, rounded to the nearest restriction site.
The actual bin pair occurs at the center of the plot and is marked by a rectangle.

The \Rcode{max.count} value controls the relative scale of the colors. 
Any pixel with a larger count will be set at the maximum color intensity.
This ensures that a few high-count regions do not dominate the plot.
A smaller \Rcode{max.count} is necessary for smaller libraries so that the intensity of the colors is comparable.
The actual color can be set by specifying \Rcode{col}.

In the example above, the differential interaction is driven mainly by the smaller bin pairs.
Changes in intensities are particularly prevalent at the top left and bottom right corners of the rectangle.
By comparison, the fold change for the entire bin pair is a little less than 30\%.
This highlights the usefulness of including analyses with smaller bin sizes.

Another example is shown below for the top DI detected in the analysis using only large bins.
Because the counts are ``averaged'' across the area of the interaction space, the change must be consistent throughout that area (and thus, more obvious) for detection to be successful.
Of course, any sharp changes within each of these large bin pairs will be overlooked as the smaller bin pairs are not used.

<<plaidplot2,fig.wide=TRUE,fig.asp=0.5>>=
# Setting up the interaction space to be plotted.
chosen <- o.r[1]
chosen.a1 <- anchors(data[chosen], type="first")
chosen.a2 <- anchors(data[chosen], type="second")
expanded1 <- resize(chosen.a1, fix="center", width=bin.size*5)
expanded2 <- resize(chosen.a2, fix="center", width=bin.size*5)
cap.wt <- 30
cap.ko <- cap.wt*data$totals[3]/data$totals[1]

# Plotting the WT sample.
par(mfrow=c(1,2))
plotPlaid(input[1], first=expanded1, second=expanded2, max.count=cap.wt, 
          width=5e4, param=mm.param, main="Flox")
rect(start(chosen.a1), start(chosen.a2), end(chosen.a1), end(chosen.a2))

# Plotting the KO sample.
plotPlaid(input[3], first=expanded1, second=expanded2, max.count=cap.ko, 
          width=5e4, param=mm.param, main="KO")
rect(start(chosen.a1), start(chosen.a2), end(chosen.a1), end(chosen.a2))
@

<<results='hide',echo=FALSE>>=
if (as.character(seqnames(chosen.a1))!="chr3" || start(chosen.a1)!=151001261 || end(chosen.a1)!=152002835 ||
    as.character(seqnames(chosen.a2))!="chr3" || start(chosen.a2)!=145998892 || end(chosen.a2)!=147005847) {
	warning("second plaid plot hits the wrong spot")
}
@

\subsection{Using rotated plaid plots}
Alternatively, users may prefer to use \Rfunction{rotPlaid} to generate rotated plaid plots.
These are more space-efficient and are easier to stack onto other genomic tracks, e.g., for ChIP-seq data.
However, rotated plots are only effective for local interactions within a specified region.
Some more effort is also required in interpretation.
In the example below, each colored box represents an interaction between two bins. 
The coordinates of each interacting bin can be identified by extending lines from opposite sides of the box until they intersect the $x$-axis.

<<rotplot, fig.asp=1.2>>=
# Setting up the interaction space to be plotted.
chosen <- o.b[4]
example <- anchors(boxed$interactions[chosen], type="second")
end(example) <- end(anchors(boxed$interactions[chosen], type="first"))
best.mid.a1 <- (results.b$best.start1[chosen]+results.b$best.end1[chosen])/2
best.mid.a2 <- (results.b$best.start2[chosen]+results.b$best.end2[chosen])/2
best.mid <- (best.mid.a1 + best.mid.a2)/2
best.gap <- best.mid.a1 - best.mid.a2

# Plotting the WT sample.
par(mfrow=c(2,1))
rotPlaid(input[1], mm.param, region=example, width=2.5e4, 
         main="Flox", max.count=cap.wt)
points(best.mid, best.gap, cex=7) 

# Plotting the KO sample.
rotPlaid(input[3], mm.param, region=example, width=2.5e4,
         main="KO", max.count=cap.ko)
points(best.mid, best.gap, cex=7) 
@

<<echo=FALSE,results='hide'>>=
.chosen.a1 <- anchors(boxed$interactions[chosen], type="first")
.chosen.a2 <- anchors(boxed$interactions[chosen], type="second")
if (as.character(seqnames(.chosen.a1))!="chr2" || start(.chosen.a1)!=136998216 || end(.chosen.a1)!=138001392 ||
    as.character(seqnames(.chosen.a2))!="chr2" || start(.chosen.a2)!=136001779 || end(.chosen.a2)!=136998219) {
	warning("third plaid plot hits the wrong spot")
}
@

The circle marks the area of the interaction space that corresponds to the top-ranked nested bin pair within the \Rcode{chosen} larger bin pair.
An increase in the interaction intensity is clearly observed in the KO condition.
This sharp change would not be observed with larger bin pairs, where the final count would be dominated by other (non-differential) areas.

\subsection{Using differential plaid plots}
In some cases, it may be more informative to display the magnitude of the changes across the interaction space.
This is achieved using the \Rfunction{plotDI} function, which assigns colors to bin pairs according to the size and direction of the log-fold change.
Visualization of the changes is useful as it highlights the DIs, whereas conventional plaid plots are dominated by high-abundance features like TADs.
The latter features may be constant between libraries and, thus, not of any particular interest.
The log-fold changes also incorporate normalization information, which is difficult to represent on a count-based plaid plot.

<<diplot, fig.small=TRUE>>=
chosen <- o.r[5]
chosen.a1 <- anchors(data[chosen], type="first")
chosen.a2 <- anchors(data[chosen], type="second")
expanded1 <- resize(chosen.a1, fix="center", width=5e7)
expanded2 <- resize(chosen.a2, fix="center", width=5e7)
colfun <- plotDI(data, result$table$logFC, expanded1, expanded2, diag=FALSE)
@

<<results='hide',echo=FALSE>>=
if (as.character(seqnames(chosen.a1))!="chr16" || start(chosen.a1)!=70000341 || end(chosen.a1)!=70999295 ||
    as.character(seqnames(chosen.a2))!="chr16" || start(chosen.a2)!=65999261 || end(chosen.a2)!=67001066) {
	warning("first DI plot hits the wrong spot")
}
@

The example above uses red and blue for positive and negative log-fold changes, respectively.
White and near-white regions correspond to those with log-fold change close to zero.
Grey regions mark the parts of the space where no bin pairs are present in \Robject{data}, possibly due to filtering on abundance.
As a result, this approach tends to be less useful when high-abundance bin pairs are more sparsely distributed, e.g., for long-range interactions.
A rotated DI plot can be similarly constructed using the \Rfunction{rotDI} function.

Both \Rfunction{rotDI} and \Rfunction{plotDI} will invisibly return another function that maps log-fold changes to colors.
This can be used to generate a custom color bar, as shown below.
A similar function that maps counts to colors is also returned by \Rfunction{plotPlaid} and \Rfunction{rotPlaid}.

<<colorbar,fig.asp=0.3>>=
logfc <- -20:20/10
plot(0,0,type="n", axes=FALSE, xlab="", ylab="", xlim=range(logfc), ylim=c(0,1))
rect(logfc - 0.05, 0, logfc + 0.05, 1, col=colfun(logfc), border=NA)
axis(1, cex.axis=1.2)
mtext("logFC", side=1, line=3, cex=1.4)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Detecting differential domain boundaries}
\label{chap:domains}
\begin{combox}
This chapter, like the cheese, stands alone, so there's not much we require from other chapters.
We'll only need the \Robject{input} and \Robject{mm.param} objects from Chapter~\ref{chap:counting}.
\end{combox}

\section{Overview}
High intensity triangles are often observed on the diagonal of the intra-chromosomal interaction space.
These correspond to topologically associating domains (TADs) where loci within each domain interact more frequently than those between domains.
Such domains are proposed to control genomic behavior by limiting the scope for interactions and restraining the spread of chromatin marks \cite{nora2013segmental}.
At higher resolutions, small domains may also be formed due to looping between specific genomic elements \cite{rao2014kilobase}.

To identify these domains, Dixon \emph{et al.} defined a ``directionality index'' for each genomic locus \cite{dixon2012topological}.
This was computed by counting the number of read pairs mapped between the target locus and an ``upstream'' interval with higher genomic coordinates;
    counting the number of read pairs mapped between the target locus and a ``downstream'' interval with lower genomic coordinates;
    and taking the normalized difference of the counts.
A region at the start of a domain will interact preferentially with upstream regions in the same domain, compared to downstream regions in a different domain.
Conversely, the region at the end of each domain will interact preferentially with the downstream regions compared to upstream regions.
This results in a characteristic pattern of positive directionality indices at the start of the domain, followed by negative values at the end.
These patterns can be extracted with hidden Markov models (HMMs) to explicitly identify the domain boundaries.

An alternative analysis is to identify loci where the size or direction of the directionality statistic changes between conditions.
This complements the standard DI analysis by focusing on regions where the domain boundary changes in strength or orientation.
Thus, changes in domain definition between conditions can be quantified in a statistically rigorous manner.

\section{Loading up- and downstream counts}
Up- and downstream counts for each 100 kbp genomic bin are loaded using the \Rfunction{domainDirections} function.
This returns a \Rclass{RangedSummarizedExperiment} object, in which each row corresponds to a bin and each column corresponds to a library.
The two matrices \Rcode{"up"} and \Rcode{"down"} contain the counts to the up- and downstream regions, respectively.

<<>>=
finder <- domainDirections(input, mm.param, width=1e5, span=10)
finder
@

The two matrices are combined into a single matrix, and the total counts for each library are also loaded.
These counts are used to construct a \Rclass{DGEList} for analysis with \Biocpkg{edgeR}. 

<<>>=
all.counts <- cbind(assay(finder, "up"), assay(finder, "down"))
totals <- totalCounts(input, mm.param)
ydom <- DGEList(all.counts, lib.size=rep(totals, 2))
@

As an aside, the same counts can be used to compute the directionality index \cite{dixon2012topological}. 
A Gaussian HMM can then be fitted to explicitly define domain boundaries, using packages like \CRANpkg{depmixS4}.
However, this analysis is largely outside the scope of \Biocpkg{diffHic} and will not be discussed here.

\section{Constructing the design matrix}
A design matrix is constructed with condition-specific coefficients for the log-fold change between up- and downstream counts.
It also contains sample-specific blocking factors to avoid incorporating unnecessary variability in the dispersion estimate due to sample-specific effects.
Recall that each library contributes two sets of counts to the final matrix, so the vectors containing condition and sample information must be doubled appropriately.

<<>>=
Condition <- factor(c("flox", "flox", "ko", "ko"))
Sample <- factor(seq_along(input))
Condition <- rep(Condition, 2)
Sample <- rep(Sample, 2)
Direction <- rep(c("Up", "Down"), each=length(input))
design <- model.matrix(~0 + Sample + Direction:Condition)
@

Some further manipulation is performed to clean up the design matrix prior to downstream analysis.
Redundant coefficients that cause linear dependencies are removed, and the remaining columns are given simpler names.
For convenience, the condition-specific up/down log-fold changes will be referred to as directionality statistics in the text below.

<<>>=
design <- design[,!grepl("DirectionDown", colnames(design))]
colnames(design) <- sub("DirectionUp:", "", colnames(design)) 
design
@

\section{Pre-processing of the count data}
Filtering is performed to remove low-abundance bins that do not contain enough counts to reject the null hypothesis.
In general, this should have little effect as most of the counts should be very large.
This is because \Rfunction{domainDirections} effectively counts read pairs for highly local interactions (as the up- or downstream intervals are adjacent to each bin).

<<>>=
ab <- aveLogCPM(ydom)
keep <- ab > 0
ydom <- ydom[keep,]
summary(keep)
@

The same filter is applied to the set of genomic bins to match with the entries of \Robject{ydom}. 

<<>>=
cur.regions <- rowRanges(finder)[keep,]
@

Normalization is not performed here as it is not required. 
The model contains sample-specific blocking factors, which largely negates the need for normalization between samples.
Additionally, the differential test for each bin will compare directionality values between conditions.
This means that any biases between the up- and downstream counts within each library or condition (e.g., due to mappability of different regions) should cancel out.

\section{Testing for significant differences with \Biocpkg{edgeR}}
The analysis proceeds directly to dispersion estimation and GLM fitting.
The mean-dependent trend in the NB dispersions is modelled using \Rfunction{estimateDisp} as previously described.

<<bcvdomplot, fig.small=TRUE>>=
ydom <- estimateDisp(ydom, design)
plotBCV(ydom)
@

A GLM is fitted to the counts for each bin using the specified model and the trended NB dispersion, and QL dispersions are estimated for all bins using robust EB shrinkage.

<<qldomplot, fig.small=TRUE>>=
fitdom <- glmQLFit(ydom, design, robust=TRUE)
plotQLDisp(fitdom)
@

The aim of the differential analysis is to identify bins where the directionality statistics change between conditions.
The \Rfunction{makeContrasts} function is used to construct an appropriate contrast vector, and the QL F-test is applied to detect significant differences.

<<>>=
con <- makeContrasts(Conditionko - Conditionflox, levels=design)
resdom <- glmQLFTest(fitdom, contrast=con)
topTags(resdom)
@

The interpretation of these changes requires some care.
If the directionality statistic was positive in the WT cells, a positive \Rcode{logFC} would represent a strengthening of the domain boundary.
However, if it was negative, a positive \Rcode{logFC} would represent a weakening of the domain boundary or a reversal of the domain orientation.
To assist interpretation, users are advised to report the condition-specific directionality statistics in the output.

<<>>=
output <- data.frame(as.data.frame(cur.regions)[,1:3],
                     KO=fitdom$coefficients[,"Conditionko"]/log(2),
                     Flox=fitdom$coefficients[,"Conditionflox"]/log(2),
                     resdom$table)
output$FDR <- p.adjust(resdom$table$PValue, method="BH")
o <- order(output$PValue)
output <- output[o,]
head(output[,1:5], 10) # not showing test statistics for brevity
@

Users can also examine whether the absolute values of the directionalities increase or decrease between conditions.
This provides an indication of whether the domains become more or less well-defined.
Here, domain definition seems to weaken in the KO condition.

<<>>=
is.sig <- output$FDR <= 0.05
summary(is.sig)
change.type <- abs(output$KO) > abs(output$Flox)
summary(change.type[is.sig])
@

Finally, the region with the top change in directionality is visualized.
Each rotated plaid plot shows the surrounding interaction space, with the region itself highlighted in red.

<<domplot, fig.asp=1.2>>=
tophit <- cur.regions[o[1]]
expanded <- resize(tophit, fix="center", width=width(tophit)*10)

par(mfrow=c(2,1))
rotPlaid(input[1], mm.param, region=expanded, width=2.5e4, main="Flox")
segments(start(tophit), 0, end(tophit), 0, col="red", lwd=10)
rotPlaid(input[3], mm.param, region=expanded, width=2.5e4, main="KO")
segments(start(tophit), 0, end(tophit), 0, col="red", lwd=10)
@

<<results='hide',echo=FALSE>>=
if (as.character(seqnames(tophit))!="chr2" || start(tophit)!=48799546 || end(tophit)!=48901284) {
	warning("first domain change hits the wrong spot")
}
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Epilogue}

\begin{combox}
Congratulations on getting to the end.
As a reward for your efforts, here is a poem:
\begin{quote}
I once had a friend named Bj\"ork, \\
With him I would always talk, \\
But he was a pig, \\
So when he got big, \\
We killed him and ate his pork.
\end{quote}
\end{combox}

\section{Data sources}
All datasets are publicly available from the NCBI Gene Expression Omnibus (GEO).
The K562/GM06990 dataset \cite{lieberman2009comprehensive} was obtained using the GEO accession number GSE18199.
The neural stem cell data set \cite{sofueva2013cohesin} was obtained using the accession GSE49017.
Finally, the RWPE1 dataset  \cite{rickman2012oncogene} was obtained using the accession GSE37752.
All libraries were processed as described in Chapter~\ref{chap:prep}.
For some datasets, multiple technical replicates are available for each library.
These were merged together prior to read pair counting.

All libraries were downloaded in the Sequence Read Archive format (SRA).
The SRA files were aligned to the reference genome using the \file{sra2bam.sh} script at \url{https://github.com/LTLA/diffHicUsersGuide}.
Some software packages are also required to run this script, such as various tools from the \software{Picard} suite -- read the source code in \file{sra2bam.sh} for more details.
The resulting BAM files were then converted into index files using the \file{bam2hdf.R} script.


\section{Session information}
<<>>=
sessionInfo()
@

\bibliography{refhic}

\end{document}