---
title: "BgeeDB, an R package for retrieval of curated expression datasets and for gene list enrichment tests"
author: "Andrea Komljenovic, Julien Roux, Marc Robinson-Rechavi, Frederic Bastian"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{BgeeDB, an R package for retrieval of curated expression datasets and for gene list enrichment tests}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```BgeeDB``` is a collection of functions to import data from the Bgee database (<https://bgee.org/>) directly into R, and to facilitate downstream analyses, such as gene set enrichment test based on expression of genes in anatomical structures. Bgee provides annotated and processed expression data and expression calls from curated wild-type healthy samples, from human and many other animal species.
 
The package retrieves the annotation of RNA-seq, single cell full length RNA-seq and Affymetrix experiments integrated into the Bgee database, and downloads into R the quantitative data and expression calls produced by the Bgee pipeline. The package also allows to run GO-like enrichment analyses based on anatomical terms, where genes are mapped to anatomical terms by expression patterns, based on the ```topGO``` package. This is the same as the TopAnat web-service available at (<https://www.bgee.org/analysis/top-anat/>), but with more flexibility in the choice of parameters and developmental stages.

In summary, the BgeeDB package allows to: 
* 1. List annotation of bulk RNA-seq, single cell RNA-seq and microarray data available the Bgee database
* 2. Download the processed gene expression data available in the Bgee database
* 3. Download the gene expression calls and use them to perform TopAnat analyses 

The pipeline used to generate Bgee expression data is documented and publicly available at (<https://github.com/BgeeDB/bgee_pipeline>)

If you find a bug or have any issues to use ```BgeeDB``` please write a bug report in our own GitHub issues manager available at (<https://github.com/BgeeDB/BgeeDB_R/issues>)

## Installation

In R:
``` {r}
#if (!requireNamespace("BiocManager", quietly=TRUE))
    #install.packages("BiocManager")
#BiocManager::install("BgeeDB")
```

In case BgeeDB is run on Windows please be sure your OS has curl installed. It is installed by default on Windows 10, version 1803 or later. If Git for Windows is installed on the OS then curl is already installed. If not installed please install it before running BgeeDB in order to avoid potential timeout errors when downloading files. 

## How to use BgeeDB package

### Load the package
``` {r, message = FALSE, warning = FALSE}
library(BgeeDB)
```

### Running example: downloading and formatting processed RNA-seq data

#### List available species in Bgee
The ```listBgeeSpecies()``` function allows to retrieve available species in the Bgee database, and which data types are available for each species. 

``` {r}
listBgeeSpecies()
```

It is possible to list all species from a specific release of Bgee with the ```release``` argument (see ```listBgeeRelease()``` function), and order the species according to a specific columns with the ```ordering``` argument. For example:

``` {r}
listBgeeSpecies(release = "13.2", order = 2)
```

#### Create a new Bgee object

In the following example we will choose to focus on zebrafish ("Danio\_rerio") RNA-seq. Species can be specified using their name or their NCBI taxonomic IDs. To specify that RNA-seq data want to be downloaded, the ```dataType``` argument is set to "rna\_seq". To download Affymetrix microarray data, set this argument to "affymetrix". To download single cell full length RNA-seq data, set this argument to "sc_full_length". To download droplet based single cell RNA-seq data, set this argument to "sc_droplet_based".   

``` {r}
bgee <- Bgee$new(species = "Danio_rerio", dataType = "rna_seq")
```

*Note 1*: It is possible to work with data from a specific release of Bgee by specifying the ```release``` argument, see ```listBgeeRelease()``` function for available releases. 

*Note 2*: The functions of the package will store the downloaded files in a versioned folder created by default in the working directory. These cache files allow faster re-access to the data. The directory where data are stored can be changed with the ```pathToData``` argument.

#### Retrieve the annotation of zebrafish RNA-seq datasets

The ```getAnnotation()``` function will output the list of bulk RNA-seq experiments and libraries available in Bgee for zebrafish. 

``` {r}
annotation_bgee_zebrafish <- getAnnotation(bgee)
# list the first experiments and libraries
lapply(annotation_bgee_zebrafish, head)
```

#### Download the processed zebrafish bulk RNA-seq data

The ```getSampleProcessedData()``` function will download processed bulk RNA-seq data from all zebrafish experiments in Bgee as a data frame.

``` {r}
# download all RNA-seq experiments from zebrafish
data_bgee_zebrafish <- getSampleProcessedData(bgee)
# number of experiments downloaded
length(data_bgee_zebrafish)
# check the downloaded data
lapply(data_bgee_zebrafish, head)
# isolate the first experiment
data_bgee_experiment1 <- data_bgee_zebrafish[[1]]
```

The result of the ```getSampleProcessedData()``` function is a data frame. Each row is a gene and the expression levels are displayed as raw read counts, RPKMs (up to Bgee 13.2), TPMs (from Bgee 14.0), or FPKMs (from Bgee 14.0 to Bgee 14.2). A detection flag indicates if the gene is significantly expressed above background level of expression. From Bgee 15.0 a pValue allows to have a precise metric indicating how much the gene is significantly expressed above background level of expression (the detection flag is still available and a gene is considered as present if pValue < 0.05).

*Note*: If microarray data are downloaded, rows corresponding to probesets and log2 of expression intensities are available instead of read counts/RPKMs/TPMs/FPKMs.

Alternatively, you can choose to download data from one RNA-seq experiment from Bgee with the `experimentId` parameter:

``` {r}
# download data for GSE68666
data_bgee_zebrafish_gse68666 <- getSampleProcessedData(bgee, experimentId = "GSE68666")
```

It is possible to download data by combining filters : 
* experimentId : one or more experimentId,
* sampleId : one or more sampleId (i.e libraryId for RNA-Seq and ChipId for Affymetrix),
* anatEntityId : one or more anatomical entity ID from the UBERON ontology (https://uberon.github.io/),
* withDescendantAnatEntities : filter on the selected anatEntityId and all its descendants (from Bgee 15.0),
* stageId : one or more developmental stage ID from the developmental stage ontologies (https://github.com/obophenotype/developmental-stage-ontologies),
* withDescendantStages : filter on the selected stageId and all its descendants (from Bgee 15.0)
* cellTypeId : one or more cell type, only for single cell datatype (from Bgee 15.0),
* withDescendantCellTypes : filter on the selected cellTypeId and all its descendants (from Bgee 15.0)
* sex : one or more sex (from Bgee 15.0),
* strain : one or more strain (from Bgee 15.0),
``` {r eval=FALSE}
# Examples of data downloading using different filtering combination
# retrieve zebrafish RNA-Seq data for heart (UBERON:0000955) or brain (UBERON:0000948)
data_bgee_zebrafish_filters <- getSampleProcessedData(bgee, anatEntityId = c("UBERON:0000955","UBERON:0000948"))
# retrieve zebrafish RNA-Seq data for embryo (UBERON:0000922) part of the experiment GSE68666
data_bgee_zebrafish_filters <- getSampleProcessedData(bgee, experimentId = "GSE68666", anatEntityId = "UBERON:0000922")
# retrieve zebrafish RNA-Seq data for head kidney (UBERON:0007132) or swim bladder (UBERON:0006860) from post-juvenile adult stage (UBERON:0000113)
data_bgee_zebrafish_filters <- getSampleProcessedData(bgee, stageId = "UBERON:0000113", anatEntityId = c("UBERON:0007132","UBERON:0006860"))
# retrieve zebrafish RNA-Seq data for brain (UBERON:0000948) and all substructures of brain from post-juvenile adult stage (UBERON:0000113)
data_bgee_zebrafish_filters <- getSampleProcessedData(bgee, stageId = "UBERON:0000113", anatEntityId = "UBERON:0000948", withDescendantAnatEntities = TRUE)
```

#### Format the sample level bulk RNA-seq data

It is sometimes easier to work with data organized as a matrix, where rows represent genes or probesets and columns represent different samples. The ```formatData()``` function reformats the data into an ExpressionSet object including:
* An expression data matrix, with genes or probesets as rows, and samples as columns (```assayData``` slot). The ```stats``` argument allows to choose if the matrix should be filled with read counts, RPKMs (up to Bgee 13.2), FPKMs (from Bgee 14.0 to Bgee 14.2), or TPMs (from Bgee 14.0) for RNA-seq data. For microarray data the matrix is filled with log2 expression intensities.
* A data frame listing the samples and their anatomical structure and developmental stage annotation (```phenoData``` slot)
* For microarray data, the mapping from probesets to Ensembl genes (```featureData``` slot)

The ```callType``` argument allows to retain only actively expressed genes or probesets, if set to "present" or "present high quality". Genes or probesets that are absent in a given sample are given ```NA``` values.

```{r}
# use only present calls and fill expression matrix with TPM values
gene.expression.zebrafish.tpm <- formatData(bgee, data_bgee_zebrafish_gse68666, callType = "present", stats = "tpm")
gene.expression.zebrafish.tpm 
```
#### Download cell level single-cell RNA-seq data

The ```getCellProcessedData()``` function will download processed single-cell RNA-seq data at the cell level. In the following example we will focus on downloading UMI counts at cell level for the Gallus gallus experiment ERP132576. The data are stored as h5ad files and can be opened using either the package \"zellkonverter\" (default) or the package \"anndata\" with the attribute ```package = "anndata```. Both \"zellkonverter\" and \"anndata\" packages require a configuration step when used the first time.

``` {r eval=FALSE}
#create a bgee object to download droplet based data from Gallus gallus
bgee <- Bgee$new(species = "Gallus_gallus", dataType = "sc_droplet_based")
# download cell data for one RNA-seq experiment
cell_data_bgee_gallus_gallus <- getCellProcessedData(bgee, experimentId = "ERP132576")
```

The result of the ```getCellProcessedData()``` depends on the package used to load the data. With \"zellkonverter\" it is an object of the class SingleCellExperiment. With \"anndata\" it is an object of the class AnnData. The expression levels are displayed as raw UMI counts.

### Running example: TopAnat gene expression enrichment analysis

For some documentation on the TopAnat analysis, please refer to our publications, or to the web-tool page (<https://www.bgee.org/analysis/top-anat/>).

#### Create a new Bgee object
Similarly to the quantitative data download example above, the first step of a topAnat analysis is to built an object from the Bgee class. For this example, we will focus on zebrafish:

```{r}
# Creating new Bgee class object
bgee <- Bgee$new(species = "Danio_rerio", release = "15.2")
```

*Note* : We are free to specify any data type of interest using the ```dataType``` argument among `rna_seq`, `sc_full_length`, `sc_droplet_based`, `affymetrix`, `est` or `in_situ`, or even a combination of data types. If nothing is specified, as in the above example, all data types available for the targeted species are used. This equivalent to specifying `dataType=c("rna_seq","sc_full_length","sc_droplet_based","affymetrix","est","in_situ")`.

#### Download the data allowing to perform TopAnat analysis

The ```loadTopAnatData()``` function loads a mapping from genes to anatomical structures based on calls of expression in anatomical structures. It also loads the structure of the anatomical ontology and the names of anatomical structures. 

```{r}
# Loading calls of expression
myTopAnatData <- loadTopAnatData(bgee)
# Look at the data
## str(myTopAnatData)
```

The stringency on the quality of expression calls can be changed with the ```confidence``` argument. Finally, if you are interested in expression data coming from a particular developmental stage or a group of stages, please specify the a Uberon stage Id in the ```stage``` argument. 

```{r, eval=FALSE}
## Loading silver and gold expression calls from affymetrix data made on embryonic samples only 
## This is just given as an example, but is not run in this vignette because only few data are returned
bgee <- Bgee$new(species = "Danio_rerio", dataType="affymetrix", release = "15.2")
myTopAnatData <- loadTopAnatData(bgee, stage="UBERON:0000068", confidence="silver")
```

*Note 1*: As mentioned above, the downloaded data files are stored in a versioned folder that can be set with the ```pathToData``` argument when creating the Bgee class object (default is the working directory). If you query again Bgee with the exact same parameters, these cached files will be read instead of querying the web-service again. **It is thus important, if you plan to reuse the same data for multiple parallel topAnat analyses, to plan to make use of these cached files instead of redownloading them for each analysis.** The cached files also give the possibility to repeat analyses offline.

*Note 2*: In releases up to Bgee 13.2 allowed ```confidence`` values were `low_quality` or or `high_quality`. Starting from Bgee 14.0 ```confidence``` values are `gold` or `silver`.

#### Prepare a topAnatData object allowing to perform TopAnat analysis with topGO

First we need to prepare a list of genes in the foreground and in the background. The input format is the same as the gene list required to build a ```topGOdata``` object in the ```topGO``` package: a vector with background genes as names, and 0 or 1 values depending if a gene is in the foreground or not. In this example we will look at genes with an annotated phenotype related to "pectoral fin" . We use the ```biomaRt``` package to retrieve this list of genes. We expect them to be enriched for expression in male tissues, notably the testes. The background list of genes is set to all genes annotated to at least one Gene Ontology term, allowing to account for biases in which types of genes are more likely to receive Gene Ontology annotation.

```{r, eval=FALSE}
# if (!requireNamespace("BiocManager", quietly=TRUE))
    # install.packages("BiocManager")
# BiocManager::install("biomaRt")
library(biomaRt)
ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset="drerio_gene_ensembl", host="mar2016.archive.ensembl.org")

# get the mapping of Ensembl genes to phenotypes. It will corresponds to the background genes
universe <- getBM(filters=c("phenotype_source"), value=c("ZFIN"), attributes=c("ensembl_gene_id","phenotype_description"), mart=ensembl)

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

# Foreground genes are those with an annotated phenotype related to "pectoral fin" 
myGenes <- unique(universe$ensembl_gene_id[universe$phenotype_description %in% phenotypes])

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

# Prepare the topGO object
myTopAnatObject <-  topAnat(myTopAnatData, geneList)
```

The above code using the `biomaRt` package is not executed in this vignette to prevent building issues of our package in case of biomaRt downtime. Instead we use a `geneList` object saved in the `data/` folder that we built using pre-downloaded data. 

```{r}
data(geneList)
myTopAnatObject <-  topAnat(myTopAnatData, geneList)
```

*Warning*: This can be long, especially if the gene list is large, since the Uberon anatomical ontology is large and expression calls will be propagated through the whole ontology (e.g., expression in the forebrain will also be counted as expression in parent structures such as the brain, nervous system, etc). Consider running a script in batch mode if you have multiple analyses to do.

#### Launch the enrichment test

For this step, see the vignette of the ```topGO``` package for more details, as you have to directly use the tests implemented in the ```topGO``` package, as shown in this example:

```{r}
results <- runTest(myTopAnatObject, algorithm = 'weight', statistic = 'fisher')
```

*Warning*: This can be long because of the size of the ontology. Consider running scripts in batch mode if you have multiple analyses to do.

#### Format the table of results after an enrichment test for anatomical terms

The ```makeTable``` function allows to filter and format the test results, and calculate FDR values. 

```{r}
# Display results sigificant at a 1% FDR threshold
tableOver <- makeTable(myTopAnatData, myTopAnatObject, results, cutoff = 0.01)
head(tableOver)
```

At the time of building this vignette (Sept. 2024), there was 9 significant anatomical structures. The first term is “pectoral fin”, and the second “pectoral appendage bud”. Other terms in the list, especially those with high enrichment folds, are clearly related to pectoral fins or substructures of fins. This analysis shows that genes with phenotypic effects on pectoral fins are specifically expressed in or next to these structures

By default results are sorted by p-value, but this can be changed with the ```ordering``` parameter by specifying which column should be used to order the results (preceded by a "-" sign to indicate that ordering should be made in decreasing order). For example, it is often convenient to sort  significant structures by decreasing enrichment fold, using `ordering = -6`. The full table of results can be obtained using `cutoff = 1`.

Of note, it is possible to retrieve for a particular tissue the significant genes that were mapped to it.

```{r}
# In order to retrieve significant genes mapped to the term " paired limb/fin bud"
term <- "UBERON:0004357"
termStat(myTopAnatObject, term) 

# 172 genes mapped to this term for Bgee 15.2
genesInTerm(myTopAnatObject, term)
# 37 significant genes mapped to this term for Bgee 15.2
annotated <- genesInTerm(myTopAnatObject, term)[["UBERON:0004357"]]
annotated[annotated %in% sigGenes(myTopAnatObject)]
```

*Warning*: it is debated whether FDR correction is appropriate on enrichment test results, since tests on different terms of the ontologies are not independent. A nice discussion can be found in the vignette of the ```topGO``` package.

#### Store expression data localy

Since version 2.14.0 (Bioconductor 3.11) BgeeDB store downloaded expression data in a local RSQLite database.
The advantages of this approach compared to the one used in the previous BgeeDB versions are:
* do not anymore need a server with lot of memory to access to subset of huge dataset (e.g GTeX dataset)
* more granular filtering using arguments in the getSampleProcessedData() function
* do not download twice the same data
* fast access to data once integrated in the local database

This approach comes with some drawbacks :
* the SQLite local database use more disk space than the previously conpressed .rds approach
* first access to a dataset takes more time (integration to SQLite local database is time consuming)

It is possible to remove .rds files generated in previous versions of BgeeDB and not used anymore since version
2.14.0. The function below delete all .rds files for the selected species and for all datatype.

```{r eval = FALSE}
bgee <- Bgee$new(species="Mus_musculus", release = "14.1")
# delete all old .rds files of species Mus musculus
deleteOldData(bgee)
```

As the new SQLite approach use more disk space it is now possible to delete all local data of one species from one release of Bgee.

```{r eval = FALSE}
bgee <- Bgee$new(species="Mus_musculus", release = "14.1")
# delete local SQLite database of species Mus musculus from Bgee 14.1
deleteLocalData(bgee)
```