Dear Günter,

Rencently, I am trying to analyze plant genome CNVs across multiple samples, and notice that your cn.MOPS is very suitable for this. We focused on a de novo plant genome, which has many scaffolds (with no chromosome-level sequences). Unfortunately we got some error message when running it.

What we did is as follows:

The first part of the code is loading library:

options(width=75)
set.seed(0)
library(Biobase)
library(GenomicRanges)
library(GenomeInfoDb)
library(cn.mops)

The second part of the code is loading data, which seems to work fine with no error report:

BAMFiles <- list.files(pattern=".bam$")
bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName="scaffolds_0",WL=2000)
bamdata <- as.data.frame(bamDataRanges)

The third part is caculating the different depth of CNV regions. This works fine for some scaffolds: 

res <- cn.mops(bamDataRanges,I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4),classes = c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6", "CN7", "CN8"))
data(cn.mops)
ls()
result <- calcIntegerCopyNumbers(res)
segm <- as.data.frame(segmentation(result))
CNVs <- as.data.frame(cnvs(result))
CNVRegions <- as.data.frame(cnvr(result))
write.table(segm,file="segmentation.tbl",sep="\t")
write.table(CNVs,file="cnvs.tbl",sep="\t")
write.table(CNVRegions,file="cnvr.tbl",sep="\t")
write.table(bamdata,file="bamdata.tbl",sep="\t")

But, when I change "refSeqName" and call cn.mops function:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName="scaffolds_133",WL=2000)
bamdata <- as.data.frame(bamDataRanges)

res <- cn.mops(bamDataRanges,I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4),classes = c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6", "CN7", "CN8"))

it comes up with the following error message:

Normalizing...
Starting local modeling, please be patient...
Reference sequence:  scaffolds_133
Starting segmentation algorithm...
Using "fastseg" for segmentation.
No CNVs detected. Try changing "normalization", "priorImpact" or "thresholds".

I have tried to change the parameter of  "norm", "priorImpact", "segAlgorithm", "lowerThreshold" or "upperThreshold". Sometime it works, and sometime it doesn't.

Could you please give us some suggestion for this? Any comments will be very much appreciated.

Best Regards

Ponyo

Hello Ponyo,

Thanks for your question. This is actually not an error message, this just means that there were no CNVs detected. The fact that even after changing the hyperparameters, such as "priorImpact", there are not detections means that cn.MOPS provided robust results. In other words, there should be no CNVs on "scaffolds_133".

A few suggestions:

- you can run cn.MOPS on all scaffolds at once using the following command (and assuming you have scaffolds numbered from 1 to 150):

bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName=paste0("scaffolds_",1:150), WL=2000)

- you can use "DNAcopy" as segmentation algorithm:

res <- cn.mops(bamDataRanges,I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4),classes = c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6", "CN7", "CN8"),segAlgorithm="DNAcopy")

- you can try to decrease the window length which would make the method more sensitive

bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName="scaffolds_0",WL=250)

I hope this helps!

Regards,

Günter

Hello Günter,

When I run cn.MOPS on all scaffolds at once using the following command:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName=paste0("scaffolds_",501:1000),WL=2000)

I found new error message:

Identified the following reference sequences:  scaffolds_0,scaffolds_1,scaffolds_2,scaffolds_3,scaffolds_4,scaffolds_5, ......
Using scaffolds_501, scaffolds_502, scaffolds_503, scaffolds_504, scaffolds_505, scaffolds_506, ............
Using indexed BAM files.
Reading file: H1.sort.dedup.realn.bam
  scaffolds_501

  scaffolds_1000
Error in .Call2("solve_user_SEW0", start, end, width, PACKAGE = "IRanges") :
  solving row 1: range cannot be determined from the supplied arguments (too many NAs)
Calls: getReadCountsFromBAM ... newGRanges -> <Anonymous> -> solveUserSEW0 -> .Call2 -> .Call
In addition: There were 50 or more warnings (use warnings() to see the first 50)
Execution halted

My BAM files do have alignments in them, so I'm not sure what's going on or how to debug this problem. I'd appreciate insights!

Best Regards

Ponyo

I find this new error message only with refSeqName=paste0("scaffolds_",751:760):

command:

bamDataRanges <- getReadCountsFromBAM(BAMFiles,mode="paired",refSeqName=paste0("scaffolds_",751:760),WL=2000)

error:

Error in .Call2("solve_user_SEW0", start, end, width, PACKAGE = "IRanges") :
  solving row 1: range cannot be determined from the supplied arguments (too many NAs)
In addition: Warning messages:
1: In FUN(c("H1.sort.dedup.realn.bam", "H10.sort.dedup.realn.bam",  :
  Some read positions are greater than length of reference sequence! File:  H1.sort.dedup.realn.bam

2: In FUN(c("H1.sort.dedup.realn.bam", "H10.sort.dedup.realn.bam",  :
  Some read positions are greater than length of reference sequence! File:  H1.sort.dedup.realn.bam

My BAM files are from bwa && samtools && GATK, and I am sure no read positions are greater than length of reference sequence in File:  H1.sort.dedup.realn.bam.

Hello,

Can you please run "traceback()" right after the error happens and post the result. I suspect that this is not a cn.MOPS error...

Regards,

Günter