---
title: "Consensus regions across replicated experiments"
author:
    -   name: "Sebastian Gregoricchio, Ph.D."
        affiliation:
            -   >-
                Netherlands Cancer Institute, Oncogenomics division,
                Amsterdam - Netherlands
        email: 's.gregoricchio@nki.nl'
date: "Last update: `r format(Sys.time(), '%d<sup>th</sup> %B, %Y')`"
output:
    BiocStyle::html_document:
        toc: true
        toc_float: true
        toc_depth: 3
        number_sections: true
        fig_caption: true
vignette: >
    %\VignetteIndexEntry{Consensus regions across replicated experiments}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

<style>
body {
text-align: justify}
</style>

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = ">",
    warning = FALSE, message = FALSE, fig.align = "center",
    dev = "png", dpi = 96, fig.width = 7, fig.height = 4.5)
```



----------------------------------------------------------

# **Introduction** {#intro}
Replicates disagree. A region can be convincing in two samples and just
below threshold in the third, and the usual response, taking the
intersection of the peak sets, throws it away. Taking the union instead
keeps it along with every piece of noise the worst replicate produced.

`consensusRegions` does neither. Each peak is examined together with the
peaks that overlap it in the other replicates, and their p-values are
combined into one. Repeated weak evidence then adds up to something that
survives, while a strong peak nobody else saw does not. The idea comes
from [MSPC](https://genometric.github.io/MSPC/docs/)
([Jalili *et al.*, 2015](https://doi.org/10.1093/bioinformatics/btv293));
what is added here is the ability to weight replicates, to calibrate the
decision threshold rather than fix it, and to work with peak files that
have lost their statistics.

```{r load}
library(consensusRegions)
```


<br>



# **Reading peaks** {#read_peaks}
The example data is three replicates over two chromosomes, with a
deliberately shallow third sample.

```{r read}
peakFiles <- system.file("extdata",
    c("rep1.narrowPeak",
        "rep2.narrowPeak",
        "rep3.narrowPeak"),
    package = "consensusRegions")

peaks <- readPeakSets(peakFiles, sampleNames = c("rep1", "rep2", "rep3"))
peaks
```

`readPeakSets()` inspects the files and reports which statistic it can
use. narrowPeak keeps -log10 p-values in its eighth column, so those go
straight in. A BED with only a score gets rank-transformed. A bare BED3
has nothing at all, and the analysis will fall back to counting
replicates instead of combining evidence.

<br>

## Chromosome naming
Peak files rarely agree on whether a chromosome is called `chr1` or `1`.
Two replicates that name the same chromosome differently share no overlap
at all, so every peak is discarded for lack of support and the result
looks like a biological finding rather than a clerical one.

`readPeakSets()` therefore harmonises the names on the way in, defaulting
to UCSC style, in the same way as `RegionSetDE::loadRegions()`:

```{r seqlevels}
withChr <- GenomicRanges::GRanges(
    "chr1", IRanges::IRanges(start = c(1000, 20000), width = 500))
withoutChr <- GenomicRanges::GRanges(
    "1", IRanges::IRanges(start = c(1100, 20100), width = 500))

mixed <- readPeakSets(list(a = withChr, b = withoutChr),
    seqlevelsStyle = "UCSC")
