rotPlaid error "Error in if (x.min >= x.max) { : missing value where TRUE/FALSE needed"
Entering edit mode
shopnil99 • 0
Last seen 2.1 years ago

I am following chapter 8 of the user guide to detect differential domain boundaries. My analysis pipeline is as follows:

I used the HICUP digester script to generate the restriction digest file for the Arima restriction enzyme cocktail and used this file to generate the hg19_param file.

I mapped my fastq files to hg19 genome using the Arima Hi-C mapping pipeline, which generates bam files for my two conditions each with 2 biological replicates. I name sorted the bam file using samtools sort -n and used these files to generate the corresponding index files using preparePairs.

After following all the steps in chapter 8, the summary(change.type[is.sig]) is as follows; Mode FALSE TRUE logical 1686 5534

Everything seems to work great up until this next step when I try to generate the rotated plots:

tophit <- cur.regions[o[1]] expanded <- resize(tophit, fix="center", width=width(tophit)*10) par(mfrow=c(2,1)) rotPlaid(input[1], mm.param, region=expanded, width=2.5e4, main="iLuc") segments(start(tophit), 0, end(tophit), 0, col="red", lwd=10) rotPlaid(input[3], mm.param, region=expanded, width=2.5e4, main="iEF") segments(start(tophit), 0, end(tophit), 0, col="red", lwd=10)

I get the following error at the first rotPlaid:

Error in if (x.min >= x.max) { : missing value where TRUE/FALSE needed

Any idea what this error might be and how to go about solving it?

diffhic rotPlaid • 283 views
Entering edit mode
Aaron Lun ★ 27k
Last seen 1 hour ago
The city by the bay

If I had to guess, your seqlengths(mm.param$fragments) is not properly set up. Make sure that the fragments object you supply to pairParam has proper seqlengths, i.e., chromosome lengths. This should be the case if they came out of cutGenome - the BioC-devel version also supports multiple restriction enzymes.

Entering edit mode

Thank you for your answer Aaron. Indeed, using the fragment output from HICUP Digester to create the param object leads to missing "seqlenghts" and eventually the above error. I used the devtools package in R/3.6.0 to turn on dev_mode and was able to use cutGenome for multiple restriction sites, which solved the issue.


Login before adding your answer.

Traffic: 317 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6