Blood Clotting Genes Linked to Thyroid Cancer Outlook, Immune Response, Stage, and Drug Sensitivity

We gathered and screened data from patients admitted to the First Affiliated Hospital of Zhengzhou University from March 1, 2024, to May 1, 2024. The inclusion criteria were: (1) first-time thyroid surgery; (2) histologically confirmed papillary thyroid cancer (PTC); (3) age >18 years. The exclusion criteria were: (1) presence of other malignancies; (2) incomplete clinical data; (3) presence of diseases related to abnormal coagulation levels (including venous thromboembolism, disseminated intravascular coagulation, etc.). Finally, 408 patients were involved in our research. Age, gender, presence of Hashimoto’s thyroiditis (HT), presence of extrathyroidal extension (ETE), multifocality of PTC, primary tumor size, and presence of lymph node metastasis (N stage, including central lymph node metastasis (CLNM) and lateral lymph node metastasis (LLNM)) were extracted from medical records. All patients underwent coagulation marker tests within two weeks before surgery. The antibody positivity often precedes clinical manifestations of thyroid dysfunction or sonographic changes in patients with HT. Studies have shown that elevated anti-TPO or anti-Tg antibodies can be present for years before the development of overt hypothyroidism or characteristic ultrasound changes, making antibody testing a valuable early diagnostic tool (9). HT was diagnosed by postoperative sectioning and examination of paraffin-embedded thyroid tissue specimens. Additionally, Serum antithyroglobulin and antithyroid peroxidase levels were measured within 30 days before surgery using the immuno-electrochemiluminescence method, and the patients were diagnosed with HT when these levels exceeded 115 IU/ml and 34 IU/ml, respectively. ETE referred to breaking through the thyroid capsule and invading adjacent soft tissues, muscles, trachea, oesophagus, nerves or blood vessels. Bilateral disease and multifocal disease were considered together for statistical analysis. Multifocality was defined as the presence of more than one lesion observed (10). Primary tumor size greater than 2 cm has been identified as an important risk factor for recurrence and lymph node metastasis in papillary thyroid carcinoma (3).

