BioXpress: an integrated RNA-seq-derived gene expression database for pan-cancer analysis

BioXpress is a gene expression and cancer association database in which the expression levels are mapped to genes using RNA-seq data obtained from The Cancer Genome Atlas, International Cancer Genome Consortium, Expression Atlas and publications. The BioXpress database includes expression data from 64 cancer types, 6361 patients and 17 469 genes with 9513 of the genes displaying differential expression between tumor and normal samples. In addition to data directly retrieved from RNA-seq data repositories, manual biocuration of publications supplements the available cancer association annotations in the database. All cancer types are mapped to Disease Ontology terms to facilitate a uniform pan-cancer analysis. The BioXpress database is easily searched using HUGO Gene Nomenclature Committee gene symbol, UniProtKB/RefSeq accession or, alternatively, can be queried by cancer type with specified significance filters. This interface along with availability of pre-computed downloadable files containing differentially expressed genes in multiple cancers enables straightforward retrieval and display of a broad set of cancer-related genes. Database URL: http://hive.biochemistry.gwu.edu/tools/bioxpress

More than a decade ago, gene expression analysis was proposed as a method to complement classification schemes based on tumor morphology because it was well known that tumors with similar histopathological appearance can have considerably different clinical outcomes (6,9).
These efforts provided the framework by which linking gene expression with cancer research could be realized (10). Hanahan and Weinberg (11) in their seminal paper 'The Hallmarks of Cancer' discussed the role of over-and under-expression of key genes in several cancers. The conjectures that both diagnosis of somatically acquired lesions in tumors and genome-wide expression profiling of tumors would become routine (11) have not yet been realized, but we anticipate that this will likely change within the next decade. With advances in next-generation sequencing (NGS) technologies, several national and international projects are underway that aim to capture and analyze the expression profiles of thousands of tumors (12)(13)(14). Additionally, there are already thousands of publications that describe over-and under-expression of specific genes in cancer. Currently, to the best of our knowledge, there is no integrated view of the expression profiles of the human genes obtained from NGS technology such as RNA sequencing (RNA-seq). Moreover, no singular effort is underway to manually curate data from publications on cancer-related gene expression, enabling easy comparison of expression data and knowledge from both small publications and large-scale studies like The Cancer Genome Atlas (TCGA: http://cancergenome.nih.gov/) and International Cancer Genome Consortium (ICGC: https:// icgc.org/). Lack of such efforts prevents us from tracking our knowledge of expression profiles of genes in different cancer types as technology improves and more data and information accumulate. Furthermore, as we move toward the translation of expression analysis through genomic or proteomic technologies to the clinic, there is no easy way to compare a patient's expression data with that extant data. BioXpress has been developed as the first step toward the provision of easy access to gene expression data from tumor and normal samples, which will be useful for clinical research, diagnostics and prognostics of cancer.
The specific technology used to measure gene expression significantly affects the cost, comprehensiveness and the time consumed to perform expression analysis. DNA microarray and quantitative polymerase chain reaction (q-PCR) are powerful approaches for measuring gene expression and have been used for many years. DNA microarray technology is efficient and cost-effective at the gene expression level, while q-PCR is considered more sensitive. However, neither of these approaches can meet the sensitivity and comprehensiveness of the newer RNA-seq technology (15,16). Despite the benefits of RNA-seq, microarrays are often preferentially used due to the higher cost and lack of standardization of pipelines using the RNA-seq technology. Once these obstacles are overcome, it is clear that RNA-seq will become the predominant tool for expression analysis (17). In addition to expression analysis, RNA-seq provides a number of other benefits. A single RNA-seq experiment output can aid in the discovery of novel and unannotated transcripts (18), single nucleotide variation (SNV) identification (19) and more (20). As RNA-seq technology and the corresponding analytical approaches grow, the application of this method is becoming indispensable for many scientific disciplines (21)(22)(23). To address this growing presence of RNA-seq data, we currently focus on large-scale integration of RNAseq-based expression data in BioXpress complemented by manual curation of information from publications reporting gene expression associated with cancer. The manual curation process allows us to collect valuable expressionrelated information from peer-reviewed publications from diverse platforms. Integration of information from both large-scale studies and publications allows users to easily compare and contrast expression profiles of their gene(s) of interest.
The advancement of expression analysis technology has led to the development of corresponding databases and standards. For example, the Minimum Information About a Microarray Dataset initiative (24) provides standards for microarray data, while databases like the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (25) and Array Express (26) have significant amounts of microarray data. Secondary databases that store and provide results and analysis of microarray and other gene expression data related to cancer such as CGED (Cancer Gene Expression Database) (27), GENT (Gene Expression across Normal and Tumor tissue) (28) and Oncomine (29) are also available. Finally, TCGA and ICGC data portals and databases, such as Expression Atlas (30), provide RNA-seq-generated data. All the above-mentioned databases provide mechanisms to retrieve gene-specific information, but, to the best of our knowledge, none of them allows integrated pan-cancer analysis across multiple projects. NCBI GEO and European Bioinformatics Institute (EBI) ArrayExpress, e.g. are public repositories for high-throughput microarray and NGS functional genomic datasets. A gene symbol-based search can result in thousands of profiles from GEO Profiles Database. CGED, on the other hand, provides data specifically obtained through collaborative efforts of Nara Institute of Science and Technology, Osaka University Medical School, Kyoto University Medical School and Osaka Medical Center for Cancer and Cardiovascular Diseases. GENT provides Affymetrix microarray data from tumor and normal samples, while Expression Atlas at EBI provides differential and baseline expression from several organisms. Similar to other public repository data, a single search can retrieve data from many experiments in these resources. Furthermore, although tools do exist which aim to analyze the same scope of data, the tools and databases of which we are aware do not facilitate the expression analysis on RNA-seq desired here. cBioPortal (31) is a widely popular resource with an emphasis on mutation analysis. Currently, Oncomine (29) does provide the means to analyze expression for microarray data, but not for RNA-seq. Thus, there is no single tool/resource available which integrates RNA-seq information that allows expression analysis to identify, e.g. cancer relatedness. Although hundreds of cancer RNA-seq studies are published each year, a cancer-centric RNA-seq expression database that integrates all cancer-related RNA-seq-based expression data from databases and publications is not available to the community. Portals like TCGA data portal and ICGC data portal, which provide RNA-seq-based expression data, only provide access to raw read counts and normalized counts: such data cannot be easily used for comparative analysis across several cancer types and existing experimental results in publications. In addition, different normalization methods are employed by different data providers, making comparison and cross-type analysis even more challenging.
The BioXpress database project collects RNA-seq data from several publicly available sources such as TCGA (http://cancergenome.nih.gov/), ICGC (12) and Gene Expression Atlas (30), and uses a standardized method to identify the expression levels of the genes. Expression levels of genes are also manually extracted from publications to supplement information gathered from large-scale studies. Additionally, all cancer types are mapped to Disease Ontology (32) terms to facilitate pan-cancer analysis. Finally, all genes are linked to a comprehensive cancerrelated non-synonymous SNV database, BioMuta (33). Together, BioMuta and BioXpress provide a detailed view of the expression and mutations of genes in cancer and therefore can be used for pan-cancer studies like the one performed by our group recently (34) and described in this manuscript.

Data Source and Metrics
The majority of RNA-seq databases provide data either in FASTQ format (sequence reads) and/or raw read count data. Read count data are calculated by analyzing the mapping file where the reads have already been aligned to a reference genome. As shown in Figure 1, BioXpress processes data based on the availability of expression data from paired data that have both normal and tumor samples  Figure 1. Flow chart of the workflow used to create BioXpress. BioXpress processes short reads and read count data through distinct pipelines. Data are further divided into two groups: paired data that have both normal and tumor samples from the same patient, and non-paired, tumor-only data. Output in BioXpress is split into three different types: differential expression (stacked bar chart), tumor-only expression (box plot) and baseline expression data (heatmap). In addition to the data integration approaches shown in the figure, gene expression information is also extracted from publications.
from the same patient, and the non-paired data from just tumor and also from just normal tissue. The data sources and statistics in terms of number of patients from each data source are shown in Table 1. To achieve comprehensiveness, data are collected from TCGA, the Curated Short Read archive (CSR) (35), ICGC (12), Gene Expression Atlas (30) and publications (Table 1). It is important to note that ICGC, at the time of writing this article, did not contain any data from paired normal and tumor samples which are not from TCGA. Therefore, the data in BioXpress are split into three different types: differential expression, tumor-only expression and baseline expression data from Illumina Human Body Map project (http://www.ebi.ac.uk/gxa/experiments/E-MTAB-513).

Data Processing
TCGA data portal TCGA-Assembler was used to download RNA-seq data from TCGA data portal. Raw counts data with paired samples (tumor and normal) were extracted and analyzed using DEseq R package with default parameters: method ¼ 'blind', sharingMode ¼ 'fit-only', fitType ¼ 'local' (36). DEseq normalization method has been reported to outperform other normalization methods (37). Fold changes, not absolute expression values, are displayed based on analysis described above (38). False discovery rates are not defined due to the low number of replicates for samples. This approach allows the user to determine the significance of differentially expressed genes on an individual basis.

ICGC data portal
ICGC contains tumor-only data (normal samples are not sequenced by the consortium currently). Gene expression data from tumor samples was downloaded from ICGC data portal (12) and analyzed using DEseq R package with default parameters (36).

Expression atlas
Normalized baseline expression was downloaded via Expression Atlas (http://www.ebi.ac.uk/gxa/download. html) (30). Because raw read counts are not available for all data retrieved from Expression Atlas, no additional normalization was performed in BioXpress.

Manual curation from publications
Decades of research on differential expression in tumor and normal samples has led to thousands of publications. Although many of these studies are based on samples from modest numbers of patients, there is value in the systematic capture and presentation of this information alongside largescale studies such as those presented by TCGA and ICGC. Although it is possible that studies may exhibit discordance, it is equally possible for the consideration of such additional experiments to contribute to the 'big picture' of differential expression between tumor and normal samples. We leave it to the discretion of individual users to decide the significance of curated publications in application to their studies.
For manual curation of expression data, genes identified in our previous pan-cancer study were prioritized (34). In addition to this prioritization, genes annotated by UniProtKB/Swiss-Prot as associated with cancer and Cancer Gene Census (http://www.sanger.ac.uk/genetics/ CGP/Census/) (39) were also targeted for manual curation. This UniProtKB/Swiss-Prot gene list was obtained using the following search string: organism: 'Homo sapiens [9606]' AND reviewed:yes AND annotation:(type:disease cancer). Briefly, the manual curation protocol involved searching PubMed (40) using the gene name (including synonyms) with accompanying text 'cancer' and 'expression'. The curator then reviewed the title to shortlist articles which appear to contain gene expression information related to cancer and have full text available. Abstracts were then read to identify potential true positive articles. All such articles were downloaded and read to extract key information such as cancer type and expression information. All cancer types were then mapped to Disease Ontology terms (32) and added to the BioXpress database. To date, 536 papers have been filtered to maintain only those focusing on human cancer after reading the 'Abstract' and 'Introduction'. Among this subset, only papers including direct evidence reflecting gene expression  Table S2).

Data Normalization and Analysis
DEseq method is regarded as one of the most robust RNAseq normalization methods (37). In the BioXpress pipeline, raw counts data were normalized by DEseq method followed by differentially expressed gene analysis. To compare non-paired samples with normalized results from DEseq pipeline, the DEseq normalization method was used [Parameters: library('DESeq'), cds ¼ newCountDataSet (data,condition), cds ¼ estimateSizeFactors(cds), result ¼ counts(cds,normalized ¼ TRUE)]. For differential expression analysis, gene expression was normalized based on each patient, and case and control were considered together. For tumor expression, all samples were collectively analyzed across different cancer types and then normalized. Heat map and clustering analysis were performed using the 'heatmap' function from the R package (http:// www.R-project.org/).

Usage and Utility
Scientists can find querying datasets useful to identify expression levels between disease and normal pairs to discover differential expression for a gene. They may also want to research on potential biomarkers or pathways that lead to tumor formation or want to explore the overall expression of specific genes across multiple cancer types. Users can search BioXpress using HGNCapproved gene symbols (HUGO Gene Nomenclature Committee), UniProtKB/Swiss-Prot accessions or RefSeq accessions. Differentially expressed genes for a specific cancer type can also be retrieved. Additionally, all data in BioXpress, including lists of genes significantly differentially expressed in two or more cancer types, can be downloaded.
Searching using gene name (gene/protein-centric search) A search using the HGNC-approved gene symbol or UniProt/RefSeq accession retrieves differential expression information (cancer vs. normal), tumor-only expression data (where normal samples are not available) and baseline expression information from normal human tissues (Illumina Human Body Map Project). The example below provides an overview of a gene/protein-centric search.

Differential expression
The abnormal spindle-like microcephaly-associated (ASPM) gene is highly expressed in several tumor cell lines (41) and cancers (42,43). Searching the BioXpress database using the gene ASPM users can retrieve the differential expression profile of this gene in different cancers. For ASPM gene, we can clearly see that this gene appears to be over-expressed in almost all cancers. Figure 2 provides a view of the BioXpress interface where the Differential Expression tab on the top menu bar is selected, and below it 'ASPM Expression Profile' is shown. The default view provides expression frequency (over-or under-expression) in the patients. The number of patients for a particular cancer type, P value and a variety of additional information is available in the table below which can be downloaded. Full cancer names are available on clicking the cancer abbreviations in figure and additional details about the data can be viewed by clicking the ' Table column description' link. All columns can be sorted and users can send an email to the help desk with comments about a specific data element by clicking on the envelope link available from each row. The tab at the top of the stacked bar chart provides an alternate view where users can see the frequency (number of patients) of significantly over-or under-expressed genes (based on a P value cutoff of 0.05). For ASPM, on clicking the Significant/Freq tab, we can see that this gene is significantly over-expressed in more than 25% of the patients in several cancers. For example, ASPM is over-expressed in breast invasive carcinoma (DOID:3459; 113 patients), lung adenocarcinoma (DOID:3907; 50 patients) and others. Combining the stacked bar frequency expression (Regulation/Freq) and the Significant/Freq, users can get a complete overview of the differential expression of a gene in all cancer types in the database.

Tumor expression
Clicking on the Tumor Expression tab on the top menu bar shows the expression profile for the ASPM gene from all patient samples without paired normal data. Although ICGC does not currently collect any paired data, tumor-only expression data can provide an overview of the expression of a specific gene in different cancer types and can be used in conjunction with differential and baseline expression data to better understand the comprehensive expression profile of a gene. The box plot provides the minimum lower quartile, median upper quartile and maximum expression value, and therefore provides a snapshot of the distribution of expression of a gene in all patients with a specific cancer. For the ASPM gene, we see that for cervical squamous cell carcinoma (CESC), the minimum, maximum and the lower and upper quartile are above the theoretical mean for all cancer types which could indicate that for CESC this gene has less variability in terms of expression in the patients and is expressed at a higher level compared with other cancers. Therefore, the box plot allows the user to identify cancer types where the lower and the upper quartile are short, signifying homogeneity in the expression of the gene for that specific cancer. The table below the box plot provides details such as UniProtKB accession, RefSeq accession and number of samples.

Baseline expression
Clicking the Baseline Expression tab for ASPM gene shows the heatmap with testis being the only tissue with increased expression of ASPM. It has been known for some time that ASPM is over-expressed in testis (41,44), although the precise function of this gene in testis development is still unknown (45).

Searching using cancer type (cancer type centric search)
Users may want to retrieve a list of genes that are significantly differentially expressed in a specific cancer. From the Home page, clicking on the Search by cancer type tab allows users to select the cancer type of interest and then retrieve genes which are either over-or under-expressed. For example, selecting lung adenocarcinoma and the default settings (over-expressed; adjusted P value and P > 0.1) retrieves the 2089 genes, out of which the top expressed gene is FAM83A (Protein FAM83A; also called Tumor antigen BJ-TSA-9). It is interesting to note that FAM83A is considered a promising tumor biomarker of lung cancer (41). Similarly, the second highly expressed gene GREM1 (Gremlin) is also known to be over-expressed in lung cancer (46).

Pan-cancer analysis
The ability to sort, filter and further analyze the gene expression data collected in BioXpress allows users to compare and contrast expression of genes across many patients and cancer types. In addition to listing the genes that are significantly differentially expressed in multiple cancers (as described in the previous paragraph), Figure 3 provides an overview of the types of analysis that users can perform A B C using the downloaded data. Figure 3A heatmap and clustering were performed based on the percent of patients who have significantly differentially expressed genes. Clustering of samples or datasets across multiple cancer types, known as one type of pan-cancer analysis, is widely conducted by the community, especially by TCGA Research Network (47,48), and is of great interest from the aspect of personalized and translational medicine. To select genes that have strong association with transcriptomic changes of tumors, we picked the top 50 genes that are differentially expressed in the highest percent of samples. The darker colors in the figure show that several cancer types have genes which are differentially expressed in a majority of the patients (red boxes). The clustering based on the heatmap indicates that several cancer types have similar patterns [kidney renal clear cell carcinoma (KIRC) and kidney renal papillary cell carcinoma (KIRP); head and neck squamous cell carcinoma (HNSC) and stomach adenocarcinoma (STAD); lung squamous cell carcinoma (LUSC) and pancreatic adenocarcinoma (PAAD); thyroid carcinoma (THCA) and lung adenocarcinoma (LUAD)]. Figure 3B shows analysis results of expression data where no normal samples are available. The figure provides a view of cancer types that cluster together based on gene expression from cancer samples only. On the basis of the color distribution, it can be seen that several cancers have similar expression patterns and hence cluster together: breast cancer (BRCA) and lymphoma (Lymph); ovarian cancer (OV) and endometrial cancer (Endca); close to them are endocrine pancreas cancer (PAEN), prostate adenocarcinoma (PRAD), lung squamous cell carcinoma (LUSC), leukemia (Leuke) and brain cancer (Braca); KIRC and THCA; colon adenocarcinoma (COAD), PAAD and rectum adenocarcinoma (READ) are also clustered. Liver cancer (Livca) shows a distinct gene expression profile with all other cancer types listed based on the selected genes. Figure 3C provides a view of tissues which have similar expression patterns. Collection of expression data from multiple cancers as presented in supplementary Table S1 allow us to identify genes that are differentially expressed in more than one cancer type. For example, from this table we can see that nine genes are differentially expressed in all cancer types ( Table 2). It is important to note that in this particular case we do not consider the number of patients who have these genes over-or under-expressed. Therefore, each gene and its expression in a cancer type needs to be carefully evaluated on a case-by-case basis if one is interested in identifying genes which are differentially expressed in majority of the patients (please see examples in the next paragraph). It is interesting to note that five of the nine proteins are glycoproteins, two are phosphoproteins, six of them are secreted and seven are involved in biological process regulation (based on UniProtKB keyword and Gene Ontology annotation). This type of filtering and sorting can reveal ideal candidates for further evaluation as diagnostic or therapeutic targets. Furthermore, literature evidence reveals that eight of the 9 genes in Table 2 are genes known to be associated with cancer. For example, the first gene listed in Table 2, CCL21, participates in leukocytes and cancer cell migration through the CCR7/CCL19 (CCL21) axis to promote the growth and metastasis of various tumors such as breast cancer, melanoma, non-small cell lung cancer, head and neck, gastrointestinal and hematologic cancer (49). Second, c-glutamyltransferase is involved in cellular glutathione homeostasis, its expression is often significantly increased in human tumors and its role in tumor progression, invasion and drug resistance has been repeatedly suggested (50). Third, alterations in the ubiquitin system have direct or indirect roles in the genesis of various tumors due to defects in the ubiquitin-dependent proteolysis of critical house-keeping genes or cell-cycle elements-p53 is a good example (51). The next genes, Matrilysin (MMP7), are frequently over-expressed in human cancer tissues and are associated with cancer progression (52) and NCAM1 has been demonstrated to be one of the immunohistochemical markers for lung neuroendocrine tumors diagnosis (53), its expression level is upregulated in large cell lung tumor cell line H460-M (54). CHRDL1 is down-regulated (79-89% of 19) in follicular thyroid carcinoma (55) and the gene, WFDC2 (HE4), contains dispersed evidence: it has been demonstrated to be a biomarker for ovarian carcinoma (56) and it is known to be over-expressed in a range of different cell lines including ovarian, renal, lung, colon and breast lines, and cancers such as endometrial adenocarcinomas (57,58) and lung adenocarcinoma (59). The next gene, LCN2, has a wide range of functions in different types of cancers (thyroid, pancreatic, breast and colon cancer), and it is a potential diagnostic and prognostic marker in both benign and malignant human diseases (60). Finally, KRT80 and its role in cancer is not well studied although there is some evidence that this gene is differentially expressed in certain types of cancer (61,62). In addition to this list, a separate, precomputed table which lists all genes and their normalized expression values in tumors across all cancer types is also provided for download. This table can be used to identify genes which have, e.g. high variability in expression in certain cancers or low variability (possible house-keeping genes).
As mentioned above, one of the key questions in pancancer analysis of gene expression is-are there any genes which are significantly over-or under-expressed in multiple cancers in a large number of the patients. Supplementary Tables S3 and S4 provide the list of genes that are significantly differentially expressed in greater than 30% and 50% of the patients. Table 3 lists the top 5 genes (sorted based on the number of cancer types it is differentially expressed in) that are significantly over-and under-expressed in more than 50% of the patients. The first gene COL11A1 is known to be over-expressed in various epithelial cancers and is prominently correlated with invasion and metastasis (63). Its over-expression is associated with colorectal cancer (64), non-small cell lung cancer (65) and several other cancers (66). The next gene MMP11 over-expression is correlated with the aggression and invasion status of various types of carcinoma and is almost absent in normal adult organs and can be considered as a biomarker for diagnosis and prognosis (67,68). TMPRSS4 is highly expressed in pancreatic, colon, lung and gastric cancers, and is also expressed in a wide range of human cancer cell lines and has been demonstrated to facilitate the invasion, migration and metastasis of tumor cells (69,70). MMP1 is highly expressed in gastric carcinoma, breast cancer, lung and other cancers (71)(72)(73)(74)(75)(76)(77)(78). ADH1B is the first gene in the table that is known to be under-expressed in multiple cancers such as oral tongue squamous cell carcinoma (79) and intrahepatic cholangiocarcinoma (80). MT1H is under-expressed in adenoid cystic carcinoma of salivary gland, prostate and liver cancer due to hypermethylation of its promoter (81,82). In the next gene MT1G, the promoter is hypermethylated which results in its down-regulation in hepatoblastoma and prostate cancer (83,84). CHRDL1 interestingly is under-expressed in colorectal cancer (85) while over-expressed in pancreatic cancer (86) and for CA4 there is currently no publication associated with expression of these gene in cancers. We believe that filtering and sorting of data in BioXpress will help researchers to focus on expression profiles of genes which currently have very little published information. Another gene SFRP1 which is also found to be under-expressed in our dataset in five cancers (>50% of the patients) is known to be under-expressed in nine cancer types: cancers of the kidney, stomach, small intestine, pancreas, parathyroid, adrenal gland, gall bladder, endometrium, renal cell carcinoma and testis (87).

Downloadable files
Websites are ideal for performing gene and cancer-centric searches as described above. Some users may wish to perform large-scale analysis or filter the data based on additional parameters. To accommodate such users, all data can be downloaded in tab-delimited format. Additionally, a table of significantly under-or over-expressed genes in one or more patients is provided that has the following columns: gene name, UniProtKB accession, protein name, cancer types where the gene is expressed and count of the number of cancer types (supplementary Table S1). This table can be used to quickly identify genes that are differentially expressed in multiple cancer types in one or more patients. Additional downloads include PubMed Identifiers (PMIDs) and accessions that were manually curated (supplementary Table S2) and all data associated with differential and tumor-only expression. Future plans include addition of additional tables based on user requests.

Future Perspective
BioXpress will be updated every 6 months and detailed statistics for each release will be provided. Such statistics will allow users to track changes in the database over time. We will also integrate BioXpress in the High-performance Integrated Virtual Environment (HIVE) NGS and proteomics analysis platform. This integration will allow users to upload RNA-seq data, map reads to the reference genome using HIVE Hexagon (88), perform expression analysis and directly compare results with those available from BioXpress. As proteomic data become available for different cancer types through programs similar to the Clinical Proteomic Tumor Analysis Consortium (CPTAC) (89), we will map such data to the genes. We also plan to augment both data and function based on input from our users. Some potential new features include the following: addition of cancer subtypes; linking BioXpress to BioMuta (33) to obtain comprehensive view of expression as it may relate to mutation; integration of clinical annotations; inclusion of additional graphical elements and more. Our preliminary results show that there is a correlation between mutation density of a gene and its expression in certain types of cancer. We intend to explore this further in our future studies.

Supplementary Data
Supplementary data are available at Database Online.