In this research, 20 PE placental tissues, among which 12 cases with severe features, were collected following cesarean sections, with diagnosis conforming to the guidelines established by the Task Force on Hypertension in Pregnancy. Correspondingly, control placental tissues (n=20), matched for age and body mass index (BMI), were also obtained. All placental tissues were sourced from pregnant women who delivered at the Affiliated Hospital of Qingdao University. Because gestational diabetes mellitus (GDM) is associated with an increased incidence of PE (30), exclusion criteria for the research were as follows: twin or multiple pregnancies; fetal structural abnormalities or chromosomal anomalies; the presence of comorbidities or complications including GDM, pre-pregnancy diabetes mellitus, chronic hypertension, cardiac, renal, or liver diseases, infectious diseases, or autoimmune disorders; history of blood transfusion, organ transplantation, or immunotherapy; and any history of smoking, alcohol consumption, or substance abuse. Basic clinical data were collected for this study, encompassing age, BMI, gestational week at delivery, parity, systolic and diastolic blood pressure, newborn weight, and one-minute Apgar score. A sample of maternal placental tissue, approximately 1 cm in diameter, was collected within ten minutes of placental delivery, specifically avoiding areas with infarcts or calcification. These samples were then placed in freezing tubes containing RNA preservation solution and kept at -80°C for subsequent analysis using RT-qPCR. All the participants of the study provided a written informed consent. The investigation was approved by the Ethics Committee of the Affiliated Hospital of Qingdao University (Approval No: QYFY WZLL 28705). It was carried out in rigorous adherence to the guidelines established by the committee.

