Integrated genomic analysis of the stemness index signature of mRNA expression predicts lung adenocarcinoma prognosis and immune landscape

Data collection and preprocessing

(1) Data collection and preprocessing

The latest clinical follow-up information or expression data from 500 LUAD patient tissues with expression profiles of RNA-seq were sourced from TCGA database (Liu et al., 2020a). GSE30219 and GSE19188 data sets containing patients’ clinical survival data were obtained from GEO database (Barrett et al., 2013; Song et al., 2023). The mRNAsi of LUAD was acquired from a previous article (Malta et al., 2018).

(2) Data preprocessing

TCGA-LUAD samples without survival time and status and clinical follow-up information were deleted. Ensembl was transformed to the gene symbol, and the mean expression value of multiple Gene symbols corresponding to one gene was taken.

Samples from GSE30219 and GSE19188 were primarily processed as those in TCGA-LUAD. The probe was transformed into a gene symbol. Only the expression of multiple probes corresponding to a gene was averaged.

MRNAsi

Firstly, mRNAsi was acquired by using OCLR algorithm to analyze stemness characteristics of stem cells and their differentiated progenitors. MRNAsi for LUAD was acquired from a previous study to explore the association between clinical features and mRNAsi (Malta et al., 2018). High and low mRNAsi groups of LUAD patients were divided by the optimal mRNAsi value.

WGCNA

Limma package was used to draw volcano plot to analyze the differentially expressed genes (DEGs) between high and low mRNAsi under the threshold of |log 2 FC| > log 2 (1.5) and FDR <  0.05 (Sui et al., 2023). According to a total of 494 expression profiles of LUAD samples in TCGA database, the expression profiles of 1,929 DEGs were extracted, and the distance between each gene was computed using Pearson correlation coefficient. The Weighted gene co-expression network analysis (WGCNA) analysis was performed using R software to generate scale-free networks (Nangraj et al., 2020). Next, topological matrix was converted from expression-based adjacency matrix. Average-linkage hierarchical clustering method was performed with 30 gene per network as the minimum gene number when using mixed dynamic shear tree. Eigengenes of each module were calculated after determining gene modules by dynamic shearing method, followed by performing cluster analysis with deepSplit = 2, mergeCutHeight = 0.25, minModuleSize = 60. Close modules were combined into a new one. The correlation between mRNAsi and each module was further analyzed.

Construction of an mRNAsi-related risk model

The glmnet in R package was employed in univariate cox regression analysis and lasso regression analysis to reduce gene number in the model and extract key genes (Guo et al., 2021).

The following formula was utilized to compute the RiskScore for each patient: ${\displaystyle \text{RiskScore}=\sum \beta i\times \text{Exp}i}$ where i is the expression of a key gene and β is the Cox regression coefficient of the corresponding gene.

Low and high risk groups were divided by the median risk score. Kaplan–Meier (KM) method was used to draw survival curves. The timeROC in R package was employed to perform receiver operating characteristic (ROC) analysis on the risk score for prognostic prediction (Zhang et al., 2022b). Prognostic prediction for 1, 3 and 5 year(s) was assessed using the training dataset.

Screening and enrichment of DEGs between the two risk groups

The DEGs between the two groups were filtered using limma package under |log 2 FC| > 1 and FDR < 0.05. The R software WebGestaltR (V0.4.4) package was used for KEGG analysis and GO enrichment analysis between the two risk groups (Liao et al., 2019).

Immune microenvironment differences

According to immune cell genes in a previous study (Charoentong et al., 2017), ssGSEA algorithm was used to score the relative enrichment of 28 immune cell subtypes between the two risk groups. The overall immunity scores were calculated by ESTIMATE algorithm for of the TCGA data.

Immunotherapy and drug sensitivity analysis

TIDE software was used to access the potential clinical effect of immunotherapy on the two risk groups (Jiang et al., 2018), with a higher TIDE score suggesting greater immune escape and less immunotherapy benefit. The correlation between risk score and TIDE, Dysfunction and Exclusion was also compared.

