correctedContact {diffHic} | R Documentation

## Iterative correction of Hi-C counts

### Description

Perform iterative correction on counts for Hi-C interactions to correct for biases between fragments.

### Usage

    correctedContact(data, iterations=50, exclude.local=1, ignore.low=0.02,
        winsor.high=0.02, average=TRUE, dist.correct=FALSE, assay=1)

### Arguments

`data` | an InteractionSet object produced by `squareCounts`
`iterations` | an integer scalar specifying the number of correction iterations
`exclude.local` | an integer scalar, indicating the distance off the diagonal under which bin pairs are excluded
`ignore.low` | a numeric scalar, indicating the proportion of low-abundance bins to ignore
`winsor.high` | a numeric scalar indicating the proportion of high-abundance bin pairs to winsorize
`average` | a logical scalar specifying whether counts should be averaged across libraries
`dist.correct` | a logical scalar indicating whether to correct for distance effects
`assay` | a string or integer scalar specifying the matrix to use from `data`

### Details

This function implements the iterative correction procedure described by Imakaev et al. in their 2012 paper. Briefly, this aims to factorize the count for each bin pair into the biases for each of the two anchor bins and the true interaction probability. The bias represents the ease of sequencing/mapping/other for the genome sequence in each bin.

The `data` argument should be generated by taking the output of `squareCounts` after setting `filter=1`. Filtering should be avoided as counts in low-abundance bin pairs may be informative upon summation for each bin. For example, a large count sum for a bin may be formed from many bin pairs with low counts. Removal of those bin pairs would result in loss of information.

For `average=TRUE`, if multiple libraries are used to generate `data`, an average count will be computed for each bin pairs across all libraries using `mglmOneGroup`. The average count will then be used for correction. Otherwise, correction will be performed on the counts for each library separately.

The maximum step size in the output can be used as a measure of convergence. Ideally, the step size should approach 1 as iterations pass. This indicates that the correction procedure is converging to a single solution, as the maximum change to the computed biases is decreasing.

### Value

A list with several components.

`truth`:

  a numeric vector containing the true interaction probabilities for each bin pair

`bias`:

  a numeric vector of biases for all bins

`max`:

  a numeric vector containing the maximum fold-change change in biases at each iteration

`trend`:

  a numeric vector specifying the fitted value for the distance-dependent trend, if `dist.correct=TRUE`

If `average=FALSE`, each component is a numeric matrix instead. Each column of the matrix contains the specified information for each library in `data`.

### Additional parameter settings

Some robustness is provided by winsorizing out strong interactions with `winsor.high` to ensure that they do not overly influence the computed biases. This is useful for removing spikes around repeat regions or due to PCR duplication. Low-abundance bins can also be removed with `ignore.low` to avoid instability during correction, though this will result in `NA` values in the output.

Local bin pairs can be excluded as these are typically irrelevant to long-range interactions. They are also typically very high-abundance and may have excessive weight during correction, if not removed. This can be done by removing all bin pairs where the difference between the first and second anchor indices is less than `exclude.local`. Setting `exclude.local=NA` will only use inter-chromosomal bin pairs for correction.

If `dist.correct=TRUE`, abundances will be adjusted for distance-dependent effects. This is done by computing residuals from the fitted distance-abundance trend, using the `filterTrended` function. These residuals are then used for iterative correction, such that local interactions will not always have higher contact probabilities.

Ideally, the probability sums to unity across all bin pairs for a given bin (ignoring `NA` entries). This is complicated by winsorizing of high-abundance interactions and removal of local interactions. These interactions are not involved in correction, but are still reported in the output `truth`. As a result, the sum may not equal unity, i.e., values are not strictly interpretable as probabilities.

### Author(s)

Aaron Lun

### References

Imakaev M et al. (2012). Iterative correction of Hi-C data reveals hallmarks of chromosome organization. Nat. Methods 9, 999-1003.

### See Also

`squareCounts`, `mglmOneGroup`

### Examples

    # Dummying up some data.
    set.seed(3423746)
    npts <- 100
    npairs <- 5000
    nlibs <- 4
    anchor1 <- sample(npts, npairs, replace=TRUE)
    anchor2 <- sample(npts, npairs, replace=TRUE)
    data <- InteractionSet(
        list(counts=matrix(rpois(npairs*nlibs, runif(npairs, 10, 100)), nrow=npairs)),
        GInteractions(anchor1=anchor1, anchor2=anchor2,
            regions=GRanges("chrA", IRanges(1:npts, 1:npts)), mode="reverse"),
    	colData=DataFrame(totals=runif(nlibs, 1e6, 2e6)))

    # Correcting.
    stuff <- correctedContact(data)
    head(stuff$truth)
    head(stuff$bias)
    plot(stuff$max)

    # Different behavior with average=FALSE.
    stuff <- correctedContact(data, average=FALSE)
    head(stuff$truth)
    head(stuff$bias)
    head(stuff$max)

    # Creating an offset matrix, for use in glmFit.
    anchor1.bias <- stuff$bias[anchors(data, type="first", id=TRUE),]
    anchor2.bias <- stuff$bias[anchors(data, type="second", id=TRUE),]
    offsets <- log(anchor1.bias * anchor2.bias)

    # Adjusting for distance, and computing offsets with trend correction.
    stuff <- correctedContact(data, average=FALSE, dist.correct=TRUE)
    head(stuff$truth)
    head(stuff$trend)
    offsets <- log(stuff$bias[anchors(data, type="first", id=TRUE),]) +
        log(stuff$bias[anchors(data, type="second", id=TRUE),])

[Package diffHic version 1.16.0 Index]