Possible pathways used to predict different stages of lung adenocarcinoma.

We aimed to find some specific pathways that can be used to predict the stage of lung adenocarcinoma.RNA-Seq expression profile data and clinical data of lung adenocarcinoma (stage I [37], stage II 161], stage III [75], and stage IV [45]) were obtained from the TCGA dataset. The differentially expressed genes were merged, correlation coefficient matrix between genes was constructed with correlation analysis, and unsupervised clustering was carried out with hierarchical clustering method. The specific coexpression network in every stage was constructed with cytoscape software. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was performed with KOBAS database and Fisher exact test. Euclidean distance algorithm was used to calculate total deviation score. The diagnostic model was constructed with SVM algorithm.Eighteen specific genes were obtained by getting intersection of 4 group differentially expressed genes. Ten significantly enriched pathways were obtained. In the distribution map of 10 pathways score in different groups, degrees that sample groups deviated from the normal level were as follows: stage I < stage II < stage III < stage IV. The pathway score of 4 stages exhibited linear change in some pathways, and the score of 1 or 2 stages were significantly different from the rest stages in some pathways. There was significant difference between dead and alive for these pathways except thyroid hormone signaling pathway.Those 10 pathways are associated with the development of lung adenocarcinoma and may be able to predict different stages of it. Furthermore, these pathways except thyroid hormone signaling pathway may be able to predict the prognosis.

Introduction

