BgeeDB, an R package for retrieval of curated expression datasets and for gene list expression localization enrichment tests

Requirements

Package installation

source("https://bioconductor.org/biocLite.R")
biocLite("BgeeDB")
# load the library
library(BgeeDB)

Data download and import of normalized expression levels

The first step of data retrieval is to initialize a new Bgee reference class object, for a targeted species and data type. Normalized expression levels are currently available in the BgeeDB package for two data types: Affymetrix microarrays and Illumina RNA-seq. The list of species available in the Bgee database for each data type, along with their NCBI taxonomy IDs and common names can be obtained with the listBgeeSpecies() function. By default, data will be downloaded from the latest Bgee release, but this can be changed with the release argument.

Next, the functions getAnnotation(), getData(), and formatData() can be called to respectively download the annotations of datasets, download the actual expression data, and reformat the expression data for more convenient use. Of note, BgeeDB creates a directory to store the downloaded annotation files and datasets, by default in the user’s R working directory, but this can be changed with the pathToData argument. These versioned cached files make it faster for the user to return to previously used data and allow for offline work.

Microarray dataset retrieval. In the following example, we look for a microarray dataset in mouse (Mus musculus), spanning multiple early developmental stages, including zygote. At the time of publication the latest Bgee release is 13.2, so if one needs to strictly reproduce the output of the code below in the future, the release="13.2" argument needs to be added when creating the Bgee object (see Supplementary file S1 and Supplementary file S2).

# specify species and data type
bgee.affymetrix <- Bgee$new(species="Mus_musculus", dataType="affymetrix")

# retrieve annotation of all mouse affymetrix datasets in Bgee
annotation.bgee.mouse.affymetrix <- getAnnotation(bgee.affymetrix)
str(annotation.bgee.mouse.affymetrix)

This creates a list of two data frames, one including the annotation of experiments, and one including the annotation of each individual sample, i.e., hybridized microarray chip. For mouse, there are 694 Affymetrix experiments and 6,077 samples available in Bgee release 13. Anatomical structures and developmental stages are annotated using the Uberon ontology^29,30. Below, we are selecting the experiments for which at least one sample is annotated to the zygote stage (UBERON:0000106).

# retrieve annotations of samples and experiments
sample.annotation <- annotation.bgee.mouse.affymetrix$sample.annotation
experiment.annotation <- annotation.bgee.mouse.affymetrix$experiment.annotation

# list experiments including a zygote sample
selected.experiments <- unique(sample.annotation$Experiment.ID[sample.annotation$Stage.ID == "UBERON:0000106"])
experiment.annotation[experiment.annotation$Experiment.ID %in% selected.experiments,]

# stages sampled in each of these experiments
unique(sample.annotation[sample.annotation$Experiment.ID %in% selected.experiments, c(1,6)])

This yields three microarray experiments, with accessions GSE1749, E-MEXP-51 and GSE18290. Among these, the accession E-MEXP-51, submitted to ArrayExpress by Wang and colleagues³¹, includes samples from more developmental stages than the other two, so we use this in the next steps. For this experiment, raw data were available from ArrayExpress, so samples were fully normalized with gcRMA³² version 2.40.0 through the Bgee pipeline.

# List all samples from E-MEXP-51 in Bgee
sample.annotation[sample.annotation$Experiment.ID == "E-MEXP-51",]

The experiment includes 35 samples that passed Bgee quality controls. They originate from 12 developmental stages: primary and secondary oocyte, zygote, early, mid and late 2-cells embryo, 4-cells embryo, 8-cells embryo, 16-cells embryo, early, mid and late blastocyst, although the developmental stages ontology used is not precise enough yet to differentiate some of these conditions: the early, mid and late 2-cells stages are annotated as Theiler stage 2 embryo, and the 4-cells and 8-cells stages are annotated as Theiler stage 3 embryo. All samples were hybridized to the Affymetrix GeneChip Murine Genome U74Av2 microarray. Let us download the normalized probesets intensities measured for all samples.

data.E.MEXP.51 <- getData(bgee.affymetrix, experimentId="E-MEXP-51")
head(data.E.MEXP.51)

The resulting data frame lists for each sample (column “Chip.ID”), the 9,017 probesets on the microarray (column “Probeset.ID”), their mapping to Ensembl gene IDs³³ (column “Gene.ID”), their logged normalized intensities (column “Log.of.normalized.signal.intensity”), and a presence/absence call and quality (columns “Detection.flag” and “Detection.quality”).

