search

ay-lab / HiCtrans

HiCtrans is a pipeline to call translocations from Hi-C data
15 stars 5 forks source link

Test for changepoint Q value #9

Closed johnstonmj closed 3 years ago

johnstonmj commented 3 years ago

Hello @abhijitcbio Thanks for the HiCTrans package!

I was encountering the error:

Error in BINSEG(sumstat, pen = pen.value, cost_func = costfunc, minseglen = minseglen,  : 
  Q is larger than the maximum number of segments 94.5

I am running 25-kb resolution data for the human genome. I only encounter this error when working with chr 21.

My guess is that this is an interaction between changepoint binseg, data resolution provided, and the default '--resolutions' argument [2,4,5,10].

chr21: 46709983 bp 25000 (data resolution) 10 (largest resolution multiplier) 2 (changepioint paired comparisons) = 47500000

For the smallest chromosome, there aren't enough bins to run changepoint with Q = 100, as in:

  seg <- ifelse(resolution > 250000, 50, 100) 
  r.cpt <- suppressWarnings(cpt.meanvar(r.center,Q=seg,method="BinSeg"))
  c.cpt <- suppressWarnings(cpt.meanvar(c.center,Q=seg,method="BinSeg"))

It worked after I changed the comparison to be >= instead of > as in: seg <- ifelse(resolution >= 250000, 50, 100)

Alternatively, could a test be added to reduce Q in the case of few data points? Something like seg = min(data rows / 2, 100)

Thanks!

ay-lab commented 3 years ago

Hi, Thanks for reporting this. I have updated the code to check the segment size.