cn.mops for WES, what am I doing wrong?

Question

[R.Blues]

Hello everyone,

I am analysing several samples through Whole Exome Sequencing with paired-reads, without having any control sample.

In any case, I am trying to find CNVs through the next code (I only have a BAM file in the Working Directory, "T.bam"). I am not sure if I should use singlecn.mops or exomecn.mops, as I do not have any control sample:

library(cn.mops)
library(Repitools)
BAMFiles <- list.files(pattern=".bam$")

BGR = BAM2GRanges(BAMFiles)
bamDataRanges <- getSegmentReadCountsFromBAM(BAMFiles, GR= BGR, mode = "paired", BAIFiles = "T.bai")

ras <- exomecn.mops(bamDataRanges)
plot(ras, which=1)
segplot(ras)

I get the following error, after a HUGE time of analysis in getSegmentReadCountsFromBAM(which I suppose is normal, as I am running it locally and not on a server, but the BAM file is only 2-3 Gb):

Error in value[[3L]](cond) : 'countBam' failed:
  record: -1168125164
  error: 0
  file: T.bam
  index: T

Any ideas or suggestions? I get a similar error if I use the parallel option. Maybe it is due to the use of BAM2GRanges, but, to be honest, I do not know how to make then the GRanges object, or how to decide its size :/

Thank you very much for your help.

[Devon Ryan]

This error ends up coming from Rsamtools, actually. I don't know what the underlying cause is, but you might be able to subset the BAM file and run the function on that to determine exactly what read is causing the error. My guess would be that either the index or part of the BAM file is corrupt.

[R.Blues]

Thank you very much. I will try with another .BAM and I will let you know

In this case, even without having a control sample, I should use exomecn, right? Do you know if it will work properly, or if I should search for another alternative for finding CNV in tumoral exome data without control sample? In that case, do you know any?

Thanks a lot again, and sorry for all these questions, but I am new to this world! :)

[Devon Ryan]

I won't pretend to be particularly knowledgeable with cn.mops, so hopefully someone else will know the exact answer to that and be able to better help there. I would say that finding CNVs without a control sample is a rather difficult thing to do with whole exome sequencing, so I'd be rather hesitant about believing the results from any package.

Answer 1

[brentp]

This may not solve your issue, but since you're using exome data, you should give it the capture file as BED and then your code would look like this:

segments = read.delim(capture_regions, header=FALSE) 
gr = GRanges(segments[,1], IRanges(segments[,2], segments[,3])) 
bam_counts = getSegmentReadCountsFromBAM(bams, GR=gr, mode="paired", parallel=4) 
res = exomecn.mops(bam_counts, parallel=4, priorImpact=20)

[R.Blues]

Thank you VERY much!! It solved the issue, moreover!