GenomeInfoDb::seqlevels(mixed[["b"]])
```

`"Ensembl"` and `"NCBI"` go the other way, stripping the prefix. Setting
`seqlevelsStyle = NULL` leaves the names untouched and stops if the
replicates disagree, which is the safe choice when you would rather see
the mismatch than have it quietly repaired.

<br>

## Call peaks permissively
This is a critical issue that can substantially affect the validity of
the entire analysis. If you call peaks at `q < 0.05` and run the
consensus on what survives, there is nothing left to rescue and you have
written an intersection with extra steps. Call at something like
`p < 0.001` and let `buildConsensus()` do the thresholding. The two-tier
design assumes exactly this.

<br>




# **Quick start** {#quick_start}
```{r basic}
result <- buildConsensus(peaks, verbose = FALSE)
result
```

The per-replicate summary is where to look first.

```{r stats}
consensusStats(result)
```

`nRescued` counts peaks that were only weak on their own and were kept
because the others agreed. That number is the reason to use this method,
and also the thing to be suspicious of. If it dwarfs everything else,
`weakThreshold` is too permissive and the combination is confirming
noise.

```{r plotRescue}
plotRescue(result)
```

The consensus regions carry the number of contributing replicates and a
significance recomputed from their member peaks.

```{r consensus}
head(consensusRanges(result), 3)
```

<br>



# **Weighting the replicates** {#weighting_reps}
Fisher's method treats every replicate as equally informative. Yours
probably are not: one is shallower, or has a worse antibody, or was
prepared on a bad day. The choices are usually to drop it, which wastes
data, or to keep it and let it drag the combined statistic around.
Weighting is the third option.

When the BAM files are available, FRiP and library size can be measured
directly:

```{r weightsBam, eval = FALSE}
weights <- computeReplicateWeights(
    peaks,
    method = "frip",
    bamFiles = c("rep1.bam", "rep2.bam", "rep3.bam"))
```

More often the alignments are long gone and only the numbers survive in
an old QC report. Pass them in:

```{r weightsFrip}
weights <- computeReplicateWeights(peaks,
    method = "frip",
    frip = c(0.21, 0.19, 0.06),
    verbose = FALSE)
weights
```

And when there is nothing but the peaks themselves, each replicate can be
scored by how well it agrees with the others:

```{r weightsIntrinsic}
computeReplicateWeights(peaks,
    method = "intrinsic",
    verbose = FALSE)
```

That last measure is circular, since agreement is also what the
consensus is testing. Use it when there is no alternative, not in
preference to FRiP.

Weights only reach the combination through Stouffer's Z or Lancaster's
method. Fisher and the rank product ignore them.

```{r weighted}
weighted <- buildConsensus(peaks,
    weights = weights,
    combinationMethod = "stouffer",
    verbose = FALSE)
consensusStats(weighted)
```

```{r jaccard}
plotJaccard(weighted)
```

<br>

## Which weighting method to use
The five settings differ in what they score a replicate on, and in what
they cost you to obtain.

| `method`        | Scored on | Needs                    | Weight is      |
|-----------------|-----------|--------------------------|----------------|
| `"equal"`       | nothing   | nothing                  | `1` throughout |
| `"frip"`        | signal    | `bamFiles` or `frip`     | the fraction   |
| `"librarySize"` | depth     | `bamFiles`/`librarySize` | sqrt of depth  |
| `"intrinsic"`   | agreement | the peak sets alone      | mean Jaccard   |
| `"custom"`      | your call | a numeric vector         | what you pass  |

`"frip"` takes the fraction straight from an old QC report when the
alignments are gone, and `"intrinsic"` scores a replicate on its mean
pairwise Jaccard index with the others.

Whichever you pick, the weights are rescaled to average one, so the
combined statistics stay on the scale they would have had unweighted and
a threshold carries across runs.

A few extra details:

*   `"equal"` is the default and the reproducible choice. Nothing about the
    data steers the result, which is what you want unless you have a reason
    to think one replicate deserves less say than another.
*   `"frip"` is the one to reach for when the alignments, or the numbers
    someone once computed from them, still exist. It measures the thing you
    care about: how much of the sequencing landed in peaks.
*   `"librarySize"` counts depth twice to some extent, because peak-caller
    p-values already respond to it. Prefer `"frip"` when both are at hand.
*   `"intrinsic"` is circular. It scores a replicate on the agreement that
    the consensus then goes on to test, so strong agreement earns a large
    weight which in turn makes the agreement look stronger. It is there for
    the case where all that survives is a folder of BED files, not as an
    alternative to FRiP.
*   There is no `"drop"`. Removing a bad replicate is something you do
    before `readPeakSets()` ever sees it, and weighting exists so that you
    usually do not have to.

<br>





# **Choosing the combination method** {#choosing_combination}

| Method        | Weights | Use it when                                     |
|---------------|---------|-------------------------------------------------|
| `fisher`      | no      | Reproducing MSPC, or replicates are comparable  |
| `stouffer`    | yes     | Default. Replicate quality varies               |
| `lancaster`   | yes     | Weights, with Fisher's pull from one strong peak |
| `rankProduct` | no      | Only scores, or p-values you do not trust       |

One caveat about Fisher that matters more than it usually gets credit
for: the combined statistic keeps accumulating as replicates are added.
With three samples this is fine. With twenty donors, a region detected
faintly in all of them starts to look overwhelmingly significant, which
is not what anyone means by reproducible. Stouffer and the rank product
degrade more gracefully at that scale.

The threshold does not carry across methods. Fisher, Stouffer and
Lancaster all sit on roughly the same scale, so a cut chosen for one is
a reasonable starting point for the others. The rank product does not.
It works on relative ranks, so the best a peak can possibly do is come
first in every replicate, and that bound depends on how many peaks each
replicate holds. With the few hundred peaks in this example the smallest
attainable combined p-value is around `1e-4`, and asking for `1e-8`
would return nothing at all. `buildConsensus()` checks for this and
refuses rather than handing back an empty result:

```{r rankProductScale, error = TRUE}
buildConsensus(peaks,
    combinationMethod = "rankProduct",
    verbose = FALSE)
