[Bioc-devel] little tiny bug in CDSID annotations from predictCoding()

[Robert Castelo]

hi,

there is a little tiny bug in the current devel version of 
VariantAnnotation::predictCoding(), and more concretely within
VariantAnnotation:::.localCoordinates(), that precludes the correct 
annotation of the CDSID column:

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

vcf <- readVcf(system.file("extdata", "CEUtrio.vcf.bgz", 
package="VariantFiltering"), genome="hg19")
seqlevelsStyle(vcf) <- seqlevelsStyle(txdb)
vcf <- dropSeqlevels(vcf, "chrM")
coding1 <- predictCoding(vcf, txdb, Hsapiens)
head(coding1$CDSID)
IntegerList of length 6
[[1]] integer(0)
[[2]] integer(0)
[[3]] integer(0)
[[4]] integer(0)
[[5]] integer(0)
[[6]] integer(0)
table(elementLengths(coding1$CDSID))

    0
6038

my sessionInfo() is at the end of the message.

here is the patch, just replacing 'map2$trancriptHits' by 
'map2$transcriptsHits':

--- R/AllUtilities.R	(revision 100756)
+++ R/AllUtilities.R	(working copy)
@@ -284,7 +284,7 @@
      cdsid <- IntegerList(integer(0))
      map2 <- mapToTranscripts(unname(from)[xHits], to,
                               ignore.strand=ignore.strand)
-    cds <- mcols(unlist(to, use.names=FALSE))$cds_id[map2$trancriptsHits]
+    cds <- mcols(unlist(to, use.names=FALSE))$cds_id[map2$transcriptsHits]
      if (length(cds)) {
          cdslst <- unique(splitAsList(cds, map2$xHits))
          cdsid <- cdslst

with this fix then things seem to work again:

coding1 <- predictCoding(vcf, txdb, Hsapiens)
 > head(coding1$CDSID)
IntegerList of length 6
[["1"]] 21771
[["2"]] 21771
[["3"]] 21771
[["4"]] 21771
[["5"]] 21428
[["6"]] 21428
table(elementLengths(coding1$CDSID))

    1    2    3    4    5    6    7    8    9   10   12   13   14   16   19
  873 1229 1024  993  615  524  324  168   82   21   12   15   42   76   40

while investigating this bug i used VariantAnnotation::locateVariants() 
which also annotates the CDSID column, and it seemed to be working. 
however, i noticed that both, predictCoding() and locateVariants(), do 
not give an identical annotation for the CDSID column in coding variants:

coding2 <- locateVariants(vcf, txdb, CodingVariants())
head(loc_all$CDSID)
IntegerList of length 6
[["1"]] 210777
[["2"]] 210777
[["3"]] 210777
[["4"]] 210778
[["5"]] 208140
[["6"]] 208141
table(elementLengths(coding2$CDSID))

    1    2    3    4
4987  901  138   12

in principle, it seems that both are annotating valid CDSID keys:

allcdsinfo <- select(txdb, keys=keys(txdb, keytype="CDSID"), 
columns="CDSID", keytype="CDSID")
sum(!as.character(unlist(coding1$CDSID, use.names=FALSE)) %in% 
allcdsinfo$CDSID)
[1] 0
sum(!as.character(unlist(coding2$CDSID, use.names=FALSE)) %in% 
allcdsinfo$CDSID)
[1] 0

but predictCoding() annotates CDSID values that are not present in 
locateVariants() annotations and viceversa:

sum(!as.character(unlist(coding1$CDSID, use.names=FALSE)) %in% 
as.character(unlist(coding2$CDSID, use.names=FALSE)))
[1] 24057
sum(!as.character(unlist(coding2$CDSID, use.names=FALSE)) %in% 
as.character(unlist(coding1$CDSID, use.names=FALSE)))
[1] 7251

length(unique(intersect(as.character(unlist(coding2$CDSID, 
use.names=FALSE)), as.character(unlist(coding1$CDSID, use.names=FALSE)))))
[1] 0

should not both annotate the same CDSID values on coding variants?


thanks!
robert.
ps: 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
  [3] LC_TIME=en_US.UTF8        LC_COLLATE=en_US.UTF8
  [5] LC_MONETARY=en_US.UTF8    LC_MESSAGES=en_US.UTF8
  [7] LC_PAPER=en_US.UTF8       LC_NAME=C
  [9] LC_ADDRESS=C              LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF8 LC_IDENTIFICATION=C

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

other attached packages:
  [1] TxDb.Hsapiens.UCSC.hg19.knownGene_3.1.1
  [2] GenomicFeatures_1.19.31
  [3] AnnotationDbi_1.29.17
  [4] Biobase_2.27.2
  [5] BSgenome.Hsapiens.UCSC.hg19_1.4.0
  [6] BSgenome_1.35.17
  [7] rtracklayer_1.27.8
  [8] VariantAnnotation_1.13.40
  [9] Rsamtools_1.19.44
[10] Biostrings_2.35.11
[11] XVector_0.7.4
[12] GenomicRanges_1.19.46
[13] GenomeInfoDb_1.3.13
[14] IRanges_2.1.43
[15] S4Vectors_0.5.22
[16] BiocGenerics_0.13.6
[17] vimcom_1.0-0
[18] setwidth_1.0-3
[19] colorout_1.0-3

loaded via a namespace (and not attached):
  [1] BiocParallel_1.1.15      biomaRt_2.23.5
  [3] bitops_1.0-6             DBI_0.3.1
  [5] GenomicAlignments_1.3.31 RCurl_1.95-4.5
  [7] RSQLite_1.0.0            tools_3.2.0
  [9] XML_3.98-1.1             zlibbioc_1.13.2