RNA sequencing data for uterine UCEC were obtained from the TCGA-UCEC cohort within TCGA database (https://portal.gdc.cancer.gov/). Raw count data were converted into transcripts per million (TPM) to construct the expression matrix. Following quality control and filtering, 584 samples (including tumor and adjacent normal tissues) and 60,616 genes were retained for downstream analyses. To normalize expression levels and reduce data skewness, the TPM matrix was log-transformed using the formula log2 (TPM +1) (Zhang et al., 2020). Genes with low-expression (mean TPM <1 across all samples) were excluded, yielding a final dataset suitable for subsequent analyses. To validate transcriptomic findings at the protein level, clinical and proteomic data for UCEC were retrieved from the CPTAC via the UALCAN portal (http://ualcan.path.uab.edu/) (Chandrashekar et al., 2017; Chandrashekar et al., 2022). In addition, transcriptomic data from the Genotype-Tissue Expression (GTEx) database were integrated to compare CCDC138 expression across 33 cancer types (Coorens et al., 2025).

ssGSEA was performed using the GSVA package (version 1.38.2) in R (version 4.2.1) to quantify the enrichment levels of specific biological pathways (Hänzelmann et al., 2013). Gene sets related to disulfidptosis, cuproptosis, and ferroptosis were curated from the Molecular Signature Database (MSigDB, version 7.5.1). Enrichment scores were computed for each sample, and hierarchical clustering was conducted to identify differential enrichment patterns between tumor and normal tissues. Heatmaps were generated using the pheatmap package (version 1.0.12) for visualization. Statistical significance was assessed using the Wilcoxon rank-sum test, with p < 0.05 considered significant.

To investigate gene co-expression patterns in TCGA-UCEC data, WGCNA was performed using the WGCNA package in R (Langfelder and Horvath, 2008). The log-transformed TPM expression matrix was used as input after filtering for genes with a mean TPM >1. Sample outliers were removed based on a WGCNA of 137,000. An appropriate soft-thresholding power (β) was determined using the pickSoftThreshold function to approximate a scale-free network topology. Gene modules were identified using the dynamic tree cut algorithm (cutreeDynamic), with the minimum module size set to 100. Module similarity was assessed through eigengene clustering, and correlations between gene modules and clinical traits were evaluated via correlation analysis. The blue module, which showed the strongest correlation with tumor traits, was selected for further analysis. Heatmaps and scatter plots were generated for data visualization.

DEGs between tumor and normal samples were identified based on the gene expression data derived from WGCNA. Gene expression values were normalized using the trimmed mean of M-values method and converted to counts per million. Genes with low expression were filtered out before analysis (Robinson and Oshlack, 2010). For each gene, fold change and false discovery rate (FDR) were computed, and those meeting the criteria of |log₂ fold change (log₂FC)| > 1 and FDR <0.05 were considered significantly differentially expressed (Zhao et al., 2021). Volcano plots were generated to visualize differential expression patterns, with significantly upregulated and downregulated genes annotated. Boxplots were used to compare the expression levels of selected DEGs between tumor and normal groups. The resulting DEG expression matrix was retained for subsequent analyses. DEG results were saved in text format, and the resulting DEG expression matrix was retained for downstream analyses.

To identify key feature genes among the 46 DEGs, five machine learning algorithms were employed: gradient boosting machine (GBM), random forest (RF), support vector machine-recursive feature elimination (SVM-RFE), least absolute shrinkage and selection operator (Lasso), and extreme gradient boosting (XGBoost) (Sanz et al., 2018; Tibshirani, 1996) (Becker et al., 2023; Chen and Guestrin, 2016). These analyses were conducted using the caret (version 6.0–90), randomForest (version 4.6–14), e1071 (version 1.7–9), glmnet (version 4.1–2), and XGBoost (version 1.5.0.2) packages in R. Each algorithm ranked genes based on their contribution to classification accuracy. Venn diagrams were generated using the VennDiagram package (version 1.7.1) to identify overlapping feature genes across all methods. The expression profiles of the consensus genes were retained for downstream analyses.

Overall survival (OS) analysis was performed using the survival (version 3.2–13) and survminer (version 0.4.9) packages in R. The prognostic relevance of candidate genes was assessed via Kaplan–Meier survival curves, with statistical significance determined by the log-rank test (p < 0.05). For CCDC138, patients were divided into high and low-expression groups based on the median expression level. Differential expression of CCDC138 between tumor and normal tissues was validated using boxplots generated with the ggplot2 package.

The transcriptional expression of CCDC138 across 33 cancer types was analyzed by integrating data from TCGA and GTEx (Tang et al., 2017). Boxplots were generated using the ggplot2 package in R to visualize differential expression. Protein level expression in UCEC was evaluated using CPTAC data assessed via the UALCAN portal. Expression differences across tumor stages, histological grades, and molecular subtypes were compared, with statistical significance determined using the Wilcoxon rank-sum test (p < 0.05). Subcellular localization of CCDC138 was assessed using immunofluorescence data from The Human Protein Atlas (https://www.proteinatlas.org/). Staining patterns in A-431 (epidermoid carcinoma), U-251MG (glioma), and U-2 OS (osteosarcoma) cell lines were examined to evaluate co-localization with cellular components, including nuclei, endoplasmic reticulum, and microtubules.

To elucidate signaling pathways associated with CCDC138, proteomic data from the CPTAC-UCEC cohort (99 tumor samples and 31 adjacent normal samples) were analyzed. Samples were grouped based on pathway activity status (“pathway-altered” vs. “other”), and differences in CDKN2A protein expression were assessed using the Wilcoxon rank-sum test. Key pathways analyzed included Hippo, NRF2, RTK, WNT, chromatin remodeling, mTOR, p53/Rb, and MYC/MYCN (Sanchez-Vega et al., 2018; Dou et al., 2020). Boxplots were generated using the ggplot2 package in R to visualize pathway-specific alterations.

The association between CCDC138 expression and immune cell infiltration was assessed using the CIBERSORT algorithm (version 1.03) and immune profiling data from the TIMER database (https://cistrome.shinyapps.io/timer/) (Newman et al., 2015; Li et al., 2016). Pearson correlation coefficients were computed for 68 immune cell types, with significance set at p < 0.05. Correlation patterns were visualized using bubble plots generated with ggplot2. TME components, including stromal and immune scores, were estimated using the ESTIMATE algorithm in R (version 1.0.13). Tumor purity was also calculated. Differences in TME metrics between high and low CCDC138 expression groups were evaluated using the Wilcoxon rank-sum test (p < 0.05).

To assess potential drug responses in TCGA-UCEC samples, TPM expression data were integrated with drug sensitivity profiles from the Genomics of Drug Sensitivity in Cancer database (Yang et al., 2013). Normal tissue samples and genes with low expression were excluded. The half-maximal inhibitory concentration (IC50) values for various therapeutic agents were predicted using a drug response modeling approach (Geeleher et al., 2014). Samples were categorized into high and low CCDC138 expression groups based on the median expression level. IC50 values were then compared between the two groups. Drugs showing significantly increased sensitivity in the high-expression groups (p < 0.001 and logFC <0) were visualized using boxplots.

The PPI network of CCDC138 was constructed using the GeneMANIA database (https://genemania.org/). Interacting proteins were identified along with their interaction types (e.g., physical, genetic) and interaction strengths. The network was visualized in Cytoscape software (version 3.9.1), with nodes representing proteins and edges denoting the type and strength of the interactions (Shannon et al., 2003).

Ishikawa cells were cultured in RPMI-1640 medium (Hyclone, Cat.No.SH30809.01B) supplemented with 10% fetal bovine serum (FBS, Hyclone, Cat.No.SH30087.01) and penicillin-streptomycin (Hyclone, Cat.No.SH30010) at 37°C with 5% CO₂. Cells were seeded at 5 × 10⁴ cells/well in 6-well plates. At 40% confluence, cells were transfected with siRNA targeting CCDC138 (si-CCDC138: ctcgactatgacatcaacattgadTdT) or negative control (si-NC: tcccgcgagacaacaccacctcadTdT) using Lipofectamine™ RNAiMAX (Invitrogen) in Opti-MEM. After 4–6 h, the transfection medium was replaced with complete medium. Transfection efficiency was assessed via RT-PCR after 24 h.

For CCK-8 assay, cells (1 × 10⁴/well) were seeded in 96-well plates. At 0, 24, 48, and 72 h, CellTiter96 AQueous One Solution (Promega, Cat.No.G3582) was added (10 µL/100 µL medium), incubated for 4 h, and absorbance measured at 490 nm using a microplate reader (Thermo Fisher Scientific, Multiscan MK3). Proliferation and inhibition rates were calculated.

For EdU assay, log-phase cells were seeded in 96-well plates, incubated with EdU (10 µM) for 2 h, fixed with 4% paraformaldehyde, and stained using an EdU kit. Nuclei were counterstained with DAPI. Images were captured using a Leica DMI6000B fluorescence microscope, and EdU-positive cells were quantified.

Total RNA was extracted using TRIzol (Invitrogen), treated with DNase I (Promega), and reverse-transcribed. RT-PCR was performed using SYBR Green qPCR SuperMix (Invitrogen) on an ABI PRISM® 7500 system with primers for CCDC138, S6K1, mTOR, p21, and GAPDH (internal control). Relative expression was calculated using the 2^−ΔΔCT method.

All statistical analyses were performed in R (version 4.2.1), with parallel processing enabled via the future.apply and doParallel packages. Continuous variables were analyzed using the Wilcoxon rank-sum test or t-test, as appropriate. Survival analysis was conducted using the log-rank test. For multiple comparisons, p-values were adjusted using the FDR method. Unless otherwise specified, p-values <0.05 were considered statistically significant.

RNA sequencing data from the TCGA-UCEC cohort were preprocessed to generate normalized expression matrices. After quality control and filtering, 584 samples and 60,616 genes were retained. Expression values were converted to the TPM matrix and subsequently log-transformed using the formula log2 (TPM +1). Genes with low-expression (mean TPM <1) were excluded, establishing a high-quality dataset for downstream analyses.

ssGSEA was performed to calculate enrichment scores for disulfidptosis, cuproptosis, and ferroptosis pathways. The resulting heatmap (Figure 1A) demonstrated distinct separation between tumor and normal samples, indicating pronounced molecular differences.

WGCNA identified six distinct co-expression modules using a soft-thresholding power of β = 16 (R ² > 0.9). After removing outliers, 584 samples were retained for analysis. The blue module, containing 1,800 genes, showed the strongest correlation with tumor traits and was selected for subsequent analysis (Figure 1B).

A total of 46 DEGs were identified within the blue module (|log2FC| > 1, FDR <0.05), comprising 28 upregulated and 18 downregulated genes. Significant DEGs were visualized using a volcano plot (Figure 1C), and expression differences were further validated through boxplots (Supplementary Material, containing boxplots of DEG expression). The resulting DEG expression matrix was retained for downstream analyses.

Five machine learning algorithms were applied to identify feature genes, yielding 19 from GBM, 18 from RF, eight from SVM-RFE, 20 from Lasso, and 16 from XGBoost. A Venn diagram revealed five overlapping genes—MAP6, CCDC138, DNAAF3, STX2, and FABP6 —which were selected for downstream analyses based on their consensus across models expression profiles (Figure 1D).

Survival analysis identified four genes significantly associated with OS (p < 0.05). Among them, high expression of CCDC138 was significantly correlated with poorer OS (p = 0.003, Figure 2A). Differential expression analysis further revealed that CCDC138 expression was significantly elevated in UCEC tumor tissues compared to that in normal tissues (Figure 2B), highlighting its potential as a prognostic biomarker.

To assess the biomarker potential of CCDC138 in UCEC, its transcriptional expression was analyzed across 33 cancer types using data from TCGA and GTEx databases (Figure 2C). CCDC138 was significantly upregulated in several cancers, including UCEC. Protein-level analysis using the CPTAC dataset via the UALCAN portal further confirmed elevated CCDC138 expression in UCEC tumor tissues compared to that in normal controls (31 normal vs. 100 tumor samples, Figure 3A). Stratified analysis revealed significantly higher protein expression across different UCEC stages (Figure 3B), histological grades (Figure 3C), and subtypes (Figure 3D), with marked elevation in early-stage disease and specific subtypes.

Immunofluorescence analysis of A-431 (epidermoid carcinoma), U-251MG (glioma), and U-2 OS (osteosarcoma) cell lines, as reported in the THPA, was used to investigate the subcellular distribution of CCDC138 (Figure 4A). CCDC138 exhibited co-localization with DAPI-stained nuclei and was detected in both the nuclear and cytoplasmic compartments in all examined cell lines.

To investigate the involvement of CCDC138 in UCEC pathogenesis, proteomic data from the CPTAC-UCEC cohort (99 tumor samples and 31 adjacent normal tissues) were analyzed. Samples were categorized into “pathway-altered” and “other” groups for key oncogenic pathways. Comparative analysis revealed significantly elevated CCDC138 protein expression in samples with alterations in the Hippo, NRF2, RTK, WNT, and chromatin remodeling pathways. In addition, CCDC138 expression was strongly associated with dysregulations of the mTOR, p53/Rb, and MYC/MYCN pathways (Figures 4B–J).

CCDC138 expression was significantly correlated with 68 immune cell types (p < 0.05), including a strong positive correlation with plasmacytoid dendritic cells (r = 0.484, p = 1.20e-33) and a negative correlation with natural killer T (NKT) cells (r = −0.378, p = 4.43e-20). A bubble plot illustrated the overall correlation patterns (Figure 5A), suggesting a potential role for CCDC138 in modulating the immune microenvironment. Further analysis using the TIMER database confirmed negative associations between CCDC138 expression and immune cell infiltration (Figure 5B), significant links with OS in patients with UCEC (Figure 5C), and the highest CCDC138 mutation frequency among TCGA cancer types (Figure 5D).

The high CCDC138 expression group exhibited significantly lower stromal scores (p = 3.92e-07) and immune scores (p = 1.41e-08; Figure 5E), with higher tumor purity (p = 1.88e-10; Figure 5F), indicating that CCDC138 may influence the TME and contribute to tumor progression.

Drug sensitivity analysis identified 81 compounds with significantly lower IC50 values in the high CCDC138 expression group (p < 0.001, logFC <0), indicating enhanced drug sensitivity. Boxplots for agents, such as 5-fluorouracil, acetalax, alpelisib, AZ6102, and AZD3759 (Supplementary Material, containing drug sensitivity boxplots) highlighted the potential of CCDC138 as a predictive biomarker for treatment response.

PPI analysis using the GeneMania database identified multiple CCDC138-interacting partners, including DCTN2, SSX2IP, CEP72, CEP290, CEP131, CEP120, OFD1, STX12, PMS1, PAK5, IKBKG, ANKRD26, RARS2, HYOU1, ZNF45, HEATR1, PLD5, ZBED5, EPHB6, and NLGN4X. Interaction types and strengths are summarized in Table 1 and illustrated in a network diagram (Figure 6).

To validate the bioinformatics findings, in vitro experiments were conducted using Ishikawa cells, a well-established UCEC model.