2.2.1 Extracting data from .hic files

Hi-C data is available from several sources and in many formats. multiHiCcompare is built to work with the sparse upper triangular matrix format popularized by the lab of Erez Lieberman-Aiden http://aidenlab.org/data.html. If you already have Hi-C data either in the form of a sparse upper triangular matrix or a full contact matrix you can skip to the creating the hicexp object section. If you obtain data from the Aiden Lab in the .hic format you will need to first extract the matrices that you wish to compare.

1. Download the straw software from https://github.com/theaidenlab/straw/wiki and install it.
2. Use straw to extract a Hi-C sparse upper triangular matrix. An example is below:

Say we downloaded and uncompressed the GSE63525_K562_combined_30.hic file from GEO https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525, direct link ftp, http

To extract the raw matrix corresponding to chromosome 22 at the 500kb resolution we would use the following command within the terminal

./straw NONE GSE63525_K562_combined_30.hic 22 22 BP 500000 > K562.chHCT116_r22.500kb.txt

This will extract the matrix from the .hic file and save it to the K562.chHCT116_r22.500kb.txt text file, in the sparse upper triangular matrix format. See more examples on how to use straw at https://github.com/theaidenlab/straw/wiki/CPP#running. Straw requires several inputs for the extraction of data from a .hic file.

<NONE/VC/VC_SQRT/KR> <hicFile(s)> <chr1>[:x1:x2] <chr2>[:y1:y2] <BP/FRAG> <binsize>

The first argument is the normalization method. For use in multiHiCcompare you want the raw data so you should select NONE. The second argument is the .hic file name. Next is the chromosome numbers of the matrix you want. For an intrachromosomal contact map, both should be the same as in the above example. If you want a matrix of interchromosomal interactions, you can use different chromosomes, i.e. interactions between chromosome 1 and chromosome 2 (Note that HiCcompare is only meant to be used on intrachromosomal interactions at this point in development). The next argument is whether you want basepair or fragment files. For multiHiCcompare use BP. The final argument is the bin size of the matrix (the resolution). To extract a matrix at a resolution of 1MB enter 10000000. Typical bin sizes are 1MB, 500KB, 100KB, 50KB, 5KB, 1KB. Note that most matrices with resolutions higher than 100KB (i.e., matrices with resolutions of 1KB - 50KB) are typically too sparse (due to insufficient sequencing coverage) for analysis in multiHiCcompare.

From here we can import the matrix into R as you would normally for any tab-delimited file.

3. Import the data into R K562.chr22 <- read.table('K562.chr22.500kb.txt', header=FALSE)
4. Repeat these steps for any other Hi-C dataset that you wish to compare to the first dataset using multiHiCcompare.

2.2.2 Extracting data from .cool files

The cooler software, http://cooler.readthedocs.io/en/latest/index.html, allows access to a large collection of Hi-C data. The cooler index ftp://cooler.csail.mit.edu/coolers contains Hi-C data for hg19 and mm9 from many different sources. To use data in the .cool format in HiCcompare follow these steps:

1. Download and install cooler from http://cooler.readthedocs.io/en/latest/index.html
2. Download a .cool file from the cooler index ftp://cooler.csail.mit.edu/coolers.
3. Say we downloaded the Dixon2012-H1hESC-HindIII-allreps-filtered.1000kb.cool file. See cooler dump --help for data extraction options. To extract the contact matrix we use the following commands in the terminal:
   cooler dump --join Dixon2012-H1hESC-HindIII-allreps-filtered.1000kb.cool > dixon.hESC.1000kb.txt
4. Read in the text file as you would any tab-delimited file in R
   hesc1000kb <- read.table("dixon.hESC.1000kb.txt", header = FALSE)
5. Convert to a sparse upper triangular matrix using the HiCcompare::cooler2sparse function.
   sparse <- cooler2sparse(hesc1000kb)
6. Repeat the steps for another Hi-C dataset that you wish to compare to the first dataset.

2.2.3 Using data from HiC-Pro