```

Give it a threshold on its own scale, or let the calibration below pick
one:

```{r rankProductWorking}
ranked <- buildConsensus(peaks,
    combinationMethod = "rankProduct",
    combinedThreshold = 0.05,
    verbose = FALSE)
length(ranked)
```

Nothing survives, and the reason is the size of the example rather than
any weakness in the method. With only a few hundred peaks per replicate
the finest relative rank available is about 1/430, which is not a small
number, so no peak accumulates much evidence however well it ranks. What
little it does accumulate is then spent on the multiple testing
correction described in the next section. Relax that correction to the
one MSPC uses and the peaks reappear:

```{r rankProductRelaxed}
relaxed <- buildConsensus(peaks,
    combinationMethod = "rankProduct",
    combinedThreshold = 0.05,
    adjustmentFamily = "confirmed",
    verbose = FALSE)
length(relaxed)
```

On a real peak set of 10^5 calls the ranks are three orders of magnitude
finer and the rank product survives the stricter correction comfortably.
Judge it on your own data, not here.

<br>



# **How many replicates have to agree** {#how_many_replicates}
The requirement can be stated either way round, because both readings are
natural and mixing them up is easy.

`minSupport` counts the *other* replicates that must hold an overlapping
peak, which is what the confirmation step actually works with.
`minReplicates` counts all of them including the one the peak came from,
which is how the requirement is usually spoken about, and it matches
MSPC's `-c`. Give one or the other, never both.

```{r minReplicates}
c(byTotal = length(buildConsensus(peaks,
            minReplicates = 3,
            verbose = FALSE)),
    bySupport = length(buildConsensus(peaks,
            minSupport = 2,
            verbose = FALSE)))
```

`minReplicates` also takes a proportion or a percentage, rounded up, so
the requirement can be written the way it is usually described:

```{r minReplicatesPercent}
length(buildConsensus(peaks,
        minReplicates = "60%",
        verbose = FALSE))
```

| Requirement | MSPC     | `minReplicates`   | `minSupport` |
|-------------|----------|-------------------|--------------|
| 2 of 3      | `-c 2`   | `2`               | `1`          |
| 3 of 3      | `-c 3`   | `3` or `"100%"`   | `2`          |
| 5 of 8      | `-c 5`   | `5`               | `4`          |
| 75% of 8    | `-c 75%` | `"75%"` or `0.75` | `5`          |

If you are coming from MSPC, use `minReplicates` and pass the same value
you passed to `-c`. Setting `minSupport` to that value instead asks for
one replicate more than you meant, which is a quiet way to lose peaks.

Note that `replicateType = "technical"` overrides both and demands every
replicate, on the grounds that technical replicates have no reason to
disagree.

<br>

# **Multiple testing** {#multiple_testing}
Peaks that clear the combined threshold are corrected within each
replicate by the Benjamini-Hochberg procedure. The question is which
peaks make up the family being corrected, and it matters more than it
first appears.

MSPC corrects across the peaks that already cleared the threshold. That
family was chosen for being significant, so the usual Benjamini-Hochberg
guarantee does not carry over to it, and the resulting values can be
optimistic. `adjustmentFamily = "tested"`, the default here, instead
corrects across every peak that entered the analysis.

```{r adjustmentFamily}
tested <- buildConsensus(peaks, verbose = FALSE)