As this format might not be the most convenient for downstream processing of an expression dataset, we offer the formatData() function, which creates an ExpressionSet object including the expression data matrix, the probesets annotation to Ensembl genes and the samples' anatomical structure and stage annotation into (assayData, featureData and phenoData slots respectively). This object class is of standard use in numerous Bioconductor packages.

data.E.MEXP.51.formatted <- formatData(bgee.affymetrix, data.E.MEXP.51, callType="all", stats="intensities")
data.E.MEXP.51.formatted
# matrix of expression intensities
head(exprs(data.E.MEXP.51.formatted))
# annotation of samples
pData(data.E.MEXP.51.formatted)
# annotation of probesets
head(fData(data.E.MEXP.51.formatted))

The callType option of the formatData() function could alternatively be set to present or present high quality to display only the intensities of probesets detected as actively expressed.

The result is a nicely formatted Bioconductor object including expression data and their annotations, ready to be used for downstream analysis with other Bioconductor packages.

RNA-seq dataset retrieval. We now search Bgee for a RNA-seq dataset sampling brain and liver tissues (Uberon Ids UBERON:0000955 and UBERON:0002107 respectively) in macaque (Macaca mulatta), and including multiple biological replicates for each tissue.

# specify species and data type
bgee.rnaseq <- Bgee$new(species="Macaca_mulatta", dataType="rna_seq")

# retrieve annotations of RNA-seq samples and experiments
annotation.bgee.macaque.rna.seq <- getAnnotation(bgee.rnaseq)
sample.annotation <- annotation.bgee.macaque.rna.seq$sample.annotation
experiment.annotation <- annotation.bgee.macaque.rna.seq$experiment.annotation

# list experiments including both brain and liver samples
selected.experiments <- intersect(unique(sample.annotation$Experiment.ID[sample.annotation$Anatomical.entity.ID == "UBERON:0000955"]),
unique(sample.annotation$Experiment.ID[sample.annotation$Anatomical.entity.ID == "UBERON:0002107"]))
experiment.annotation[experiment.annotation$Experiment.ID %in% selected.experiments,]

# check whether experiments include biological replicates
sample.annotation[sample.annotation$Experiment.ID %in%
selected.experiments & (sample.annotation$Anatomical.entity.ID == "UBERON:0000955" | 
sample.annotation$Anatomical.entity.ID == "UBERON:0002107"), 1:6]

Accessions GSE41637³⁴ and GSE30352³⁵ both include biological replicates for brain and liver. We focus on GSE41637 for the next steps since it includes three replicates of each tissue, vs. only two for GSE30352. We download the dataset and reformat it to obtain an ExpressionSet including counts of mapped reads on each Ensembl gene for each sample.

data.GSE41637 <- getData(bgee.rnaseq, experimentId="GSE41637")
data.GSE41637.formatted <- formatData(bgee.rnaseq, data.GSE41637, callType="all", stats="counts")
data.GSE41637.formatted

Instead of mapped read counts, it is also possible to fill the data matrix with expression levels in the RPKM unit (reads per kilobase per million reads), using the option stats="rpkm". In the next Bgee release (release 14), it will be possible to obtain expression levels in the TPM unit (transcript per million)^36,37 from pseudo-mapping of reads computed in Bgee using the Kallisto software³⁸.

Presence/absence calls retrieval. It is often difficult to compare expression levels across species³⁹, and even within species, across datasets generated by different experimenters or laboratories^40–42. Batch effects have indeed been shown to impact extensively gene expression levels, confounding biological signal differences.

Encoding gene expression as present or absent in a sample allows a more robust comparison across such conditions. In addition to retrieving RNA-seq and Affymetrix quantitative expression levels, BgeeDB also allows to retrieve calls of presence or absence of expression computed in the Bgee database for each gene (RNA-seq) or probeset (Affymetrix), in the column “Detection.flag” of the data.E.MEXP.51 and data.GSE41637 objects created above. And interestingly, expression calls are also available in Bgee for ESTs and in situ hybridization data, as well as for the consensus of the four data types for each triplet “gene / tissue / developmental stage”.