HiC-Pro is another tool for processing raw Hi-C data into usable matrix files. HiC-Pro will produce a .matrix file and a .bed file for the data. These .matrix files are in a sparse upper triangular format similar to the results of Juicer and the dumped contents of a .hic file, however instead of using the genomic start coordinates for the first two columns of the sparse matrix they use an ID number. The .bed file contains the mappings for each of these IDs to their genomic coordinates. The original HiCcompare package includes a function to convert the results of HiC-Pro into a usable format for analysis in multiHiCcompare. When using data from HiC-Pro, it is important to use the raw .matrix files and NOT the iced .matrix files. The iced .matrix files have already had ICE normalization applied to them and are not suitable for entry into multiHiCcompare. Here we convert HiC-Pro data for input into multiHiCcompare:

# read in files
mat <- read.table("hic_1000000.matrix")
bed <- read.table("hic_1000000_abs.bed")
# convert to BEDPE
dat <- HiCcompare::hicpro2bedpe(mat, bed)
# NOTE: hicpro2bedpe returns a list of lists. 
#   The first list, dat$cis, contains the intrachromosomal contact matrices
# NOTE: dat$trans contains the interchromosomal 
#   contact matrix which is not used in multiHiCcompare.

See the help using ?HiCcompare::hicpro2bedpe for more details.

2.4.1 Sparse upper triangular format

A sparse matrix format represents a relatively compact and human-readable way to store pair-wise interactions. It is a tab-delimited text format containing three columns: “region1” - a start coordinate (in bp) of the first region, “region2” a start coordinate of the second region, and “IF” - the interaction frequency between them (IFs). Zero IFs are dropped (hence, the sparse format). Since the full matrix of chromatin interactions is symmetric, only the upper triangular portion, including the diagonal, is stored. Typically matrices in this format are stored in separate text files for each chromosome. For use in multiHiCcompare, you will need to add a column for the chromosome number. The chromosome number should be entered as just the number. Chromosomes such as X, Y, etc. should be entered as 23, 24, etc. If you are planning to analyze data for more than a single chromosome, you will need to concatenate these matrices together. A sparse Hi-C matrix ready to be input into multiHiCcompare should look like the following:

data("HCT116_r1") # load example sparse matrix
head(HCT116_r1)
#>   "22"       V1       V2 V3
#> 1   22 16000000 16000000 11
#> 2   22 16100000 16100000  1
#> 3   22 16200000 16200000  3
#> 4   22 16300000 16300000 15
#> 5   22 16400000 16400000  3
#> 6   22 16400000 16500000  1
colnames(HCT116_r1) <- c('chr', 'region1', 'region2', 'IF') # rename columns
head(HCT116_r1) # matrix ready to be input into multiHiCcompare
#>   chr  region1  region2 IF
#> 1  22 16000000 16000000 11
#> 2  22 16100000 16100000  1
#> 3  22 16200000 16200000  3
#> 4  22 16300000 16300000 15
#> 5  22 16400000 16400000  3
#> 6  22 16400000 16500000  1

If you have full Hi-C contact matrices you can convert them to sparse upper triangular matrices using the HiCcompare::full2sparse function and then add a column indicating the chromosome.

Say we have data from 2 experimental conditions with 2 samples each. We can make a hicexp object by doing the following.

data("HCT116_r1", "HCT116_r2", "HCT116_r3", "HCT116_r4")
hicexp1 <- make_hicexp(HCT116_r1, HCT116_r2, HCT116_r3, HCT116_r4, 
                       groups = c(0, 0, 1, 1), 
                       zero.p = 0.8, A.min = 5, filter = TRUE,
                       remove.regions = hg19_cyto)
hicexp1
#> Hi-C Experiment Object 
#> 2 experimental groups 
#> Group 1 has 2 samples 
#> Group 2 has 2 samples

The groups option specifies the experimental groups. You must enter a vector the length of the number of Hi-C matrices with indicators for which group each matrix belongs to. An optional covariate data.frame with rows corresponding the Hi-C matrices and columns for each additional covariate can be provided with the covariates option.

