---
title: "Searching NCBI and Ensembl"
---
This lab guides you through a routine session with online biological databases.
The first sections give an overview of some of the existing internet databases
for sequence and general information search; the last part uses **R** to gain
programmatic access to annotation datasets through **Bioconductor** libraries.
By the end of this laboratory you will be able to:
- navigate the main **NCBI** resources (Entrez, Gene, Genome, PubMed) and extract
genome, gene, and literature information
- retrieve gene and transcript information from **Ensembl**
- look up protein sequences and structures in **UniProtKB/Swiss-Prot** and the **RCSB PDB**
- query genome-wide annotation from R using the Bioconductor `org.Hs.eg.db` package
- query NCBI live from R with **`rentrez`** — and know when to prefer a local
snapshot over a live query
::: {.callout-note title="Installing Bioconductor packages"}
This lab uses packages from **Bioconductor**. If your R installation is missing a
package that belongs to Bioconductor (rather than CRAN), install it like this:
```r
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("<packageName>") # a specific package
```
:::
```{r}
#| label: ncbi-setup
#| include: false
source("R/setup.R")
options(width = 75)
```
::: {.callout-note appearance="minimal" title="Lab at a glance"}
**Goal:** query biological databases by hand and, reproducibly, from R.
**Prerequisites:** [Introduction to R](intro-r.qmd). **Time:** ~2 h.
**Packages:** `org.Hs.eg.db`, `dplyr`, `rentrez` (and `biomaRt` for live Ensembl
access, shown not run). **Deliverable:** a lab report.
:::
## NCBI {#sec-ncbi}
NCBI stands for the *National Center for Biotechnology Information*. It is a
division of the National Library of Medicine, which belongs to the National
Institutes of Health (NIH). Created in 1988 as a US national resource for
molecular-biology information, all of its databases are public. Besides
molecular-biology information, NCBI provides tools for analysing that data and is
active in basic research, databases, and education.
NCBI links many databases through the **Entrez** search system and provides
tools for sequence search, genome browsing, data download, 3D structure
visualisation, and literature discovery.
Important resources include **Bookshelf**, **Gene**, **Nucleotide**, **Protein**,
**Assembly**, **Taxonomy**, **PubMed**, **GEO**, **SRA**, **dbGaP**, **dbSNP**,
**ClinVar**, and the **Conserved Domain Database**. Useful analysis and download
tools include **BLAST**, **Genome Data Viewer (GDV)**, **NCBI Datasets**, **ORF
Finder**, and the browser-based **iCn3D** structure viewer.
### Genomes {#sec-genomes}
NCBI Datasets is the current entry point for finding assemblies and downloading
genome, transcript, protein, and annotation files; Genome Data Viewer (GDV)
provides the interactive chromosome view (the older Map Viewer was retired in
2018). A good first exercise is to compare the *reference assemblies* of three
very different organisms.
::: {.callout-tip title="✅ Exercise 2.1 — Genomes"}
1. Open the [NCBI Datasets genome page](https://www.ncbi.nlm.nih.gov/datasets/genome/).
2. Search for the reference genomes of *Caenorhabditis elegans*, *Severe acute
respiratory syndrome coronavirus 2*, and *Mus musculus* and fill in a
comparison table like the one below (report it in your lab report):
| organism | accession | assembly level | genome size | annotated genes |
|----------|----------|----------------|-------------|-----------------|
| *C. elegans* | | | | |
| SARS-CoV-2 | | | | |
| *M. musculus* | | | | |
Explain why "number of chromosomes" is not a useful field for every viral or
draft assembly.
3. Open the *C. elegans* or mouse reference assembly in
[Genome Data Viewer](https://www.ncbi.nlm.nih.gov/gdv/). Choose a chromosome,
navigate to a gene, zoom from chromosome to exon scale, and identify the
assembly and RefSeq annotation release shown.
4. From NCBI Datasets, download a small data package containing the genome data
report and GFF3 annotation. Describe the files and identify which one contains
feature coordinates.
:::
### Gene search {#sec-gene-search}
From the NCBI entry point <https://www.ncbi.nlm.nih.gov/>, open the **Gene**
database and look for the human gene that encodes a receptor involved in
*familial hypercholesterolemia* (FH).
::: {.callout-tip title="✅ Exercise 2.2 — Gene search"}
1. What is the full name of the gene? What is its official *Symbol*? Which
organisation standardises the official symbols?
2. Describe the *gene type*. What other possibilities exist for the *Gene Type*
field? (Hint: look at the help.)
3. On which chromosome is this gene? Give the *location* of the gene.
4. Search for and download the gene sequence. Compute its length in base pairs.
5. Search for other *loci* related to familial hypercholesterolemia
(see <https://www.ncbi.nlm.nih.gov/omim/143890>); write down their position
and length.
6. Describe what the *phenotypes* and *homology* fields mean on a gene entry.
7. To which databases do the links in the *Pathway* section connect? Describe the
site you reach and what it contains.
:::
### PubMed {#sec-pubmed}
PubMed indexes biomedical literature. Its field tags make searches reproducible
and help separate words in titles/abstracts from author affiliations.
::: {.callout-tip title="✅ Exercise 2.3 — PubMed"}
1. Search for what *COVID-19* is in PubMed / *Entrez* (Publications),
<https://pubmed.ncbi.nlm.nih.gov/>. Explain what this database is.
2. Compare the broad query *COVID-19 Wang* with the author-field query
`COVID-19 AND Wang[Author]`.
3. Entrez allows control over the search by specifying fields. Search for
`COVID-19[Title/Abstract] AND China[Affiliation]` in PubMed. Explain what the fields mean and what
other search options exist (the *help* button is your friend).
:::
## Ensembl {#sec-ensembl}
The Ensembl Project is based at the **European Molecular Biology Laboratory's
European Bioinformatics Institute (EMBL-EBI)**. It produces and integrates
genome annotations for vertebrates and other eukaryotes. Its website and APIs
expose genes, transcripts, comparative genomics, variation, and regulation data.
::: {.callout-tip title="✅ Exercise 2.4 — Ensembl"}
1. Point your browser to <https://www.ensembl.org>.
2. Look for the gene *LDLR* in *Homo sapiens*.
3. Within the gene entry, navigate to *Sequence* in the left bar; observe the
sequence and the marked exonic regions.
4. Open the transcript Ensembl marks as **canonical** and record its stable
transcript and protein IDs and lengths. Do not use "the first transcript":
display order can change between releases.
5. Export one exon in FASTA format and record the Ensembl release and human
assembly so your answer remains reproducible.
:::
## Protein databases {#sec-protein-databases}
PubMed can retrieve protein information, but dedicated protein databases are also
useful. A good example is UniProt (Swiss-Prot).
### UniProtKB (Swiss-Prot) {#sec-uniprot}
We will inspect the reviewed UniProtKB/Swiss-Prot record for the SARS-CoV-2 spike
glycoprotein. The protein is viral; *Homo sapiens* is its host, not its source
organism.
::: {.callout-tip title="✅ Exercise 2.5 — UniProtKB"}
1. Point a browser to <https://www.uniprot.org/>.
2. Search for accession **P0DTC2** (*SPIKE_SARS2*) and confirm that the organism
is SARS-CoV-2. (P59594 is the related spike protein from the 2003 SARS-CoV.)
3. Explain the function of the SARS-CoV-2 spike protein and distinguish the S1
receptor-binding and S2 membrane-fusion roles.
4. Retrieve the amino-acid protein sequence in FASTA format.
5. What is the molecular weight (mass) of the protein, in kDa?
:::
### Research Collaboratory for Structural Bioinformatics (RCSB) {#sec-rcsb}
RCSB PDB stores experimentally determined macromolecular structures and links to
computed structure models. Connect at <https://www.rcsb.org>.
::: {.callout-tip title="✅ Exercise 2.6 — RCSB"}
1. What does RCSB mean?
2. Search for PDB entry **6M0J**, the SARS-CoV-2 spike receptor-binding domain
bound to human ACE2.
3. Describe the differences between UniProtKB and the PDB. What type of information
can you obtain from RCSB, apart from the protein sequence?
4. Identify the chains corresponding to viral spike and human ACE2 (UniProt
**Q9BYF1**), the experimental method, resolution, and biological assembly.
Hint: the complex is a 1:1 construct (one RBD chain bound to one ACE2
chain), which helps you tell spike apart from ACE2 when more than two
chains are present in the asymmetric unit.
5. Inspect the interface in the browser-based Mol* viewer. Download the mmCIF
file (the archival standard); optionally compare it with the legacy PDB format
in [PyMOL](https://pymol.org).
:::
::: {.callout-note title="About the PDB file format"}
The primary information stored in the PDB archive consists of coordinate files
for biological molecules, listing the atoms in each protein and their 3D
locations. The PDB archive's standard text format is **PDBx/mmCIF**; the older
fixed-column PDB format remains available for compatible structures. A typical
legacy PDB-formatted file has a large *header* section summarising the protein,
citation information, and the details of the structure solution, followed by the
sequence and a long list of atoms and their coordinates. The archive also
contains the experimental observations used to determine the coordinates. Further
information about the format is available
[here](https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/introduction).
:::
## Programmatic access from R {#sec-programmatic-access}
Everything so far happened in a browser. In this section we move the same
queries into R, which is what makes them reproducible and automatable.
Through Bioconductor, R offers many local annotation packages (see the
[AnnotationData BiocViews](https://www.bioconductor.org/packages/release/BiocViews.html#___AnnotationData)),
while services such as `biomaRt` and `rentrez` provide version-aware access to
live databases. We start with a local snapshot and finish with live queries.
### `org.Hs.eg.db` {#sec-org-hs-eg-db}
`org.Hs.eg.db` is a Bioconductor package providing genome-wide annotation for
*Homo sapiens*, primarily based on mappings using Entrez Gene identifiers.
Common package families include organism-level mappings such as `org.Hs.eg.db`
and `org.Mm.eg.db`, assay annotation such as `hgu133plus2.db`, ontology metadata
in `GO.db`, transcript models in `TxDb.*`, and Ensembl-based models in `EnsDb.*`.
Always record the package version and genome assembly; identifier mappings and
annotations change between releases.
Most annotation packages use four standard access methods: `columns()`,
`keytypes()`, `keys()`, and `select()`.
Load `org.Hs.eg.db` and inspect its contents.[^suppress-note]
[^suppress-note]: `suppressPackageStartupMessages()` hides the package start-up
messages, which is convenient when generating a document.
```{r}
#| label: ls
suppressPackageStartupMessages({
library(org.Hs.eg.db)
})
ls("package:org.Hs.eg.db")
```
To learn which kinds of data are retrievable via `select()`, list the available
columns:
```{r}
#| label: columns
columns(org.Hs.eg.db)
```
You can retrieve help for a particular column with, e.g., `help("GO")`.
Query keys are available through `keytypes()`:
```{r}
#| label: keytypes
keytypes(org.Hs.eg.db)
```
Get keys of a particular type with `keys()` — for example the first gene symbols,
OMIM entries, and GO entries:
```{r}
#| label: keys-first
head(keys(org.Hs.eg.db, keytype = "SYMBOL"))
head(keys(org.Hs.eg.db, keytype = "OMIM"))
head(keys(org.Hs.eg.db, keytype = "GO"))
```
Now define a general query on the `org.Hs.eg.db` object through
`select()`.[^mapids-note] You specify: **keys** (what to look up), **columns**
(what to return), and **keytype** (the type of key you pass in). For example,
query all *Gene Ontology* entries for coagulation Factor V, whose gene symbol is
`F5`:
[^mapids-note]: Use `mapIds()` if you need only a single column.
```{r}
#| label: f5-call
select(org.Hs.eg.db, keys = "F5",
columns = c("GO", "GENENAME"), keytype = "SYMBOL")
```
`select()` scales to the whole genome, and the *shape* of the result teaches
you something about annotation itself. As a worked example, count how many
OMIM records each UniProt accession links to, and look at the distribution
(@fig-omim-prot): most proteins carry none or a few cross-references, while a
long tail of heavily curated proteins links to dozens of OMIM entries.
```{r}
#| label: fig-omim-prot
#| fig-cap: "Number of distinct OMIM identifiers cross-referenced per UniProt accession (log scale)."
#| fig-alt: "Histogram with a tall spike at zero cross-references and a long right tail up to 18 OMIM entries per protein."
#| fig-width: 6
#| fig-height: 3.4
#| message: false
library(dplyr)
omim_uniprot <- AnnotationDbi::select(org.Hs.eg.db,
keys = keys(org.Hs.eg.db, keytype = "OMIM"),
columns = "UNIPROT", keytype = "OMIM") |>
filter(!is.na(UNIPROT)) |>
distinct(OMIM, UNIPROT) |>
count(UNIPROT, name = "n_omim")
all_uniprot <- data.frame(UNIPROT = keys(org.Hs.eg.db, keytype = "UNIPROT")) |>
left_join(omim_uniprot, by = "UNIPROT")
all_uniprot$n_omim[is.na(all_uniprot$n_omim)] <- 0
ggplot(all_uniprot, aes(n_omim)) +
geom_histogram(bins = 19, fill = "#0e8d7c") +
scale_y_log10() +
labs(x = "OMIM cross-references per UniProt accession",
y = "proteins (log scale)")
```
::: {.callout-note title="A note on name clashes"}
`dplyr` also has a `select()` function, and whichever package was loaded *last*
wins. Once `library(dplyr)` has run, call the annotation query with its full
name `AnnotationDbi::select(...)` (or `org.Hs.eg.db`'s own `select`) to be
unambiguous — this is the standard way R resolves name conflicts between
packages.
:::
::: {.callout-tip title="✅ Exercise 2.7 — org.Hs.eg.db"}
Using the *genome-wide annotation for Human* package (`org.Hs.eg.db`), work
through these questions in order — each builds on the previous one. As you go,
keep a single tidy data frame of OMIM × UniProt (or OMIM × SYMBOL) pairs and
reuse it rather than re-querying the package:
1. How many distinct OMIM identifiers does the annotation package contain?
(Note that OMIM identifiers can represent gene *and* phenotype records.)
2. What fraction of the UniProt accessions in `org.Hs.eg.db` have at least one
OMIM cross-reference? And what fraction of gene symbols (SYMBOL) do?
3. List the five proteins linked to the largest number of distinct OMIM
identifiers, with their gene names. What do these proteins have in common?
4. Compute the number of genes (SYMBOL) and the number of proteins (UNIPROT)
in the package. Why do the two counts differ — and why is neither count
"the number of genes in the human genome"?
5. Build a table of the gene–protein correspondence. Which are the 10 genes
with the largest number of related proteins? Find all UniProt IDs for the
first gene in that list, and explain why it has so many.
:::
::: {.content-hidden when-format="pdf" unless-profile="solutions"}
::: {.callout-tip collapse="true" title="Solution 2.7 — org.Hs.eg.db"}
Start from all OMIM keys, and build the OMIM × UniProt cross-reference once:
```{r}
#| label: sol-omim
#| message: false
library(dplyr)
omim <- keys(org.Hs.eg.db, keytype = "OMIM")
length(omim) # 1. number of OMIM entries
uniprot <- keys(org.Hs.eg.db, keytype = "UNIPROT")
g <- AnnotationDbi::select(org.Hs.eg.db, keys = omim,
columns = "UNIPROT", keytype = "OMIM") |>
filter(!is.na(UNIPROT)) |>
distinct(OMIM, UNIPROT)
# 2. fraction of annotated proteins with at least one OMIM cross-reference
n_distinct(g$UNIPROT) / length(uniprot)
```
So about `r round(100 * n_distinct(g$UNIPROT) / length(uniprot))` % of UniProt
accessions carry at least one OMIM link — OMIM coverage of *reviewed, disease
curated* proteins is extensive, but a 10 % tail has none (housekeeping and
poorly characterised proteins).
```{r}
#| label: sol-top-prot
# 3. proteins associated with the most OMIM identifiers
top <- g |>
count(UNIPROT, name = "n_omim", sort = TRUE) |>
slice_head(n = 5)
top |>
left_join(AnnotationDbi::select(org.Hs.eg.db, keys = top$UNIPROT,
columns = "SYMBOL", keytype = "UNIPROT"),
by = "UNIPROT")
```
They are all **COL2A1** isoforms — collagen II is heavily cross-referenced
because a large number of skeletal-dysplasia phenotypes map to the same
protein product.
```{r}
#| label: sol-gene-prot
# 4. genes vs proteins
n_sym <- length(keys(org.Hs.eg.db, keytype = "SYMBOL"))
n_up <- length(uniprot)
c(genes = n_sym, proteins = n_up)
# 5. gene-protein correspondence, top-10 genes by number of proteins
gp <- AnnotationDbi::select(org.Hs.eg.db,
keys = keys(org.Hs.eg.db, keytype = "ENTREZID"),
columns = c("SYMBOL", "UNIPROT"),
keytype = "ENTREZID") |>
filter(!is.na(SYMBOL), !is.na(UNIPROT)) |>
distinct()
gp |>
count(SYMBOL, name = "n_uniprot", sort = TRUE) |>
slice_head(n = 10)
# all UniProt accessions for the top gene
gp |> filter(SYMBOL == "HLA-DRB1") |> pull(UNIPROT)
```
The counts differ for two reasons: one gene can map to several UniProt
accessions (isoforms, reviewed vs unreviewed entries) and — the bigger effect
in the other direction — many `SYMBOL` keys correspond to pseudogenes, RNA
genes, and other loci that have no protein product at all. Neither number is
"the number of genes in the human genome": the package aggregates several
databases, each with its own scope. The top-10 list is dominated by the **HLA**
locus: `HLA-DRB1` alone has hundreds of accessions because each serotype and
allele variant of the highly polymorphic HLA proteins gets its own reviewed
UniProt entry.
:::
:::
### Live database access with `rentrez` {#sec-ncbi-programmatic}
`org.Hs.eg.db` is a **local snapshot**: it installs once, works offline, and is
perfectly reproducible — but it only knows what was true at its release date,
and NCBI and Ensembl keep moving. For **live** queries use the packages that
drive the same web APIs you clicked through above: `rentrez` speaks the NCBI
E-utilities behind Entrez, and `biomaRt` speaks to the Ensembl BioMart behind
the Ensembl browser. A pragmatic workflow is to prototype live, then pin a
snapshot (package version or API date) for anything you publish.
`rentrez` follows the web workflow you did by hand: `entrez_search()` finds
records, `entrez_summary()` fetches their metadata, `entrez_fetch()` downloads
full records. Let us revisit the LDLR gene from Exercise 2.2 — this time from
R:
```{r}
#| label: rentrez-ldlr
#| cache: true
library(rentrez)
ldlr <- entrez_search(db = "gene", term = "LDLR[Gene Name] AND Homo sapiens[Organism]")
ldlr$ids
```
```{r}
#| label: rentrez-summary
#| cache: true
entrez_summary("gene", id = ldlr$ids)
```
The summary list behaves like a data frame: `extract_from_esummary()` pulls
the fields you need. Note how the *live* answer — description, chromosome,
coordinates, last update — matches what you found on the web page:
```{r}
#| label: rentrez-extract
#| cache: true
ldlr_sum <- entrez_summary("gene", id = ldlr$ids)
extract_from_esummary(ldlr_sum, c("name", "description", "chromosome"))
```
The full record (the equivalent of the web page's *Genomic regions, transcripts,
and products* table) comes from `entrez_fetch()`:
```{r}
#| label: rentrez-fetch
#| cache: true
ldlr_tab <- entrez_fetch("gene", id = ldlr$ids, rettype = "gene_table", retmode = "text")
substr(ldlr_tab, 1, 400) # first lines of the gene table
```
The same search pattern drives PubMed. The query from Exercise 2.3 becomes an
`entrez_search()` call, and `retmax` caps how many records come back:
```{r}
#| label: rentrez-pubmed
#| cache: true
es <- entrez_search(db = "pubmed",
term = "COVID-19[Title/Abstract] AND China[Affiliation]",
retmax = 3)
es$count # total number of matching publications
es$ids # the first 3 PMIDs
```
::: {.callout-tip title="✅ Exercise 2.8 — rentrez"}
1. Fetch the *titles* and *first authors* of the three publications found
above (`extract_from_esummary()` on the PubMed summaries; the fields are
`title` and `sortfirstauthor`).
2. The first `rentrez` query in this section used a *search*. Repeat it the
other way round: fetch the full **gene_table** record for Entrez Gene ID
**3949** directly with `entrez_fetch()` (hint: `id = 3949`).
3. `es$count` changes daily. Which of the two numbers you computed in this
section would you report in a paper, and how would you make it
reproducible?
:::
::: {.content-hidden when-format="pdf" unless-profile="solutions"}
::: {.callout-tip collapse="true" title="Solution 2.8 — rentrez"}
```{r}
#| label: sol-rentrez
#| cache: true
pm <- entrez_summary("pubmed", id = es$ids)
extract_from_esummary(pm, c("sortfirstauthor", "title"))
# 2. direct fetch by known ID — no search step
entrez_fetch("gene", id = 3949, rettype = "gene_table", retmode = "text") |>
substr(1, 200)
```
3. Never report `es$count` bare: it grows with every indexed paper. Report the
query string, the database, the **date of the query** (and ideally the API
version), or snapshot the retrieved records so the analysis is fixed.
:::
:::
`biomaRt` is the Ensembl counterpart (it needs network access and an Ensembl
mirror, so it is shown here rather than run). It queries the same Ensembl
database you browsed in Exercise 2.4:
```r
library(biomaRt)
mart <- useEnsembl("genes", dataset = "hsapiens_gene_ensembl")
getBM(attributes = c("hgnc_symbol", "chromosome_name", "start_position", "end_position"),
filters = "hgnc_symbol", values = "LDLR", mart = mart)
```
## Further reading
- [NCBI Datasets documentation](https://www.ncbi.nlm.nih.gov/datasets/docs/) and the [BLAST topics guide](https://blast.ncbi.nlm.nih.gov/doc/blast-topics/).
- [Ensembl Help](https://www.ensembl.org/info/) — gene/transcript pages and BioMart.
- The `AnnotationDbi` vignette (`vignette("IntroToAnnotationPackages")`) for `select()`/`mapIds()`.
::: {.content-visible when-format="html" unless-profile="solutions"}
Worked solutions for the exercises are in the collapsed callouts below each exercise.
:::
::: {.content-visible when-format="html" when-profile="solutions"}
Solutions are collected in the [Solutions appendix](../appendix/solutions.qmd) (the org.Hs.eg.db answer is expanded above).
:::
::: {.content-visible when-format="pdf" unless-profile="solutions"}
Worked solutions for this lab are omitted from this PDF — find them as collapsed callouts in the HTML edition.
:::
## Report {#sec-report-ncbi}
Write a laboratory report following the structure of this text. It should contain
an introduction, the goals of each section, and a description of the results of
each part. Include any graphics relevant to making your work understandable.
**The lab should be delivered within two weeks.**
