[Bioc-devel] contradictions in the return value of predictCoding() from VariantAnnotation
Dear Valerie and colleagues, I have recently started using the VariantAnnotation package, which is a huge asset for much of my work with sequence variants - thank you! I have a question regarding the following example (it makes use of two small text files, "annot.gtf" and "variants.vcf", attached to this message): library(VariantAnnotation) library(GenomicFeatures) library(BSgenome.Celegans.UCSC.ce6) # build a TranscriptDb from "annot.gtf" chrominfo <- data.frame(chrom=as.character(seqnames(Celegans)), length=seqlengths(Celegans), is_circular=FALSE) txdb <- makeTranscriptDbFromGFF(file="annot.gtf", format="gtf", chrominfo=chrominfo, dataSource="WormBase v. WS190", species="Caenorhabditis elegans") txdb #TranscriptDb object: #| Db type: TranscriptDb #| Supporting package: GenomicFeatures #| Data source: WormBase v. WS190 #| Organism: Caenorhabditis elegans #| miRBase build ID: NA #| transcript_nrow: 7 #| exon_nrow: 42 #| cds_nrow: 42 #| Db created by: GenomicFeatures package from Bioconductor #| Creation time: 2013-03-07 14:21:08 +0100 (Thu, 07 Mar 2013) #| GenomicFeatures version at creation time: 1.11.14 #| RSQLite version at creation time: 0.11.2 #| DBSCHEMAVERSION: 1.0 # read in sequence variants vcf <- readVcf("variants.vcf", genome="ce6") dim(vcf) #[1] 3 1 # predict the impact on coding sequences aa <- predictCoding(query=vcf, subject=txdb, seqSource=Celegans) length(aa) #[1] 7 # My question is related to the impact of variant "chrV:15822727": # This SNV is overlapping two transcripts (same base, plus strand), # yet the reported VARCODON values are different (GAT -> TAT/AAT): aa[names(aa)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] #GRanges with 2 ranges and 5 metadata columns: #seqnames ranges strand | # | # chrV:15822727 chrV [15822727, 15822727] + | # chrV:15822727 chrV [15822727, 15822727] + | # REFALT varAllele # # chrV:15822727 G A A # chrV:15822727 G A A # REFCODON VARCODON # # chrV:15822727GATTAT # chrV:15822727GATAAT # --- # seqlengths: # chrI chrV # NA NA # interestingly, when altering the annotation or the set of variants, # this contradictions disapears (GAT -> AAT): vcf2 <- vcf[-1] # remove the first SNV aa2 <- predictCoding(query=vcf2, subject=txdb, seqSource=Celegans) length(aa2) #[1] 5 aa2[names(aa2)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] #GRanges with 2 ranges and 5 metadata columns: #seqnames ranges strand | # | # chrV:15822727 chrV [15822727, 15822727] + | # chrV:15822727 chrV [15822727, 15822727] + | # REFALT varAllele # # chrV:15822727 G A A # chrV:15822727 G A A # REFCODON VARCODON # # chrV:15822727GATAAT # chrV:15822727GATAAT # --- # seqlengths: # chrI chrV # NA NA Do you know why this could be the case? Regards, Michael Here is my envirnoment: sessionInfo() R Under development (unstable) (2013-02-25 r62061) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] C attached base packages: [1] parallel stats graphics grDevices utils datasets [7] methods base other attached packages: [1] BSgenome.Celegans.UCSC.ce6_1.3.19 [2] BSgenome_1.27.1 [3] GenomicFeatures_1.11.14 [4] AnnotationDbi_1.21.12 [5] Biobase_2.19.3 [6] VariantAnnotation_1.5.42 [7] Rsamtools_1.11.21 [8] Biostrings_2.27.11 [9] GenomicRanges_1.11.35 [10] IRanges_1.17.35 [11] BiocGenerics_0.5.6 [12] RColorBrewer_1.0-5 loaded via a namespace (and not attached): [1] DBI_0.2-5 RCurl_1.95-4.1 RSQLite_0.11.2 [4] XML_3.95-0.1 biomaRt_2.15.0 bitops_1.0-5 [7] rtracklayer_1.19.9 stats4_3.0.0 tools_3.0.0 [10] zlibbioc_1.5.0 chrVCoding_transcript exon6512065 6512138 . - 2 gene_id "WBGene1073"; transcript_id "F46E10.9.1"; chrVCoding_transcript exon6512189 6512279 . - 0 gene_id "WBGene1073"; transcript_id "F46E10.9.1"; chrVCoding_transcript exon6512325 6512472 . - 1 gene_id "WBGene1073"; transcript_id "F46E10.9.1"; chrVCoding_transcript exon6512521 6512649 . - 1 gene_id "WBGene1073"; transcript_id "F46E10.9.1"; chrVCoding_t
Re: [Bioc-devel] contradictions in the return value of predictCoding() from VariantAnnotation
Here is the second attachment ("variants.vcf") that was missing from the first message. Michael On 03/07/2013 02:40 PM, Michael Stadler wrote: > Dear Valerie and colleagues, > > I have recently started using the VariantAnnotation package, which is a > huge asset for much of my work with sequence variants - thank you! > > I have a question regarding the following example (it makes use of two > small text files, "annot.gtf" and "variants.vcf", attached to this message): > > library(VariantAnnotation) > library(GenomicFeatures) > library(BSgenome.Celegans.UCSC.ce6) > > # build a TranscriptDb from "annot.gtf" > chrominfo <- data.frame(chrom=as.character(seqnames(Celegans)), > length=seqlengths(Celegans), > is_circular=FALSE) > txdb <- makeTranscriptDbFromGFF(file="annot.gtf", > format="gtf", > chrominfo=chrominfo, > dataSource="WormBase v. WS190", > species="Caenorhabditis elegans") > txdb > #TranscriptDb object: > #| Db type: TranscriptDb > #| Supporting package: GenomicFeatures > #| Data source: WormBase v. WS190 > #| Organism: Caenorhabditis elegans > #| miRBase build ID: NA > #| transcript_nrow: 7 > #| exon_nrow: 42 > #| cds_nrow: 42 > #| Db created by: GenomicFeatures package from Bioconductor > #| Creation time: 2013-03-07 14:21:08 +0100 (Thu, 07 Mar 2013) > #| GenomicFeatures version at creation time: 1.11.14 > #| RSQLite version at creation time: 0.11.2 > #| DBSCHEMAVERSION: 1.0 > > # read in sequence variants > vcf <- readVcf("variants.vcf", genome="ce6") > dim(vcf) > #[1] 3 1 > > # predict the impact on coding sequences > aa <- predictCoding(query=vcf, subject=txdb, seqSource=Celegans) > length(aa) > #[1] 7 > > > # My question is related to the impact of variant "chrV:15822727": > # This SNV is overlapping two transcripts (same base, plus strand), > # yet the reported VARCODON values are different (GAT -> TAT/AAT): > aa[names(aa)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] > > #GRanges with 2 ranges and 5 metadata columns: > #seqnames ranges strand | > # | > # chrV:15822727 chrV [15822727, 15822727] + | > # chrV:15822727 chrV [15822727, 15822727] + | > # REFALT varAllele > # > # chrV:15822727 G A A > # chrV:15822727 G A A > # REFCODON VARCODON > # > # chrV:15822727GATTAT > # chrV:15822727GATAAT > # --- > # seqlengths: > # chrI chrV > # NA NA > > > # interestingly, when altering the annotation or the set of variants, > # this contradictions disapears (GAT -> AAT): > vcf2 <- vcf[-1] # remove the first SNV > aa2 <- predictCoding(query=vcf2, subject=txdb, seqSource=Celegans) > length(aa2) > #[1] 5 > aa2[names(aa2)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] > > #GRanges with 2 ranges and 5 metadata columns: > #seqnames ranges strand | > # | > # chrV:15822727 chrV [15822727, 15822727] + | > # chrV:15822727 chrV [15822727, 15822727] + | > # REFALT varAllele > # > # chrV:15822727 G A A > # chrV:15822727 G A A > # REFCODON VARCODON > # > # chrV:15822727GATAAT > # chrV:15822727GATAAT > # --- > # seqlengths: > # chrI chrV > # NA NA > > Do you know why this could be the case? > Regards, > Michael > > > Here is my envirnoment: > sessionInfo() > R Under development (unstable) (2013-02-25 r62061) > Platform: x86_64-unknown-linux-gnu (64-bit) > > locale: > [1] C > > attached base packages: > [1] parallel stats graphics grDevices utils datasets > [7] methods base > > other attached packages: > [1] BSgenome.Celegans.UCSC.ce6_1.3.19 > [2] BSgenome_1.27.1 > [3] GenomicFeatures_1.11.14 > [4] AnnotationDbi_1.21.12 > [5] Biobase_2.19.3 > [6] VariantAnnotation_1.5.42 > [7] Rsamtools_1.11.21 > [8] Biostrings_2.27.11 > [9] GenomicRanges_1.11.35 > [10] IRanges_1.17.35 > [11] BiocGenerics_0.5.6 > [12] RColorBrewer_1.0-5 > > loaded via a namespace (and not attached): > [1] DBI_0.2-5 RCurl_1.95-4.1 RSQLite_0.11.2 > [4] XML_3.95-0.1 biomaRt_2.15.0 bitops_1.0-5 > [7] rtracklayer_1.19.9 stats4_3.0.0 tools_3.0.0 > [10] zlibbioc_1.5.0 > > > > ___ > Bioc-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/bioc-
Re: [Bioc-devel] contradictions in the return value of predictCoding() from VariantAnnotation
Hi Michael, Thanks for reporting the bug and sending a reproducible example. Sorry it took a few days to get to this. A fix has been checked into versions 1.4.12 (release) and 1.5.44 (devel). Valerie On 03/07/13 05:43, Michael Stadler wrote: Here is the second attachment ("variants.vcf") that was missing from the first message. Michael On 03/07/2013 02:40 PM, Michael Stadler wrote: Dear Valerie and colleagues, I have recently started using the VariantAnnotation package, which is a huge asset for much of my work with sequence variants - thank you! I have a question regarding the following example (it makes use of two small text files, "annot.gtf" and "variants.vcf", attached to this message): library(VariantAnnotation) library(GenomicFeatures) library(BSgenome.Celegans.UCSC.ce6) # build a TranscriptDb from "annot.gtf" chrominfo <- data.frame(chrom=as.character(seqnames(Celegans)), length=seqlengths(Celegans), is_circular=FALSE) txdb <- makeTranscriptDbFromGFF(file="annot.gtf", format="gtf", chrominfo=chrominfo, dataSource="WormBase v. WS190", species="Caenorhabditis elegans") txdb #TranscriptDb object: #| Db type: TranscriptDb #| Supporting package: GenomicFeatures #| Data source: WormBase v. WS190 #| Organism: Caenorhabditis elegans #| miRBase build ID: NA #| transcript_nrow: 7 #| exon_nrow: 42 #| cds_nrow: 42 #| Db created by: GenomicFeatures package from Bioconductor #| Creation time: 2013-03-07 14:21:08 +0100 (Thu, 07 Mar 2013) #| GenomicFeatures version at creation time: 1.11.14 #| RSQLite version at creation time: 0.11.2 #| DBSCHEMAVERSION: 1.0 # read in sequence variants vcf <- readVcf("variants.vcf", genome="ce6") dim(vcf) #[1] 3 1 # predict the impact on coding sequences aa <- predictCoding(query=vcf, subject=txdb, seqSource=Celegans) length(aa) #[1] 7 # My question is related to the impact of variant "chrV:15822727": # This SNV is overlapping two transcripts (same base, plus strand), # yet the reported VARCODON values are different (GAT -> TAT/AAT): aa[names(aa)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] #GRanges with 2 ranges and 5 metadata columns: #seqnames ranges strand | # | # chrV:15822727 chrV [15822727, 15822727] + | # chrV:15822727 chrV [15822727, 15822727] + | # REFALT varAllele # # chrV:15822727 G A A # chrV:15822727 G A A # REFCODON VARCODON # # chrV:15822727GATTAT # chrV:15822727GATAAT # --- # seqlengths: # chrI chrV # NA NA # interestingly, when altering the annotation or the set of variants, # this contradictions disapears (GAT -> AAT): vcf2 <- vcf[-1] # remove the first SNV aa2 <- predictCoding(query=vcf2, subject=txdb, seqSource=Celegans) length(aa2) #[1] 5 aa2[names(aa2)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] #GRanges with 2 ranges and 5 metadata columns: #seqnames ranges strand | # | # chrV:15822727 chrV [15822727, 15822727] + | # chrV:15822727 chrV [15822727, 15822727] + | # REFALT varAllele # # chrV:15822727 G A A # chrV:15822727 G A A # REFCODON VARCODON # # chrV:15822727GATAAT # chrV:15822727GATAAT # --- # seqlengths: # chrI chrV # NA NA Do you know why this could be the case? Regards, Michael Here is my envirnoment: sessionInfo() R Under development (unstable) (2013-02-25 r62061) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] C attached base packages: [1] parallel stats graphics grDevices utils datasets [7] methods base other attached packages: [1] BSgenome.Celegans.UCSC.ce6_1.3.19 [2] BSgenome_1.27.1 [3] GenomicFeatures_1.11.14 [4] AnnotationDbi_1.21.12 [5] Biobase_2.19.3 [6] VariantAnnotation_1.5.42 [7] Rsamtools_1.11.21 [8] Biostrings_2.27.11 [9] GenomicRanges_1.11.35 [10] IRanges_1.17.35 [11] BiocGenerics_0.5.6 [12] RColorBrewer_1.0-5 loaded via a namespace (and not attached): [1] DBI_0.2-5 RCurl_1.95-4.1 RSQLite_0.11.2 [4] XML_3.95-0.1 biomaRt_2.15.0 bitops_1.0-5 [7] rtracklayer_1.19.9 stats4_3.0.0 tools_3.0.0 [10] zlibbioc_1.5.0 ___ Bioc-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
Re: [Bioc-devel] contradictions in the return value of predictCoding() from VariantAnnotation
Hi Valerie, Thank you for the quick fix - it works like a charm! Michael On 03/10/2013 02:45 AM, Valerie Obenchain wrote: > Hi Michael, > > Thanks for reporting the bug and sending a reproducible example. Sorry > it took a few days to get to this. > > A fix has been checked into versions 1.4.12 (release) and 1.5.44 (devel). > > Valerie > > > On 03/07/13 05:43, Michael Stadler wrote: >> Here is the second attachment ("variants.vcf") that was missing from the >> first message. >> >> Michael >> >> On 03/07/2013 02:40 PM, Michael Stadler wrote: >>> Dear Valerie and colleagues, >>> >>> I have recently started using the VariantAnnotation package, which is a >>> huge asset for much of my work with sequence variants - thank you! >>> >>> I have a question regarding the following example (it makes use of two >>> small text files, "annot.gtf" and "variants.vcf", attached to this >>> message): >>> >>> library(VariantAnnotation) >>> library(GenomicFeatures) >>> library(BSgenome.Celegans.UCSC.ce6) >>> >>> # build a TranscriptDb from "annot.gtf" >>> chrominfo <- data.frame(chrom=as.character(seqnames(Celegans)), >>> length=seqlengths(Celegans), >>> is_circular=FALSE) >>> txdb <- makeTranscriptDbFromGFF(file="annot.gtf", >>> format="gtf", >>> chrominfo=chrominfo, >>> dataSource="WormBase v. WS190", >>> species="Caenorhabditis elegans") >>> txdb >>> #TranscriptDb object: >>> #| Db type: TranscriptDb >>> #| Supporting package: GenomicFeatures >>> #| Data source: WormBase v. WS190 >>> #| Organism: Caenorhabditis elegans >>> #| miRBase build ID: NA >>> #| transcript_nrow: 7 >>> #| exon_nrow: 42 >>> #| cds_nrow: 42 >>> #| Db created by: GenomicFeatures package from Bioconductor >>> #| Creation time: 2013-03-07 14:21:08 +0100 (Thu, 07 Mar 2013) >>> #| GenomicFeatures version at creation time: 1.11.14 >>> #| RSQLite version at creation time: 0.11.2 >>> #| DBSCHEMAVERSION: 1.0 >>> >>> # read in sequence variants >>> vcf <- readVcf("variants.vcf", genome="ce6") >>> dim(vcf) >>> #[1] 3 1 >>> >>> # predict the impact on coding sequences >>> aa <- predictCoding(query=vcf, subject=txdb, seqSource=Celegans) >>> length(aa) >>> #[1] 7 >>> >>> >>> # My question is related to the impact of variant "chrV:15822727": >>> # This SNV is overlapping two transcripts (same base, plus strand), >>> # yet the reported VARCODON values are different (GAT -> TAT/AAT): >>> aa[names(aa)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] >>> >>> >>> #GRanges with 2 ranges and 5 metadata columns: >>> #seqnames ranges strand | >>> # | >>> # chrV:15822727 chrV [15822727, 15822727] + | >>> # chrV:15822727 chrV [15822727, 15822727] + | >>> # REFALT varAllele >>> # >>> # chrV:15822727 G A A >>> # chrV:15822727 G A A >>> # REFCODON VARCODON >>> # >>> # chrV:15822727GATTAT >>> # chrV:15822727GATAAT >>> # --- >>> # seqlengths: >>> # chrI chrV >>> # NA NA >>> >>> >>> # interestingly, when altering the annotation or the set of variants, >>> # this contradictions disapears (GAT -> AAT): >>> vcf2 <- vcf[-1] # remove the first SNV >>> aa2 <- predictCoding(query=vcf2, subject=txdb, seqSource=Celegans) >>> length(aa2) >>> #[1] 5 >>> aa2[names(aa2)=="chrV:15822727",c("REF","ALT","varAllele","REFCODON","VARCODON")] >>> >>> >>> #GRanges with 2 ranges and 5 metadata columns: >>> #seqnames ranges strand | >>> # | >>> # chrV:15822727 chrV [15822727, 15822727] + | >>> # chrV:15822727 chrV [15822727, 15822727] + | >>> # REFALT varAllele >>> # >>> # chrV:15822727 G A A >>> # chrV:15822727 G A A >>> # REFCODON VARCODON >>> # >>> # chrV:15822727GATAAT >>> # chrV:15822727GATAAT >>> # --- >>> # seqlengths: >>> # chrI chrV >>> # NA NA >>> >>> Do you know why this could be the case? >>> Regards, >>> Michael >>> >>> >>> Here is my envirnoment: >>> sessionInfo() >>> R Under development (unstable) (2013-02-25 r62061) >>> Platform: x86_64-unknown-linux-gnu (64-bit) >>> >>> locale: >>> [1] C >>> >>> attached base packages: >>> [1] parallel stats graphics grDevices utils datasets >>> [7] methods base >>> >>> other attached packages: >>> [1] BSgenome.Celegans.UCSC.ce6_1.3.19 >>> [2] BSgenome_1.27.1 >>> [3]