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

BgeeDB is a collection of functions to import into R re-annotated, quality-controlled and re-processed expression data available in the Bgee database. This includes data from thousands of wild-type healthy samples of multiple animal species, generated with different gene expression technologies (RNA-seq, Affymetrix microarrays, expressed sequence tags, and in situ hybridizations). BgeeDB facilitates downstream analyses, such as gene expression analyses with other Bioconductor packages. Moreover, BgeeDB includes a new gene set enrichment test for preferred localization of expression of genes in anatomical structures (“TopAnat”). Along with the classical Gene Ontology enrichment test, this test provides a complementary way to interpret gene lists. Availability: https://www.bioconductor.org/packages/BgeeDB/


Introduction
Gene expression levels influence the behavior of cells, the functionality of tissues, and a wide range of processes from development and aging to physiology or behavior. It is of particular importance that researchers are able to take advantage of the vast amounts of publicly available gene expression datasets to reproduce and validate results, or to investigate new research questions [1][2][3] .
To that purpose, one should be able to easily query and import gene expression datasets generated using different technologies, and their associated metadata. The R environment 4 has now become a standard for bioinformatics and statistical analysis of gene expression data, through the Bioconductor framework and its many open source packages 5,6 . It is thus desirable to provide an access to gene expression datasets programmatically and directly into R. For example, the Bioconductor packages ArrayExpress 7 , GEOquery 8 and SRAdb 9 provide access to the reference databases ArrayExpress 10 , GEO 11 and SRA 12 respectively. However, such databases are primary archives aiming at comprehensiveness. They include gene expression datasets and other functional genomics data, generated from diverse experimental conditions, of diverse quality. The data provided are heterogeneous, with some datasets including only unprocessed raw data, and others including only data processed using specific analysis pipelines. For instance, over the 44,481 RNA array assay experiments stored in ArrayExpress with processed data available as of June 2018, 7,544 do not include the raw data. Metadata are often provided as free-text information that is difficult to query. For instance, the GEO database encourages submitters of high-throughput sequencing experiments to provide MINSEQE elements, but does not enforce this practice (see, e.g., GEO submission guidelines, and GEO Excel template for submissions). Unless the user needs to retrieve a specific known dataset from its accession number, it can be difficult to identify relevant available datasets. This can ultimately constitute an obstacle to data reuse.
One response to this diversity of primary archives is topical databases 1 . They can be useful for researchers of specialized fields, and even more so if they propose an R package for data access. For example, the BrainStars Bioconductor package allows access to microarray data of mouse brain regions samples from the BrainStars project 13,14 . The ImmuneSpaceR Bioconductor package allows access to the gene expression data generated by the Human Immunology Project Consortium 15 . Such efforts allow a better control of the data and annotation quality, but by nature they include a limited number of conditions, which only fit the needs of specialized projects. Similarly, numerous "ExperimentData" packages are available on the Bioconductor repository, which each include a single curated and well-formatted expression dataset (see https://www.bioconductor.org/packages/ release/BiocViews.html#___ExpressionData). But these packages are rarely updated and are mostly meant to be used as examples in software packages vignettes, for teaching, or to provide supplementary data of publications. The package ExperimentHub 16 also provides access to a central location where this type of single datasets can be retrieved, but it does not address the difficulty of integrating datasets annotated and processed in different ways.
Finally, added-value databases aim at filtering, annotating, and possibly reprocessing all or some of the datasets available from the primary archives 1 . For example, there is a Bioconductor package to access the Expression Atlas, which includes a selection of microarray and RNA-seq datasets from ArrayExpress that are re-annotated and

Amendments from Version 1
We thank the reviewers for their work, and we feel that, thanks to their comments, the manuscript has been greatly improved. We have updated considerably the Bgee database, and the corresponding documentation. This addresses the comments made by the reviewers about a lack of transparency of our data processing steps. Moreover, we have set up the documentation is such a manner as to insure that it remains updated with the progress of the database in the future. We have added information about the processing of gene expression data performed at the Bgee database. We notably now link to the source code of our pipeline for data processing. We have added some information about similar tools allowing to perform gene list expression localization enrichment analyses. We also have updated the examples and results based on the use of the latest Bgee release and BgeeDB package version. The code examples have been made more robust to potential future changes to the format of the files used by the package. We felt that it was necessary to link the revised publication of the BgeeDB package to the updated documentation. We have added an author, Julien Wollbrett, to the author list. Since this manuscript was first submitted, Julien Wollbrett made significant contributions to the development of the package described in this paper, and notably towards the aim of submitting this revised manuscript. Please also note that we have updated all figures and supplementary files, so that they are based on the latest releases of our database and R package.
The Bgee database (https://bgee.org/) 22 is another added-value database, which currently offers access to reprocessed gene expression datasets from 29 animal species. Bgee aims at comparisons of gene expression patterns across tissues, developmental stages, ages and species. It provides manually curated annotations to ontology terms, describing precisely the experimental conditions used. It integrates expression data generated with multiple technologies: RNA-Seq, Affymetrix microarrays, in situ hybridization, and expressed sequence tags (ESTs) in release 14. An important characteristic of Bgee is that all datasets are manually curated to retain only "normal" healthy wild-type samples, i.e., excluding gene knock-out, treatments, or diseases. Finally, Bgee datasets are carefully checked for quality issues, and reprocessed to produce normalized expression level, calls of presence/absence of expression, and of differential expression. Bgee thus provides a reference of high-quality and reusable gene expression datasets that are relevant for biological insights into normal conditions of gene expression. The release 14.0 of Bgee covers 29 animal species, and includes 5,745 RNA-seq libraries, 12,996 Affymetrix chips, 360,653 results from 49,241 in situ hybridization experiments, and 3,335 EST libraries. This includes 4,860 human RNA-Seq libraries from the GTEx project 23,24 .
Until 2016 the Bgee database lacked a programmatic access to data through a R package, a shortcoming that we have addressed with the release of the BgeeDB Bioconductor package, available at https://www.bioconductor.org/packages/BgeeDB/. The package provides functions for fast extraction of data and metadata. The data structures used in the package can be easily incorporated with other Bioconductor packages, offering a wide range of possibilities for downstream analyses.
Moreover, we introduce in BgeeDB the possibility to run TopAnat analyses, i.e., anatomical expression enrichment tests on gene lists provided by the user. This functionality is based on the topGO package 25,26 , modified to use Bgee data (A. Alexa, personal communication). TopAnat is similar to the widely used Gene Ontology enrichment test [27][28][29] . But in our case, the enrichment test is applied to terms from an anatomical ontology, mapped to genes by expression patterns. The reference set of genes in a given species consists of all genes for which at least one "present" expression call is available in Bgee. The expression calls are propagated to parent anatomical structures by part_of and is_a relations, using the Uberon anatomical ontology 30,31 (e.g., a gene expressed in the "hindbrain" is also considered expressed in the parent structure "brain"). Different algorithms, from TopGO, are available in TopAnat to account for the non-independence of anatomical structures, and avoid the overrepresentation of lowly-informative top-level terms. Enrichment of expression is tested for each anatomical structure independently with a Fisher exact test, and the resulting p-values for all anatomical structures are then corrected using a FDR correction 32 . As a result, TopAnat allows to discover the tissues where a set of genes is preferentially expressed. This feature is available as a web-tool at https://bgee.org/?page=top_anat, but the R package offers more flexibility in the choice of input data and analysis parameters, and possibilities of inclusion within programs or pipelines.
There exist few other tools allowing to perform anatomical expression enrichment tests. For instance, the webapplication Tissue Specific Expression Analysis (TSEA 33 based on methods from refs 34,35 ) allows to perform such analyses, but only in human and mouse, while TopAnat can be used for any species integrated in Bgee (29 species as of Bgee release 14.0). For human, TSEA is based on the Genotype-Tissue Expression (GTEx) RNA-Seq dataset, while Bgee integrates GTEx data, but also other RNA-Seq datasets, and datasets from different data types, providing a higher diversity of anatomical structures. TSEA was last updated on March 2014. The database wormbase also proposes a similar tool, but for analyses only in C. elegans 36 . There exists another application for expression enrichment analyses, but focused on analyzing gene regulatory networks in human 37 .
The pipeline to process the data accessible through the BgeeDB package is documented in detail at https://github. com/BgeeDB/bgee_pipeline. In brief, for RNA-seq experiments: data present in SRA 12 are selected and annotated using information from GEO 11 or from papers, or provided by the Model Organism Database Wormbase 38 . GTF annotation files and genome sequence fasta files are retrieved from Ensembl and Ensembl Genomes Metazoa 39,40 . After quality control steps, the Kallisto software is used to perform a pseudo-mapping of the reads to the transcriptome 41 . TMM normalization 42 is used to normalize TPM and F/RPKM values within each experiment independently. Present/absent expression calls are produced for each library by comparing the level of expression of each gene to the background transcriptional noise in the library (estimated by using the level of expression of intergenic regions; Roux J., Rosikiewicz M., Wollbrett J., Robinson-Rechavi M., Bastian F.B.; in preparation). In brief, for Affymetrix experiments: data present in ArrayExpress and GEO are selected and annotated using the information available in these repositories, or in papers, or provided by the Model Organism Database Wormbase. Mappings of probesets to genes are retrieved from Ensembl and Ensembl Genomes Metazoa. Quality controls are performed to remove low quality chips and redundant chips 43,44 . When raw data are available, they are normalized using gcRMA (using version 2.42.0 for Bgee release 14.0) within each experiment independently 45 . Present/absent expression calls are generated either from the MAS5 processed data 46 , based on the perfect match/mismatch probesets, or using the raw data when available, by comparing the signal of a probeset to a subset of weakly expressed probesets 47 .
The BgeeDB package information is available on the Bioconductor website at https://bioconductor.org/packages/BgeeDB/. The source code is available at https://github.com/BgeeDB/BgeeDB_R. The preferred location for filing bug reports and suggestions is the issue tracker on GitHub.
In the following sections we provide some typical examples of usage of the BgeeDB package. • Working internet connection Please note that an earlier version of Bioconductor and of R (>= 3.3) could be used, but would require to clone our GitHub repository and to use the R CMD BUILD command to build the package.

Use cases
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 to work offline.
Microarray dataset retrieval. In the following example, we will look for a microarray dataset in mouse (Mus musculus), spanning multiple early developmental stages, including zygote.
# specify species and data type # the examples in this paper are based on Bgee release 14.0 # the following line targets the latest Bgee release. # In order to target specifically the release 14.0, # add the parameter 'release="14.0"' 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) This creates a list of two data frames, one including the annotation of experiments, and the other including the annotation of each individual sample, i.e., hybridized microarray chip. For mouse, there are 698 Affymetrix experiments and 6,095 samples available in Bgee release 14.0. Anatomical structures and developmental stages are annotated using the Uberon ontology 30,31 . Sex and strain information is also provided. Below, we are selecting the experiments for which at least one sample is annotated to the zygote stage (UBERON:0000106). 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 48 , includes samples from more developmental stages than the other two, so we will choose this one for the next steps. For this experiment, raw data were available from ArrayExpress, so samples were fully normalized with gcRMA 49 trough 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. The developmental stage 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.
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.
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 will 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. As for Affymetrix data, Bgee RNA-seq annotations provide information about anatomical structure, developmental stage, sex, and strain.

Presence/absence calls retrieval.
It is often difficult to compare expression levels across species 54 , and even within species, across datasets generated by different experimenters or laboratories 55-57 . 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 combination "gene / tissue / developmental stage / sex / strain".
A powerful use of these expression calls is the anatomical expression enrichment test "TopAnat". TopAnat uses a similar approach to Gene Ontology enrichment tests 27 , but genes are associated to the anatomical structures where they display expression, instead of to their functional classification. These tests allow discovering 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.; in preparation). 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 58,59 .

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

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 (in Bgee 14, 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 silver and gold quality are used, which can be changed with the confidence argument of the loadTopAnatData() function (in releases of Bgee up to 13, "high" and "low" confidence were used). 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 will use expression calls for zebrafish genes using all sources of expression data.
# the examples in this paper are based on Bgee # release 14.0 # the following line targets the latest Bgee release. In order to target # specifically the release 14.0, add the parameter 'release="14.0"' bgee.topanat <-Bgee$new(species="Danio_rerio") top.anat.data <-loadTopAnatData(bgee.topanat) We will 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 64 and integrated into the Ensembl database 39 . We will thus retrieve the targeted genes using the biomaRt 65 Bioconductor package (version 2.36.1 for this paper).
biocLite("biomaRt") library(biomaRt) # zebrafish data in Ensembl 84, that Bgee 14.0 uses (stable link) ensembl <-useMart("ENSEMBL_MART_ENSEMBL", dataset="drerio_gene_ensembl", host="mar2016.archive.ensembl.org") # get the mapping of Ensembl genes to phenotypes genes.to.phenotypes <-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(genes.to.phenotypes$phenotype_description), value=T) # select the genes annotated to select phenotypes genes <-unique(genes.to.phenotypes$ensembl_gene_id[genes.to.phenotypes$phenotype_ description %in% phenotypes]) This gives a list of 147 zebrafish genes implicated in the development and function of pectoral fin. The next step of the analysis relies on the topGO Bioconductor package. We will 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 66 .
Here we choose as background all zebrafish Ensembl genes with an annotated phenotype from ZFIN.
# prepare the gene list vector gene.list <-factor(as.integer(unique(genes.to.phenotypes$ensembl_gene_id) %in% genes)) names(gene.list) <-unique(genes.to.phenotypes$ensembl_gene_id) summary(gene.list) # prepare the topAnat object based on topGO top.anat.object <-topAnat(top.anat.data, gene.list) top.anat.object 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 can 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 26 . Here we will use a Fisher test, coupled with the "weight" decorrelation algorithm.
results <-runTest(top.anat.object, algorithm='weight', statistic='fisher') results Finally, we implemented 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 32 ) 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 table.Over <-makeTable(top.anat.data, top.anat.object, results, cutoff=0.01) head(table.over) The 27 anatomical structures displaying a significant enrichment at a FDR threshold of 1% are show in Table 1. The first term is "pectoral fin", and the second "paired limb/fin bud". Other terms in the list, especially those with high enrichment folds, are clearly related to pectoral fins (e.g., "pectoral appendage field"), or substructures of fins (e.g., "fin bone"). This analysis shows that genes with phenotypic effects on pectoral fins are specifically expressed in or next to these structures. More generally, it demonstrates the relevance of TopAnat analysis for the characterization of lists of genes.
Of note, it is possible to retrieve for a particular tissue the significant genes that were mapped to it.

Conclusion
In summary, the BgeeDB package serves as a bridge between curated data from the Bgee database and the R/Bioconductor environment, facilitating the access to high-quality curated and re-analyzed gene expression datasets, and significantly reducing time for downstream analyses of the datasets. Moreover, it provides access to TopAnat, a new enrichment tool allowing to make sense of lists of genes, by uncovering their preferential localization of expression in anatomical structures. The TopAnat workflow is straightforward; for users already using topGO in their analysis pipelines, performing a TopAnat analysis on the same gene list only requires 6 additional lines of code. Author contributions AK and JR contributed equally to this work. AK developed the initial BgeeDB R package and made it available in Bioconductor. JR implemented the enrichment analyses, and refined the data download part. JW corrected some bugs, updated the package and the files used by the package, to use the latest Bgee release and for better compatibility between operating systems. FBB developed the server-side responses. MRR and FBB tested and commented on the package development. AK and JR wrote the manuscript. All authors discussed the results and implications and commented on the manuscript at all stages.

Competing interests
No competing interests were disclosed.

Grant information
This work was supported by SIB Swiss Institute of Bioinformatics project Bgee, Swiss National Science Foundation grant 31003A_153341, SystemsX.ch project AgingX, and Etat de Vaud.

Supplementary material
File S1. R markdown file including code from the paper.
Click here to access the data.
File S2. PDF file including the results of execution of the code from File S1.
Click here to access the data. I want to highlight all the work the authors did ensuring that their paper is reproducible (they mention versions of packages throughout the text) and in clarifying their work (mostly updates on the introduction). I was able to download supplementary file 1 (an .Rmd file) and execute it without any problems or modifications at all using R 3.5.1 with Bioconductor 3.7 on a Mac (BgeeDB 2.6.2, MFuzz 2.40.0, biomaRt 2.36.1). The authors did a great job with https://github.com/BgeeDB/bgee_pipeline which I believe includes all the information needed for anyone who is interested in the Bgee pipeline. If not, the authors seem responsive on https://github.com/BgeeDB/bgee_pipeline/issues and https://github.com/BgeeDB/BgeeDB_R/issues which is a great sign. Thanks to their changes in the introduction I now have a better understanding of their anatomical expression enrichment test.

Open Peer Review
I look forward to seeing how the community uses BgeeDB in the future and the use cases they apply it to.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. In this manuscript the authors describe the BgeeDB Bioconductor package and show how to use it (as of Bioconductor 3.4) to interact with Bgee in order to get the data from Bgee into your R session. This allows users to then perform differential expression analyses and integrate Bgee with other data sets such as unpublished data. The manuscript includes code that shows how to use BgeeDB and showcases it's 1 as unpublished data. The manuscript includes code that shows how to use BgeeDB and showcases it's different features including their unique anatomical expression enrichment analysis method. I find interesting that you can use BgeeDB to get data from different platforms and from different organisms. Most of this can be done using other packages such as GEOquery, but BgeeDB makes it so the user doesn't have to do all the processing of the data and standarization over multiple projects.
My main concern with the manuscript in its current form and the BgeeDB package itself is the lack of clarity on how the data has been processed and how the anatomical expression test works. That is, it could potentially become a black box that produces interesting output but hides information that could be important.
For example, I'm sure some of the Affymetrix data could be downloaded with other packages and I do not know what would be the differences between the raw data and the data downloaded via BgeeDB. Is the data in BgeeDB normalized? If so, how? The help pages of BgeeDB, the package vignette, the original Bgee publication and did not help me fully answer these questions (I might http://bgee.org/?page=doc have missed the information). Maybe the functions in BgeeDB could print a message describing the main steps of how a given data set was processed or this could be added to the help pages. I currently ignore if all data sets were treated the same. For instance, is all the Affymetrix data normalized with the same method and same parameters? I assume that the answer is yes but I don't know. I suggest that the authors describe in more detail the data available in Bgee. The authors might want to consider making the processing code public at or citable via . https://github.com/BgeeDB figshare With the anatomical expression test it's not clear to me how to interpret the results from BgeeDB::makeTable(). I understand that the authors will describe the details of how their anatomical test works in a future publication, which they did before with Bgee and Homolanto . Ideally, the anatomical expression test would have been described first followed by BgeeDB. Without hindering the current plan, I believe that the authors could provide a summary of how TopAnat works. Then they can explain it fully in the planned future TopAnat publication. I am also curious on how users could use their own data to improve the TopAnat results, although that could be work for the TopAnat paper or future work.
I think that the manuscript is overall well written and will be more appealing if the data and main features (TopAnat) are described in more detail. I hope that the authors are not discouraged by my report.

Best, Leonardo
Minor comments: I'm an author of recount which is incorrectly cited here. The pre-print version of had data from 2040 different projects which together https://jhubiostatistics.shinyapps.io/recount/ made up more than 60,000 RNA-seq samples. The current version has data from over 70,000 Illumina human RNA-seq samples from SRA, GTEx and TCGA. I don't think that it makes sense to include the str() calls in the paper. They do make sense in the supplementary material (the html and pdf rendered versions of the paper code) since those include the output. Also, while str() shows all the details of an object, it can encourage users to write code that depends on the internal structure of the object. You might want to consider adding accessor 1 1 2 3 that depends on the internal structure of the object. You might want to consider adding accessor functions.
If you added indentation the code that runs over multiple lines would be easier to read. You can use the Bioconductor standard of using 4 spaces at the start of the line. Also make sure that object names don't get split into multiple lines. For example check the line after the "list experiments including both brain and liver samples" comment where "Anatomical.entity.ID" gets split into "Anatomical.e" and "ntity.ID" in the html version of the paper. Copy pasting works fine, but if someone prints the paper they might introduce errors can be avoided with better formatting. F1000Research's team should be able to tell you what is the character limit per line to use so that the PDF and HTML versions look great. The package might be useful here. formatR I would not use numerical indexes in the code since the results could change with time in such a way that the current code would not work in the future or worse, it might run without error but change the results in a way a new user would not notice. For example, change the code on the line after the "order developmental stages" comment which currently reads: data.E.MEXP.51.formatted <-data.E.MEXP.51.formatted [, c(5,8,9,3,2,1,4,7,6)] The comment that reads with "retrieve anatomical structures enriched at a 1% FDR threshold" is mixed with the code. That is, you are missing a new line character.
The package's vignette is missing a title as currently shown at . http://bioconductor.org/packages/release/bioc/html/BgeeDB.html I recommend adding internal R links to your manual pages. For example, ?topAnat mentions loadTopAnatData(). Those links make it easier for a user to browse the help pages.
I was able to run all the code without any edits (beyond that new line issue I already mentioned) using Regarding Virag Sharma's peer review report , I assume that Virag was using an earlier R version (and thus an earlier Bioconductor version). The current development version of BgeeDB uses "dataType" and not "datatype", just like the release version. Check . Hopefully the https://github.com/Bioconductor-mirror/BgeeDB/search?utf8=%E2%9C%93&q=datatype authors won't change the spelling of arguments in the future since that's confusing for users, although that's certainly doable following the deprecated/defunct code cycle.
Regarding Daniel S. Himmelstein's peer review report , there is no need to add a license file when the license is specified in the DESCRIPTION file of an R package. See 4 5 license is specified in the DESCRIPTION file of an R package. See where they state that https://github.com/Bioconductor-mirror/BgeeDB/blob/master/DESCRIPTION#L14 the license is GPL-2. Although the authors should make sure that they correctly specify which license their software is released on: GPL-2 or GPLv3 as Daniel mentioned. Regarding where to place bug reports, the authors could resolve this by specifying the "BugReports" field in their DESCRIPTION file. For example see . I also agree with https://github.com/Bioconductor-mirror/recount/blob/master/DESCRIPTION#L63 Daniel that currently BgeeDB has a bit of a messy download structure. I would prefer if the files were downloaded in a single directory (say "bgee_downloads") instead of the current working directory. My main concern with the manuscript in its current form and the BgeeDB package itself is the lack of clarity on how the data has been processed and how the anatomical expression test works. That is, it could potentially become a black box that produces interesting output but hides information that could be important. For example, I'm sure some of the Affymetrix data could be downloaded with other packages and I do not know what would be the differences between the raw data and the data downloaded via BgeeDB. Is the data in BgeeDB normalized? If so, how? The help pages of BgeeDB, the package vignette, the original Bgee publication1 and http://bgee.org/?page=doc did not help me fully answer these questions (I might have missed the information).
We have made public the Bgee pipeline source code at https://github.com/BgeeDB/bgee_pipeline. We also have added a paragraph at the end of the "Introduction" section, pointing to the relevant part of the documentation for RNA-Seq and Affymetrix analyses, and describing them in brief.

---
Maybe the functions in BgeeDB could print a message describing the main steps of how a given data set was processed or this could be added to the help pages. data set was processed or this could be added to the help pages.
We have opened an issue on our tracker related to this point, see https://github.com/BgeeDB/BgeeDB_R/issues/22. We will add a function pointing to the relevant documentation in a future release.
---I currently ignore if all data sets were treated the same. For instance, is all the Affymetrix data normalized with the same method and same parameters? I assume that the answer is yes but I don't know.
The Affymetrix data are not treated in the same way depending on whether the raw data were available, or only the data processed by using the MAS5 software. This is clarified at the end of the "Introduction" section. Also, in the package, this information about raw data availability can be retrieved in the annotation data frame.
---I suggest that the authors describe in more detail the data available in Bgee. The authors might want to consider making the processing code public at https://github.com/BgeeDB or citable via figshare.

---
With the anatomical expression test it's not clear to me how to interpret the results from BgeeDB::makeTable(). I understand that the authors will describe the details of how their anatomical test works in a future publication, which they did before with Bgee and Homolanto. Ideally, the anatomical expression test would have been described first followed by BgeeDB. Without hindering the current plan, I believe that the authors could provide a summary of how TopAnat works. Then they can explain it fully in the planned future TopAnat publication.
We have added a brief description of how TopAnat works in the "Introduction" section.
---I am also curious on how users could use their own data to improve the TopAnat results, although that could be work for the TopAnat paper or future work.
This represents an advanced use of TopAnat that we don't find suitable for the paper. But users can override the association file, mapping genes to anatomical structures in the BgeeDB directory, to use their own data. Also, since the source code of the package is public, users can also modify the mapping files used by modifying the source code.
---I'm an author of recount which is incorrectly cited here. The pre-print version of https://jhubiostatistics.shinyapps.io/recount/ had data from 2040 different projects which together https://jhubiostatistics.shinyapps.io/recount/ had data from 2040 different projects which together made up more than 60,000 RNA-seq samples. The current version has data from over 70,000 Illumina human RNA-seq samples from SRA, GTEx and TCGA.
We have updated the number in our paper. We apologize for the mistake.
---I don't think that it makes sense to include the str() calls in the paper. They do make sense in the supplementary material (the html and pdf rendered versions of the paper code) since those include the output. Also, while str() shows all the details of an object, it can encourage users to write code that depends on the internal structure of the object. You might want to consider adding accessor functions.
We have removed str() calls from the paper. For the future, we will think of adding accessor functions, although several are already available thanks to the topGO package.

---
If you added indentation the code that runs over multiple lines would be easier to read. You can use the Bioconductor standard of using 4 spaces at the start of the line. Also make sure that object names don't get split into multiple lines. For example check the line after the "list experiments including both brain and liver samples" comment where "Anatomical.entity.ID" gets split into "Anatomical.e" and "ntity.ID" in the html version of the paper. Copy pasting works fine, but if someone prints the paper they might introduce errors can be avoided with better formatting. F1000Research's team should be able to tell you what is the character limit per line to use so that the PDF and HTML versions look great. The formatR package might be useful here.
We didn't know about the formatR package and will have a look at it. In the meantime, we have split such offending lines, as identified by the reviewer.
---I would not use numerical indexes in the code since the results could change with time in such a way that the current code would not work in the future or worse, it might run without error but change the results in a way a new user would not notice. For example, change the code on the line after the "order developmental stages" comment which currently reads: data.E. c(5,8,9,3,2,1,4,7,6)] We have replaced all lines using numerical indexes, with use of column names.

