Two predicted models based on ceRNAs and immune cells in lung adenocarcinoma

Identification of differentially expressed ceRNAs

In the present study, expression profiles of LUAD and the corresponding clinical information were derived from The Cancer Genome Atlas (TCGA) database, including lncRNA,miRNA, and mRNA. The R package “DEseq2” was used to screen the differentially expressed mRNAs,miRNAs and lncRNAs (DEmRNA, DEmiRNA and DElncRNA) between normal and tumor samples (Love, Huber & Anders, 2014). Specifically, a false discovery rate (FDR) <0.05 and log2 —fold change —> 1 were set as cutoff criteria.

ceRNA network construction

Both lncRNA-miRNA and miRNA-mRNA interactions were predicted using, miRcode and starBase (Chou et al., 2018; Jeggari, Marks & Larsson, 2012; Li et al., 2014). DEmiRNAs, DEmRNAs, and DE lncRNAs, which presented a significant difference in hypergeometric distribution detection and correlation analysis (P < 0.05), were selected to construct a lncRNA-miRNA-mRNA regulatory network and then was visualized by Cytoscape 3.80 (Shannon et al., 2003).

Construction of nomogram based on key ceRNAs

Univariate Cox analysis was employed to screen the survival-associated ceRNAs. Lasso regression analysis was applied to ensure that the Cox model was not overfitted. Eventually, the candidate ceRNAs were integrated to establish the multivariate cox regression models. Risk scores were described as following: risk score=(Exp RNA1* β1) + (Exp RNA2* β2 ) + ... + (Exp RNAn* βn). It is worth noting that Exp stands for the expression level, and β denotes the regression coefficient. Based on the risk scores, patients with LUAD were divided into high-risk and low-risk groups. The corresponding Kaplan–Meier survival curves were carried out to reveal the overall survival (OS) within different groups. In addition, we also performed the ROC curve to evaluate the specificity and sensitivity of the model. Finally, according to the multivariate model result, we constructed a nomogram to predict patients’ prognostic values. The calibration curve was applied to evaluate the accuracy of the nomogram.

Construction of nomogram based on immune cells

With a good deconvolution performance in gene expression profiles, CIBERSORT method could estimate the proportion of particular cell types (Newman et al., 2015). In this study, CIBERSORT algorithm was utilized to estimate 22 tumour-infiltrating immune cells in LUAD by collating and calculating the genes expression in the tumor samples in the TCGA database. Samples with a CIBERSORT output of p < 0.05 presented that the predicted proportion of immune cells were correct.

The Wilcoxon test was conducted to examine the significant difference of immune cells between the tumor samples and healthy samples. We also performed multivariate Cox regression and produced the corresponding risk scores. Based on the multivariate Cox model results, we finally built a predicted nomogram. Similar to ceRNAs, patients were also stratified into high-risk and low-risk groups on the basis of the mean risk score. Subsequently, the prognostic model’s sensitivity and specificity were examined by the ROC curve, and the calibration curve investigated the accuracy of the nomogram.

Validation of prognostic equations with GEO

To better elucidate the accuracy of both multivariate Cox models, an external validation dataset GSE72094 (n = 442) was applied in the present study. We selected this dataset because it contained most LUAD patients with clinical data in the Gene Expression Omnibus (GEO) database. Based on the median risk score generated from TCGA, we also classified patients with high- and low-risk patients. K-M and ROC curves were applied to evaluate the prognostic efficacy of multivariate models in TCGA.

Estimation of tumor stemness and microenvironment

