Identification of iron metabolism-related genes as prognostic indicators for papillary thyroid carcinoma: a retrospective study

Abstract

Background: The thyroid cancer subtype that occurs more frequently is papillary thyroid carcinoma (PTC). Despite a good surgical outcome, treatment with traditional antitumor therapy does not offer ideal results for patients with radioiodine resistance, recurrence, and metastasis. The evidence for the connection between iron metabolism imbalance and cancer development and oncogenesis is growing. Nevertheless, the iron metabolism impact on PTC prognosis is still indefinite.

Methods: Herein, we acquired the medical data and gene expression of individuals with PTC from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database. Typically, three predictive iron metabolism-related genes (IMRGs) were examined and employed to build a risk score (RS) model via the least absolute shrinkage and selection operator (LASSO) regression, univariate Cox, and differential gene expression analyses. Then we analyzed somatic mutation and immune cell infiltration among RS groups. We also validated the prognostic value of two IMRGs (SFXN3 and TFR2) by verifying their biological function through in vitro experiments.

Results: Based on RS, all patients with PTC were stratified into low- and high-risk groups, where Kaplan-Meier analysis indicated that disease-free survival (DFS) in the high-risk group was much lower than in the low-risk group (P < 0.0001). According to ROC analysis, the RS model successfully predicted the 1-, 3-, and 5-year DFS of individuals with PTC. Additionally, in the TCGA cohort, a nomogram model with RS was developed and exhibited a strong capability to anticipate PTC patients' DFS. In the high-risk group, the enriched pathological processes and signaling mechanisms were detected utilizing the gene set enrichment analysis (GSEA). Moreover, the high-risk group had a significantly higher level of BRAF mutation frequency, tumor mutation burden, and immune cell infiltration than the low-risk group. In vitro experiments found that silencing SFXN3 or TFR2 significantly reduced cell viability.

Conclusion: Collectively, our predictive model depended on IMRGs in PTC, which could be potentially utilized to predict the PTC patients' prognosis, schedule follow-up plans, and provide potential targets against PTC.

Keywords: Iron metabolism; Nomogram; Papillary thyroid cancer; Prognostic signature; Risk score.

Figure 1. Flow chart.

Figure 2. Expression and correlation analysis of IMRGs.

(A) The volcano plot of comparison of IMRGs expression levels between the tumor and the control groups. The green and orange dots represented the down-regulated and up-regulated genes, respectively. (B) The IMR-DEGs expression level in different samples is shown in a heatmap. (C) The correlation matrix of IMR-DEG expression levels. Blue revealed a negative association, and red revealed a positive association. The darker the color, the greater the degree of correlation. Those without statistical significance were indicated in black.

Figure 3. Pearson’s correlation analysis point line diagram of IMRGs.

(A–I) The significant correlation analysis outcomes of IMRGs expression level, respectively, and the point line diagram of results were shown with correlation coefficients more than 0.45.

Figure 4. Interaction network analysis of IMR-DEGs.

(A) The PPI network of 12 IMR-DEGs, blue for down-regulated genes, red for up-regulated genes, and the bigger the circle, the greater the network weight. (B) The interaction network of IMR-DEGs and transcription factors, the core network was displayed after the source network was simplified by the minimum network algorithm. (C) The interaction network of IMR-DEGs and miRNA, the core network was displayed after the source network was simplified by the minimum network algorithm. (D) The interaction network of IMR-DEGs and small molecular compounds, the core network was displayed after the source network was simplified by the minimum network algorithm. (E–G) The interaction network of IMR-DEGs and drugs.

Figure 5. Survival analysis of IMRGs and LASSO regression was employed to examine predictive markers.

(A) Cox regression analysis of IMRGs expression level impact on the patient’s prognosis, displayed in the forest map. (B) The fitting result of LASSO regression used the ROC curve to assess the predictive ability of prognosis, with AUC as the area under the curve. (C and D) The LASSO Cox regression model was employed to screen predictive markers, and the incomplete likelihood deviation with 10 times cross-validation was utilized to calculate the best λ. (E–G) The high-/low- expression groups’ survival curves (grouped by median) with significant results of Cox analysis, respectively.

Figure 6. Evaluation of the prognostic ability of RS for prognosis survival time of PTC patients.

