2.2 Gene expression analysis

In order to define expression-based concepts, we developed a gene expression analysis pipeline that uses a carefully chosen, statistical method for each step. The gene expression concept type is populated with human Affymetrix experiments in GEO. Details of the analysis pipeline are provided in the Supplementary Methods section, and outlined briefly here. The pipeline downloads the raw data, pre-processes it with a Entrez ID centered CDF package (Dai et al., 2005) and normalizes it, outputs quality control data, and tests for differentially expressed genes using an empirical Bayes method (Sartor et al., 2006). The tests are set up manually through an interface, and following testing, gene sets (concepts) are defined by the top ranked genes.

2.3 Enrichment testing

For public concepts, all pairs of concepts from all concept types were tested for whether there exists a larger number of overlapping genes than is expected by chance. We use a slightly modified Fisher's exact test, identical to the ‘Ease score’, which helps to stabilize results by penalizing tests with small numbers of overlap (Hosack et al., 2003). P-values are adjusted for multiple testing by calculating q-values using the FDR method (Benjamini and Hochberg, 1995). The default display is those concepts with q-value <0.05, but the user may choose a different q- or p-value cutoff.

Gene lists uploaded by the user are considered concepts in their private, Experimental concept type. The gene lists are converted to human Entrez Gene IDs, if necessary, and stored and tested in a private concept type. Users have the option of converting from a list of mouse or rat genes using NCBI Homologene homolog families with our conversion tool, or from a set of compounds/metabolites. Users also have the option to upload a background gene set consisting of all the genes that were interrogated in their study (e.g. all genes on a microarray platform). If no background set is provided, ConceptGen uses all Entrez IDs as default. The modified Fisher's exact test described above is then implemented, and q-values are calculated for the experimental list within each concept type separately. Users can filter results based on concept type(s) and/or significance levels (p- or q-values) for export in spreadsheet format or visualizations.

A background set is defined as all genes (Entrez IDs) that were interrogated in creating the concept type. For example, for GO concept types, the background set is all genes that are assigned to at least one ontology term. It is important to use the correct background gene set for each enrichment test, and for that we use the intersection of the background gene sets for the two concept types of the concepts being tested. Thus, for example, if we are testing a GO term versus a microRNA target list, we use all genes that are in both the GO background set and the microRNA target background set.

2.4 Gene set relation mapping—graph network visualization

Users can explore concepts from any of the public sources, or load their own sets of genes to define private concepts. Once a concept is selected, the concepts that are paired with it (those whose enrichment scores are significant) can be selected by category, significance or individually, to participate in a concept-to-concept graph (see the figs in Section 3.2), with nodes representing concepts, and edges representing significant enrichment of overlapping genes. The concept networks, nodes and edges are displayed in an interactive web application (coded in Adobe^® Flex/flash). The graph is laid out on the display using a standard layout from Adobe's open source code, which implements a force directed layout algorithm (Fruchterman and Reingold, 1991). This layout algorithm results in highly interconnected groups of concepts clustering together. Within that layout, the concept type of each concept node is shown by the color, the size of the concept node is based on the number of genes in the concept, and the thickness of the edge lines is based on the number of overlapping genes. These graphical network displays can be further explored to find the genes that concepts have in common, filter the graph based on sets of genes, display the statistics associated with a concept or edge, or explore protein interactions within a node or edge. By moving among the results panel, the graphics display, and the protein interaction networks, the user can narrow in on a set of concepts of interest and conceptualize results.

2.5 Gene set relation mapping—heatmap view

The heatmap view offers an alternative view to the network graph, convenient for visualizing large numbers of concepts. It also allows one to see at a glance which genes are responsible for the enrichment of which concepts, and which gene groups co-occur in the same concepts most often (see Supplementary Figures S8 and S9 for examples). The heatmap plots genes (columns) versus enriched concepts (rows), and values used are 0 when a gene does not belong to the concept, or the number of enriched concepts to which the gene belongs otherwise. Genes and concepts are clustered using the complete linkage hierarchical clustering method with the Euclidean distance measure. The color of columns ranges from black (gene does not belong to the enriched concept) to bright red (genes belonging to the most enriched concepts.) Users can toggle between the heatmap and graph network views, and can use a ‘draw tool’ to choose a cluster or section of the heatmap to explore in the graph network display.

