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 2506265. 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 | 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.
To begin, install the {org.Hs.eg.db} package.
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")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"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 pseudogeneWithin 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 3Using 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).
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 MFFor 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 18400993We 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 columnsa 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 columnshead(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 40369259One 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] org.Hs.eg.db_3.19.1
[7] AnnotationDbi_1.66.0
[8] IRanges_2.38.1
[9] S4Vectors_0.42.1
[10] Biobase_2.64.0
[11] BiocGenerics_0.50.0
[12] 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