Expression profiles, RNA sequencing data, mutation information and clinical features of samples were downloaded from the open TCGA database (https://portal.gdc.cancer.gov/). Samples with incomplete data were excluded. Furthermore, samples with survival times <20 days or >8 years were removed from the present study since their survival results were not relevant to EC. For survival time <20 days, the death was considered to be associated with postoperative complications; for survival time >8 years, the present study considered the cause of death to be unrelated to EC as patients with EC with a survival of ≥5 years are defined clinically cured. Thus, 518 EC samples and 23 normal samples were enrolled for further study. Copy number variation (CNV) data was acquired from the UCSC Xena database (http://xena.ucsc.edu/). The MSI data of the samples was downloaded from the Cancer Imaging Archive database (https://tcia.at/home).

FAM-associated genes were obtained from two online gene sets and one published gene set. The HALLMARK_FATTY_ACID_METABOLISM gene set, consisting of 158 genes, was downloaded from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Using ‘fatty acid metabolism’ as the keyword, a total of 543 genes with relevance scores >20 were downloaded from GeneCards (https://www.genecards.org/). A total of 92 FAM-associated genes were obtained from a previous study (29), which introduced metabolism-associated signatures based on fatty acid degradation, elongation and biosynthesis. The feasibility of this genetic source has been established (30). A total of 669 genes were obtained after excluding duplicate genes (Table SI).

The ssGSEA algorithm was used to evaluate the immune cell content and immune function of each sample. Based on the evaluation results, the samples were classified into high- and low-immune groups using the R package ‘limma’ (version 3.26.9) with a false discovery rate (FDR) >0.05 and fold change >2. The CIBERSORT algorithm was used to quantify 22 types of tumor-infiltrating immune cells in high- and low-immune samples using the R package ‘CIBERSORT’ (version 1.04) as previously described (31). Microenvironment Cell Populations-Counter (MCPcounter) analysis is a robust quantification method (32), which is based on the R packages ‘MCPcounter’ and ‘limma’. CIBERSORT, MCPcounter and ssGSEA were used to compare the differences in immune cell types and immune-associated functions between the EC and normal samples. The intersection of the three algorithms was considered more reliable, strengthening the confidence in the evidence. The R packages ‘limma’ (version 3.26.9) and ‘ggpubr’ (version 0.6.0) were used to visualize the differential results. The immune cell estimation of each EC sample was also downloaded from Tumor Immune Estimation Resource (TIMER; version 2.0, http://cistrome.shinyapps.io/timer/).

The common genes between the union of FAM-associated genes and differential genes of the two immune groups were defined as differentially expressed genes (DEGs). A total of 50 genes were identified. The association between prognosis and these 50 genes was analyzed by univariate Cox regression analysis. A total of 22 genes with prognostic values were obtained. Least absolute shrinkage and selection operator (LASSO) Cox regression analysis of prognosis-associated genes was conducted using the R package ‘glmnet’ (version 4.1–8). The risk score for each sample is the sum of the multiplication of the expression values of genes and correlation coefficients, which can be further explained as follows:

Where Expk indicates the expression level of genes and Coefk denotes the risk coefficient of genes. The R package ‘survival’ (version 3.7.0) was used for Kaplan-Meier survival analysis, and univariate and multivariate Cox regression analyses were conducted. Additionally, to determine the accuracy and reliability of the risk assessment model, the R package ‘timeROC’ (version 0.4) was used to plot the receiver operating characteristic (ROC) curve and calculate the area under the curve (AUC). The R packages ‘timeROC’ (version 0.4) and ‘rms’ (version 6.9–0) were used to develop a predictive model that merged gene expression values and clinical features. The R package ‘ggDCA’ (version 1.2) was used for the decision curve analysis to evaluate the predictive ability of the nomogram.

The R package ‘scatterplot3d’ (version 0.3–44) was used to conduct principal component analysis (PCA) to improve data visibility. ‘Rtsne’ R package (version 0.17) was employed to conducted t-distributed stochastic neighbor embedding dimension reduction (t-SNE).

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses were conducted using the R package ‘clusterProfiler’ (version 4.14.4) for the differential genes between the high- and low-immune groups and differential genes between high- and low-risk groups. The R package ‘GSVA’ (version 2.0.2) was used for enrichment analysis of biological functions. The enrichment analysis results were visualized using the R packages ‘enrichplot’ (version 1.13.1.994) and ‘ggplot2’ (version 0.6.0). Statistical significance was defined as an FDR <0.05.

Human EC tumor samples and normal endometrium tissues were obtained from female patients aged 25–70 who were treated surgically from 2022 to 2023 at the Department of Obstetrics and Gynecology, Qilu Hospital of Shandong University (Jinan, China). The patients diagnosed with EC were the experimental group and patients who did not have EC or a endometrial precancerous lesion but received hysterectomy for other reasons were set as the normal controls. All methods were conducted in accordance with the relevant guidelines and regulations and the experiment was approved by the Ethics Committee of Qilu Hospital of Shandong University (Jinan, China). Clinical specimens were fixed in 4% paraformaldehyde for 12 h at room temperature, then embedded in paraffin. The specimens were cut into 4-µm paraffin sections, dewaxed with xylene, hydrated with a descending alcohol series, and endogenous peroxidase activity was blocked according to the instruction (ZSGB-BIO; cat. no. PV-9000). Subsequently, antigen retrieval was carried at 95°C using pH 6.0 sodium citrate and non-specific binding sites were blocked with goat serum (Beyotime Institute of Biotechnology; cat. no. C0265) for 30 min at room temperature. The samples were then washed with PBS. Rabbit anti-PPARG coactivator 1α antibody (cat. no. A20995; ABclonal Biotech Co., Ltd.) was administered to the slides at a 1:2,000 dilution and slides were incubated at 4°C overnight. The samples were then incubated with horseradish peroxidase-labeled secondary antibody, according to the instructions of the kit (Mouse/Rabbit Polymer Detection System; ZSGB-BIO; cat. no. PV-9000), at room temperature for 30 min, and treated with diaminobenzidine for 1 min, followed by hematoxylin staining for 5 min at 37°C. Finally, images were captured with an Olympus light microscope (IX71; Olympus Corporation). The IHC images were analyzed using Image-Pro Plus 6.0 (Media Cybernetics, Inc.) to determine the average optical density (AOD), which can be expressed as AOD=integrated optical density/area.

Bioinformatic analyses were conducted using R (https://www.r-project.org/; version 4.4.2) and the RStudio (https://posit.co/products/open-source/rstudio/; version 2023.06.2+561.pro5) software. χ² test was used to analyze differences in clinical information. Kaplan-Meier analysis and the weight of ‘aalen’ calculated by the R package ‘timeROC’ (version 0.4) were used to conduct survival analyses and draw the ROC curves. Log-rank test was used to calculate the P-value in the survival analysis. The ‘cph’ function from the ‘rms’ package (version 6.9–0) was used for survival calibration in risk model. Analysis of paired data was performed using the paired t-test. Analysis of group data was performed using the Wilcoxon rank sum test and Kruskal-Wallis test (non-parametric) using the ‘stat_compare_means’ algorithm in the ‘ggpubr’ R package (version 0.6.0). An unpaired Student's t-test was conducted for the statistical analysis of the IHC staining. Linear regression and Spearman's rank correlation analysis was performed. All statistical P-values were two-sided, and P<0.05 was considered to indicate a statistically significant difference. The waterfall function of the R package ‘maftools’ (version 3.20) was used to compare the TMB in high- and low-risk groups. Heat maps were created using the R package ‘ComplexHeatmap’ (version 2.22.0).

The immune cell content in each EC sample was calculated using the ssGSEA algorithm. Using this evaluation, samples were clustered into high- and low-immune groups (Fig. 1A). t-SNE further indicated that the clustering was effective and the difference between groups was notable (Fig. 1B). The TME score of each sample and the content of immune cells computed by ssGSEA in the high- and low-immune groups are displayed in Fig. 1C. The analysis indicated that low-immune groups contained lower immune cell infiltration and response, resulting in lower immune scores in the TME. There was also a lower stromal score in the low-immune group, which led to a lower ESTIMATE score since this is the sum of the immune and stromal scores. Therefore, the low-immune group had a higher tumor purity. The opposite was true in the high-immune group. For the estimated TME score, quantitative analysis demonstrated the difference in TME between the high- and low-immune groups in immune, stromal and ESTIMATE scores (P<0.001; Fig. 1D).

Enrichment analysis between groups revealed that the enriched GO molecular function terms were associated with receptor binding. The enriched GO cellular component terms were associated with the plasma membrane, and the enriched biological process terms were associated with ‘positive regulation of cell activation’ and immunity (Fig. 1E). KEGG enrichment analysis of DEGs was also conducted, and ‘cytokine-cytokine receptor interaction’ and ‘cell adhesion molecules’ were detected (Fig. 1F).

Between the high- and low-immune groups, 2,609 DEGs were identified (Fig. 2A). A total of 669 FAM-associated genes were obtained after deleting duplicate genes from the two online gene sets and one published gene set. Differentially expressed FAM-associated genes were obtained by intersecting the total FAM-associated genes and DEGs from the aforementioned immune cluster, leaving 50 genes (Fig. 2B). Based on Cox regression analysis, 22 prognostic genes were selected (Fig. 2C). Tumor samples from TCGA were randomly divided into training and test cohorts, with an equal number of samples in both groups. Subsequently, clinical manifestation analysis was conducted for patients with clinical manifestations in the two cohorts. Any patients in the unknown group were not analyzed. The analysis revealed that there was no difference in clinical manifestations between the two cohorts (Table I). After LASSO analysis, three genes with the best λ value were chosen to construct the optimum signature (Fig. 2D and E). The risk score of each sample was calculated according to the following formula: Risk score=(expression level of PPARGC1A × 0.653085171388467) + (expression level of RBP4 × 0.323298560719871) + (expression level of CRABP2 × 0.647075550554985).

A series of analyses on the constructor genes, CRABP2, RBP4 and PPARGC1A, were conducted. According to TCGA, CRABP2 had lower expression levels and RBP4 had higher expression levels in normal samples, but the P-value did not reach significance (Fig. 2F and G). PPARGC1A had significantly higher expression levels in normal tissue compared with tumor tissue (Fig. 2H). The mRNA expression levels of PPARGC1A were found to be significantly reduced in both carcinoma and adjacent tissue within the same patient (Fig. 2I). The difference in expression levels of PPARGC1A among different grades of EC is displayed in Fig. 2J, and the difference in expression levels of PPARGC1A among different stages of EC is shown in Fig. 2K. Clinical samples were used to confirm the difference in expression levels of PPARGC1A between EC and normal tissue. IHC analysis of 103 EC and 17 endometrium samples revealed higher PPARGC1A levels in normal samples (unpaired Student's t-test; Fig. 2L and M). Together, these results suggested that PPARGC1A may be an indicator of EC.

Based on the median risk score of 0.93, samples were divided into high-risk and low-risk groups in the training cohort. A significant difference was shown between risk groups in terms of overall survival (OS) and progression-free survival; the low-risk group had an increased survival rate (Fig. 3A and B). The risk curve suggested that as the risk increased, the number of mortalities increased (Fig. 3C). A significant but weak correlation was found between the OS and the risk score based on correlation analysis (Fig. 3D). The time-dependent ROC curve of the prognostic signature in the training cohort indicated a reliable prediction of prognosis in terms of 1-, 3- and 5-year survival (Fig. 3E). An evaluation of the calibration curve also confirmed the effectiveness of the signature (Fig. 3F). The test cohort and the whole tumor samples were split into high- and low-risk groups based on the same cut-off value as the training cohort. The same analysis was applied in the test cohort (Fig. S1A-F) and all samples (Fig. S2A-F), and the results were consistent with the aforementioned results.

Disruption of FAM can also lead to changes in body weight (6). Therefore, the BMI of each sample in the high- and low-risk groups was investigated. There was no significant difference between the high- and low-risk groups of the training cohort (data not shown); however, a significant difference was observed in the test cohort (Fig. S1G) and all samples (Fig. S2G), indicating that patients in the low-risk group had a higher BMI score. The relationship between abnormal FAM and BMI and EC needs to be further verified.

PCA was conducted on all EC samples. Compared with other classification methods (Fig. 3G-I), the gene signature established displayed differences between groups as the distribution of dots were separated (Fig. 3J), confirming that grouping and demarcation points of the risk score had been selected effectively. Based on univariate logistic analysis, the signature was associated with a poor prognosis (Fig. 3K). According to multivariate logistic analysis, grade and International Federation of Gynecology and Obstetrics (FIGO) stage were independent prognostic factors, while the signature did not have statistical validity (Fig. 3L). Thus, the evidence was insufficient to indicate that the signature was an independent prognostic factor.

Gene set variation analysis (GSVA) of GO and KEGG terms was conducted in high- and low-risk groups. GSVA of GO terms between groups was mainly associated with dendritic regulation and metabolic processes (Fig. S3A; Table SII). The top three pathways of GSVA enrichment on KEGG between groups were ‘fatty acid metabolism’, ‘alpha linolenic acid metabolism’ and ‘maturity onset diabetes of the young’ (Fig. S3B; Table SIII). All enrichment analyses indicated that metabolism and diabetes mellitus were closely associated with EC groups, suggesting that these factors may have an impact on EC progression.

An immune-associated analysis was conducted to determine whether high- and low-risk groups had different immune compositions. The MCPcounter method was used to assess the immune cell contents of each group. The levels of three types of immune cells were found to be statistically different between the two groups. Levels of B lineage cells and cytotoxic lymphocytes were higher in the low-risk group, while the level of T cells was higher in the high-risk group (Fig. S3C-E). The immune cell estimation of each EC sample was downloaded from TIMER version 2.0 (https://cistrome.shinyapps.io/timer/) (33). The results indicated that B cells, CD4 T cells, CD8 T cells, dendritic cells and macrophages were weakly negatively associated with the risk score, and a significant difference was observed (Fig. S3F-J). CIBERSORT was then used to analyze the differences in immune cell infiltration between high- and low-risk groups. The levels of six immune cell types were different between the groups. The levels of five cell types, including resting dendritic cells, were higher in the low-risk group and the level of activated dendritic cells was higher in the high-risk group (Fig. S3K). The immune cell infiltration of each sample in the high- and low-risk groups is displayed in Fig. S3L. Furthermore, ssGSEA was conducted to compare immune cell types between high- and low-risk groups (Fig. 4A). A total of 16 immune cell types were found to be significantly different between the low-risk and high-risk groups, with increased clustering in the low-risk group. The results of immune infiltration analysis assessed by four independent algorithms are summarized in Table SIV. The majority of the correlations were negative, which indicated that immune cells were gathered in low-risk groups. T cells counted using MCPcounter and dendritic cells assessed by CIBERSORT exhibited positive correlations. The positive correlation of activated dendritic cells was consistent with the negative correlation of resting dendritic cells. However, the positive correlation of T cells calculated by MCPcounter was contradictory to the finding that all types of T cells exhibited a negative correlation according to the four algorithms and this warrants further investigation.

In addition to immune cells, the association between immune-associated functions and the risk score was also examined. Based on the ssGSEA algorithm, each sample was scored for immune-associated functions, and the scores of samples were displayed based on the high- and low-risk groups (Fig. 4B). Distinct differences in the score were observed between the two groups. The qualitative differences in scores between groups are shown in Fig. 4C. Based on the results, it was concluded that several immune-associated functions were more concentrated in the low-risk group, with the exception of antibody-drug conjugates and type I IFN response. Further investigation of the expression level difference of human leukocyte antigen subtypes between high- and low-risk groups was carried out (Fig. 4D). Of the 24 subtypes studied, 13 were expressed differently between groups and all 13 subtypes were highly expressed in the high-risk group.

The aforementioned results indicated differences in immune infiltration between the high- and low-risk groups. The low-risk group had higher immune content, which indicated that immunotherapy may be beneficial in these patients with EC. The findings also suggested that the created signature could discern between immunity levels, and the cut-off points that were set could identify these differences.

To determine the effect of age, samples were divided into two groups, <60 and ≥60 years, as a study has shown that the incidence of EC increases significantly over the age of 60 (4). The risk scores were significantly different between the two groups; the ≥60 years group had higher risk scores compared with the <60 years group (Fig. 4E). The differences between grades 1 and 3, and grades 2 and 3 were significant, while there was no difference between grades 1 and 2. Therefore, if grades 1 and 2 are joined into a low-grade group and grade 3 is denoted as a high-grade group, the difference between the high- and low-grade groups is significant (Fig. 4F). Stage I of the FIGO stage was classed as the early stage, and stages II–IV were classed as the late stage. A difference in risk score was observed between the early and late stages; the higher the stage, the higher the risk score (Fig. 4G). The aforementioned results indicated the relationship between the risk score and clinical traits, and they also redefined and divided clinical traits according to the prognostic model to improve alignment of clinical traits with the FAM prognostic signature.

Based on the aforementioned analysis, a nomogram, a fusion model of gene expression levels and clinical traits, including age, grade and stage information, was established (Fig. 5A). ROC curves demonstrated considerable prognostic effects for 1, 3 and 5 years with an AUC of 0.772, 0.780 and 0.814, respectively (Fig. 5B). The calibration curve also confirmed the validity of the prediction (Fig. 5C). Numerous prognostic models have been developed to predict factors that may affect prognosis from different perspectives. To compare the effectiveness of the models created in the present study with previously established models, five established models were chosen at random. These were from the studies of Chen et al (34), Liu et al (35), Yin et al (36), Zhang and Yang (37) and Zhang et al (38). In terms of 1-, 3- and 5-year predictions, the nomogram of the present study ranked highest and in terms of 1- and 5-year predictions, the signature created in the present study ranked second (Fig. 5D-F). Both the FAM gene-based signature and nomogram model performed well in prognosis prediction, suggesting that FAM serves an important role in EC.

According to the median number of total scores from the nomogram, the samples were evenly divided into high- and low-nomo-score groups. Genes were mutated at different frequencies in samples. For the 20 genes with the highest mutation rates, the mutation percentage of the low-nomo-score group was higher than that of the high-nomo-score group, except for TP53 (Fig. 5G). The TP53 mutation occurred at a rate of 14% in the low-nomo-score samples, whereas the mutation rate was 56% in the high-nomo-score samples, with missense mutations accounting for the majority of mutations. The three genes with the highest mutation rates were PTEN, PIK3CA and ARID1A in both high and low nomo-score groups. The TMB of samples from different nomo groups was then calculated. In general, the low-nomo-score group had a higher TMB, with a statistically significant difference (Fig. 5H). Samples were divided into high- and low-TMB groups according to the median value of TMB and survival analysis demonstrated that the high-TMB group had a higher survival probability (Fig. 5I), which was consistent with the established conclusion that a high TMB is associated with prolonged survival after immune checkpoint inhibitor (ICI) treatment in several types of cancer (39). The combination of TMB classification and nomo-risk showed a higher survival rate for patients with a high TMB and a low nomo-risk score. Furthermore, patients with low TMB and high nomo-risk score had the least favorable outcomes (Fig. 5J). In summary, the aforementioned results indicated that low-nomo-score groups have a higher TMB, thus predicting an improved outcome of immunotherapy.

As ICI therapy strategies have been discussed previously (17), the relationship between the nomo-score and ICI genes was investigated. A total of 34 ICI-associated genes were identified and the results of correlation analysis revealed that ICI genes were correlated with the nomogram (Fig. 5K). Furthermore, the expression levels of ICI genes in the high- and low-nomo-risk groups were calculated. There were 17 genes with different expression levels between groups (Fig. 5L). Thus, a prevalent relationship between the ICI genes and the nomogram was observed. The MSI of EC also affected the efficacy of immunotherapy. The MSI-H group had a lower nomo-score, indicating improved immunotherapy efficacy for the low-risk group (Fig. 5M and N). All results indicated that low-nomo-score groups had improved survival rates and responded effectively to immunotherapy. A prediction of drug sensitivity was conducted based on gene signatures and nomo-score groups using the R package ‘pRRophetic’. There were 35 predicted sensitive drugs in both categories, and 29 of them were identified in both the gene signature and nomo-score groups (Table SV). These results demonstrated that the nomo-score evaluation system could predict the survival status and therapy effect, and that a low nomo-score indicated improved outcomes.