[Bioc-devel] Gene annotation: TxDb vs ENSEMBL/NCBI inconsistency
On 06/08/2015 11:43 PM, Rainer Johannes wrote:
dear Robert and Ludwig, the EnsDb packages provide all the gene/transcript etc annotations for all genes defined in the Ensembl database (for a given species and Ensembl release). Except the column/attribute "entrezid" that is stored in the internal database there is however no link to NCBI or UCSC annotations. So, basically, if you want to use "pure" Ensembl based annotations: use EnsDb, if you want to have the UCSC annotations: use the TxDb packages. In case you need EnsDbs of other species or Ensembl versions, the ensembldb package provides functionality to generate such packages either using the Ensembl Perl API or using GTF files provided by Ensembl. If you have problems building the packages, just drop me a line and I'll do that.
Two other sources of Ensembl TxDb's are GenomicFeatures::makeTxDbFromBiomart() and AnnotationHub. For the latter, I'll add a variant of the following to the AnnotationHub HOWTO vignette http://bioconductor.org/packages/devel/bioc/html/AnnotationHub.html later today. ## Gene models _Bioconductor_ represents gene models using 'transcript' databases. These are available via packages such as [TxDb.Hsapiens.UCSC.hg38.knownGene](http://bioconductor.org/packages/TxDb.Hsapiens.UCSC.knownGene.html), or can be constructed using functions such as `[GenomicFeatures](http://bioconductor.org/packages/GenomicFeatures.html)::makeTxDbFromBiomart()` or `GenomicFeatures::makeTxDbFromGRanges()`. _AnnotationHub_ provides an easy way to work with gene models published by Ensembl. Here we discover the Ensemble release 80 r esources for pufferfish,_Takifugu rubripes_ ```{r takifugu-gene-models} query(ah, c("Takifugu", "release-80")) ``` We see that there is a GTF file, as well as various DNA sequences. Let's retrieve the GTF and top-level sequence files. The GTF file is imported as a _GRanges_ instance, the DNA sequence as a compressed, indexed Fasta file ```{r takifugi-data} gtf <- ah[["AH47101"]] dna <- ah[["AH47477"]] head(gtf, 3) dna head(seqlevels(dna)) ``` It is trivial to make a TxDb instance ```{r takifugi-txdb} library(GenomicFeatures) txdb <- makeTxDbFromGRanges(gtf) ```` and to use that in conjunction with the DNA sequence, e.g., to find exon sequences of all annotated genes. ```{r takifugi-exons} library(Rsamtools) # for getSeq,FaFile-method exons <- exons(txdb) getSeq(dna, exons) ``` Some difficulties arise when working with this partly assembled genome that require more advanced GenomicRanges skills, see the [GenomicRanges](http://bioconductor.org/packages/GenomicRanges.html) vignettes, especially "GenomicRanges HOWTOs" and "An Introduction to GenomicRanges".
cheers, jo
On 03 Jun 2015, at 15:56, Robert M. Flight <rflight79 at gmail.com> wrote: Ludwig, If you do this search on the UCSC genome browser (which this annotation package is built from), you will see that the longest variant is what is shown http://genome.ucsc.edu/cgi-bin/hgTracks?clade=mammal&org=Human&db=hg38&position=brca1&hgt.positionInput=brca1&hgt.suggestTrack=knownGene&Submit=submit&hgsid=429339723_8sd4QD2jSAnAsa6cVCevtoOy4GAz&pix=1885
If instead of "genes" you do "transcripts", you will see 20 different
transcripts for this gene, including the one listed by NCBI. I havent tried it yet (haven't upgraded R or bioconductor to latest version), but there is now an Ensembl based annotation package as well, that may work better?? http://bioconductor.org/packages/release/data/annotation/html/EnsDb.Hsapiens.v79.html
-Robert
On Wed, Jun 3, 2015 at 7:04 AM Ludwig Geistlinger < Ludwig.Geistlinger at bio.ifi.lmu.de> wrote:
Dear Bioc annotation team, Querying TxDb.Hsapiens.UCSC.hg38.knownGene for gene coordinates, e.g. for BRCA1; ENSG00000012048; entrez:672 via
genes(TxDb.Hsapiens.UCSC.hg38.knownGene, vals=list(gene_id="672"))
gives me: GRanges object with 1 range and 1 metadata column: seqnames ranges strand | gene_id <Rle> <IRanges> <Rle> | <character> 672 chr17 [43044295, 43170403] - | 672 ------- seqinfo: 455 sequences (1 circular) from hg38 genome However, querying Ensembl and NCBI Gene http://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000012048
http://www.ncbi.nlm.nih.gov/gene/672
the gene is located at (note the difference in the end position) Chromosome 17: 43,044,295-43,125,483 reverse strand How is the inconsistency explained and how to extract an ENSEMBL/NCBI conform annotation from the TxDb object? (I am aware of biomaRt, but I want to explicitely use the Bioc annotation functionality). Thanks! Ludwig -- Dipl.-Bioinf. Ludwig Geistlinger Lehr- und Forschungseinheit f?r Bioinformatik Institut f?r Informatik Ludwig-Maximilians-Universit?t M?nchen Amalienstrasse 17, 2. Stock, B?ro A201 80333 M?nchen Tel.: 089-2180-4067 eMail: Ludwig.Geistlinger at bio.ifi.lmu.de
_______________________________________________ Bioc-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
[[alternative HTML version deleted]]
_______________________________________________ Bioc-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
_______________________________________________ Bioc-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793