---
The comment that reads with "retrieve anatomical structures enriched at a 1% FDR threshold" is mixed with the code. That is, you are missing a new line character.

Reference 46 is incorrect. It's edgeR, not eedgeR.
This was fixed.
---I recommend adding internal R links to your manual pages. For example, ?topAnat mentions loadTopAnatData(). Those links make it easier for a user to browse the help pages.
We thank the reviewer for the suggestion, and we will implement this in a future release. This is indeed a change that we introduced in an earlier version, in an effort to name all our arguments in a consistent manner. We will try not to change this in the future.
---I also agree with Daniel that currently BgeeDB has a bit of a messy download structure. I would prefer if the files were downloaded in a single directory (say "bgee_downloads") instead of the current working directory.
While another directory can be specified by using the "pathToData" argument, it is true that the solution proposed by the reviewer would be convenient, and we will try to update the package accordingly in the future.
No competing interests were disclosed. This study describes the BgeeDB R package, which provides a programmatic interface for accessing Bgee gene expression data. Bgee is a valuable resource because it integrates gene expression results across many experiments.
for its presence/absence of expression calls and its Previously, I've used Bgee differential expression calls.
In my opinion, Bgee's ability to provide a genome-wide profile of expression for a given species, developmental stage, and anatomical structure is its most powerful capability. It was not clear to me whether BgeeDB provides this functionality. For example, can the user retrieve the normalized expression level across several experiments for the same species-stage-anatomy combination? In general, I think users will be more interested in this high-level functionality than the low level access BgeeDB currently provides. An example here would likely clear things up for me.
Is it possible to integrate expression levels across Affymetrix and RNA-Seq experiments?
The of the source code specifies as the license. This is great, but it's ideal to also In my opinion, Bgee's ability to provide a genome-wide profile of expression for a given species, developmental stage, and anatomical structure is its most powerful capability. It was not clear to me whether BgeeDB provides this functionality. For example, can the user retrieve the normalized expression level across several experiments for the same species-stage-anatomy combination? In general, I think users will be more interested in this high-level functionality than the low level access BgeeDB currently provides. An example here would likely clear things up for me.
Indeed there is currently no easy way to do this. As mentioned in , it would be nice to have a getDataByCondition https://github.com/BgeeDB/BgeeDB_R/issues/7 function that would return all processed data for chips / libraries matching a queried organ/stage/(sex)/(strain). But it was hard to set priorities for the initial development (should the package complement the web interface, or be orthogonal to it?), and we will likely implement it in the near future.
If the reviewer means to integrate levels of expression, it is then not the aim of Bgee: Bgee integrate different data types and different experiments, processed and normalized independently (but in a consistent manner).

