Identification of disulfidptosis-associated genes and characterization of immune cell infiltration in thyroid carcinoma

Source and processing of data

The genomic sequencing data for THCA were obtained from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). This study employed expression data from TCGA-THCA-FPKM, alongside pertinent clinical information from 507 THCA cases, encompassing demographics such as age, gender, stage, TNM classification, survival status, and duration. The characteristics of the THCA patients within the TCGA database are outlined in Supplementary Table 1. Before analyzing the expression profile data, transcripts per million (TPM) from the exon model per thousand bases underwent standardization. Drawing from prior literature [3, 4, 10, 11], a curated set of 106 DRGs was compiled and detailed in Supplementary Table 2.

Additionally, the THCA dataset GSE33630, comprising 49 THCA samples and 45 normal tissues, was accessed from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Leveraging the labeling information of the GPL570 platform, probes were converted into corresponding gene symbols. This dataset served as the validation cohort for our study.

Consensus clustering

Consensus clustering serves as a methodological cornerstone for discerning the optimal number of clusters across diverse clustering algorithms [12]. In our study, leveraging the expression profiles of DRGs, we employed the ConsensusClusterPlus R package to effectively stratify THCA samples into discernible subtypes. To ensure robustness, the ConsensusClusterPlus command was iterated 1000 times, thereby fortifying classification stability. The determination of the cluster number (k) was facilitated by the cumulative distribution function (CDF) of uniform clustering. Subsequently, survival analysis was meticulously conducted, juxtaposing Cluster C1 and Cluster C2 based on survival status and duration, utilizing the survival R package. Moreover, we undertook a comprehensive assessment of survival disparities concerning gender, disease stage, and TNM classification.

Differential analysis

The expression levels of DRGs were extracted from the TCGA-THCA dataset utilizing the limma package. Subsequently, the transcriptomic data of TCGA-THCA-DRGs underwent rigorous analysis employing the DESeq2 R package. Employing stringent criteria, wherein |log2FC| > 1.5 and a corrected P < 0.05 were utilized as screening thresholds, the DEDRGs were meticulously identified. Simultaneously, employing the same rigorous methodology, we endeavored to identify DEDRGs within the validation cohort. This comprehensive approach ensured consistency and reliability across datasets, facilitating robust comparative analyses.

Verification of DEDRGs in the different database

