text: user298: 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? 
  
  user298: 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.
  
original_text: 298: <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>
   298: <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>
  