The genomics of drug sensitivity in cancer (GDSC) stores the data on drug sensitivity (http://www.cancerrxgene.org). Spearman correlation analysis was employed to compute the correlation between drug sensitivity and risk sore with AUC in LUAD cell lines.

Nomogram

Risk subgroups were classified using a decision tree in terms of T. stage, N. stage, M. stage, Stage, risk score, age, and gender in the TCGA-LUAD cohort. Risk score and clinicopathological features were subjected to univariate and multifactorial cox regression analysis to determine significant prognostic factors. A nomogram was designed combining the risk score and other clinicopathological features for quantifying risk assessment and survival of LUAD patients. Calibration curve and Decision curve (DCA) were applied to assess the prediction performance and stability of the model.

Cell culture and drug treatment

In this study, human bronchial epithelial cells BEAS-2B (BNCC359274) as well as human LUAD cells NCI-H2009 (BNCC358037) and NCI-H1975 (BNCC340345) were purchased from BeiNa Culture Bio (Xinyang, China) and were used for the experiments. High-glucose DMEM (Cat. #31600, Solarbio, Beijing, China), Eagle’s Minimal Essential Medium (Cat. #LA10099, BALAB, Beijing, China), and RPMI-1640 media (Cat. #90023, Solarbio, China) were respectively used for culturing these cells. The culture media were additionally supplemented with 10% FBS (S9020, Solarbio, China). All cells were cultured in an environment with 5% CO₂ at 37 °C.

Fresh medium containing 10 µM GSK2256098 (TargetMol, Cat. #T2281) was used to replace the existing medium for GSK2256098 treatment (Zhang et al., 2022a), during which the cells were incubated for 3 h.

Transient transfection

The EIF5A-specific small interfering RNA (si-EIF5A) and the corresponding negative control (si-NC) were obtained from GenePharma Co. located in Shanghai, China. These plasmids were introduced into H1975 and H2009 cells using Lipofectamine 3,000 reagents (Cat. #L3000-001) from Invitrogen (Carlsbad, CA, USA) in the USA. AGCAAGTTTTCTGAAACGTTTGG (si-EIF5A) was the target sequence for EIF5A siRNA. After a period of 48 h, the efficiency of the transfection was evaluated via qRT-PCR.

QRT-PCR

Total RNA was extracted from the cells (including LUAD cells and BEAS-2B cells) using TRIzol reagent (Cat. #15596-026, Thermo Fisher, Waltham, MA, USA). The quality of extracted RNA was determined by NanoDrop (Thermo Fisher, USA) to ensure the purity of RNA. To synthesize cDNA, we used HiScript II SuperMix (Cat. #R223-01, Vazyme, Nanjing, China) for the reverse transcription reaction. Each reaction system contained 500 ng of RNA, 2 µL of 5 × HiScript II SuperMix, RNAase-free water, and was adjusted to a final volume of 10 µL. The conditions of the reverse transcription reaction were as follows: pre-incubation at 25 °C for 10 min, reaction at 42 °C for 1 h, and termination of the reaction at 85 °C. The cDNA after reverse transcription was stored at −20 °C for backup. The qRT-PCR experiments were performed using SYBR Green fluorescent dye (Roche, Basel, Switzerland) for real-time quantitative PCR, and the experiments were performed in an ABI 7500 real-time quantitative PCR system (Thermo Fisher, Waltham, MA, USA). PCR amplification was performed under the following conditions: initial denaturation at 94 °C for 10 min; followed by 45 cycles: 94 °C for 10 s and 60 °C for 45 s. ACTB was acted as the internal reference. Experiments were performed in triplicate. Data results were normalized by the 2^−ΔΔCt method and all data are expressed as mean ± standard deviation (SD). Statistical analysis was performed using GraphPad Prism 9.0 software. Table 1 displayed the sequence list of primer pairs.

Transwell assay

Migration and invasion of H1975 and H2009 cell lines were detected by performing Transwell assay. Cells (5 × 10⁴) were inoculated onto chambers coated (for invasion) or uncoated with Matrigel (BD Biosciences, Franklink Lakes, NJ, USA) (for migration). The upper layer and the lower layer were added with serum-free medium and complete DMEM medium, respectively. After incubation for 24 h, migrated or invaded cells were fixed using 4% paraformaldehyde (Cat. # P1110, Solarbio, China) and stained by 0.1% crystalline violet (Cat. #G1063, Solarbio, China). A light microscope (Thermo Fisher Scientific, Waltham, MA, USA) was employed for counting cell numbers.

Cell viability

CCK-8 kit (Cat. #CK-04, Dojindo Laboratories, Kumamoto, Japan) was conducted for measuring cell viability. Differentially treated cells were cultured at 1 × 10³ cells/well in 96-well plates and added with CCK-8 solution at specific time points for 2-h incubation at 37 °C. Cell viability was indicated by the O.D 450 value of each well measured by a microplate reader.

Statistic analysis

All the statistical data were analyzed using R language (version 4.3.1) or GraphPad Prism 9.0. Specific analysis methods and software packages were mentioned in the above contents. Sangerbox (http://sangerbox.com/) offered auxiliary analysis in this article.

The correlation between clinical features and mRNAsi in LUAD patients

We first looked at the correlation between mRNAsi and clinical characteristics of LUAD patients. There was no significant correlation between mRNAsi and stage (P = 0.064), N. stage (P = 0.29) and M. stage (P = 0.076) (Figs. 1B, 1D and 1E) but mRNAsi was significantly related to gender (P = 0.00027) and T. stage (P = 0.019) (Figs. 1A and 1C). In addition, compared to those with a high mRNAsi, patients with low mRNAsi tended to show a better prognosis (Fig. 1F), which provides important support for the application of mRNAsi-based prognostic prediction models in LUAD.

Constructing a weighted correlation network using mRNAsi-related DEGs

A total of 1,929 DEGs (1,402 down-regulated and 527 were up-regulated) were visualized in the volcano plot (Fig. 2A). β = 5 was selected to analyze the co-expression networks of different soft threshold capabilities (Figs. 2B and 2C). Gene dendrogram and module colors distinguished six modules, among which grey module was a gene set that cannot be classified into other modules (Fig. 2D).

Analysis on the correlation between each module and mRNAsi showed that the turquoise module was the one the most significantly correlated with mRNAsi (Fig. 2E). A total of 451 genes with high MM and GS were screened from the 775 module genes significantly correlated with mRNAsi in the turquoise module (Fig. 2F, MM > 0.4, GS > 0.4). These 451 genes were then used for subsequent studies.

Construction of an mRNAsi-related RiskScore model

A total of 345 genes showing greater prognostic impact were filtered from the 451 gene modules using univariate cox regression analysis. These genes were reduced by LASSO, and a mutual increase between lambda and the number of independent variable coefficients approaching 0 could be observed (Fig. 3A). A model was created by 10-fold cross-validation. The model was the optimal when lambda = 0.0458, under which a total of eight genes (GNG7, EIF5A, ANLN, FKBP4, GAPDH, GNPNAT1, E2F7 and CISH) were considered as the target genes in the following study (Figs. 3A–3B). The eight genes were subjected to multivariate Cox analysis and the risk coefficient for each gene was calculated (Fig. 3C).

Finally, our risk model was designed: Risk score = −0.108*GNG7 + 0.154*EIF5A + 0.051*ANLN + 0.08*FKBP4 + 0.116*GAPDH + 0.229*GNPNAT1 + 0.06*E2F7 −0.042*CISH.

The KM survival curve showed a significantly lower survival probability of the high RiskScore group in TCGA than patients with a low risk score (P < 0.0001). The model had a high AUC value, as observed from 1-, 3- and 5-year prognostic ROC curves of the TCGA dataset (Fig. 3D). In TCGA dataset, the cutoff of risk score was used to divide high and low risk types. We found that survival was significantly better in high-risk patients than in low-risk patients. The expression of GNG7 and CISH showed a negative correlation with risk score, while that of six genes (EIF5A, ANLN, FKBP4, GAPDH, GNPNAT1, E2F7) was positively correlated with risk score (Fig. 3E).

Validation was further performed on LUAD datasets of GSE30219 and GSE19188 to assess the model’s robustness. Similar results in the two validation cohorts as in the training set were observed, with the ROC curves of both datasets having high AUC values. High risk score could predict a worse prognosis than a better prognosis of low risk score (Figs. 3F–3G).

Identified and enrichment of DEGs between high and low risk groups

In the TCGA dataset, a total of 400 DEGs (179 down-regulated and 221 up-regulated genes) between the two groups were found (Fig. 4A).

The DEGs showing high and low risk score were subjected to KEGG and GO analyses. In GO analysis, we found that the top 10 biological processes and cellular components were all related to metabolism, cell cycle, and chromosome activity (Figs. 4B–4C).

Immune microenvironment differences between patients with high and low risk score

In the TCGA dataset, 13 types of immune cells in the low group with a better prognosis showed higher scores, while six types of immune cells in the high risk group with a worse prognosis had higher scores (Fig. 5A). Correlation analysis revealed a positive correlation between the risk score and infiltration of most immune cells (Fig. 5B). Innate and acquired immunity was better in the low risk group (Fig. 5C). We used ESTIMATE to assess the overall immunity score of the TCGA data, with the low risk group showing higher immunity scores (Fig. 5D). The MCP-counter immune score indicated that most of the immune cells with differences also scored higher in the low risk group (Fig. 5E).

Patients’ sensitivity to chemotherapy and immunotherapy assessed by risk score

The low risk group in the TCGA cohort showed lower TIDE scores, which indicated that those patients could benefit from taking immunotherapy (Fig. 6A). The dysfunction score was lower in the high risk group, while the exclusion score was lower in the low risk group (Figs. 6B–6C). At the same time, we computed the correlation between risk score and TIDE, dysfunction and exclusion, and the results indicated that our risk score was significantly positively correlated with TIDE and exclusion but significantly negatively correlated with dysfunction (Figs. 6D–6F).

Following the analysis on the correlation between drug sensitivity and risk score, GSK2256098C, Carmustine, Dacarbazine, PCI-34051 and Elephantin were selected as the sensitive drugs for LUAD patients (Fig. 6G). Patients in high risk group were more sensitive to Elephantin and GSK2256098C (Fig. 6H).

Development and verification of a clinical prediction nomogram based on mRNAsi-related risk model

As shown in Fig. 7A, we constructed a survival prediction decision tree based on the risk model combined with clinical features (Stage) and evaluated the overall survival of LUAD patients by different risk groupings (C1, C2, and C3). Overall survival varied significantly in the three risk subgroups, with the highest survival probability in C1 and the lowest probability in C3 (Fig. 7B). Among them, all the patients in risk subgroup C2 were high risk patients, while all patients in group C1 were low risk patients (Fig. 7C). Patients’ survival status varied significantly in different risk subgroups, with the most favorable survival status in subgroup C1 and the worst survival status in subgroup C3 (Fig. 7D). The risk score was validated as the most significant prognosis factor (Figs. 7E–7F). Moreover, a nomogram integrating risk score and other clinicopathological characteristics was developed, and here risk score once again showed the strongest influence on evaluation LUAD survival (Fig. 7G).

Calibration curve reflected the prediction accuracy of the model. The predicted calibration curves for 1, 3, and 5 year(s) were close to the standard curve, which indicated that the nomogram could predict correctly (Fig. 7H). Additionally, significantly higher risk score and nomogram benefits than the extreme curve was observed, as shown in the DCA. Both the risk score and nomogram showed the strongest survival prediction (Fig. 7I). This further supports that our constructed risk model combined with clinicopathological features can effectively predict the survival prognosis of LUAD patients and provide potential clinical guidance value for individualized treatment.

Validation of the prognosis model using wet experiments

We experimentally validated the makers involved in the risk model using PCR experiments. The results of qRT-PCR indicated that the expression of GNG7, GNPNAT1 and CISH was downregulated in LUAD cell lines (H1975 and H2009). EIF5A, ANLN, FKBP4, E2F7 were upregulated in H1975 and H2009 cell lines, while GAPDH was not differentially expressed in BEAS-2B, H1975 and H2009 cells (Fig. 8A). We further explored the functions of the EIF5A due to its highest expression level in LUAD cells. The invasive and migratory abilities of LUAD cells were significantly decreased after suppressing EIF5A expression in H1975 and H2009 cells, indicating that EIF5A may promote LUAD progression (Figs. 8B–8C). Finally, we preliminary explored the effect of therapeutic drug GSK2256098. The mRNA expression of EIF5A was noticeably decreased in H1975 and H2009 cells after treatment with GSK2256098 at the concentration of 10 µM (Fig. 8D). Also, cell viability was significantly affected following the treatment of GSK2256098 at 10 µM (Fig. 8E). These results further indicate that the key genes identified based on mRNAsi characterization may have important biological functions in LUAD and also provide a rationale for their use as potential therapeutic targets.