2.1Installation

BiocManager::install("TADCompare")

library(dplyr)

## ## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats': ## ## filter, lag

## The following objects are masked from 'package:base': ## ## intersect, setdiff, setequal, union

library(SpectralTAD) library(TADCompare)

3.1Working with\(n \times n\)matrices

\(n \times n\)contact matrices are most commonly associated with data coming from the Bing Ren lab (http://chromosome.sdsc.edu/mouse/hi-c/download.html). These contact matrices are square and symmetric with entry\(ij\)corresponding to the number of contacts between region\(i\)and region\(j\). Below is an example of a\(5 \times 5\)region of an\(n \times n\)contact matrix derived from Rao et al. 2014 data, GM12878 cell line(Rao et al. 2014), chromosome 22, 50kb resolution. Note the symmetry around the diagonal - the typical shape of chromatin interaction matrix. The figure was created using thepheatmappackage.

3.2Working with\(n \times (n+3)\)matrices

\(n \times (n+3)\)matrices are commonly associated with theTopDomTAD caller (http://zhoulab.usc.edu/TopDom/). These matrices consist of an\(n \times n\)matrix but with three additional leading columns containing the chromosome, the start of the region and the end of the region. Regions in this case are determined by the resolution of the data. The subset of a typical\(n \times (n+3)\)matrix is shown below.

## chr start end X18500000 X18550000 X18600000 X18650000 ## 1 chr22 18500000 18550000 13313 4817 1664 96 ## 2 chr22 18550000 18600000 4817 15500 5120 178 ## 3 chr22 18600000 18650000 1664 5120 11242 316 ## 4 chr22 18650000 18700000 96 178 316 162

3.3Working with sparse 3-column matrices

Sparse 3-column matrices are matrices where the first and second columns refer to region\(i\)and region\(j\)of the chromosome, and the third column is the number of contacts between them. This style is becoming increasingly popular and is associated with raw data from Lieberman-Aiden lab (e.g.,https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525), and is the data output produced by the Juicer tool(Durand et al. 2016). 3-column matrices are handled internally in the package by converting them to\(n \times n\)matrices using theHiCcomparepackage’ssparse2full()function. The first 5 rows of a typical sparse 3-column matrix are shown below.

## region1 region2 IF ## 1: 16050000 16050000 12 ## 2: 16200000 16200000 4 ## 3: 16150000 16300000 1 ## 4: 16200000 16300000 1 ## 5: 16250000 16300000 1 ## 6: 16300000 16300000 10

3.4Working with other data types

3.5Working with .hic files

.hic files are a common form of files generally associated with the lab of Erez Lieberman-Aiden (http://aidenlab.org/data.html). To use .hic files you must use the following steps.

1. Downloadstrawfromhttps://github.com/aidenlab/straw/and follow instalation instructions.
2. Download .hic data files. Here, we use data from Rao 2014 and download them using the following commands:

wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE63nnn/GSE63525/suppl/GSE63525_GM12878_insitu_primary_30.hic

wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE63nnn/GSE63525/suppl/GSE63525_GM12878_insitu_replicate.hic

3. Extract chromosome 22 at 50kb resolution with no normalization:

./straw NONE GSE63525_GM12878_insitu_primary_30.hic 22 22 BP 50000 > primary.chr22.50kb.txt

./straw NONE GSE63525_GM12878_insitu_replicate_30.hic 22 22 BP 50000 > replicate.chr22.50kb.txt

4. Analyze normally:

#Read in data primary = read.table('primary.chr22.50kb.txt', header = FALSE) replicate = read.table('replicate.chr22.50kb.txt', header = FALSE) #Run TADCompare tad_diff=TADCompare(primary, replicate, resolution=50000)

3.6Working with .cool files

Users can also find TADs from data output bycooler(http://cooler.readthedocs.io/en/latest/index.html) and HiC-Pro (https://github.com/nservant/HiC-Pro) with minor pre-processing using theHiCcomparepackage.

The cooler software can be downloaded fromhttps://mirnylab.github.io/cooler/. A catalog of popular HiC datasets can be found atftp://cooler.csail.mit.edu/coolers. We can extract chromatin interaction data from .cool files using the following steps:

1. Follow instructions to install the cooler software,https://mirnylab.github.io/cooler/
2. Download the first contact matrix wgetftp://cooler.csail.mit.edu/coolers/hg19/Zuin2014-HEK293CtcfControl-HindIII-allreps-filtered.50kb.cool
3. Convert the first matrix to a text file usingcooler dump --join Zuin2014-HEK293CtcfControl-HindIII-allreps-filtered.50kb.cool > Zuin.HEK293.50kb.Control.txt
4. Download the second contact matrix wgetftp://cooler.csail.mit.edu/coolers/hg19/Zuin2014-HEK293CtcfDepleted-HindIII-allreps-filtered.50kb.cool
5. Convert the matrix to a text file usingcooler dump --join Zuin2014-HEK293CtcfDepleted-HindIII-allreps-filtered.50kb.cool > Zuin.HEK293.50kb.Depleted.txt
6. Run the code below

# Read in data cool_mat1 <- read.table("Zuin.HEK293.50kb.Control.txt") cool_mat2 <- read.table("Zuin.HEK293.50kb.Depleted.txt") # Convert to sparse 3-column matrix using cooler2sparse from HiCcompare sparse_mat1 <- HiCcompare::cooler2sparse(cool_mat1) sparse_mat2 <- HiCcompare::cooler2sparse(cool_mat2) # Run TADCompare diff_tads = lapply(names(sparse_mat1), function(x) { TADCompare(sparse_mat1[[x]], sparse_mat2[[x]], resolution = 50000) })

3.7Working with HiC-Pro files

HiC-Pro data is represented as two files, the.matrixfile and the.bedfile. The.bedfile contains four columns (chromosome, start, end, ID). The.matrixfile is a three-column matrix where the 1^stand 2^ndcolumns contain region IDs that map back to the coordinates in the bed file, and the third column contains the number of contacts between the two regions. In this example we analyze two matrix filessample1_100000.matrixandsample2_100000.matrixand their corresponding bed filessample1_100000_abs.bedandsample2_100000_abs.bed. We do not include HiC-Pro data in the package, so these serve as placeholders for the traditional files output by HiC-Pro. The steps for analyzing these files is shown below:

# Read in both files mat1 <- read.table("sample1_100000.matrix") bed1 <- read.table("sample1_100000_abs.bed") # Matrix 2 mat2 <- read.table("sample2_100000.matrix") bed2 <- read.table("sample2_100000_abs.bed") # Convert to modified bed format sparse_mats1 <- HiCcompare::hicpro2bedpe(mat1,bed1) sparse_mats2 <- HiCcompare::hicpro2bedpe(mat2,bed2) # Remove empty matrices if necessary # sparse_mats$cis = sparse_mats$cis[sapply(sparse_mats, nrow) != 0] # Go through all pairwise chromosomes and run TADCompare sparse_tads = lapply(1:length(sparse_mats1$cis), function(z) { x <- sparse_mats1$cis[[z]] y <- sparse_mats2$cis[[z]] #Pull out chromosome chr <- x[, 1][1] #Subset to make three column matrix x <- x[, c(2, 5, 7)] y <- y[, c(2, 5, 7)] #Run SpectralTAD comp <- TADCompare(x, y, resolution = 100000) return(list(comp, chr)) }) # Pull out differential TAD results diff_res <- lapply(sparse_tads, function(x) x$comp) # Pull out chromosomes chr <- lapply(sparse_tads, function(x) x$chr) # Name list by corresponding chr names(diff_res) <- chr

3.8Effect of matrix type on runtime

The type of matrix input into the algorithm can affect runtimes for the algorithm.\(n \times n\)matrices require no conversion and are the fastest. Meanwhile,\(n \times (n+3)\)matrices take slightly longer to run due to the need to remove the first 3 columns. Sparse 3-column matrices have the highest runtimes due to the complexity of converting them to an\(n \times n\)matrix. The times are summarized below, holding all other parameters constant.

库(微基准测试)#阅读第二马trix data("rao_chr22_rep") # Converting to sparse prim_sparse <- HiCcompare::full2sparse(rao_chr22_prim) rep_sparse <- HiCcompare::full2sparse(rao_chr22_rep) # Converting to nxn+3 # Primary prim_n_n_3 <- data.frame(chr = "chr22", start = as.numeric(colnames(rao_chr22_prim)), end = as.numeric(colnames(rao_chr22_prim))+50000, rao_chr22_prim) # Replicate rep_n_n_3 <- data.frame(chr = "chr22", start = as.numeric(colnames(rao_chr22_rep)), end = as.numeric(colnames(rao_chr22_rep))+50000, rao_chr22_rep) # Defining each function # Sparse sparse <- TADCompare(cont_mat1 = prim_sparse, cont_mat2 = rep_sparse, resolution = 50000) # NxN n_by_n <- TADCompare(cont_mat1 = prim_sparse, cont_mat2 = rep_sparse, resolution = 50000) # Nx(N+3) n_by_n_3 <- TADCompare(cont_mat1 = prim_n_n_3, cont_mat2 = rep_n_n_3, resolution = 50000) # Benchmarking different parameters bench <- microbenchmark( # Sparse sparse <- TADCompare(cont_mat1 = prim_sparse, cont_mat2 = rep_sparse, resolution = 50000), # NxN n_by_n <- TADCompare(cont_mat1 = rao_chr22_prim, cont_mat2 = rao_chr22_rep, resolution = 50000), # Nx(N+3) n_by_n_3 <- TADCompare(cont_mat1 = prim_n_n_3, cont_mat2 = rep_n_n_3, resolution = 50000), times = 5, unit = "s" ) summary_bench <- summary(bench) %>% dplyr::select(mean, median) rownames(summary_bench) <- c("sparse", "n_by_n", "n_by_n_3") summary_bench

## mean median ## sparse 0.3675062 0.2253169 ## n_by_n 0.1239759 0.1239285 ## n_by_n_3 0.1470194 0.1423968

The table above shows the mean and median of runtimes for different types of contact matrices measured in seconds. As we see,TADCompareis extremely fast irrespectively of the parameters. However, sparse matrix inputs will slow down the algorithm. This can become more apparent as the size of the contact matrices increase.