Ferroptosis-related gene transferrin receptor protein 1 expression correlates with the prognosis and tumor immune microenvironment in cervical cancer

TCGA, GTEx, and GEO datasets

Gene expression quantification data and clinical information on females with CC were downloaded from the TCGA database (Zhou et al., 2023). Gene expression data from healthy cervical tissues were obtained from the GTEx database (Chen et al., 2022). The validation cohorts, consisting of complete expression profile data (GSE63514 and GSE63678), were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/gds).

Identification and analysis of differentially expressed genes

Prior to conducting the differential expression analysis, we conducted background correction and quantile normalization utilizing the robust multi-array analysis (RMA) method implemented in the limma R package, resulting in the generation of the normalized gene expression matrix. P-value < 0.05 and log₂ fold-change > |1| were set as the cut-off criteria to identify the statistically significant differentially expressed genes (DEGs) (Shen et al., 2020; Yi et al., 2021; Zhen et al., 2020).

Identification of ferroptosis modules using weighted gene co-expression network analysis

We performed gene co-expression network analysis using the R package weighted gene co-expression network analysis (WGCNA) (Bhatia et al., 2021; Shi et al., 2023). Briefly, we converted the expression levels of individual transcripts into a similarity matrix based on Pearson’s correlation values between paired genes. To detect the outliers, we constructed a hierarchical clustering tree based on the expression matrix. To construct a scale-free network, we determined an appropriate soft thresholding power (β) as 6 by plotting R² (scale-free topology fitting index) against the soft thresholds. Next, we calculated the intramodular connectivity between genes exhibiting similar expression using the topological overlap dissimilarity measure. Finally, we utilized the dynamic hybrid cutting method to establish a hierarchical clustering tree and identify the co-expressed gene modules. Ferroptosis scores were obtained by applying Gene Set Variation Analysis to the ferroptosis genes present in the samples.

Identification of ferroptosis-related prognostic genes using least absolute shrinkage and selection operator Cox regression

In the next part of the study, we used the LASSO logistic regression method to identify the prognostic genes associated with ferroptosis in CC. The ‘glmnet’ R package was employed for the LASSO feature selection (Zhang, Pei & Zhu, 2023). Initially, 10-fold cross-validation was performed to determine the optimal value of the tuning parameter lambda, which influenced the shrinkage penalty applied to the regression coefficients. The risk score signature was defined as the combination of Coefi (representing the coefficient obtained from LASSO Cox regression) and Expri (Gao et al., 2022; Wang et al., 2021).

Functional annotation of FRGs

To investigate potential biological processes and signaling pathways associated with the identified FRGs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the R package ‘clusterProfiler’ (Ding et al., 2022). The crucial intersecting genes were subjected to this analysis. GO enrichment analysis elucidated three aspects: biological processes, cellular components, and molecular functions.

Multi-omics analyses of the FRGs

Mutation data for patients with CC were downloaded from the TCGA database. The relationship between the risk score and gene expression was visualized using the Kaplan–Meier (KM) curve, which was implemented using the KM R package ‘survival’ (Zengin & Önal-Süzek, 2021). Immunohistochemical images obtained from The Human Protein Atlas were used to determine the protein expression levels of the five key ferroptosis-related prognostic genes in both CC and normal cervical tissues (Xu et al., 2020b). Immunohistochemistry staining outcomes for five key genes were derived from the Human Protein Atlas (HPA) database.

Construction and validation of an FRGs-based prognostic signature for CC cohorts

We employed the KM method to evaluate the disparity in survival outcomes between the high-risk and low-risk groups. From the FRGs, we carefully selected those essential for constructing a prognostic risk score model to predict the overall survival (OS) of patients with CC. To confirm the predictive capability of the FRGs-based signature independently, we conducted univariate and multivariate Cox regression analyses using the “survival” package in R. These analyses involved assessing risk scores along with other clinical factors such as T stage, N stage, M stage, and age.

TFRC expression analysis of prognosis, model development, and assessment

