Last updated: 2024-10-23

Checks: 7 0

Knit directory: muse/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20200712) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version b3485f3. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/pbmc3k.csv
    Ignored:    data/pbmc3k.csv.gz
    Ignored:    data/pbmc3k/
    Ignored:    r_packages_4.4.0/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/annotation_dbi.Rmd) and HTML (docs/annotation_dbi.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd b3485f3 Dave Tang 2024-10-23 reactome.db
html de5bf16 Dave Tang 2024-10-23 Build site.
Rmd 2506265 Dave Tang 2024-10-23 Using AnnotationDbi

The Bioconductor annotation packages are an extensive collection of annotations. For this post I simply illustrate the basics of probing these annotation packages. For the first example I will use the org.Hs.eg.db package, which provides genome wide annotations for the human genome.

Installation

To begin, install the {org.Hs.eg.db} package.

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("org.Hs.eg.db")

Package

Load package.

packageVersion("org.Hs.eg.db")
[1] '3.19.1'
suppressPackageStartupMessages(library(org.Hs.eg.db))
class(org.Hs.eg.db)
[1] "OrgDb"
attr(,"package")
[1] "AnnotationDbi"

Getting started

We can query the package by using the AnnotationDbi::select() function; to find out what we can select and return we can use the keys(), columns(), and keytypes() functions.

head(keys(org.Hs.eg.db))
[1] "1"  "2"  "3"  "9"  "10" "11"
#store the first six keys
my_keys <- head(keys(org.Hs.eg.db))
 
keytypes(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
[21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
[26] "UNIPROT"     
columns(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GENETYPE"     "GO"           "GOALL"        "IPI"          "MAP"         
[16] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"        
[21] "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
[26] "UNIPROT"     
AnnotationDbi::select(
  org.Hs.eg.db,
  keys = my_keys,
  columns=c("ENTREZID","SYMBOL","GENENAME"),
  keytype="ENTREZID"
)
'select()' returned 1:1 mapping between keys and columns
  ENTREZID SYMBOL                           GENENAME
1        1   A1BG             alpha-1-B glycoprotein
2        2    A2M              alpha-2-macroglobulin
3        3  A2MP1 alpha-2-macroglobulin pseudogene 1
4        9   NAT1              N-acetyltransferase 1
5       10   NAT2              N-acetyltransferase 2
6       11   NATP     N-acetyltransferase pseudogene

Within the annotation package are many objects, which we can also probe.

ls("package:org.Hs.eg.db")
 [1] "org.Hs.eg"                "org.Hs.eg_dbconn"        
 [3] "org.Hs.eg_dbfile"         "org.Hs.eg_dbInfo"        
 [5] "org.Hs.eg_dbschema"       "org.Hs.eg.db"            
 [7] "org.Hs.egACCNUM"          "org.Hs.egACCNUM2EG"      
 [9] "org.Hs.egALIAS2EG"        "org.Hs.egCHR"            
[11] "org.Hs.egCHRLENGTHS"      "org.Hs.egCHRLOC"         
[13] "org.Hs.egCHRLOCEND"       "org.Hs.egENSEMBL"        
[15] "org.Hs.egENSEMBL2EG"      "org.Hs.egENSEMBLPROT"    
[17] "org.Hs.egENSEMBLPROT2EG"  "org.Hs.egENSEMBLTRANS"   
[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"         
[21] "org.Hs.egENZYME2EG"       "org.Hs.egGENENAME"       
[23] "org.Hs.egGENETYPE"        "org.Hs.egGO"             
[25] "org.Hs.egGO2ALLEGS"       "org.Hs.egGO2EG"          
[27] "org.Hs.egMAP"             "org.Hs.egMAP2EG"         
[29] "org.Hs.egMAPCOUNTS"       "org.Hs.egOMIM"           
[31] "org.Hs.egOMIM2EG"         "org.Hs.egORGANISM"       
[33] "org.Hs.egPATH"            "org.Hs.egPATH2EG"        
[35] "org.Hs.egPFAM"            "org.Hs.egPMID"           
[37] "org.Hs.egPMID2EG"         "org.Hs.egPROSITE"        
[39] "org.Hs.egREFSEQ"          "org.Hs.egREFSEQ2EG"      
[41] "org.Hs.egSYMBOL"          "org.Hs.egSYMBOL2EG"      
[43] "org.Hs.egUCSCKG"          "org.Hs.egUNIPROT"        

Ensembl Gene IDs to Entrez Gene IDs.

class(org.Hs.egENSEMBL2EG)
[1] "AnnDbBimap"
attr(,"package")
[1] "AnnotationDbi"
head(keys(org.Hs.egENSEMBL2EG))
[1] "ENSG00000121410" "ENSG00000175899" "ENSG00000291190" "ENSG00000171428"
[5] "ENSG00000156006" "ENSG00000196136"

Create a lookup table.

my_ensg_keys <- head(keys(org.Hs.egENSEMBL2EG))
 
AnnotationDbi::select(
  org.Hs.eg.db,
  keys = my_ensg_keys,
  columns=c("ENSEMBL","ENTREZID","SYMBOL","GENENAME"),
  keytype="ENSEMBL"
)
'select()' returned 1:1 mapping between keys and columns
          ENSEMBL ENTREZID   SYMBOL                           GENENAME
1 ENSG00000121410        1     A1BG             alpha-1-B glycoprotein
2 ENSG00000175899        2      A2M              alpha-2-macroglobulin
3 ENSG00000291190        3    A2MP1 alpha-2-macroglobulin pseudogene 1
4 ENSG00000171428        9     NAT1              N-acetyltransferase 1
5 ENSG00000156006       10     NAT2              N-acetyltransferase 2
6 ENSG00000196136       12 SERPINA3           serpin family A member 3

Using the annotation package {org.Hs.eg.db} we can easily convert different gene identifiers, obtain their gene symbols, and descriptions (as well as all other keytypes).

reactome.db

reactome.db is a set of annotation maps for Reactome, assembled using data from Reactome. This package has been created by a third-party developer, and is not affiliated with the Reactome team.

To begin, install the {reactome.db} package.

BiocManager::install("reactome.db")

Load the package.

packageVersion('reactome.db')
[1] '1.88.0'
suppressPackageStartupMessages(library(reactome.db))
class(reactome.db)
[1] "ReactomeDb"
attr(,"package")
[1] "AnnotationDbi"

Demo.

keytypes(reactome.db)
[1] "ENTREZID"   "GO"         "PATHID"     "PATHNAME"   "REACTOMEID"
columns(reactome.db)
[1] "ENTREZID"   "GO"         "PATHID"     "PATHNAME"   "REACTOMEID"
head(keys(reactome.db))
[1] "1"         "10"        "100"       "1000"      "10000"     "100000160"
my_keys <- head(keys(reactome.db))

AnnotationDbi::select(
  reactome.db,
  keys = my_keys,
  columns="REACTOMEID",
  keytype="ENTREZID"
) |>
  head()
'select()' returned 1:many mapping between keys and columns
  ENTREZID    REACTOMEID
1        1  R-HSA-109582
2        1  R-HSA-114608
3        1  R-HSA-168249
4        1  R-HSA-168256
5        1 R-HSA-6798695
6        1   R-HSA-76002

GO.db

There are also other annotation packages, such as GO.db, which contains a set of annotation maps describing the entire Gene Ontology that we can probe in the same manner.

To begin, install the {GO.db} package.

BiocManager::install("GO.db")

Load the package.

packageVersion('GO.db')
[1] '3.19.1'
suppressPackageStartupMessages(library(GO.db))
class(GO.db)
[1] "GODb"
attr(,"package")
[1] "AnnotationDbi"

Use the same functions as before on the AnnotationDbi object.

keytypes(GO.db)
[1] "DEFINITION" "GOID"       "ONTOLOGY"   "TERM"      
my_go_keys <- head(keys(GO.db))
my_go_keys
[1] "GO:0000001" "GO:0000002" "GO:0000003" "GO:0000006" "GO:0000007"
[6] "GO:0000009"
AnnotationDbi::select(
  GO.db,
  keys = my_go_keys,
  columns=c("GOID", "TERM", "ONTOLOGY"),
  keytype="GOID"
)
'select()' returned 1:1 mapping between keys and columns
        GOID                                                     TERM ONTOLOGY
1 GO:0000001                                mitochondrion inheritance       BP
2 GO:0000002                         mitochondrial genome maintenance       BP
3 GO:0000003                                             reproduction       BP
4 GO:0000006    high-affinity zinc transmembrane transporter activity       MF
5 GO:0000007 low-affinity zinc ion transmembrane transporter activity       MF
6 GO:0000009                   alpha-1,6-mannosyltransferase activity       MF

GenomicFeatures

For genome annotation purposes, one may be interested in the gene models, which are provided in the TxDb.Hsapiens.UCSC.hg38.knownGene package.

BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")

Load the package.

packageVersion('TxDb.Hsapiens.UCSC.hg38.knownGene')
[1] '3.18.0'
suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene))
class(TxDb.Hsapiens.UCSC.hg38.knownGene)
[1] "TxDb"
attr(,"package")
[1] "GenomicFeatures"

To see what columns can be returned use the columns() function

columns(TxDb.Hsapiens.UCSC.hg38.knownGene)
 [1] "CDSCHROM"   "CDSEND"     "CDSID"      "CDSNAME"    "CDSPHASE"  
 [6] "CDSSTART"   "CDSSTRAND"  "EXONCHROM"  "EXONEND"    "EXONID"    
[11] "EXONNAME"   "EXONRANK"   "EXONSTART"  "EXONSTRAND" "GENEID"    
[16] "TXCHROM"    "TXEND"      "TXID"       "TXNAME"     "TXSTART"   
[21] "TXSTRAND"   "TXTYPE"    

Selecting.

AnnotationDbi::select(
  TxDb.Hsapiens.UCSC.hg38.knownGene,
  keys = head(keys(TxDb.Hsapiens.UCSC.hg38.knownGene), 2),
  columns=c('GENEID', 'TXCHROM', 'TXSTART', 'TXEND', 'TXID'),
  keytype="GENEID"
)
'select()' returned 1:many mapping between keys and columns
   GENEID   TXID TXCHROM  TXSTART    TXEND
1       1 228947   chr19 58345178 58347634
2       1 228948   chr19 58345183 58353492
3       1 228949   chr19 58346854 58356225
4       1 228950   chr19 58346858 58353491
5       1 228951   chr19 58346860 58347657
6       1 228952   chr19 58348466 58362751
7       1 228953   chr19 58350594 58353129
8       1 228954   chr19 58353021 58356083
9      10 104459    chr8 18386311 18388323
10     10 104461    chr8 18391282 18401218
11     10 104462    chr8 18391287 18400993

Combining databases

We can combine different annotation packages too; say I have two genes LRRK2 and PINK1, and I want their genomic locations as well as OMIM ids.

my_genes <- c("LRRK2","PINK1")
 
AnnotationDbi::select(
  org.Hs.eg.db,
  keys = my_genes,
  columns=c("ENTREZID", "SYMBOL","OMIM"),
  keytype="SYMBOL"
) -> a
'select()' returned 1:many mapping between keys and columns
a 
  SYMBOL ENTREZID   OMIM
1  LRRK2   120892 607060
2  LRRK2   120892 609007
3  PINK1    65018 605909
4  PINK1    65018 608309
# symbol and transcript location
AnnotationDbi::select(
  TxDb.Hsapiens.UCSC.hg38.knownGene,
  keys = a$ENTREZID,
  columns=c('GENEID', 'TXCHROM', 'TXSTART', 'TXEND', 'TXID'),
  keytype="GENEID"
) -> b
'select()' returned many:many mapping between keys and columns
head(b)
  GENEID   TXID TXCHROM  TXSTART    TXEND
1 120892 148742   chr12 40196744 40283952
2 120892 148744   chr12 40224977 40242782
3 120892 148745   chr12 40224997 40367077
4 120892 148746   chr12 40224997 40369285
5 120892 148747   chr12 40225011 40304976
6 120892 148748   chr12 40225080 40369259
# change GENEID column to ENTREZID for consistency and for merging
names(b) <- c('ENTREZID', 'TXID', 'TXCHROM', 'TXSTART', 'TXEND')
c <- merge(a, b, 'ENTREZID')
head(c)
  ENTREZID SYMBOL   OMIM   TXID TXCHROM  TXSTART    TXEND
1   120892  LRRK2 607060 148742   chr12 40196744 40283952
2   120892  LRRK2 607060 148744   chr12 40224977 40242782
3   120892  LRRK2 607060 148745   chr12 40224997 40367077
4   120892  LRRK2 607060 148746   chr12 40224997 40369285
5   120892  LRRK2 607060 148747   chr12 40225011 40304976
6   120892  LRRK2 607060 148748   chr12 40225080 40369259

Conclusions

One important point to be aware of when using the Bioconductor annotation packages is to appreciate the lag between updates to the packages. From the FAQ on the Bioconductor website:

Different sources take different approaches to managing annotations. The annotation packages in Bioconductor are based on downloads obtained shortly before each Bioconductor release, and so can lag by six months (at the end of the release cycle) compared to on-line resources. The advantage of this approach is that the annotations do not change unexpectedly during development of an analysis, while the disadvantage is that the resource is not quite up-to-date with current understanding.


sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

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       

time zone: Etc/UTC
tzcode source: system (glibc)

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

other attached packages:
 [1] TxDb.Hsapiens.UCSC.hg38.knownGene_3.18.0
 [2] GenomicFeatures_1.56.0                  
 [3] GenomicRanges_1.56.2                    
 [4] GenomeInfoDb_1.40.1                     
 [5] GO.db_3.19.1                            
 [6] reactome.db_1.88.0                      
 [7] org.Hs.eg.db_3.19.1                     
 [8] AnnotationDbi_1.66.0                    
 [9] IRanges_2.38.1                          
[10] S4Vectors_0.42.1                        
[11] Biobase_2.64.0                          
[12] BiocGenerics_0.50.0                     
[13] workflowr_1.7.1                         

loaded via a namespace (and not attached):
 [1] blob_1.2.4                  Biostrings_2.72.1          
 [3] bitops_1.0-7                fastmap_1.2.0              
 [5] RCurl_1.98-1.14             GenomicAlignments_1.40.0   
 [7] promises_1.3.0              XML_3.99-0.16.1            
 [9] digest_0.6.37               lifecycle_1.0.4            
[11] processx_3.8.4              KEGGREST_1.44.1            
[13] RSQLite_2.3.7               magrittr_2.0.3             
[15] compiler_4.4.0              rlang_1.1.4                
[17] sass_0.4.9                  tools_4.4.0                
[19] utf8_1.2.4                  yaml_2.3.8                 
[21] rtracklayer_1.64.0          knitr_1.47                 
[23] S4Arrays_1.4.1              bit_4.0.5                  
[25] curl_5.2.1                  DelayedArray_0.30.1        
[27] abind_1.4-5                 BiocParallel_1.38.0        
[29] grid_4.4.0                  fansi_1.0.6                
[31] git2r_0.33.0                SummarizedExperiment_1.34.0
[33] cli_3.6.3                   rmarkdown_2.27             
[35] crayon_1.5.2                rstudioapi_0.16.0          
[37] httr_1.4.7                  rjson_0.2.21               
[39] DBI_1.2.3                   cachem_1.1.0               
[41] stringr_1.5.1               zlibbioc_1.50.0            
[43] parallel_4.4.0              XVector_0.44.0             
[45] restfulr_0.0.15             matrixStats_1.3.0          
[47] vctrs_0.6.5                 Matrix_1.7-0               
[49] jsonlite_1.8.8              callr_3.7.6                
[51] bit64_4.0.5                 jquerylib_0.1.4            
[53] glue_1.7.0                  codetools_0.2-20           
[55] ps_1.7.6                    stringi_1.8.4              
[57] later_1.3.2                 BiocIO_1.14.0              
[59] UCSC.utils_1.0.0            tibble_3.2.1               
[61] pillar_1.9.0                htmltools_0.5.8.1          
[63] GenomeInfoDbData_1.2.12     R6_2.5.1                   
[65] rprojroot_2.0.4             evaluate_0.24.0            
[67] lattice_0.22-6              png_0.1-8                  
[69] Rsamtools_2.20.0            memoise_2.0.1              
[71] httpuv_1.6.15               bslib_0.7.0                
[73] Rcpp_1.0.12                 SparseArray_1.4.8          
[75] whisker_0.4.1               xfun_0.44                  
[77] fs_1.6.4                    MatrixGenerics_1.16.0      
[79] getPass_0.2-4               pkgconfig_2.0.3