Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-06T17:01:53.317Z
Answers.LastActivityDate: 2019-05-31T08:35:59.723Z
Answers.LastEditDate: 2019-05-31T08:35:59.723Z
Answers.UserDisplayName: terdon
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:42:01.447Z
Comments.Text: 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`?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:43:49.293Z
Comments.Text: @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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:44:48.520Z
Comments.Text: 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`?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:54:12.417Z
Comments.Text: 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
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:55:58.023Z
Comments.Text: 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
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:57:27.323Z
Comments.Text: @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`.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-06T11:58:23.407Z
Comments.Text: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/90663/discussion-between-terdon-and-llrs).
Comments.UserDisplayName: terdon
Tags: <r><bioconductor><cnv>
Title: cn.MOPS fails with 'missing value where TRUE/FALSE needed'
ThreadText: # 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](https://www.bioinf.jku.at/software/cnmops/cnmops.html) 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**: @username_1 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**: @username_1 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-username_0-and-username_1).  
  
  ---
  
  ## 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.
  
  
  ---
  
  