Immune- and Stemness-Related Genes Revealed by Comprehensive Analysis and Validation for Cancer Immunity and Prognosis and Its Nomogram in Lung Adenocarcinoma

We downloaded the data from The Cancer Genome Atlas database (TCGA, https://portal.gdc.cancer.gov/↗), which contained transcript data of 535 tumor tissues with LUAD and 59 normal tissues (TCGA-LUAD) and clinical data of 522 patients with LUAD. We also download GSE68465↗ from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/↗). Transcriptome profiling data of 443 tumor tissues with LUAD and 19 normal tissues in the GSE68465↗ dataset were used for further analysis.

Using mRNA stemness index (mRNAsi) as a variable, the R packages “survminer” was applied to analyze the correlation of mRNAsi with clinical parameters. Then, according to the median of mRNAsi, the tumor components were separate into two groups (high–mRNAsi level and low–mRNAsi level group) for survival analysis. The mRNAsi index of LUAD was acquired from the supplemental information in the study of Malta et al. (14).

After the analysis above, we first assessed the difference between tumor and normal group according to mRNAsi, and then we used the R package “limma” (the Wilcoxon test) to screen the differential expression genes (DEGs) related to LUAD. The DEGs were next used to construct a coexpression module using a weighted gene coexpression network analysis (WGCNA). The construction process includes the following main steps: (1) give a definition for similarity matrix; (2) use the function pickSoftThreshold to select the soft threshold powerβ; (3) convert the adjacency matrix into a topological overlap matrix (TOM); (4) execute hierarchical aggregation of dissTOM derived from TOM; (5) from the hierarchical clustering tree, use the dynamic tree cutting method to distinguish modules with identical expression profiles; (6) quantify the coexpression similarity of the entire modules and compute their characteristic genes, etc. (19). At last, we selected two modules with the highest absolute value associated with mRNAsi (including positive and negative correlations) for the following analysis.

For better understanding the functions and pathways of the two mRNAsi-related modules above in LUAD, each of them was analyzed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEG)G enrichment, respectively. R package “colorspace”, “stringi”, and “ggplot2” were used. The GO enrichment analysis included three components: molecular function (MF), cellular component (CC), and biological process (BP). Choose the threshold as p-value <0.05 and q-value <0.05.

