This article is mentioned in:

Introduction

Lung cancer is the type of cancer with the highest incidence rate worldwide (1). Lung adenocarcinoma is the most common histological type, accounting for ~50% of all lung cancer cases (2–4). The prognosis for lung adenocarcinoma is generally poor, since its histopathological sub-classification is unclear and it tends to metastasize widely at an early stage of the disease (5,6). Therefore, the majority of patients with metastatic adenocarcinoma receive empirical chemotherapy regimens (7). Different genetic or molecular backgrounds may affect the response of an individual to cancer treatment and finally lead to a survival discrepancy (8). Thus, improved understanding of the mechanisms underlying the occurrence and development of lung adenocarcinoma is needed, and prognostic factors that identify patients at high risk for recurrence or metastasis are urgently required.

Among the numerous pathogenic factors of lung adenocarcinoma, long non-coding RNAs (lncRNAs), which are ncRNAs of >200 nucleotides in length, have attracted great interest. Lung adenocarcinoma-associated lncRNAs are increasingly being identified. For example, colon cancer associated transcript 2, a lung adenocarcinoma-specific lncRNA, promotes the invasion of non-small cell lung cancer and may be a potential biomarker for lymph node metastasis (9). HNF1A antisense RNA 1 is overexpressed in lung adenocarcinoma compared with corresponding non-tumor tissues, and its expression level is significantly associated with Tumor-Node-Metastasis (TNM) stage, tumor size and lymph node metastasis, leading to worse overall survival (10). The expression of the lncRNAs DKFZP434 L187 and LOC285548 have been demonstrated to be negatively associated with the overall survival of patients with lung adenocarcinoma (3). A total of three differentially expressed lncRNAs [maternally expressed 3 (MEG3), myocardial infarction associated transcript (MIAT) and MIR4697 host gene) are reportedly associated with clinical features of lung adenocarcinoma, while nine differentially expressed lncRNAs [long intergenic non-protein coding RNA 115 (LINC00115), long intergenic non-protein coding RNA 265, long intergenic non-protein coding RNA 1001, long intergenic non-protein coding RNA 1002, MIR22 host gene, NFYC antisense RNA 1, small nucleolar RNA host gene 10, THUMPD3 antisense RNA 1 and TMPO antisense RNA 1] were revealed as prognostic biomarkers for lung adenocarcinoma (11). MEG3 and MIAT may be important in the development of lung adenocarcinoma via their interactions with microRNA (miR)-106, with the consequent regulation of mitogen-activated protein kinase (MAPK)9 in MAPK signaling pathways, while LINC00115 may interact with miR-7 to regulate fibroblast growth factor 2 and thus be involved in ‘pathways in cancer’ (11).

Despite the marked progress that has been made, the prognostic roles of lncRNAs in lung adenocarcinoma and the underlying mechanisms remain poorly characterized. Further functional studies are required to identify additional lung adenocarcinoma associated lncRNAs and verify their functional mechanisms in lung adenocarcinoma.

In the present study, extensive RNA sequencing (RNA-seq) data and clinical survival prognosis information from patients with lung adenocarcinoma was downloaded from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases at the National Center for Biotechnology Information (NCBI), in order to construct a co-expression network and mine network modules with particular biological functions. Our goal was to construct a prognostic risk assessment model based on the expression of these lncRNAs. The reliability of the prognostic risk assessment model was further validated in two independent datasets. Target genes regulated by the prognostic lncRNAs were investigated, lncRNA-mRNA networks were built, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed to decipher the possible biological function of the prognostic lncRNAs in the pathogenesis of lung adenocarcinoma.

Materials and methods

Data