A powerful use of these expression calls is the anatomical expression enrichment test “TopAnat”. TopAnat uses a similar approach to Gene Ontology enrichment tests²⁶, but genes are associated to the anatomical structures where they display expression, instead of to their functional classification. These tests allow detecting where a set of genes is preferentially expressed as compared to a background universe (Roux J., Seppey M., Sanjeev K., Rech de Laval V., Moret P., Artimo P., Duvaud S., Ioannidis V., Stockinger H., Robinson-Rechavi M., Bastian F.B.; unpublished report). We show an example of such an analysis in the section “Anatomical expression enrichment analysis” below.

Of note, the expression calls imported from BgeeDB can also be used for other downstream analyses. For example, when studying protein-protein interaction datasets, it might be biologically relevant to retain only interactions for which both members are expressed in the same tissues^43,44.

Downstream analysis examples

Clustering analysis. A variety of downstream analyses can be performed on the imported expression data. Below we detail an example of gene expression clustering analysis on the developmental time-series microarray experiment imported above. The analysis, performed with the Mfuzz package^45,46 (version 2.34.0 for this paper), aims at uncovering genes with similar expression profiles across development. We can readily start with the ExpressionSet object previously created.

# for simplicity, keep only one sample per condition
data.E.MEXP.51.formatted <- data.E.MEXP.51.formatted[, !duplicated(pData(data.E.MEXP.51.formatted)[2:5])]

# order developmental stages
data.E.MEXP.51.formatted <- data.E.MEXP.51.formatted[, c(5,8,9,3,2,1,4,7,6)]

# filter out rows with no variance
data.E.MEXP.51.formatted <- data.E.MEXP.51.formatted[apply(exprs(data.E.MEXP.51.formatted), 1, sd) != 0, ]

# Mfuzz clustering
biocLite("Mfuzz")
library(Mfuzz)
# standardize matric of expression data
z.mat <- standardise(data.E.MEXP.51.formatted)
# cluster data into 16 clusters
clusters <- mfuzz(z.mat, centers=16, m=1.25)

# visualizing clusters
mfuzz.plot2(z.mat, cl=clusters, mfrow=c(4,4), colo="fancy",
time.labels=row.names(pData(z.mat)), las=2, xlab="", ylab="Standardized expression level", x11=FALSE)

The resulting plot can be seen in Figure 1.

Figure 1.Standardized expression levels of 16 groups of microarray probesets, clustered according to their expression during mouse early development.

The x-axis displays sample names (column “Chip.ID” of the data.E.MEXP.51 object).

Differential expression analysis. Below, we detail a differential expression analysis, with the package edgeR^47,48 (version 3.16.1 for this paper), on the previously imported RNA-seq dataset of macaque tissues. We aim at isolating genes differentially expressed between brain and liver.

# differential expression analysis with edgeR
biocLite("edgeR")
library(edgeR)

# subset the dataset to brain and liver
brain.liver <- data.GSE41637.formatted[, pData(data.GSE41637.formatted)$Anatomical.entity.name %in% c("brain", "liver")]  

# filter out very lowly expressed genes
brain.liver.filtered <- brain.liver[rowSums(cpm(brain.liver) > 1) > 3, ]

# create edgeR DGElist object
dge <- DGEList(counts=brain.liver.filtered,
group=pData(brain.liver.filtered)$Anatomical.entity.name)
dge <- calcNormFactors(dge)
dge <- estimateCommonDisp(dge)
dge <- estimateTagwiseDisp(dge)
de <- exactTest(dge, pair=c("brain","liver"))
de.genes <- topTags(de, n=nrow(de))$table

# MA plot with DE genes highlighted
plotSmear(dge, de.tags=rownames(de.genes)[de.genes$FDR < 0.01], cex=0.3)

The resulting plot can be seen in Figure 2.

Figure 2.Mean-average (MA) plot of differential gene expression between brain and liver in macaque based on RNA-seq data.

Significantly differentially expressed genes (FDR < 1%) are highlighted in red.

Anatomical expression enrichment analysis

The loadTopAnatData() function loads the names of anatomical structures, and relationships between them, from the Uberon anatomical ontology (based on parent-child “is_a” and “part_of” relationships). It also loads a mapping from genes to anatomical structures, based on the presence calls of the genes in the targeted species. These calls come from a consensus of all data types specified in the input Bgee class object. We recommend to use all available data types (RNA-seq, Affymetrix, EST and in situ hybridization) for both genomic coverage and anatomical precision, which is the default behavior if no dataType argument is specified when the Bgee class object is created.