To discover the immune-related genes (IRGs) in LUAD, we first downloaded the IRG data from the immunology database and analysis website (ImmPort, https://www.immport.org/↗). By taking the intersection with the DEGs that we screened before, we extracted the LUAD immune-related DEGs for the next step. We further analyzed the intersection of mRNAsi-related DEGs and immune-related DEGs via Venn diagram.

Using the R package “survival”, we further screened for prognosis-related hub genes by univariate Cox regression analysis. We selected the genes with P < 0.05 and HR ≠ 1 from the univariate Cox analysis. The Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis is a popular algorithm, which was extensively utilized in medical studies (20, 21). Using the R package “glmnet” and “survival”, the optimal model based on prognosis-related stem cell index–related differential genes and immune-related differential genes (SCIRGs) was subsequently identified utilizing LASSO regression analysis (22). The model formula is

where Coef refers to the regression coefficient of SCIRGs in LASSO Cox regression analysis, “Ni” is the expression value of the gene, and “n” is the number of SCIRGs.

To verify the predictive ability of the model, we assessed the model through the training set (TCGA-LUAD) and the validation set (GSE68465↗), respectively. Using R package “survival” and “survminer” for survival analysis, we drew a Kaplan–Meier curve in TCGA and GEO datasets, separately. To explore high- and low-risk hub genes in the model and the risk score distribution in LUAD, we used the “pheatmap” package to depict risk curves, survival status maps, and risk heat maps. In addition, using R package “survival” for an independent prognostic analysis of the training and validation set, these helped us to understand whether the risk score can be used as a prognostic factor independent of clinical parameters. We used R package “survivalROC” to draw a multi-index ROC curve to assess prediction accuracy of the model. Last, we further take risk score with clinical parameters to draw a nomogram. The clinical parameters included age, gender, TNM (TNM Classification of Malignant Tumors, UICC 8th edition), and stage. The nomogram was used to evaluate the 1-, 2-, and 3-year survival rates of patients. The predictive capability of the model was assessed by calculating the C-index and plotting the calibration curves.

Gene set enrichment analysis (GSEA) and the frequency of gene mutations were analyzed in high- and low-risk groups by utilizing the Maftools package of R. Furthermore, the association of high- and low-risk groups with TIME was also validated. CIBERSORT (https://cibersort.stanford.edu/↗) was used to input the data and perform 1,000 iterations to explore the 22 immune cells’ proportions. In addition to this, we compared the difference of 22 immune cells’ related function between the two subgroups. The correlation between risk score and common oncogene (EGFR, ALK, ROS1, KRAS, and TP53), CD274, and stem cell index (DNAss and RNAss) were also explored. Then, correlation analysis was performed between tumor mutation burden (TMB) and risk score, and the difference analysis between TMB and the subgroups was also explored. We also explored the distribution of every samples classified by clinical parameters between the subgroups, and the difference of stage and immunophenotyping was further demonstrated. Finally, the drug sensitivity analysis of every hub gene was demonstrated using CellMiner.

To gain a deep understanding of the key genes in the model, we use Oncomine (https://www.oncomine.org↗), UALCAN (http://ualcan.path.uab.edu↗), Kaplan–Meier (http://kmplot.com↗), TIMER (https://cistrome.shinyapps.io/timer/↗), GEPIA (https://gepia.cancer-pku.cn↗), and other web-based bioinformation tools to perform differential analysis, survival analysis, immune infiltration analysis, and correlation analysis in LUAD. In the correlation analysis, so as to represent the strength of the interrelationship between gene expression and tumor immune infiltration in TIMER, we categorized it as follows: 0.00–0.19, “very weak”; 0.20–0.39, “weak”; 0.40–0.59, “moderate”; 0.60–0.79, “strong”; and 0.80–1.0, “very strong”.

The infiltration proportion of every immune cell in the two risk groups is shown, respectively (Figure 7A). Plasma cells, CD8 T cells, activated memory CD4 T cells, M0 and M1 macrophages, and activated mast cells were more abundant in the high-risk subgroup; correspondingly, immune-related function like inflammation-promoting, MHC class I, NK cells, and Tfh were more frequency in the high-risk subgroup. Whereas memory B cells, memory resting CD4 T cells, monocytes, M2 macrophages cell, resting dendritic cells, resting mast cells, etc., were more abundant in the low-risk subgroup; correspondingly aDCs, B cells, DCs, HLA, mast cells, etc., were more common in the low-risk group (Figures 7B, C). The association of risk score with TIME shows that both immune score and stromal score were negatively relevant to risk score (Figures 7D, E).

Furthermore, to prove the relevance of these genes to immunity, we compared the correlation between hub genes and immune cells through TIMER. ADRB2 has moderate correlation with dendritic cell; CX3CR1 has moderate correlation with macrophage and neutrophil; IL3RA also has moderate correlation with neutrophil and dendritic cell. Apart from these genes, other genes also have weak correlation with immune cell (Figures 7F–I; Supplementary Figure 3).

The relationship between risk score and TMB was further probed. The results exhibited that TMB was markedly higher in the high-risk group than in the low-risk group, and the higher the risk score, the larger the TMB (R = 0.31, p = 4e−12; Figures 8A, B). To explore the difference of every samples classified by clinical parameters between those two groups, clinical relevance heat map was used. In Figure 8C, age and T staging have markedly difference between the two groups. We also found that the proportion of stage IV samples has almost equal distributions between the two groups, and there were more samples in the high-risk subgroup and fewer samples in the low-risk subgroup in stages II–III, but there was an opposite result in stage I (p = 0.003, chi-square test) (Figure 8D). Then, 446 TCGA samples were further classified according to immune subtype. As shown in Figure 8E, there were more C3 subtypes in the SCIRG-low subgroup, whereas more C1 and C2 subtypes in SCIRG-high subgroup (p = 0.001, chi-square test).

Finally, the SCIRG gene was analyzed in combination with drug sensitivity, and the first 16 drugs with statistically significant differences were selected. The results uncovered that the expression level of CX3CR1 was positively relevant to the sensitivity of Alectinib, LDK-378, Denileukin Diftitox Ontak, Estramustine, Nelfinavir, PF-06463922, and Carmustine. This indicated that the higher the expression of CX3CR1, the stronger the sensitivity to the abovementioned drugs. We also ascertained that the CX3CR1 expression level was negatively relevant to the sensitivity of Irofulven. In addition, CRABP1 expression level was positively relevant to the sensitivity of Bendamustine and Dexrazoxane. The expression of MET was negatively relevant to the sensitivity of Bendamustine and Dexrazoxane. The expression of LIFR and BDNF was negatively relevant to the sensitivity of Tamoxifen, whereas the expression of GPER1 was positively relevant to Procarbazine (Supplementary Figure 4).

We then implement GSEA analysis to find out in which function the two subgroups of genes were up- or downregulated. For example, the genes of the high-risk group were upregulated in chromosome segregation, cornification, DNA dependent, DNA replication, and epidermal cell differentiation, whereas they also upregulated in cell cycle, DNA replication, proteasome, pyrimidine metabolism, and spliceosome in KEGG. This tells us that the high-risk group was mainly correlated with proliferation in LUAD. On the other hand, the genes of the low-risk group were downregulated in cilium movement, rDNA heterochromatin assembly, ciliary plasma, cilium, and DNA packaging complex in GO, whereas they also have the same performance in asthma, Fc epsilon RI signaling pathway, long-term depression, systemic lupus erythematosus, and vascular smooth muscle contraction in KEGG (Supplementary Figures 5A–D, p < 0.05). Next, we analyzed gene mutations to gain further insight in the charicteristics of the subgroups. We found 96.03% samples were altered in high-risk groups, whereas 80.58% samples were altered in low-risk groups. Missense variations were the most common mutation type (Supplementary Figures 5E, F). In addition to this, the risk score was also markedly relevant to common oncogenes expression, such as ALK, ROS1, KRAS, and TP53, but no markedly relevant to CD274 (Supplementary Figure 6).

Here, to find out potential biomarkers in LUAD, we explored the 16 genes in the model through bioinformatics. First, Oncomine was used to explore the overall difference of the above genes in lung cancer, and UALCAN was utilized to seek every gene differential expression between LUAD and normal. As shown in Supplementary Figures 7, 8, whether in Oncomine or in UALCAN, ADRB2, ARRB1, BDNF, etc., had a lower expression in tumor group than normal tissues, whereas CBLC, GPI, and PAK1 had a higher expression in tumor group than normal tissues. There were no studies of ANGPTL4, CRABP1, MET, and SEMA3A in Oncomine, but in UALCAN, they were frequently expressed in tumor group than normal tissues. Second, the survival analysis of key genes in LUAD was analyzed by the Kaplan–Meier method. The results exhibited that the high expression groups of ANGPTL4, CBLC, CRABP1, etc., have a poorer prognosis in LUAD than the low expression group, whereas the survival analysis of ADRB2, CX3CR1, GPI, etc., in LUAD exhibited that the prognosis was better in the high expression group (Supplementary Figure 9). To further understand the correlation between these key genes and common oncogenes such as TP53, EGFR, and CD274, we explored the correlation through GEPIA (spearman, P value < 0.05 and R > ± 0.1). Genes related to TP53 include ANGPTL4, ARRB1, CBLC, etc.; genes related to EGFR include ADRB2, ANGPTL4, ARRB1, etc; genes related to CD274 include ADRB2, BDNF, CBLC, etc. (Supplementary Figures 10–12).

We respectively compared the expression of the SCIRGs in normal lung epithelial cells, lung cancer cells, and lung CSCs and repeatedly compared them in different cell lines (A549 and HCC827) (Figure 9; Supplementary Figure 13). The results showed that in the A549 cell line, the expression results of nine genes in the 16 genes were consistent with those in UALCAN (Figure 9; Supplementary Figures 8, 13): ADRB2, ANGPTL4, BDNF, CBLC, CRABP1, CX3CR1, GPI, IL3RA, and SCH3; in the HCC827 cell line, the expression results of nine genes among the 16 genes were consistent with those in UALCAN: ADRB2, ANGPTL4, BDNF, CBLC, CX3CR1, IL3RA, LIFR, and MET. Therefore, half or more of the genes in our model were consistent with the gene expression results of external data.