Filtering can be performed when creating a hicexp object using the zero.p and A.min options in the make_hicexp function. The zero.p option allows for filtering by the proportion of zero IFs for an interaction. The A.min allows for filtering by a minimum average IF value. These options can be used together or individually to filter your data. Filtering is important to remove interactions with lots of 0 IFs and low average expression. These interactions tend to not be very interesting and can easily become a false positive during difference detection. Additionally, removing these interactions will increase the computational speed of multiHiCcompare. If for some reason you do not want to filter the data simply set filter = FALSE.

Additionally, you can filter out specific genomic regions such as centromeres or blacklisted regions. multiHiCcompare comes with built-in regions to be filtered for hg19 and hg38 which can be accessed like so.

data("hg19_cyto")
data("hg38_cyto")

hg19_cyto
#> GRanges object with 70 ranges and 2 metadata columns:
#>        seqnames              ranges strand |      arm  feature
#>           <Rle>           <IRanges>  <Rle> | <factor> <factor>
#>    [1]        1 121500000-125000000      * |    p11.1     acen
#>    [2]        1 125000000-128900000      * |      q11     acen
#>    [3]        1 128900000-142600000      * |      q12     gvar
#>    [4]       10   38000000-40200000      * |    p11.1     acen
#>    [5]       10   40200000-42300000      * |    q11.1     acen
#>    ...      ...                 ...    ... .      ...      ...
#>   [66]       23   58100000-60600000      * |    p11.1     acen
#>   [67]       23   60600000-63000000      * |    q11.1     acen
#>   [68]       24   11600000-12500000      * |    p11.1     acen
#>   [69]       24   12500000-13400000      * |    q11.1     acen
#>   [70]       24   28800000-59373566      * |      q12     gvar
#>   -------
#>   seqinfo: 24 sequences from an unspecified genome; no seqlengths

By default, the make_hicexp object will have the remove.regions option set to use the hg19_cyto object. If your data was not aligned to hg19 or you want other regions to be removed, you can create a GenomicRanges object containing the ranges to be removed and the remove.regions option to this object.

2.4.2 The hicexp object

The hicexp S4 class has several slots which can be accessed with the accessor functions hic_table(), results(), and meta(). The hic_table slot contains the Hi-C matrix in sparse format. The first four columns are the chromosome, region1 start location, region2 start location, and unit distance. All following chromosomes represent the IFs for that interacting pair from each sample. The comparison slot is empty at creation but will be filled following use of one of the comparison functions. It contains the same first four columns as the hic_table slot, but also has the logFC - log fold change between conditions, logCPM - log counts per million, p.value, and p.adj - multiple testing corrected p-value columns which indicate the significance of the difference for each interacting pair of regions between the conditions. Access the comparison slot using results(). The metadata slot contains the data.frame of covariates for the experiment. Access the metadata slot by using meta(). The other slots are mainly for internal use, and the typical user does not need to be concerned with them.

2.5.1 Library scaling

The simplest form of normalization to compare Hi-C data is library scaling. multiHiCcompare provides the hic_scale() function to scale the Hi-C libraries from each sample to the size of the smallest library. If you believe that any trends present in your data are important differences and not due to bias, then you can use library scaling for normalizing your data as follows.

data("hicexp2")
hicexp2 <- hic_scale(hicexp2)

Note that you need to use either simple scaling or loess normalization method. It is recommended to use either cyclic loess or fast loess that will implicitly rescale the libraries and remove unwanted trends.

2.5.2 Cyclic Loess Normalization

multiHiCcompare provides a cyclic loess method for the joint normalization of multiple Hi-C datasets. The method is based on representing the data on an MD plot. The MD plot is similar to the MA plot (Bland-Altman plot) which is commonly used for the visualization of gene expression differences. \(M\) is defined as the log difference between the two data sets \(M = log_2(IF_2/IF_1)\), where \(IF_1\) and \(IF_2\) are interaction frequencies of the first and the second Hi-C datasets, respectively. \(D\) is defined as the distance between two interacting regions, expressed in unit-length of the \(X\) resolution of the Hi-C data. A loess regression curve is fit through the MD plot and used to remove global biases by centering the \(M\) differences around \(M=0\) baseline.

