[Bioc-devel] subtle error with VariantAnnotation::locateVariants()

Hi Robert,

Thanks for reporting this. Now fixed in 1.13.27.

I introduced this bug when switching the coordinate mapping over from
mapCoords() to the new mapToTranscripts() family. The output from the
new mapper is cleaner and allowed me to tidy some old code. In the
process I realized there could be cases with >1 CDSID per row and
changed the data type from 'integer' to 'IntegerList'. I thought I
changed the default in all places but obviously missed a couple. The
specific bug you hit was complaining that GRanges with different types
for CDSID couldn't not be combined.

Valerie

On 02/03/2015 05:57 AM, Robert Castelo wrote:

hi,

VariantAnnotation::locateVariants() is breaking in a very specific way
in its current devel version (1.13.26). Since I'm using it while working
on VariantFiltering, I have reverted my copy of VariantAnnotation to
version 1.13.24 to be able to continue working.

so at the moment this is not much of a problem for me, but I thought you
might be interested in knowing about it, just in case you were not aware
of it. This is the minimal example I've come up reproducing it:

library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

gr <- GRanges(rep("chr20", 1), IRanges(start=44501458, width=1))

loc <- locateVariants(gr, txdb,
AllVariants(intergenic=IntergenicVariants(0, 0)))
Error in .Primitive("c")(<S4 object of class "CompressedIntegerList">, :
all arguments in '...' must be CompressedList objects
traceback()
traceback()
23: stop("all arguments in '...' must be CompressedList objects")
22: .Primitive("c")(<S4 object of class "CompressedIntegerList">,
NA_integer_)
21: .Primitive("c")(<S4 object of class "CompressedIntegerList">,
NA_integer_)
20: do.call(c, unname(cols))
19: do.call(c, unname(cols))
18: FUN(c("LOCATION", "LOCSTART", "LOCEND", "QUERYID", "TXID", "CDSID",
"GENEID", "PRECEDEID", "FOLLOWID")[[6L]], ...)
17: lapply(colnames(df), function(cn) {
cols <- lapply(args, `[[`, cn)
isRle <- vapply(cols, is, logical(1L), "Rle")
if (any(isRle) && !all(isRle)) {
cols[isRle] <- lapply(cols[isRle], S4Vectors:::decodeRle)
}
isFactor <- vapply(cols, is.factor, logical(1L))
if (any(isFactor)) {
cols <- lapply(cols, as.factor)
levs <- unique(unlist(lapply(cols, levels), use.names =
FALSE))
cols <- lapply(cols, factor, levs)
}
rectangular <- length(dim(cols[[1]])) == 2L
if (rectangular) {
combined <- do.call(rbind, unname(cols))
}
else {
combined <- do.call(c, unname(cols))
}
if (any(isFactor))
combined <- structure(combined, class = "factor", levels =
levs)
combined
})
16: lapply(colnames(df), function(cn) {
cols <- lapply(args, `[[`, cn)
isRle <- vapply(cols, is, logical(1L), "Rle")
if (any(isRle) && !all(isRle)) {
cols[isRle] <- lapply(cols[isRle], S4Vectors:::decodeRle)
}
isFactor <- vapply(cols, is.factor, logical(1L))
if (any(isFactor)) {
cols <- lapply(cols, as.factor)
levs <- unique(unlist(lapply(cols, levels), use.names =
FALSE))
cols <- lapply(cols, factor, levs)
}
rectangular <- length(dim(cols[[1]])) == 2L
if (rectangular) {
combined <- do.call(rbind, unname(cols))
}
else {
combined <- do.call(c, unname(cols))
}
if (any(isFactor))
combined <- structure(combined, class = "factor", levels =
levs)
combined
})
15: .Method(..., deparse.level = deparse.level)
14: eval(expr, envir, enclos)
13: eval(.dotsCall, env)
12: eval(.dotsCall, env)
11: standardGeneric("rbind")
10: (function (..., deparse.level = 1)
standardGeneric("rbind"))(<S4 object of class "DataFrame">, <S4
object of class "DataFrame">,
<S4 object of class "DataFrame">, <S4 object of class
"DataFrame">,
<S4 object of class "DataFrame">, <S4 object of class
"DataFrame">,
<S4 object of class "DataFrame">)
9: do.call(rbind, lapply(x, mcols, FALSE))
8: do.call(rbind, lapply(x, mcols, FALSE))
7: .unlist_list_of_GenomicRanges(args, ignore.mcols = ignore.mcols)
6: .local(x, ..., recursive = recursive)
5: c(coding, intron, fiveUTR, threeUTR, splice, promoter, intergenic)
4: c(coding, intron, fiveUTR, threeUTR, splice, promoter, intergenic)
3: .local(query, subject, region, ...)
2: locateVariants(gr, txdb, AllVariants(intergenic =
IntergenicVariants(0,
0)))
1: locateVariants(gr, txdb, AllVariants(intergenic =
IntergenicVariants(0,
0)))

interestingly, if i replace the chromosome name in the GRanges object of
'chr20' by 'chr2', then it works:

gr <- GRanges(rep("chr2", 1), IRanges(start=44501458, width=1))
loc <- locateVariants(gr, txdb,
AllVariants(intergenic=IntergenicVariants(0, 0)))

or if i replace the start of the ranges of 44501458 by 1, then it also
works:

gr <- GRanges(rep("chr20", 1), IRanges(start=1, width=1))
loc <- locateVariants(gr, txdb,
AllVariants(intergenic=IntergenicVariants(0, 0)))

here is the sessionInfo():

R Under development (unstable) (2014-10-14 r66765)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF8 LC_NUMERIC=C LC_TIME=en_US.UTF8
LC_COLLATE=en_US.UTF8
[5] LC_MONETARY=en_US.UTF8 LC_MESSAGES=en_US.UTF8
LC_PAPER=en_US.UTF8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF8
LC_IDENTIFICATION=C

attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
methods base

other attached packages:
[1] TxDb.Hsapiens.UCSC.hg19.knownGene_3.0.0 GenomicFeatures_1.19.17
AnnotationDbi_1.29.17
[4] Biobase_2.27.1 VariantAnnotation_1.13.26
Rsamtools_1.19.26
[7] Biostrings_2.35.7 XVector_0.7.3
GenomicRanges_1.19.35
[10] GenomeInfoDb_1.3.12 IRanges_2.1.36
S4Vectors_0.5.17
[13] BiocGenerics_0.13.4 vimcom_1.0-0
setwidth_1.0-3
[16] colorout_1.0-3

loaded via a namespace (and not attached):
[1] base64enc_0.1-2 BatchJobs_1.5 BBmisc_1.8
BiocParallel_1.1.13 biomaRt_2.23.5
[6] bitops_1.0-6 brew_1.0-6 BSgenome_1.35.16
checkmate_1.5.1 codetools_0.2-10
[11] DBI_0.3.1 digest_0.6.8 fail_1.2
foreach_1.4.2 GenomicAlignments_1.3.27
[16] iterators_1.0.7 RCurl_1.95-4.5 RSQLite_1.0.0
rtracklayer_1.27.7 sendmailR_1.2-1
[21] stringr_0.6.2 tools_3.2.0 XML_3.98-1.1
zlibbioc_1.13.0


this problem dissapear if i revert to the older version (1.13.24),

svn merge -r 98874:98522 .

build and install it. in a new fresh R-devel session the same code
smoothly with some warnings, i believe unrelated to this problem:

[...]
loc <- locateVariants(gr, txdb,
+ AllVariants(intergenic=IntergenicVariants(0, 0)))
Warning messages:
1: '.local' is deprecated.
Use 'mapToTranscripts' instead.
See help("Deprecated")
2: '.local' is deprecated.
Use 'mapToTranscripts' instead.
See help("Deprecated")
3: '.local' is deprecated.
Use 'mapToTranscripts' instead.
See help("Deprecated")
4: '.local' is deprecated.
Use 'mapToTranscripts' instead.
See help("Deprecated")

this is the sessionInfo() of this run:

R Under development (unstable) (2014-10-14 r66765)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF8 LC_NUMERIC=C LC_TIME=en_US.UTF8
LC_COLLATE=en_US.UTF8
[5] LC_MONETARY=en_US.UTF8 LC_MESSAGES=en_US.UTF8
LC_PAPER=en_US.UTF8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF8
LC_IDENTIFICATION=C

attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
methods base

other attached packages:
[1] TxDb.Hsapiens.UCSC.hg19.knownGene_3.0.0 GenomicFeatures_1.19.17
AnnotationDbi_1.29.17
[4] Biobase_2.27.1 VariantAnnotation_1.13.24
Rsamtools_1.19.26
[7] Biostrings_2.35.7 XVector_0.7.3
GenomicRanges_1.19.35
[10] GenomeInfoDb_1.3.12 IRanges_2.1.36
S4Vectors_0.5.17
[13] BiocGenerics_0.13.4 vimcom_1.0-0
setwidth_1.0-3
[16] colorout_1.0-3

loaded via a namespace (and not attached):
[1] base64enc_0.1-2 BatchJobs_1.5 BBmisc_1.8
BiocParallel_1.1.13 biomaRt_2.23.5
[6] bitops_1.0-6 brew_1.0-6 BSgenome_1.35.16
checkmate_1.5.1 codetools_0.2-10
[11] DBI_0.3.1 digest_0.6.8 fail_1.2
foreach_1.4.2 GenomicAlignments_1.3.27
[16] iterators_1.0.7 RCurl_1.95-4.5 RSQLite_1.0.0
rtracklayer_1.27.7 sendmailR_1.2-1
[21] stringr_0.6.2 tools_3.2.0 XML_3.98-1.1
zlibbioc_1.13.0

cheers,

robert.