Question: Union of GenomicRanges problem when one range is empty
0
gravatar for WouterDeCoster
3.4 years ago by
Antwerp, Belgium
WouterDeCoster110 wrote:

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 • 399 views
ADD COMMENTlink modified 3.4 years ago • written 3.4 years ago by WouterDeCoster110

Would you please provide a reproducible example?

ADD REPLYlink written 3.3 years ago by Michael Lawrence11k
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 16.09
Traffic: 177 users visited in the last hour