Lung adenocarcinoma, the most common subtype of nonsmall cell lung cancer, is the leading cause of cancer-related death worldwide.[ It is responsible for more than 500,000 deaths annually worldwide.[ The all-stage 5-year survival rate for lung cancer is only 15% despite of some improvements in therapeutic method.[ It is reported that tumor, node, metastasis stage is the most important factor to predict the prognosis in lung cancer, including lung adenocarcinoma.[ In spite of the complete resection by lymph node dissection, the survival rates for the lung adenocarcinoma patients were 73% (stage IA) and 58% (stage IB).[ Thus, it is necessary and important to identify patients with lung adenocarcinoma in early stage, and then implement appropriate treatments.

Understanding of the molecular mechanisms underlying lung adenocarcinoma may provide good options for the identification of patients with lung adenocarcinoma. Takeuchi et al[ indicated that transcription termination factor, RNA polymerase I was a marker for lung adenocarcinoma. Yoshizawa et al[ believed that molecular testing for epidermal growth factor receptor and kirsten rat sarcoma viral oncogene homolog mutations was helpful for predicting prognosis among patients with resectable lung cancer. Furthermore, 1 study suggested that multiple therapies targeted against alterations in the RTK/RAS/RAF pathway advanced the treatment of lung adenocarcinoma.[ Inhibition of the Notch pathway may be an appropriate target for the treatment of lung adenocarcinoma.[ Besides, some drugs play roles in lung adenocarcinoma via some specific pathways. Such as, curcumin can induce the autophagy of lung adenocarcinoma via AMP-activated protein kinase signaling pathway, and angiotensin-(1–7) inhibits the migration and invasion of A549 human lung adenocarcinoma cells via inactivation of the phosphatidylinositol 3-kinase/Akt and mitogen-activated protein kinase signaling pathways.[ Some important pathways associated with lung adenocarcinoma have been researched, but that is not enough to clarify the molecular mechanisms of this disease.

In our present study, we used RNA-Seq expression profile data and clinical data of lung adenocarcinoma to perform gene correlation analysis, coexpression network analysis, and pathway enrichment analysis. Furthermore, classification model was constructed. We aimed to find some specific pathways that can be used to predict the stage of lung adenocarcinoma.

Materials and methods

Data remodeling and grouping

RNA-Seq expression profile data (LUAD_rnaseqv2__illuminahiseq_rnaseqv2.txt) and clinical data (nationwidechildrens.org_clinical_patient_luad.txt) of lung adenocarcinoma were obtained from the TCGA dataset. The platform for expression profile data was illuminahiseq_rnaseqv2 Level 3. The expression profile data samples were grouped into stage I, stage II, stage III, and stage IV based on clinical information. The sample size was 37, 161, 75, and 45 for stage I, stage II, stage III, and stage IV, respectively.

Our present study was a secondary analysis of RNA-Seq expression profile data downloaded from TCGA database. Thus, the ethical approval was not necessary.

Data normalization

The stage I was used as control group, and the mean and standard deviation of each gene in the control group were calculated. Then, Z-score normalization[ for all samples was performed, and the expressions of genes in the samples were subject to standardized normal distribution with mean value = 0 and the variance = 1.

Gene correlation analysis

We performed difference test for any 2 different stage lung adenocarcinoma samples by using Student t test, and the significant P value is set to .01. The genes differentially expressed compared with other stage were set as specific genes of this stage. The interaction between genes was changed in the process of lung adenocarcinoma development from one stage to next stage. The 2 interrelated genes have common biological effect in a nondisease state, cooperate with each other in biological terms, and show coexpression relationship on the expression level.[ In the disease state, the function of genes is abnormal, and subsequently this coexpression relationship will be changed. Thus, we examined the expression correlation of genes in every disease state with Pearson correlation coefficient. It was considered as positive correlation when correlation coefficient above 0.5 and negative correlation when correlation coefficient below −0.5.

Then, all the differentially expressed genes were merged, correlation coefficient matrix between genes was constructed with correlation analysis, and unsupervised clustering of samples and genes was carried out with hierarchical clustering method (hclust function in R package).[

Coexpression network analysis

The specific network in every stage was constructed based on the coexpression relationship between genes with cytoscape software, and there were edges between genes if coexpression relationship existed between 2 genes. We analyzed coexpression network by average shortest path length (ASLP), betweenness coefficient (BC), clustering coefficient (CC), degree, eccentricity (EC), and topological coefficient (TC). ASLP measured the distance between any 2 nodes in the network, and the mean value of distance set was obtained. The shorter the ASLP was, more closer the network was, and the higher the efficiency was. BC was used to assess the connected property of other gene pass through 1 gene. The higher the BC was the more obvious the roles of this gene as bridge. Degree was the number of edges connected to a node. The higher the degree was, the more obvious the proximity relationship of nodes was, and the higher the importance was. EC was the degree of node from the network core. The higher the EC was, the looser the network was. TC assessed the modularity of the network. Higher TC indicated more functional specificity, and lower TC indicated deficiency of functional specificity. If the edge of network was missing (the coexpression relationship was lost), ASLP, BC, CC, degree, EC, and TC would decrease, and the transmission efficiency of network signal would decrease.

Pathway enrichment analysis

Students t test was performed for any 2 groups, and genes differentially expressed in any 2 groups were regarded as specific gene of this group (P = .01). We merged 4 group specific genes and took the union to carry out the functional analysis. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was performed with KOBAS database and Fisher exact test.[

Functional pathway score

The genes that are functionally related to each other were enriched in the same pathway. We used Euclidean distance algorithm[ to calculate total deviation score.

The P pathway contains n enrichment genes, mean is expression mean of any gene (Gi) in the stage I group. Higher score (P) suggests that the pathway is obviously deviate from the normal level; lower score (P) suggests that the pathway is close to the normal level of stage I group.

Construction of classification model

Support vector machines, a supervised machine learning algorithm, can calculate the geometric hyperplane by combining characteristics in multidimensional data, and achieve the greatest degree of classification for 2 groups.[ The diagnostic model was constructed with support vector machines algorithm, and pathway deviation score of every sample was regarded as characteristic value. The stage I and stage II were combined into early benign group, and stage IV and stage III were combined into advanced malignant group, and the risks were predicted for 2 combined samples. Parameter is linear kernel, penalty coefficient = 1, and gamma = 0. Receiver-operating characteristic (ROC) curve was drawn and classification efficiency of model was assessed with precision, recall, f1-score, and 5-fold cross-validation.

Correlation between the identified 10 pathway and patients’ prognosis

In order to further verify the identified 10 pathway and predict patients’ prognosis, we compared the correlation between the pathway and patients’ prognosis.

Results

Totally 18 specific genes were obtained by getting intersection of 4 group differentially expressed genes. The heat map for stage I, stage II, stage III, and stage IV (supplement Fig. 1) showed that there were significant correlations between the 18 genes and formed significant correlation module. These positively related genes or negatively related genes were not completely consistent in different stage, which suggests that the expression correlation between genes was changed with the development of lung adenocarcinoma. The change of expression correlation might be the reason of functional difference.

The results of hierarchical clustering (supplement Fig. 2) showed that almost all of the stage III and stage IV samples were clustered together, and stage I and stage II samples were clustered together. Thus, we concluded that there was significant difference among sample groups of different risks at the molecular level, and early benign group (stage I–stage II) samples tended to have higher similarity and advanced malignant group (stage III–stage IV) samples tended to have higher similarity.

The coexpression networks for 4 stages are shown in Fig. 1A–D. The results of network analysis are shown in Fig. 2. From stage I to stage IV, ASLP gradually increased, CC and degree gradually decreased, and there was no significant difference for BC, EC, and TC. The increase of ASLP indicated that in the whole network the average distance between any 2 nodes was increased, the network connectivity was decreased, and the signal transmission efficiency was reduced. The decrease of CC and degree indicated that the degree of network dispersion was increased, and the structure was loose.

Figure 1

The coexpression networks for 4 stages (A: stage I, B: stage II, C: stage III, and D: stage IV).

Figure 2

The results of network analysis.

Totally 10 significantly enriched pathways are obtained (Table 1). It showed that genes associated with malignant degree of lung adenocarcinoma were enriched in mineral absorption, aldosterone-regulated sodium reabsorption, adipocytokine signaling pathway, and thyroid hormone signaling pathway, in addition to several cancer pathways (small cell lung cancer, endometrial cancer, nonsmall cell lung cancer, and pancreatic cancer). It indicated that lack of trace elements, sodium ion concentration, adipocyte activity, as well as changes in hormone levels were important factors affecting different malignant progression of lung adenocarcinoma.

Table 1

Significantly enriched pathways by differentially expressed genes.

Figure 3 is the distribution map of 10 pathways score in different groups. From the deviation level, degrees that sample groups deviated from the normal level were as follows: stage I < stage II < stage III < stage IV.

Figure 3

The distribution map of 10 pathways score in different groups.

The boxplot of 10 pathway score in different stages (Fig. 4) showed that the score of 4 stages exhibited linear change in some pathways, and the score of 1 or 2 stages were significantly different from the rest stages in some pathways. It demonstrated that these pathways showed different functional level in different stage, and these pathways could be used to classify lung adenocarcinoma samples of different malignancy degrees.

Figure 4

The boxplot of 10 pathway score in different stages.

The ROC curve is shown in Fig. 5. The mean ROC was 0.84. Lung adenocarcinoma with different disease stage could be distinguished effectively by using stage specific pathway of lung adenocarcinoma as characteristic. Using stage specific pathway of lung adenocarcinoma as characteristic had more advantages compared with serological markers or genes as characteristic, such as pathways could show synergistic action of multiple associated genes systematically and had more robustness. Classification reports for early stage and late stage are shown in Table 2. The average precision and recall for early stage was both higher than that for late stage, but the prediction effect for early stage alone was not as good as that for late stage. It suggested that our model had higher sensitivity for predicting patients with advanced lung adenocarcinoma. At the same time, it explained again that these10 pathways had significantly different functional levels in different lung adenocarcinoma stage. Furthermore, these pathways might provide new ideas for explaining risk and prognosis of lung adenocarcinoma at different stages and realization of personalized therapy.

Figure 5

The receiver-operating characteristic (ROC) curve.

Table 2

Classification reports for early stage (A) and late stage (B).

We analyzed the correlation between the identified 10 pathway and patients’ prognosis. There was significant difference between dead samples and alive samples for small cell lung cancer, PPAR signaling pathway, mineral absorption, endometrial cancer, nonsmall cell lung cancer, hepatitis C, pancreatic cancer, aldosterone-regulated sodium reabsorption, and adipocytokine-signaling pathway (P < .05, Fig. 6). No significant difference was found for thyroid hormone signaling pathway (P > .05).

Figure 6

Correlation between the identified pathway and patients’ prognosis. The horizontal ordinate, patients’ prognosis; longitudinal coordinates, pathscore; read sample, dead; and green sample, alive.

Discussions

Lung adenocarcinoma is a highly malignant cancer, and its cells are involved in pulmonary circulation system. The cells can be released into the circulatory system and grow in other tissues and organs, and then the spread and metastasis of cancer was achieved. Thus, it is necessary to identify the stage of lung adenocarcinoma and then block it before the malignant progression. In this study, from stage I to stage IV, ASLP gradually increased, CC and degree gradually decreased, and there was no significant difference for BC, EC, and TC in coexpression network. Totally 10 significantly enriched pathways were obtained. In the distribution map of 10 pathways score in different groups, degrees that sample groups deviated from the normal level were as follows: stage I < stage II < stage III < stage IV. The boxplot of 10 pathway score in different stages showed that the score of 4 stages exhibited linear change in some pathways, and the score of 1 or 2 stages were significantly different from the rest stages in some pathways.

In the present study, from stage I to stage IV, ASLP gradually increased, CC and degree gradually decreased, and there was no significant difference for BC, EC, and TC in coexpression network. These indicators suggested that the efficiency of the system network was gradually reduced in the progress from stage I to stage IV. In the state of disease, the decrease of network efficiency and signal transmission efficiency led to the reduction of some important function, particularly functions related with disease response or self-activated immune system. Under the stimulation of the disease, due to the decrease of the signal transmission efficiency, functional mechanism of stress and repair was suppressed, resulting in the disease progression and the increase of the malignancy degree.

Furthermore, totally 10 significantly enriched pathways were obtained from pathway enrichment analysis. They were cancer-related pathways (small cell lung cancer, endometrial cancer, nonsmall cell lung cancer, and pancreatic cancer), PPAR-signaling pathway, mineral absorption, aldosterone-regulated sodium reabsorption, adipocytokine signaling pathway, hepatitis C, and thyroid hormone signaling pathway. One study showed that the expression status of PPAR-γ correlated with differentiation status and survival in the lung cancer patients, and many PPAR-γ ligands (exert their effects through both PPAR-γ dependent and independent pathways) could inhibit tumor growth and progression in lung cancer preclinical models by modulating various cellular processes in cancer cells.[ The mineral absorption was also a significantly downregulated pathway in the article of Wang et al[ studying lung adenocarcinoma by using microarray analysis. Mijatovic et al[ indicated that the α1 subunit of the sodium pump could be a target to treat nonsmall cell lung cancers, and thus sodium reabsorption might also be important to lung adenocarcinoma. Adipocytokines play a central role in many aspects of inflammation and immunity,[ and the hepatitis C virus is also associated with chronic inflammation.[ Furthermore, lung adenocarcinoma is associated with inflammation.[ The metastatic thyroid cancer is originate from lung adenocarcinoma,[ and thyroid transcription factor-1 is almost exclusively expressed in the normal thyroid and lung, as well as in carcinomas derived from these organs.[ Therefore, our present study is in line with the former researches and suggests that these 10 pathways may play significant roles in lung adenocarcinoma.

In addition, in the distribution map of 10 pathways score in different groups, degrees that sample groups deviated from the normal level were as follows: stage I < stage II < stage III < stage IV. In the boxplot of 10 pathway score in different stages, the score of 4 stages exhibited linear change in small cell cancer pathway, mineral absorption pathway, endometrial cancer pathway, hepatitis C pathway, pancreatic cancer pathway, and aldosterone-regulated sodium reabsorption pathway, and the score of stage I was almost equal to stage II, stage III equal to stage IV for nonsmall cell lung cancer pathway and PPAR-signaling pathway. However, for the adipocytokine-signaling pathway, the score of stage I and II developed in the forward direction, then the score of stage III and IV developed toward negative direction. The score of stage I, II, and III developed in the forward direction, then the score of stage IV developed toward negative direction for thyroid hormone signaling pathway. These pathways showed significantly different functional levels in different process of lung adenocarcinoma development. On the one hand, with development of lung adenocarcinoma, those pathways exhibited linear change, nonsmall cell lung cancer pathway and PPAR signaling pathway developed linear or periodical change of function, and the functions gradually strengthen or gradually weaken. On the other hand, adipocytokine-signaling pathway and thyroid hormone signaling pathway developed functional change in the opposite direction at a certain stage with development of lung adenocarcinoma, indicating the existence of a critical point of functional variation in these 2 pathways. Furthermore, the mean ROC of the classification model was 0.84. Therefore, these 10 pathways have important roles for predicting lung adenocarcinoma in different stages, and these pathways may provide new ideas for explaining the pathogenic mechanism of lung adenocarcinoma and new treatment. But concrete prediction methods need further studies.

Besides, in our present study, the analysis for the correlation between the identified pathway and patients’ prognosis showed that there was significant difference between dead samples and alive samples for these pathways except thyroid hormone signaling pathway. Thus, these 9 pathways may be able to predict the patients’ prognosis.

In conclusion, cancer-related pathways (small cell lung cancer, endometrial cancer, nonsmall cell lung cancer, and pancreatic cancer), PPAR-signaling pathway, mineral absorption, aldosterone-regulated sodium reabsorption, adipocytokine-signaling pathway, hepatitis C, and thyroid hormone signaling pathway are associated with the development of lung adenocarcinoma, and may be able to predict different stages of lung adenocarcinoma. Furthermore, these pathways except thyroid hormone signaling pathway may be able to predict the patients’ prognosis. Although these pathways can be used to predict different stages of lung adenocarcinoma, whether these pathways are also applied to squamous cell cancer or large cell cancer need further research in the future. However, no experimental verification is the limitation of this study, and further studies are needed to verify our results. Besides, for further exact analysis of this topic, the analyses of the data with exact subgroup of the cancer in the TCGA datasheet are needed in our next research plan.