Errors in the getNAD.R script and issues with the wf

[bfosso]

Hi,

I'm trying to replicate your workflow but actually I've found some issue:

1. In the Bedgraph files normalization it is not clear what this formula means: mat$V6 <- as.numeric(mat$V6)/sum(as.numeric(mat$V6)) * 100000000. Actually, bedgraphs have only 4 columns and not 6. I supposed: mat$V6 <- as.numeric(mat$V4)/sum(as.numeric(mat$V4)) * 100000000. Is that correctly?
2. in the HMM calling NADs step I got this error: Error in getClass(Class, where = topenv(parent.frame())) : “TDist” is not a defined class. It seems you used the tileHMM package but I cannot find any import of this library in your code. Could you please better specify this step.
3. Finally, even when the tileHMM library is installed and loaded the getNAD function stucks at the viterbiEMstep.

Thank you in advance

[ChengLiLab]

Sorry for the delay.

1. The total column of the read count result depends on the hg19_1kb.bed file, if there are 3 column in this file, then the read count column is the 4th column. If there are 5 column in this file (including bin length and strand information), then the read count column is the 6th column.

2. Yes, we used the tileHMM package. Sorry, we missed the code for importing the library and we have added it now. Thanks for your reminding !

3. We have tested the code and it worked well with our datasets. We added a test dataset (mat.rds) for step5 HMM calling NADs. You can try to test with it.

[rivendele]

Hi, I also encountered a problem when running getNAD.R using the provided test dateset mat.rds.
after running getNAD(wgmat[["chr1"]]),
I got this error
`Error in length(init) == 0 || init@alpha != states: 'length = 2' in coercion to 'logical(1)'
Traceback:

1. getNAD(wgmat[["chr1"]])
2. hmm.setup.our(obs, state = c("enriched", "non-enriched"), pos.state = 1,
   . probe.region = 1000, frag.size = 3000)
3. new("contHMM", transition = transition, emission = emission,
   . init = pi)
4. initialize(value, ...)
5. initialize(value, ...)
6. .local(.Object, ...)`

it seems the problem with the line hmm.init <- new("contHMM", transition = transition, emission = emission,init=pi) in hmm.setup.our function