Identification of Mitotic Catastrophe-related Prognostic Genes in Lung Adenocarcinoma via Transcriptome Analysis

Summary

This study identifies eight mitotic catastrophe-related prognostic genes (PDGFB, LDHA, ZEB2, H2AX, FKBP4, DMD, ANXA2, S100B) in lung adenocarcinoma, establishes a risk model with acceptable predictive performance, and reveals distinct immune microenvironment characteristics between risk groups.

Abstract

Lung adenocarcinoma (LUAD), the most common lung cancer subtype, exhibits a poor prognosis. Although mitotic catastrophe can be induced in LUAD, leading to cell death, the prognostic and biological significance of mitotic catastrophe-related genes (MCRGs) in this disease remains unclear. Differentially expressed genes (DEGs) associated with LUAD were screened from the TCGA-LUAD cohort, and overlapping genes were defined by intersecting DEGs with MCRGs. Candidate prognostic genes were selected via machine learning, after which a multigene prognostic model was constructed. An eight-candidate prognostic signature comprising PDGFB, LDHA, ZEB2, H2AX, FKBP4, DMD, ANXA2, and S100B was established, demonstrating acceptable predictive reliability in training with 2-year, 3-year and 5-year area under the curves (AUCs) of 0.723, 0.715, and 0.625, respectively, in the training set (bootstrap-corrected C-index: 0.699, 95% CI: 0.657-0.742) and consistent validation in GSE31210 (AUCs: 0.812, 0.766, and 0.814). Low-risk patients exhibited significantly better survival outcomes, enhanced immune infiltration, higher immune and stromal scores, increased immunophenoscores, and more active cancer immunity cycles. Significant differences in tumor mutation burden and drug sensitivity were observed between the high- and low-risk groups. Quantitative polymerase chain reaction(qPCR) validated the transcriptomic profiles of six candidate outcome-related genes in clinical samples. A nomogram incorporating the risk score and clinical stage showed good predictive performance. These findings demonstrate that the eight MCRG-related candidate genes have strong potential to predict patient outcomes and stratify risk in LUAD, with distinct immune microenvironment characteristics between risk groups.

Introduction

Lung cancer remains among the most prevalent and deadly cancers worldwide^1,2. Non-small cell lung cancer (NSCLC) is dominated by lung adenocarcinoma (LUAD), which accounts for about 40% of all lung cancer cases^3,4. Although significant strides have been made in operative procedures, molecularly targeted therapies, immune checkpoint inhibitors, and combination regimens, the clinical outlook for LUAD patients remains unfavorable^5,6. Although pathological and imaging-based diagnostic approaches are widely used in LUAD, they have inherent limitations in capturing tumor heterogeneity and molecular characteristics. In contrast, molecular biomarker-based tumor tissue analysis can provide more comprehensive information on tumor biology and may contribute to prognostic assessment and individualized treatment strategies⁷. Therefore, identifying reliable molecular biomarkers remains important for improving prognostic prediction and biological understanding of LUAD.

Mitosis, an essential mechanism of cellular replication, guarantees the faithful partitioning of replicated chromosomes between two daughter cells⁸. Dysregulation of mitotic processes can promote uncontrolled cell proliferation, genomic instability, and tumor progression. Mitotic catastrophe (MC) defines a cell death pathway triggered upon failure of mitosis completion, typically resulting from spindle apparatus damage or defective cell cycle checkpoints⁹. At the morphological level, MC presents with cellular abnormalities including multinucleation, micronuclei, multipolar spindle configurations, and polyploidization¹⁰. As a mechanism for eliminating cells with severe mitotic defects, MC may contribute to the suppression of tumor progression⁹. The p53 deletion in A549 LUAD cells can induce mitotic catastrophe and ultimately cell death¹¹, implying a potential link between MC and LUAD-related cellular processes. Nevertheless, the prognostic significance and biological functions of MC-related genes (MCRGs) in LUAD remain poorly elucidated.

In this research, the prognostic importance and putative biological relevance of MCRGs in LUAD were systematically interrogated utilizing data from TCGA and GEO. Candidate MCRGs were first identified through differential expression and intersection analyses. Multiple machine learning approaches screened candidate genes and established a prognostic risk signature, which was subsequently externally validated. A nomogram incorporating the risk score in conjunction with clinical characteristics was further developed to improve survival prediction. In addition, immune microenvironment characterization, somatic mutation assessment, and chemosensitivity evaluation served to examine the biological implications and clinical utility of the identified genes. Collectively, these findings demonstrate the prognostic potential of MCRGs in LUAD and may inform biomarker development and individualized management.

Protocol

