Comparative gene co-expression network analysis of epithelial to mesenchymal transition reveals lung cancer progression stages

Abstract

Background: The epithelial to mesenchymal transition (EMT) plays a key role in lung cancer progression and drug resistance. The dynamics and stability of gene expression patterns as cancer cells transition from E to M at a systems level and relevance to patient outcomes are unknown.

Methods: Using comparative network and clustering analysis, we systematically analyzed time-series gene expression data from lung cancer cell lines H358 and A549 that were induced to undergo EMT. We also predicted the putative regulatory networks controlling EMT expression dynamics, especially for the EMT-dynamic genes and related these patterns to patient outcomes using data from TCGA. Example EMT hub regulatory genes were validated using RNAi.

Results: We identified several novel genes distinct from the static states of E or M that exhibited temporal expression patterns or 'periods' during the EMT process that were shared in different lung cancer cell lines. For example, cell cycle and metabolic genes were found to be similarly down-regulated where immune-associated genes were up-regulated after middle EMT stages. The presence of EMT-dynamic gene expression patterns supports the presence of differential activation and repression timings at the transcriptional level for various pathways and functions during EMT that are not detected in pure E or M cells. Importantly, the cell line identified EMT-dynamic genes were found to be present in lung cancer patient tissues and associated with patient outcomes.

Conclusions: Our study suggests that in vitro identified EMT-dynamic genes capture elements of gene EMT expression dynamics at the patient level. Measurement of EMT dynamic genes, as opposed to E or M only, is potentially useful in future efforts aimed at classifying patient's responses to treatments based on the EMT dynamics in the tissue.

Keywords: Cancer progression; Comparative network analysis; EMT-dynamic signature genes; Epithelial to mesenchymal transition; Gene regulatory network; Lung cancer.

Fig. 1

Previously identified EMT signature genes have distinct temporal expression dynamics during epithelial to mesenchymal transition in lung cancer. a The heatmaps show the normalized gene expression levels of 76 known EMT genes across H358’s ten EMT stages (left, 0 h, 1 h, 2 h, 4 h, 6 h, 8 h, 16 h, 24 h, 72 h, 168 h) and A549’s eight EMT stages (right, 0 h, 6 h, 12 h, 24 h, 36 h, 48 h, 72, 96 h) [16, 17]. These EMT genes were predicted according to their fold changes between epithelial and mesenchymal states only. Red: highly expressed. Green: lowly expressed. b PCA of 76 known EMT genes using their gene expression data in H358 EMT. The dots are genes. The x-axis is the PC1 coefficient, and the y-axis is the PC2 coefficient. The four gene groups have been clustered by K-means. The embedded boxplots display the gene expression level distributions across H358 EMT stages for four groups. The cyan group represents genes with an increasing expression pattern at middle EMT stages (~72 h and continuing) that includes the EMT associated EGFR resistance oncogene AXL [10]. The red group consists of EMT genes including TGFB1 having an increasing expression pattern at ~ 16 h which decays after 168 h. The gene expression in the green group increases slowly from 16 h but dramatically decreases after 168 h. The blue group includes genes that are decreasing in expression during EMT (from 24 h on)

Fig. 2

The eigengenes and enriched pathways of gene co-expression modules in H358 EMT. a The heatmap shows the eigengene expression levels across H358 EMT stages for 55 gene co-expression modules. Red: high expression level; Green: low expression level. These eigengenes represent the gene expression dynamic patterns at the system level in H358 EMT. The gene co-expression modules are identified using WGCNA [18]. b The enriched KEGG pathways of gene co-expression modules, which are found by clusterProfiler [20]. The rows are the enriched pathways, and the columns are modules. The dot size is proportional to the modular gene fraction involved in the pathway (i.e., number of pathway genes in the module over number of total pathway genes). The darkness of color is proportional to the enrichment score (adjusted p value)

Fig. 3

The gene coefficients over module eigenegenes, enriched functions and pathways of EMT-dynamic genes. a The heatmap shows the correlation matrix between ~18,000 genes and the eigengenes of 55 modules in H358 EMT. b The barplot shows the enriched pathways and functions (y-axis) of 254 EMT-dynamic genes. The x-axis is the number of EMT-dynamic genes for the enriched pathway/function. The colors correspond to the –log10(enrichment p value)

Fig. 4

The gene expression dynamics and regulatory networks for cell cycle modules in H358 EMT. a The heatmap shows the eigengene expression level across H358 EMT stages for all 7 cell cycle modules. Red: high expression level; Blue: low expression level. b The predicted gene regulatory network controlling the cell cycle modules. Nodes are the transcription factors (TFs). The TFs in the network have significantly large numbers of target genes in the cell cycle modules (p < 0.05). The orange TFs have highly positive correlated expression with cell cycle eigengenes (Pearson correlation coefficient > 0.7), and the light-blue TFs have negatively correlation (Pearson correlation coefficient < −0.7). c The gene expression fold changes by RNAi depletion in M cells relative to E cells. The dashed line highlights the 1-fold of down-regulation

Fig. 5

EMT-dynamic genes reveal EMT periods of TCGA lung cancer patients. a The PCA plots show the PC1 and PC2 coefficients of 522 TCGA LUAD patients using their TCGA gene expression data for three gene sets: 254 EMT-dynamic genes (left), 76 previously identified ETM genes (top right) and all ~18,000 genes (top bottom). The dots are patients. The colors denote the patients personalized EMT periods (PEPs) in H358 (Methods): green (h0-8), brown (h16-24), purple (h72), magenta (h168) and yellow (>h500). b To evaluate how well TCGA LUAD patients are clustered on Panel A using three gene sets, the boxplots show the silhouette value distributions. The 254 EMT-dynamic genes have significantly higher silhouette values than others. c The personalized EMT periods (PEPs) of 522 TCGA LUAD patients significantly classify the patient survival rates by the Kaplan–Meier analysis. Three major PEP groups are h16 (red), h168 (green) and h72 (blue). The reference group of hazard ratios consists of the patients whose PEPs is h0-8