The cyclic loess algorithm proceeds through the following steps.

1. Choose two out of the \(N\) total samples then generate an MD plot.
2. Fit a loess curve \(f(d)\) to the MD plot.
3. Subtract \(f(d)/2\) from the first dataset and add \(f(d)/2\) to the second.
4. Repeat until all unique pairs have been compared.
5. Repeat until convergence.

To perform cyclic loess on your Hi-C data you will need to use the cyclic_loess() function as shown below:

hicexp1 <- cyclic_loess(hicexp1, verbose = FALSE, 
                        parallel = FALSE, span = 0.2)
# make MD plot
MD_hicexp(hicexp1)

As can be seen in the above MD plots, the data for each sample has been jointly normalized with all other samples. Note that the user can set the span option. A user-set span will run quicker than the default option of automatically calculating the span. It is best to use the automatic span calculation if you have not worked with the data before, but if you are familiar with it then setting the span is a way to speed up processing. The hic_table slot in the hicexp object has also been updated with the normalized IFs.

hic_table(hicexp1)
#>        chr  region1  region2 D        IF1        IF2        IF3        IF4
#>     1:  22 18000000 18000000 0 4848.93586 5117.56104 5757.67702 4330.37698
#>     2:  22 18000000 18100000 1 1307.53629 1128.57032 1317.91263  893.55501
#>     3:  22 18000000 18200000 2  715.99832  734.05089  729.32691  743.67190
#>     4:  22 18000000 18300000 3  352.17410  409.94247  420.22430  387.50225
#>     5:  22 18000000 18400000 4  296.77137  274.80602  306.21630  314.82774
#>    ---                                                                    
#> 43740:  22 51000000 51100000 1 1622.03802 1762.87853 1664.49400 1271.16626
#> 43741:  22 51000000 51200000 2   30.51641   35.06727   31.06704   22.70308
#> 43742:  22 51100000 51100000 0 3681.72778 4131.14202 4126.33930 3604.34388
#> 43743:  22 51100000 51200000 1   79.15003   63.50592   77.51615   57.56003
#> 43744:  22 51200000 51200000 0   17.19701   25.33298   28.78896   36.23247

The runtime of cyclic loess can be decreased when multiple processors are available by setting the parallel option to TRUE. This option splits up the data by chromosome and sends each chromosome’s data to a parallel processor.

2.5.3 Fast Loess Normalization (Fastlo)

In addition to the standard cyclic loess method, multiHiCcompare also implements the Fast Loess (Fastlo) joint normalization algorithm. Our implementation of fastlo is adapted to Hi-C data on a per-distance basis. To perform “fastlo” on Hi-C data we first split the data into \(p\) pooled matrices. The “progressive pooling” is used to split up the Hi-C matrix by unit distance such that distance 0 is its own pool, distances 1 and 2 are pooled, distance 3, 4, 5 are pooled, and so on until all unit distances belong to one of \(p\) pools. Each matrix will have an \(IF_{gj}\) value with \(g\) interacting pairs for each of the \(j\) samples. These \(p\) matrices can then be input into the “fastlo” algorithm using the following steps.

1. Create the vector \(\hat{IF}_{pgj}\), the row means of the \(p^{th}\) matrix. This is the equivalent of creating an average IF at distance pool \(p\).
2. Plot \(\hat{IF}_{p}\) versus \((IF_{pg} - \hat{IF_p})\) for each sample \(j\). This is equivalent to an MA plot at a genomic distance pool \(p\).
3. Fit a loess curve \(f(x)\) to the plot.
4. Subtract \(f(x)\) from sample \(j\).
5. Repeat for all remaining replicates.
6. Repeat until algorithm converges.

