rotPlaid error "Error in if (x.min >= x.max) { : missing value where TRUE/FALSE needed"
1
0
Entering edit mode
shopnil99 • 0
@shopnil99-20966
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
ADD COMMENT
1
Entering edit mode
Aaron Lun ★ 27k
@alun
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.

ADD COMMENT
0
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.

ADD REPLY

Login before adding your answer.

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

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

Powered by the version 2.3.6