Introduction

The purpose of this package is to detect meso-scale changes in chromatin conformation from 3C data. Typical 3C data might look like the following:

The segments shown are restriction enzyme digestion fragments. The counts show the numbers of experimentally captured interactions captured interactions between each of these fragments and the Paupar viewpoint. This viewpoint lies in the region 105661048 -105661864 and we see higher numbers of captures for digestion fragments near the viewpoint. Packages like r3Cseq attempt to determine p-values for the numbers of captures for each fragment. In order to do this, it estimates a ‘background’ level of interaction as a function of distance from the viewpoint. Comparing the count for each fragment to this background level allows it to estimate significance for each fragment.

This method makes only limited use of positional information. While it considers distance from the viewpoint in estimating background levels it ignores the proximity of individual fragments. But suppose we have ten consecutive fragments, each with counts only slightly above background. None may rise to statistical significance on its own, but their co-location may make the group highly significant. We will exploit this observation to look for statistically significant changes in counts between treatments or cell types.

Consider 3C data for the same viewpoint for two replicates each for two different cell types.

We wish to discover regions where there is a statisitcally significant change in the interaction levels between the two cell types.

Difference in mean normalized counts:

The first task is to normalize the values for each of the cell types. This will allow us to find the mean count for each of the digestion fragments in each of the cell types and consequently find the difference in the mean normalized counts for the two cell types. Normalizing depends on estimating a size factor for each replicate. Here we are relying on DESeq2::estimateSizeFactorsForMatrix(). For further details see the DESeq2 vignette or Love, M.I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 Genome Biology 15(12):550 (2014)

We have included a miniature version of the summarized experiment miniSE in order to demonstrate the process of extracting the delta mean normalized value miniDeltaSE.

{r}
library(deltaCaptureC)
miniDeltaSE = getDeltaSE(miniSE)

The actual deltaSE comes from a SummarizedExperiment representing the complete data set for these four replicates and is subsequently trimmed to a region on chromosome2 surrounding our region of interest. After binning to a bin size of 1kb it looks like this:

Our interface carries out the binning inside the call to the function getSignificantRegions(). The delta data (here deltaSE) is then mined for significant regions.

The Algorithm

The algorithm proceeds as follows:

1. The data is binned to smallBinSize (here 1kb) and trimmed to the region of interest. The region of interest here is 500kb up and downstream of the mid-point of the viewpoint.
2. We now rebin this data to bigBinSize (here 10kb). This value should be an integer multiple of smallBinSize.
3. The interaction count data are highest near the viewpoint. We therefore distinguish a viewpoint region.
   • In the viewpoint region, we compute the lopsidedness of the data, i.e., the absolute value of the difference between the sum of the bins to the left of the viewpoint and the sum of the bins to the right of the viewpoint.
   • Outside the viewpoint region we define runs to be maximal consecutive sequences of bins where the delta value does not change sign. We compute the totals of the bins in each run. We call these run totals.
4. We now perform random permutations of the small-bin delta data.
   • Within the viewpoint region, these permutations respect the distance from the viewpoint. Thus, for each distance k from the viewpoint, the two bins at distance k either are or are not swapped.
   • We perform arbitrary permutations on the non-viewpoint bins.
5. For each permutation, we rebin the resulting scrambled data to the large bin size and then compute the lopsidedness and run totals for that permutation.
6. For each permutation we record the lopsidedness and maximum and minimum of the run totals. Across all permutations this produces
   • a distribution of lopsidedness values
   • a distribution of maximum run totals
   • a distribution of minimum run totals
7. Significance in the original (unpermuted) data is then assessed by comparing the lopsidedness and run totals in the original data to these distributions.

Invoking the algorithm and plotting the results

significantRegions = getSignificantRegions(deltaSE,
                                           regionOfInterest,
                                           viewpointRegion,
                                           smallBinSize,
                                           bigBinSize,
                                           numPermutations,
                                           pValue)

which are then ploted using

library(deltaCaptureC)
significantRegionsPlot = plotSignificantRegions(significantRegions,
                                                significanceType,
                                                plotTitle)

This gives the following result:

Under the hood

library(SummarizedExperiment)
library(deltaCaptureC)
se = miniSE
counts = assays(se)[['counts']]
sizeFactors = DESeq2::estimateSizeFactorsForMatrix(counts)
colData(se)$sizeFactors = sizeFactors
assays(se)[['normalizedCounts']] = counts
for(i in seq_len(ncol(assays(se)[['normalizedCounts']])))
{
    assays(se)[['normalizedCounts']][,i] =
            assays(se)[['normalizedCounts']][,i] /
        colData(se)$sizeFactors[i]
}

library(SummarizedExperiment)
library(deltaCaptureC)
meanNormalizedCountsSE = getMeanNormalizedCountsSE(miniSE)
meanCounts = assay(meanNormalizedCountsSE)
delta = matrix(meanCounts[,1] - meanCounts[,2],ncol=1)
colData = data.frame(delta=sprintf('%s - %s',
                                    as.character(colData(meanNormalizedCountsSE)$treatment[1]),
                                    as.character(colData(meanNormalizedCountsSE)$treatment[2])),
                                    stringsAsFactors=FALSE)
deltaSE = SummarizedExperiment(assay=list(delta=delta),
                                          colData=colData)
rowRanges(deltaSE) = rowRanges(meanNormalizedCountsSE)

library(deltaCaptureC)
print(length(smallSetOfSmallBins))
#> [1] 201
print(length(smallerDeltaSE))
#> [1] 401
tictoc::tic('binning into small bins')
binnedSummarizedExperiment = binSummarizedExperiment(smallSetOfSmallBins,smallerDeltaSE)
tictoc::toc()
#> binning into small bins: 73.06 sec elapsed

library(deltaCaptureC)
tictoc::tic('rebinning to larger bin size')
rebinnedSummarizedExperiment =
      rebinToMultiple(binnedSummarizedExperiment,10)
tictoc::toc()
#> rebinning to larger bin size: 0.14 sec elapsed