The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Anhui Chest Hospital (K2025-007) on 22 April 2025. Informed consent was obtained from all subjects involved in the study.

Data extraction and normalization 

Transcriptomic profiles and corresponding clinical datasets for LUAD were sourced from the TCGA and GEO cohorts. The TCGA-LUAD dataset was designated as the training set, with GSE72094, GSE31210, and GSE26939 serving as cohorts for external validation (Table 1). In addition, 900 MCRGs were collected from a previous study¹²(Supplementary Table 1). Transcriptome data were annotated using GENCODE v36 or the corresponding GPL platform annotation files. Probe IDs were converted into gene symbols, duplicated genes were merged using the avereps function, and only protein-coding genes were retained to generate gene-level expression matrices. For the TCGA-LUAD training set, genes with Fragments per kilobase of exon model per million mapped fragments (FPKM) < 1 in more than 50% of samples were filtered out, and the remaining expression values were log2 transformed (log2[FPKM+1]). For the GEO validation cohorts, raw expression data were downloaded, probe IDs were mapped to gene symbols using the respective platform annotation files, and multiple probes corresponding to the same gene collapsed by averaging their expression values. These datasets were log2-transformed when necessary. No cross-platform batch effect correction was applied between TCGA and GEO, as we adopted a per-cohort standardization strategy to ensure relative comparability. Specifically, for both the training and validation cohorts, gene expression values were centered and scaled (z-score transformation) using the mean and standard deviation of each dataset individually. The same Cox regression coefficients derived from the training set were then used to calculate risk scores for all cohorts. To maintain clinical applicability and avoid overfitting to any validation set, the median risk score of the training cohort was employed as a fixed cut-off to stratify patients into high and low-risk groups across all external validation cohorts. Clinical information, including age, sex, pathological stage, Tumor-Node-Metastasis (TNM) stage, histological type, survival time, survival status, and tissue type, was extracted when available. The endpoint was overall survival (OS). Samples with incomplete survival information or survival time < 30 days were excluded. Survival time was converted into years, and survival status was coded as 0 for alive and 1 for dead.

Identification and functional analyses of candidate genes

The Limma package identified differentially expressed genes (DEGs) between LUAD tumor and normal samples in the training set¹³. The following criteria defined DEGs: |log₂FC| > 0.5 and adjusted p-value < 0.05. Subsequently, the mfuzz fuzzy clustering algorithm in the R package ClusterGVis was utilized to partition DEGs into distinct expression clusters. Gene Ontology–Biological Process (GO-BP) analysis was completed on the top five representative genes in each cluster based on their membership scores. A set of shared genes was derived by taking the intersection of DEGs with MCRGs. Functional enrichment analysis using Gene Ontology/Kyoto Encyclopedia of Genes and Genomes (GO/KEGG) assessed the biological relevance of overlapping genes. Protein–protein interaction (PPI) networks originated from the STRING database¹⁴. Only interactions with confidence scores > 0.7 remained to improve network reliability.

Prognostic genes screening

The Survival package was used to conduct an univariate Cox regression analysis to identify probable genes linked to overall survival in LUAD¹⁵. Genes with p-value < 0.05 were considered potential prognostic indicators. The TCGA-LUAD training cohort included 500 patients with complete survival data, among whom 216 (43.2%) experienced death events during follow-up. The ratio of candidate genes (n = 108) to events (n = 216) was approximately 1:2, which is acceptable for Cox regression analysis. Subsequently, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis and Extreme Gradient Boosting (XGBoost) model further selected features. Cox proportional hazards models were built with family = "cox" via the cv.glmnet function of the glmnet package. The optimal regularization parameter was identified using 10-fold cross-validation, with λ.min value representing the minimum cross-validation error, which was selected as optimal λ value. Genes with non-zero regression coefficients were extracted as candidate features. For the XGBoost model, survival time and survival status were combined as the outcome variable, with positive values assigned to death events and negative values assigned to censored cases. The parameters were set as objective = "survival: cox" and eval_metric = "cox-nloglik", with 100 iterations and a learning rate of 0.1. After model training, gene importance scores were calculated using feature gain values. The top 20 genes were retained after sorting the importance scores in descending order to reduce feature dimensionality and model complexity. Genes overlapping between the LASSO and XGBoost results were identified as candidate prognostic genes.

Construction and assessment of a prognostic model

A prognostic model was developed using multivariate Cox regression analysis of the identified candidate genes. Risk scores were individually computed as follows:

.