---
The Zenodo archive of the source code specifies GPLv3 as the license. This is great, but it's ideal to also add a LICENSE file to the GitHub.
We have added the LICENSE file to GitHub (GPL 3.0).

---
It looks like there are at least two potential places where bug reports should be filed: on Bioconductor Support and GitHub Issues. It would be nice to clarify the preferred location for filing bug reports go and opening pull requests.
We have added the preferred location for filing bug reports at the end of the "Introduction" section (GitHub), and in the DESCRIPTION file of the source code.
---Currently, the GitHub repository BgeeDB/BgeeDB_R mentioned in the manuscript is forked from wirawara/BgeeDB. I expect this may cause some confusion, as BgeeDB/BgeeDB_R should be the upstream repository that users fork and contribute back to. If you make wirawara/BgeeDB private, this should break the relationship. @wirawara can then fork BgeeDB/BgeeDB_R to continue contributions if desired.
We thank the reviewer for the suggestion, we have now made wirawara/BgeeDB private.
---Finally, I created some GitHub issues as part of this review: Sample annotation variable names https://github.com/BgeeDB/BgeeDB_R/issues/5 We have replied on the issue. Our answer was that is a bit of a controversial topic. For example Google's R Style Guide (https://google.github.io/styleguide/Rguide.xml#identifiers) advise against the use of underscores (although they do not justify why, and we agree that the "words separated 1.

