Coagulation-related genes for thyroid cancer prognosis, immune infltration, staging, and drug sensitivity

Article type: Research Article

Keywords: thyroid cancer, coagulation, prognosis signature, immune infiltration, clinical relevance, drug sensitivity

Authors: Yuxiao Sun, Yifei Zhang, Yuchuan Yang, Weihao Liu, Detao Yin

Affiliations: Department of Thyroid Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China; Department of Radiotherapy, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China; Engineering Research Center of Multidisciplinary Diagnosis and Treatment of Thyroid Cancer of Henan Province, Zhengzhou, China; Key Medicine Laboratory of Thyroid Cancer of Henan Province, Zhengzhou, China

License: Copyright © 2024 Sun, Zhang, Yang, Liu and Yin CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Article links: DOI: 10.3389/fimmu.2024.1462755  |  PubMed: 39497824  |  PMC: PMC11532168

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (10.2 MB)

Introduction

Thyroid carcinoma (THCA) is the most prevalent endocrine tumor (ref. 1), and the incidence rate of THCA in China has significantly increased since 2000, with an estimated 22,000 new cases reported in 2022 (ref. 2). THCA is classified into three pathological types: differentiated thyroid cancer (DTC), medullary thyroid cancer (MTC), and anaplastic thyroid cancer (ATC) (ref. 3). Despite the generally good prognosis for most thyroid cancer patients, 10%-15% of patients experience recurrence, and 5% develop distant metastases (ref. 4). To further improve the prognosis of THCA patients, it is important to identify suitable prognostic biomarkers.

Coagulation, one of the hallmarks of tumor, could be a consequence of increasing plasma extravasation and vascular permeability which leads to extravascular coagulation, or be activated by disruption of vessels which leads to intravascular coagulation (ref. 5). A study has shown that tumor cells can release procoagulant factors, such as tissue factor, which may trigger the coagulation cascade both in vitro and in vivo (ref. 6). On the other hand, tumor coagulum, a cancer-driven network of molecular effectors favoring bleeding or thrombosis, could interact with the tumor microenvironment (TME) to orchestrate cancer inhibition or progression (ref. 7). Additionally, thyroid hormones directly regulate the transcription of genes encoding coagulation proteins in hepatocytes and endothelial cells (ref. 8). Therefore, the coagulation pathway may interact with the occurrence and development of THCA, but the role of coagulation in THCA patients remains unclear.

In this study, using The Cancer Genome Atlas (TCGA) cohort, we identified coagulation-related molecular subtypes via consensus clustering analysis and compared their prognoses. We used gene expression data and clinical information from TCGA’s THCA to identify differentially expressed Coagulation-related genes (CRGs) in THCA and analyzed their correlation with prognosis. Through LASSO regression, we further filtered out 8 prognostic CRGs to construct a prognostic model and examined the relationships between risk scores and clinicopathological features, molecular functions, pathways, outcomes, immune infiltration, and immunotherapy. Additionally, we assessed the predictive performance of coagulation indicators for LLNM.

Methods