The expression profile datasets GSE24129 (31), GSE30186 (32), GSE54618 (33) and GSE75010 (34) of patients with PE were obtained from GEO database (35) utilizing the GEOquery package (36). The GSE24129 dataset included 16 placental samples, evenly distributed between 8 PE cases and 8 Control group. The GSE30186 dataset comprised 12 placental samples, with an equal number of PE cases and Control group. The GSE54618 dataset consisted of 17 placental samples, including 5 from PE cases and 12 from Control group. Lastly, the GSE75010 dataset contained 80 PE cases and 77 Control group. The dataset GSE24129 and GSE75010 utilized the GPL6244 [HuGene-1_0-st-v1] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]. For the datasets GSE30186 and GSE54618, the associated platform was the GPL10558 Illumina HumanHT-12 V4.0 expression beadchip. The microarray GPL platform files facilitated related annotation for the probe names across these datasets. Specific information for each dataset is depicted in Supplementary Table 1 . Utilizing “mitochondrial energy metabolism” as the search keyword and focusing solely on protein-coding genes, we extracted 219 MMRGs from the datebase of GeneCards (https://www.genecards.org/),which offers extensive data on the human genes (37). Additionally, we derived 188 MMRGs from the published literature (38). By integrating these datasets and removing duplicates, we compiled a consolidated list of 384 MMRGs. The specific names of these genes are listed in Supplementary Table 2 .

We integrated the GSE24129, GSE30186 and GSE54618 datasets and then eliminated batch effects by applying the ComBat algorithm from the R package (39), followed by normalization using the normalize Between Arrays function. Thus, the Combined dataset (including 19 PE cases and 26 Control group) was obtained. Subsequently, we obtained DEGs by utilizing R’s limma package to carry out a differential analysis of the expression of all genes among the PE and the control cohort samples of the combined dataset. To make sure to capture all changes in expression levels, whether up-regulated or down-regulated, we made a screening standard of p < 0.05 plus | logFC | > 0 to further study the DEGs (40). The findings of variance analysis through the R package ggplot2 map volcano to display. Then, we took MMRGs and DEGs intersection to obtain the MMRDEGs.

The GO (41) approach is frequently employed in large-scale functional enrichment investigations for the purpose of categorizing genes into groups that are associated with biological process (BP), molecular function (MF), and cellular component (CC). The KEGG (42) serves as a crucial repository for genomic information, diseases, drug-related data and biological pathways. We conducted GO and KEGG annotation analyses of MMRDEGs by employing the R package clusterProfiler (43). We set a marked threshold (p < 0.05) for pathway selection, ensuring that only statistically significant pathways were considered in our analysis.

GSEA (44) is a widely utilized method for assessing variations in pathway activity and biological process involvement across different sample groups within an expression dataset. In this research, we initially carried out a differential gene expression analysis between various groups (PE/Control and High/Low Risk score) within the combined dataset. Subsequently, genes were categorized into two cohorts based on their logFC values: those with positive and those with negative logFC. For the enrichment analysis of these categorized genes, we utilized the clusterProfiler package. The GSEA configuration for this analysis utilized the following specifications: a seed of 2022, 1000 permutations, and a gene set size ranging from a minimum of 10 to a maximum of 500 genes. We retrieved the gene set “c2.Cp.All.V2022.1.Hs.Symbols.GMT [All Canonical Pathways]” containing 3050 entries from the Molecular Signatures Database (MSigDB) (45). Pathways which got a significant enrichment level (p < 0.05) were deemed markedly enriched.

GSVA (46) was designed to assess gene set enrichment within microarray and nuclear transcriptome data. This technique enables the conversion from diverse samples into a sample-specific gene expression matrix and can evaluate the pathway enrichment across multiple specimens. In this study, we also employed the gene set used earlier when GSEA analysis was performed. GSVA was carried out on gene expression matrices derived from distinct groups (PE/Control or High/Low Risk score) within the Combined dataset, utilizing this reference gene set. The analysis revealed functional disparities in enriched pathways between sample cohorts within the Combined dataset. Pathways with a significance level (p < 0.05) were further scrutinized; specifically, we selected and examined the 10 pathways exhibiting both the largest and smallest log fold change (logFC).

The RandomForest (RF) (47) technique is a collective learning approach that integrates numerous decision tree models. It belongs to the bagging (bootstrap aggregation) ensemble algorithm, which consists of multiple algorithms. RF is a commonly used approach for model building. By constructing multiple decision trees, the prediction results of each tree in the forest are aggregated using a voting method to obtain the final prediction result for a given sample. In this study, we utilized the MMRDEGs expression levels in the Combined dataset’s expression matrix to build a model using the RandomForest package with parameter set.seed (2023) and ntree = 1000.

We conducted a logistic regression analysis on MMRDEGs to construct a Logistic diagnostic model of the Combined dataset. Moreover, we employed the Logistic regression to analyze the association of the independent variables and dependent variables, when considering dependent variables as binary variables(PE cases and Control group). p < 0.05 was a significance level as criteria for identifying MMRDEGs and constructing the Logistic diagnostic model. The molecular expressions of MMRDEGs incorporated in this logistic regression model were visualized through Forest Plot.

Furthermore, we conducted the Least Absolute Shrinkage and Selection Operator (LASSO, the seed number is 2022) by R package glmnet (48) to process the MMRDEGs, which were screened out by utilizing our logistic regression model, to obtain the Logistic-LASSO regression model. LASSO regression analysis reduces overfitting incorporating a penalization factor (lambda × absolute value of slope), thereby improving its capacity for generalization while maintaining interpretability. The results obtained from LASSO analysis were depicted through variable trajectory plot techniques and diagnostic model plot.

Subsequently, we identified the common MMRDEGs by intersecting the MMRDEGs derived from both the RF model and the Logistic-LASSO regression model, which were then visualized using a Venn diagram. The expression levels of the common MMRDEGs in the Combined dataset were combined with the coefficients of these genes in the regression model of Logistic-LASSO to establish an MMRDEGs diagnostic model and to calculate corresponding Risk-scores. A Nomogram (49), a visual depiction of interrelations among several independent variables on a rectangular plane-coordinate system, was constructed based on the gene expression levels derived from the MMRDEGs diagnostic model generated through Logistic LASSO regression analysis in the Combined dataset. To examine the precision and distinguishing capability of our MMRDEGs diagnostic models, Decision Curve Analysis (DCA) (50), a straightforward approach for appraising molecular markers, diagnostic tests and clinical prediction models, was performed using the ggDCA R package.

The relative abundance of a variety of immune cell infiltrates within every sample was quantified utilizing the single-sample gene-set enrichment analysis (ssGSEA) algorithm. The method facilitated to label various immune cell types. For instance, regulatory T cells, CD8⁺ T cells, dendritic cells and macrophages. We represented the relative abundance of each immune cell type across the samples by enrichment scores, which were calculated utilizing ssGSEA methodology (51, 52). Using the ssGSEA algorithm from the GSVA R package (version 1.46.0), we calculated the enrichment scores of groups within high and low risk cohorts according to the MMRDEGs diagnostic model from the Combined dataset. These scores depicted the extent of immune cell infiltrations in individual specimen, thereby illustrating disparities of the abundance of immune cell infiltration among the different (High and Low) risk cohorts through box plots. Additionally, we examined the correlation of immune cell abundances among the high and low risk cohorts utilizing scatter plots. The association among immune cells and commonly altered MMRDEGs across these groups was analyzed using the Spearman statistical method and depicted in correlation dot plots, increasing our understandings of the immune landscape in relation to preeclampsia risk stratification.

CIBERSORT (53) is a kind of immune infiltration algorithm, that deconvoluted transcriptome expression matrices based on linear support vector regression, to assess the abundance and composition of immune cells within different samples. For this analysis, we input the expression matrix data of samples of the High and the Low risk groups defined by the MMRDEGs diagnostic model in the Combined dataset to CIBERSORT. Using the feature gene matrix of LM22, we refined the results by retaining solely those data points with immune cell enrichment scores >0, thus obtaining and visualizing the comprehensive findings of the immune cell infiltration abundance matrix. Those disparities in immune cell infiltration between the high risk and low risk cohorts were depicted utilizing stacked bar charts. We employed the Spearman statistical method to analyze the correlations among immune cells within the Combined dataset and utilized the R package ggplot2 to visualize the results. Moreover, the interactions among immune cells and commonly altered MMRDEGs were depicted using correlation dot plots, providing insights into the immune dynamics associated with different risk stratifications in PE.

The protein-protein interaction (PPI) network consists of individual proteins that engage with one another. In this study, we constructed the common MMRDEGs PPI network(minimum required interaction score: low confidence (0.150)) using the database of STRING (54). The network was visualized using Cytoscape, which allowed us to identify densely interconnected clusters within the PPI network. These clusters potentially signify molecular assemblies with unique biological roles, offering insights into the molecular mechanisms underlying PE.

ENCORI database (55) (https://starbase.sysu.edu.cn/) facilitates the exploration of interactions among various RNA types, including microRNAs-ncRNA, microRNAs-mRNA, ncRNA-RNA, and RNA-RNA, as well as interactions among RNA-binding proteins (RBPs) and ncRNAs or mRNAs. These interactions are curated utilizing degradome sequencing data and CLIP-seq, supporting comprehensive visual tools of investigating miRNA targets. In our study, we utilized the ENCORI database to forecast RBPs interacting with the commonly altered MMRDEGs. We established “pancancerNum> 27” as the threshold for selecting significant interactions. and the mRNA-RBP interaction network was rendered utilizing Cytoscape.

HTFtarget database (56) (http://bioinfo.life.hust.edu.cn/hTFtarget) integrates human transcription factors (TFs) and their corresponding control targets data. The CHIPBase database (https://rna.sysu.edu.cn/chipbase/) predicted transcriptional regulatory relationships among millions of TFs and genes. Utilizing both HTFtarget databases and CHIPBase, the TFs that link to common MMRDEGs were identified. We applied the screening criteria of having an upstream and downstream sample count greater than zero. Subsequently, the mRNA-TF interactive network was rendered visually utilizing Cytoscape software.

Placental tissues were lysed using FreeZol reagent (Vazyme, R711) following the manufacturer’s instructions. RNA concentration and purity were measured using a spectrophotometer. The isolated RNA was reverse transcribed into cDNA with a reverse transcription kit (Agbio, AG11705). Real-time polymerase chain reaction (qPCR) was then performed using the SYBR Green Pre-Mix Pro Taq HS qPCR Kit (Agbio, AG11701). Relative gene expression levels were normalized to β-actin and calculated using the 2^^(-ΔΔCt) method. The Supplementary Table 3 lists the primer sequences of mRNAs and internal control.

The entirety of data manipulation and statistical evaluation in this investigation was executed utilizing R software (Release 4.1.2). The independent Student’s t-test was used to compare continuous variables (fit normal distribution). We employed the Mann-Whitney U test (Wilcoxon rank sum test) for variables lacking normally distributed. Use Spearman correlation analysis to computer the findings unless otherwise specified. The p-values for statistical tests are two-tailed, and a threshold of 0.05 is deemed indicative of statistically meaningful results.

According to the technical roadmap of this experiment ( Figure 1 ), we first combined the three datasets (GSE24129, GSE30186, and GSE54618), then batched the data using the ComBat function from R’s sva package, and then utilized the Normalize Between Arrays function of the limma package to perform standardization procedures. The dataset of 19 PE cases and 26 Control group, that was combined, was obtained.

The before and after data processing boxplots and PCA plots of the combined dataset, according to the sample source, were showed in the Supplementary Figures 1A–D , respectively. The findings demonstrated that the expression profiles of samples from the Combined dataset exhibited a high degree of consistency, indicating successful removal of batch effects through data processing. The Combined dataset utilized for subsequent analyses represented the batch effect-corrected and normalized data.

We further utilized the limma package to standardize the GSE75010 dataset and compared the pre- and post-processing states of the dataset using boxplots ( Supplementary Figures 1E, F ). The boxplot analysis demonstrated that the expression levels of samples in the GSE75010 dataset became significantly more consistent after data processing.

The placenta serves as a critical organ facilitating material transport between mother and fetus, performing multiple functions during pregnancy such as immune protection, endocrine regulation, and serving as a conduit for nutrient and oxygen delivery. Its condition is closely associated with the health of both mother and child during gestation. Torbergsen T et al. first described a high incidence of preeclampsia in a family with mitochondrial disorder (57). Recent research into the mechanisms underlying preeclampsia has revealed mitochondrial dysfunction in both patients with preeclampsia and animal models (58).

We utilized limma package to explore the Combined dataset of PE cases and Control group. And we got 1073 differentially expressed genes using the threshold of | logFC | > 0 and p < 0.05, including 603 highly expressed genes in PE cases(the Control group of samples low expressed, logFC is positive, raised genes), and 470 genes low expressed in PE cases(the Control group of samples increased, logFC is negative, cut genes). And then, we presented the outcomes of differential expression analysis between the two groups in the Combined dataset using the volcano plot ( Figure 2A ). We then intersected these 1073 differently expressed genes (DEGs) with 384 mitochondrial energy metabolism-related genes (MMRGs). And then, 24 mitochondrial energy metabolism-related differentially expressed genes (MMRDEGs) were identified. The 24 MMRDEGs were ACSL3, ALDH16A1, ALDH1A3, ALDH4A1, ATG7, BTD, FBXL4, FOXO1, GAPDH, GLS, HK2, KCNJ2, LDHA, MFN2, NDUFS6, OCRL, PC, PGK1, PPARG, RARS2, SOD1, TPI1, VDAC1, XBP1 ( Figure 2B ).

We also generated a comparative map to analyze the differential expression of 24 MMRDEGs between the two cohorts ( Figure 2C ). The analysis suggested that 16 MMRDEGs exhibited significant differences, with ALDH16A1, ALDH4A1, ATG7, GLS and SOD1 significantly down-regulated while ALDH1A3, BTD, FOXO1, GAPDH, HK2, KCNJ2, LDHA, OCRL, PGK1, TPI1, and VDAC1 significantly up-regulated.

Then, we drew a simple numerical heat map derived from the expression matrix of these 16 MMRDEGs above ( Figure 2D ), and the visualization revealed substantial disparities in the expression patterns of the 16 MMRDEGs between the two sample groups. Additionally, we annotated the positions of these 16 MMRDEGs and draw a chromosome localization map ( Figure 2E ) by employing the RCircos package, from which the specific distribution of the 16 MMRDEGs on each chromosome can be obtained.

The biological processes (BP), molecular functions (MF), cellular components (CC), relationships between biological pathways, pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene function enrichment analysis based on Gene Oncology (GO) were carried out to analyze the 16 MMRDEGs. Pathways that below the P threshold of 0.05, were considered to be statistically significant. The outcomes showed that the 16 MMRDEGs main enriched in those BPs, such as pyruvate metabolic process, glycolytic process, ATP generation from ADP, generation of precursor metabolites and energy, ATP metabolic process. And in the CCs of the mitochondrial matrix. It was enriched in acting on the aldehyde or oxo group of donors, oxidoreductase activity, NAD or NADP as acceptor, aldehyde dehydrogenase (NAD+) activity, aldehyde dehydrogenase [NAD(P)+] activity, protein phosphatase binding, oxidoreductase activity, acting on the CH-NH group of donors, NAD or NADP as acceptor and other MFs ( Figure 3A ). It was also enriched in Glycolysis/Gluconeogenesis, HIF-1 signaling pathway, Carbon metabolism, Alanine, aspartate, glutamate metabolism, Inositol phosphate metabolism ( Figure 3B ) and other KEGG pathways ( Supplementary Table 4 ). In addition, the enrichment consequences of the BP pathways ( Figure 3C ), CC pathways ( Figure 3D ), MF pathways ( Figure 3E ), and KEGG pathways ( Figure 3F ) of GO analysis were presented utilizing ring network diagrams.

We utilized the Pathview R package for pathway mapping to illustrate the KEGG enrichment results of Glycolysis/Gluconeogenesis ( Supplementary Figure 2A ), Carbon metabolism ( Supplementary Figure 2B ), Alanine, aspartate and glutamate metabolism ( Supplementary Figure 2C ), Inositol phosphate metabolism ( Supplementary Figure 2D ), and HIF-1 signaling pathway ( Supplementary Figure 2E ).

To appraise the influence of gene expression levels of genes from PE and Control groups of Combined dataset on PE, we examined the relationships between the expression levels of all genes in different groups (PE/Control) of the Combined dataset and the BPs, CCs, and MFs they played, by employing the Gene Set Enrichment Analysis (GSEA). p < 0.05 was set as the significant enrichment criterion. The results demonstrated a significant enrichment of genes from different (PE/Control) groups in the Combined dataset, specifically in the vascular smooth muscle contraction pathway ( Figure 4B ), IL9 signaling pathway ( Figure 4C ), Notch signaling pathway ( Figure 4D ), IL2 signaling pathway ( Figure 4E ), IL6/7 signaling pathway ( Figure 4F ), cell surface interactions at the vascular wall ( Figure 4G ), and other pathways ( Supplementary Table 5 ). In addition, the outcomes of GSEA analyzing genes between distinct cohorts (PE/Control) of the Combined dataset were depicted by mountain plot ( Figure 4A ).

To investigate the distinctions between disease and controls from the Combined dataset, we then performed Gene Set Variation Analysis (GSVA). From the pathways with p < 0.05, we identified 10 pathways with the highest and lowest logFC for further examination ( Supplementary Table 6 ), respectively.

The results of GSVA analysis on all the genes of the Combined dataset revealed significant differences among PE and Control groups. Specifically, IKEDA Mir133 targets DN, hyaluronan biosynthesis and export, RHOT1 GTPASE cycle, neurofascin interactions, Irinotecan pathway, Aflatoxin B1 metabolism, Sulindac metabolic pathway, weber methylated LCN in SPERM DN, Tomlins metastasis upregulation of steroid biosynthesis. Additionally, activated NTRK2 signals through FYN and PI3K pathways were observed along with NTRK2 activation of RAC1. Furthermore, HIF1A and PPARG were found to regulate glycolysis. Calvet Rinotecan sensitive vs resistant upregulation was also identified as well as Korkola choriocarcinoma involvement. Lastly erythrocytes demonstrated oxygen uptake and carbon dioxide release while Tesar Alk targeted human es 4D and 5D DN along with JAK targeting mouse es D4 DN. Utilizing the outcomes derived from GSVA, we carried out a differential expression analysis of 20 pathways among PE and control cohorts of the Combined dataset. Subsequently, we created a heatmap illustrating the particular differential analysis outcomes ( Supplementary Figure 3A ) employing the R package. Furthermore, we assessed the extent of group divergence for these 20 pathways across various cohorts from the Combined dataset, utilizing the Mann-Whitney U test, and we use a group comparison plot to illustrate the outcomes ( Supplementary Figure 3B ). The findings demonstrated marked differences in pathway expression among disease control cohorts within the Combined dataset.

Based on the Combined dataset, we examined the expression levels of the 16 MMRDEGs using Random Forest algorithm (RF) to evaluate the values in diagnosis of the 16 MMRDEGs ( Figure 5A ). IncNodePurity (Increase in NodePurity) indicates the enhancement in node purity. The higher the node purity, the less impurities it contains (that is, the smaller the Gini coefficient). We applied an IncNodePurity threshold of > 0.5 to filter the specific analysis outcomes. The findings ( Figure 5B ) revealed that 7 diagnostic markers were obtained by RF algorithm. They are: OCRL, GAPDH, TPI1, LDHA, SOD1, HK2 and PGK1.

Logistic regression was performed utilizing the expression levels of 16 MMRDEGs in the Combined dataset, with the screening criterion of p< 0.05 ( Figure 5C ). The Logistic regression model included a total of 8 MMRDEGs (ALDH16A1, ATG7, BTD, GAPDH, HK2, LDHA, OCRL and TPI1), and the diagnostic model was developed by the expression relative quantities of the 8 genes in the combined dataset (the expression levels were evaluated by Least Absolute Shrinkage and Selection Operator (LASSO) analysis). And the findings of the LASSO analysis were illustrated via the LASSO regression model diagram ( Figure 5D ) and the LASSO variable trajectory plot ( Figure 5E ). The findings indicated that the diagnostic model comprised 6 MMRDEGs, which were: ALDH16A1, ATG7, GAPDH, LDHA, OCRL and TPI1.

Then we interposed the MMRDEGs from the RF model and the MMRDEGs from the Logistic-LASSO regression model ( Figure 5F ), and 4 Common MMRDEGs (p < 0.05) were obtained, which were OCRL, GAPDH, TPI1 and LDHA.

Next, we examined the differential expression of the 4 Common MMRDEGs in the placental tissues of preeclamptic and normal mothers using RT-qPCR. The demographic characteristics of the PE patients are presented in Supplementary Table 7 . The findings indicated that the mRNA expressions of the 4 common genes were notably elevated in the placental tissues of the PE cases relative to the Control group (p < 0.05, Figure 5G ). These four Common MMRDEGs (OCRL, GAPDH, TPI1, and LDHA) were identified for the first time in a study of PE. This novel discovery offers fresh insights into the role of mitochondrial metabolism in preeclampsia and may establish a foundation for the development of future biomarkers and therapeutic targets.

And then, utilizing the expression level of the four Common MMRDEGs in Combined dataset and corresponding coefficients established by applying LASSO analysis, we obtained the MMRDEGs diagnostic model of 4 Common MMRDEGs.

The diagnostic model for MMRDEGs included four Common MMRDEGs. We used combined logistic regression analysis to process the dataset’s expression levels to construct a logistic regression model for MMRDEGs. Additionally, we generated a nomogram depicting the impact of these four common MMRDEGs on the logistic regression model ( Figure 6A ). Our findings revealed that among all variables, OCRL exhibited notably superior effectiveness within the MMRDEGs logistic regression model.

The diagnostic model’s clinical value was appraised through decision curve analysis (DCA), and the findings were presented in Figure 6B . In the DCA graph, a model’s line consistently surpasses those of “All negative” and “All positive” within a specific range, greater net benefits can be obtained, indicating a stronger model performance. Our findings demonstrate that our constructed model exhibits considerable accuracy in diagnosing PE.

To further substantiate the value of the MMRDEGs diagnostic model, we drew ROC curves utilizing the Risk Scores of the diagnostic model of MMRDEGs and the information for grouping (PE/Control) of the Combined dataset and displayed the outcomes ( Figure 6C ). The MMRDEGs diagnostic model exhibited substantial precision in the diagnosis of the two groups (PE/Control) (AUC = 0.970, CI=0.930-1.000, Figure 6C ).

We further validated the diagnostic performance of the MMRDEGs diagnostic model using the GSE75010 dataset. Specifically, we calculated the risk scores by applying the formula derived from the MMRDEGs diagnostic model and the gene expression profiles in GSE75010. Subsequently, we incorporated the grouping information to construct the ROC curve. The results indicated that the MMRDEGs diagnostic model exhibited satisfactory accuracy in distinguishing the PE and Control groups within the GSE75010 dataset (AUC = 0.877, CI 0.823-0.932, Figure 6D ).

We also performed functional similarity analysis for four Common MMRDEGs and displayed them using a boxplot. We calculated the semantic similarity of sets of GO terms, GO terms, gene products and gene clusters through the R package GOSemSim. Similarity analysis was performed only on genes that were annotated to pathways in MF, BP, and CC. Finally, functional similarity analysis results between four Common MMRDEGs were obtained and visualized by Boxplot ( Figure 6E ). The findings indicated that LDHA exhibited the greatest functional similarity score in comparison to other Common MMRDEGs (the X-axis of D graph is the similarity score, with higher values indicating increased functional similarity to other genes).

Initially, we categorized those disease samples from the Combined dataset into the Low-Risk Score group and the High-Risk Score group utilizing the median Risk-Score of the previous MMRDEGs diagnostic model and performed a differential analysis between the two groups utilizing the limma package ( Figure 7A ). Based on the results of the differential analysis, we conducted GSEA to explore the relationship among the MFs, the CCs, the BPs and the expression of all genes involved between the different groups (Low/High Risk-Score group) in the Combined dataset, using the threshold of p < 0.05 for enrichment selection. The findings demonstrated a significant enrichment of genes linked to the citric acid TCA cycle and respiratory electron transport ( Figure 7C ), IL5 signaling pathway ( Figure 7D ), IL7 signaling pathway ( Figure 7E ), IL6 signaling pathway ( Figure 7F ), energy metabolism ( Figure 7G ), electron transport chain Oxphos system in mitochondria ( Figure 7H ) as well as other pathways, indicating their association with High and Low Risk cohorts ( Supplementary Table 8 ). Furthermore, the GSEA outcomes of genes among the High Risk-Score cohort and the Low Risk-Score cohort in the Combined dataset were presented by mountain plot ( Figure 7B ).

To investigate the disparities among the High Risk-Score cohort and the Low Risk-Score cohort in the Combined dataset, we subsequently conducted GSVA. From the pathways with a p < 0.05, we identified 10 pathways with the highest and lowest logFC for further analysis (refer to Supplementary Table 9 for detailed information), respectively. The GSVA results of all genes revealed significant differences in 20 pathways among the High and Low Risk-Score cohorts in the Combined dataset. These pathways include defective CSF2RB causes SMDP5, IPS LCP with H3K4ME3 and H3K27ME3, Korkola choriocarcinoma DN, FGFR3B ligand binding and activation, Aml methylation Cluster 7 DN, Turashvili breast carcinoma Ductal vs Lobular DN, FTO obesity variant mechanism, miscellaneous substrates, PEPI pathway, ES LCP with H3K4ME3 and H3K27ME3 angiogenic targets of VHL HIF2A up regulation Biocarta Myosin pathway OPN targets Cluster 3 Myc targets DN CTNNB1 pathway and proliferation mesothelioma survival up schavolt targets of TP53 and TP63 MAPK11 targets Pujana breast cancer with BRCA1 mutated DN regulation of PTEN localization. According to the GSVA outcomes, we analyzed the differential expression of 20 pathways among the Low-Risk cohort and the High-Risk cohort in the Combined dataset, and the specific differential analysis findings ( Supplementary Figure 4A ) was showed as a heatmap by the R package. Furthermore, we employed the Mann-Whitney U test to examine the group distinction level of 20 pathways between diverse cohorts in the Combined dataset and displayed the findings by group comparison plot ( Supplementary Figure 4B ). The findings suggested that all the expressions of the 20 pathways were markedly different among the Low-Risk cohort and High-Risk cohort in the Combined dataset (p<0.05).

We categorized PE samples in the Combined dataset into the Low Risk-Score and the High Risk-Score cohorts by the median Risk-Score of the MMRDEGs diagnostic model.

To study the difference of immune infiltration between the Low/High Risk-Score groups of the Combined dataset, we applied ssGSEA algorithm to computer the abundance of 28 immune cell infiltration in the two risk-score sample groups. And then, we used Mann-Whitney U test to analyze the differences of the two abundances of the Low and High Risk-Score groups, using group comparison plot to exhibit the results ( Figure 8A ). The findings suggested that there were two immune cells, namely Neutrophil and Plasmacytoid Dendritic cell, showing statistically differences in the abundance between the Low and High Risk-Score groups (p < 0.05).

We plotted the correlation scatter plots showing the relationship among Neutrophil and Plasmacytoid Dendritic cells in the Low-Risk cohort ( Figure 8B ) and the High-Risk cohort ( Figure 8C ) from the Combined dataset. The outcomes showed that, in the Low Risk-Score group, there was a marked inverse association among Neutrophil and Plasmacytoid Dendritic cells ( Figure 8B , R = -0.709, p = 0.028). However, there was no association between the two immune cells in the High Risk-Score group ( Figure 8C ).

We used Spearman’s statistical algorithm to calculate the association between the infiltrating abundances of the Neutrophil, Plasmacytoid Dendritic cells in the Low and High Risk-Score cohorts, and the expression of the four Common MMRDEGs in the Combined dataset data group ( Figures 8D, E ). The findings suggested that Neutrophil was positively correlated with the four Common MMRDEGs in the Low Risk-Score cohort of the Combined dataset ( Figure 8D ), moreover Neutrophil and OCRL had the strongest association. In the High Risk-Score cohort of the Combined dataset data, Plasmacytoid Dendritic cells had the strongest correlation with OCRL ( Figure 8E ).

The CIBERSORT method was utilized to estimate the abundance of 22 immune cell infiltrations in both the Low and High Score cohorts. A stacked bar chart was employed to graphically depict the distribution of immune cells across the dataset samples ( Supplementary Figure 5A ). There were 22 immune cells with non-zero infiltration abundances within the Combined dataset according to the results.

We used Spearman’s statistical algorithm to assess the relationships among the 22 immune cells ( Supplementary Figure 5B ), and the findings suggested that the number of positive and negative associations between the 22 immune cells was basically equal, among which Mast cells activated and B cells memory had the strongest correlation.

We subsequently computed the association among immune cells and the four Common MMRDEGs using Spearman’s statistical algorithm ( Supplementary Figure 5C ). The results showed that T cells CD4 memory activated, Dendritic cells resting, and T cells gamma delta were moderately positively correlated with the four Common MMRDEGs in the Combined dataset. Among all the associations examined, the most pronounced relationship was detected among naïve B cells and GAPDH.

Since these four Common MMRDEGs (OCRL, GAPDH, TPI1, LDHA) are the most potentially valuable biomarker genes identified during model construction, they are suspected to play a crucial role in related disease processes. Therefore, conducting an in-depth study on their interactions can aid in understanding their biological mechanisms and clinical applications. Therefore, we utilized the STRING database (PPI network, minimum required interaction score: low confidence (0.150)) to perform the PPI analysis of the 4 Common MMRDEGs (treated as hub genes) and visualized by Cytoscape software ( Figure 9A ).

And then, we utilized the ENCORI database to forecast RNA binding proteins (RBPs) that interacted with four Common MMRDEGs and subsequently visualized the mRNA-RBP interaction network by Cytoscape software ( Figure 9B ). The mRNA-RBP interaction network, in which green quadrilateral blocks presenting RBPs and the blue quadrilateral blocks presenting mRNAs, was composed of 4 Common MMRDEGs (OCRL, GAPDH, TPI1 and LDHA) and 51 RBP molecules, which constituted 58 pairs of mRNA-RBP interaction relationships. The specific mRNA-RBP interaction relationships are depicted in Supplementary Table 10 .

We utilized the CHIPBase database (version 3.0) and hTFtarget database to identify transcription factors (TFs) that bound to the four Common MMRDEGs. Then we screened by “Number of samples found (downstream)>0” and “Number of samples found (upstream) >0”, and finally got 3 Common MMRDEGs (OCRL, GAPDH, TPI1) and 39 pairs of interaction data of 29 TFs were graphically represented utilizing Cytoscape software ( Figure 9C ). In the mRNA-TF interaction network, those blue quadrilateral blocks represent mRNAs, and the green quadrilateral blocks are TFs. The detailed mRNA-TF interactions are depicted in the Supplementary Table 11 .

We employed the Comparative Toxicogenomics Database (CTD) to identify small molecule compounds or potential drugs that interact with four commonly observed MMRDEGs. The selection criterion for mRNA-Drugs interaction pairs was set as “Reference Count” > 1. To render the mRNA-Drug interaction network ( Figure 9D ), we employed Cytoscape software. Within the mRNA-Drugs interaction network, the blue quadrilateral blocks signify mRNAs, while the green quadrilateral blocks denote drugs. Our analysis revealed that our mRNA-Drugs interaction network consisted of three common MMRDEGs (OCRL, LDHA, and TPI1) and twenty-four drug molecules, forming thirty mRNA-Drugs interaction associations. Detailed information regarding these specific interactions can be found in Supplementary Table 12 .

The AlphaFold Protein Structure Database (https://www.alphafold.ebi.ac.uk/) encompasses approximately 350,000 protein structure predictions generated by the AlphaFold AI system. This comprehensive database includes predictions for humans and 20 widely studied model organisms in biological research, such as E. coli, Drosophila, zebrafish, and mice. Remarkably, AlphaFold has successfully predicted the structures of 98.5% of human proteins within the human proteome. By combining AlphaFold’s structural prediction, we can more comprehensively construct and understand complex interaction networks, revealing the important roles of these genes in cellular metabolic regulation. This, in turn, provides a molecular basis for exploring the mechanisms of related diseases. To investigate the protein structures of four common MMRDEGs, we leveraged the resources provided by the AlphaFold website and presented our findings in Figures 10A–D .