By default, presence calls of both high and low quality are used, which can be changed with the confidence argument of the loadTopAnatData() function. Finally, it is possible to specify the developmental stage under consideration, with the stage argument. By default expression calls generated from samples of all developmental stages are used, which is equivalent to specifying stage="UBERON:0000104" (“life cycle”, the root of the stage ontology). Data are stored in versioned tab-separated cached files that will be read again if a query with the exact same parameters is launched later, to save time and server resources, and to work offline.

In this example, we use expression calls for zebrafish genes using all sources of expression data.

bgee.topanat <- Bgee$new(species="Danio_rerio")
myTopAnatData <- loadTopAnatData(bgee.topanat)
str(myTopAnatData)

We look at the expression localization of the genes with an annotated phenotype related to pectoral fin (i.e., genes which upon knock-out or knock-down led to abnormal phenotypes of pectoral fin or its components). Zebrafish phenotypic data are available from the ZFIN database⁴⁹ and integrated into the Ensembl database⁵⁰. We thus retrieve the targeted genes using the biomaRt⁵¹ Bioconductor package (version 2.30.0 for this paper).

biocLite("biomaRt")
library(biomaRt)

# zebrafish data in Ensembl 85 (stable link)
ensembl <- useMart("ENSEMBL_MART_ENSEMBL",
dataset="drerio_gene_ensembl", host="jul2016.archive.ensembl.org")

# get the mapping of Ensembl genes to phenotypes
genesToPhenotypes <- getBM(filters=c("phenotype_source"), value=c("ZFIN"),
attributes=c("ensembl_gene_id","phenotype_description"), mart=ensembl)

# select phenotypes related to pectoral fin
myPhenotypes <- grep("pectoral fin", unique(genesToPhenotypes$phenotype_description), value=T)

# select the genes annotated to select phenotypes
myGenes <- unique(genesToPhenotypes$ensembl_gene_id[genesToPhenotypes$phenotype_description %in% myPhenotypes])

This gives a list of 150 zebrafish genes implicated in the development and function of pectoral fin. The next step of the analysis relies on the topGO Bioconductor package. We prepare a modified topGOdata object allowing to handle the Uberon anatomical ontology instead of the Gene Ontology, and perform a GO-like enrichment test for anatomical terms. As for a classical topGO analysis, we need to prepare a vector including all background genes, and with values 0 or 1 depending if genes are part of the foreground or not. The choice of background is very important since the wrong background can lead to spurious results in enrichment tests⁵². Here we choose as background all zebrafish Ensembl genes with an annotated phenotype from ZFIN.

# prepare the gene list vector 
geneList <- factor(as.integer(unique(genesToPhenotypes$ensembl_gene_id) %in% myGenes))
names(geneList) <- unique(genesToPhenotypes$ensembl_gene_id)
summary(geneList)

# prepare the topAnat object based on topGO
myTopAnatObject <- topAnat(myTopAnatData, geneList)

At this step, expression calls are propagated through the whole ontology (e.g., expression in the forebrain will also be counted as expression in the brain, the nervous system, etc). This can take some time, especially if the gene list is large.

Finally, we launch an enrichment test for anatomical terms. The functions of the topGO package can directly be used at this step. See the vignette of this package for more details²⁵. Here we use a Fisher test, coupled with the “weight” decorrelation algorithm.

results <- runTest(myTopAnatObject, algorithm='weight', statistic='fisher')

Finally, we implement a function to display results in a formatted table. By default anatomical structures are sorted by their test p-value, which is displayed along with the associated false discovery rate (FDR⁵³) and the enrichment fold. Sorting on other columns of the table (e.g., on decreasing enrichment folds) is possible with the ordering argument. Of note, it is debated whether a FDR correction is relevant on such enrichment test results, since tests on different terms of the ontologies are not independent. An interesting discussion can be found in the vignette of the topGO package.

# retrieve anatomical structures enriched at a 1% FDR threshold tableOver <- makeTable(myTopAnatData, myTopAnatObject, results, cutoff=0.01)

The 22 anatomical structures displaying a significant enrichment at a FDR threshold of 1% are show in Table 1. The first term is “paired limb/fin bud”, and the second “pectoral fin”. Other terms in the list, especially those with high enrichment folds, are clearly related to pectoral fins (e.g., “pectoral appendage cartilage tissue”), substructures of fins (e.g., “fin bone”), or located next to them (e.g., “ceratohyal cartilage”). This analysis shows that genes with phenotypic effects on pectoral fins are specifically expressed in or next to these structures. More generally, it proves the pertinence of TopAnat analysis for the characterization of lists of genes.