Union of GenomicRanges problem when one range is empty
Entering edit mode
Last seen 2.7 years ago
Antwerp, Belgium


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.

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 • 579 views
Entering edit mode

Would you please provide a reproducible example?


Login before adding your answer.

Traffic: 356 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6