Multi‑omics identification of a signature based on malignant cell-associated ligand-receptor genes for lung adenocarcinoma

Abstract

Purpose: Lung adenocarcinoma (LUAD) significantly contributes to cancer-related mortality worldwide. The heterogeneity of the tumor immune microenvironment in LUAD results in varied prognoses and responses to immunotherapy among patients. Consequently, a clinical stratification algorithm is necessary and inevitable to effectively differentiate molecular features and tumor microenvironments, facilitating personalized treatment approaches.

Methods: We constructed a comprehensive single-cell transcriptional atlas using single-cell RNA sequencing data to reveal the cellular diversity of malignant epithelial cells of LUAD and identified a novel signature through a computational framework coupled with 10 machine learning algorithms. Our study further investigates the immunological characteristics and therapeutic responses associated with this prognostic signature and validates the predictive efficacy of the model across multiple independent cohorts.

Results: We developed a six-gene prognostic model (MYO1E, FEN1, NMI, ZNF506, ALDOA, and MLLT6) using the TCGA-LUAD dataset, categorizing patients into high- and low-risk groups. This model demonstrates robust performance in predicting survival across various LUAD cohorts. We observed distinct molecular patterns and biological processes in different risk groups. Additionally, analysis of two immunotherapy cohorts (N = 317) showed that patients with a high-risk signature responded more favorably to immunotherapy compared to those in the low-risk group. Experimental validation further confirmed that MYO1E enhances the proliferation and migration of LUAD cells.

Conclusion: We have identified malignant cell-associated ligand-receptor subtypes in LUAD cells and developed a robust prognostic signature by thoroughly analyzing genomic, transcriptomic, and immunologic data. This study presents a novel method to assess the prognosis of patients with LUAD and provides insights into developing more effective immunotherapies.

Keywords: Immunotherapy; Lung adenocarcinoma; Prognostic model; Single-cell; Tumor microenvironment.

Fig. 1

The workflow illustrating the schematic overview of single-cell sequencing and GSE171145 dataset analysis (upper) and prognostic model establishment (lower)

Fig. 2

Definition of cell clusters. (A) The t-distributed stochastic neighbor embedding (t-SNE) plot of nine samples in the GSE171145 dataset, colored to indicate sample names. (B) The t-SNE plot of the distribution of 27 clusters, colored to indicate cell clusters. (C) The t-SNE plot of eight cell types after cell annotation, colored to indicate cell types. (D) Dot plots of the top five marker genes contributing to the clusters, x-axis: cell types, y-axis: marker genes, dot colors: average expression (blue represents low expression and red represents high expression), and dot size: percent expressed cells in the cluster. (E and F) Numbers and proportions of cell types in each sample after annotation, x-axis: cell numbers and proportions and y-axis: cell types. (G) The t-SNE plot of aneuploid and diploid cells, colored to indicate aneuploid and diploid cells

Fig. 3

Transcription factor regulatory networks in malignant tumor subpopulations and trajectory analysis of malignant cells in lung adenocarcinoma. (A) Heatmap with regulon area under the curve (AUC) matrix of scaled AUC values (columns) detected in different cell types (rows). Blue represents low expression, yellow represents moderate expression, and red represents high expression. (B) Density map of steady-state cells. The darker color represents more steadiness. (C) Heatmap of transcriptional regulatory activity (columns) of nine cell types (rows). Blue represents low expression and red represents high expression. (D) Point plots of the top five regulon specificity scores. X-axis: rank and y-axis: regulon specificity scores. (E and F) Monocle 2 trajectory plots showing state dynamics and pseudotime curves. Each dot represents a singlet and the color gradient represents the pseudotemporal order. States 1–3 are labeled in the same topology. (G) Heatmap hierarchical clustering of differentially expressed transcription factor genes (columns) along the pseudotime curve (rows). Blue represents low expression, gradient represents moderate expression, and red represents high expression

Fig. 4

Cell–cell communication analysis and identification of molecular subtypes. (A and B) Circle plots showing the number and strength of cell type interactions. The ligand–receptor expressed by each cell type, the thicker the lines, the greater the number/intensity of ligand–receptor. Dot size represents the number of cells in the subpopulation. (C) Enrichment of tumor microenvironment-related pathways inputs and outputs among cell types. (D) Hazard ratio distribution plot for univariate Cox analysis of malignant cell ligand–receptor-related gene sets. X-axis: cox coefficient and y-axis: −log10(p-value), colored to indicate cell states. (E) Cumulative distribution function (CDF). X-axis: consensus index and y-axis: CDF, colored to indicate clustering number. (F) Delta area curve for The Cancer Genome Atlas cohort samples. X-axis: k and y-axis: relative change in area under CDF curve. (G) Heatmap of sample clustering when k = 2. (H) Kaplan–Meier survival analysis comparing the prognosis of two subtypes. X-axis: years and y-axis: survival probability