3.1 Properties of the data in ConceptGen

The general aspects involved in development of a gene set enrichment or gene set relation mapping tool are: (i) the test used to assess significance of enrichment, (ii) the concepts and concept types used, (iii) usability of the software tool, (iv) visualizations offered by the tool, (v) quality of statistical methods used to build concepts from experimental data (if applicable); and (vi) adjusting the significance levels for multiple comparisons. The quality of the software application will be a function of the above aspects. Thus in the development of ConceptGen, we have carefully chosen the implementation of each of these aspects to create an easily accessible, powerful, stream-lined and accurate enrichment testing and gene set relation mapping tool.

Table 1 lists the types of biological knowledge incorporated in ConceptGen and the names of the concept type(s) used to represent each of them. These concept types encompass pathways, molecular processes, chromosome location, biological targets, protein–protein interactions, diseases, and experimental and literature-derived data. Each biological knowledge type provides the ability to generate a different type of hypothesis; which concept types are of greatest interest will depend on the context of the application. For example, testing the protein interactions concept type will identify proteins that interact with a significant number of user-supplied gene products, possibly generating a hypothesis for what other proteins may be activated or play a central role in the process under study. For protein interactions, we used the Michigan molecular interactions (MiMI) database, which deep-merges several sources of interactions, resulting in a comprehensive database of human protein-protein interactions and thus greater power to detect significant enrichments than using any one data source alone. A notable member of our concept types is MeSH [derived from Medical Subject Headings (MeSH^®)]. Medical subject headings is the National Library of Medicine's controlled vocabulary built in a hierarchical structure, which covers subjects ranging from anatomy, drugs, and diseases to social phenomena, geographical locations and proteins. Therefore, the MeSH concept type in itself contains concepts relating to various types of biological knowledge. Use of the MeSH concept type will help users identify relationships with concepts in literature that may otherwise be overlooked.

Table 2 compares features of ConceptGen with four related software programs performing gene set enrichment and/or gene set relation mapping: DAVID/EASE (Dennis et al., 2003), GSEA and MSigDB (Subramanian et al., 2005), Oncomine's molecular concept mapping (Rhodes et al., 2007) and GeneDecks (Stelzer et al., 2005). Similar to other programs, such as DAVID, certain concept types tend to have larger or smaller concept sizes (number of genes populating each concept) than others, and different distributions (Supplementary Fig. S1). For example, gene sets based on gene expression experiments and microRNA targets had the largest concepts on average. Others, such as GO and MeSH have a broad range of concept sizes due to their hierarchical nature.

3.1.1 Inter-connectivity among concept types

To determine how closely-related the concept types are amongst themselves, we calculated a measure of connectivity. Connectivity is defined as s/n where s = No. of tests between concept types with q-value <0.05 and n = total number of possible connections. As seen in Figure 1, most although not all concept types have highest connectivity with self, and the pathway databases KEGG, Biocarta and Panther form a block of high inter-connectivity. In addition to these expected results, we observe that the pathways have high inter-connectivity with protein interactions and GO. Gene expression data has relatively high connectivity with KEGG Pathway, GO cellular component, GO biological process, MiMI, MeSH, and interestingly, microRNA targets (miRBase). We can speculate that the higher connectivity of gene expression data with KEGG compared with the other pathway databases may be an indicator of how well KEGG approximates the reality of biological systems, although more research is needed to test this prediction. The clustering divides the concept types into two groups, which can be interpreted as (i) those that assess the pathways, molecular processes and cellular networks, and (ii) all others.

Open in a separate window

Fig. 1.

Inter-connectivity among concept types. Connectivity is defined as s/n where s = No. of tests between concept types with q-value <0.05 and n = total number of possible connections. Values were clustered using complete linkage hierarchical clustering in R.

3.2 Application to TGFβ-induced EMT in lung adenocarcinoma cells