mspcStyle <- buildConsensus(peaks,
    adjustmentFamily = "confirmed",
    verbose = FALSE)

c(tested = length(tested), confirmed = length(mspcStyle))
```

On this data the two agree, and that is the normal case: when the
combined evidence is strong the correction is a rounding error against
it. The choice only bites where the evidence was marginal, which is
where you would want it to bite. Use `"confirmed"` when you are
deliberately reproducing MSPC.

Neither setting delivers exact false discovery rate control, and the
package does not claim it does. Two things stand in the way. The
supporting peak for each replicate is chosen as the most significant of
those overlapping, and a maximum over a selected set is not distributed
the way an unselected p-value would be. Overlapping peaks also share
supporters, so their combined p-values are not independent. If you need
a defensible error rate rather than a ranking, use the empirical route
in the next section.

<br>



# **Calibrating the threshold** {#calibrating_threshold}
The cut on the combined p-value is conventionally `1e-8`. That number is
a guess, and the right one depends on the number of replicates, how the
peaks were called, and how densely they sit.

`calibrateThreshold()` estimates it. Peak positions are shuffled within
each chromosome, keeping widths, counts and statistics, so any overlap
between the shuffled replicates arises by chance. Running the combination
on those gives a null, and the threshold is read off where the expected
number of null peaks falls to a chosen fraction of the observed count.

```{r calibrate}
calibration <- calibrateThreshold(peaks,
    nPermutations = 10,
    targetFDR = 0.05,
    seed = 42,
    verbose = FALSE)

calibration$threshold
```

```{r plotCalibration}
plotCalibration(calibration)
```

```{r calibrated}
calibrated <- buildConsensus(peaks,
    combinedThreshold = calibration$threshold,
    verbose = FALSE)
length(calibrated)
```

Ten permutations is enough for a demonstration; fifty is a reasonable
working number. The positions are drawn at random, and `seed` is what
brings the same threshold back. The draws happen inside `BiocParallel`,
which keeps its random numbers apart from those of the session, so a
`set.seed()` before the call does not reach them. The seed is kept with
the result, in `calibration$seed`, and the random numbers of the session
are left as they were.

The permutations are independent of one another and are where nearly all
the time goes, so they are handed to `BiocParallel`. On a genome-scale
peak set a single round takes tens of seconds, which makes fifty of them
most of an afternoon in series:

```{r parallel, eval = FALSE}
calibration <- calibrateThreshold(peaks, nPermutations = 50, BPPARAM = 8,
    seed = 42)
```

`BPPARAM` takes the number of cores directly, which is how the question
usually gets asked, and the same seed returns the same threshold on one
core or on eight. A `BiocParallelParam` object is accepted too, for finer
control of the workers. It is used as it is, so the seed then travels with
the backend, as its `RNGseed`, and `seed` is ignored:

```{r parallelSeed, eval = FALSE}
calibration <- calibrateThreshold(
    peaks,
    nPermutations = 50,
    BPPARAM = BiocParallel::MulticoreParam(workers = 8, RNGseed = 42))
```

Passing a blacklist through `excludeRegions` makes the
null more honest, because shuffled peaks otherwise land in artefact
regions where real peaks cluster too.


<br>



# **Input without statistics** {#input_no_stats}
Sometimes all that survives of an old analysis is a folder of BED3 files.
There is no p-value to combine, so the package says so and switches to a
weighted presence rule: a region is kept when the weights of the
replicates supporting it sum to at least `minSupport`. With equal weights
that is the familiar k-of-n rule.

```{r bed3}
minimalFile <- system.file("extdata", "rep1_minimal.bed",
    package = "consensusRegions")

minimal <- readPeakSets(rep(minimalFile, 3),
    sampleNames = c("a", "b", "c"))