Stemness scores of “RNAss” and “DNAss” referred to the result based on mRNA and DNA-methylation respectively. The larger stemness scores were, the higher likelihood more stem cells were infiltrating. ESTIMATE algorithm could speculate the level of stromal cell, immune cells and tumor purity, respectively (Yoshihara et al., 2013). Correspondingly, stromal/immune/Estimate scores were significantly associated with stromal cells, immune cells and tumor purity, respectively. We utilized the “ESTIMATE” package to assess immune, stromal scores, and the sum of both in individual patients. The higher scores, the higher proportion of corresponding cells were. Additionally, the correlation of hub ceRNAs with above scores was examined by Spearman analysis. The ssGSEA method could apply expressed traits of immune cell population to individual tumor samples (Barbie et al., 2009). Based upon the results of ssGSEA, we calculated infiltrating level of immune cells and immune-associated functions in LUAD samples by using the “gsva” package (Hänzelmann, Castelo & Guinney, 2013). Furthermore, we compared the differences in these immune data sets between the high-risk group and the low-risk group.

Statistical analysis

All statistical analyses were accomplished with R software 4.0.2. The Wilcoxon rank sum tests were conducted for comparisons between two groups. Only a two-sided P-value <0.05 was recognized as statistically significant.

Identification of DEmRNAs, DEmiRNAs and DElncRNAs

Figure 1 presents our workflow for the bioinformatic analysis. Among the whole expression profiles, totally 1645 upregulated DEmRNAs, and 1,344 downregulated DEmRNA were screened out (Figs. 2A and 2B). Furthermore, we also identified 163 upregulated DElncRNA, 52 downregulated DElncRNA (Figs. 2C and 2D), 111 upregulated DEmiRNA, and 87 downregulated DEmiRNA (Fig. 3B). Specifically, the top 20 significantly upregulated- and downregulated- miRNAs were shown in Fig. 3A.

Construct Nomogram of ceRNAs

To investigate the endogenously competitive relationships between the lncRNA and mRNA, we constructed an intersected ceRNA network. Seven lncRNAs (H19, AC074117.1, FBXL19-AS1, MAGI2-AS3, SNHG3, PVT1 and SNHG1), 94mRNA, and 15miRNA were involved in the network (Table S2). Hsa-miR-130b-3p and hsa-miR-29b-3p were the top 2 regulated miRNA in the ceRNA network (Fig. 4). In addition, the initial univariate Cox regression (Fig. S1) and Lasso regression analysis were further conducted to investigate the key biomarkers (Figs. 5B and 5C). Seven survival-associated biomarkers were finally involved in a new multiple Cox regression analysis (Fig. 5A). According to the result, all biomarkers belonged to protein-coding RNAs. Consequently, the corresponding risk score was generated (Table 3): risk score = [Exp CDC14A * (-0.281888)]+(Exp LOXL2*0.114969) + (Exp CCT6A * 0.174244)+(Exp E2F7 *0.187164)+(Exp GPR37*0.099737)+(Exp H1F0 *0.197590)+(Exp SMOC1 * 0.068385).

As shown in Fig. 5D, the prognostic score could be a promising indicator to distinguish the LUAD patients according to the individual prognosis (P = 5.701e−8). The area under ROC reflected the predictive efficacy of the prognostic signature (AUC of 3-year survival:0.772) (Fig. 5E). By the GEO validation set, the AUC of gene signature at 3-year was 0.711 (Fig. S2A). Additionally, K-M analysis further confirmed the difference between high- and low-risk groups (P = 9.963e−05) (Fig. S2B). Eventually, we assessed the prognosis of patients with LUAD by a nomogram (Fig. 6A), and the accuracy of the nomogram was confirmed by the calibration curve (Fig. 6C). As shown in Fig. 7A, based on this prognostic model, the risk score in the group with fewer smoking years ( ≤20) was lower than the group with more smoking years (>20) (P = 0.025).

Construct nomogram of immune cells