TCGA database (https://gdc-portal.nci.nih.gov) was used to download mRNA-seq expression profiles of lung adenocarcinoma, which included 513 lung adenocarcinoma samples detected in the Illumina HiSeq 2000 RNA Sequencing platform. The NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo) was also used for a search based on the key word ‘lung adenocarcinoma’ to retrieve all publicly uploaded expression profiles. According to the different data screening criteria, dataset I and dataset II were obtained. The screening criteria of dataset I were gene expression data, solid tissue lung cancer specimens (non-blood, non-cell lines etc.), human expression profiles and a total number of samples ≥80. A total of three chip datasets [GSE30219 (12), GSE37745 (13) and GSE50081 (14)] were included in the analysis, which contained 84, 106 and 127 samples, respectively. The detection platform was the Affymetrix-GPL570 (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA).

The screening criteria of dataset II were gene expression data, lung adenocarcinoma tumor tissue specimens (non-blood and non-cell lines) patients, presence of control tissues, human expression profiles and a total number of samples ≥100. Finally, a total of four sets was included in the analysis of the microarray datasets [GSE32863 (15), GSE75037 (16), GSE10072 (17) and GSE43458 (18)]. They included 116 (58 tumor and 58 normal control), 166 (83 tumor and 83 normal control), 107 (58 tumor and 49 normal control) and 110 (80 tumor and 30 normal control) samples, respectively. The detection platforms for the samples were the Illumina-GPL6884 (Illumina, Inc., San Diego, CA, USA), Affymetrix-GPL96 (Affymetrix; Thermo Fisher Scientific, Inc.) and Agilent-GPL6244 (Agilent Technologies, Inc., Santa Clara, CA, USA) systems.

Data preprocessing

The original format of the downloaded expression spectrum data was divided into three different formats, according to the different detection platforms. The original expression data were preprocessed in the following three different ways. In the first approach, the expression profile data files were downloaded from the Affy platform in the original CEL format, and their formats were transformed using the oligo package, version 1.41.1, in R language 3.4.1 (19) (http://www.bioconductor.org/packages/release/bioc/html/oligo.html). The median method was supplemented using the missing values, and the data were standardized using MAS 5.0 algorithm (19) and the quantile method.

In the second approach, the expression profile data file was downloaded from the Illumina and Agilent platforms in the original TXT format. The log2 logarithm was performed using the Limma package, version 3.34.0, in R language 3.4.1 (20) (https://bioconductor.org/packages/release/bioc/html/limma.html) to transform the expressed data from a partial distribution to an approximate normal distribution. The data normalization was performed using the median method.

In the third approach, fragments per kilobase of exon per million fragments mapped expression profile data downloaded from TCGA dataset were normalized using the quantile standardization method in the preprocessCore package, version 1.40.0, of R language 3.4.1 (21) (http://bioconductor.org/packages/release/bioc/html/preprocessCore.html).

Ref_seq and Transcript_ID provided by the annotation platforms were used to annotate the lncRNAs. Comparisons between the detection sequence provided by the platform and the human genome sequence (version GRCh38) were performed using Clustal 2 (22) (http://www.clustal.org/clustal2). Finally, multiple annotation results were combined to identify each lncRNA and its corresponding expression information.

Selection of stable modules using weighted gene co-expression network analysis (WGCNA)

As a bioinformatics algorithm, WGCNA is commonly used to construct co-expression networks to identify modules associated with diseases to screen important pathogenic mechanisms or potential therapeutic targets (23). In this study, the stable modules associated with lung adenocarcinoma were selected from the co-expression networks built using the WGCNA package, version 1.61, (https://cran.r-project.org/web/packages/WGCNA/index.html) (24) in R language 3.4.1, on the basis of the GSE50081 dataset. The GSE30219 and GSE37745 datasets in this section were used as validation sets. The WGCNA algorithm was executed in three steps: i) Calculation of the correlated expression between two of the three aforementioned datasets selected randomly; ii) definition of adjacency function; and iii) module partition. A module containing ≥100 RNAs with cutHeight = 0.99 was set as the screening threshold for module partition. The functions of the significantly stable modules were annotated using the userListEnrichment function, and lncRNAs in the stable module were defined as those with an association with lung adenocarcinoma.

MetaDE analysis of expression differences in integrated multi-data

The meta-analysis of the GSE32863, GSE75037, GSE10072 and GSE43458 datasets was performed using MetaDE.ES of the MetaDE package (https://cran.r-project.org/web/packages/MetaDE) (25,26) in R language 3.4.1. The aim was a comprehensive screen of significantly differently expressed RNAs with consistent expression across these four datasets between lung adenocarcinoma and control samples. Thresholds of tau2 = 0 and Qpval >0.05 were set as consistency selection parameters. P-values and false discovery rate values <0.05 were selected as the significant difference screening parameters to identify expression differences in lncRNAs.

Construction of risk prediction model

Based on the important lncRNAs that were validated as being associated with lung adenocarcinoma through the aforementioned analyses, the optimal prognostic lncRNA combination was screened using the Cox proportional hazards (Cox-PH) (27) model, according to the L1-penalized regular regression algorithm. The screening model was derived from the penalized package (28) (http://bioconductor.org/packages/penalized/) in R language 3.4.1, in which the optimized parameter of ‘lambda’ was obtained using the 1,000 cycle calculation of cross-validation likelihood (CVL). The sample risk assessment system was established based on the lncRNA expression, which was weighted by the regression coefficient of each lncRNA in the optimal combination, and the risk-score of each sample was obtained. The risk assessment score was calculated as follows:

Risk-score = β lncRNA1 × exprlncRNA1 + ··· + β lncRNAn × exprlncRNAn.

The GSE37745 dataset was used as a validation dataset, which was performed using WGCNA analysis along with GSE50081, and TCGA dataset was used as an independent validation dataset to evaluate the risk of concentrating samples. The viability of the risk-score model was judged by the significance of the prognosis difference between the two groups, with high and low risks identified by this model.

Analysis of important lncRNA-associated pathways

KEGG pathways were enriched using the Gene Set Enrichment Analysis (GSEA) method (http://software.broadinstitute.org/gsea/index.jsp) (29). Signaling pathways that significantly correlated with lncRNAs were screened out when the gene sets had been extracted from the corresponding important lncRNA modules. The basic principle of GSEA is to use a predefined gene set that usually comes from functional annotations, such as the KEGG pathway, or previous experimental results, to rank genes according to their expression levels in samples, and to test whether the predefined gene set is enriched in the top or the bottom of the ranking table. GSEA analysis is able to evaluate gene sets rather than single gene expression changes; thus, it may contain these subtle expression changes, with better results expected.

Results

Screening of stable modules significantly correlated with lung adenocarcinoma using WGCNA algorithm

RNA of each dataset was annotated when the expression profiles of the downloaded datasets had been standardized. There were 15,988 mRNAs and 851 lncRNAs shared in the GSE30219, GSE37745 and GSE50081 datasets. Using GSE50081 expression data as the training set and the remaining data as the validation set, RNA modules significantly associated with lung adenocarcinoma were selected using the WGCNA algorithm. The steps and results are described as follows.

Correlation analysis of RNA expression among the three datasets

The consistency of expression values of the common RNAs was first checked in the three datasets to ensure the comparability of RNA expression levels in each dataset. The results indicated that the correlation distribution of RNA expression levels among the three datasets always exceeded 0.9, while the P-values were all <1×10−200 (Fig. 1). The results suggested a significantly positive correlation between each paired combination of the datasets.

Figure 1. Correlation analysis of expression levels between each pair of datasets including GSE30219, GSE37745 and GSE50081. (A) GSE50081-GSE30219. (B) GSE50081-GSE37745. (C) GSE30219-GSE37745.

Selection of β parameter value in the adjacency function

A prerequisite of WCGNA analysis is a scale-free network distribution. Thus, the appropriate adjacency matrix weight parameter β (power) should be chosen to ensure that the constructed co-expression network is as close to a scale-free distribution as possible. With a β value ranging between 1 and 20, a linear model was established using the logarithms of the adjacency degree of a node (log k) and the probability of the node appearance [log P (k)], respectively. The parameter β is the square value of coefficient R. A higher R2 value indicates a closer association between the network and the scale-free distribution (Fig. 2A). In the present study, the β (power) of 5 when the R2 value was ~0.9 for the first time was finally chosen. This selection ensured that the network connection was close to a scale-free distribution, and made the curve smooth as the minimum threshold. When the β value was equal to 5, the average connection degree of RNAs in the network was 5 (Fig. 2B), consistent with the small world network property in the scale-free network.

Figure 2. Determination of the β parameter value in the adjacency function for construction of a scale-free network. (A) The selection diagram of adjacency matrix weighting parameter β (power). The horizontal axis depicts the weight parameter power, and the longitudinal axis depicts the square value R2 of the correlation coefficient between log k and log P (k) in the corresponding network. The red line indicates the standard line with an R2 value of 0.9. (B) Schematic diagram of the RNA average connectivity under the conditions of different power parameters.

Screening of stable RNA function modules associated with lung adenocarcinoma

The RNA adjacency matrix, which is the correlation matrix among RNAs, was constructed on the basis of the GSE50081 data set as the training set, with the β value set at 5. The hierarchical clustering tree was built based on this matrix, in which each RNA module contained at least 100 RNAs, according to the standard of the hybrid dynamic shear tree, with the pruning height set at cutHeight = 0.99. A total of eight modules were obtained: M1-black, M2-blue, M3-brown, M4-green, M5-gray, M6-red, M7-turquoise, and M8-yellow (Fig. 3A). The stability of the modules in training set GSE50081 was evaluated by partitioning the corresponding module in the other validation datasets, GSE30219 and GSE37745 (Fig. 3B and C).

Figure 3. Corresponding tree diagrams of RNA module partition on the analyzed datasets. (A) GSE50081. (B) GSE30219. (C) GSE37745. Each color represents an individual module.

The results of module partition and correlation in the GSE50081 dataset are displayed in Fig. 4. RNAs in the same module (denoted by the same color dots) aggregated together (Fig. 4A), indicating that RNAs in the same module had similar expression. The brown-red module and blue-green-yellow module exhibited the characteristics of independent branches (Fig. 4B).

Figure 4. Partition and correlation of modules in the GSE50081 dataset. (A) The MDS graph for each gene in the corresponding module. The X and Y axes represent the first and second principal components, respectively. (B) The system tree diagram of the modules. MDS, multidimensional scaling; ME, module eigengene.

The eight modules were partitioned and their stabilities were analyzed. Except for the gray, turquoise and black modules, the stability scores (preservation Z-scores) of the remaining six modules exceeded 5 (Table I), indicating that they were all stable. The top five modules with high stability were the red, blue, yellow, green and brown modules, in which RNAs were likely to be relevant to the pathogenesis of lung adenocarcinoma. Table I lists the functional annotation information of these eight modules, of which the top three were red, blue and brown. A total of six lncRNAs in the red module, 37 in the blue module, and 24 in the brown module were associated with cellular immune responses, the cell cycle and focal adhesion, respectively.

Table I. Details of preservation Z-scores and annotation in the eight stable modules of the GSE50081, GSE30219 and GSE37745 datasets.

Table I.

Details of preservation Z-scores and annotation in the eight stable modules of the GSE50081, GSE30219 and GSE37745 datasets.

Module	Color	Module size	mRNAs	lncRNAs	Preservation z-score	Module GO annotation
Module 7	Turquoise	911	887	24	3.1526	Mitotic nuclear division
Module 5	Grey	792	771	21	4.8675	Cell-cell signaling
Module 1	Black	129	121	8	5.7128	Epidermis development
Module 4	Green	290	283	7	19.2659	Epithelial cell proliferation
Module 8	Yellow	306	288	18	19.5524	Oxidation-reduction process
Module 3	Brown	318	294	24	22.9983	Focal adhesion
Module 2	Blue	413	376	37	47.6821	Cell cycle
Module 6	Red	246	240	6	50.8970	Immune response

[i] A higher Z value indicates greater module stability: Stable modules, 5<Z<10; and highly stable modules, Z>10. GO, Gene Ontology; lncRNA, long non-coding RNA.

The correlation between each module and the various clinical factors was calculated by integrating the clinical information of samples in the GSE50081 dataset. The results presented in Fig. 5 indicated that the top three modules in terms of stability (red, blue, and brown) were significantly correlated with the T and N staging (T, the range of the primary tumors; N, the metastasis of peripheral lymph nodes) and recurrence of lung adenocarcinoma. The remaining two modules with high stability (green and yellow) were significantly associated with smoking. Therefore, lncRNAs in these five modules, which were considered to be important factors associated with lung adenocarcinoma, were focused upon.

Figure 5. Heat map of the correlation between the RNA functional modules and their clinical attributes, based on the GSE50081 dataset. Red asterisks represent a significant positive correlation. ME, module eigengene; T, tumor; N, node.

Screening of important lncRNAs associated with lung adenocarcinoma

Based on the integration of the GSE32863, GSE75037, GSE10072 and GSE43458 datasets, a total of 1,691 RNAs with consistently significant expression were screened through the comprehensive analysis of the MetaDE package, including 58 differentially expressed lncRNAs and 1,633 differentially expressed mRNAs. Considering the previous results of module stability, annotation and correlation between modules and clinical factors, lncRNAs in the brown, red, blue, green and yellow modules were regarded as the factors significantly associated with the incidence and prognosis of lung adenocarcinoma. In addition, their differential expression was considered to be associated with cancer pathogenesis. Therefore, combining the two aforementioned aspects, 32 differentially expressed lncRNAs were selected from the brown, red, blue, green and yellow modules for subsequent survival analysis in relation to lung adenocarcinoma. The distributions in the eight modules are presented in Fig. 6.

Figure 6. Distribution number and proportion of the long non-coding RNAs with significantly different expression in each color module.

Establishment and evaluation of risk assessment model

Selection of optimal lncRNAs

According to the results of the phylogenetic tree (Fig. 4B), 17 and 15 differentially expressed lncRNAs in the blue-green-yellow and brown-red modules, respectively, were divided into two groups. The third group was defined as the combined group, which included these 32 lncRNAs. The lncRNAs in each group were further optimized and screened using the Cox-PH model based on the L1-penalized regularization regression algorithm in the penalized package to obtain the optimal prognosis-associated lncRNAs. As presented in Fig. 7A, the 1,000-cycle calculation of the CVL algorithm was performed to obtain the optimal lambda values in the Cox-PH model. The lambda value corresponding to the maximum CVL value was selected as the input value in the Cox-PH model. Using the Cox-PH model, the best combinations composed of 5, 5 and 6 lncRNAs emerged in the blue-green-yellow, brown-red and combined modules, respectively. The corresponding coefficients of the three selected sets of the optimal combination of lncRNAs are presented in Fig. 7B (the algorithm yielded different coefficients to the optimal lncRNAs, and the coefficients of the excluded lncRNAs were all reduced to zero).

Figure 7. Lambda parameter filtration and coefficient distribution of the best prognosis-associated lncRNAs. (A) Curve graphs of lambda parameters screened by CVL in the blue-green-yellow, brown-red and combined groups, respectively. The horizontal and vertical axes represent the different values of lambda and CVL, respectively. The intersection of the dotted lines indicates the value of the lambda parameter when the CVL is maximal. (B) Distribution charts of the optimal prognosis-associated lncRNAs selected by the Cox proportional hazards model when the value of lambda parameter corresponding to the maximum CVL value was taken. Each point corresponds to each lncRNA. lncRNA, long non-coding RNA; CVL, cross-validation likelihood; CPS1-IT1, CPS1 intronic transcript 1; LINC01140, long intergenic non-protein coding RNA 1140; LINC01711, long intergenic non-protein coding RNA 1711; SNHG14, small nucleolar RNA host gene 14; TP73-AS1, TP73 antisense RNA 1; DIAPH2-AS1, DIAPH2 antisense RNA 1; LINC00652, long intergenic non-protein coding RNA 652; FOXN3-AS2, FOXN3 antisense RNA 2; MEG3, maternally expressed 3; RHPN1-AS1, RHPN1 antisense RNA 1 (head to head).

Establishment, evaluation, and selection of risk prediction model

Based on the corresponding Cox-PH regression coefficients of the three optimal lncRNAs obtained in the previous step, three sets of prediction models of sample risk scoring were constructed as follows:

Risk-score (blue-green-yellow) = (0.182) × ExpCPS1-IT1 + (−0.060) × ExpLINC01140 + (0.244) × ExpLINC01711 + (−0.259) × ExpSNHG14 + (−0.480) × ExpTP73-AS1;

Risk-score (brown-red) = (0.381) × ExpDIAPH2-AS1 + (−0.097) × ExpFOXN3-AS2 + (0.738) × ExpLINC00652 + (0.041) × ExpMEG3 + (−0.0003) × ExpRHPN1-AS1;

Risk-score (combined) = (0.122) × ExpCPS1-IT1 + (−0.065) × ExpLINC01140 + (0.134) × ExpLINC01711 + (−0.262) xExpSNHG14 + (−0.445) × ExpTP73-AS1 + (0.017) × ExpMEG3.

Based on these three models, the risk-score of each sample in training dataset GSE50081 was calculated, and these samples were divided into high and low-risk groups according to the median score. Kaplan-Meier survival curve analysis was used to assess the significant differences in the survival time between the two groups. All three constructed models, including blue-green-yellow, brown-red and combined, had a significantly different effect on the high and low-risk groups (P<0.05 for all models) of the GSE50081 training set (Fig. 8A).

Figure 8. Kaplan-Meier survival curve of samples in the analyzed datasets. (A) Training set GSE50081. (B) Verification set GSE37745. (C) Verification set TCGA. The black and red curves represent the low-risk and high-risk groups, respectively. TCGA, The Cancer Genome Atlas.

GSE50081 together with the GSE37745 dataset in the WGCNA analysis were used as the validation datasets to evaluate the efficiency of the prediction models. Only the lncRNAs in the brown-red model were able to better separate the low-risk and high-risk prognostic samples in the GSE37745 dataset (Fig. 8B; P=0.015). Furthermore, the effectiveness of the three models was further verified by the TCGA independent validation dataset. Only the brown-red model produced a clear discrimination result (Fig. 8C; P=0.0037), suggesting that the brown-red model was very robust, in which five lncRNAs [DIAPH2 antisense RNA 1 (DIAPH-AS1), FOXN3 antisense RNA 2 (FOXN3-AS2), long intergenic non-protein coding RNA 652 (LINC00652), maternally expressed 3 (MEG3) and RHPN1 antisense RNA 1 (head to head) (RHPN1-AS1)] were significantly correlated with prognosis, were stable in distinguishing between the prognosis samples, and could be universally applied to other new samples.

Network construction and analysis of functional pathways

A total of five important lung adenocarcinoma-associated lncRNAs, including DIAPH2-AS1, FOXN3-AS2, LINC00652, MEG3 and RHPN1-AS1, were obtained according to the previous analysis. Combined with their associated mRNAs selected in the brown and red modules using WGCNA, the important lncRNA-mRNA networks were built, as presented in Fig. 9.

Figure 9. lncRNA-mRNA networks corresponding to the red and brown modules. (A) Red module. (B) Brown module. The circular nodes represent the significantly differentially expressed genes (up- and downregulated) screened by MetaDE. The square nodes represent lncRNAs; the green and red edges represent positive and negative expression, respectively. lncRNA, long non-coding RNA.

GSEA-based pathway enrichment analysis was performed on the genes in two lncRNA-mRNA networks, in which three key statistical values were used. The first value was the enrichment score (ES) that was the original result of the GSEA analysis, reflecting the enrichment degree value of a functional gene set at the anterior or posterior of an array when all the hybridization data have been sorted. The basic principle of the calculation is to scan the sequence of sequences. When a gene of this set is prepared, the ES value increases; otherwise, it decreases. The second value is the normalized enrichment score (NES, which is the standardized processing of ES values. The final value is the nominal P-value, which describes the statistical significance of the ES of a functional gene subset. A smaller P-value indicates better gene enrichment. When the absolute value of NES is greater, the P-value will be spontaneously smaller, suggesting a higher degree of enrichment and greater confidence in the analysis result. In the present study, P<0.05 was chosen as the threshold for screening KEGG pathways in which the relevant module genes were significantly enriched.

As illustrated in Fig. 10A and Table II, six pathways were significantly associated with the RHPN1-AS1 network in the red module. In the brown module, three pathways were simultaneously associated with the DIAPH2-AS1, FOXN3-AS2, LINC00652, and MEG3 networks (Fig. 10B and Table III).

Figure 10. Pathways enriched in the red and brown network modules, and heat maps of the associated genes. (A) Red module. (B) Brown module. The x-axes indicate the genes involved in the pathways; the y-axes indicate the enriched pathways. A deeper red color indicates a stronger positive correlation. MEG3, maternally expressed 3; FOXN3-AS2, FOXN3 antisense RNA 2; LINC00652, long intergenic non-protein coding RNA 652; DIAPH2-AS1, DIAPH2 antisense RNA 1; PRKG1, protein kinase cGMP-dependent 1; MYL9, myosin light chain 9; KCNMB1, potassium calcium-activated channel subfamily M regulatory β-subunit 1; ACTG2, actin, γ2, smooth muscle, enteric; ADCY2, adenylate cyclase 2; CCL11, C-C motif chemokine ligand 11; TNFSF4, TNF superfamily member 4.

Table II. KEGG pathways significantly associated with the red module.

Table II.

KEGG pathways significantly associated with the red module.

KEGG pathways	ES	NES	NOM p-val
T_CELL_RECEPTOR_SIGNALING_PATHWAY	0.3936	1.3243	0.0001
JAK_STAT_SIGNALING_PATHWAY	0.2300	1.1451	0.0158
NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY	0.4680	1.0535	0.0285
TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY	0.6984	1.0162	0.0300
B_CELL_RECEPTOR_SIGNALING_PATHWAY	−0.2269	−1.0968	0.0333
CELL_ADHESION_MOLECULES_CAMS	0.3157	1.0298	0.0394

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; ES, enrichment score; NES, normalized enrichment score; NOM p-val, nominal P-value.

Table III. KEGG pathways significantly associated with the brown module.

Table III.

KEGG pathways significantly associated with the brown module.

A, MEG3
KEGG pathways	ES	NES	NOM p-val
VASCULAR_SMOOTH_MUSCLE_CONTRACTION	−0.4164	−1.0637	0.0344
CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	0.3187	1.3840	0.0468
CHEMOKINE_SIGNALING_PATHWAY	0.2642	1.0044	0.0448
B, LINC00652
KEGG pathways	ES	NES	NOM p-val
VASCULAR_SMOOTH_MUSCLE_CONTRACTION	−0.3550	−0.9555	0.0066
CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	0.2316	1.0048	0.0476
CHEMOKINE_SIGNALING_PATHWAY	−0.2200	−0.8889	0.0482
C, FOXN3-AS2
KEGG pathways	ES	NES	NOM p-val
VASCULAR_SMOOTH_MUSCLE_CONTRACTION	−0.4266	−1.1554	0.0180
CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	0.2684	1.1074	0.0800
CHEMOKINE_SIGNALING_PATHWAY	0.2802	1.0802	0.0860
D, DIAPH2-AS1
KEGG pathways	ES	NES	NOM p-val
VASCULAR_SMOOTH_MUSCLE_CONTRACTION	−0.3925	−1.0452	0.3050
CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	0.2130	0.9520	0.6348
CHEMOKINE_SIGNALING_PATHWAY	0.2426	1.0289	0.3527

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; ES, enrichment score; NES, normalized enrichment score; NOM p-val, nominal P-value; MEG3, maternally expressed 3; FOXN3-AS2, FOXN3 antisense RNA 2; LINC00652, long intergenic non-protein coding RNA 652; DIAPH2-AS1, DIAPH2 antisense RNA 1.

The pathway enrichment analysis demonstrated that the RHPN1-AS1 network was significantly associated with cell receptor-mediated signaling and cell adhesion pathways (Table II). A total of 26 genes, which included C-X-C motif chemokine ligand (CXCL)1, interleukin (IL)6, CXCL5, IL15, IL2 inducible T cell kinase (ITK), RAS guanyl releasing protein 1 (RASGRP1), protein tyrosine phosphatase, receptor type C (PTPRC), cluster of differentiation (CD)274, CD19, CD72 and selectin E (SELE), were jointly involved in the six pathways and positively associated with them. Furthermore, the DIAPH2-AS1, FOXN3-AS2, LINC00652 and MEG3 networks were significantly associated with ‘cytokine cytokine receptor signaling pathway’, ‘chemokine signaling pathway’ and ‘vascular smooth muscle contraction pathway’ (Table III). A total of seven genes, which included C-C motif chemokine ligand 11 (CCL11), TNF superfamily member 4 (TNFSF4), adenylate cyclase 2 (ADCY2), actin, γ2, smooth muscle, enteric (ACTG2), myosin light chain 9 (MYL9) and protein kinase cGMP-dependent 1 (PRKG1), were commonly involved in these three functional enrichment pathways, with significant positive correlations between genes and pathways.

Discussion

It has been considered that lncRNAs serve important roles in the development of lung adenocarcinoma, and may be used as diagnostic markers for tumors (3,9–11). To elucidate an lncRNA-based signature that predicts the prognosis of patients with lung adenocarcinoma, a large amount of of RNA-seq data and clinical survival prognosis information from patients with adenocarcinoma was initially downloaded from TCGA and the GEO database, and a co-expression network was constructed to mine network modules with particular biological functions. A total of eight modules were obtained in the GSE50081 training set. Of these, five (red, blue, yellow, green and brown modules) featured high stability and appeared likely to be associated with the pathogenesis of lung adenocarcinoma. Integrated module analyses that included stability, functional annotations, correlation analysis with clinical factors and the differential expression of lncRNAs, identified 32 lncRNAs in the brown, red, blue, green and yellow modules. These were used for follow-up survival analysis. Based on the results of the system clustering tree, Cox-PH model and Kaplan-Meier survival curve, three risk scoring systems were constructed, among which five lncRNAs (DIAPH2-AS1, FOXN3-AS2, LINC00652, MEG3 and RHPN1-AS1) comprised a signature-based risk scoring system that successfully classified the low-risk and the high-risk samples in the training set.

Considering that MEG3 has been related to the pathogenesis of lung adenocarcinoma (11,30,31), it is possible that the other four lncRNAs may be potential novel prognostic factors for lung adenocarcinoma. DIAPH2 is a mammalian homolog of Drosophila diaphanous. The latter is a protein required for cytokinesis, and belongs to a family of formin-associated proteins containing repetitive polyproline stretches (32). DIAPH2 interacts with RhoD (33). However, there is no information on the association between DIAPH2-AS1 and human cancer. FOXN3 is an important member of the FOX family of transcription factors that is essential in organ differentiation, development, cell growth and cancer (34). FOXN3-AS2 may serve an important role in esophageal cancer (35). The function of FOXN3-AS2 in lung adenocarcinoma is not well understood. The present study provides the first clue that FOXN3-AS2 may be a potential diagnostic marker for lung adenocarcinoma, to the best of our knowledge. RHPN1 was originally identified as a RhoA GTPase-interacting partner (36). It was recently implicated as being essential for the integrity of the glomerular filtration barrier and is a key determinant of podocyte cytoskeleton architecture (37). RHPN1-AS1 has a potential role in the progression of uveal melanoma and might be an attractive biomarker and therapeutic target in uveal melanoma (38). LINC00652 is located in human chromosome 20; no physiological function has been previously reported, to the best of our knowledge. The risk stratification capability of the five-lncRNAs signature was confirmed in the validation set and independent validation set. The present data reveal a potential biomarker that may predict the prognosis of lung adenocarcinoma patients. The data may be helpful in future explorations of the pathogenesis of lung adenocarcinoma.

It has been demonstrated that lncRNAs serve important roles in a variety of biological processes by regulating target genes at the transcriptional, posttranscriptional, and epigenetic levels (39,40). Therefore, the present study attempted to investigate the target genes regulated by the five prognostic lncRNAs to decipher their possible biological functions in the pathogenesis of lung adenocarcinoma. Since the results of the module functional annotation demonstrated that lncRNAs in the red and brown modules were primarily associated with the cellular immune response and focal adhesion, respectively, these two modules were used to select the associated mRNAs and to build the important lncRNA-mRNA networks. Further KEGG pathway enrichment analysis revealed that six pathways were significantly associated with the RHPN1-AS1 network in the red module. These included cell receptor-mediated signaling and cell adhesion pathways, which have notable roles during the multistage progression of human carcinogenesis (41). Therefore, it was hypothesized that RHPN1-AS1 has the same directional correlation with the six important pathways by positively regulating the expression level of CXCL1, IL6, CXCL5, IL15, ITK, PRSGRP1, PTPRC, CD274, CD19, CD72, SELE and other genes. CXCL1 has been associated with lung cancer (42), upregulation of CD274 is associated with the poor prognosis of lung adenocarcinoma (43,44), and IL6 (45) and IL15 (46) are relevant to the pathogenesis of lung adenocarcinoma. Another three pathways, including ‘cytokine cytokine receptor signaling pathway’, ‘chemokine signaling pathway’ and ‘vascular smooth muscle contraction pathway’, were significantly associated with the DIAPH2-AS1, FOXN3-AS2, LINC00652 and MEG3 networks in the brown module. A possible explanation is that DIAPH2-AS1, FOXN3-AS2, LINC00652 and MEG3 have the same directional correlation with the three important pathways by positively regulating the expression levels of CCL11, TNFSF4, ADCY2, ACTG2, MYL9 and RPKG1. It was reported that the differential expression of ACTG2 is useful in distinguishing lung adenoma from normal lung samples (47), and low MYL9 expression might be associated with the development and metastasis of non-small cell lung cancer (48). Moreover, CCL11 is reported to have a direct and selective profibrogenic effect on lung and bronchial fibroblasts (49). The present data revealed the possible pathogenesis of lung adenocarcinoma, in which five lncRNAs (DIAPH2-AS1, FOXN3-AS2, LINC00652, MEG3 and RHPN1-AS1) may affect the occurrence and development of lung adenocarcinoma by regulating the expression levels of target genes.

It is noteworthy that this study was an extensive bioinformatics study based on published data. The results require further validation by in vitro or in vivo models. These useful clues may help other researchers to perform relevant research. In a further study, the expression levels of the five lncRNAs in clinical samples with prognostic information may be measured using experimental methods to validate the predictive value of this five-lncRNA signature.

In conclusion, the present study identified and validated a five-lncRNA signature predicting the prognosis of patients with lung adenocarcinoma. The predictive ability of this signature may be exploited as a promising prognostic biomarker for lung adenocarcinoma. The target genes of these five prognostic lncRNAs are associated with a number of cellular processes and signaling pathways, including the cell receptor-mediated signaling pathway and cell adhesion pathway. Therefore, these five prognostic lncRNAs may be potential diagnostic markers. The present results may be helpful to elucidate the possible pathogenesis of lung adenocarcinoma.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

JL performed data analyses and wrote the manuscript. XY, SO, KL, YK, HL, YO and QL contributed significantly to the data analyses and manuscript revision. RX conceived and designed the study. All authors read and approved the final manuscript.

Ethics approval and consent to participate

In the original article of the datasets, the trials were approved by the local institutional review boards of all participating centers, and informed consent was obtained from all patients.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References