where Coef_i refers to the coefficient for gene i, and Exp_i indicates the respective gene expression value. Individuals were subsequently divided into two groups: high-risk and low-risk, using the median risk score as the cutoff. Then, time-dependent receiver operating characteristic (ROC) curves were created. To assess the potential for overfitting, bootstrap internal validation with 1,000 resampling iterations was performed to calculate bias-corrected C-index and time-dependent AUCs with 95% confidence intervals. Calibration curves were generated to evaluate the agreement between predicted and observed survival probabilities at 2 years, 3 years, and 5 years. Furthermore, decision curve analysis (DCA) was performed using the ggDCA package in R to evaluate the clinical net benefit of the model at 2, 3 and 5 year time points, quantifying the potential value of the risk score in clinical decision-making across different threshold probabilities. Survival disparities between risk-stratified groups and across other clinical categories, were compared using Kaplan–Meier(KM) survival curves with log-rank testing. Furthermore, elucidating individual gene contributions to the model performance employed Shapley Additive exPlanations (SHAP) analysis for post-hoc explanatory interpretation.

Nomogram development and external validation

The relationships between the calculated risk scores and various clinical features (including gender, age, and TNM stage) were examined using Wilcoxon rank-sum or Kruskal-Wallis tests to evaluate the model's clinical applicability. To assess whether the risk score functioned as an independent prognostic factor, clinical variables, along with the risk score, were incorporated into multivariate Cox regression modeling. Then, a prognostic nomogram combining the independent clinical risk factors (e.g., Stage) and the genetic risk score constructed through the regplot R package to individualize survival probability predictions. Calibration curves were used to assess agreement of the nomogram-predicted survival probability with actual survival outcomes. Finally, the ultimate predictive capacity and generalizability of the integrated nomogram system were rigorously validated through time-dependent ROC curves and comprehensive KM clinical subgroup analyses across the cohorts.

Immune infiltration and immune subtype analyses

CIBERSORT, using the leukocyte gene (LM22) signature matrix, was used to estimate the relative proportions of 22 immune cell types to assess immune cell infiltration in LUAD patients. The relationships between prognostic gene expression levels and immunological infiltration were evaluated by Spearman correlation analysis. Immune, stromal, tumor purity, and ESTIMATE scores were derived via the ESTIMATE algorithm, and the Wilcoxon test assessed differences across risk groups. LUAD patients were assigned to six immune subtypes using the ImmuneSubtypeClassifier package¹⁶. Wilcoxon test was additionally used to compare immune subtype distributions across risk groups.

Immune checkpoint, immunophenoscore, and cancer immunity cycle analyses

In this study, the Wilcoxon rank-sum test served to assess 21 immune checkpoint genes¹⁷ across risk-stratified groups, aiming to characterize the LUAD immune landscape. Spearman correlation linked candidate prognostic genes to immune checkpoint genes. To assess differences in response to immune checkpoint inhibitors (ICIs) in LUAD patients at different risk levels, immunophenoscore (IPS) data for anti-PD-1 and anti-CTLA-4 treatments were obtained from The Cancer Immunome Atlas (TCIA)¹⁸, and the Tracking Tumor Immunophenotype (TIP)¹⁹ database was used to evaluate cancer-immunity cycle activity, comparing corresponding scores between risk groups.

Somatic mutation and drug sensitivity analysis

The TCGA mutations tool retrieved somatic mutation profiles for TCGA-LUAD cases to investigate variation in mutation patterns across risk groups. The maftools package processed and visualized mutation data. Tumor mutation burden (TMB) levels were determined for each specimen and contrasted between the two risk categories. Pharmacogenomic sensitivity analysis was conducted using the pRRophetic package according to the Genomics of Drug Sensitivity in Cancer(GDSC) database²⁰. Half-maximal inhibitory concentration (IC₅₀) values for anticancer drugs were predicted for each LUAD patient, and risk-group differences were quantified using the Wilcoxon rank-sum test.

Evaluation of expression levels for prognostic genes

Each dataset served to assess transcript levels of selected candidate outcome-related genes. To link gene expression with patient prognosis, optimal cutoffs were derived via the surv_cutpoint function in the survminer R package. Based on these thresholds, LUAD cases were categorized into high- and low-expression subsets for subsequent survival analysis.

Furthermore, five matched LUAD tumors and adjacent normal tissue pairs were obtained from Anhui Chest Hospital and qPCR validation was subsequently performed. Every participant provided written informed permission. Six candidate prognostic genes (PDGFB, LDHA, ZEB2, FKBP4, DMD, and S100B) were selected for qPCR validation. RNA was extracted from homogenized tissue samples using RNA extraction reagent, followed by chloroform extraction and isopropanol precipitation. The spectrophotometer enabled measurement of RNA concentration and purity. qPCR validation of the six candidate prognostic genes was performed by SYBR Green-based PCR master mix on a real-time PCR system: initial denaturation at 95 °C for 30 s, followed by 40 cycles of 95 °C for 20 s, 55 °C for 20 s, and 72 °C for 20 s. Relative expression was computed and standardized to Glyceraldehyde-3-Phosphate Dehydrogenase(GAPDH) with the 2^-ΔΔCt technique. Details of all reagents and instruments are provided in the Table of Materials.