We employed the “Expression Boxplot” module within the GEPIA database (http://gepia.cancer-pku.cn/) [13] to evaluate the differential expression patterns of DEDRGs between normal tissues and THCA samples. This facilitated a comprehensive assessment of the transcriptional alterations of DEDRGs across distinct biological contexts. Additionally, we utilized the “Survival Plots” module within the same platform to delve into the prognostic implications associated with DEDRGs. This enabled a detailed exploration of the potential prognostic significance of DEDRGs in the context of THCA patient outcomes. Furthermore, to elucidate the precise cellular distribution and immunohistochemical (IHC) characteristics of DEDRGs, we conducted verification using the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) [14]. This comprehensive analysis provided valuable insights into the spatial localization and potential functional roles of DEDRGs within cellular frameworks. Moreover, the TIMER database (https://cistrome.shinyapps.io/timer/) [15] served as an essential resource for examining the intricate relationship between DEDRGs and immune cell infiltration. This facilitated a thorough interrogation of the immune microenvironment within the context of THCA, shedding light on potential immunomodulatory roles of DEDRGs [14].

Gene regulation analysis of DEDRGs

GenDoma represents a pioneering effort in China, constituting the premier comprehensive database encompassing human genes and diseases. By furnishing a multi-faceted analysis, GenDoma is poised to significantly bolster the realm of precision medicine. A central tenet of this study revolves around the Gene Regulatory Network (GRN), which serves as a compendium of molecular regulators governing gene expression at both mRNA and protein levels. In the pursuit of elucidating the regulatory landscape of DEDRGs, we harnessed the resources of the GenDoma. Through this platform, we accessed the GRN pertaining to DEDRGs, thereby unraveling the intricate web of regulatory interactions governing their expression dynamics.

GO and KEGG analysis

Utilizing the clusterProfiler R package [16], we conducted an in-depth analysis to elucidate the potential biological functions [17, 18] of the DEDRGs. Concurrently, we endeavored to explore their underlying molecular mechanisms. A significance threshold of P < 0.05 was applied to discern statistically significant findings. To present the results in a visually intuitive manner, we employed a diverse array of visualization techniques, including bubble charts, circle charts, and lollipop charts. These visualizations were meticulously crafted using the ggplot2 R package, facilitating the elucidation and interpretation of complex biological phenomena encoded within the DEDRGs.

Analysis of mutation landscape

Somatic mutation data pertinent to THCA were meticulously retrieved from the esteemed TCGA database in the form of Mutation Annotation Format (MAF) files. Subsequently, to offer a comprehensive visualization of the mutational landscape, a Waterfall Map was meticulously generated leveraging the Maftools R Package.

Correlation analysis of immune cell infiltration

In order to infer the abundance of 22 invasive immune cell types within THCA samples and normal tissues, we employed the CIBERSORT algorithm. CIBERSORT operates through a deconvolution analysis rooted in the linear support vector regression principle, thereby facilitating the prediction of the proportion of diverse cell types within gene expression profiles. A significance threshold of P < 0.05 was established as the standard for statistical significance [19]. The resultant data were visualized utilizing the barplot R package. Furthermore, the immune score and estimate score of THCA were meticulously computed utilizing the ESTIMATE algorithm. This allowed for a nuanced evaluation of the tumor microenvironment. To delve deeper into the correlation between DEDRGs and immune checkpoints within distinct subtypes, we employed the Pearson correlation coefficient to rigorously examine their association.

Construction and verification of prognostic model and nomogram

The glmnet R package served as a pivotal tool in conducting LASSO regression analysis, enabling the identification of pertinent variables associated with dependent variables (P < 0.05), while mitigating concerns related to variable collinearity [20]. This rigorous analysis facilitated the construction of a prognostic model by integrating survival time and survival status data for each sample. Subsequently, a corresponding risk score was meticulously generated, encapsulating the prognostic potential of the model. To comprehensively evaluate the predictive performance of the model, area under the curve (AUC) and receiver operating characteristic (ROC) analysis were employed, providing valuable insights into its efficacy.

Furthermore, the risk score, alongside age, gender, TNM, and stage variables, underwent rigorous assessment through univariate and multivariate Cox regression analyses. Variables exhibiting a significance level of P < 0.05 were identified as independent prognostic factors for THCA patient outcomes. Leveraging these independent prognostic factors, a nomogram was meticulously constructed utilizing the regplot R package. This nomogram facilitates the prediction of 1, 3, and 5-year survival rates of THCA patients, thereby offering a valuable clinical tool for prognostic assessment.

Bayesian co-localization analysis

Bayesian co-localization analysis plays a pivotal role in evaluating the presence of shared causal variation within two phenotypic regions. In our study, we employed the “coloc” R package to meticulously scrutinize the co-localization of seven DEDRGs with THCA, excluding FLNA due to the absence of GWAS data. This comprehensive analysis provides insights through four hypothetical posterior probabilities, elucidating the likelihood of a single variable being shared by both traits. To discern gene co-localization, we imposed a stringent criterion, stipulating that the PPH4 score exceed 0.75 [21]. The GWAS data for THCA (ebi-a-GCST90013867) were meticulously obtained from the IEU OpenGWAS platform (https://gwas.mrcieu.ac.uk/).

Experimental verification in vitro

Landscape of DRGs

To unravel the intricate interaction among various DRGs, we meticulously constructed a PPI network of DRGs (Figure 2A) utilizing the STRING database, a renowned resource for elucidating molecular interactions. Furthermore, the expression profiles pertinent to DRGs were meticulously visualized as a heatmap employing the pheatmap R package (Figure 2B), providing a comprehensive overview of the expression patterns across different DRGs.

Consensus clustering

Employing the Consensus clustering method, THCA patients were meticulously stratified into two distinct subtypes (Figure 3A, Supplementary Figure 1A–1D). The heatmap representing the two subtypes distinctly portrayed that cluster C1 exhibited a heightened expression level of DRGs, while cluster C2 showcased a comparatively lower expression level. This observation suggests a clear dichotomy, with cluster C1 corresponding to the high subtype of DRGs, and cluster C2 denoting the low subtype of DRGs (Figure 3B).

Survival analysis unequivocally demonstrated that cluster C1, characterized by a high DRGs level, boasted a notably favorable prognosis, while cluster C2, marked by a lower DRGs level, exhibited a considerably poorer prognosis (Figure 4A). Furthermore, survival curves based on clinical characteristics elucidated that females exhibited a superior prognosis compared to males (Figure 4B), while earlier TNM staging was distinctly associated with a more favorable prognosis for THCA patients (Figure 4C–4F), underscoring the prognostic significance of clinical staging in THCA management.

Differential analysis

The comparative analysis of TCGA-THCA-DRGs and GSE33630-DRGs datasets meticulously revealed statistically significant differences in DRGs expression between normal and THCA samples (P < 0.05). From the TCGA-THCA-DRGs dataset, a total of eight DEDRGs were discerned, wherein seven genes exhibited up-regulation while one gene displayed down-regulation (Figure 5A). Notably, Cluster C1 exhibited a distinctive hallmark of heightened expression of DEDRGs, while Cluster C2 was distinctly associated with a lower expression level of DEDRGs (Supplementary Figure 2A). The correlation analysis of TCGA-DEDRGs distinctly unveiled predominantly positive correlations among the DEDRGs (Supplementary Figure 2B).

In parallel, nine DEDRGs were successfully identified from the GSE33630-DRGs dataset, comprising five up-regulated genes and four down-regulated genes (Figure 5B). Remarkably, DEDRGs exhibited a markedly elevated expression level in the THCA group compared to the normal cohort (Supplementary Figure 2C). Moreover, the correlation analysis of GEO-DEDRGs underscored predominantly positive correlations among the DEDRGs (Supplementary Figure 2D). Notably, the intersection of DEDRGs from both datasets facilitated the identification of three co-DEDRGs (Figure 5C), highlighting potential candidates with robust implications across different datasets.

Verification of DEDRGs in the different database

The findings from the GEPIA database revealed significant differences in the expression levels of MYH10 and ME1 in THCA, with MYH10 exhibiting notably high expression and ME1 showing low expression (P < 0.05) (Supplementary Figure 3). Concurrently, insights gleaned from the HPA database illuminated the IHC staining patterns of the DEDRGs, showcasing distinct expression profiles between THCA samples and normal tissues (Figure 6). Specifically, MYH10 emerged as a potential risk factor for THCA, while ME1 and ZHX2 exhibited characteristics suggestive of protective factors against THCA development.

Moreover, Figure 7 provides a visual depiction of the specific spatial distribution of the DEDRGs within cellular contexts, unraveling deeper insights into their intracellular localization. Notably, FLNA, IQGAP1, and ACTN1 predominantly localize within the plasma membrane, hinting at their potential involvement in plasma membrane signaling and regulation. Conversely, ZHX2 and ME1 exhibit predominant localization within the nucleoplasm, suggesting roles in nuclear processes. Furthermore, MYH10 and PDLIM1 are primarily distributed within the actin filaments, indicative of their pivotal roles in actin filament signaling and receptor activation. Additionally, INF2 predominantly localizes within the endoplasmic reticulum, implicating potential roles in endoplasmic reticulum signaling. These distinctive spatial distribution patterns underscore the intricate regulatory mechanisms and diverse functional roles of the DEDRGs within THCA cells, offering valuable insights into their cellular localization and potential physiological functions.

Furthermore, analysis from the TIMER database unveiled significant correlations between DEDRGs and immune cell infiltration in THCA (Supplementary Figure 4A–4H), shedding light on the intricate interplay between DEDRGs and immune cell infiltration dynamics within the THCA microenvironment.

Gene regulation analysis of DEDRGs

Analysis of the gene regulation mechanisms associated with DEDRGs revealed intricate interactions with various miRNAs, transcription factors, and proteins, exerting regulatory control over the onset and advancement of THCA. Notably, the DEDRGs collectively modulate the activity of key transcription factors, including KLF1, JUN, and FLI1 (Figure 8A–8D). Furthermore, predictive modeling enabled the identification of potential targeted drugs and compounds, presenting promising prospects for future drug development and molecular therapy interventions aimed at disrupting the dysregulated pathways associated with THCA progression.

GO and KEGG analysis

The Gene Ontology (GO) function analysis, dissected into biological processes (BP), molecular functions (MF), and cellular components (CC) [22], unveiled significant enrichments, providing valuable insights into the functional roles of DEDRGs. Specifically, in BP, enrichment was evident in processes related to actin crosslink formation, cell junction organization, and exocytosis, highlighting their involvement in crucial cellular activities (Figure 9A). In CC, notable enrichment was observed in components associated with the actin filament and actin cytoskeleton, indicative of their pivotal roles in cytoskeletal organization and dynamics (Figure 9B). Meanwhile, MF analysis revealed significant enrichments in functions such as ADP binding, Ras GTPase binding, and actin filament binding, underscoring the diverse molecular functions of DEDRGs (Figure 9C). Furthermore, KEGG pathway analysis delineated the functional implications of DEDRGs in THCA pathogenesis. The analysis revealed significant enrichments in pathways including the MAPK, PPAR signaling pathway, and Proteoglycans in cancer, implicating their involvement in diverse signaling cascades and cellular processes critical for THCA progression (Figure 9D). These findings shed light on the intricate molecular mechanisms underlying THCA pathogenesis and highlight potential therapeutic targets for further investigation and intervention.

Immune cell infiltration analysis and mutation landscape

The results of immune cell infiltration analysis revealed that Cluster C2 exhibited higher immune cell infiltration compared to Cluster C1, indicative of distinct immune microenvironments within these subtypes (Figure 10A–10C). Notably, although Cluster C2 displayed elevated immune cell infiltration, it was associated with a poorer prognosis, while Cluster C1, characterized by lower immune cell infiltration, demonstrated a better prognosis. This observation led us to postulate that Cluster C2 represents an immune rejection type, whereas Cluster C1 corresponds to an immune inflammation type, suggesting divergent immunological responses in THCA (Figure 10A–10C). Furthermore, correlation analysis of immune cells unveiled intriguing relationships among different immune cell populations (Figure 10D).

To further elucidate the interplay between DEDRGs and immune checkpoints in THCA, a Wilcoxon signed-rank sum test confirmed significant upregulation of TIGIT, PD-L1, CTLA4, and CD274 in Cluster C1 (P < 0.0001), while PD-1 was highly expressed in Cluster C2 (P < 0.05). Interestingly, no statistical difference in the expression of LAG3 was observed between the two subtypes, suggesting differential regulation of immune checkpoint molecules in distinct THCA subtypes (Figure 11A–11F). The correlation heatmap of immune checkpoints revealed intricate relationships among these molecules (Figure 11G, 11H).

Moreover, analysis of the mutation landscape provided insights into susceptible genes in THCA. The waterfall map of gene mutations highlighted missense mutations as the predominant classification. Remarkably, the BRAF mutation rate (76.2%) emerged as the highest in THCA, with Cluster C1 exhibiting a higher BRAF mutation frequency compared to Cluster C2, whereas NRAS mutation frequency was higher in Cluster C2 than in Cluster C1 (Figure 11I).

Construction and verification of prognostic model and nomogram

LASSO regression analysis was performed to ascertain the optimal λ value through cross-validation employing 10 folds. Subsequently, five genes demonstrating significant prognostic effects were identified (Figure 12A–12D). Cox multivariate regression analysis was then employed to construct a risk model for calculating the risk score, formulated as follows: RiskScore = (0.319775943653634 × ANGPTL7) + (1.61175649349686 × FIRRE) + (5.33230222379885 × ODAPH) + (1.05543854580178 × PROKR1) + (0.196685509128214 × SFRP5). Based on the median risk score, all samples were stratified into high and low-risk groups. Kaplan-Meier survival curve analysis unveiled that the survival rate of the high-risk group was significantly lower compared to the low-risk group (Figure 12E). Furthermore, the prognostic model exhibited robust forecasting ability, as evidenced by the high AUC values for predicting the OS of THCA patients at 1, 3, and 5 years, which were 0.99, 0.89, and 0.89, respectively (Figure 12F).

The box diagram illustrated that prognosis genes ANGPTL7, FIRRE, and ODAPH were highly expressed in the high-risk group (P < 0.05), while PROKR1 and SFRP5 were highly expressed in the low-risk group (P < 0.0001) (Figure 12G).

Through univariate and multivariate Cox analysis, it was ascertained that age, stage, gender, and risk score served as independent prognostic factors for THCA patients. This underscores the efficacy of the prognosis model in assessing THCA risk. Additionally, a nomogram model was developed by integrating clinical factors to evaluate the prognosis of THCA patients (Figure 13A, 13B). Calibration curve analysis revealed a close fit between the calibration curve and the ideal curve, indicating high prediction accuracy of the model (Figure 13C). Similarly, consistent conclusions were obtained in the validation cohort, where significant differences in OS were observed between the high and low-risk groups, with patients in the high-risk group exhibiting poorer prognosis (Figure 13D).

Furthermore, we evaluated the relationship between risk scores and immune cell infiltration. The results indicated that patients with increased risk scores were negatively correlated with CD8+ T cells, CD4 memory cells, and M2 macrophages (Figure 14A–14I). This suggests that patients with low-risk scores could potentially benefit more from immunotherapy.

Bayesian co-localization analysis

In the GWAS co-location analysis between DEDRGs and THCA, our investigation revealed PDLIM1 (PPH4 = 0.78) and INF2 (PPH4 = 0.99) to be collocated with THCA. Notably, our analysis did not unveil any GWAS co-location findings for FLNA with THCA (Figure 15A–15G, Supplementary Table 4). This observation suggests that the remaining five DEDRGs (IQGAP1, MYH10, ACTN1, ZHX2, ME1) do not exhibit a genetic correlation with THCA, despite their documented involvement in its pathogenesis as elucidated in our study. We postulate that observational studies may be susceptible to acquired confounding factors, such as environmental influences and acquired genetic variations.

Experimental verification in vitro

qRT-PCR analysis was utilized to discern the disparity in PDLIM1 and INF2 mRNA expression levels between the normal and control groups. Remarkably, a discernible decrease was observed in the expression level of INF2 mRNA within the control group (FTC-133) compared to the normal group. Conversely, the expression level of PDLIM1 mRNA within the control group (FTC-133) exhibited a marked increase relative to the normal group. These disparities yielded statistically significant results (P < 0.05), as depicted in Figure 15H. Western blot results indicated a significant decrease in the expression of INF2 protein in the control group, while the expression of PDLIM1 protein in the control group showed a significant increase (Figure 15I), which is consistent with the qRT-PCR validation. These results further corroborate the findings of our study.