# The bnbc User’s Guide

Kipper Fletez-Brant and Kasper Daniel Hansen

#### 27 October 2020

#### Abstract

A comprehensive guide to using the bnbc package for normalizing Hi-C replicates.

#### Package

bnbc 1.12.0

# 1 Introduction

The bnbc package provides functionality to perform normalization and batch correction across samples on data obtained from Hi-C (Lieberman-Aiden et al. 2009) 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.

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).

## 1.1 Citing bnbc

If you use this package, please cite (Fletez-Brant et al. 2017).

## 1.2 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. (Yang et al. 2017)), 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”.

## 1.3 Dependencies

This document has the following dependencies

    library(bnbc)

# 2 The ContactGroup class from bnbc

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:

- `rowData`: a `GRanges` object that has 1 element for each bin of the partitioned genome.
- `colData`: a `DataFrame` object that contains information on each sample (i.e. gender).
- `contacts`: a `list` of contact matrices.

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

    ## Object of class `ContactGroup` representing contact matrices with
    ##  1282 bins
    ##  40 kb average width per bin
    ##  14 samples

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

    ## Object of class `ContactGroup` representing contact matrices with
    ##  1282 bins
    ##  40 kb average width per bin
    ##  14 samples

## 2.1 Alternatives

The 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 Rcpp-based routines for getting and setting bands of normal matrices, which means `ContactGroup` is a pretty fast for our purposes.

# 3 Getting your data into bnbc

