1
$\begingroup$

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?

$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$ – llrs Mar 6 '19 at 11:42
  • $\begingroup$ @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. $\endgroup$ – terdon Mar 6 '19 at 11:43
  • $\begingroup$ 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? $\endgroup$ – terdon Mar 6 '19 at 11:44
  • $\begingroup$ 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 $\endgroup$ – llrs Mar 6 '19 at 11:54
  • $\begingroup$ 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 $\endgroup$ – llrs Mar 6 '19 at 11:55
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.