minimalResult <- buildConsensus(minimal, verbose = FALSE)
analysisParameters(minimalResult)$presenceOnly
```

This is a fallback, not an equivalent. Without per-peak statistics
nothing can be rescued, and the result is reproducibility filtering
rather than combined evidence.



<br>


# **Merging, and relative issues** {#merging}
Peak A overlaps B, B overlaps C, and A never touches C. Transitive
merging puts all three in one region, and on permissive input this
occasionally produces regions tens of kilobases long that mean nothing.

Two defences. `maxConsensusWidth` catches oversized regions and rebuilds
only those. `mergeMethod = "iterative"` avoids the problem entirely by
seeding on the most significant peak, claiming what it overlaps, and
repeating, which is the approach the ATAC-seq peak atlases use.

```{r merge}
seeded <- buildConsensus(peaks, mergeMethod = "iterative",
    verbose = FALSE)

summary(GenomicRanges::width(consensusRanges(result)))
summary(GenomicRanges::width(consensusRanges(seeded)))
```

<br>



# **Fixed-width regions** {#fixed_width}
Consensus regions inherit their edges from whichever peaks were merged,
so their widths vary. That variation biases everything downstream: a
wider region collects more reads and more motif matches for reasons that
have nothing to do with biology. For transcription factors and ATAC data,
recentre on the summit and use a common width.

```{r recentre}
fixed <- recentrePeaks(result, width = 400, summitSource = "best")
head(fixed, 3)
```

Do not do this to broad histone domains. Collapsing a domain to a point
discards the thing being measured.


<br>



# **Handing over to a differential analysis**
If the consensus is a step on the way to differential binding rather than
the endpoint, convert it and count reads over the row ranges with csaw,
DiffBind or edgeR.

```{r se}
se <- asSummarizedExperiment(result)
se
```

The assay holds per-replicate significance, not counts. It is useful for
spotting a replicate that disagrees with the others; it is not a
substitute for counting reads.

And a word on where to spend your effort. If the consensus is an
intermediate, do not over-engineer it. Aggressive filtering upfront
mostly discards peaks that a proper count-based model would have handled
correctly anyway. The stringency of this package earns its keep when the
peak set is the deliverable.

<br>



# **Blacklisting**
Filter after the consensus, not before. Removing regions upfront distorts
the rank distributions the combination relies on, which is why
`excludeRegions` is applied at the end of `buildConsensus()`.

```{r blacklist, eval = FALSE}
result <- buildConsensus(peaks, excludeRegions = blacklistGRanges)
```

<br>



# **Checking the result**
Nothing in this package can tell you whether your consensus set is any
good. Validate it against something orthogonal: central motif enrichment
for a transcription factor, TSS enrichment and cCRE overlap for ATAC,
correlation with expression for activating marks. If the rescued peaks
are motif-poor and TSS-depleted, the settings are too permissive
regardless of what the statistics say.

<br>

----------------------------------------------------

# **Everything in one call** {#wrapper}
The steps above almost always run in the same order, so `runConsensus()`
chains them: read, weight, optionally calibrate, build, optionally
recentre, optionally write out.

```{r wrapper}
result <- runConsensus(peakFiles,
    sampleNames = c("rep1", "rep2", "rep3"),
    minReplicates = 2,
    verbose = FALSE)
length(result)
```

Anything it does not name itself is passed through to
`buildConsensus()`, so `combinationMethod`, `mergeMethod`,
`adjustmentFamily` and the rest all work here too:

```{r wrapperOptions, eval = FALSE}
result <- runConsensus(peakFiles,
    sampleNames = c("rep1", "rep2", "rep3"),
    weightMethod = "frip",
    bamFiles = c("rep1.bam", "rep2.bam", "rep3.bam"),
    calibrate = TRUE,
    nPermutations = 50,
    combinationMethod = "stouffer",
    minReplicates = "75%",
    recentre = TRUE, width = 400,
    excludeRegions = blacklist,
    outputFile = "consensus.bed",
    BPPARAM = 8)
```

Use the individual functions when a step needs looking at before the next
one runs. Calibration in particular deserves a glance at
`plotCalibration()` rather than blind trust, and the wrapper gives you no
chance to take it.

<br>

----------------------------------------------------

# **Session info**
```{r sessionInfo, echo=FALSE}
sessionInfo()
```