(A) ROC curve and calculated AUC value of RS for predicting 1-, 3- and 5-year survival. (B) The calculation dot plot of the best cutoff value of RS, with the cutoff value marked by a dotted line. (C) Survival curve for high/low expression groups of RS (K–M method). (D) Cox regression analysis of the RS, grouping, and other medical characteristics effects on the prognosis of patients, presented as a forest map.

Figure 7. Assessment and validation of the predictive ability of candidate prognostic markers and RS for clinical prognosis of subjects with PTC.

(A) The three prognostic markers’ expression differences were analyzed in the TCGA validation set between PTC and normal tissues. Kruskal-Wallis test was employed to examine the variations between groups, and asterisks (*P < 0.05, ***P < 0.001, ****P < 0.0001) indicated the statistically significant differences. (B) Differential expression analysis of three prognostic markers between PTC and healthy tissues in the GEO dataset GSE33630. (C) Survival curve of high-/low- expression groups of RS in TCGA validation set (K–M method).

Figure 8. Comparison of clinical characteristics of candidate prognostic markers.

(A–C) The boxplots of LTF, SFXN3, and TFR2 expression in different subgroups with different clinical characteristics, respectively. Wilcoxon or Kruskal-Wallis test was employed to evaluate variations between groups, and asterisks (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001) indicated statistically significant differences.

Figure 9. Multivariate survival analysis of RS and construction of the prognostic model.

(A) The Cox model of T and N stages, in addition to RS, is presented in a forest map. (B) The Cox model of tumor stage and RS is presented in a forest map. (C–F) The nomogram and calibration curve of a multivariate Cox regression model for predicting patient survival with RS.

Figure 10. DEGs analysis and functional enrichment analysis in RS grouping.

(A) Differential analysis volcano plot, blue dots for down-regulated, red dots for up-regulated. (B) Heatmap of differential analysis. (C) The bubble chart of KEGG analysis. The closer the color was to red, the smaller the P was, and the larger the bubble was, the more DEGs were enriched in this pathway. (D) The bar chart of KEGG analysis. The horizontal axis exhibited the number of genes enriched by the pathway, and the smaller the P, the closer the color was to red. (E) The network diagram of KEGG analysis, displaying the relevant genes in the neuroactive ligand-receptor interaction, mineral absorption, and tyrosine metabolism pathways. (F) The network diagram of KEGG analysis, displaying the relevant genes in the neuroactive ligand-receptor interaction, mineral absorption, and tyrosine metabolism pathways and representing the fold change values of differential analysis in colors. Among them, red represented up-regulation, purple represented down-regulation, and darker colors indicated larger values. (G) The enrichment analysis of BP, CC, and MF by GO analysis is displayed in a bubble chart. More DEGs were enriched in this pathway. The closer the color was to red, the lower the P, and the bigger the bubble. (H) The enrichment analysis of BP, CC, and MF by GO analysis is displayed in a bar chart. The horizontal axis revealed the number of genes enriched by the pathway, and the lesser the P, the closer the color was to red.

Figure 11. Correlation analysis of RS grouping and somatic mutation and immune cell infiltration of candidate genes.

(A) Genes with different somatic mutation proportions between high- and low-risk scores groups. (B) Comparison of TMB between high- and low-risk scores groups. Differences between groups were assessed by the Wilcoxon test. (C) Comparison of immune cell infiltration degree between high- and low-risk scores groups. Differences between groups were assessed by the Wilcoxon test, and asterisks (*P < 0.05, **P < 0.01) indicated statistically significant differences. (D) Correlation matrix between candidate gene expression and degree of immune cell infiltration. Positive correlations were represented by red, negative correlations by blue, and darker colors indicated a stronger correlation. A grid within the matrix displayed the correlation coefficient and P-value.

Figure 12. IMRGs expression verification and functional analysis.

(A) In normal and PTC tissues, IHC staining was employed to analyze SFXN3 and TFR2 expression. (B) Verification of SFXN3-siRNA and TFR2-siRNA silencing efficacy at the mRNA level in TPC-1 cells. (C) Cell viability assay was measured after SFXN3 or TFR2 silencing. (D-E) mRNA and protein levels of iron-related proteins expression after SFXN3 or TFR2 silencing. DMT1 (−) IRE, divalent metal transporter 1 without iron responsive element; DMT1 (+) IRE, divalent metal transporter 1 with iron responsive element. *P < 0.05, **P < 0.01, ***P < 0.001.