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

Andrea Komljenovic, Julien Roux, Marc Robinson-Rechavi, Frederic Bastian

2016-10-17

BgeeDB is a collection of functions to import data from the Bgee database (http://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 or 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 (http://bgee.org/?page=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 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

Installation

In R:

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

How to use BgeeDB package

Load the package

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.

listBgeeSpecies()
## 
## Querying Bgee to get release information...
## 
## Building URL to query species in Bgee release 13_2...
## 
## Submitting URL to Bgee webservice... (http://r13_2.bgee.org/?page=species&display_type=tsv)
## 
## Query to Bgee webservice successful!
##       ID           GENUS SPECIES_NAME COMMON_NAME AFFYMETRIX   EST IN_SITU
## 1   6239  Caenorhabditis      elegans   c.elegans       TRUE FALSE    TRUE
## 2   7227      Drosophila melanogaster    fruitfly       TRUE  TRUE    TRUE
## 3   7955           Danio        rerio   zebrafish       TRUE  TRUE    TRUE
## 4   8364         Xenopus   tropicalis     xenopus      FALSE  TRUE    TRUE
## 5   9031          Gallus       gallus     chicken      FALSE FALSE   FALSE
## 6   9258 Ornithorhynchus     anatinus    platypus      FALSE FALSE   FALSE
## 7   9544          Macaca      mulatta     macaque      FALSE FALSE   FALSE
## 8   9593         Gorilla      gorilla     gorilla      FALSE FALSE   FALSE
## 9   9597             Pan     paniscus      bonobo      FALSE FALSE   FALSE
## 10  9598             Pan  troglodytes  chimpanzee      FALSE FALSE   FALSE
## 11  9606            Homo      sapiens       human       TRUE  TRUE   FALSE
## 12  9823             Sus       scrofa         pig      FALSE FALSE   FALSE
## 13  9913             Bos       taurus         cow      FALSE FALSE   FALSE
## 14 10090             Mus     musculus       mouse       TRUE  TRUE    TRUE
## 15 10116          Rattus   norvegicus         rat      FALSE FALSE   FALSE
## 16 13616     Monodelphis    domestica     opossum      FALSE FALSE   FALSE
## 17 28377          Anolis carolinensis      anolis      FALSE FALSE   FALSE
##    RNA_SEQ
## 1     TRUE
## 2    FALSE
## 3    FALSE
## 4     TRUE
## 5     TRUE
## 6     TRUE
## 7     TRUE
## 8     TRUE
## 9     TRUE
## 10    TRUE
## 11    TRUE
## 12    TRUE
## 13    TRUE
## 14    TRUE
## 15    TRUE
## 16    TRUE
## 17    TRUE

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:

listBgeeSpecies(release = "13.2", order = 2)
## 
## Querying Bgee to get release information...
## 
## Building URL to query species in Bgee release 13_2...
## 
## Submitting URL to Bgee webservice... (http://r13_2.bgee.org/?page=species&display_type=tsv)
## 
## Query to Bgee webservice successful!
##       ID           GENUS SPECIES_NAME COMMON_NAME AFFYMETRIX   EST IN_SITU
## 17 28377          Anolis carolinensis      anolis      FALSE FALSE   FALSE
## 13  9913             Bos       taurus         cow      FALSE FALSE   FALSE
## 1   6239  Caenorhabditis      elegans   c.elegans       TRUE FALSE    TRUE
## 3   7955           Danio        rerio   zebrafish       TRUE  TRUE    TRUE
## 2   7227      Drosophila melanogaster    fruitfly       TRUE  TRUE    TRUE
## 5   9031          Gallus       gallus     chicken      FALSE FALSE   FALSE
## 8   9593         Gorilla      gorilla     gorilla      FALSE FALSE   FALSE
## 11  9606            Homo      sapiens       human       TRUE  TRUE   FALSE
## 7   9544          Macaca      mulatta     macaque      FALSE FALSE   FALSE
## 16 13616     Monodelphis    domestica     opossum      FALSE FALSE   FALSE
## 14 10090             Mus     musculus       mouse       TRUE  TRUE    TRUE
## 6   9258 Ornithorhynchus     anatinus    platypus      FALSE FALSE   FALSE
## 9   9597             Pan     paniscus      bonobo      FALSE FALSE   FALSE
## 10  9598             Pan  troglodytes  chimpanzee      FALSE FALSE   FALSE
## 15 10116          Rattus   norvegicus         rat      FALSE FALSE   FALSE
## 12  9823             Sus       scrofa         pig      FALSE FALSE   FALSE
## 4   8364         Xenopus   tropicalis     xenopus      FALSE  TRUE    TRUE
##    RNA_SEQ
## 17    TRUE
## 13    TRUE
## 1     TRUE
## 3    FALSE
## 2    FALSE
## 5     TRUE
## 8     TRUE
## 11    TRUE
## 7     TRUE
## 16    TRUE
## 14    TRUE
## 6     TRUE
## 9     TRUE
## 10    TRUE
## 15    TRUE
## 12    TRUE
## 4     TRUE

Create a new Bgee object

In the following example we will choose to focus on mouse (“Mus_musculus”) 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”.

bgee <- Bgee$new(species = "Mus_musculus", dataType = "rna_seq")
## 
## Querying Bgee to get release information...
## 
## Building URL to query species in Bgee release 13_2...
## 
## Submitting URL to Bgee webservice... (http://r13_2.bgee.org/?page=species&display_type=tsv)
## 
## Query to Bgee webservice successful!
## 
## API key built: 5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad7

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

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 mouse RNA-seq datasets

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

annotation_bgee_mouse <- getAnnotation(bgee)
## 
## Downloading annotation files...
## 
## Saved annotation files in /tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Mus_musculus_Bgee_13_2 folder.
# list the first experiments and libraries
lapply(annotation_bgee_mouse, head)
## $sample.annotation
##   Experiment.ID Library.ID Library.secondary.ID Anatomical.entity.ID
## 1      GSE30617  GSM759583            ERX012363       UBERON:0000948
## 2      GSE30617  GSM759584            ERX012348       UBERON:0000948
## 3      GSE30617  GSM759585            ERX012344       UBERON:0000948
## 4      GSE30617  GSM759586            ERX012362       UBERON:0000948
## 5      GSE30617  GSM759587            ERX012378       UBERON:0000948
## 6      GSE30617  GSM759588            ERX012374       UBERON:0000948
##   Anatomical.entity.name       Stage.ID      Stage.name
## 1                  heart MmusDv:0000052 8 weeks (mouse)
## 2                  heart MmusDv:0000052 8 weeks (mouse)
## 3                  heart MmusDv:0000052 8 weeks (mouse)
## 4                  heart MmusDv:0000052 8 weeks (mouse)
## 5                  heart MmusDv:0000052 8 weeks (mouse)
## 6                  heart MmusDv:0000052 8 weeks (mouse)
##                   Platform.ID Library.type Library.orientation Read.count
## 1 Illumina Genome Analyzer II       paired          unstranded   31000737
## 2 Illumina Genome Analyzer II       paired          unstranded    8605668
## 3 Illumina Genome Analyzer II       paired          unstranded   30075234
## 4 Illumina Genome Analyzer II       paired          unstranded   29498377
## 5 Illumina Genome Analyzer II       paired          unstranded   27824366
## 6 Illumina Genome Analyzer II       paired          unstranded   31160789
##   Left.part.mapped.read.count.Mapped.read.count
## 1                                      23674754
## 2                                       7099912
## 3                                      24308959
## 4                                      24792064
## 5                                      19632932
## 6                                      21342634
##   Right.part.mapped.read.count Min..read.length Max..read.length
## 1                     23809225               76               76
## 2                      7083684               76               76
## 3                     23884114               76               76
## 4                     24122931               76               76
## 5                     20578789               76               76
## 6                     23529915               76               76
##   All.genes.percent.present Protein.coding.genes.percent.present
## 1                     50.73                                31.84
## 2                     41.94                                25.69
## 3                     50.70                                31.66
## 4                     49.12                                30.51
## 5                     50.24                                31.51
## 6                     50.32                                31.72
##   Intergenic.regions.percent.present   Run.IDs Discarded.run.IDs
## 1                               0.98 ERR032227                NA
## 2                               1.17 ERR032228                NA
## 3                               0.95 ERR032229                NA
## 4                               0.87 ERR032238                NA
## 5                               0.97 ERR032230                NA
## 6                               1.04 ERR032231                NA
##   Data.source                                             Data.source.URL
## 1         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759583
## 2         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759584
## 3         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759585
## 4         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759586
## 5         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759587
## 6         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM759588
##                                                                                                                                             Bgee.normalized.data.URL
## 1 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 2 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 3 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 4 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 5 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 6 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
##                                                Raw.file.URL
## 1 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759583
## 2 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759584
## 3 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759585
## 4 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759586
## 5 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759587
## 6 http://trace.ncbi.nlm.nih.gov/Traces/study/?acc=GSM759588
## 
## $experiment.annotation
##   Experiment.ID Library.count Condition.count Organ.count Stage.count
## 1      GSE30617            36               6           6           1
## 2      GSE41637            26               9           9           1
## 3      GSE30352            17               6           6           1
## 4      GSE36026            12              12          12           3
## 5      GSE43520             9               5           4           2
## 6      GSE41338             6               5           5           1
##   Data.source                                            Data.source.URL
## 1         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30617
## 2         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41637
## 3         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30352
## 4         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36026
## 5         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE43520
## 6         GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41338
##                                                                                                                                             Bgee.normalized.data.URL
## 1 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip
## 2 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE41637.tsv.zip
## 3 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30352.tsv.zip
## 4 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE36026.tsv.zip
## 5 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE43520.tsv.zip
## 6 ftp://lausanne.isb-sib.ch/pub/databases/Bgee/bgee_v13_1/download/processed_expr_values/rna_seq/Mus_musculus/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE41338.tsv.zip
##                                                             Experiment.name
## 1                                          [E-MTAB-599] Mouse Transcriptome
## 2 Evolutionary dynamics of gene and isoform regulation in mammalian tissues
## 3               The evolution of gene expression levels in mammalian organs
## 4                                                  RNA-seq from ENCODE/LICR
## 5  The evolution of lncRNA repertoires and expression patterns in tetrapods
## 6  The evolutionary landscape of alternative splicing in vertebrate species
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                Experiment.description
## 1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 Sequencing the transcriptome of DBAxC57BL/6J mice. To study the regulation of transcription, splicing and RNA turnover we have sequenced the transcriptomes of tissues collected DBAxC57BL/6J mice.
## 2                                                                                                                                                                                                                                                                                                                                                      Most mammalian genes produce multiple distinct mRNAs through alternative splicing, but the extent of splicing conservation is not clear. To assess tissue-specific transcriptome variation across mammals, we sequenced cDNA from 9 tissues from 4 mammals and one bird in biological triplicate, at unprecedented depth. We find that while tissue-specific gene expression programs are largely conserved, alternative splicing is well conserved in only a subset of tissues and is frequently lineage-specific. Thousands of novel, lineage-specific and conserved alternative exons were identified; widely conserved alternative exons had signatures of binding by MBNL, PTB, RBFOX, STAR and TIA family splicing factors, implicating them as ancestral mammalian splicing regulators. Our data also indicates that alternative splicing is often used to alter protein phosphorylatability, delimiting the scope of kinase signaling.
## 3 Changes in gene expression are thought to underlie many of the phenotypic differences between species. However, large-scale analyses of gene expression evolution were until recently prevented by technological limitations. Here we report the sequencing of polyadenylated RNA from six organs across ten species that represent all major mammalian lineages (placentals, marsupials and monotremes) and birds (the evolutionary outgroup), with the goal of understanding the dynamics of mammalian transcriptome evolution. We show that the rate of gene expression evolution varies among organs, lineages and chromosomes, owing to differences in selective pressures: transcriptome change was slow in nervous tissues and rapid in testes, slower in rodents than in apes and monotremes, and rapid for the X chromosome right after its formation. Although gene expression evolution in mammals was strongly shaped by purifying selection, we identify numerous potentially selectively driven expression switches, which occurred at different rates across lineages and tissues and which probably contributed to the specific organ biology of various mammals. Our transcriptome data provide a valuable resource for functional and evolutionary analyses of mammalian genomes.
## 4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       Using RNA-Seq (Mortazavi et al., 2008), high-resolution genome-wide maps of the mouse transcriptome across multiple mouse (C57Bl/6) tissues and primary cells were generated.
## 5                                                                                                                                                                                                                                                                                                                                                                                                                                To broaden our understanding of lncRNA evolution, we used an extensive RNA-seq dataset to establish lncRNA repertoires and homologous gene families in 11 tetrapod species. We analyzed the poly- adenylated transcriptomes of 8 organs (cortex/whole brain without cerebellum, cerebellum, heart, kidney, liver, placenta, ovary and testis) and 11 species (human, chimpanzee, bonobo, gorilla, orangutan, macaque, mouse, opossum, platypus, chicken and the frog Xenopus tropicalis), which shared a common ancestor ~370 millions of years (MY) ago. Our dataset included 47 strand-specific samples, which allowed us to confirm the orientation of gene predictions and to address the evolution of sense-antisense transcripts. See also GSE43721 (Soumillon et al, Cell Reports, 2013) for three strand-specific samples for mouse brain, liver and testis.
## 6                                                                                       How species with similar repertoires of protein coding genes differ so dramatically at the phenotypic level is poorly understood. From comparing the transcriptomes of multiple organs from vertebrate species spanning ~350 million years of evolution, we observe significant differences in alternative splicing complexity between the main vertebrate lineages, with the highest complexity in the primate lineage. Moreover, within as little as six million years, the splicing profiles of physiologically-equivalent organs have diverged to the extent that they are more strongly related to the identity of a species than they are to organ type. Most vertebrate species-specific splicing patterns are governed by the highly variable use of a largely conserved cis-regulatory code. However, a smaller number of pronounced species-dependent splicing changes are predicted to remodel interactions involving factors acting at multiple steps in gene regulation. These events are expected to further contribute to the dramatic diversification of alternative splicing as well as to other gene regulatory changes that contribute to phenotypic differences among vertebrate species.

Download the processed mouse RNA-seq data

The getData() function will download processed RNA-seq data from all mouse experiments in Bgee as a list.

# download all RNA-seq experiments from mouse
data_bgee_mouse <- getData(bgee)
## 
## The experiment is not defined. Hence taking all rna_seq experiments available for Mus_musculus.
## 
## Downloading expression data...
## 
## Saved expression data file in /tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Mus_musculus_Bgee_13_2 folder. Now unzipping file...
## 
## Saving all data in .rds file...
# number of experiments downloaded
length(data_bgee_mouse)
## [1] 7
# check the downloaded data
lapply(data_bgee_mouse, head)
## [[1]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE36026  GSM929703       single ENSMUSG00000000001
## 2      GSE36026  GSM929703       single ENSMUSG00000000003
## 3      GSE36026  GSM929703       single ENSMUSG00000000028
## 4      GSE36026  GSM929703       single ENSMUSG00000000031
## 5      GSE36026  GSM929703       single ENSMUSG00000000037
## 6      GSE36026  GSM929703       single ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0001348   brown adipose tissue MmusDv:0000074
## 2       UBERON:0001348   brown adipose tissue MmusDv:0000074
## 3       UBERON:0001348   brown adipose tissue MmusDv:0000074
## 4       UBERON:0001348   brown adipose tissue MmusDv:0000074
## 5       UBERON:0001348   brown adipose tissue MmusDv:0000074
## 6       UBERON:0001348   brown adipose tissue MmusDv:0000074
##         Stage.name Read.count    RPKM Detection.flag Detection.quality
## 1 24 weeks (mouse)        412 5.49619        present      high quality
## 2 24 weeks (mouse)          0 0.00000         absent      high quality
## 3 24 weeks (mouse)         53 1.46527        present      high quality
## 4 24 weeks (mouse)         94 3.22440        present      high quality
## 5 24 weeks (mouse)          4 0.05278         absent      high quality
## 6 24 weeks (mouse)         10 0.65387         absent      high quality
##                            State.in.Bgee
## 1                Used in expression call
## 2 Result excluded, reason: pre-filtering
## 3                Used in expression call
## 4                Used in expression call
## 5             Used in no-expression call
## 6             Used in no-expression call
## 
## [[2]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE41637 GSM1020640       paired ENSMUSG00000000001
## 2      GSE41637 GSM1020640       paired ENSMUSG00000000003
## 3      GSE41637 GSM1020640       paired ENSMUSG00000000028
## 4      GSE41637 GSM1020640       paired ENSMUSG00000000031
## 5      GSE41637 GSM1020640       paired ENSMUSG00000000037
## 6      GSE41637 GSM1020640       paired ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000955                  brain UBERON:0000113
## 2       UBERON:0000955                  brain UBERON:0000113
## 3       UBERON:0000955                  brain UBERON:0000113
## 4       UBERON:0000955                  brain UBERON:0000113
## 5       UBERON:0000955                  brain UBERON:0000113
## 6       UBERON:0000955                  brain UBERON:0000113
##                  Stage.name Read.count     RPKM Detection.flag
## 1 post-juvenile adult stage       2677 12.41761        present
## 2 post-juvenile adult stage          0  0.00000         absent
## 3 post-juvenile adult stage         87  0.83635         absent
## 4 post-juvenile adult stage         27  0.32204         absent
## 5 post-juvenile adult stage        140  0.64249         absent
## 6 post-juvenile adult stage         15  0.34104         absent
##   Detection.quality                                  State.in.Bgee
## 1      high quality                        Used in expression call
## 2      high quality         Result excluded, reason: pre-filtering
## 3      high quality                        Used in expression call
## 4      high quality Result excluded, reason: noExpression conflict
## 5      high quality                     Used in no-expression call
## 6      high quality                        Used in expression call
## 
## [[3]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE30617  GSM759583       paired ENSMUSG00000000001
## 2      GSE30617  GSM759583       paired ENSMUSG00000000003
## 3      GSE30617  GSM759583       paired ENSMUSG00000000028
## 4      GSE30617  GSM759583       paired ENSMUSG00000000031
## 5      GSE30617  GSM759583       paired ENSMUSG00000000037
## 6      GSE30617  GSM759583       paired ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000948                  heart MmusDv:0000052
## 2       UBERON:0000948                  heart MmusDv:0000052
## 3       UBERON:0000948                  heart MmusDv:0000052
## 4       UBERON:0000948                  heart MmusDv:0000052
## 5       UBERON:0000948                  heart MmusDv:0000052
## 6       UBERON:0000948                  heart MmusDv:0000052
##        Stage.name Read.count     RPKM Detection.flag Detection.quality
## 1 8 weeks (mouse)        941 11.57664        present      high quality
## 2 8 weeks (mouse)          0  0.00000         absent      high quality
## 3 8 weeks (mouse)         72  1.83570        present      high quality
## 4 8 weeks (mouse)        225  7.11763        present      high quality
## 5 8 weeks (mouse)         15  0.18255         absent      high quality
## 6 8 weeks (mouse)          2  0.12059         absent      high quality
##                            State.in.Bgee
## 1                Used in expression call
## 2 Result excluded, reason: pre-filtering
## 3                Used in expression call
## 4                Used in expression call
## 5             Used in no-expression call
## 6                Used in expression call
## 
## [[4]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE41338 GSM1015150       paired ENSMUSG00000000001
## 2      GSE41338 GSM1015150       paired ENSMUSG00000000003
## 3      GSE41338 GSM1015150       paired ENSMUSG00000000028
## 4      GSE41338 GSM1015150       paired ENSMUSG00000000031
## 5      GSE41338 GSM1015150       paired ENSMUSG00000000037
## 6      GSE41338 GSM1015150       paired ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000955                  brain UBERON:0000113
## 2       UBERON:0000955                  brain UBERON:0000113
## 3       UBERON:0000955                  brain UBERON:0000113
## 4       UBERON:0000955                  brain UBERON:0000113
## 5       UBERON:0000955                  brain UBERON:0000113
## 6       UBERON:0000955                  brain UBERON:0000113
##                  Stage.name Read.count    RPKM Detection.flag
## 1 post-juvenile adult stage        722 4.17829        present
## 2 post-juvenile adult stage          0 0.00000         absent
## 3 post-juvenile adult stage         59 0.70761         absent
## 4 post-juvenile adult stage         26 0.38690         absent
## 5 post-juvenile adult stage         17 0.09733         absent
## 6 post-juvenile adult stage         46 1.30483        present
##   Detection.quality                                  State.in.Bgee
## 1      high quality                        Used in expression call
## 2      high quality         Result excluded, reason: pre-filtering
## 3      high quality                        Used in expression call
## 4      high quality Result excluded, reason: noExpression conflict
## 5      high quality                     Used in no-expression call
## 6      high quality                        Used in expression call
## 
## [[5]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE30352  GSM752614       single ENSMUSG00000000001
## 2      GSE30352  GSM752614       single ENSMUSG00000000003
## 3      GSE30352  GSM752614       single ENSMUSG00000000028
## 4      GSE30352  GSM752614       single ENSMUSG00000000031
## 5      GSE30352  GSM752614       single ENSMUSG00000000037
## 6      GSE30352  GSM752614       single ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000955                  brain UBERON:0000113
## 2       UBERON:0000955                  brain UBERON:0000113
## 3       UBERON:0000955                  brain UBERON:0000113
## 4       UBERON:0000955                  brain UBERON:0000113
## 5       UBERON:0000955                  brain UBERON:0000113
## 6       UBERON:0000955                  brain UBERON:0000113
##                  Stage.name Read.count     RPKM Detection.flag
## 1 post-juvenile adult stage        550 11.92921        present
## 2 post-juvenile adult stage          0  0.00000         absent
## 3 post-juvenile adult stage         12  0.53941         absent
## 4 post-juvenile adult stage          2  0.11157         absent
## 5 post-juvenile adult stage         13  0.27897         absent
## 6 post-juvenile adult stage          7  0.74416         absent
##   Detection.quality                                  State.in.Bgee
## 1      high quality                        Used in expression call
## 2      high quality         Result excluded, reason: pre-filtering
## 3      high quality                        Used in expression call
## 4      high quality Result excluded, reason: noExpression conflict
## 5      high quality                     Used in no-expression call
## 6      high quality                        Used in expression call
## 
## [[6]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE43721 GSM1069683       single ENSMUSG00000000001
## 2      GSE43721 GSM1069683       single ENSMUSG00000000003
## 3      GSE43721 GSM1069683       single ENSMUSG00000000028
## 4      GSE43721 GSM1069683       single ENSMUSG00000000031
## 5      GSE43721 GSM1069683       single ENSMUSG00000000037
## 6      GSE43721 GSM1069683       single ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000955                  brain MmusDv:0000055
## 2       UBERON:0000955                  brain MmusDv:0000055
## 3       UBERON:0000955                  brain MmusDv:0000055
## 4       UBERON:0000955                  brain MmusDv:0000055
## 5       UBERON:0000955                  brain MmusDv:0000055
## 6       UBERON:0000955                  brain MmusDv:0000055
##         Stage.name Read.count    RPKM Detection.flag Detection.quality
## 1 11 weeks (mouse)       1024 9.36418        present      high quality
## 2 11 weeks (mouse)          0 0.00000         absent      high quality
## 3 11 weeks (mouse)         61 1.15606        present      high quality
## 4 11 weeks (mouse)         19 0.44676         absent      high quality
## 5 11 weeks (mouse)         45 0.40712         absent      high quality
## 6 11 weeks (mouse)         35 1.56882        present      high quality
##                                    State.in.Bgee
## 1                        Used in expression call
## 2         Result excluded, reason: pre-filtering
## 3                        Used in expression call
## 4 Result excluded, reason: noExpression conflict
## 5                     Used in no-expression call
## 6                        Used in expression call
## 
## [[7]]
##   Experiment.ID Library.ID Library.type            Gene.ID
## 1      GSE43520 GSM1064841       single ENSMUSG00000000001
## 2      GSE43520 GSM1064841       single ENSMUSG00000000003
## 3      GSE43520 GSM1064841       single ENSMUSG00000000028
## 4      GSE43520 GSM1064841       single ENSMUSG00000000031
## 5      GSE43520 GSM1064841       single ENSMUSG00000000037
## 6      GSE43520 GSM1064841       single ENSMUSG00000000049
##   Anatomical.entity.ID Anatomical.entity.name       Stage.ID
## 1       UBERON:0000955                  brain UBERON:0000113
## 2       UBERON:0000955                  brain UBERON:0000113
## 3       UBERON:0000955                  brain UBERON:0000113
## 4       UBERON:0000955                  brain UBERON:0000113
## 5       UBERON:0000955                  brain UBERON:0000113
## 6       UBERON:0000955                  brain UBERON:0000113
##                  Stage.name Read.count    RPKM Detection.flag
## 1 post-juvenile adult stage        852 9.45702        present
## 2 post-juvenile adult stage          0 0.00000         absent
## 3 post-juvenile adult stage         36 0.82814         absent
## 4 post-juvenile adult stage          3 0.08563         absent
## 5 post-juvenile adult stage         31 0.34042         absent
## 6 post-juvenile adult stage          3 0.16322         absent
##   Detection.quality                                  State.in.Bgee
## 1      high quality                        Used in expression call
## 2      high quality         Result excluded, reason: pre-filtering
## 3      high quality                        Used in expression call
## 4      high quality Result excluded, reason: noExpression conflict
## 5      high quality                     Used in no-expression call
## 6      high quality                        Used in expression call
# isolate the first experiment
data_bgee_experiment1 <- data_bgee_mouse[[1]]

The result of the getData() function is, for each experiment, a data frame with the different samples listed in rows, one after the other. Each row is a gene and the expression levels are displayed as raw read counts or RPKMs. A detection flag indicates is the gene is significantly expressed above background level of expression.

Note 1: An additional column in the data frame including expression in the TPM unit will be available from Bgee release 14 (planned for the end of 2016).

Note 2: If microarray data are downloaded, rows correspond to probesets and log2 of expression intensities are available instead of read counts/RPKMs.

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

# download data for GSE30617
data_bgee_mouse_gse30617 <- getData(bgee, experimentId = "GSE30617")
## 
## Downloading expression data for the experiment GSE30617 ...
## 
## Saved expression data file in /tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Mus_musculus_Bgee_13_2 folder. Now unzipping /tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Mus_musculus_Bgee_13_2/Mus_musculus_RNA-Seq_read_counts_RPKM_GSE30617.tsv.zip file...
## 
## Saving all data in .rds file...

Format the 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 (and soon TPMs) 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.

# use only present calls and fill expression matric with RPKM values
gene.expression.mouse.rpkm <- formatData(bgee, data_bgee_mouse_gse30617, callType = "present", stats = "rpkm")
## 
## Extracting expression data matrix...
##   Keeping only present genes.
## 
## Extracting features information...
## 
## Extracting samples information...
gene.expression.mouse.rpkm 
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 39179 features, 36 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: GSM759583 GSM759584 ... GSM759618 (36 total)
##   varLabels: Library.ID Anatomical.entity.ID ... Stage.name (5
##     total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ENSMUSG00000000001 ENSMUSG00000000003 ...
##     ENSMUSG00000099334 (39179 total)
##   fvarLabels: Gene.ID
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

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 (http://bgee.org/?page=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:

# Creating new Bgee class object
bgee <- Bgee$new(species = "Danio_rerio")
## 
## NOTE: You did not specify any data type. The argument dataType will be set to c("rna_seq","affymetrix","est","in_situ") for the next steps.
## 
## Querying Bgee to get release information...
## 
## NOTE: the file describing Bgee species information for release 13_2 was found in the download directory /tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes. Data will not be redownloaded.
## 
## API key built: 5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad7

Note : We are free to specify any data type of interest using the dataType argument among rna_seq, 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","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.

# Loading calls of expression
myTopAnatData <- loadTopAnatData(bgee)
## 
## Building URLs to retrieve organ relationships from Bgee.........
##    URL successfully built (http://r13_2.bgee.org/?page=dao&action=org.bgee.model.dao.api.ontologycommon.RelationDAO.getAnatEntityRelations&display_type=tsv&species_list=7955&attr_list=SOURCE_ID&attr_list=TARGET_ID&api_key=5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad7&source=BgeeDB_R_package&source_version=2.0.0)
##    Submitting URL to Bgee webservice (can be long)
##    Got results from Bgee webservice. Files are written in "/tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Danio_rerio_Bgee_13_2"
## 
## Building URLs to retrieve organ names from Bgee.................
##    URL successfully built (http://r13_2.bgee.org/?page=dao&action=org.bgee.model.dao.api.anatdev.AnatEntityDAO.getAnatEntities&display_type=tsv&species_list=7955&attr_list=ID&attr_list=NAME&api_key=5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad7&source=BgeeDB_R_package&source_version=2.0.0)
##    Submitting URL to Bgee webservice (can be long)
##    Got results from Bgee webservice. Files are written in "/tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Danio_rerio_Bgee_13_2"
## 
## Building URLs to retrieve mapping of gene to organs from Bgee...
##    URL successfully built (http://r13_2.bgee.org/?page=dao&action=org.bgee.model.dao.api.expressiondata.ExpressionCallDAO.getExpressionCalls&display_type=tsv&species_list=7955&attr_list=GENE_ID&attr_list=ANAT_ENTITY_ID&api_key=5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad7&source=BgeeDB_R_package&source_version=2.0.0)
##    Submitting URL to Bgee webservice (can be long)
##    Got results from Bgee webservice. Files are written in "/tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Danio_rerio_Bgee_13_2"
## 
## Parsing the results.............................................
## 
## Adding BGEE:0 as unique root of all terms of the ontology.......
## 
## Done.
# Look at the data
str(myTopAnatData)
## List of 4
##  $ gene2anatomy       :List of 18715
##   ..$ ENSDARG00000000001: chr [1:3] "UBERON:0000468" "UBERON:0001997" "ZFA:0001093"
##   ..$ ENSDARG00000000002: chr [1:11] "UBERON:0000019" "UBERON:0000468" "UBERON:0000965" "UBERON:0002329" ...
##   ..$ ENSDARG00000000018: chr [1:28] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000019: chr [1:23] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000068: chr [1:53] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000069: chr [1:23] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000102: chr [1:7] "UBERON:0000151" "UBERON:0000468" "UBERON:0001774" "UBERON:0002329" ...
##   ..$ ENSDARG00000000103: chr [1:20] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000142: chr [1:4] "UBERON:0000080" "UBERON:0000468" "UBERON:0000992" "UBERON:0001555"
##   ..$ ENSDARG00000000151: chr [1:13] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000161: chr "UBERON:0000955"
##   ..$ ENSDARG00000000175: chr [1:36] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000183: chr [1:10] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000189: chr [1:6] "UBERON:0000965" "UBERON:0001891" "UBERON:0001893" "UBERON:0001894" ...
##   ..$ ENSDARG00000000212: chr [1:20] "UBERON:0000004" "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" ...
##   ..$ ENSDARG00000000229: chr [1:4] "UBERON:0000307" "UBERON:0000468" "UBERON:0000992" "UBERON:0007010"
##   ..$ ENSDARG00000000241: chr [1:22] "UBERON:0000080" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000324: chr "UBERON:0000473"
##   ..$ ENSDARG00000000369: chr [1:2] "UBERON:0000080" "UBERON:0000468"
##   ..$ ENSDARG00000000370: chr "UBERON:0000468"
##   ..$ ENSDARG00000000380: chr [1:9] "UBERON:0000019" "UBERON:0000468" "UBERON:0000922" "UBERON:0000955" ...
##   ..$ ENSDARG00000000384: chr [1:25] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000394: chr [1:17] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000423: chr [1:3] "UBERON:0000468" "UBERON:0000473" "UBERON:0000992"
##   ..$ ENSDARG00000000442: chr [1:20] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000460: chr "UBERON:0000948"
##   ..$ ENSDARG00000000472: chr [1:24] "UBERON:0000468" "UBERON:0000955" "UBERON:0001675" "UBERON:0001700" ...
##   ..$ ENSDARG00000000476: chr [1:16] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000000489: chr [1:18] "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" "UBERON:0000922" ...
##   ..$ ENSDARG00000000503: chr [1:14] "UBERON:0000019" "UBERON:0000468" "UBERON:0000473" "UBERON:0000922" ...
##   ..$ ENSDARG00000000516: chr [1:3] "UBERON:0000468" "UBERON:0000473" "UBERON:0004734"
##   ..$ ENSDARG00000000529: chr [1:22] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000000540: chr [1:20] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000551: chr [1:22] "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" "UBERON:0000955" ...
##   ..$ ENSDARG00000000563: chr [1:31] "UBERON:0000019" "UBERON:0000080" "UBERON:0000151" "UBERON:0000468" ...
##   ..$ ENSDARG00000000567: chr [1:2] "UBERON:0000019" "UBERON:0000468"
##   ..$ ENSDARG00000000568: chr [1:20] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000588: chr [1:2] "UBERON:0000468" "UBERON:0001997"
##   ..$ ENSDARG00000000606: chr [1:2] "UBERON:0000019" "UBERON:0000468"
##   ..$ ENSDARG00000000631: chr [1:21] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000638: chr [1:3] "UBERON:0000019" "UBERON:0000468" "UBERON:0000966"
##   ..$ ENSDARG00000000651: chr [1:14] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000656: chr [1:12] "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000000666: chr [1:6] "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000000690: chr [1:24] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000729: chr [1:19] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000730: chr [1:4] "UBERON:0000468" "UBERON:0002107" "UBERON:0008897" "ZFA:0001093"
##   ..$ ENSDARG00000000747: chr [1:3] "UBERON:0000019" "UBERON:0000468" "UBERON:0001555"
##   ..$ ENSDARG00000000760: chr [1:20] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000000767: chr [1:30] "UBERON:0000019" "UBERON:0000033" "UBERON:0000080" "UBERON:0000178" ...
##   ..$ ENSDARG00000000779: chr [1:24] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000796: chr [1:31] "UBERON:0000004" "UBERON:0000019" "UBERON:0000080" "UBERON:0000151" ...
##   ..$ ENSDARG00000000804: chr [1:3] "UBERON:0000468" "UBERON:0001997" "ZFA:0001093"
##   ..$ ENSDARG00000000830: chr [1:5] "UBERON:0000080" "UBERON:0000468" "UBERON:0000992" "UBERON:0002107" ...
##   ..$ ENSDARG00000000837: chr [1:3] "UBERON:0000178" "UBERON:0000468" "UBERON:0001997"
##   ..$ ENSDARG00000000853: chr [1:22] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000000857: chr [1:34] "UBERON:0000004" "UBERON:0000080" "UBERON:0000151" "UBERON:0000307" ...
##   ..$ ENSDARG00000000861: chr [1:13] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000000887: chr [1:29] "UBERON:0000019" "UBERON:0000080" "UBERON:0000151" "UBERON:0000307" ...
##   ..$ ENSDARG00000000906: chr [1:6] "UBERON:0000080" "UBERON:0000468" "UBERON:0001555" "UBERON:0001997" ...
##   ..$ ENSDARG00000000935: chr [1:18] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000943: chr [1:3] "UBERON:0000080" "UBERON:0000468" "ZFA:0001093"
##   ..$ ENSDARG00000000966: chr [1:26] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000000972: chr [1:23] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001014: chr [1:4] "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" "UBERON:0008897"
##   ..$ ENSDARG00000001015: chr [1:3] "UBERON:0000019" "UBERON:0000468" "UBERON:0000955"
##   ..$ ENSDARG00000001057: chr [1:33] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000001127: chr [1:5] "CL:0005000" "UBERON:0000468" "UBERON:0000955" "UBERON:0001898" ...
##   ..$ ENSDARG00000001129: chr [1:20] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000001162: chr [1:17] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000001169: chr [1:14] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000001210: chr [1:7] "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000001220: chr [1:23] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001234: chr [1:3] "UBERON:0000468" "UBERON:0000955" "UBERON:0000966"
##   ..$ ENSDARG00000001241: chr [1:27] "UBERON:0000019" "UBERON:0000080" "UBERON:0000151" "UBERON:0000307" ...
##   ..$ ENSDARG00000001244: chr [1:22] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001249: chr [1:4] "UBERON:0000080" "UBERON:0000468" "UBERON:0000473" "UBERON:0000992"
##   ..$ ENSDARG00000001285: chr [1:2] "UBERON:0000468" "UBERON:0001997"
##   ..$ ENSDARG00000001303: chr [1:16] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001313: chr [1:25] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000473" ...
##   ..$ ENSDARG00000001354: chr "UBERON:0000468"
##   ..$ ENSDARG00000001414: chr [1:7] "UBERON:0000468" "UBERON:0000473" "UBERON:0000922" "UBERON:0001997" ...
##   ..$ ENSDARG00000001431: chr [1:28] "UBERON:0000019" "UBERON:0000080" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000001437: chr [1:8] "UBERON:0000033" "UBERON:0000468" "UBERON:0000966" "UBERON:0001017" ...
##   ..$ ENSDARG00000001442: chr [1:13] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000001447: chr [1:3] "UBERON:0000307" "UBERON:0000468" "UBERON:0007010"
##   ..$ ENSDARG00000001452: chr [1:13] "UBERON:0000080" "UBERON:0000468" "UBERON:0000473" "UBERON:0000926" ...
##   ..$ ENSDARG00000001463: chr [1:4] "UBERON:0000019" "UBERON:0000468" "UBERON:0000966" "UBERON:0001997"
##   ..$ ENSDARG00000001498: chr [1:4] "UBERON:0000080" "UBERON:0000468" "UBERON:0000992" "UBERON:0001846"
##   ..$ ENSDARG00000001549: chr [1:13] "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000001557: chr [1:25] "UBERON:0000019" "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" ...
##   ..$ ENSDARG00000001558: chr [1:22] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001559: chr "ZFA:0001093"
##   ..$ ENSDARG00000001572: chr "UBERON:0000080"
##   ..$ ENSDARG00000001578: chr [1:28] "UBERON:0000019" "UBERON:0000080" "UBERON:0000307" "UBERON:0000468" ...
##   ..$ ENSDARG00000001621: chr "UBERON:0000468"
##   ..$ ENSDARG00000001634: chr [1:13] "UBERON:0000019" "UBERON:0000080" "UBERON:0000151" "UBERON:0000468" ...
##   ..$ ENSDARG00000001637: chr [1:16] "UBERON:0000307" "UBERON:0000468" "UBERON:0000922" "UBERON:0000948" ...
##   ..$ ENSDARG00000001656: chr [1:13] "UBERON:0000019" "UBERON:0000468" "UBERON:0000473" "UBERON:0000948" ...
##   .. [list output truncated]
##   ..- attr(*, "dim")= int 18715
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr [1:18715] "ENSDARG00000000001" "ENSDARG00000000002" "ENSDARG00000000018" "ENSDARG00000000019" ...
##  $ organ.relationships:List of 12587
##   ..$ AEO:0000013   : chr "UBERON:0000479"
##   ..$ AEO:0000127   : chr "UBERON:0005423"
##   ..$ AEO:0000129   : chr "UBERON:0005423"
##   ..$ AEO:0000131   : chr "UBERON:0005423"
##   ..$ AEO:0000153   : chr "AEO:0000187"
##   ..$ AEO:0000162   : chr "UBERON:0006984"
##   ..$ AEO:0000173   : chr [1:2] "UBERON:0002416" "UBERON:0000020"
##   ..$ AEO:0000180   : chr "UBERON:0002418"
##   ..$ AEO:0000185   : chr [1:2] "UBERON:0004288" "UBERON:0007475"
##   ..$ AEO:0000187   : chr "UBERON:0000481"
##   ..$ AEO:0000188   : chr [1:2] "UBERON:0000463" "UBERON:0000476"
##   ..$ AEO:0000189   : chr "AEO:0000187"
##   ..$ AEO:0000196   : chr "UBERON:0002548"
##   ..$ AEO:0000210   : chr "UBERON:0000480"
##   ..$ AEO:0000211   : chr [1:2] "UBERON:0000477" "UBERON:0005423"
##   ..$ AEO:0000215   : chr "UBERON:0010064"
##   ..$ AEO:0000217   : chr "UBERON:0005423"
##   ..$ AEO:0001003   : chr "UBERON:0000476"
##   ..$ AEO:0001006   : chr "UBERON:0000463"
##   ..$ AEO:0001007   : chr [1:2] "UBERON:0000483" "UBERON:0005423"
##   ..$ AEO:0001009   : chr "UBERON:0004022"
##   ..$ AEO:0001010   : chr "UBERON:0004022"
##   ..$ AEO:0001013   : chr [1:2] "UBERON:0010314" "UBERON:0002020"
##   ..$ AEO:0001014   : chr "UBERON:0005423"
##   ..$ AEO:0001019   : chr "UBERON:0005423"
##   ..$ AEO:0001020   : chr "UBERON:0000479"
##   ..$ AEO:0001021   : chr "UBERON:0000479"
##   ..$ CARO:0000049  : chr "UBERON:0010895"
##   ..$ CARO:0000050  : chr "UBERON:0010895"
##   ..$ CARO:0000064  : chr "CL:0000003"
##   ..$ CL:0000000    : chr "GO:0005575"
##   ..$ CL:0000001    : chr "CL:0000010"
##   ..$ CL:0000003    : chr [1:3] "GO:0005623" "CL:0000000" "UBERON:0000061"
##   ..$ CL:0000005    : chr [1:2] "UBERON:0010313" "CL:0000057"
##   ..$ CL:0000006    : chr [1:2] "CL:0000197" "CL:0000101"
##   ..$ CL:0000007    : chr "CL:0002321"
##   ..$ CL:0000008    : chr "CL:0000333"
##   ..$ CL:0000010    : chr "CL:0000578"
##   ..$ CL:0000011    : chr [1:2] "UBERON:0003083" "CL:0000333"
##   ..$ CL:0000014    : chr [1:2] "CL:0000039" "CL:0000034"
##   ..$ CL:0000015    : chr "CL:0000586"
##   ..$ CL:0000016    : chr [1:2] "CL:0000014" "CL:0000015"
##   ..$ CL:0000017    : chr [1:2] "UBERON:0000473" "CL:0000015"
##   ..$ CL:0000018    : chr [1:2] "CL:0000015" "CL:0000413"
##   ..$ CL:0000019    : chr [1:3] "CL:0000064" "CL:0000219" "CL:0000408"
##   ..$ CL:0000020    : chr [1:2] "UBERON:0000473" "CL:0000015"
##   ..$ CL:0000021    : chr "CL:0000586"
##   ..$ CL:0000022    : chr [1:2] "CL:0000014" "CL:0000021"
##   ..$ CL:0000023    : chr [1:2] "UBERON:0001305" "CL:0000021"
##   ..$ CL:0000024    : chr "CL:0000021"
##   ..$ CL:0000025    : chr "CL:0000675"
##   ..$ CL:0000026    : chr [1:2] "CL:0000039" "CL:0000412"
##   ..$ CL:0000027    : chr [1:2] "UBERON:0010313" "CL:0000192"
##   ..$ CL:0000028    : chr "CL:2000029"
##   ..$ CL:0000029    : chr [1:2] "UBERON:0010313" "CL:0000540"
##   ..$ CL:0000030    : chr [1:2] "CL:0000055" "UBERON:0004121"
##   ..$ CL:0000031    : chr [1:2] "CL:0000047" "CL:0000055"
##   ..$ CL:0000032    : chr [1:2] "CL:0002321" "UBERON:0004121"
##   ..$ CL:0000033    : chr "CL:0000151"
##   ..$ CL:0000034    : chr "CL:0000003"
##   ..$ CL:0000035    : chr "CL:0000723"
##   ..$ CL:0000036    : chr "CL:0000035"
##   ..$ CL:0000037    : chr [1:3] "UBERON:0002390" "CL:0000048" "CL:0008001"
##   ..$ CL:0000038    : chr [1:2] "CL:0000764" "CL:0000839"
##   ..$ CL:0000039    : chr "CL:0000548"
##   ..$ CL:0000040    : chr [1:3] "CL:0000763" "CL:0002194" "UBERON:0004120"
##   ..$ CL:0000041    : chr [1:2] "CL:0000771" "CL:0002242"
##   ..$ CL:0000042    : chr [1:2] "CL:0000835" "CL:0000834"
##   ..$ CL:0000043    : chr [1:3] "CL:0000767" "CL:0002274" "CL:0002242"
##   ..$ CL:0000047    : chr [1:3] "CL:0000048" "CL:0002319" "UBERON:0001017"
##   ..$ CL:0000048    : chr "CL:0000723"
##   ..$ CL:0000049    : chr [1:2] "UBERON:0004120" "CL:0002032"
##   ..$ CL:0000050    : chr [1:3] "CL:0000839" "CL:0000763" "CL:0002032"
##   ..$ CL:0000051    : chr [1:2] "UBERON:0004120" "CL:0002032"
##   ..$ CL:0000052    : chr "CL:0000723"
##   ..$ CL:0000054    : chr "CL:0000151"
##   ..$ CL:0000055    : chr "CL:0000003"
##   ..$ CL:0000056    : chr "CL:0000680"
##   ..$ CL:0000057    : chr [1:2] "UBERON:0004120" "CL:0002320"
##   ..$ CL:0000058    : chr [1:2] "CL:0000667" "CL:0000055"
##   ..$ CL:0000059    : chr [1:3] "CL:0000151" "CL:0000146" "UBERON:0005176"
##   ..$ CL:0000060    : chr [1:4] "UBERON:0001754" "CL:0000055" "CL:0002159" "UBERON:0005176"
##   ..$ CL:0000061    : chr [1:4] "UBERON:0000078" "CL:0000062" "CL:0000151" "UBERON:0005176"
##   ..$ CL:0000062    : chr [1:4] "UBERON:0010313" "CL:0000055" "CL:0002320" "UBERON:0004120"
##   ..$ CL:0000064    : chr "CL:0000003"
##   ..$ CL:0000065    : chr [1:3] "UBERON:0001016" "CL:0000067" "CL:0000710"
##   ..$ CL:0000066    : chr [1:2] "UBERON:0000483" "CL:0000548"
##   ..$ CL:0000067    : chr [1:2] "CL:0000075" "CL:0000064"
##   ..$ CL:0000068    : chr [1:3] "CL:0000066" "CL:0002371" "UBERON:0000058"
##   ..$ CL:0000069    : chr "CL:0000068"
##   ..$ CL:0000071    : chr [1:3] "CL:0000076" "CL:0002139" "UBERON:0004638"
##   ..$ CL:0000072    : chr [1:2] "CL:0000068" "UBERON:0004120"
##   ..$ CL:0000073    : chr [1:2] "CL:0000066" "CL:0002371"
##   ..$ CL:0000075    : chr [1:2] "CL:0000066" "CL:0002371"
##   ..$ CL:0000076    : chr "CL:0000066"
##   ..$ CL:0000077    : chr [1:3] "CL:0000213" "CL:0000076" "UBERON:0004120"
##   ..$ CL:0000078    : chr [1:2] "CL:0000076" "CL:0002371"
##   ..$ CL:0000079    : chr [1:2] "CL:0000066" "CL:0002371"
##   ..$ CL:0000080    : chr [1:2] "CL:0000003" "UBERON:0000179"
##   .. [list output truncated]
##  $ organ.names        :'data.frame': 12588 obs. of  2 variables:
##   ..$ ID  : chr [1:12588] "AEO:0001009" "AEO:0001010" "AEO:0001013" "CL:0000005" ...
##   ..$ NAME: chr [1:12588] "proliferating neuroepithelium" "differentiating neuroepithelium" "neuronal column" "fibroblast neural crest derived" ...
##  $ bgee.object        :Reference class 'Bgee' [package "BgeeDB"] with 13 fields
##   ..$ species          : chr "Danio_rerio"
##   ..$ speciesName      : chr "Danio_rerio"
##   ..$ speciesId        : num 7955
##   ..$ dataType         : chr [1:4] "rna_seq" "affymetrix" "est" "in_situ"
##   ..$ pathToData       : chr "/tmp/Rtmpy2umFm/Rbuild63a2eeaa29/BgeeDB/vignettes/Danio_rerio_Bgee_13_2"
##   ..$ release          : chr "13_2"
##   ..$ annotationUrl    : chr(0) 
##   ..$ experimentUrl    : chr(0) 
##   ..$ allExperimentsUrl: chr(0) 
##   ..$ topAnatUrl       : chr "http://r13_2.bgee.org/"
##   ..$ sendStats        : logi TRUE
##   ..$ quantitativeData : logi FALSE
##   ..$ apiKey           : chr "5b534b6e616600e6ea5445ae289801ef455ba79e2d5d9da478fab78551afbab3d942fa64e3346172580b8967c6ecd9da65294a3abd55a43d7b9814dbb0357ad"| __truncated__
##   ..and 18 methods, of which 4 are  possibly relevant:
##   ..  format_data, get_annotation, get_data, initialize

The strigency 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.

## Loading only high-quality 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")
## myTopAnatData <- loadTopAnatData(bgee, stage="UBERON:0000068", confidence="high_quality")

Note: 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.

Prepare a topGOdata-like object allowing to perform TopAnat analysis

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, annotated with “spermatogenesis” in the Gene Ontology (using the biomaRt package). We expect these genes to be enriched for expression in male tissues, notablythe 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.

# source("https://bioconductor.org/biocLite.R")
# biocLite("biomaRt")
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("drerio_gene_ensembl", mart=ensembl)

# Foreground genes are those with GO annotation "spermatogenesis"
myGenes <- getBM(attributes= "ensembl_gene_id", filters=c("go_id"), values=list(c("GO:0007283")), mart=ensembl)

# Background are all genes with GO annotation
universe <- getBM(attributes= "ensembl_gene_id", filters=c("with_go_go"), values=list(c(TRUE)), mart=ensembl)

# Prepare the gene list vector 
geneList <- factor(as.integer(universe[,1] %in% myGenes[,1]))
names(geneList) <- universe[,1]
head(geneList)
## ENSDARG00000063895 ENSDARG00000063899 ENSDARG00000063905 
##                  0                  0                  0 
## ENSDARG00000063908 ENSDARG00000063910 ENSDARG00000063911 
##                  0                  0                  0 
## Levels: 0 1
summary(geneList == 1)
##    Mode   FALSE    TRUE    NA's 
## logical   19766      31       0
# Prepare the topGO object
myTopAnatObject <-  topAnat(myTopAnatData, geneList)
## 
## Checking topAnatData object.........
## 
## Checking gene list..................
## 
## WARNING: Some genes in your gene list have no expression data in Bgee, and will not be included in the analysis. 13684 genes in background will be kept.
## 
## Building 'most specific' Terms......  ( 1051 Terms found. )
## 
## Build DAG topology..................  ( 1841 terms and 3704 relations. )
## 
## Annotating nodes (Can be long)......  ( 13684 genes annotated to the nodes. )

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:

results <- runTest(myTopAnatObject, algorithm = 'classic', statistic = 'fisher')
## 
##           -- Classic Algorithm -- 
## 
##       the algorithm is scoring 314 nontrivial nodes
##       parameters: 
##           test statistic: fisher

You can also choose one of the topGO decorrelation methods, for example the “weight” method, allowing to avoid redundant results induced by the structure of the ontology.

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.

# Display results sigificant at a 5% FDR threshold
makeTable(myTopAnatData, myTopAnatObject, results, cutoff = 0.05)
## 
## Building the results table for the 17 significant terms at FDR threshold of 0.05...
## Ordering results by pValue column in increasing order...
## Done
##                       organId                  organName annotated
## UBERON:0000079 UBERON:0000079   male reproductive system      4135
## UBERON:0000473 UBERON:0000473                     testis      4135
## UBERON:0003135 UBERON:0003135    male reproductive organ      4135
## UBERON:0003101 UBERON:0003101              male organism      4138
## CL:0000255         CL:0000255            eukaryotic cell       170
## CL:0000548         CL:0000548                animal cell       170
## CL:0000003         CL:0000003                       cell       171
## GO:0005623         GO:0005623                       cell       171
## CL:0000000         CL:0000000                       cell       173
## GO:0005575         GO:0005575         cellular_component       173
## CL:0000007         CL:0000007       early embryonic cell        84
## CL:0000353         CL:0000353            blastoderm cell        84
## CL:0002321         CL:0002321             embryonic cell        84
## UBERON:0000474 UBERON:0000474 female reproductive system      4049
## UBERON:0000992 UBERON:0000992               female gonad      4049
## UBERON:0003134 UBERON:0003134  female reproductive organ      4049
## UBERON:0003100 UBERON:0003100            female organism      4051
##                significant expected foldEnrichment       pValue        FDR
## UBERON:0000079          18     7.86       2.290076 4.625880e-05 0.01233300
## UBERON:0000473          18     7.86       2.290076 4.625880e-05 0.01233300
## UBERON:0003135          18     7.86       2.290076 4.625880e-05 0.01233300
## UBERON:0003101          18     7.86       2.290076 4.676018e-05 0.01233300
## CL:0000255               4     0.32      12.500000 2.776418e-04 0.03131243
## CL:0000548               4     0.32      12.500000 2.776418e-04 0.03131243
## CL:0000003               4     0.32      12.500000 2.839259e-04 0.03131243
## GO:0005623               4     0.32      12.500000 2.839259e-04 0.03131243
## CL:0000000               4     0.33      12.121212 2.968003e-04 0.03131243
## GO:0005575               4     0.33      12.121212 2.968003e-04 0.03131243
## CL:0000007               3     0.16      18.750000 5.239205e-04 0.04251817
## CL:0000353               3     0.16      18.750000 5.239205e-04 0.04251817
## CL:0002321               3     0.16      18.750000 5.239205e-04 0.04251817
## UBERON:0000474          16     7.69       2.080624 7.112340e-04 0.04440597
## UBERON:0000992          16     7.69       2.080624 7.112340e-04 0.04440597
## UBERON:0003134          16     7.69       2.080624 7.112340e-04 0.04440597
## UBERON:0003100          16     7.70       2.077922 7.155465e-04 0.04440597

There is only one significant term, testis, which makes biological sense: there is an expression bias for testis of genes involved in spermatogenesis.

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.

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.