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?

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.