Identification of tumor heterogeneity associated with KRAS/TP53 co-mutation status in lung adenocarcinoma based on single-cell RNA sequencing

Abstract

Lung cancer stands as the predominant cause of cancer-related mortality globally. Lung adenocarcinoma (LUAD), being the most prevalent subtype, garners extensive attention due to its notable heterogeneity, which significantly influences tumor development and treatment approaches. This research leverages single-cell RNA sequencing (scRNA-seq) datasets to delve into the impact of KRAS/TP53 co-mutation status on LUAD. Moreover, utilizing the TCGA-LUAD dataset, we formulated a novel predictive risk model, comprising seven prognostic genes, through LASSO regression, and subjected it to both internal and external validation sets. The study underscores the profound impact of KRAS/TP53 co-mutational status on the tumor microenvironment (TME) of LUAD. Crucially, KRAS/TP53 co-mutation markedly influences the extent of B cell infiltration and various immune-related pathways within the TME. The newly developed predictive risk model exhibited robust performance across both internal and external validation sets, establishing itself as a viable independent prognostic factor. Additionally, in vitro experiments indicate that MELTF and PLEK2 can modulate the invasion and proliferation of human non-small cell lung cancer cells. In conclusion, we elucidated that KRAS/TP53 co-mutations may modulate TME and patient prognosis by orchestrating B cells and affiliated pathways. Furthermore, we spotlight that MELTF and PLEK2 not only function as prognostic indicators for LUAD, but also lay the foundation for the exploration of innovative therapeutic approaches.

Keywords: B cell; KRAS; Lung adenocarcinoma; TP53; mutation; prognosis; single-cell RNA sequenceing.

Figure 1

Cell type analysis and classification in KRAS/TP53 mutation. A. The t-SNE diagram shows the results of the dimensionality reduction analysis. B. Dot diagram displays the average expression of marker genes in 26 cell clusters. C. The t-SNE diagram shows ten cell types, namely T cells, epithelial cells, macrophage, monocytes, DC, neutrophils, B cells, endothelial cells, NK cells and Tissue stem cells. D. Cell proportions of ten cell types in KRAS/TP53 MUT group and WT group.

Figure 2

B cell subtypes and interactions in KRAS/TP53 mutant context. A. The UMAP diagram shows the situation of ten B cell subtypes in the KRAS/TP53 MUT group and WT group. B. The proportion of ten B cell subtypes in the KRAS/TP53 MUT group and WT group. C and D. Number of interactions and interaction weights/strength for each cell type. E. The heat map shows the outgoing and incoming signal pattern of each cell type recognition cell in the KRAS/TP53 MUT group.

Figure 3

Gene module analysis and correlation in B cells with KRAS/TP53 mutations. A. hdWGCNA lineage dendrogram of B cells. B. Optimal soft thresholds were selected and maximum, median and average connectivity are shown. C. Heat map showing correlation analysis between three gene modules. D. Three gene modules were obtained according to the standard process and the top 10 huh genes were presented accordingly. E. Distribution of the three gene modules in B cells in the KRAS/TP53 MUT group.

Figure 4

Univariate regression analysis and survival outcomes in LUAD patient clusters. A. Volcano plot shows the results of univariate regression analysis. Red represents up-regulated genes, blue represents down-regulated genes and gray represents genes with no significant difference. B. Matrix heatmap shows LUAD patients divided into two clusters (k = 2). C and D. CDF is displayed in delta plot and cumulative distribution curve plot. E. Kaplan-Meier curve show OS for two LUAD clusters.

Figure 5

Comparative analysis of immune infiltration and functionality between clusters. A. Degree of immune cell infiltration between clusters. B. ImmuneScore, StromalScore, and ESTIMATEScore between clusters. C. Expression of immune-related genes between clusters. D-I. Differences in TIDE, IFNG, dysfunction score, exclusion score, TAM M2 and MDSC proportions between clusters (*P < 0.05, **P < 0.01 and ***P < 0.001).