2.
3. the use of underscores (although they do not justify why, and we agree that the "words separated with dots" convention can be disturbing for python users).

A less messy default download directory
This point was discussed in https://github.com/BgeeDB/BgeeDB_R/issues/4. We notably mention that another directory can be specified by using the "pathToData" argument. This parameter is mentioned at the end of the section "Data download and import of normalized expression levels". We agree that a default directory should be used in future releases.
No competing interests were disclosed. In the manuscript, Komljenovic et al. present BgeeDB which is an R package for retrieval of expression datasets which have been curated. Additionally, they also provide a method (TopAnat) to determine tissue-specific enrichments for a given list of genes and species.
The former is a very useful resource because there is clearly a need for a database that provides gene expression datasets which are homogenous in nature and are of comparable quality. The BgeeDB database contains gene expression datasets from 17 species across different tissues and developmental stages, which is impressive. The fact that the database can be queried via a Bioconductor package should ensure that the database will be used by both -wet-lab biologists and computational scientists. Similarly, the TopAnat method also provides a useful functionality to determine anatomical expression enrichment on a user specified list.
I have a few minor comments regarding the manuscript: The authors should include some details about how they have reprocessed the gene expression datasets that are a part of BgeeDB. At the moment, it is rather unclear how this was achieved. I assume that the authors have an automated pipeline in place but it would be beneficial for readers to know how this was done.
The authors state that "TopAnat allows for discovery of tissues where a set of genes is preferentially expressed". Is TopAnat the only tool that offers such a functionality? A brief background of similar tools that are currently available will be useful for the readers.
I was not able to run the workflow that the authors have included in the Supplementary material: We have added this example at the end of the "Anatomical expression enrichment analysis" section.
No competing interests were disclosed.

Competing Interests:
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias You can publish traditional articles, null/negative results, case reports, data notes and more The peer review process is transparent and collaborative Your article is indexed in PubMed after passing peer review Dedicated customer support at every stage For pre-submission enquiries, contact research@f1000.com