EMT is a highly conserved embryonic and developmental process that facilitates the dispersion of cells. During EMT, cells lose their epithelial properties, while acquiring mesenchymal properties which enable them to migrate to a predetermined destination in order to generate distinct tissue types (Kalluri and Weinberg, 2009). A similar process is reactivated in cancer cells as an early event during tumor metastasis and confers certain fundamental abilities essential for the process of metastasis. These include the ability to migrate, invade, resist apoptosis and evade immune surveillance (Kalluri and Weinberg, 2009). Multifunctional cytokine TGF-β is a potent inducer of EMT (Thiery and Sleeman, 2006). TGF-β-induced EMT in A549 lung adenocarcinoma cell line provides an excellent in vitro correlate for mechanistic investigations of EMT in the context of tumor progression (Keshamouni et al., 2006, 2009). Here we demonstrate the functionalities of ConceptGen in enabling visualization of data and contextualizing prior knowledge to generate new hypotheses, by utilizing time course gene expression data (Control versus 0.5, 1, 2, 4, 8, 16, 24 and 72 h) obtained from A549 cells undergoing TGF-β-induced EMT.

Data were processed and analyzed as previously described (Keshamouni et al., 2009). The data have been deposited in NCBI's GEO (Barrett et al., 2007) and are accessible as GSE17708. The set of differentially expressed genes at each time (defined by p-value <0.001 and >2-fold change) was uploaded to ConceptGen for analysis, and network graphs and heatmaps of results were visualized for each time point. All concept types were used for testing, although in other situations one may wish to limit testing to a subset of concept types. Supplementary Figure 4 illustrates the main results page of ConceptGen and Figure 2 displays the network graph for the 30 min time point. In this figure, colors represent different knowledge types. Results for each time point were then exported from ConceptGen in spreadsheet format and the enrichment significance profiles for concepts with q-value <0.05 for at least one time point were clustered for an overall visualization of enriched concept profiles. This allowed us to easily visualize which processes were being affected throughout the time course. Model-based clustering (Medvedovic and Sivaganesan, 2002) of the −log₁₀(p-values) for enrichment was performed separately for GO terms (Supplementary Fig. S5), Gene expression enrichment (Supplementary Fig. S6), and other concept types (Fig. 3) for ease of interpretation.

Open in a separate window

Fig. 2.

ConceptGen Network Graph after 30 min TGFβ-induction of A549 cells shows TGF-β receptors, SMAD proteins, Myod1 interactions, and transcription factor activity to be enriched. Shown are all concepts with q-value <0.05 after filtering for Gene Expression.

Open in a separate window

Fig. 3.

Heatmap of enriched Concept profiles throughout TGF-β-induced EMT transition in the A549 cell line for concept types other than GO and Gene Expression: rows are concepts from KEGG, Biocarta, and Panther Pathways, MiMI protein interactions, MeSH, Pfam and Transfac. All concepts with q-value <0.05 for ≥1 time point were clustered using −log₁₀(p-values) to create a ‘bird's eye view’ of which processes were turning on and off throughout the time course. For presentation, enriched concepts were summarized for each main cluster. Heatmaps for GO and gene expression concepts are presented in Supplementary Material.

TGF-β initiates signaling by binding to type II receptors on the cell surface, triggering the activation of type I receptors. Activated type I receptors phosphorylate SMAD proteins. Phosphorylated/activated SMADs move into the nucleus where they bind to transcriptional co-activators or co-repressors to regulate target gene expression (Massague, 2000). Enrichment analysis for 30 min to 6 h time points indicated effects of a robust transcriptional reprogramming with the modulation of several transcriptional factors (summarized in Fig. 3 and Supplementary Fig. S5), including SMADs, MyoD and induction of type I TGF-β receptors (Fig. 2), reflecting the primary mode of action of TGF-β. Consistent with growth inhibitory affects of TGF-β (Massague, 2000), we observed enrichment of concepts that are reflective of negative growth regulation in the middle time points (Fig. 3 and Supplementary Fig. S5). Similarly, around the same time points we also observed the enrichment of GO terms reflecting inhibition of apoptosis (Supplementary Fig. S5), consistent with the EMT induced apoptotic resistance demonstrated in earlier studies (Lee et al., 2006). Enrichment of concepts of cell movement, cell adhesion, extracellular matrix and matrix metalloproteases (Figs 3 and ​and4)4) are reflective of migratory and invasive phenotypes acquired during EMT (Keshamouni et al., 2006), whereas the concepts such as regulation of cell size, actin binding and cytoskeletal protein binding are consistent with the dramatic change in morphology and robust cytoskeletal reorganization that occurs during EMT (Thiery and Sleeman, 2006). On the whole, functional enrichment analysis of multiple biological knowledge types with ConceptGen coupled with a heatmap viewer provided a reliable overview of the entire data set, accurately representing the underlying biology across time points.