The flow chart is shown in Figure 1. The RNA-seq data, clinical information and mutation data of THCA were acquired from the TCGA-THCA. After removing data with missing or less than 30 days of follow-up, missing outcome time occurrence status, and duplicates, the TCGA-THCA set remained 500 THCA samples. Based on the GeneCards website (https://www.genecards.org/↗), we retrieved 378 coagulation-related genes (CRGs, Relevance Score≥3) via searching the term “coagulation” (7, 11).

The “limma” R package was conducted to figure out differentially expressed genes (DEGs) between normal and tumor tissues (FDR < 0.05 and log2FC >=1) in TCGA-THCA sets. Then, we took the intersection of these DEGs and CRGs to obtain differentially expressed CRGs (DECRGs). To explore the prognostic value of DECRGs in THCA, univariate Cox analysis was conducted in TCGA-THCA set (p < 0.05). Te co-expression network of the prognostic DECRGs was explored by the “igraph” R package. Using the STRING (https://string-db.database↗, the protein-protein interaction (PPI) network org/) of proteins coded by the mitochondrial dynamic prognostic DECRGs was constructed and visualized. Finally, we validated these prognostic DECRGs using survival analyses and log-rank tests.

Based on these prognostic DECRGs, we used the “ConsensusClusterPlus” package to perform consensus clustering on THCA samples from the TCGA-THCA dataset. To determine the optimal number of clusters, we used the Consensus Cumulative Distribution Function (CDF) to depict the CDF distribution patterns at different cluster numbers (k), and by plotting the Delta area plot to visually show the rate of change in the area under the CDF curve as the number of clusters increases from k to k+1, using this as a basis to select the optimal number of clusters and further subdivide TC patients into different subtypes. Then, we used the “survminer” package to plot Kaplan-Meier survival curves to analyze the differences in progression-free interval (PFI) between different subtypes.

We randomly split the samples from the whole set into a train set and a test set according to a 6:4 ratio using the R package “caret”. To further compress these prognostic DECRGs, we performed LASSO analysis using 10-fold cross-validation (12). Finally, we constructed a prognostic model using stepwise multivariate Cox regression analysis. The scores for each sample were calculated according to the following formula: Risk score = coefficient1 * gene 1 expression +…+ coefficientN * gene N expression, and the samples were classified into high or low risk groups based on the median score.

Kaplan-Meier analysis with chi-squared test was applied to estimate survival differences between risk groups in the model. The test set, and whole set were applied to validate the internal stability of the model. The time-dependent receiver operating characteristic (ROC) curves were employed to analyze the predictive performance of the model and clinical characteristics through the R package “timeROC” (13). To assess the applicability of the model to patients with diverse clinical characteristics, we compared the survival differences between different risk score groups within each subgroup. To assess whether the model was an independent predictor for predicting patient prognosis, we performed univariate and multivariate Cox analyses including risk scores and clinical characteristics. To assess the clinical relevance of the model, the study examined the relationship between groups and clinical features. To better apply the model to clinical work, we constructed a nomogram combining the model and clinical features to predict the 3-, 5-, and 7-year survival probabilities of patients.

To explore the differences in somatic mutations between different risk groups, we performed analysis and visualization using the R package “maftools” (14). To find differences in molecular mechanisms and relevant pathways between risk groups, we recognized differentially expressed risk genes (DERGs) between risk groups (|logFC > 1| and FDR < 0.05) and conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (p < 0.05) (15–18).

Furthermore, we employed the EPIC, MCP-COUNTER, TIMER, XCELL, QUANTISEQ, CIBERSORT, and CIBERSORT-ABS algorithms to examine immune cell infiltration (19–25). Subsequently, we scored the immune function of different risk groups using single-sample gene set enrichment analysis (ssGSEA) and compared the scores using the Wilcoxon test (26, 27). We also explored the differences in immune checkpoint gene expression levels between groups using the wilcox test. To predict the response to immunotherapy in each risk group, we calculated the tumor immune dysfunction and exclusion (TIDE) scores for each group and compared them (28). Finally, we used ESTIMATE to compare TME scores (StromalScore, ImmuneScore, and EstimateScore) across different risk groups.

In this study, the calcPhenotype function in the “oncoPredict” R package was utilized to predict drug sensitivity in the Genomics of Drug Sensitivity in Cancer (GDSC) database (29).

Using qPCR technology to detect the differential expression of genes used to construct a prognostic model in human papillary thyroid carcinoma cell lines (IHH4, KTC-1, TPC-1) and normal thyroid cell lines (Nthy ori-3-1). First, TRIzol reagent (Invitrogen, USA) was used to extract total RNA from the cells. Subsequently, the extracted RNA was reverse-transcribed into cDNA using the Prime Script RT reagent kit (Takara, Japan). After obtaining the cDNA template, quantitative PCR was carried out according to the instructions of the SYBR Green Quantitative PCR Detection Kit (Takara, Japan). All the reactions were repeated for at least 3 times. GAPDH was used as an internal control gene, and the 2−△△Ct method was used for the relative quantification of the target genes. The primer sequences used in this study were: AZU1, Forward: 5′-AGAACCTGAACGACCTGATGC-3′ and Reverse: 5′-CCTGGGAAAACGGGAGAGA-3′; COL3A1, Forward: 5′-TTGCTGTGGTGGTGTTGGAG-3′ and Reverse: 5′-TTCTAGCGGGGTTTTTACGA-3′; CP, Forward: 5′-AGGAGATTCGGTCGTGTGGT-3′and Reverse: 5′-TTGAGGGAAGAGGTTTGCTG-3′; CSF2, Forward: 5′-AGAGACACTGCTGCTGAGATG-3′and Reverse: 5′-CAGGAAGTTTCCGGGGTT-3′;F12, Forward: 5′-AGGACCAGCGATGGGGATA-3′ and Reverse: 5′-TGTGGAAAAACCGGAGAAGC-3′; GNA14, Forward: 5′-TGTTACGACAGGAGGAGGGA-3′ and Reverse: 5′-CGAAGCACATCTTGTTGGGT-3′; IL1RN, 5′-AACAGAAAGCAGGACAAGCG-3′ and Reverse: 5′-CCTTCGTCAGGCATATTGGT-3′; SERPIND1, Forward: 5′-ATGGGTATGATTTCCTTAGGTCTG-3′ and Reverse: 5′-GGAAGAGATTATGAATGGTCGTG-3′; GAPDH, Forward: 5′-GGCAAATTCCATGGCACCG -3′ and Reverse: 5′-TCGCCCCACTTGATTTTGGA-3′.

All data analyses were completed using R software (version 4.3.9) and GraphPad Prism (version 8.0.2). ROC curves were applied to evaluate the area under the curve (AUC) for each coagulation index. We used univariate and multivariate logistic regression analyses to determine the association between clinical characteristics and LLNM. Spearman correlation analysis was used to evaluate the correlation between two continuous variables. The Wilcoxon rank-sum test or Student’s T-test was used for the statistical analysis of two groups of continuous variables. Kaplan-Meier survival curves and log-rank tests were used to assess differences in DFS between different groups. P values less than 0.05 were considered to indicate statistical significance. ns indicates P>0.5, * indicates P<0.05, ** indicates P<0.01, *** indicates P<0.001.

We collected 408 cases of adult PTC, and their clinicopathological characteristics are shown in Figure 2A. The mean of coagulation profiles between LLNM and non-LLNM patients was compared in Figure 2B. The ROC curves revealed that D-dimer had superior predictive performance compared to prothrombin time (PT), prothrombin activity (PTA), international normalized ratio (INR), activated partial thromboplastin time (APTT), fibrinogen, and thrombin time (TT). The area under the curve (AUC) was 0.656 (95% CI 0.580-0.733), with a cut-off value of 0.065 mg/l (Figure 2C). The cut-off value for PT was 11.15s; for PTA, it was 130%; for INR, the cut-off value was 1.025; for APTT, it was 35.8s; for fibrinogen, the cut-off value was 3.175g/l; and for TT, it was 14.95s. Subsequently, through univariate and multivariate analyses, we identified the factors affecting LLNM, including extrathyroidal extension, multifocality, and D-dimer>0.065mg/l (Figure 2D).

By comparing gene expression levels in tumor and normal tissue samples, we identified 1374 differentially expressed genes. The heatmap displaying differentially expressed genes is shown in Figure 3A. For deeper analysis, we obtained 377 CRGs from the Genecards website (https://www.genecards.org↗), and 64 of these CRGs were differentially expressed between THCA tissues and normal controls as shown in the Venn plot (Figure 3B). We used univariate COX regression analysis to assess the relationship between the 64 CRGs and PFI, identifying 23 prognostic CRGs (Figure 3C). The correlations of these genes are shown in Figure 3D. Additionally, we constructed a PPI network of the 23 prognostic CRGs (Figure 3E). The Kaplan-Meier survival curves of 23 prognostic CRGs in THCA are shown in Figure 4.

We obtained transcriptomic data and corresponding clinical characteristics of patients from the TCGA-THCA cohort. Based on the 23 prognostic CRGs, we classified THCA patients into two subtypes using the unsupervised clustering method, including coagulation-related cluster 1 and cluster 2 (Figures 5A, B). We conducted K-M survival analysis on the TCGA-THCA cohort, and the results indicated that cluster 1 had a better prognosis than cluster 2 (p=0.014) (Figure 5C).

We divided the patients in the TCGA-THCA dataset into training and validation sets in a 6:4 ratio. In the training cohort of TCGA-THCA, we conducted LASSO regression analysis (Figures 6A,B), screening out 8 key genes (AZU1, COL3A1, CP, CSF2, F12, GNA14, IL1RN, and SERPIND1). Using stepwise multivariate Cox regression analysis, we determined the coefficients associated with each gene (Figure 6C): Riskscore = 0.23947 * AZU1 + 0.00155 * COL3A1 + 0.08232 * CP + 0.03633 * CSF2 + 0.61566 * F12 + 0.09073 * GNA14 + 0.03401 * IL1RN + 0.03270 * SERPIND1.

We classified THCA patients in the training and validation sets into high-risk and low-risk groups based on the median risk score. Figures 6G–I provide detailed information on the risk score distribution and survival status of each THCA patient. Kaplan-Meier survival analysis curves were drawn to compare PFI differences between high- and low-risk groups in the training, validation, and overall sets. The results indicated that PFI was shorter in high-risk patients compared to low-risk patients (all P<0.05, Figures 6D–F), suggesting that the prognostic model has clinical utility in predicting THCA prognosis. Our study found that AZU1, COL3A1, CP, CSF2, F12, GNA14, IL1RN, and SERPIND1 were overexpressed in the high-risk group (Figures 6J–L).

The ROC curves further showcased the excellent performance of the prognostic model in predicting PFI. We found that the AUC for the prognostic features in the training set reached 1.000, 0.901, and 0.924 at 1, 3, and 5 years, respectively. Likewise, in the test set, the AUC were 0.790, 0.777, and 0.828, respectively. For the entire cohort, the AUC reached 0.865, 0.860, and 0.877, respectively (Figures 7A, C, E). Figures 7B, D, F illustrate the ROC curves for the risk score and different clinical characteristics predicting 5-year PFI.

To assess the applicability of the model for patients with various clinical characteristics, we compared survival differences between different risk score groups within each subgroup (Figure 8). We found that the prognostic model exhibited good performance in predicting PFI across subgroups including male, female, N0, N1, Stage II-IV, M0, T1-2, and T3-4 (Figure 8).

A heatmap illustrated the distribution of clinical characteristics and gene expression in the high-risk and low-risk groups (Figure 9B). To assess the independence of the risk score as a prognostic indicator for THCA, univariate and multivariate Cox regression analyses were conducted in the TCGA-THCA cohort. The analyses consistently affirmed the risk score’s status as an independent prognostic predictor (all p < 0.05) (Figure 9A). We created a nomogram that integrates the risk score and clinical characteristics to enhance the clinical application of our results (Figure 9C). Figures 10A, B depict the gene mutation occurrences in the low-risk and high-risk groups, respectively. In the low-risk group, the somatic mutation frequency rankings were BRAF (75%), NRAS (4%), TG (4%), with mutation frequencies of other genes below 4%. In the high-risk group, the somatic mutation frequency rankings were BRAF (43%), NRAS (13%), HRAS (6%), TG (4%), with mutation frequencies of other genes below 4%. It is noteworthy that the primary type of somatic mutation in both high- and low-risk groups was missense mutation, indicating that missense mutations may play a crucial role in the common mechanisms of tumorigenesis across different risk levels.

We examined the biological functions and pathways of the high-risk and low-risk groups using GO and KEGG pathway analyses. GO analysis indicated that the Biological process, cellular component, and Molecular Function affected in the low-risk group mainly included antimicrobial peptide production, extracellular matrix disassembly, lamellar body, neuron to neuron synapse, hippocampal mossy fiber to CA3 synapse, serine-type endopeptidase activity, and serine-type peptidase activity (Figure 10C). KEGG analysis revealed that the significantly enriched pathways in the low-risk group were the Ras signaling pathway, Taste transduction, Pancreatic secretion, and Neuroactive ligand-receptor interaction (Figure 10D). GO analysis indicated that the Biological process, cellular component, and Molecular Function affected in the high-risk group mainly included ossification, osteoblast differentiation, collagen-containing extracellular matrix, collagen trimer, extracellular matrix structural constituent, and extracellular matrix structural constituent conferring tensile strength (Figure 10E). KEGG analysis revealed that the significantly enriched pathways in the high-risk group were Protein digestion and absorption, Cytoskeleton in muscle cells, and ECM-receptor interaction (Figure 10F).

We further analyzed the differences in immune infiltration levels between the high-risk and low-risk groups. The results indicated that the risk score was negatively correlated with the infiltration levels of macrophage M1 and M2 subtypes, CD8+ T cells, and NK cells, with correlation coefficients all below -0.2 (Figure 11A), suggesting that the role of these cells may weaken as the risk score increases. The immune function analysis of different risk groups is depicted in Figure 11B. The differential expression analysis of immune checkpoints between different risk groups showed that immune checkpoints significantly expressed in the low-risk group compared to the high-risk group included CD274, TMIGD2, and TNFRSF18, while immune checkpoints significantly expressed in the high-risk group compared to the low-risk group included CD27, IDO2, and NRP1 (Figure 12A). We also calculated the Tumor Immune Dysfunction and Exclusion (TIDE) score for each group and compared them. The results indicated that the high-risk group had higher TIDE scores, and higher TIDE scores were associated with poorer prognosis. The group with high TIDE scores and high-risk scores had a relatively worse prognosis compared to other groups (Figures 12B–D). TME scores also suggested better prognosis for the low-risk group (Figure 12E).

Drug sensitivity analysis results indicated that drugs such as Cisplatin, Dactinomycin.1, Docetaxel, Erlotinib, and Fludarabine had higher sensitivity in the high-risk group, implying they might have greater potential in the treatment of patients in the high-risk group (Figure 13).

The expression levels of 10 SHMRGs were detected in human PTC cell lines (IHH4, KTC-1, TPC-1) and human normal thyroid cell line (Nthy ori-3-1) using qPCR technology (Figure 14).