Fig. 5

Immune infiltration analysis in molecular subtypes and differential expression of malignant cell-associated ligand–receptor genes (A) Relative abundance of immune cells infiltrating the tumor microenvironment between molecular subtypes, x-axis: infiltrating immune cells and y-axis: score, colored to indicate different cell clusters, red, cluster1; green, cluster2. (B) Differences in stromal, immune, and “ESTIMATE” scores in molecular subtypes. X-axis: immune scores and y-axis: score, colored to indicate different cell clusters, red, cluster1; green, cluster2. (C) Expression levels of 47 immune checkpoints between molecular subtypes. X-axis: genes and y-axis: expression, colored to indicate different cell clusters, red, cluster1; green, cluster2. (D) Differences of TIDE, IFNG, MDSC, Exclusion, Dysfunction, and TAM.M2 in molecular subtypes. X-axis: cluster and y-axis: immune suppressive score, colored to indicate different cell clusters, red, cluster1; green, cluster2. (E) The volcano plot of differentially expressed genes was identified between cluster1 and cluster2 (false discovery rate [FDR] < 0.05). X-axis: log2(FoldChange), y-axis: −log10(FDR), color of bubbles: red, considerably upregulated, and blue, considerably downregulated. (F) Bar chart of top five terms showing pathway enriched in biological process, cellular component, and molecular function. X-axis: gene counts in the enriched pathway and y-axis: pathway, colored to indicate enriched − log10(p-value). (G) Top 15 terms of the Kyoto Encyclopedia of genes and genomes (KEGG) pathways enrichment visualized via a bubble chart. X-axis: gene ratio in the enriched pathway and y-axis: pathway, colored to indicate enriched − log10(p-value), and the bubble size indicates the count of enriched genes

Fig. 6

Identification of prognostic gene signature. (A) 95 predictive models using diverse machine learning techniques, employing a tenfold cross-validation method. The C-index for each model was computed, covering both the TCGA-LUAD and GSE72094 cohorts. (B) Lambda trajectory of differentially expressed genes. X-axis: −In (lambda) and y-axis: coefficients, colored to indicate genes. (C) Confidence interval under lambda. X-axis: In (lambda) and y-axis: partial likelihood deviance, colored to indicate genes. (D) Kaplan–Meier survival analysis in The Cancer Genome Atlas (TCGA) dataset. X-axis: years and y-axis: survival probability. (E) Receiver operator characteristic (ROC) curve analysis-based evaluation of the prediction performance of gene signature in TCGA. X-axis: false positive fraction and y-axis: true positive fraction, colored to show time site. (F) Kaplan–Meier survival analysis in GSE31210. X-axis: years and, y-axis: survival probability. (G) ROC curve analysis-based evaluation of the prediction performance of gene signature in GSE31210. X-axis: false positive fraction and y-axis: true positive fraction, colored to show time site. (H) Pie plot of the difference in clinical characteristics between high- and low-risk groups (Wilcox test, *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001)

Fig. 7

Immune infiltration analysis and drug sensitivity analysis in high- and low-risk groups. (A) Comparison of 28 immune cell scores in high- and low-risk groups. (B) Comparison of immune checkpoint expression in high-and low-risk groups. (C) Analyzing the association between IC50 values and the risk scores in patients with lung adenocarcinoma. (D-G) Analysis of correlation and differences in sensitivity to drugs among potential medications derived from the CTRP and PRISM datasets

Fig. 8

Construction of nomogram. (A) Univariate Cox regression analysis of LUAD patients. (B) Multivariate Cox regression analysis of LUAD patients. (C) AUC analysis of risk score, nomogram, stage, T stage, and N stage. (D) Calibration curve of the nomogram. (E) Decision curves of “risk score”, “nomogram”, “T stage”, “N stage”, “stage”, “all”, and “None” models. (F) Nomogram for predicting the 1-, 3-, and 5-year survival rates based on the risk score

Fig. 9

MYO1E promotes proliferation and migration of LUAD cells. (A) RT-qPCR analyse confirmed MYO1E knockdown in A549 and H1299 cells using two siRNAs (B) Colony formation of A549 cells and H1299 cells transfected with control or si-MYO1E was measured by ImageJ. (C) Edu assay assessed the cell proliferation of control cells compared to MYO1E knockdown cells. (D) Transwell assay demonstrated the cell migration of control cells compared to MYO1E knockdown cells. (E) Wound healing assay showed the cell migration of control cells compared to MYO1E knockdown cells