Union of GenomicRanges problem when one range is empty
0
0
Entering edit mode
@wouterdecoster-10838
Last seen 2.7 years 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 • 579 views
0
Entering edit mode

Would you please provide a reproducible example?