Union of GenomicRanges problem when one range is empty
0
0
Entering edit mode
@wouterdecoster-10838
Last seen 5 months ago
Antwerp, Belgium

Hello,

I'm integrating the calls of different CNV prediction tools for whole exome sequencing (conifer, exomeCopy and CODEX).
In order to combine their predictions per sample I produce genomicranges objects (split by duplications and deletions, which I obviously don't want to merge), perform pairwise intersection and then get the union of those objects. As such I select regions which are at least supported by two tools, as a form of quality control without being too stringent.

I added my main code below (the part before it in which I format the output of all tools is not included). This works fine most of the time, but occasionally an intersection is empty (0 ranges). When creating the union of an empty Granges object with another I get the following error: `Error in validObject(.Object) :  invalid class "GRanges" object: 'seqlevels(seqinfo(x))' and 'levels(seqnames(x))' are not identical`

However, all intersection objects are (gr12, gr13 and gr23) are valid (tested by e.g. `validObject(gr12)`

Do you have suggestions how to fix this? Thanks in advance.

library(GenomicRanges)
res = data.frame()

for (indiv in unique(conifer$sample)) {
    for (copy in c('dup', 'del')) {
        t1 = subset(conifer, conifer$sample == indiv & conifer$state == copy)
        t2 = subset(exomecopy, exomecopy$sample == indiv & exomecopy$state == copy)
        t3 = subset(codex , codex $sample == indiv & codex $state == copy)
        if (dim(t1)[[1]] > 0 & dim(t2)[[1]] > 0  & dim(t3)[[1]] > 0 ) {
            gr1 = with(t1, GRanges(chromosome, IRanges(start=begin, end=end)))
            gr2 = with(t2, GRanges(chromosome, IRanges(start=begin, end=end)))
            gr3 = with(t3, GRanges(chromosome, IRanges(start=begin, end=end)))
            gr12 = intersect(gr1, gr2)
            gr13 = intersect(gr1, gr3)
            gr23 = intersect(gr2, gr3)
            consensus = as.data.frame(union(union(gr12, gr13), gr23))
                out = cbind(consensus, state = copy, sample = indiv)
                res = rbind(res, out)
                }
            }
    }

 

genomicranges R cnv • 1.1k views
ADD COMMENT
0
Entering edit mode

Would you please provide a reproducible example?

ADD REPLY

Login before adding your answer.

Traffic: 617 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6