Open in a separate window

Fig. 4.

ConceptGen Network Graph after 8hrs of TGFβ induction in A549 cells shows extracellular matrix, cell adhesion, cell movement and metallothioneins are enriched, consistent with the migratory and invasive phenotypes acquired during EMT. Shown are the top significant concepts from Pfam, MeSH, and MiMI. GO biological processes indicated that vascular and blood vessel development were also enriched.

ConceptGen also allowed us to contextualize the data in hand with prior knowledge in the public domain which includes previously published gene expression data sets (Supplementary Fig. S5) and protein interaction maps. Enrichment analysis of EMT time course data with publicaly available gene expression data sets has returned interesting results. Similar to the analysis of other concept types, these results were reflective of TGF-β biology and the complexity of EMT process in a time depedent manner, with specific examples provided in supplemental material. Overlap with other expression sets allows identification of common sets of responses that in combination are not well represented in any other knowledge domain. For example, one could detect other conditions that resulted in a certin mix of DNA repair, response to stress and apoptosis. Using the MiMI NetBrowser tool in ConceptGen, we can also visualize a network of protein–protein interactions for the genes belonging to any node or edge, and can identify potential regulatory hubs at each time point. For example, at the 1 h time point, we identified JUN and TGFBR1 as hubs, consistent with the induction and activation of these two proteins respectively in response to TGF-β (Fig. 5).

Open in a separate window

Fig. 5.

Protein interactions among differentially expressed genes after 1hr of TGFβ induction in A549 cells, visualized by clicking on the gear icon for the uploaded concept in the network browser view, which links to the MiMI NetBrowser view of known protein interactions.

3.3 Application with compounds

A novel feature of ConceptGen is its ability to take one or more compounds as input (in the Upload Gene Set page) and find the genes encoding metabolic enzymes related to those compounds. As an example, we input six metabolites that were found to significantly change in LnCAP cells in the disease progression from benign to localized to metastatic prostate cancer (Sreekumar et al., 2009). This list of metabolites (glycerol-3-phosphate, kynurenine, leucine, proline, sarcosine and uracil) was generated by unbiased metabolomics experiments, and was linked to a total of 37 genes using the methods described for building the Metabolite concept type (see Supplemental Methods section). The first step in characterizing a set of metabolites is identification of the biochemical pathways to which they belong. ConceptGen offers a convenient way to do this, without having to manually map the compounds individually onto KEGG or a similar database.

The conclusions from this simple analysis are consistent with the findings reported in Sreekumar et al. (2009), namely that amino acid metabolism (q-value = 7.4 × 10⁻¹²), nitrogen breakdown (q = 2.3 × 10⁻⁴) and amino methyltransferase activity (q = 9.6 × 10⁻³) from GO are enriched (Supplementary Fig. S7). Upon inspection with the heatmap view, we saw that the above-mentioned concepts all clustered with the compound sarcosine (Supplementary Fig. S8), indicating that mainly enzymes related to sarcosine drove the enrichment of these concepts. Additional processes identified as enriched were aminoacyl tRNA biosynthesis KEGG pathway [recently identified as disregulated by androgen in prostate cancer (Vellaichamy et al., 2009)], and procollagen–proline dioxygenase activity from GO molecular function, which is involved in collagen biosynthesis and folding. This is interesting because collagen production is known to be involved in cancer progression, and because the prolyl hydroxylase gene family was recently noted as being a novel class of tumor suppressors, specifically in breast cancer (Shah et al., 2009).

We would like to thank all members of the National Center for Integrative Biomedical Informatics (NCIBI) who contributed to database support, gave us valuable input on improving the user interface, or other discussion, particularly Xiaosong Wang. We thank Jyoti Athanikar for a thorough reading of the manuscript.

Funding: This work was funded by the NIGMS (grant U54 DA021519-01A1) (NCIBI), NIH R01 (CA132571), American Cancer Society (CSM-116801), and MTTC GR687.

Conflict of interest: none declared