Application of Mendelian randomization and bioinformatic analysis to construct a prognostic model for thyroid cancer and perform pan-cancer analysis

Abstract

Objective: This study aimed to identify causal effects and potential molecular mechanisms of genes associated with THCA development.

Methods: Bioinformatic analyses were performed to identify differentially expressed genes (DEGs) associated with THCA. Subsequently, Mendelian randomization (MR) analysis was conducted using large-scale eQTL data and THCA GWAS summary statistics to screen for candidate genes. The intersection of DEGs and MR-derived candidate genes was used to determine DEGs with potential causal associations with thyroid carcinogenesis. Functional enrichment analysis, pathway analysis, and immune cell infiltration profiling were performed. External datasets were used for validation. Additionally, prognostic modeling and pan-cancer analyses of the candidate genes were conducted.

Results: IVW-based MR analysis revealed that elevated expression levels of ALOX15B [OR = 1.647, 95% CI (1.120-2.420), P < 0.05], TIAM1 [OR = 1.270, 95% CI (1.001-1.611), P < 0.05], and TMC6 [OR = 1.250, 95% CI (1.021-1.530), P < 0.05] were associated with an increased risk of THCA. Conversely, elevated expression of JUN [OR = 0.795, 95% CI (0.653-0.967), P < 0.05], PAPSS2 [OR = 0.779, 95% CI (0.608-1.000), P < 0.05], and RAP1GAP [OR = 0.895, 95% CI (0.810-0.989), P < 0.05] was associated with a reduced risk. Gene set enrichment analysis (GSEA) indicated that risk genes were enriched in proliferation- and metastasis-related pathways, such as extracellular matrix (ECM)-receptor interaction and cell adhesion molecules (CAMs). Findings from the training set were further validated experimentally and via external datasets. Additionally, candidate risk genes demonstrated associations with the development and progression of multiple tumor types.

Conclusion: This study identified ALOX15B, TIAM1, and TMC6 as potential risk genes and JUN, PAPSS2, and RAP1GAP as protective genes in THCA. These genes may serve as promising biomarkers and therapeutic targets for THCA, offering novel insights into precision oncology.

Keywords: Bioinformatics; Immune infiltration; Mendelian randomization; Prognostic model; Thyroid cancer.

Fig. 1

The flowchart graph of this study

Fig. 2

A: Heatmap of DEGs between THCA and control tissues; B: Volcano plot of DEGs; red represents up-regulated DEGs and blue represents down-regulated DEGs; C: Venn diagram showing the intersection of 420 up-regulated DEGs with 106 low-risk candidate genes, identifying three THCA suppressor genes; D: Venn diagram showing the intersection of 303 down-regulated DEGs with 83 high-risk candidate genes, identifying three THCA-promoting genes

Fig. 3

A₁: Summary results of the MR analysis; A₂: Circular plot illustrating the chromosomal distribution of the six core genes; B_1–B₆: Scatter plots of causal effects between six candidate genes and THCA in the MR analysis: B₁: JUN; B₂: PAPSS2; B₃: RAP1GAP; B₄: ALOX15B; B₅: TIAM1; B₆: TMC6; C₁–C₆: Forest plots showing the causal effects of SNPs associated with each gene on THCA: C₁: JUN-associated SNPs; C₂: PAPSS2-associated SNPs; C₃: RAP1GAP-associated SNPs; C₄: ALOX15B-associated SNPs; C₅: TIAM1-associated SNPs; C₆: TMC6-associated SNPs; D₁–D₆: Leave-one-out sensitivity analyses evaluating the stability of MR results for each gene: D₁: JUN; D₂: PAPSS2; D₃: RAP1GAP; D₄: ALOX15B; D₅: TIAM1; D₆: TMC6

Fig. 4

A, B: GO enrichment analysis of core genes showing involvement in BP, CC, and MF; C: KEGG pathway enrichment histogram—Y-axis indicates pathway names; X-axis indicates the number of core genes involved. D–I: GSEA plots of six core genes: D: JUN; E: PAPSS2; F: RAP1GAP; G: ALOX15B; H: TIAM1; I: TMC6

Fig. 5

Immune infiltration analysis and correlation with core gene expression. A: Proportion of 22 immune cell types in control versus THCA samples; B: Boxplot comparing immune cell infiltration between control (blue) and THCA (red) samples (* represents < 0.05, ** represents < 0.01, *** represents < 0.001); C: Immune correlation network. The lower left panel shows correlations among immune cell types; red indicates positive and blue indicates negative correlations. The upper right panel illustrates correlations between core genes and immune cell types; red lines represent positive correlations, green lines indicate negative correlations, and line thickness reflects the magnitude of the correlation coefficient

Fig. 6

Prognostic signature and immune cell correlation in THCA. A: Stratification of patients into high- and low-risk groups based on risk scores; B: Scatterplot of survival status and time versus risk score; C: KM curve comparing OS between high-risk and low-risk groups; D: Time-dependent ROC curve analysis for 1-, 3-, and 5-year OS; E: Correlation between risk scores and immune cell infiltration levels

Fig. 7

Differential Expression and Correlation of Core Genes in the Validation Dataset. A: Boxplot showing expression levels of intersecting genes in THCA (red) versus control tissues (blue). B: Heatmap illustrating pairwise correlations among the six core genes

Fig. 8

Validation of core gene expression in THCA and normal thyroid cell lines. A–F: Relative mRNA expression levels of ALOX15B, TIAM1, TMC6, JUN, PAPSS2, and RAP1GAP measured by qRT-PCR. G–I: Western blot analysis and quantification of TIAM1 and PAPSS2 protein expression in OR2, TPC, and IHH4 cells. TPC and IHH4 are THCA cell lines, while OR2 is a normal thyroid epithelial cell line. (* indicates comparison with OR2 p < 0.05; *** indicates comparison with OR2, p < 0.001)

Fig. 9

Pan-cancer expression analysis of THCA risk genes. A: ALOX15B expression across cancer types; B: TIAM1 expression across cancer types; C: TMC6 expression across cancer types