The present study examined prognostic parameters, namely Disease Specific Survival (DSS), and Progress Free Interval (PFI) by analyzing patient data obtained from TCGA. The analysis was conducted within the clinical meaning module of the Xiantao platform (https://www.xiantao.love/), employing Cox regression and Kaplan-Meier methods. The threshold value for categorizing TFRC gene expression into low and high groups was determined based on the median value. To establish the association between clinical-pathological characteristics and TFRC gene expression, we employed the Wilcoxon signed-rank sum test in conjunction with logistic regression. The findings derived from the Cox regression model were then combined with the independent prognostic variables obtained from the univariate analysis. Subsequently, survival probabilities for 1, 3, and 5 years were projected using these integrated data. The accuracy of these projections was evaluated by comparing them to the actual occurrences through calibration curves.

Analysis of the infiltration of immune cells

We applied the ssGSEA method to investigate tumor infiltration through the analysis of 24 different immune cell types. The Spearman correlation algorithm was employed to compare immune cell infiltration levels between subgroups with high and low TFRC gene expression and to evaluate the strength of association between TFRC gene expression and the concentrations of the 24 distinct immune cell types. The relationship between TFRC gene expression and immune infiltration, as well as the association between infiltrating levels of immune cells and the values obtained in different subgroups of TFRC gene expression, were analyzed using the module provided by the “Xiantao tool” based on the findings associated with immune infiltration.

Patients and tissue samples

The CC tissues and the adjacent paracancerous tissues were obtained from the Department of Pathology, Second Affiliated Hospital, School of Medicine, Zhejiang University. The present study adhered to the ethical principles outlined in The Declaration of Helsinki. Study procedures were approved by the Ethics Committee of the Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou, China (approval no: 2023-1138). The Clinical Research Ethics Committee of the Second Affiliated Hospital, Zhejiang University School of Medicine, provided ethical approval and waived informed consent. A total of 33 patients with CC who had undergone surgical resection at the Second Affiliated Hospital of Zhejiang University School of Medicine were included in the present study. CC samples and the corresponding medical information were collected from all patients.

Immunohistochemistry staining

The protocol for IHC was essentially as we described previously (Li et al., 2022).

Quantitative real-time PCR

Total RNA was isolated utilizing the TRIzol reagent manufactured by Invitrogen (United States). The synthesis of complementary DNA (cDNA) was performed using the PrimeScript RT kit provided by Takara. Quantification of mRNA levels was performed using qRT-PCR. The primers were designed by Shanghai Bio-Tech, as follows:

TFRC: forward, 5′-AGGTCAAAGACAGCGCTCAA-3′ and reverse, 5′-GCCACATAACCCCCAGGATT-3′.

Actin: forward, 5′-AGCCTTCCTTCCTGGGCAT-3′ and reverse, 5′-CTGTGTTGGCGTACAGGTCT-3′.

Statistical analysis

Statistical and bioinformatics analyses were conducted using R software (version 4.2.0). We performed the Wilcoxon rank-sum test to examine the differential gene expression of TFRC in CC and normal tissues. Additionally, the Kruskal-Wallis test, logistic regression analysis, and Wilcoxon rank-sum test were used to investigate the relationship between TFRC gene expression and clinicopathological characteristics. The statistical significance of the observed variations was assessed through appropriate methods including the unpaired Student’s t-test, Spearman’s correlation analysis, Chi-square test, and Yates’ correction, depending on the nature of comparisons conducted. Finally, Table 1 presents the complete list of abbreviations in English.

Identification of ferroptosis-related modules in CC using WGCNA

To construct a WGCNA network, we initially calculated β and estimated the co-expression similarity to compute the adjacency. Network topology analysis was performed using the WGCNA PickSoft threshold function. Subsequently, β was set to 6 as the scale independence reached 0.9 and exhibited relatively high average connectivity (Fig. 1A). The gene network and distinct modules were established using the one-step network construction function in the ‘WGCNA’ R package. Specifically, a WGCNA network was constructed, resulting in the generation of eight gene co-expression modules (Fig. 1B). The heatmap in Fig. 1C shows the topological overlap matrix of all the analyzed genes. There was a substantial degree of independence between these modules, implying relative independence in gene expression. Furthermore, the yellow modulus significantly correlated with ferroptosis (Fig. 1D).

Identification of prognostic genes using LASSO analysis

We conducted an overlap analysis between the identified DEGs and the genes corresponding to the ferroptosis-related modules generated. A total of 3,514 DEGs and 644 ferroptosis-related module genes were selected for further analysis (Fig. 2A). Subsequently, using LASSO Cox regression analysis, we identified five key genes with the best prognostic value–KIFC2, TFRC, SMYD2, OPN3, and ASMTL. This method helped in reducing the dimensionality and calculating the correlation coefficients among genes (Fig. 2B).

Functional enrichment analysis of the ferroptosis-related genes identified in the multigene signature

To investigate the biological functions of DEGs related to ferroptosis, we conducted GO annotation and KEGG enrichment analyses. GO analysis revealed significant enrichment in several biological processes, including chromosome segregation, organelle fission, nuclear division, nuclear chromosome segregation, and sister chromatid segregation. Among the cellular components, the enriched terms mainly included spindle, chromosomal region, and condensed chromosomal structures. Among the molecular functions, there was significant enrichment of the DNA-binding transcription activator activity, tubulin binding, microtubule binding, and cyclin-dependent protein serine kinase regulator activity (Figs. 3A, 3C and Tables 2, S1). Additionally, the GO analysis suggested that changes in the chromatin endoreduplication and segmentation complex were enriched in nuclear division, chromatin binding, and condensed chromosomes. Furthermore, KEGG enrichment analysis indicated the involvement of seven main pathways: cell cycle, cellular senescence, the p53 signaling pathway, oocyte meiosis, carbon metabolism, amino acid biosynthesis, and glutathione metabolism (Fig. 3B and Table S2). These results align with the established functions associated with ferroptosis in CC, thereby corroborating our findings.

FRGs in the multigene signature displayed an association with prognostic value for CC

As seen in the volcano map given in Fig. 4A, the tumor samples displayed significant upregulation of KIFC2, TFRC, SMYD2, and OPN3 compared to the normal samples. Conversely, ASMTL was significantly downregulated in the tumor samples when compared to the normal samples. Additionally, we investigated the mutation landscape of the FRGs involved in CC. Among the 289 samples analyzed, 3.81% exhibited mutations in at least one key molecule. ASMTL, KIFC2, and TFRC displayed low mutation frequencies, whereas SMYD2 and OPN3 showed no mutations in the CC samples. The waterfall plot shown in Fig. 4B, which represents the mutation landscape of the five signature molecules, revealed that missense mutations were predominant. Furthermore, patients with CC who displayed low expression levels of KIFC2, TFRC, SMYD2, and OPN3 had a significant survival advantage over those with high expression levels of the same (Fig. 5A). Interestingly, high ASMTL expression conferred a significant survival advantage. Moreover, the expression levels of KIFC2, TFRC, SMYD2, and OPN3 were significantly higher in the CC tissues than in the normal tissues (Fig. 5B). Immunohistochemical staining results for the five genes obtained from the Human Protein Atlas database revealed that the tumor group displayed notably higher protein expression levels of KIFC2, TFRC, SMYD2, and OPN3 than the normal group (http://www.proteinatlas.org/) (Fig. 5C); conversely, ASMTL was downregulated in the tumor group when compared to that in the normal group. The mRNA expression levels of the five genes in the external validation datasets (GSE63514 and GSE63678) demonstrated similar findings to those in TCGA. With the exception of ASMTL, we observed upregulated mRNA levels of KIFC2, SMYD2, TFRC, and OPN3 in tumor tissues within these datasets (Figs. 5D and 5E).

The FRGs-based prognostic model for CC

The patients with CC were stratified into high-risk and low-risk groups. Kaplan-Meier survival analysis demonstrated that the high-risk group exhibited poorer prognostic outcomes compared to the low-risk group within the TCGA dataset (Fig. 6A). Consistent with the findings from the TCGA training set, the high-risk groups demonstrated poor prognostic outcomes in the validation sets. The patients were stratified into high-risk and low-risk cohorts using a carefully determined cutoff value for the risk score, which was efficiently established through the MaxStat R package (Fig. 6B). Additionally, we examined the distribution of expression patterns of five FRGs, risk scores, survival time, and survival status in the TCGA dataset. Our results demonstrated the significant prognostic value of the ferroptosis-related signature (Fig. 6C).

Univariate and multivariate Cox analysis of the VMTRG signature

Univariate and multivariate Cox regression analyses were performed to establish if the five-gene signature is an independent predictor for OS of CC patients. Univariate Cox regression analysis of TCGA data revealed that T, N, M, riskscore and FRGs are independent prognostic variables for CC patients (Fig. 6D). Multivariate Cox regression analysis also showed that riskscore and FRGs are independent prognostic factors for CC patients in the TCGA dataset (Fig. 6E).

Correlation between TFRC gene expression and whole gene expression patterns

We conducted an analysis of the TFRC gene expression profile to enhance our comprehension of the biological significance of the TFRC gene in CC. The gene expression heat map presented the top nine genes exhibiting aberrant expression levels (with a logFC > |1| and P < 0.05) (Fig. 7A). To investigate the biological and functional pathways associated with differential TFRC gene expression, we employed GSEA using data from TCGA. This analysis allowed us to determine the pathways that distinguish between high- and low-TFRC gene expression groups. The enrichment signaling pathway that exhibited the highest relevance to TFRC gene expression was selected based on the calculated normalized enrichment scores. The GSEA analysis revealed a predominant concentration of the TFRC gene expression phenotype in the ferroptosis, diseases of immune system, jak stat signaling pathway, ECM affiliated, RHO GTPases activate NADPH oxidases, JNK pathway, and condensation of chromosomes (Fig. 7B). Our findings indicate that TFRC likely plays a significant role in the progression of CC.

Prognostic relevance of TFRC expression in CC

The TNM stage has long been recognized as an independent prognostic factor for CC survival. The T refers to the primary tumor, N refers to regional lymph node involvement, and M refers to distant metastasis. The term “OS event” denotes the duration spanning from treatment commencement to the patient’s demise. The relationship between TFRC expression and clinical parameters was assessed using the Kruskal-Wallis and Wilcoxon signed-rank tests. Elevated levels of TFRC expression were found to be associated with higher T stage, and N stage, as well as occurrence of OS events (Figs. 8A–8C). Table 3 displays the correlation between the clinicopathological characteristics of 306 CC patients and their TFRC protein levels. Patients exhibiting high TFRC expression demonstrated poorer survival outcomes. These findings provide evidence suggesting an association between increased TFRC expression and advanced tumor T stage, thus indicating a crucial role of TFRC in CC prognosis. We developed a clinical prognostic risk score for CC incorporating T stage, N stage, M stage, clinical stage, histological grade, age, and TFRC expression (Fig. 8D). Additionally, a calibration chart was utilized to evaluate the accuracy of the model’s predictions (Fig. 8E). TFRC expression has the potential to offer improved accuracy in predicting patients’ survival probabilities for 3- and 5-year intervals. In general, there was a demonstrated correlation between TFRC expression and the prognosis of patients with CC.

Overall, TFRC expression was shown to correlate with the prognosis of patients with CC. The figures illustrate the associations between TFRC expression and prognosis indicators using data from the TCGA database, including DSS and PFI. Increased TFRC expression correlated with poor outcomes in terms of DSS (HR = 2.10 (1.20–3.66), P = 0.009, Fig. 8F) and PFI (HR = 1.93 (1.19–3.13), P = 0.007, Fig. 8G).

The association between TFRC gene expression and immune cell infiltration

This study investigated and analyzed the relationship between TFRC gene expression and 24 distinct immune cell subtypes in CC. There was a significant positive correlation between TFRC gene expression and infiltration of Th2 cells, Tgd cells, and Tcm cells, while a strong inverse correlation was observed with T cells, pDCs, and cytotoxic cells, among others (Fig. 9A). Subsequent analysis revealed significant variations in TFRC gene expression levels across various infiltrating immune cell types, including neutrophils, NK CD56bright cells, pDCs, T cells, iDCs, DCs, cytotoxic cells, CD8 T cells, B cells, aDCs, Treg cells, and TH1 cells, among others (Figs. 9B–9D). Moreover, this study examined various functional subsets of T cells, including Tregs, exhausted T cells, Th1, Th2, Th9, Th17, Th22, and Tfh cells.

Experimental verification of the expression of TFRC

The transferrin receptor TFRC serves as a crucial component of the channel responsible for facilitating the uptake of iron ions into cells. It plays a pivotal role in regulating cellular iron metabolism and maintaining iron balance. To investigate the variations in the mRNA and protein expression of TFRC between paracancerous tissue and CC tissues, we examined the expression levels of TFRC in tissues by using qRT-PCR and IHC (Figs. 10A–10C). All procedures for cDNA synthesis, and qPCR were designed to follow the MIQE guidelines. The findings indicate that TFRC showed significantly higher expression levels in CC tissues compared to the adjacent paracancerous tissues.

Prognostic significance of TFRC in CC

Patient follow-up commenced from the resection date and extended until October 2023, with survival or death serving as the definitive endpoint for overall survival assessment. The prognostic significance of TFRC was assessed using Kaplan-Meier (KM) survival analysis. Patients with CC who displayed low TFRC expression levels exhibited a significant survival benefit in comparison to those with high TFRC expression levels (Fig. 10D). This implies that TFRC holds promise as a prognostic biomarker in CC and could potentially be targeted for therapeutic interventions in CC treatment.