<p>I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome <code>chrUn_gl000228</code>. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the <code>exomecn.mops(X)</code> was failing. </p>

<p>So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on <code>chrUn_gl000228</code> and confirmed that the program worked as expected and didn't throw any errors.</p>

CC BY-SA 4.0
2019-03-06T17:01:53.317Z
2019-05-31T08:35:59.723Z
2019-05-31T08:35:59.723Z
terdon
<p>I am trying to use <a href="https://www.bioinf.jku.at/software/cnmops/cnmops.html" rel="nofollow noreferrer">cn.MOPS</a> to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:</p>

<pre><code>source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles &lt;- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments &lt;- read.table("region.bed",sep="\t",as.is=TRUE)
gr &lt;- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X &lt;- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS &lt;- exomecn.mops(X)
</code></pre>

<p>The last command, <code>resCNMOPS &lt;- exomecn.mops(X)</code> fails with:</p>

<pre><code>Error in if (all(segMedianT == 0)) { : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments
</code></pre>

<p>If I run <code>traceback</code> immediately afterwards, I get:</p>

<pre><code>&gt; traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth, 
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact, 
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu, 
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor, 
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold, 
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount, 
       useMedian = useMedian, returnPosterior = returnPosterior, 
       ...)
1: exomecn.mops(X)
</code></pre>

<p>I am using R version <code>3.4.4</code> on an Ubuntu 16.04.5 LTS machine.</p>

<p>My bed file only has one region which is on <code>chrUn_gl000228</code>:</p>

<pre><code>$ grep chrUn_gl000228 region.bed 
chrUn_gl000228  112604  113879
</code></pre>

<p>If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?</p>

<p>Presumably, once I have built these two objects:</p>

<pre><code>segments &lt;- read.table("region.bed",sep="\t",as.is=TRUE)
gr &lt;- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X &lt;- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
</code></pre>

<p>I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the <code>GRanges</code> object. But how can I do this? </p>

CC BY-SA 4.0
2019-03-06T11:42:01.447Z
You start your question with `biocLite`, which indicates that this is not the latest version of R and Bioconductor. Which versions are you using? Have you tried with the most up to date version? Why does it complain that you have low number of segments? Could you post how many segments do you have on `X`?
llrs
CC BY-SA 4.0
2019-03-06T11:43:49.293Z
@llrs no, I haven't used the latest, but upgrading is not an option (this is a large server). I would be happy to post that, but can you give me a command for it? My R knowledge is extremely limited. That sort of command is exactly the kind of thing I am asking about. I have no idea how to inspect and get statistics about the contents of a GRanges object.
terdon
CC BY-SA 4.0
2019-03-06T11:44:48.520Z
Also, hang on, why are you saying it isn't the latest version? I got that code from the documentation of cn.MOPS. Is there a better way to load `biocLite`?
terdon
CC BY-SA 4.0
2019-03-06T11:54:12.417Z
Well you will need to use `sessionInfo()` and paste here the result to see the version of R and Bioconductor. The most up to date documentation of cn.MOPS is on the bioconductor page AFAIK (if packages are not up to the standards they are kicked out of the project). Try with `summary(X)` and `dim(X)` or simply print the object and post the first lines
llrs
CC BY-SA 4.0
2019-03-06T11:55:58.023Z
I said you weren't using the latest version because in the latest version the recommended way to install packages changed to `BiocManager::install("cn.Mops")` you need a package in CRAN to install packages of Bioconductor
llrs
CC BY-SA 4.0
2019-03-06T11:57:27.323Z
@llrs yes, I added the R version (3.4.4). I had already tried `summary(X)` but that only prints `"GRanges object with 191743 ranges and 7 metadata columns"`.  `dim(X)` gives `NULL`.
terdon
CC BY-SA 4.0
2019-03-06T11:58:23.407Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/90663/discussion-between-terdon-and-llrs).
terdon
<r><bioconductor><cnv>
cn.MOPS fails with 'missing value where TRUE/FALSE needed'
# cn.MOPS fails with 'missing value where TRUE/FALSE needed'

(Asked by: username_0 on 2019-03-06)

I am trying to use cn.MOPS to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:
source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles <- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS <- exomecn.mops(X)

The last command, resCNMOPS <- exomecn.mops(X) fails with:
Error in if (all(segMedianT == 0)) { : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments

If I run traceback immediately afterwards, I get:
> traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth, 
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact, 
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu, 
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor, 
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold, 
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount, 
       useMedian = useMedian, returnPosterior = returnPosterior, 
       ...)
1: exomecn.mops(X)

I am using R version 3.4.4 on an Ubuntu 16.04.5 LTS machine.
My bed file only has one region which is on chrUn_gl000228:
$ grep chrUn_gl000228 region.bed 
chrUn_gl000228  112604  113879

If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?
Presumably, once I have built these two objects:
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)

I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the GRanges object. But how can I do this?

**Tags:** r, bioconductor, cnv
> **Comments:**
> **username_1**: You start your question with `biocLite`, which indicates that this is not the latest version of R and Bioconductor. Which versions are you using? Have you tried with the most up to date version? Why does it complain that you have low number of segments? Could you post how many segments do you have on `X`?  
> **username_0**: @llrs no, I haven't used the latest, but upgrading is not an option (this is a large server). I would be happy to post that, but can you give me a command for it? My R knowledge is extremely limited. That sort of command is exactly the kind of thing I am asking about. I have no idea how to inspect and get statistics about the contents of a GRanges object.  
> **username_0**: Also, hang on, why are you saying it isn't the latest version? I got that code from the documentation of cn.MOPS. Is there a better way to load `biocLite`?  
> **username_1**: Well you will need to use `sessionInfo()` and paste here the result to see the version of R and Bioconductor. The most up to date documentation of cn.MOPS is on the bioconductor page AFAIK (if packages are not up to the standards they are kicked out of the project). Try with `summary(X)` and `dim(X)` or simply print the object and post the first lines  
> **username_1**: I said you weren't using the latest version because in the latest version the recommended way to install packages changed to `BiocManager::install("cn.Mops")` you need a package in CRAN to install packages of Bioconductor  
> **username_0**: @llrs yes, I added the R version (3.4.4). I had already tried `summary(X)` but that only prints `"GRanges object with 191743 ranges and 7 metadata columns"`.  `dim(X)` gives `NULL`.  
> **username_0**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/90663/discussion-between-terdon-and-llrs).  

---

## 1 answers

### Answer by username_0 on 2019-03-06. Score: 2 | Accepted answer

I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome `chrUn_gl000228`. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the `exomecn.mops(X)` was failing.

So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on `chrUn_gl000228` and confirmed that the program worked as expected and didn't throw any errors.


---