Figure 6

Prognostic gene identification and risk assessment in LUAD via regression and LASSO model. A. Volcano plot showing the results of univariate Cox regression analysis. B and C. Coefficient distribution diagram of each gene and optimal lambda for constructing LASSO model. D. Cox coefficient of seven prognostic genes. E. The Kaplan-Meier analysis results of high- and low-risk groups based on TCGA-LUAD dataset. F. The ROC analysis based on TCGA-LUAD dataset.

Figure 7

Comprehensive evaluation and validation of prognostic nomogram. A and B. Univariate and multivariate Cox analyzes assessed the independence of the predictive risk models. C. Construct the prognostic nomogram combining risk scores and clinical characteristics. D. Calibration curves show the correlation between the nomogram’s 1-, 3-, and 5-year predictions and actual observations. E. Standardized net benefit decision curve graph. F. ROC analyzes the AUC of the nomogram, the risk score and each clinical feature.

Figure 8

Correlation of prognostic genes with immune and stromal components. A. Correlation between prognostic gene expression levels and ImmuneScore. B. The heatmap shows correlation of prognostic genes with ImmuneScore, StromalScore and ESTIMATEScore. C. The heatmap shows the correlation of prognostic genes with 22 types of immune-related cells (*P < 0.05, **P < 0.01, ***P < 0.001). D. The heatmap shows the correlation of prognostic genes with 10 types of immune-related cells.

Figure 9

Evaluating a risk model’s predictive accuracy for treatment response and survival. A. Kaplan-Meier analysis based on IMvigor210 dataset. B. Differences in risk scores for different clinical response groups based on the IMvigor210 dataset. C. Proportions of high- and low-risk in different clinical response groups based on the IMvigor210 dataset. D. Based on the IMvigor210 dataset, Kaplan-Meier analysis was performed on the risk scores of clinical stages I and stages II patients. E. Based on the IMvigor210 dataset, Kaplan-Meier analysis was performed on the risk scores of clinical stages III and stages IV patients. F. Kaplan-Meier analysis based on GSE78220 dataset. G. Differences in risk scores for different clinical response groups based on the GSE78220 dataset. H. Proportions of high- and low-risk in different clinical response groups based on the GSE78220 dataset.

Figure 10

Cell type scoring, interaction, and trajectory analysis in KRAS/TP53 MUT and WT groups. A. The Violin plot shows the prognostic genesets scored for ten cell types in the KRAS/TP53 MUT group and WT group based on different algorithms. B. The Bubble chart shows the scoring cluster markers of ten types of cells by different algorithms. C. The Violin plot shows the scoring of ten cell types on the prognostic geneset between the KRAS/TP53 MUT group and the WT group (*P < 0.05, **P < 0.001, ***P < 0.0001 and ****P < 0.00001). D and E. The number of interactions between cells and the proportion of interactions. F. The heatmap shows the signal flow pattern of mutual recognition cells between cells in the KRAS/TP53 MUT group. G and H. Trajectories showing pseudo-time-dependent cellular states of epithelial cells in the IKRAS/TP53 MUT group. I. The heatmap shows the expression of seven prognostic genes over pseudo-time.

Figure 11

Impact of MELTF and PLEK2 knockdown on protein levels, proliferation, invasion, and migration in A549 and H1299 cells. A. Protein levels of MELTF and PLEK2 in A549 and H1299 cells after siRNA interference respectively. B. Proliferation rate of A549 and H1299 cells post MELTF and PLEK2 interference using the CCK8 assay. C. Invasive capability of A549 and H1299 cells after MELTF knockdown via wound-healing assay. D. Invasive capability of A549 and H1299 cells after PLEK2 knockdown via wound-healing assay. E. Migration ability of A549 and H1299 cells post MELTF interference in the transwell assay. F. Statistical analysis of the migration ability of A549 and H1299 cells post PLEK2 interference in the transwell assay.