Currently I am in the process of analyzing HiC data from two different conditions.  When I talked to some colleagues about best way of analyzing and comparing two HiC data-sets, I was asked to look at the `diffHiC`.   I read the user-manual and the paper.  Congratulations on the paper.  I have a question.  The paper starts with reads or `bam` files.  What if I already have to HiC interaction matrices from two different conditions.  Can I still use `diffHiC`?  I did not see a simple path to do that.  If I am missing something please point me in the right direction.  Any help in this regard is highly appreciated.

What I have used is `HiCPro` (the first author of this paper also suggested that I use `diffHiC`) to process the data.  Each data has ~5 B reads.  I have already aligned and generated ICE normalized contact matrices.  There is quite a bit of difference between the two interms of read-depth.  When I plotted them using `HiCPlotter`, I can see that there are differences, but I would rather have some kind of statistical measure/FDR for those differences so that I can explain to the collaborator how to interpret that data.

Yes, you can use interaction matrices, provided that the entries are counts. That is, not normalized intensities, but the raw read pair counts. If you only have intensities, that's too bad; diffHic (and, in general, most statistical software for handling sequencing data) needs to know the original count size in order to do a proper analysis. I'm not sure what else HiC-Pro produces, but I bet the relevant information is there somewhere.

Assuming you have or can get interaction matrices containing counts, you can load them into memory and convert them into ContactMatrix objects using the code below. You'll also need GRanges objects describing the genomic coordinates of the region represented by each row/column.

cm1 <- ContactMatrix(first.matrix, row.regions, column.regions)
cm2 <- ContactMatrix(second.matrix, row.regions, column.regions)

This can be converted into an InteractionSet object for entry into diffHic, like so:

to.keep <- as.matrix(cm1)!=0 | as.matrix(cm2)!=0
iset1 <- deflate(cm1, extract=to.keep)
iset2 <- deflate(cm2, extract=to.keep)
data <- cbind(iset1, iset2)

The to.keep setting basically specifies that we want to keep all entries of the interaction matrix that are non-zero for one or more of the libraries, to avoid storing all possible region pairs. The data object needs a little bit of work before it's fully ready:

interactions(data) <- as(interactions(data), "ReverseStrictGInteractions")

... and then you can plug it into the downstream diffHic functions. I should note, though, that a proper analysis needs replicates for each group. It sounds as if you only have one matrix for each condition, which would be a struggle to analyse rigorously. In any case, the code above will work with multiple ContactMatrix objects; just expand it to cm3, cm4, iset3, iset4, etc. as necessary.

Hi @Aaron, thank you.  This is useful.  Can you give me an idea of what columns should be in the `GRanges` object.  I do have raw matrices generated as part of the pipeline.  The `bin` information is as follows:

    chr1    0    2500    1
    chr1    2500    5000    2

Is this information enough?  Additionally, I do have a total of two replicates per condition.  But at 1M resolution.  I will try this analysis with that resolution first, and go from there.