Hi,
I'm trying to extract all snps on TNF transcript from SNPlocs.Hsapiens.dbSNP144.GRCh38, but only one is returned (rs281865419). In the NCBI variation viewer, there are 255 snps.
Is my code right? And is there a better approach?
My code and sessionInfo are below.
Many thanks,
Yusuke
----R-code----
library(SNPlocs.Hsapiens.dbSNP144.GRCh38); library(TxDb.Hsapiens.UCSC.hg38.knownGene); library(VariantAnnotation); txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene; tx.by.gene <- transcriptsBy(txdb, "gene"); snps <- SNPlocs.Hsapiens.dbSNP144.GRCh38; chr6 <- snpsBySeqname(snps, "ch6"); chr6 <- renameSeqlevels(chr6, sub("^chrMT", "chrM", sub("^ch", "chr", seqlevels(chr6)))); seqinfo(chr6)@genome <- rep("hg38", 25); hit <- findOverlaps(chr6, tx.by.gene["7124"]); # EntrezID of TNF chr6[hit@queryHits]; GRanges object with 1 range and 2 metadata columns: seqnames ranges strand | RefSNP_id alleles_as_ambig <Rle> <IRanges> <Rle> | <character> <character> [1] chr6 [31577157, 31577157] + | rs281865419 Y ------- seqinfo: 25 sequences (1 circular) from hg38 genome sessionInfo() R version 3.2.0 (2015-04-16) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets [8] methods base other attached packages: [1] VariantAnnotation_1.14.1 [2] Rsamtools_1.20.4 [3] TxDb.Hsapiens.UCSC.hg38.knownGene_3.1.2 [4] GenomicFeatures_1.20.6 [5] AnnotationDbi_1.30.1 [6] Biobase_2.28.0 [7] SNPlocs.Hsapiens.dbSNP144.GRCh38_0.99.11 [8] BSgenome_1.36.3 [9] rtracklayer_1.28.4 [10] Biostrings_2.36.4 [11] XVector_0.8.0 [12] GenomicRanges_1.20.8 [13] GenomeInfoDb_1.4.3 [14] IRanges_2.2.9 [15] S4Vectors_0.6.6 [16] BiocGenerics_0.14.0
Hi Herve,
Thanks for your response.
It would be very helpful and I could find the SNPs in ds_flat_chMulti.flat.gz file.
I hope you will update SNPlocs packages.
Thank you for your help,
Yusuke
I'll take a shot at this. That won't happen before January though...
H.