Statistical analysis

Statistical analyses were conducted using statistical computing and graphing software. The protein-protein interaction network was visualized using network-analysis software. Following normality assessment, Student’s t-test was used for normally distributed continuous variables, and the Mann-Whitney U test for non-normally distributed ones.

Results

Candidate gene identification and functional analysis

A total of 4,915 DEGs were detected, including 2,485 down-DEGs and 2,430 up-DEGs (Figure 1A). Functional module clustering analysis classified these DEGs into six expression modules (C1-C6), among which the C1, C2, C3, and C5 modules showed relatively higher expression in tumor samples (Figure 1B). Functional annotation indicated that these modules were linked to b...

Discussion

LUAD is among the leading forms of lung cancer and continues to pose a considerable public health burden²¹. MC is a cell death-related process that may contribute to clearing cells with severe mitotic abnormalities and the suppression of tumor progression²². Therefore, exploring the potential role of MCRGs in LUAD may help identify prognostic biomarkers and improve understanding of LUAD biology. The study analyzed the prognostic relevance and potential functional significan...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
ClusterGVis	CRAN	version 0.1.2	Gene Expression Clustering and Visualization
Cytoscape	Cytoscape Consortium	version 3.9.1	Visualization of Protein Interaction Networks RRID:SCR_003032
CFX384 Real-Time Quantitative Fluorescent PCR System	Bio-Rad	CFX384 Touch	qPCR reagent
DMD-primer	Tsingke	N/A	Forward: 5’-GCTCAACCATCGATTTGCAGCC-3’ Reverse: 5’-TTCAGCCTCCAGTGGTTCAAGC-3’
FKBP4-primer	Tsingke	N/A	Forward: 5’-TGACTCCAGTCTGGATCGCAAG-3’ Reverse: 5’-CTGGTTTGCAGGTGATGTGGCA-3’
GEO	NCBI	Online Tool	Obtaining the LUAD External Validation Dataset RRID:SCR_005012
GAPDH-primer	Tsingke	N/A	Forward: 5’-GTCTCCTCTGACTTCAACAGCG-3’ Reverse: 5’-ACCACCCTGTTGCTGTAGCCAA-3’
glmnet	CRAN	version 4.1-8	Feature Selection in LASSO Regression RRID:SCR_015505
LDHA-primer	Tsingke	N/A	Forward: 5’-GGATCTCCAACATGGCAGCCTT-3’ Reverse: 5’-AGACGGCTTTCTCCCTCTTGCT-3’
limma	Bioconductor	version 3.60.6	DEG screening RRID:SCR_010943
maftools	Bioconductor	version 2.20.0	Visualization of Somatic Mutation Data RRID:SCR_024519
pRRophetic	CRAN	0.5	Prediction of Drug IC50 Values RRID:SCR_024417
PDGFB-primer	Tsingke	N/A	Forward: 5’-GAGATGCTGAGTGACCACTCGA-3’ Reverse: 5’-GTCATGTTCAGGTCCAACTCGG-3’
qPCR SYBR Green Master Mix	Vazyme	Q111-02	qPCR reagent
R software	R Foundation	version 4.2.2	Bioinformatics Analysis
regplot	CRAN	version 1.1	Creating Line Charts
survival	CRAN	version 3.8-3	Survival Analysis and Cox Regression RRID:SCR_021137
STRING	STRING Consortium	Online Tool	Protein Interaction Network Analysis RRID:SCR_005223
S100B-primer	Tsingke	N/A	Forward: 5’-GAAGAAATCCGAACTGAAGGAGC-3’ Reverse: 5’-TCCTGGAAGTCACATTCGCCGT-3’
survminer	CRAN	version 0.5.0	Survival Curve Plotting and Optimization RRID:SCR_021094
TRIzol	Tiangen	DP424	qPCR reagent
TCGA	NCI	Online Tool	LUAD Transcriptome and Clinical Data Acquisition RRID:SCR_003193
ZEB2-primer	Tsingke	N/A	Forward: 5’-AATGCACAGAGTGTGGCAAGGC-3’ Reverse: 5’-CTGCTGATGTGCGAACTGTAGG-3’