---
title: "`HiCcompare` Vignette"
author: "John Stansfield, Mikhail Dozmorov"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{HiCcompare Usage Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r set-options, echo=FALSE, cache=FALSE}
options(width = 400)
```
# Introduction

`HiCcompare` provides functions for the joint normalization and detection of differential chromatin interactions between two or multiple Hi-C datasets. `HiCcompare` operates on processed Hi-C data in the form of chromosome-specific chromatin interaction matrices. It accepts three-column tab-separated text files storing chromatin interaction matrices in a sparse matrix format (see [Creating the hic.table object](#hic_table)). Functions to convert popular Hi-C data formats (`.hic`, `.cool`) to sparse format are available (see ?cooler2sparse). `HiCcompare` differs from other packages that attempt to compare Hi-C data in that it works on processed data in chromatin interaction matrix format instead of pre-processed sequencing data. In addition, `HiCcompare` provides a non-parametric method for the joint normalization and removal of biases between two Hi-C datasets for the purpose of comparative analysis. `HiCcompare` also provides a simple yet robust permutation method for detecting differences between Hi-C datasets. 

The `hic_loess` function outputs normalized chromatin interactions for both matrices ready for the comparative analysis. The `hic_compare` function performs the comparative analysis and outputs genomic coordinates of pairs of regions detected as differentially interacting, interaction frequencies, the difference and the corresponding permutation p-value.

# How to use `HiCcompare`

## Installation

Install `HiCcompare` from Bioconductor.

```{r, eval = FALSE}
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("HiCcompare")
library(HiCcompare)
```


## Getting Hi-C Data

You will need processed Hi-C data in the form of sparse upper triangular matrices or BEDPE files in order to use `HiCcompare`. Data is available from several sources and two examples for downloading and extracting data are listed below. If you have full Hi-C contact matrices you can convert them to sparse upper triangular format using the full `full2sparse` function as shown in [additional functions](#addfunc)

### Extracting data from `.hic` files

Hi-C data is available from several sources and in many formats. `HiCcompare` is built to work with the sparse upper triangular matrix format popularized by the lab of Erez Lieberman-Aiden [http://aidenlab.org/data.html](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 hic.table object](#hic_table) 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](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 the `GSE63525_K562_combined_30.hic ` file from GEO [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE63525)  

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

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

This will extract the matrix from the `.hic` file and save it to the `K562.chr22.500kb.txt` text file, in the sparse upper triangular matrix format. See more examples on on how to use `straw` at [https://github.com/theaidenlab/straw/wiki/CPP#running](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 `HiCcompare` you want the raw data so you should selected `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 `HiCcompare` use `BP`. The final argument is the binsize 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 `HiCcompare`. 

