kundajelab / 3DChromatin_ReplicateQC

Software to compute reproducibility and quality scores for Hi-C data
MIT License
43 stars 16 forks source link

Summary tables not created when HiRep gives NaN back #7

Closed cgirardot closed 5 years ago

cgirardot commented 5 years ago

Hi,

for some reasons, HiC-rep fails for one of my resolution (50K) while everything goes fine at 10K i.e. here is the result file SMPL1.vs.SMPL2.txt in the reproducibility/HiCRep folder:

SMPL1   SMPL2   chr2L        NA
SMPL1   SMPL2   chr2R       0.961
SMPL1   SMPL2   chr3L        0.975
SMPL1   SMPL2   chr3R        0.964
SMPL1   SMPL2   chrX          0.938

Weirdly, I have 6 comparisons and it constantly fails for chr2L for all 6 comparisons (Dm data). The result above is gained using corrected counts. But I have the same with raw counts but here all car failed :

SMPL1   SMPL2   chr2L       NA
SMPL1   SMPL2   chr2R       NA
SMPL1   SMPL2   chr3L       NA
SMPL1   SMPL2   chr3R       NA
SMPL1   SMPL2   chrX         NA

I manually executed your HiCRep_wrapper.R code and this is what the get_scc() returns :

> get.scc(m1_big[,-c(1:3)], m2_big[,-c(1:3)], resol, h, 0, 100000)
$corr
logical(0)

$wei
logical(0)

$scc
     [,1]
[1,]  NaN

$std
[1] NaN

which ends up as NA and unfortunately make the generation of the summary results in the scores folder fails.

cgirardot commented 5 years ago

I just tried to run the same code with HiC-rep 1.4, which solves the issue i.e. :

Replace SCC.out = get.scc(m1_big[,-c(1:3)], m2_big[,-c(1:3)], resol, h, 0, maxdist) with

processed <- prep(m1_big, m2_big, resol, h, maxdist)
SCC.out = get.scc(processed, resol, maxdist)

and I now get (truncated for clarity) :

$corr
        1         2         3         4   
0.9657219 0.9644333 0.9631286 0.9616232
$wei
       1        2        3        4        
78.25000 78.08333 77.91667 77.75000 

$scc
          [,1]
[1,] 0.9283221

$std
[1] 0.002273394
oursu commented 5 years ago

Thanks very much for the detailed fix. I updated the install scripts to pull in the latest version of HiCRep from Bioconductor, and used your updated command.