Question: cn.mops Error in referencecn.mops Genomic Ranges must be unique
0
gravatar for wes3985
2.5 years ago by
wes39850
wes39850 wrote:

Hi All, I have come across an error in cn.mops that is not clear how to solve. I haven't used the software before so it may be a quick fix for anyone who has. The function ' referencecn.mops' is returning the following error:

    Error in referencecn.mops(X[, 1], X[, 2], norm = 0, I = c(0.025, 0.5,  :
    Genomic Ranges must be unique. Check "all(isUnique(input))" and remove identical segments.
    Execution halted

As I understand, a 'genomic range object' is a read count matrix that the algorithm converts a BAM file into [1]. So it seems to be suggesting that the BAM files are identical? I have done a head(X[,1]) and head(X[,2]) just to show the inputs if anyone can spot a problem here:

> head(X[,1])
GRanges object with 6 ranges and 1 metadata column:
      seqnames         ranges strand | g14_picard_dedup.sorted.bam
         <Rle>      <IRanges>  <Rle> |                   <numeric>
  [1]     chr1 [12098, 12258]      * |            392.216485975478
  [2]     chr1 [12553, 12721]      * |            542.181612966102
  [3]     chr1 [13331, 13701]      * |            1445.20239302964
  [4]     chr1 [30334, 30503]      * |            267.168641623296
  [5]     chr1 [35045, 35544]      * |            2001.22694079488
  [6]     chr1 [35618, 35778]      * |            488.655598409448

> head(X[,2])
GRanges object with 6 ranges and 1 metadata column:
      seqnames         ranges strand | fc14_picard_dedup.sorted.bam
         <Rle>      <IRanges>  <Rle> |                    <numeric>
  [1]     chr1 [12098, 12258]      * |             408.547928776133
  [2]     chr1 [12553, 12721]      * |             617.262631520462
  [3]     chr1 [13331, 13701]      * |             1989.45078360552
  [4]     chr1 [30334, 30503]      * |             386.747940481754
  [5]     chr1 [35045, 35544]      * |              2282.1358116319
  [6]     chr1 [35618, 35778]      * |             756.136631025392

The row identifiers on the left seem to indicate a genomic range which seems to be identical in both files, though this is what I would expect as the algorithm is supposed to compare local regions for variations in depth to look for cnvs. The numeric columns are a bit variable but again I would expect this for regions of similar coverage from different samples. Am I missing something?

Here is the full code up to where the error occurs:

library(cn.mops)
BAMFiles <- c("somatic.bam", "tumour.bam")
BaiFiles <- c("somatic.bai", "tumour.bai")
segments <- read.table("truseq_hg38.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X1 <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,mode="paired", BAIFiles=BaiFiles, parallel=12)
X <- normalizeGenome(X1, normType="poisson")
ref_analysis_14 <- referencecn.mops(X[,1], X[,2],
norm=0,
I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64),
classes = paste0("CN", c(0:8, 16, 32, 64, 128)),
segAlgorithm="DNAcopy")
cn.mops next-gen cnv calling • 923 views
ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by wes39850
1
gravatar for wes3985
2.5 years ago by
wes39850
wes39850 wrote:

For anyone else who comes across this error: I think it is caused by the segments input file having some duplicate lines in it which causes errors downstream in the genomic range objects.

It can be solved as follows:

library(cn.mops)
BAMFiles <- c("somatic.bam", "tumour.bam")
BaiFiles <- c("somatic.bai", "tumour.bai")
segments <- read.table("truseq_hg38.bed",sep="\t",as.is=TRUE)
segments <- unique(segments)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X1 <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,mode="paired", BAIFiles=BaiFiles, parallel=12)
X <- normalizeGenome(X1, normType="poisson")
ref_analysis_14 <- referencecn.mops(X[,1], X[,2],
norm=0,
I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64),
classes = paste0("CN", c(0:8, 16, 32, 64, 128)),
segAlgorithm="DNAcopy")
ADD COMMENTlink written 2.5 years ago by wes39850
Please log in to add an answer.

Help
Access

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