By applying the CIBERSORT, we estimated the percentages of 22 immune cell types in LUAD (Figs. S3B and S3C) (Chen et al., 2018). The violin plot indicated that 13 immune cell types showed a significant difference in percentages of immune cells between tumor samples and healthy samples by the Wilcoxon test (Fig. S3D). The correlation analysis of immune cells could be found in Fig. S3A. Subsequently, after the analysis of univariate cox regression and lasso regression, seven immune cells were taken as prognostic biomarkers and eventually integrated into a new multivariate model (Figs. 8A, 8B and 8C). The formula was presented as following (Table 4): risk score = [Exp CD8 *(-3.425649 )]+(Exp Tregs*9.175062 )+[Exp Monocytes*(-5.318492)]+(Exp M1 *6.134551)+(Exp Dendritic cells activated*5.539130 )+(Exp Mast cells activated *24.611110 )+(Exp Eosinophils*37.351722). The K-M analysis results showed a significant difference in overall survival between high- and low- risk groups (P = 3.342e−05) (Fig. 8D). The ROC curve (AUC of 3-years: 0.766) demonstrated that the multivariate model was recognized as a fair potential to monitor the prognostic efficacy (Fig. 8E). In the external validation set, the model also carried out adequate predicted capacity. The AUC of immune cell signature was 0.738 at 3-years (Fig. S2C). K-M curves reflected that the low-risk group had longer overall survival than the high-risk group (P = 2.922e−04) (Fig. S2D). A predicted nomogram was also performed, and the discrimination was conducted to test the calibration quality (Figs. 6B and 6D). Similarly, we could also found that the longer smoking history was related to the higher the risk scores based on the predicted model of immune cells in smoking groups (P = 0.0046) (Fig. 7B). The correlations of key prognostic biomarkers in the two models illustrate that the expression of E2F7 is positively related to Macrohages M1 (R = 0.4) (Fig. S4).

Besides, the Univariate and multivariate survival analyses revealed that risk scores in the two models were both independent prognostic factors (Fig. S5).

Correlation of hub ceRNAs expression with tumor stemness and microenvironment

As shown in Fig. 9, CCT6A, H1F0 and E2F7 were positively related to RNAss and DNAss (R = 0.22 to 0.53). A significant negative association was detected between CDC14A and RNAss (R =  − 0.25, P = 7.5e−8). Nevertheless, we did not find a significant correlation between LOXL2, GPR37 and SMOC1 with tumor stemness. In addition, H1F0 displayed a negative association with stromal scores, immune scores and estimate scores (R =  − 0.34 to −0.27). However, CDC14A showed a positive association with the microenvironment-related scores (R = 0.31 to 0.35).

Multiple databases validation

To validate CDC14A and H1F0 expression in LUAD, we utilized multiple databases. Based on the result of UCLAN, GEPIA and TIMER, it could be concluded that H1F0 was highly expressed while CDC14A was weakly expressed in LUAD (Fig. S6). Furthermore, the CDC14A and H1F0 protein expression were obtained from The Human Protein Atlas data (HPA), indicating that H1F0 was high staining in tumor cells, while it was medium staining in normal tissues (Figs. S7A and S7B). The opposite trend was observed in CDC14A. Besides, the association of H1F0 and CDC14A with immune cells were validated in TISIDB. CDC14A expression was positively associated with Immature B cell (R = 0.345), NK cell (R = 0.308), Eosinophils (R = 0.357), and Mast cell (R = 0.332) (Fig. S8A–S8E). Controversially, H1F0 expression was negatively related to Eosinophil (R =  − 0.34), Macrophage (R =  − 0.359), NKT (R =  − 0.323), Tem_CD8 (R =  − 0.323) and Th1 (R =  − 0.407) (Figs. S8F–S8J).

Differences in ssGSEA scores in TCGA/GEO cohort

To further elucidate the immune status in risk- and low groups of ceRNA signatures, we performed the ssGSEA algorithm to estimate the immune infiltrating between two groups (Fig. S9). Of note, low-risk group were associated with higher scores of aDCs, B_cells, CD8+_T_cells, iDCs, Mast_cells, Neutrophils, pDCs, T_helper_cells and TIL (Figs. 10A and 10B). Similarly, in the lower-risk group, the scores of Check-point, HLA, T_cell_co-stimulation and Type_II_IFN_Reponse in still higher in the low-risk group.Additionally, the differences above were confirmed by the GEO cohort (Figs 10C and 10D).