PulmonDB: a curated lung disease gene expression database

Chronic Obstructive Pulmonary Disease (COPD) and Idiopathic Pulmonary Fibrosis (IPF) have contrasting clinical and pathological characteristics and interesting whole-genome transcriptomic profiles. However, data from public repositories are difficult to reprocess and reanalyze. Here, we present PulmonDB, a web-based database (http://pulmondb.liigh.unam.mx/) and R library that facilitates exploration of gene expression profiles for these diseases by integrating transcriptomic data and curated annotation from different sources. We demonstrated the value of this resource by presenting the expression of already well-known genes of COPD and IPF across multiple experiments and the results of two differential expression analyses in which we successfully identified differences and similarities. With this first version of PulmonDB, we create a new hypothesis and compare the two diseases from a transcriptomics perspective.


Results
PulmonDB is a relational database implemented in MySQL with lung disease transcriptome measurements, re-annotated platform probes, and manually curated data with a controlled vocabulary designed for lung diseases ( Fig. 1). Tables were created to describe each feature and to connect the information across experiments, samples, measurements, platforms, genes, and annotated information. The full database scheme is provided in Supplementary Fig. 1. pulmonDB a curated gene expression lung disease database. PulmonDB is a curated gene expression database of human lung diseases, with RNA-seq and microarray data from different platforms that have been uniformly preprocessed and manually curated to add sample and experiment information. In addition, we developed a website to access and visualize homogenized data (http://pulmondb.liigh.unam.mx/), and we also developed an R package (https://github.com/AnaBVA/pulmondb) to download curated annotation and preprocessed data that can be used for further analysis in the R environment.
Our database has a total of 76 GSEs, corresponding to 4481 unique preprocessed GSM contrasts that used 26 different platforms or GPLs (platform ID from GEO) (Fig. 2C). PulmonDB contains different sample types, we searched for human gene expression experiments related to COPD and IPF without any restriction. Lung biopsies account for 37.8% of samples, and 33.2% are blood samples. However, different cell types can be found in PulmonDB: some of them are primary cells (e.i. alveolar macrophages, fibroblasts, alveolar epithelial cells, etc.), and others are cell lines (e.i. A549) ( Fig. 2A). Of the samples, 34.9% correspond to COPD, 40.5% to control samples (30.9% healthy plus 9.6% match tissue), 17.2% to IPF, and 1.5% to other diseases (Fig. 2B and Supplementary  Table 2). We separated control tissues into two groups, "healthy" individuals, as far as the authors are aware and "match_tissue_controls" which refers to tissue samples from a phenotypically healthy region of a patient who had a tumor removed (e.i. non-tumor tissue from a cancer patient).
Although other resources reuse and reanalyze GEO data using web interfaces 9 , those tools are not specialized for lung diseases. Their limitations include the need for previous manual curation in each analysis, and they consider a small number of COPD and IPF experiments due to the fact that only curated GEO data are used. We designed a web interface that enables data exploration and visualization to facilitate lung disease analysis. This interface uses Clustergrammer 13 to visualize gene expression values and the creation of interactive heatmaps that allow data exploration. A valuable feature of Clustergrammer is to be connected to EnrichR 14 , which provides pathway enrichment analysis. All these features together should help to generate new hypotheses about the pathologies of lung diseases to perform exploratory analyses, to visualize specific gene expression across public experiments for comparing results, and to generate new insights based on different data sets. pulmonDB can recapitulate gene expression patterns expected in copD and ipf. To show that PulmonDB can be used to recapitulate previously reported knowledge regarding COPD and IPF biology, we performed a literature search and manually selected relevant genes for each disease. We selected 19 genes related to IPF (not necessarily associated with gene expression in lung tissues) to visualize their gene expression: CCL18 15 27 . Then, we selected eight IPF experiments performed with lung tissue biopsy samples (GSE32537, GSE21369, GSE24206, GSE94060, GSE72073, GSE35145, GSE31934), and using the PulmonDB website, we created a heatmap with the gene expression patterns and observed that the hierarchical clustering of these data separates IPF and control data sets (Fig. 3A, green and gray clusters at the bottom). For COPD, we curated 16 genes from the literature that were deemed relevant to this disease: HHIP 28 32 , and GATA2 32 . We selected five experiments (GSE27597, GSE37768, GSE57148, GSE8581, GSE1122) performed on lung tissue biopsy samples from COPD patients and controls. Our hierarchical clustering analysis of the expression profiles using the PulmonDB interface allowed us to cluster patients and controls into two different groups (Fig. 3B), similar to the case of IPF. In conclusion, PulmonDB not only helps to recapitulate previously published work (Supplementary PulmonDB was created using COMMAND by downloading, parsing and storing COPD and IPF public transcriptomic data into a MySQL database. Then, we remapped microarray probes to establish a uniform gene annotation, and we also created a controlled vocabulary for clinical and biological annotations for each sample. We created contrasts based on the original hypothesis, selecting a sample as the reference. Finally, the data were homogenized and subjected to a quality check. www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 3) but also helps to verify gene expression stability across experiments. This may help to analyze concordance in different experiments, contrast study results, show implications of using different control groups, etc. We believe this resource can be used to drive, make decisions, and support new hypotheses in experimental laboratories for studying molecular or cellular disease mechanisms.

Differences and similarities in COPD and IPF.
PulmonDB can be used not only to replicate previous knowledge but also to provide a framework to test new hypotheses. In this context, we set out to investigate the differences and similarities between COPD and IPF in lung tissue when compared to samples from healthy individuals (Fig. 4A). Using PulmonDB in the R environment, we selected contrasts where the sample was annotated as lung biopsy and the reference status as HEALTHY/CONTROLs (GSE52463, GSE63073, GSE1122, GSE72073, GSE24206, GSE27597, GSE29133, GSE31934, GSE37768) (Fig. 4B), and then using limma 38 we assessed differential gene expression between COPD and IPF. We identified 1781 differentially expressed genes ( Supplementary  Fig. 4). To have a visual representation of the differences between COPD and IPF, we selected the top 20 differentially expressed genes and visualized their expression using the PulmonDB website tool (Fig. 4C). We observed that data sets tend to cluster by test status; Fig. 4C shows IPF contrasts on the left (turquoise), control contrasts in the middle (blue), and COPD contrasts on the right (red). Genes are clustered in two groups (left panel, y-axis); the first gene group (I) is overexpressed in IPF while it is barely expressed or underexpressed in COPD contrasts. By comparison, the second gene cluster (group II) is overexpressed in COPD contrasts and underexpressed in IPF. To correlate similarities among samples, the 20 top differentially expressed genes were used (Fig. 4C, right panel); samples from the same disease group showed higher correlations and tended to have a null or negative correlation with the HEALTHY/CONTROL and the opposite disease (Fig. 4C). For example, FOSB and CXCL2 have opposite behaviors, as both genes are overexpressed in COPD and underexpressed in IPF. FOSB is part of the family of Fos genes that can dimerize with JUN family proteins to form the transcription factor complex AP-1, which is related to COPD 39 . CXCL2 is a chemokine secreted in inflammation that induces chemotaxis in www.nature.com/scientificreports www.nature.com/scientificreports/ neutrophils 40,41 ; these cells are predominant in COPD, and they are key mediators in tissue damage 42 . While neutrophils are also important in IPF, we observed their underexpression in this disease.
We also asked the opposite question, i.e., whether we could identify which genes that are shared between these two diseases. We assigned a weight to COPD and IPF expression to perform limma contrasts (Fig. 4D), which enabled us to identify when both diseases drove a differential expression profile. We selected the 20 top differentially expressed genes and visualized their expression patterns using PulmonDB website tool, and we could see that a set of genes was consistently overexpressed or underexpressed in both COPD and IPF. In particular, VCAM1 and FCN3 are differentially expressed in COPD and IPF, with a similar trend in both diseases when compared with HEALTHY/CONTROLs. VCAM1 is the vascular cell adhesion molecule 1, and it is important in the immune response for mediating cellular adhesion in leukocytes 43 ; it is overexpressed in these two diseases, suggesting infiltration of immune cells in both pathologies 44,45 . In contrast, FCN3 (or ficolin 3) is underexpressed in both diseases: this gene is a collagen-like protein associated with the innate immune defense, as it activates the lectin complement pathway 46 , which has been shown to be important in pulmonary pathologies 47,48 . www.nature.com/scientificreports www.nature.com/scientificreports/ As a result, PulmonDB assisted our analysis of COPD and IPF analogous and antagonist genes and can thus be used to dissect common molecular mechanisms, because both lung diseases are present under heterogeneous conditions with progressive and irreversible phenotypes mainly caused by smoking and by aging, plus both diseases entail cellular matrix remodeling. Furthermore, the differential gene signatures between COPD and IPF might explain the particularities of each disease. In both (C, D) columns are sample contrasts, rows are genes, the first covariate is colored by each corresponding experiment, the second covariate is the sample type (in this case, lung tissue is shown in light brown), the third row is the test status, and the fourth is the reference status. Columns are ordered by test status and genes by hierarchical clusterization. The right heatmap is the correlation among sample contrasts, and the covariates are the same.

Discussion
The present methodology had been previously applied for the study of bacterial and grapevine gene expression in different experiments and conditions, allowing for the integration of data from a diverse origin. Here we prove this methodology can also be applied to human data to exploit publicly available resources better, we hope these methods will be taken by other teams to create databases to help understand relevant diseases in other tissues.
PulmonDB can help the scientific community to study which genes have a distinct expression profile in COPD and IPF, explore experiments across technologies and platforms, identify interesting expression patterns across different diseases, generate new hypotheses, and find relationships among clinical or experimental variables. This database also enables comparisons of an updated collection of expression profiles already homogenized for their analyses of specific diseases. Additionally, having different lung diseases (COPD and IPF) in the same database creates the opportunity to observe their similarities and differences. In the future, we aim for PulmonDB to grow and include more diseases. To our knowledge, there is no other resource for transcriptomic analysis focused on the same lung diseases; for this reason, we believe researchers of different backgrounds can use and benefit from the information contained in PulmonDB, by using the web interface and the R package.
An integrated comparable collection of homogenized values with controlled vocabulary describing biological and technical characteristics will facilitate further comparative analyses, such as the study of profiles in COPD and IPF, exploration of experiments across technologies and platforms, identification of interesting coexpression patterns across different diseases, the generation of new hypotheses, and determination of relationships among clinical or experimental variables.
This project sets the foundation to integrate transcriptomics data of other respiratory diseases or related phenotypes and thus facilitates the identification of common and divergent pathways that lead to a pathological state. PulmonDB platform will be expanded in the future to include other lung diseases.
Methods platform and metadata. Most of the metadata was obtained from GEO. For specific cases, the platform information (.cdf file) was obtained from the Affymetrix website (http://www.affymetrix.com/site/mainPage. affx). Additional information (e.g., clinical data, source of the biological sample), was obtained either from metadata or manually curated from the original papers.
inclusion criteria for transcriptome data. The experiments currently included in PulmonDB are listed in Supplementary Table 2.
We used two main resources to download raw data and preprocessed counts, GEO and Recount2.
Gene expression omnibus. Using GEO 2,3 , we searched datasets related to COPD and IPF for gene expression data. The following queries were used to retrieve the experiments: GEO experiments were manually curated, abstracts and related articles were revised, and only datasets confirmed as having COPD and/or IPF samples were considered. In order for an experiment to be included in PulmonDB we used the following criteria: The data set had to be original, samples had to be unique, raw data had to be public and available, platform information must had sequence probes, and custom platforms must have had information to link raw expression signal with the probe sequence. Otherwise, data sets were not taken into account for PulmonDB.
Recount2. Recount2 is an online resource with RNA-seq human experiments already preprocessed using Rail-RNA alignment and summarized by gene and exon counts 49 . We used the keywords "IPF" and "COPD" separately in Recount2 to retrieve counts form RNA-seq. compendium creation. The compendium creation process was done as previously described in COLOMBOS and VESPUCCI 10,12 . The platform was developed in bacteria and later employed in grapevine, but in this paper, we used COLOMBOS for the first time in human data. After we selected the datasets using the experiment ID from GEO (GSE), we worked on COMMAND>_ 50 .
COMMAND. COMMAND stands for COMpendia MANagement Desktop, it is a web application tool that provides a framework to facilitate and perform the following steps: (1) download data from selected experiments, (2) parse files and store data in database form, (3) probe-to-gene (re)mapping process, (4) sample curation and annotation with a controlled vocabulary, (5) selection of references and sample experiments to determine contrasts, (6) homogenization (and normalization) of data, and (7) perform data quality control (Fig. 1). This software can be used for any transcriptomic data 50 .
In more detail, each experiment with a GSE ID, also referred to as a data set, was normalized independently without performing background correction, as explained in 11 . We defined a contrast for each sample with a GSM ID (sample ID from GEO) by using a unique control reference sample per data set. The sample contrast per gene was defined as the log ratio between the expression value in the test condition (i.e., IPF, COPD) and the (2020) 10:514 | https://doi.org/10.1038/s41598-019-56339-5 www.nature.com/scientificreports www.nature.com/scientificreports/ expression value in the reference condition (i.e., healthy, untreated, smokers without COPD) (Fig. 1, step 5). This gives every comparison an interpretable biological meaning when combined with extensive manual curated annotation. The condition properties describing the contrasts were then structured in a condition-controlled vocabulary tree. Finally, all contrasts were homogenized, resulting in direct comparable log ratios across all experiments; this information later became part of the final compendium of expression data ( Supplementary Fig. 2). pulmonDB uses a controlled vocabulary to describe sample metadata. A controlled vocabulary is required to create databases with homogeneous and standard information. For PulmonDB, we created a controlled vocabulary organized in a hierarchical structure that contains terms to annotate transcriptome experiments in lung diseases. We defined classes describing the main categories and terms that can be found in experiments, with some of them as mandatory features (e.i. sample type, sample status, and platform). Some non-IPF or non-COPD diseases were included in the controlled vocabulary because the original experiments used them.
Once the controlled vocabulary was established, each article related to the experiment was manually curated, and whenever it was necessary, new terms were added, making the vocabulary flexible and allowing for the inclusion of other diseases to our database in the future. Complete definitions of the terms are provided in Supplementary Table 1. experiment annotation. Each sample was manually annotated using the controlled vocabulary; when necessary, the vocabulary was updated to include new features. The information was curated by experts who reviewed the associated articles and protocols to retrieve data such as age, sex, ancestry, stage of disease or treatment, DLCO (the diffusing capacity of the lung for carbon monoxide, a common functional test), etc., from either GEO or the associated paper.
Homogenization and quality control. As described before, data homogenization was done with COMMAND>_ 11,12 . This step was performed on raw data without background correction, as it has been shown to retrieve more errors [51][52][53] . A nonlinear model was applied to homogenize raw data. We used RMA Quantile for Affymetrix samples and loess fit for the other platforms. The next step was to summarize probes per transcript using RMA median polish summary from Affymetrix or with data averaged across replicates for the other platforms. After performing the homogenization step, low-quality microarrays were identified using MA plots and histograms of raw and homogenized data.

Website implementation.
PulmonDB has a web interface that uses Clustergrammer (https://clustergrammer.readthedocs.io/index.html) 13 to visualize gene expression contrasts. Clustergrammer has a frontend in javascript and a backend in python, supporting an interactive web application for gene expression exploration. The PulmonDB web interface requires one or several GSE identifiers and more than two gene names to generate interactive heatmaps.
In addition, Clustergrammer is connected with EnrichR (http://amp.pharm.mssm.edu/Enrichr/) 14 , an integrative web application tool for enrichment analysis that helps the user explore not only potentially differentiated genes but also enriched pathways, facilitating the discovery of transcriptomic signature patterns in lung diseases or related phenotypes. copD and ipf comparative analysis. We used limma 3.40.0 in a Rstudio environment 3.6.0 for our comparative analyses, and the GSE ID was included in the linear model. Then, two contrasts were created: (1) "COPD -IPF", for obtaining differentially expressed genes between COPD and IPF, and (2) "(COPD + IPF)/2 -CONTROL", for genes similarly expressed between COPD/IPF and CONTROL. Differential gene expression analyses were adjusted for multiple testing using the false discovery rate (FDR) method, also referred to as Benjamini & Hochberg adjustment. We applied a cutoff of the adjusted p-value < 0.05, and after sorting based on the log fold change, the top 20 genes were obtained.