From here we can just 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  `HiCcompare`, then proceed to [Creating the `hic.table` object](#hic_table).
 
### Extracting data from `.cool` files

The `cooler` software, [http://cooler.readthedocs.io/en/latest/index.html](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](ftp://cooler.csail.mit.edu/coolers) contains Hi-C data for `hg19` and `mm9` from many different sources. 

`.cool` files can be read directly into R using the `cooler2bedpe` function provided by `HiCcompare`.

```{r, eval=FALSE}
dat <- cooler2bedpe(path = "Dixon2012-H1hESC-HindIII-allreps-filtered.1000kb.cool")
```


Alternatively, data in the `.cool` format can be import into `HiCcompare` using the following steps:

  1. Download and install `cooler` from [http://cooler.readthedocs.io/en/latest/index.html](http://cooler.readthedocs.io/en/latest/index.html)
  2. Download a `.cool` file from the cooler index [ftp://cooler.csail.mit.edu/coolers](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 then proceed to [Creating the `hic.table` object](#hic_table).


### 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 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. `HiCcompare` includes a function to convert the results of HiC-Pro into a usable format for analysis in `HiCcompare`. 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 `HiCcompare`. Here we convert HiC-Pro data for input into HiCcompare:

```{r, eval=FALSE}
# read in files
mat <- read.table("hic_1000000.matrix")
bed <- read.table("hic_1000000_abs.bed")
# convert to BEDPE
dat <- 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 HiCcompare.
```


See the help using `?hicpro2bedpe` for more details.

## Detecting Copy Number Variations and Excluding regions

Before you start a `HiCcompare` analysis, you may want to detect CNV and exclude these regions along with any other regions known to exhibit changes or undesirable sequencing characteristics which could cause false positives. `HiCcompare` provides a function to perform a CNV detection analysis utilizing the `QDNAseq` R package. You will need to have your Hi-C data in `.bam` files. This can be accomplished by downloading the raw sequencing results and aligning them using one of the many Hi-C processing pipelines. Once you have `.bam` files place them in a specified folder and you can then run `get_CNV`. Make sure to specify the bin size you will be using in kilobase pairs. The bin size should be the same as the resolution of the Hi-C data that will be used in the `HiCcompare` analysis. The `CNV.level` option will allow you to choose which level of CNV you would like to exclude. The CNV calls are defined as -2 = deletion, -1 = loss, 0 = normal, 1 = gain, 2 = amplification. In order to exclude amplifications and deletions set `CNV.level = 2`. To exclude any amount of CNV set `CNV.level = 1`.

```{r, eval = FALSE}
cnv <- get_CNV(path2bam = 'path/to/bamfiles', out.file = 'path/to/bamfiles/outfile', 
               bin.size = 1000, genome = 'hg19', CNV.level = 2)
```

This function will run the CNV detection steps provided in `QDNAseq` and export a `.txt` file containing the copy number calls and a `.bed` file containing the regions detected at or above the chosen CNV level. It will also return a data.frame containing the regions at or above the chosen CNV level.  It is recommended to exclude regions with scores of -2 or 2, however this decision is left up to the user. 

You may want to exclude the regions with CNV along with blacklisted regions from any further analysis. `HiCcompare` has the ENCODE blacklists for hg19 and hg38 included (available using `data("hg19_blacklist")` or `data("hg38_blacklist")`). We can now create a data.frame (or `GenomicRanges` object) containing all of the regions we want to exclude from any further analysis.

```{r, eval = FALSE}
data('hg19_blacklist')

# combine cnv excluded regions with blacklist regions
exclude <- cbind(cnv, hg19_blacklist)
```

Now that we have a data.frame containing the regions to be excluded we just simply set the `exclude.regions` options in the `create.hic.table` function to the `exclude` data.frame.


## Creating the `hic.table` object {#hic_table}

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

If you have two full Hi-C contact matrices you can convert them to sparse upper triangular matrices using the `HiCcompare::full2sparse` function. Once you have two sparse matrices you are ready to create a `hic.table` object. This will be illustrated using the included sparse matrices at 500kb resolution for chromosome 22 from the HMEC and NHEK cell lines.

```{r, warning=FALSE, message=FALSE}
library(`HiCcompare`)
# load the data
data("HMEC.chr22")
data("NHEK.chr22")
head(HMEC.chr22)
```

Now that we have 2 sparse upper triangular matrices we can create the `hic.table` object. `create.hic.table` requires the input of 2 sparse matrices and the chromosome name.

```{r}
# create the `hic.table` object
chr22.table <- create.hic.table(HMEC.chr22, NHEK.chr22, chr = 'chr22')
head(chr22.table)
```

By default, all regions for the data entered are used as is. It may be desirable to exclude certain regions from any further analysis because they contain CNV or are overlapping blacklisted regions with known sequencing issues. To exclude the CNV and blacklisted regions we defined above, create the `hic.table` object as follows. 

```{r, eval = FALSE}
chr22.table <- create.hic.table(HMEC.chr22, NHEK.chr22, chr = 'chr22', exclude.regions = exclude, exclude.overlap = 0.2)
```

The `exclude.overlap` option controls the percentage overlap required with the regions listed in your data in order to be exclude. CNV regions should in the same resolution and bins as your data so they will overlap any of your regions by 100%. Blacklisted regions tend to be small and may only overlap your regions by a few basepairs. If you want any amount of overlap to result in exclusion set `exclude.overlap = 0`. To require 100% overlap before exclusion set `exclude.overlap = 1`. By default this option will exclude regions with 20% or more overlap.

We can also list multiple `hic.tables` in order to utilize parallel computing. First we create another `hic.table` as was done above. Then we combine these two `hic.tables` into a list.

```{r}
# create list of `hic.table` objects
data("HMEC.chr10")
data("NHEK.chr10")

# create the `hic.table` object
chr10.table <- create.hic.table(HMEC.chr10, NHEK.chr10, chr = 'chr10')
hic.list <- list(chr10.table, chr22.table)
head(hic.list)
```

The `hic.table` object contains a summary of the differences between the two matrices. "IF1" and "IF2" correspond to interaction frequencies in the first and second matrices, "D" is the unit distance (length of each unit is equivalent to the resolution of the data, e.g., 500kb), "M" is the $log_2(IF2)-log_2(IF1)$ difference.

### BEDPE format

A `hic.table` object can also be created using data in the 7 column BEDPE format. An example of BEDPE data for the HMEC dataset used above is shown below.

```{r, echo=FALSE}
HMEC.chr22_BEDPE <- chr22.table[, 1:7, with=FALSE]
NHEK.chr22_BEDPE <- chr22.table[, c(1:6, 8), with=FALSE]
head(HMEC.chr22_BEDPE)
```

To create a `hic.table` object using BEDPE data is very similar to using data in the sparse upper triangular format.

```{r}
bed.hic.tab <- create.hic.table(HMEC.chr22_BEDPE, NHEK.chr22_BEDPE)
head(bed.hic.tab)
```

### HiC-Pro format

If you are using data from HiC-Pro you can create a `hic.table` object as follows.

```{r, eval=FALSE}
# first dataset
mat1 <- read.table("hic1_1000000.matrix")
bed1 <- read.table("hic1_1000000_abs.bed")
dat1 <- hicpro2bedpe(mat, bed)
dat1 <- dat1$cis # extract intrachromosomal matrices
# second dataset
mat2 <- read.table("hic2_1000000.matrix")
bed2 <- read.table("hic2_1000000_abs.bed")
dat2 <- hicpro2bedpe(mat, bed)
dat2 <- dat2$cis # extract intrachromosomal matrices

# for chr1
hic.table <- create.hic.table(dat1[[1]], dat2[[1]])

# for all chromosomes
hic.list <- mapply(create.hic.table, dat1, dat2, SIMPLIFY = FALSE)
```


### InteractionSet/GRanges format

A `hic.table` can also be created using an `InteractionSet` object. Simply enter two `InteractionSets` representing two Hi-C matrices into the `create.hic.table` function and they will be converted to the proper format. See the `InteractionSet` vignette for creating `InteractionSet` objects [here](https://bioconductor.org/packages/release/bioc/vignettes/InteractionSet/inst/doc/interactions.html)

```{r}
data("hmec.IS")
data("nhek.IS")
head(hmec.IS)
IS.hic.tab <- create.hic.table(hmec.IS, nhek.IS)
```

To shorten computational time, or if one is only interested in a subsection of a Hi-C matrix, one may wish to use a subset of the `hic.table` object. Use the `subset.dist` or `subset.index` options, see help for the `create.hic.table` function.

## Total Sum Scaling

If you are performing your `HiCcompare` analysis on the entire genome you can perform total sum scaling on the data before performing loess normalization. This requires the use of a list of `hic.table` objects, one for each chromosome. Additionally, when you create these `hic.table` objects you must set `scale = FALSE`. Once you have your list of `hic.table` objects you can then perform total sum scaling.

```{r, eval = FALSE}
hic.list <- total_sum(hic.list)
```

It is recommended to use total sum scaling when performing a `HiCcompare` analysis on the entire genome. If you are only using data for a single chromosome then total sum scaling should be equivalent to the `scale` option in the `create.hic.table` function.

## Joint Normalization

Now that you have created a `hic.table` object you can jointly normalize your two Hi-C matrices. The `hic_loess` function has many options and can accept a single `hic.table` or a list of hic.tables. If for example you wish to perform joint normalization for every chromosome on two cell lines while utilizing parallel computing you can create a list containing the hic.tables for each chromosome. 

To change the degree of the polynomial for the loess joint normalization you can utilize the `degree` option, default is 1 (linear regression). A user-defined span, or the amount of data used to build the loess model, can also be set with the `span` option. However if `span = NA` (the default) the automatic smoothing parameter selection process will run and determine the optimal span for the data. The type of selection process can be changed using the `loess.criterion` option. Available settings are `gcv` (the default) for generalized cross-validation or `aicc` for Akaike information criterion. The loess automatic smoothing parameter selection uses a customized version of the `fANCOVA::loess.as` function. For more information on parameter selection please see the `fANCOVA` reference manual. It is recommended to use the default settings. If you have already run `hic_loess` on a dataset and know the span or would like to use a specific span then manually setting the span to a value can significantly reduce computation time especially for high resolution data.

`hic_loess` can utilize the `BiocParallel` package to perform parallel computing and lower computation time. The `parallel` option (FALSE by default) will only speed up computation if a list of hic.tables is entered into the function, i.e., it parallelizes processing of several chromosome-specific matrices. This is useful for performing joint normalization for every chromosome between two Hi-C datasets. For more information on `BiocParallel` see the reference manual [here.](http://bioconductor.org/packages/release/bioc/vignettes/BiocParallel/inst/doc/Introduction_To_BiocParallel.pdf)

The basis of `HiCcompare` rests on the novel concept termed the MD plot. The MD plot is similar to the MA plot or the Bland-Altman plot. $M$ is the log2 difference between the interaction frequencies from the two datasets. $D$ is the unit distance between the two interacting regions. Loess is performed on the data after it is represented in the MD coordinate system. To visualize the MD plot the `Plot` option can be set to TRUE. 

If you wish to detect differences between the two Hi-C datasets immediately following joint normalization you can set the `check.differences` option to TRUE. However, if you only want to jointly normalize the data for now keep this option set to FALSE. The difference detection process will be covered in the next section.  

```{r, fig.width=7, fig.height=4}
# Jointly normalize data for a single chromosome
hic.table <- hic_loess(chr22.table, Plot = TRUE, Plot.smooth = FALSE)
knitr::kable(head(hic.table))
```

```{r, message=FALSE, eval=FALSE}
# Multiple hic.tables can be processed in parallel by entering a list of hic.tables
hic.list <- hic_loess(hic.list, parallel = TRUE)
```

The `hic_loess` joint normalization function extends the `hic.table` with the adjusted interaction frequencies, adjusted "M", the "mc" correction factor, and the average expression "A".

<!-- ## Weighting of M values -->

<!-- This is an optional step of the analysis process for advanced users. This step will require the user to look at several plots to decide if weighting is necessary and if so, decide on the input parameters for the weighting. Weighting of M values may not be necessary for data at low resolutions such as 1MB or 500KB. For higher resolution data this process will probably be more necessary. However, you will need to look at the data and plots in order to determine this for yourself. If you perform a standard `HiCcompare` analysis without weighting of M values on your data and find that a large number of detected differences have low average expression it may be helpful to use M weighting. -->

<!-- M values are defined as $log2(IF2 / IF1)$. Since M values are the result of a ratio, differences on the log scale can appear identical when on the raw scale they are actually very different. For instance The M values for IF1 = 1 and IF2 = 10 and for IF1 = 10 and IF2 = 100 are both 3.32, despite the fact that the second pair of interaction frequencies has a much larger difference on the raw scale. This could potentially lead to issues in the analysis of Hi-C data, as IFs with lower average expression values are less trustworthy and more prone to display differences due solely to chance or technical biases rather than true biological differences. Therefore, M values derived from IFs with large average expression should be treated as more trustworthy than those derived from IFs exhibiting small average expression. One solution to this issue to weight M values based on their corresponding average expression before performing difference detection.  -->

<!-- To weight M values, first we calculate the average expression associated with each M value, $A = (IF1 + IF2) / 2$. Simple plotting of A values indicates that average expression between two Hi-C datasets is roughly distributed in a powerlaw fashion, with the majority of A values being rather small and a long right tail of very large A values. A values less than 1 should not be trusted at all and so they can safely be ignored for now. At this point we will need to plot the distribution of A in order to determine a minimum value of A at which to start the weighting and at which quantile of A we should stop weighting. -->

<!-- We can plot the distribution of A using the `plot_A` function on a list of hic.table objects. Here we use Hi-C data comparing two different regions of the brain for chromosomes 18 and 22 at 100KB resolution. The data has already been read into a `hic.table` list and jointly normalized with `hic_loess`.  -->

<!-- ```{r, fig.height=8, fig.width=6, warning=FALSE, message=FALSE} -->
<!-- data("brain_table") -->
<!-- aplot <- plot_A(brain_table) -->
<!-- ``` -->

<!-- From looking at the plots we can see that the density levels out by the 75th quantile so we can focus on the values of A < 75th quantile. Next we will need to determine the minimum value of A at which to start fitting the powerlaw for weighting. From the 4th plot we can see that after a value of 4 the density of A starts to more or less uniformly decrease so this should be a safe guess for the minimum A value. The peak of the log2(A) histogram looks to occur at about 1.5 so if we take the anti-log we get 2.8 further suggesting that a value between 2.8 and 4 should be okay for fitting the weights.  -->

<!-- Now we can perform weighting of the M values. -->

<!-- ```{r, warning=FALSE} -->
<!-- brain_table <- weight_M(brain_table, A.min = 4, quant = 0.75, Plot.diagnostic = TRUE, Plot.MD = TRUE) -->
<!-- ``` -->

<!-- The diagnostic plots are the histogram of raw A values between the `A.min` value and the `quant` quantile. The solid line represents the powerlaw fit to the density. You should always check these plots to make sure the fit of the weights is okay. You can also see the new MD plots after M weighting above. The diagnostic plots indicate that the parameters chosen worked for the data. You can now proceed to difference detection on your data.  -->

<!-- Weighting of the M values occurs as follows. We estimate the kernel density for the A values which fall into the range $A \le A_{min} \ and \ A < A_{quant}$. The approximate probability density function, $f_A(.)$, of the kernel density is then interpolated from the kernel density. Next, we fit a powerlaw to this subset of A. The powerlaw is defined as $P(A) = cA^{-\alpha}$. We get the best fit of the powerlaw to the kernel density by minimizing the sum of squared residuals (SSR) between the powerlaw fit and the kernel density. The weights, $w$, can be calculated as follows -->
<!-- $$w = \begin{cases} 0, \  A < A_{min} \\ -->
<!-- 1 - \frac{P(A)}{max(P(A))}, \ A_{min} \le A < A_{0.75} \\ -->
<!-- 1, \ A \ge A_{0.75} -->
<!-- \end{cases} -->
<!-- $$ -->
<!-- where $P(A)$ is the powerlaw fit based on the minimizing the SSR. This weighting scheme will produce weights on [0, 1] corresponding to the powerlaw fit to the estimated kernel density of A. The weights can then be multiplied to the normalized M values before difference detection is performed. The weights can only make the M values closer to 0 and thus the differences represented will become less significant if their corresponding A value is low. M values corresponding to high A values will be unchanged.  -->


## Difference Detection

<!-- Difference detection also makes use of the MD plot. The jointly normalized datasets are once again plotted on the MD plot. Difference detection takes place on a per-unit-distance basis using a permutation test. Permutations are broken into blocks for each unit distance. See methods section of the manuscript [here](https://doi.org/10.1101/147850) for more details. Difference detection can be performed at the same time as joint normalization if the `check.differences` option in `hic_loess` is set to TRUE. Otherwise a `hic.table` object output from the `hic_loess` function can be input into the `hic_diff` function for difference detection.  -->

<!-- `hic_diff` accepts a `hic.table` object or a list of hic.table objects output from `hic_loess`. `hic_diff` can also utilize the `BiocParallel` package the same way as the `hic_loess` loess function. The same limits and uses apply as stated above. If the `Plot` option is set to TRUE, an MD plot will be displayed showing the normalized MD plot along with coloring based on the significance of each difference. The dotted lines on this MD plot represent the `diff.thresh` threshold. This threshold is automatically calculated as 2 standard deviations of M from 0 on the MD plot by default but can be set by the user or turned off. This option helps to deal with false positives by restricting the regions that can be labelled as significant when they fall inside the threshold. The threshold can be thought of as a range in which differences do not exceed the level of the expected noise in the Hi-C data. If you used M weighting on your data you should use a `diff_thresh` otherwise you will get many false positives for longer range interactions. The function will return a hic.table or a list of hic.tables with an additional column containing the p-value for the significance of the difference between the IFs in each cell of the two Hi-C matrices.  -->

<!-- If you wish to visualize the p-values in the form of a Hi-C contact map you can use the `visualize_pvals` function to transform the results to a contact matrix with the p-values inserted instead of IFs as illustrated in the example below. The `visualize_pvals` function will multiply the -log10 p-value by 1 if the fold change for that interaction is positive or by -1 if the fold change is negative thus allowing for differentiation between regions that display either an increase or decrease in interaction frequencies. -->

<!-- ```{r, fig.width=7, fig.height=6, warning=FALSE, message = FALSE} -->
<!-- # input hic.table object into hic_diff -->
<!-- hic.table <- hic_diff(hic.table, Plot = TRUE, Plot.smooth = FALSE) -->

<!-- # input a list of hic.tables -->
<!-- hic.list <- hic_diff(hic.list) -->

<!-- # visualize results as contact map of p-values -->
<!-- visualize_pvals(hic.table) -->
<!-- ``` -->

the MD plot is used to represent the differences $M$ between two normalized datasets at a distance $D$. Normalized $M$ values are centered around 0 and the distribution at each distance remains fairly consistent. The $M$ values are approximately normally distributed. $M$ values are converted to Z-scores using the standard approach:


$$
Z_i = \frac{M_i - \bar{M}}{\sigma_M}
$$
where $\bar{M}$ is the mean value of all $M$'s on the chromosome and $\sigma_M$ is the standard deviation of all $M$ values on the chromosome and $i$ is the $i$th interacting pair on the chromosome. 

During Z-score conversion the average expression of each interacting pair is considered. Due to the nature of $M$, a difference represented by an interacting pair with IFs 1 and 10 is equivalent to an interacting pair of IFs 10 and 100 with both differences having an $M$ value of 3.32. However the average expression of these two differences is 5.5 and 55, respectively. Differences with higher average expression are more trustworthy than those with low average expression due to the fact that a relatively small difference on the raw scale can lead to a large log2 fold change. Additionally since the IFs represent Hi-C sequencing reads, lower average numbers of reads are also less trustworthy. Thus we filter out differences with low average expression before calculating Z-scores when average expression ($A$) is less than a user set quantile (typically 5-20%) or a user set specific value of $A$. The Z-scores can then be converted to p-values using the standard normal distribution.

To determine how to filter your data you can use the `filter_params()` function. This will produce a plot of the Matthews Correlation Coefficient (MCC) vs. the A minimum value filtered out. This plot should help you determine a reasonable value to set your filtering parameter in `hic_compare`. You may need to adjust the `numChanges` option and the `FC` option depending on the resolution of the data and how noisy it is.

Here we check the filtering:

```{r, message=FALSE}
filter_params(hic.table)
```


Difference detection can be performed using the `hic_compare()` function. Before difference detection you will need to decide how you would like to filter out interactions with low average expression. You can filter based on a specific value of A. Once you have a normalized `hic.table` or list of `hic.tables` detect differences with an A filter as follows:

```{r}
hic.table <- hic_compare(hic.table, A.min = 15, adjust.dist = TRUE, p.method = 'fdr', Plot = TRUE)
```

If you leave the `A.min = NA` then the 10th percentile of A will be calculated and used for filtering by default.
The resulting table will display the results as shown below:

```{r}
knitr::kable(head(hic.table))
```

Where `Z` is the Z-score calculated for the interaction (Z-scores with NA values result from filtering based on A), `p.value` is the un-adjusted p-value based on the standard normal distribution, and `p.adj` is the p-value after the multiple testing correction specified was applied.

## Converting `HiCcompare` results to `InteractionSet` objects

If after running `hic_loess` or `hic_compare` on your data you wish to perform additional analyses which require the `GRanges` object class you can convert the results using the `make_InteractionSet` function. This function will produce an `InteractionSet` object with the genomic ranges contained in the `hic.table` along with several metadata files containing the additional information produced by `hic_loess` or `hic_compare`.

```{r}
IntSet <- make_InteractionSet(hic.table)
```



## Simulating Hi-C Data

`HiCcompare` includes functions for simulating Hi-C data. The `hic_simulate` function allows you to simulate two Hi-C matrices with added bias and true differences at a specified fold change. As an example we will simulate two matrices with 250 true differences added at a fold change of 4.

```{r fig.width=7, fig.height=4}
number_of_unitdistances <- 100 # The dimensions of the square matrix to be simualted
number_of_changes       <- 250 # How many cells in the matrix will have changes
i.range <- sample(1:number_of_unitdistances, number_of_changes, replace = TRUE) # Indexes of cells to have controlled changes
j.range <- sample(1:number_of_unitdistances, number_of_changes, replace = TRUE) # Indexes of cells to have controlled changes

sim_results <- hic_simulate(nrow = number_of_unitdistances, medianIF = 50000, sdIF = 14000, powerlaw.alpha = 1.8, fold.change = 4, i.range = i.range, j.range = j.range, Plot = TRUE, alpha = 0.1)
```

The results of the simulation are saved in a list.

```{r}
names(sim_results)
```

TPR is the true positive rate, SPC is the specificity, pvals is a vector of the p-values for every cell in the matrices, hic.table is the resulting hic.table object after `hic_loess` and `hic_compare` have been applied to the simulated data, true.diff is a table for the cells that had the specified fold change applied to them, truth is a vector of 0's and 1's indicating if a cell had a true difference applied - this is useful for input into ROC packages, sim.table is the simulated data in a hic.table object before being scaled, normalized, and analyzed for differences.

The `sim_matrix` function will produce two simulated Hi-C matrices without performing any analysis.

```{r}
sims <- sim_matrix(nrow = number_of_unitdistances, medianIF = 50000, sdIF = 14000, powerlaw.alpha = 1.8, fold.change = 4, 
                   i.range = i.range, j.range = j.range)
```


## Additional Functions {#addfunc}

`HiCcompare` contains some additional functions that may be useful. 

If you do not choose to show the MD plots when you initially run `hic_loess` or `hic_compare` you can use the `MD.plot1` or `MD.plot2` functions. `MD.plot1` will create a side by side MD plot showing before and after loess normalization. Enter your original M and D vectors along with the M correction factor, `mc`, calculated by `hic_loess`. The `smooth` option controls if the plot is plotted using `smoothScatter` or as a `ggplot2` scatter plot.

```{r, fig.height=4, fig.width=7}
MD.plot1(M = hic.table$M, D = hic.table$D, mc = hic.table$mc, smooth = TRUE)
```

`MD.plot2` will create a standard MD plot with optional coloring based on p-value. Just enter an M and D vector and a p-value vector if desired.

```{r, fig.height=4, fig.width=7}
# no p-value coloring
MD.plot2(M = hic.table$adj.M, D = hic.table$D, smooth = FALSE)

# p-value coloring
MD.plot2(M = hic.table$adj.M, D = hic.table$D, hic.table$p.value, smooth = FALSE)
```

There are two matrix transformation functions included. `sparse2full` will transform a sparse upper triangular matrix to a full Hi-C matrix. `full2sparse` will transform a full Hi-C matrix to sparse upper triangular format.

```{r}
full.NHEK <- sparse2full(NHEK.chr22)
full.NHEK[1:5, 1:5]

sparse.NHEK <- full2sparse(full.NHEK)
head(sparse.NHEK)
```



`KRnorm` will perform Knight-Ruiz normalization on a Hi-C matrix. Just enter the full Hi-C matrix to be normalized.

```{r}
KR.NHEK <- KRnorm(full.NHEK)
```


`SCN` will perform Sequential Component Normalization on a Hi-C matrix. Just enter the full Hi-C matrix to be normalized.

```{r}
SCN.NHEK <- SCN(full.NHEK)
```


`MA_norm` will perform MA normalization on a `hic.table` object.

```{r, message=FALSE}
result <- MA_norm(hic.table, Plot = TRUE)
```


## Using HiCcompare on a cluster

Below is an example script for using `HiCcompare` on a computing cluster. The data was concatenated into a single sparse matrix for each group, containing information for all the chromosomes. 


```{r, eval = FALSE}
library(HiCcompare)
library(BiocParallel)

args = commandArgs(trailingOnly=TRUE)

dat1 <- read.table(args[1], header=FALSE, col.names=c("chr1", "start1", "end1", "chr2", "start2", "end2", "IF"))
dat2 <- read.table(args[2], header=FALSE, col.names=c("chr1", "start1", "end1", "chr2", "start2", "end2", "IF"))

dat1 <- dat1[dat1$chr1==dat1$chr2, ]
dat2 <- dat2[dat2$chr1==dat2$chr2, ]

dat1 <- split(dat1, dat1$chr1)
dat2 <- split(dat2, dat2$chr1)

hic.list <- mapply(create.hic.table, dat1, dat2, SIMPLIFY = FALSE, scale=FALSE)

hic.list <- total_sum(hic.list)

register(MulticoreParam(workers = 10), default = TRUE)

hic.list <- hic_loess(hic.list, Plot=TRUE, parallel=TRUE)
hic.list <- hic_compare(hic.list, A.min = NA, adjust.dist = TRUE, p.method = 'fdr', Plot = TRUE, parallel=TRUE)

hic.list <- do.call(rbind, hic.list)

hic.list <- hic.list[hic.list$p.adj<0.05,]

write.table(hic.list, args[3])
```

The above can be saved as an R script and then submitted to a cluster using the following submission script.

```{bash, eval = FALSE}
#!/bin/sh
# Sun Grid Engine options (lines prefixed with #$)
#$ -cwd
#$ -N HiCCompare
#$ -l h_rt=06:00:00
#$ -l h_vmem=32G
#$ -pe sharedmem 11
#$ -V

Rscript 01_HiCcompare.R $1 $2 $3
```


```{bash, eval = FALSE}
qsub 01_launch_HiCcompare.sh dat1.1000000.pixels.tsv dat2.1000000.pixels.tsv hiccompare_1Mb.tsv
```

Credit to Ilya Flyamer (@Phlya on GitHub) for this example of using `HiCcompare` on a cluster.