The bnbc User's Guide

Introduction

The r Biocpkg("bnbc") package provides functionality to perform normalization and batch correction across samples on data obtained from Hi-C [@HiCassay] experiments.

In this package we implement tools for general subsetting of and data extraction from sets of Hi-C contact matrices, as well as smoothing of contact matrices, cross-sample normalization and cross-sample batch effect correction methods.

r Biocpkg("bnbc") expects as input

  1. a GRanges object representing the genome assayed, with individual ranges having the width that is equal to the bin size used to partition the genome.
  2. a list of (square) contact matrices.
  3. a DataFrame object containing sample-level covariates (i.e. gender, date of processing, etc).

Citing bnbc

If you use this package, please cite [@bnbc].

Terminology

It is well appreciated that Hi-C contact matrices exhibit an exponential decay in observed number of contacts as a function of the distance between the pair of interacting loci. In this work we operate, as has recently been done (i.e. [@HiCRep]), on the set of all loci interacting at a specific distance, one chromosome at a time. For a given distance $k$, the relevant set of loci are listed in each contact matrix as the entries comprising the $k$-th matrix diagonal (with the main diagonal being referred to as the first diagonal). We refer to these diagonals as matrix "bands".

Dependencies

This document has the following dependencies

library(bnbc)

The ContactGroup class from bnbc

r Biocpkg("bnbc") uses the ContactGroup class to represent the set of contact matrices for a given set of genomic locis interactions. The class has 3 slots:

We expect each ContactGroup to represent data from 1 chromosome. We are thinking about a whole-genome container. An example dataset for chromosome 22 is supplied with the package.

data(cgEx)
cgEx

Creating a ContactGroup object requires specifying the 3 slots above:

cgEx <- ContactGroup(rowData=rowData(cgEx),
                     contacts=contacts(cgEx),
                     colData=colData(cgEx))

Note that in this example, we used the accessor methods for each of these slots; there are also corresponding 'setter' methods, such as rowData(cgEx)<-.

Printing a ContactGroup object gives the number of bins represented by the rowData slot, the width of the bin in terms of genomic distances (i.e. 100kb) and the number of samples:

cgEx

Alternatives

The r Biocpkg("InteractionSet") package contains a class called InteractionSet which is essentially an extension of the ContactGroup class. The internal storage format is different and InteractionSet is not restricted to square contact matrices like the ContactGroup class. We are interested in porting the bnbc() function to using InteractionSet, but bnbc() extensively uses band matrices and we have optimized r CRANpkg("Rcpp")-based routines for getting and setting bands of normal matrices, which means ContactGroup is a pretty fast for our purposes.

Getting your data into bnbc

To get data into r Biocpkg("bnbc") you need a list of contact matrices, one per sample. We assume the contact matrices are square, with no missing values. We do not require that data have been transformed or pre-processed by various bias correction software and provide methods for some simple pre-processing techniques.

There is currently no standard Hi-C data format. Instead, different groups produces custom formats, often in forms of text files. Because contact matrices are square, it is common to only distribute the upper or lower triangular matrix. In that case, you can use the following trick to make the matrix square:

mat <- matrix(1:9, nrow = 3, ncol = 3)
mat[lower.tri(mat)] <- 0
mat
## Now we fill in the lower triangular matrix with the upper triangular
mat[lower.tri(mat)] <- mat[upper.tri(mat)]
mat

Below, we demonstrate the steps needed to convert a set of hypothetical contact matrices into a ContactGroup object. The object upper.mats.list is supposed to be a list of contact matrices, each represented as an upper triangular matrix. We also suppose LociData to be a GenomicRanges object containing the loci of the contact matrices, and SampleData to be a DataFrame of per-sample information (i.e. cell type, sample name etc). We first convert all contact matrices to be symmetric matrices, then use the constructor method ContactGroup() to create the object.

## Example not run
## Convert upper triangles to symmetry matrix
MatsList <- lapply(upper.mats.list, function(M) {
    M[lower.tri(M)] <- M[upper.tri(M)]
})
## Use ContactGroup constructor method
cg <- ContactGroup(rowData = LociData, contacts = MatsList, colData = SampleData)

For this to work, the contacts list has to have the same names as the rownames of colData.

Working with bnbc contact matrices

We provide setter and getter methods for manipulating individual matrix bands for contact matrices as well. First, we have functions for working with bands of individual matrices (not bnbc related):

mat.1 <- contacts(cgEx)[[1]]
mat.1[1000:1005, 1000:1005]
b1 <- band(mat=mat.1, band.no=2)
band(mat=mat.1, band.no=2) <- b1 + 1
mat.1[1000:1005, 1000:1005]

In this example, the main diagonal of the contact matrix is also the main diagonal of the printed example above. Similarly, band number two, which is also the first off-diagonal, is also the first off-diagonal of the printed example. As can be seen from the printed example, updating a matrix band is a symmetric operation, and updated the first off-diagonal in both the upper and lower triangles of the matrix.

To utilize this across a list of contact matrices, we have the cgApply() function which applies the same function to each of the matrices. It supports parallelization using r CRANpkg("parallel").

Per-Sample Adjustments

To adjust for differences in depth sequencing, we first apply the logCPM transform [@voom] to each contact matrix. This transformation divides each contact matrix by the sum of the upper triangle of that matrix (adding 0.5 to each matrix cell and 1 to sum of the upper triangle), scales the resulting matrix by $10^6$ and finally takes the log of the scaled matrix (a fudge factor is added to both numerator and denominator prior to taking the logarithm)..

cgEx.cpm <- logCPM(cgEx)

Additionally, we smooth each contact matrix with a square smoothing kernel to reduce artifacts of the choice of bin width. We support both box and Gaussian smoothers.

cgEx.smooth <- boxSmoother(cgEx.cpm, h=5)
## or
## cgEx.smooth <- gaussSmoother(cgEx.cpm, radius=3, sigma=4)

Cross Sample Normalization

BNBC operates on each matrix band separately. For each matrix band $k$, we extract each sample's observation on that band and form a matrix $M$ from those bands; if band $k$ has $d$ entries, then after logCPM transformation, $M \in \mathbb{R}^{n \times d}$. For each such matrix, we first apply quantile normalization [@QN] to correct for distributional differences, and then ComBat [@ComBat] to correct for batch effects.

Here we will use r Biocpkg("bnbc") to do batch correction on the first 10 matrix bands, beginning with the second matrix band and ending on the eleventh.

cgEx.bnbc <- bnbc(cgEx.smooth, batch=colData(cgEx.smooth)$Batch,
                  threshold=1e7, step=4e4, nbands=11, verbose=FALSE)

sessionInfo()

sessionInfo()

References



Try the bnbc package in your browser

Any scripts or data that you put into this service are public.

bnbc documentation built on Nov. 1, 2018, 2:37 a.m.