Question: How to use cn.mops to detect CNVs from Targeted DNA Sequencing data?
gravatar for
2.9 years ago by
amalbouzid.ing0 wrote:


I'm using cn.mops to detect CNVs from Targeted DNA Sequencing data by "TruSight One Sequencing Panel from Illumina" investigated for one patient and one control samples, on an Illumina MiSeq sequencer with 150 bp paired-end run.
By running the algorithm with the referencecn.mops function, it results #GRanges object with 67 ranges.
Then, by verifying, the resulting CNVs with the coverage of the enriched regions in the raw data,
I found that several detected CNVs aren't included in the targeted regions enriched by the panel!!
Please, How can you explain that?
Find below my code, what am I doing wrong?


BAMFiles <- list.files(pattern=".bam$")

bamDataRanges<-getReadCountsFromBAM(BAMFiles,sampleNames=paste("Sample",2:3), refSeqName= refSeqName,mode="paired",WL=5000)

case.ctrl2<-referencecn.mops(cases=bamDataRanges[,3],controls=bamDataRanges[,2])  case.ctrl2.ICN <- calcIntegerCopyNumbers(case.ctrl2)
(cnvs(case.ctrl2.ICN ))




ADD COMMENTlink modified 2.9 years ago by Günter Klambauer540 • written 2.9 years ago by amalbouzid.ing0
Answer: How to use cn.mops to detect CNVs from Targeted DNA Sequencing data?
gravatar for Günter Klambauer
2.9 years ago by
Günter Klambauer540 wrote:

Hello Amal,

The reason for this behaviour of cn.mops is that you used the "whole genome sequencing" script. The function "getReadCountsFromBAM" partitions the genome into equally sized windows and counts how many reads are in these windows. However, you have targeted sequencing data and you should therefore restrict the analysis to the targeted segments. For this case, you can use the function "getSegmentReadCountsFromBAM" and "exomecn.mops". In the following code, I assume that you have a list of targeted regions provided as BED file (tab-separated):



BAMFiles <- list.files(pattern=".bam$")

segments <- read.table("targetRegions.bed",sep="\t",

gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))

X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,mode="unpaired")

resCNMOPS <- exomecn.mops(X)

resCNMOPS <- calcIntegerCopyNumbers(resCNMOPS)


I hope this helps!



ADD COMMENTlink modified 2.9 years ago • written 2.9 years ago by Günter Klambauer540

Thank you Günter for your reply, but I do not have a list of targeted regions provided as BED file, how can I retrieve it??


Best Regards


ADD REPLYlink written 2.9 years ago by amalbouzid.ing0

Hello Amal,

There is a workaround: Go back to the initial "getReadCountsFromBAM" and use a small segment size, e.g. "WL=50". Then you get a read count matrix with segments/windows along the whole genome. Remove rows with zero read counts across all samples (these are the DNA segments that were not targeted). Then run "exomecn.mops" on the the read count matrix with the remaining rows.




ADD REPLYlink written 2.9 years ago by Günter Klambauer540

Hi Günter, 

Thanks for your suggestion!

I tried but I think that I make a mistake somewhere?!

BAMFiles <- list.files(pattern=".bam$")

DataRanges.smallWL <- getReadCountsFromBAM(BAMFiles, sampleNames=paste("Sample",1:3), refSeqName= refSeqName,mode="paired",WL=50)

#results: GRanges object with 61913560 ranges and 3 metadata columns:

matrix <- DataRanges.smallWL[DataRanges.smallWL[,1:3] != 0, ]

#results: GRanges object with 0 ranges and 3 metadata columns!



ADD REPLYlink written 2.9 years ago by amalbouzid.ing0

I think your code for filtering the DNA segments does not work. I give you code how to do this on the data that comes with cn.MOPS:



idx <- apply(values(XRanges),1, function(x) !all(x < 70))

XRangesSmall <- XRanges[idx, ]

The code above removes DNA segments where all samples have less than 70 reads. You should use a smaller number, e.g. 5 or 10, instead of 70.



ADD REPLYlink written 2.9 years ago by Günter Klambauer540

Dear  Günter

Thanks a lot for all your suggestions that helped me cope! Just to remember, I'm working with targeted sequencing on two samples, one case and one control. I tried with your two previous proposals: 1/The first one:To generate the read counts from BAM Files comparing to whole genome use a small segment size

>BAM2Files <- list.files(pattern=".bam$")
>DataRanges.smallWL.case_ctrl<-getReadCountsFromBAM(BAM2Files, sampleNames=paste("Sample",1:2), refSeqName= refSeqName,mode="paired",WL=50)
>idx.case_ctrl <- apply(values(DataRanges.smallWL.case_ctrl),1, function(x) !all(x < 5))
> table(idx.case_ctrl)
#61436686   476874
>DataRanges.smallWL.case_ctrl.idx <- DataRanges.smallWL.case_ctrl[idx.case_ctrl, ]
#GRanges object with 176 ranges and 4 metadata columns

## I tried also with referencecn.mops:

>ref.DataRanges.smallWL.case_ctrl.idx<referencecn.mops(cases=DataRanges.smallWL.case_ctrl.idx[,2],controls= DataRanges.smallWL.case_ctrl.idx[,1])
#GRanges object with 45 ranges and 4 metadata columns:

# =>Of a first view, overlapping the two resulted GRanges objects only 5 CNVrs have the CNVs (the same starts)!

2/The second one: To generate the read counts from BAM Files basing on the list of targeted regions in the panel provided as a BED file.

>BAM2Files <- list.files(pattern=".bam$")
>targeted.segments <- read.table("TruSight_One_v1.1.bed",sep="\t",
> dim(targeted.segments) # [1] 62309     4
> exomecn.getSegmentReadCounts <- exomecn.mops(getSegmentReadCounts)
> exomecn.getSegmentReadCounts
#CNV regions:
# GRanges object with 5 ranges and 2 metadata columns:

#I tried also with referencecn.mops insted of exomecn.mops

>ref.getSegmentReadCounts<-referencecn.mops(cases=getSegmentReadCounts[,2],controls= getSegmentReadCounts[,1])
> ref.getSegmentReadCounts
#CNV regions:
#GRanges object with 23 ranges and 1 metadata column:

#=> When overlapping the two resulted GRanges objects (with SegmentReadCounts) only 1 CNV was deduced with the same region! And it’s also the only CNVr between the four previous resulted CNV regions GRanges objects.

I expected more overlapping CNVs?! Really, I lost between the different methods, on which one I have to trust?  And how to choose the best approach for my analysis?

Looking for youy reply,                                                                                                 Many thanks, most grateful!

Best, Amal

ADD REPLYlink written 2.8 years ago by amalbouzid.ing0

Please,  there is someone who can answer me for my previous comment??!

ADD REPLYlink written 2.8 years ago by amalbouzid.ing0

Hello Amal,

Any results that do not come from "referencecn.mops" do not make sense, since you only have two samples. If there is a CNV in one of the two samples, it is randomly assigned to one of the two... The assumption of cn.mops (including exomecn.mops) is that the majority of samples has copy number 2, but here -- if there is a CNV -- the "majority" can only be 50%. I would not recommend to use cn.mops below 8 samples.

To sum it up: you should trust only the results of "referencecn.mops".




ADD REPLYlink written 2.8 years ago by Günter Klambauer540

Dear Günter, 

Many thanks for your reply. 

I will take that into account!



ADD REPLYlink written 2.8 years ago by amalbouzid.ing0
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 394 users visited in the last hour