You can perform fastlo normalization on your data as follows:

data("hicexp2")
# perform fastlo normalization
hicexp2 <- fastlo(hicexp2, verbose = FALSE, parallel = FALSE)
# make MD plot
MD_hicexp(hicexp2)

Again, the above MD plots show the normalized data. fastlo() can also make use of parallelization to speed up computation speeds by setting the parallel option to TRUE. The results of fastlo() and cyclic_loess() may be slightly different, but both should result in the removal of biases between Hi-C datasets. fastlo() will have quicker run times compared to cyclic_loess(), but cyclic_loess() will likely give a slightly better normalization.

2.6.1 Exact Test

For simple Hi-C experiments the hic_exactTest() function can be used as shown below:

hicexp1 <- hic_exactTest(hicexp1, p.method = 'fdr', 
                         parallel = FALSE)

# plot results
MD_composite(hicexp1)

The above composite MD plot displays where any significant differences are detected between the two groups. This function can also be sped up by using the parallel option. The results of the comparison are then saved in the comparison slot of the hicexp object.

results(hicexp1)
#>        chr  region1  region2 D        logFC    logCPM    p.value     p.adj
#>     1:  22 18000000 18000000 0 -0.003277066 11.162876 0.98355590 0.9978994
#>     2:  22 18000000 18100000 1 -0.240244231  9.056102 0.24699655 0.6200377
#>     3:  22 18000000 18200000 2  0.021080246  8.394195 0.91509193 0.9704410
#>     4:  22 18000000 18300000 3  0.092429543 10.298587 0.60748736 0.9685625
#>     5:  22 18000000 18400000 4  0.142103770  9.903732 0.45349569 0.9206592
#>    ---                                                                    
#> 43740:  22 51000000 51100000 1 -0.322643637  9.499253 0.07028879 0.3629536
#> 43741:  22 51000000 51200000 2 -0.432487991  3.866202 0.34735616 0.6875078
#> 43742:  22 51100000 51100000 0 -0.039714461 10.798143 0.80184233 0.9212657
#> 43743:  22 51100000 51200000 1 -0.140516043  5.032406 0.69220609 0.8790205
#> 43744:  22 51200000 51200000 0  0.869666933  3.727353 0.08237680 0.3831598

In this data.table the first 3 columns represent the identity of the interaction, then followed by the unit genomic distance (D), the log fold change difference between the groups (logFC), the log counts per million for the interaction (logCPM), the p-value and finally the multiple testing correction p-value (p.adj). The type of multiple testing applied can be changed using the p.method option. To view what other adjustment methods are available look at the help for ?p.adjust.

2.6.2 GLM Methods

For more complex Hi-C experiments the hic_glm() function must be used. To use this function, a design matrix must first be created. Here use the hicexp2 object with some covariates and create the design matrix.

batch <- c(1,2,1,2)
# produce design matrix
d <- model.matrix(~factor(meta(hicexp2)$group) + factor(batch))

The design matrix should contain the covariates of interest. Any categorical variables should be entered as factors. Next, the comparison of interest will need to be specified using either the contrast or the coef option. For this example we are interested in the group difference thus we can set coef = 2 to test if the group effect is equal to 0. For more information on using contrast and coef please see the edgeR user manual. Now we are ready to perform the comparison.

hicexp2 <- hic_glm(hicexp2, design = d, coef = 2, method = "QLFTest", p.method = "fdr", parallel = FALSE)

There are three methods by which hic_glm() can be performed. The default method is to the use the QLFTest which makes use of the quasi-likelihood model. Additionally, there is the LRTest which conducts a likelihood ratio test. The final method is the Treat method which conducts a test relative to a specified fold change threshold. For this option, the M option will need to be used to specify the log2 fold change threshold.

# use Treat option
hicexp2 <- hic_glm(hicexp2, design = d, coef = 2, method = "Treat", 
                  M = 0.5, p.method = "fdr", parallel = FALSE)