Patient selection and data collection

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 (ref. 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 (ref. 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 (ref. 3).

Data download and organization

The flow chart is shown in f1. 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” (ref. 7, ref. 11).

Collection of differentially expressed CRGs

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.

Identification of coagulation subtypes and survival analysis

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.

Construction of prognostic model

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 (ref. 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.

Validation of prognostic signature

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” (ref. 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.

Enrichment analysis and somatic mutations

To explore the differences in somatic mutations between different risk groups, we performed analysis and visualization using the R package “maftools” (ref. 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) (ref. 15–ref. 18).

Immune microenvironment and immunotherapy

Furthermore, we employed the EPIC, MCP-COUNTER, TIMER, XCELL, QUANTISEQ, CIBERSORT, and CIBERSORT-ABS algorithms to examine immune cell infiltration (ref. 19–ref. 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 (ref. 26, ref. 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 (ref. 28). Finally, we used ESTIMATE to compare TME scores (StromalScore, ImmuneScore, and EstimateScore) across different risk groups.

Development of individualized anti-tumor treatment protocols

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 (ref. 29).

Reverse transcription−qPCR (RT−qPCR)

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′.

Statistical analysis

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.

Results

Predictive performance of coagulation indexes

We collected 408 cases of adult PTC, and their clinicopathological characteristics are shown in f2. The mean of coagulation profiles between LLNM and non-LLNM patients was compared in f2. 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 (f2). 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 (f2).

Dysregulated CRGs in THCA and their prognosis

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 f3. 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 (f3). We used univariate COX regression analysis to assess the relationship between the 64 CRGs and PFI, identifying 23 prognostic CRGs (f3). The correlations of these genes are shown in f3. Additionally, we constructed a PPI network of the 23 prognostic CRGs (f3). The Kaplan-Meier survival curves of 23 prognostic CRGs in THCA are shown in f4.

Identification of coagulation-related subtypes

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 (f5). 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) (f5).

Prognostic model construction and validation

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 (f6), 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 (f6): 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. f6 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, f6), 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 (f6).

Validation of the 8-Gene signature

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 (f7). f7 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 (f8). 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 (f8).

A heatmap illustrated the distribution of clinical characteristics and gene expression in the high-risk and low-risk groups (f9). 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) (f9). We created a nomogram that integrates the risk score and clinical characteristics to enhance the clinical application of our results (f9). f10 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.

Functional, immune infiltration analyses and drug sensitivity analysis in different risk groups

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 (f10). 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 (f10). 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 (f10). 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 (f10).

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 (f11), 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 f11. 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 (f12). 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 (f12). TME scores also suggested better prognosis for the low-risk group (f12).

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 (f13).

Validation of the genes expression

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 (f14).

Discussion

THCA accounts for 3% of all cancers. While it is mostly regarded as an indolent disease, its recurrence rate remains at 10%-15% (ref. 4). Generally, cancer is associated with an imbalance in the hemostatic system, but the pathogenesis of cancer-related coagulation disorders is complex (ref. 30). However, there is limited research on the relationship between THCA and coagulation. This study aims to clarify the correlation between coagulation and prognosis in THCA patients.

Several studies have indicated that PTC patients with LLNM have a higher risk of recurrence compared to those without LLNM (ref. 31, ref. 32). Generally, prophylactic central lymph node dissection is deemed necessary. However, prophylactic dissection of negative lateral lymph nodes is not recommended because it may lead to numerous surgical complications. Therefore, it is important to determine the presence of LLNM preoperatively. This study analyzed the influence of coagulation markers (PT, PTA, INR, APTT, fibrinogen, TT, and D-dimer) on LLNM. The ROC curve showed that D-dimer has a good predictive performance for LLNM. Multivariate analysis indicated that extrathyroidal extension, multifocality, and D-dimer >0.065mg/l are independent predictors of LLNM. The effects of extrathyroidal extension and multifocality on LLNM have been confirmed in numerous studies (ref. 33, ref. 34). The research has shown that elevated D-dimer is significantly associated with cancer recurrence, metastasis, and worse survival outcomes (ref. 35). This study suggests that D-dimer >0.065mg/l is significantly associated with LLNM in PTC.

In this study, we divided THCA patients into two distinct subtypes based on the expression of prognostically valuable CRGs. K-M survival analysis indicated that cluster1 had better outcomes than cluster2 (p=0.014), offering a reference for refining the classification and management of THCA patients in clinical practice. Subsequently, we used LASSO regression to screen 8 genes (AZU1, COL3A1, CP, CSF2, F12, GNA14, IL1RN, and SERPIND1) to construct the prognostic model. ROC results from both the training and validation sets suggested that the constructed model had high predictive efficacy, and the high-risk group had poorer outcomes. We further used univariate and multivariate Cox regression analyses to validate the effectiveness of the prognostic model. We also developed a nomogram integrating risk scores and clinical features to facilitate the clinical application of our research findings.

The tumor immune microenvironment (TME) plays an important role in tumor development. Recently, a pan-cancer analysis of the human tumor coagulome revealed that coagulation links to the TME (ref. 36). In our study results, the high-risk group had lower infiltration levels of macrophage M1 and M2 subtypes, CD8+ T cells, and NK cells, which may be related to poorer prognosis. Macrophages may combat cancer progression by directly engulfing cancer cells or activating antitumor immune responses, and CD8+ T cells and NK cells also have strong antitumor functions (ref. 37). To explore the biological functional differences between high and low-risk groups, we conducted GO and KEGG analyses. The results indicated significant differences in biological functions between the two groups, but the underlying mechanisms remain to be clarified. Finally, we investigated the intrinsic connection between drug sensitivity and risk groups by correlating THCA patients’ prognostic risk scores with the half-maximal inhibitory concentration (IC50) values of chemotherapeutic drugs.

Research has indicated that these genes are involved in cancer development. In gastric cancer, AZU1 was upregulated (ref. 38). COL3A1 could be an oncogene and promote drug resistance in lung cancer (ref. 39). The product expressed by CP, cancer procoagulant, is a biomarker for cancer. Elevated CP activity has been detected in pancreatic, breast, lung, digestive system, and urinary system cancers (ref. 40). In breast cancer, CSF2 could activate the Stat3 pathway in CAA via paracrine or autocrine mechanisms, leading to increased expression and secretion of CXCL3. CXCL3 binds to CXCR2, a receptor on breast cancer cells, and activates FAK, which promotes a mesenchymal phenotype, invasion, and metastasis of breast cancer cells (ref. 41) F12 expression might affect the OS of PTC patients by regulating metabolic pathways (ref. 42). GNA14 might accelerate colorectal cancer cell proliferation and malignant tumor progression through ERK and β-catenin pathways (ref. 43). IL1RN was a good prognostic and diagnostic biomarker for PTC and might promote thyroid cancer progression through immune-related pathways (ref. 44). SERPIND1 promoted the proliferation, migration, invasion, G1-to-S phase transition, and epithelial-mesenchymal transition of ovarian cancer cells and inhibited their apoptosis by promoting phosphorylation in the phosphoinositide 3-kinase/protein kinase B (PI3K/AKT) pathway (ref. 36). We also validated the high expression of these 8 genes in PTC by qRT-PCR.

There are some limitations in this study. Firstly, the TCGA-THCA cohort has a limited number of cases, and more datasets are needed to validate these findings. Otherwise, the selected genes were only validated by qRT-PCR in PTC, lacking comprehensive experimental validation. Moreover, although we found that the coagulation pathway affects the prognosis and immune microenvironment of THCA patients, its underlying mechanisms need further investigation.

Conclusion

Our study demonstrated that coagulation is related to immune infiltration and prognosis in THCA. Our study suggested the possibility of D-dimer predicting LLNM. Subsequently, we used the TCGA cohort to construct a new coagulation-related THCA risk scoring model for prognosis assessment and risk stratification in THCA patients. This risk model could provide a robust prognostic tool and promote clinical guidance for THCA patients. We also constructed a nomogram combining the model and clinical features to predict the 3-, 5-, and 7-year survival probabilities of patients. We analyzed and compared high- and low-risk groups regarding immune infiltration, somatic mutations, and other aspects, offering some guidance for treatment strategy selection. In summary, our research will aid in understanding the role and significance of the coagulation in THCA.

References

1. S Vaccarella, L Dal Maso, M Laversanne, F Bray, M Plummer, S Franceschi. The impact of diagnostic changes on the rise in thyroid cancer incidence: A population-based study in selected high-resource countries.. Thyroid. (, 2015. [DOI]
2. C Xia, X Dong, H Li, M Cao, D Sun, S He. Cancer statistics in China and United States, 2022: profiles, trends, and determinants.. Chin Med J (Engl). (, 2022. [DOI]
3. BR Haugen, EK Alexander, KC Bible, GM Doherty, SJ Mandel, YE Nikiforov. 2015 American thyroid association management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer: the american thyroid association guidelines task force on thyroid nodules and differentiated thyroid cancer.. Thyroid. (, 2016. [DOI | PubMed]
4. D Laha, N Nilubol, M Boufraqech. New therapies for advanced thyroid cancer.. Front Endocrinol (Lausanne). (, 2020. [DOI | PubMed]
5. B Furie, BC Furie. Mechanisms of thrombus formation.. N Engl J Med. (, 2008. [DOI]
6. B Kocatürk, HH Versteeg. Tissue factor isoforms in cancer and coagulation: may the best isoform win.. Thromb Res. (, 2012. [DOI]
7. J Yang, C Wang, Y Zhang, S Cheng, M Wu, S Gu. Clinical significance and immune infiltration analyses of a novel coagulation-related signature in ovarian cancer.. Cancer Cell Int. (, 2023. [DOI | PubMed]
8. A Bano, L Chaker, MPM de Maat, F Atiq, M Kavousi, OH Franco. Thyroid function and cardiovascular disease: the mediating role of coagulation factors.. J Clin Endocrinol Metab. (, 2019. [DOI]
9. MP Vanderpump, WM Tunbridge, JM French, D Appleton, D Bates, F Clark. The incidence of thyroid disorders in the community: a twenty-year follow-up of the Whickham Survey.. Clin Endocrinol (Oxf). (, 1995. [DOI | PubMed]
10. Y Sun, Y Liu, H Li, Y Tang, W Liu, Y Zhang. The significance and prognostic value of multifocal papillary thyroid carcinoma in children and adolescents.. BMC Cancer. (, 2024. [DOI | PubMed]
11. G Stelzer, N Rosen, I Plaschkes, S Zimmerman, M Twik, S Fishilevich. The geneCards suite: from gene data mining to disease genome sequence analyses.. Curr Protoc Bioinf. (, 2016. [DOI]
12. J Friedman, T Hastie, R Tibshirani. Regularization paths for generalized linear models via coordinate descent.. J Stat Software. (, 2010. [DOI]
13. P Blanche, JF Dartigues, H Jacqmin-Gadda. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks.. Stat Med. (, 2013. [DOI]
14. A Mayakonda, DC Lin, Y Assenov, C Plass, HP Koeffler. Maftools: efficient and comprehensive analysis of somatic variants in cancer.. Genome Res. (, 2018. [DOI]
15. ME Ritchie, B Phipson, D Wu, Y Hu, CW Law, W Shi. limma powers differential expression analyses for RNA-sequencing and microarray studies.. Nucleic Acids Res. (, 2015. [DOI | PubMed]
16. T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.. Innovation (Camb). (, 2021. [DOI | PubMed]
17. Z Gu, L Gu, R Eils, M Schlesner, B Brors. circlize Implements and enhances circular visualization in R.. Bioinformatics. (, 2014. [DOI]
18. Z Gu, R Eils, M Schlesner. Complex heatmaps reveal patterns and correlations in multidimensional genomic data.. Bioinformatics. (, 2016. [DOI]
19. AM Newman, CL Liu, MR Green, AJ Gentles, W Feng, Y Xu. Robust enumeration of cell subsets from tissue expression profiles.. Nat Methods. (, 2015. [DOI]
20. T Li, J Fan, B Wang, N Traugh, Q Chen, JS Liu. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells.. Cancer Res. (, 2017. [DOI]
21. R Dienstmann, G Villacampa, A Sveen, MJ Mason, D Niedzwiecki, A Nesbakken. Relative contribution of clinicopathological variables, genomic markers, transcriptomic subtyping and microenvironment features for outcome prediction in stage II/III colorectal cancer.. Ann Oncol. (, 2019. [DOI]
22. D Aran, Z Hu, AJ Butte. xCell: digitally portraying the tissue cellular heterogeneity landscape.. Genome Biol. (, 2017. [DOI | PubMed]
23. F Finotello, C Mayer, C Plattner, G Laschober, D Rieder, H Hackl. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data.. Genome Med. (, 2019. [DOI | PubMed]
24. E Schenk, J Boland, A Mansfield, MC Aubry, A Dietz. Local and systemic immunity predict survival in patients with pulmonary sarcomatoid carcinoma.. Med Oncol. (, 2017. [DOI | PubMed]
25. J Racle, K de Jonge, P Baumgaertner, DE Speiser, D Gfeller. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data.. Elife. (, 2017. [DOI]
26. S Hänzelmann, R Castelo, J Guinney. GSVA: gene set variation analysis for microarray and RNA-seq data.. BMC Bioinf. (, 2013. [DOI]
27. Z Zhang. Reshaping and aggregating data: an introduction to reshape package.. Ann Transl Med. (, 2016. [DOI | PubMed]
28. J Fu, K Li, W Zhang, C Wan, J Zhang, P Jiang. Large-scale public data reuse to model immunotherapy response and resistance.. Genome Med. (, 2020. [DOI | PubMed]
29. D Maeser, RF Gruener, RS Huang. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data.. Brief Bioinform. (, 2021. [DOI]
30. A Ordookhani, A Motazedi, KD Burman. Thrombosis in thyroid cancer.. Int J Endocrinol Metab. (, 2018. [DOI | PubMed]
31. SK Kim, I Park, N Hur, M Rayzah, JH Lee, JH Choe. Patterns, predictive factors, and prognostic impact of contralateral lateral lymph node metastasis in N1b papillary thyroid carcinoma.. Ann Surg Oncol. (, 2017. [DOI]
32. N Chéreau, C Buffet, C Trésallet, F Tissier, L Leenhardt, F Menegaux. Recurrence of papillary thyroid carcinoma with lateral cervical node metastases: Predictive factors and operative management.. Surgery. (, 2016. [DOI]
33. YK So, MJ Kim, S Kim, YI Son. Lateral lymph node metastasis in papillary thyroid carcinoma: A systematic review and meta-analysis for prevalence, risk factors, and location.. Int J Surg. (, 2018. [DOI | PubMed]
34. HJ Kim, SY Sohn, HW Jang, SW Kim, JH Chung. Multifocality, but not bilaterality, is a predictor of disease recurrence/persistence of papillary thyroid carcinoma.. World J Surg. (, 2013. [DOI]
35. W Kong, L Zhang, R An, M Yang, H Wang. Diagnostic value of serum D-dimer for detection of gallbladder carcinoma.. Cancer Manag Res. (, 2021. [DOI]
36. Z Saidak, S Soudet, M Lottin, V Salle, MA Sevestre, F Clatot. A pan-cancer analysis of the human tumor coagulome and its link to the tumor immune microenvironment.. Cancer Immunol Immunother. (, 2021. [DOI]
37. KE de Visser, JA Joyce. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth.. Cancer Cell. (, 2023. [DOI | PubMed]
38. X Ran, X Xu, Y Yang, S She, M Yang, S Li. A quantitative proteomics study on olfactomedin 4 in the development of gastric cancer.. Int J Oncol. (, 2015. [DOI]
39. L Wang, Y Sun, Z Guo, H Liu. COL3A1 overexpression associates with poor prognosis and cisplatin resistance in lung cancer.. Balkan Med J. (, 2022. [DOI | PubMed]
40. SD Szajda, B Darewicz, J Kudelski, T Werel, K Zwierz, Z Skrzydlewski. Assessment of the efficacy of cancer procoagulant (CP) activity evaluation in the oncological diagnosis of the urinary system.. Pol Merkur Lekarski. (, 2005
41. X He, L Wang, H Li, Y Liu, C Tong, C Xie. CSF2 upregulates CXCL3 expression in adipocytes to promote metastasis of breast cancer via the FAK signaling pathway.. J Mol Cell Biol. (, 2023. [DOI]
42. JH Luo, XX Zhang, WH Sun. F12 as a reliable diagnostic and prognostic biomarker associated with immune infiltration in papillary thyroid cancer.. Aging (Albany NY). (, 2022. [DOI]
43. R Park, S Lee, H Chin, AT Nguyen, D Lee. Tumor-promoting role of GNA14 in colon cancer development.. Cancers (Basel). (, 2023. [DOI]
44. Z Xie, X Li, Y He, S Wu, S Wang, J Sun. Analysis of the expression and potential molecular mechanism of interleukin-1 receptor antagonist (IL1RN) in papillary thyroid cancer via bioinformatics methods.. BMC Cancer. (, 2020. [DOI | PubMed]