To get data into 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

    ##      [,1] [,2] [,3]
    ## [1,]    1    4    7
    ## [2,]    0    5    8
    ## [3,]    0    0    9

    ## Now we fill in the lower triangular matrix with the upper triangular
    mat[lower.tri(mat)] <- mat[upper.tri(mat)]
    mat

    ##      [,1] [,2] [,3]
    ## [1,]    1    4    7
    ## [2,]    4    5    8
    ## [3,]    7    8    9

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`.

## 3.1 Getting your data out of `*.cooler` files

The `.cooler` file format is widely adopted and supported by `bnbc`. We assume a simple cooler file format (see `?getChrCGFromCools` for a full description; importantly, we assume the same interactions are observed in all samples, even if some have a value of 0) of one resolution per file, generated by the `cooler` program. Our point of entry is to catalog which interactions are stored in the cooler file. We do this by generating an index of the positions of the file, using the function `getGenomeIdx()`.

    coolerDir <- system.file("cooler", package = "bnbc")
    cools <- list.files(coolerDir, pattern="cool$", full.names=TRUE)

    step <- 4e4

    bin.ixns.list <- bnbc:::getGenomeIdx(cools[1], step)
    bin.ixns <- bin.ixns.list$bin.ixns

We have, as output, a list `bin.ixns.list`, of two elements. The first is the element `bin.ixns`, which is a `data.table` object that lists the set of interactions observed in the cooler file. The other element is a `data.frame` object `bins`, which lists the set of genomic bins and their coordinates. These are a convenience used in some functions and generally are not of interest to the end user.

With our index, we can proceed to load our data into memory, one chromosome’s data at a time (at this time our method does not handle -interactions). We emphasize that with all observations from interactions between loci on one chromosome in memory, our algorithm is extremely efficient, with custom routines for matrix updating, and requires only pass over the data.

    data(cgEx)
    cool.cg <- bnbc:::getChrCGFromCools(bin.ixns,
                                        files = cools,
                                        chr = "chr22",
                                        colData = colData(cgEx)[1:2,])
    all.equal(contacts(cgEx)[[1]], contacts(cool.cg)[[1]])

    ## [1] TRUE

In this example, we load the `ContactGroup` object `cgEx` into memory to compare with the representation of it in `cool` files generated by the `cooler` program. We then use the method `getChrCGFromCools()` to load an entire chromosome’s interaction matrices (observed on all subjects) into memory. At this point, users have a valid `ContactGroup` object, and can proceed with their analyses as described in subsequent sections.

# 4 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]

    ##      [,1] [,2] [,3] [,4] [,5] [,6]
    ## [1,]    5   30   36   32   53   20
    ## [2,]   30    7   36   38   50   16
    ## [3,]   36   36    3   44   51   15
    ## [4,]   32   38   44    4   89   39
    ## [5,]   53   50   51   89   36   55
    ## [6,]   20   16   15   39   55    5

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

    ##      [,1] [,2] [,3] [,4] [,5] [,6]
    ## [1,]    5   31   36   32   53   20
    ## [2,]   31    7   37   38   50   16
    ## [3,]   36   37    3   45   51   15
    ## [4,]   32   38   45    4   90   39
    ## [5,]   53   50   51   90   36   56
    ## [6,]   20   16   15   39   56    5

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 parallel.

# 5 Per-Sample Adjustments

To adjust for differences in depth sequencing, we first apply the `logCPM` transform (Law et al. 2014) 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)

# 6 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 (Bolstad et al. 2003) to correct for distributional differences, and then ComBat (Johnson, Li, and Rabinovic 2007) to correct for batch effects.

Here we will use 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)

    ## Found 76 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 94 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 102 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 78 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 66 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 77 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 64 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 79 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 80 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
    ## Found 74 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.

# 7 sessionInfo()

    ## R version 4.0.3 (2020-10-10)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats4    parallel  stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] bnbc_1.12.0                 SummarizedExperiment_1.20.0
    ##  [3] Biobase_2.50.0              GenomicRanges_1.42.0
    ##  [5] GenomeInfoDb_1.26.0         IRanges_2.24.0
    ##  [7] S4Vectors_0.28.0            MatrixGenerics_1.2.0
    ##  [9] matrixStats_0.57.0          BiocGenerics_0.36.0
    ## [11] BiocStyle_2.18.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_1.0.5             locfit_1.5-9.4         lattice_0.20-41
    ##  [4] fftwtools_0.9-9        png_0.1-7              digest_0.6.27
    ##  [7] R6_2.4.1               tiff_0.1-5             RSQLite_2.2.1
    ## [10] evaluate_0.14          sva_3.38.0             httr_1.4.2
    ## [13] zlibbioc_1.36.0        rlang_0.4.8            data.table_1.13.2
    ## [16] annotate_1.68.0        blob_1.2.1             Matrix_1.2-18
    ## [19] preprocessCore_1.52.0  rmarkdown_2.5          splines_4.0.3
    ## [22] BiocParallel_1.24.0    stringr_1.4.0          htmlwidgets_1.5.2
    ## [25] RCurl_1.98-1.2         bit_4.0.4              DelayedArray_0.16.0
    ## [28] compiler_4.0.3         xfun_0.18              mgcv_1.8-33
    ## [31] htmltools_0.5.0        GenomeInfoDbData_1.2.4 bookdown_0.21
    ## [34] edgeR_3.32.0           XML_3.99-0.5           bitops_1.0-6
    ## [37] rhdf5filters_1.2.0     grid_4.0.3             nlme_3.1-150
    ## [40] xtable_1.8-4           DBI_1.1.0              magrittr_1.5
    ## [43] stringi_1.5.3          XVector_0.30.0         genefilter_1.72.0
    ## [46] limma_3.46.0           vctrs_0.3.4            EBImage_4.32.0
    ## [49] Rhdf5lib_1.12.0        tools_4.0.3            bit64_4.0.5
    ## [52] jpeg_0.1-8.1           abind_1.4-5            survival_3.2-7
    ## [55] yaml_2.2.1             AnnotationDbi_1.52.0   rhdf5_2.34.0
    ## [58] BiocManager_1.30.10    memoise_1.1.0          knitr_1.30

# References

Bolstad, B M, R A Irizarry, M Astrand, and T P Speed. 2003. “A Comparison of Normalization Methods for High Density Oligonucleotide Array Data Based on Variance and Bias.” Bioinformatics 19:185–93. https://doi.org/10.1093/bioinformatics/19.2.185.

Fletez-Brant, Kipper, Yunjiang Qiu, David U Gorkin, Ming Hu, and Kasper D Hansen. 2017. “Removing Unwanted Variation Between Samples in Hi-C Experiments.” bioRxiv, 214361. https://doi.org/10.1101/214361.

Johnson, W Evan, Cheng Li, and Ariel Rabinovic. 2007. “Adjusting Batch Effects in Microarray Expression Data Using Empirical Bayes Methods” 8:118–27. https://doi.org/10.1093/biostatistics/kxj037.

Law, Charity W, Yunshun Chen, Wei Shi, and Gordon K Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.” Genome Biology 15:R29. https://doi.org/10.1186/gb-2014-15-2-r29.

Lieberman-Aiden, Erez, Nynke L van Berkum, Louise Williams, Maxim Imakaev, Tobias Ragoczy, Agnes Telling, Ido Amit, et al. 2009. “Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome.” Science 326:289–93. https://doi.org/10.1126/science.1181369.

Yang, Tao, Feipeng Zhang, Galip Gurkan Yardimci, Fan Song, Ross C Hardison, William Stafford Noble, Feng Yue, and Qunhua Li. 2017. “HiCRep: Assessing the Reproducibility of Hi-C Data Using a Stratum- Adjusted Correlation Coefficient.” Genome Research. https://doi.org/10.1101/gr.220640.117.