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

Thanks! Do you have an explanation for the apparent disagreement in
CDSID annotations that i described below the bug, between
predictCoding() and locateVariants()?

robert.


-------- Original message --------
From: Valerie Obenchain
Date:17/03/2015 19:18 (GMT+01:00)
To: bioc-devel at r-project.org
Subject: Re: [Bioc-devel] little tiny bug in CDSID annotations from
predictCoding()

Hi Robert,

Thanks for reporting the typo and bug. Now fixed in 1.13.41.

Valerie

On 03/17/2015 10:58 AM, Robert Castelo wrote:
? > in my message below, the line that it says:
? >
? > head(loc_all$CDSID)
? >
? > it should say
? >
? > head(coding2$CDSID)
? >
? > cheers,
? >
? > robert.
? >
? > ==========================================
? > 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
? >
? > _______________________________________________
? > Bioc-devel at r-project.org mailing list
? > https://stat.ethz.ch/mailman/listinfo/bioc-devel