Development of a Risk Predictive Model for Evaluating Immune Infiltration Status in Invasive Thyroid Carcinoma

Aims This study aimed to reveal the molecular characteristics and potential biomarker of immune-activated and immunosuppressive invasive thyroid carcinoma. Methods Expression and clinical data for invasive thyroid carcinoma were obtained from the TCGA database. Tumor samples were divided into immune-activated or immunosuppressive groups based on the immune enrichment score calculated by ssGSEA. Differentially expressed genes (DEGs) between tumor vs. normal groups or between immune-activated vs. immunosuppressive groups were screened, followed by functional enrichment. Immune infiltration was evaluated using the ESTIMATE, CIBERSORTx, and EPIC algorithms, respectively. A random forest algorithm and Lasso cox analysis were used to identify gene signatures for risk model construction. Results Totally 1171 DEGs were screened between tumor vs. normal groups, and multiple tumorigenesis-associated pathways were significantly activated in invasive thyroid carcinoma. Compared to immune-activated samples, immunosuppressive samples showed higher tumor purity, lower immune/stromal scores, and lower expression of immune markers, as well as lower infiltration abundance of CD4+ T cells and CD8+ T cells. A risk model based on a 12-immune signature (CCR7, CD1B, CD86, CSF2RB, HCK, HLA-DQA1, LTA, LTB, LYZ, NOD2, TNFRSF9, and TNFSF11) was developed to evaluate the immune infiltration status (AUC = 0.998; AUC of 0.958 and 0.979 in the two external validation datasets), which showed a higher clinical benefit and high accuracy. Immune-activated samples presented lower IC50 value for bortezomib, MG.132, staurosporine, and AZD8055, indicating sensitivity to these drugs. Conclusion A 12-gene-based immune signature was developed to predict the immune infiltration status for invasive thyroid carcinoma patients and then to identify the subsets of invasive thyroid carcinoma patients who might benefit from immunotherapy.


Introduction
yroid carcinoma is the most common malignant tumor of the endocrine system, originating in the thyroid follicular epithelium or parafollicular epithelium [1]. Global Cancer Statistics indicate that there are an estimated 586,202 new cases of thyroid carcinoma in 2020, accounting for 3.0% of the global cancer burden [2]. In China, thyroid carcinoma ranks the fourth most common type of cancer in females and the ninth in males, accounting for 38% of global new cases of thyroid carcinoma [3]. yroid carcinoma can be divided into four histological types based on their tumor origin and differentiation, including papillary, follicular, medullary, and anaplastic thyroid carcinoma [4]. Of which, papillary thyroid carcinoma (PTC) is the most frequent type, accounting for 85% to 90% of thyroid carcinoma. Due to relatively inert clinical biological behavior, the prognosis of PTC is generally excellent, with a 10-year-survival rate of over 90% [5]. However, about 5%-34% patients with differentiated thyroid carcinoma present extrathyroidal extension, which negatively affects the prognosis [6][7][8]. Extrathyroidal lesions may implicate crucial structures in the central neck, and gross residual lesions postoperatively may lead to recurrence [9]. Hence, it is necessary to propose effective approaches to improve the prognosis of invasive thyroid carcinoma.
Reportedly, a risk predictive model contributes to make corresponding decisions for treatment and management of patients with thyroid carcinoma by dividing patients into different risk groups [10,11]. Several risk models have been established for thyroid carcinoma [12][13][14]. For example, Shen et al. established a six-genotype genetic molecular (involves mutations in RAS, BRAF V600E, and TERT promoter) prognostic system for PTC, which helps to identify the subsets of PTC patients with the highest risk of invasion for personalized and precise treatment [15]. Wang et al. constructed a risk predictive model based on the N 6methyladenosine-related signature to predict the diseasefree survival for PTC patients [14]. However, few studies focus on immune-related genes signature in the risk prediction ofinvasive thyroid carcinoma.
Tumor cells escaping immune surveillance and destruction present a hallmark of tumors. Immune cells in the tumor microenvironment affect the development and progression of tumors by secreting soluble cytokines to mediate the proliferation of tumor cells [16]. Immunotherapy has shown survival benefits by unleashing the immune system and activating cytotoxic lymphocytes to kill cancer cells [17]. However, only part of patients benefit from immunotherapy. erefore, this study proposed a risk model based on immune-related signature to evaluate the immune activation and immunosuppression status for invasive thyroid carcinoma patients. is will help to identify the subsets of invasive thyroid carcinoma patients who might benefit from immunotherapy or not, so as to get personalized and precise treatment.

Data Acquisition and Preprocessing.
Data of genes expression, clinical phenotype, mutation, and survival for thyroid carcinoma in e Cancer Genome Atlas (TCGA) database were downloaded from UCSC Xena (https:// xenabrowser.net/datapages/). From which, invasive thyroid carcinoma samples (extrathyroid carcinoma present extension status ≥ T3) were selected for analysis. Totally 176 samples with both expression data and clinical phenotype were included, containing 154 tumor samples and 22 normal samples. Two microarray datasets, GSE35570, and GSE65074 were downloaded from the Gene Expression Omnibus (GEO) database. ese two microarray datasets were used to validate the risk-score model.

Differential Expression Analysis.
Principal component analysis (PCA) dimension reduction was performed for the expression data. en, differential expression analysis was performed using the Limma package in R to screen the differentially expressed genes (DEGs) between tumor vs. normal groups or between immune activation vs. immunosuppression groups, respectively. e DEGs were screened with the cut-off value of |logFC| > 2 and adjusted P < 0.01.

Functional Enrichment
Analysis. Enrichment analysis for gene ontology (GO) annotation terms and KEGG pathways was conducted using the clusterProfiler package [18] in R to uncover the involved functions of DEGs. e significantly enriched GO terms and pathways were selected with P < 0.05. Gene set enrichment analysis (GSEA) was also performed using the clusterProfiler package in R to reveal the differences in terms of pathways between groups.

Immune-Activated and Immunosuppressive Tumors.
Based on the 29 immune-related gene sets (Table S1), singlesample GSEA (ssGSEA) was conducted using the GSVA package [19] to calculate the enrichment score for each tumor sample. en, the 154 tumor samples were assigned into immune-activated and immunosuppressive groups using sparse hierarchical clustering. After PCA dimension reduction, the DEGs between immune activation vs. immunosuppression groups were screened with the cut-off value of |logFC| > 1.5 and adjusted P < 0.01.

Evaluation of Immune Infiltration.
e stromal-score, immune-score, and tumor-purity of each tumor sample were evaluated using the ESTIMATE package in R. Composition and infiltration abundance of immune cells of each tumor sample were evaluated using CIBERSORTx and EPIC methods provided in TIMER 2.0 database [20].

Feature Genes Screening Using Machine Learning.
Immune-related genes (IRGs , Table S2) were obtained from the ImmPort (https://immport.niaid.nih.gov) database [21]. en, IRGs were merged with DEGs that screened between immune activation vs. immunosuppression groups, and the overlapped genes were considered as DEIRGs. DEIRGs were then analyzed by a random forest algorithm provided in the Boruta package [22] to screen feature genes. e feature genes were then uploaded to the STRING database to reveal the interactions among their encoded proteins.

Construction of Risk-Score Model.
e optimal gene signature was further identified from the feature genes using the Lasso cox analysis provided in the glmSparseNet package. e gene signature identified by Lasso cox analysis was used to construct a risk-score model. A receiver operator characteristic (ROC) curve was plotted using the pROC package [23] to evaluate the predictive performance of the risk-score model. Decision curve analysis (DCA) and the clinical impact curve were plotted using the DecisionCurve package (version 1.3) to evaluate the clinical benefit. e GSE35570 and GSE65074 datasets were used as external datasets for validation. e invasive tumor samples in these two datasets were divided into immune-activated and immunosuppressive groups based on the methods mentioned above. Based on the identified gene signature, a risk-score model was established based on these two datasets, followed by evaluating the predictive performance using the ROC curve, DCA, and clinical impact curve.

Mutation
Analysis and Drug Sensitivity Analysis. In order to investigate the possible causes of differential expression for signature genes, mutation analysis was conducted using the maftools package [  Evidence-Based Complementary and Alternative Medicine 3 SNP data in the TCGA database. Based on the Genomics of Drug Sensitivity in Cancer (GDSC) database, the drug sensitivity of each sample was evaluated by calculating the IC50 value using the OncoPredict package (version 0.2) to investigate the degree of patient response to drugs.

Statistical
Analysis. e differences in stromal-score, immune-score, tumor-purity, immune markers (including T cell stimulators, T cell inhibitors, and histocompatibility complexes) expression, gene mutation, and drug sensitivity between immune-activated and immunosuppressive groups were compared using the T-test. P < 0.05 shows the statistical significance.

Genes Differentially Expressed in Invasive yroid
Carcinoma Compared to Normal Control. In order to investigate the differences between invasive thyroid carcinoma and normal control at the molecular level, we screened the DEGs between these two groups first. PCA indicated that samples from the same group were clustered together, indicating that gene expression in these two groups was significantly distinguished (Figure 1(a)). Totally 1171 DEGs were screened, of which 471 genes were upregulated while 700 genes were downregulated in invasive thyroid carcinoma compared to normal control (Figure 1(b)). e bidirectional clustering heatmap confirmed that these DEGs could clearly distinguish the samples into two groups (Figure 1(c)). Enrichment analysis was performed to investigate the involved functions of these DEGs. e results showed that DEGs were enriched in GO annotation terms, including cell junction assembly, synapse assembly, collagen-containing extracellular matrix, extracellular matrix structural constituent, and glycosaminoglycan binding (Figures S1A-C), and in KEGG pathways, such as neuroactive ligand-receptor interaction, ECM-receptor interaction, and cell adhesion molecules ( Figure S1D). e DEGs might contribute to tumor invasion by these biological functions and pathways.
GSEA analysis revealed significant differences in terms of pathways between tumor and normal groups. As shown in Figures 1(d) and 1(e), various tumorigenesis-associated pathways were significantly activated in invasive thyroid carcinoma samples, including p53 signaling pathway, cell adhesion molecules, cytokine-cytokine receptor interaction, microRNAs in cancer, and transcriptional misregulation in cancer. Multiple metabolism pathways, such as propanoate metabolism, butanoate metabolism, and thyroid hormone synthesis, were markedly inhibited in invasive thyroid carcinoma samples.

Immune-Activated and Immunosuppressive Invasive
yroid Carcinoma. e immune enrichment score of each tumor sample was calculated using ssGSEA, and then the tumor samples were assigned into immune-activated and immunosuppressive groups using sparse hierarchical clustering (Figure 2(a)). Compared to immune-activated samples, immunosuppressive samples showed lower immune, stromal scores, and expression levels of immune-related gene sets, while showing higher tumor purity (Figure 2(b)). Additionally,   Evidence-Based Complementary and Alternative Medicine expression levels of T cell stimulators, T cell inhibitors, and histocompatibility complexes were all decreased compared to immune-activated samples (Figures 2(c)-2(e)). ese results suggested that immune-activated patients might be more likely to benefit from immunotherapy.

Genes Differentially Expressed between
Immune-Activated and Immunosuppressive Groups. Immune-activated and immunosuppressive samples were distinguished in PCA, indicating the practicability and successfulness of such groups (Figure 3(a)). In total, 765 DEGs were screened between immune-activated and immunosuppressive groups, and the majority of genes (724 genes) were upregulated in immune-activated groups (Figures 3(b) and 3(c)). Functional enrichment indicated that these DEGs were involved in multiple immune related biological function, including T cell activation, leukocyte cell-cell adhesion, lymphocyte   Evidence-Based Complementary and Alternative Medicine proliferation, MHC protein complex, cytokine activity, and immune receptor activity. (Figures 3(d)-3(f )). ese DEGs also significantly enriched in cytokine-cytokine receptor interaction, chemokine signaling pathway, and cell adhesion molecules pathways (Figure 3(g)). GSEA analysis indicated that cell adhesion molecules, cytokine-cytokine receptor interaction, 17 cell differentiation, and chemokine signaling pathways were markedly activated in immune-activated groups ( Figure S2). ese findings suggested that immune-related genes and pathways were markedly upregulated in immune-activated groups.

Immune Infiltration between Immune-Activated and
Immunosuppressive Groups. Immune cells composition for all tumor samples was evaluated using CIBERSORTx. Resting memory CD4+ T cells, CD8+ T cells, naive B cells, and M0/M2 macrophages accounted for high composition (Figure 4(a)). Infiltration abundance of follicular helper T cells, resting NK cells, resting myeloid dendritic cells, and M1 macrophages in immune-activated groups significantly differed from immunosuppressive groups (Figure 4(b)). Based on the EPIC tool, we found CD4+ T cells, CD8+ T cells, and endothelial cells showed higher composition in tumor samples (Figure 4(c)). Of which, immune-activated samples showed higher infiltration abundance of CD4+ T cells and CD8+ T cells compared to that of immunosuppressive samples (Figure 4(d)), indicating high antitumor immunity in immune-activated samples. (DEIRGs). From the ImmPort database, 1793 IRGs were obtained. Of which, 183 DEIRGs were screened by filtrating with the 765 DEGs between immune-activated and immunosuppressive groups (Figure 4(e)). ese DEIRGs were significantly enriched in various biological processes associated with T cells (Figure 4(f )), including T cell activation/ proliferation/differentiation, T cell mediated immunity, and T cell mediated cytotoxicity. Various macrophage-associated biological processes were also enriched ( Figure 4(g)), including macrophage activation/chemotaxis/migration/differentiation, macrophage activation involved in immune response, and macrophage cytokine production. is was consistent with the high infiltration composition of T cells and macrophages in the tumor microenvironment.

Machine Learning and Model Construction.
Using a random forest algorithm provided in Boruta, we screened 81 feature genes from the 183 DEIRGs (Figures 5(a) and 5(b)). e protein-protein interactions were further predicted. Of which, PTPRC, CSF2, CD86, CCR7, and CTLA4 were the top five hub genes with higher degree (Figure 5(c)). Lasso cox analysis identified 12 optimal genes signature ( Figure 5(d)), and a risk-score model was constructed with the following formula: risk score � 47.    Evidence-Based Complementary and Alternative Medicine activated samples than immunosuppressive samples ( Figure 5(f )). DCA curve indicated that the risk score model showed higher clinical benefit ( Figure 5(g)). e clinical impact curve confirmed that the prediction result of riskscore model was close to clinical reality, suggesting a high accuracy ( Figure 5(h)). ese findings suggested that the risk-score model had a good performance to evaluate the immune-activated or immunosuppressive status for invasive thyroid carcinoma.

Validation of Risk-Score Model by External Datasets.
To avoid overfitting of the model in the TCGA-cohort, two external GEO datasets (GSE35570 and GSE65074) were used to validate the model. First, the tumor samples in these two datasets were clustered into immune-activated or immunosuppressive groups, respectively (Figure 6(a)). Similarly, immune-activated samples showed a high expression of immune-related gene sets and a high immune score, but it showed a low tumor purity (Figures 6(b) and 6(c)). Riskscore models were constructed as the formula mentioned above. As shown in Figure 6(d), the model had a good performance to evaluate the risk of immune-activated and immunosuppressive status, with an AUC of 0.958 in GSE35570 dataset and AUC of 0.979 in GSE65074 dataset. Consistently, all the 12 genes in model were highly expressed in immune-activated groups (Figures 6(e) and 6(f )). Furthermore, DCA curves ( Figure 6(g)) and clinical impact curves (Figure 6(h)) showed that the model showed higher clinical benefit and a high accuracy to evaluate the risk of immune-activated and immunosuppressive status.

Drug Sensitivity Analysis.
In order to screen possible drugs for the treatment of invasive thyroid carcinoma, the drugs with IC50 < 1 for all tumor samples were analyzed (Figure 7(a)). From which, four drugs, including bortezomib, MG.132, staurosporine, and AZD8055, showed significant differences between immune-activated and immunosuppressive groups (Figure 7(b)). All these four drugs showed lower IC50 values in immune-activated samples, suggesting that immune-activated patients were more sensitive to these drugs. Immunosuppressive group showed significant high-risk score than immune-activated group (Figure 7(c)). is also confirmed that immune-activated patients might be more likely to benefit from drug therapy.

Mutation Analysis and Correlation Analysis.
Since all the genes in the model were highly expressed in the immuneactivated group compared to immunosuppressive group, we analyzed the possible reasons for the gene expression differences from the perspective of gene mutation. ere were no differences on tumor mutation burden between immuneactivated and immunosuppressive groups (Figure 7(d)). Missense mutation accounted for the majority of the variants (Figure 7(e)). Among the top 20 mutation genes, mutation in BRAF gene accounted for 78% of tumor samples (Figure 7(f )), suggesting that BRAF mutation might be correlated with the development of invasive thyroid carcinoma. Genes in the model were not found in the top 20 mutation genes, suggesting that the expression differences of these genes might not be caused by gene mutation. We  OLR1  CCR5  HLA-DRB1  MARCO  CCL20  IL1RL1  CXCL13  HLA-G  VAV1  TNFSF18  ITGAL  IL32  TUBB3  PNOC  SH2D1A  C3  DUOX2  TNFRSF18  HLA-DRA  CD3E  HLA-DPB1  IL2RA  SEMA3C  HLA-DOA  GREM1  TNFSF13B  MMP12  IL21R  CD3D  CXCR2  S100A8  TNFSF11  PTGER2  CCL23  CCR8  S100A9  HLA-DPA1  FLT3  IL12RB1  PTPRC  IL7R  CCL17  HLA-DQA2  IL22RA2  TLR8  FCGR2B  CCL13  RNASE2  CSF2  COLEC12  CD48  SECTM1  BTK  CD28  ITGB2  HLA-DOB  EBI3  CCL19  TYROBP  CSF2RA  CYBB  CD1A  ICOS  LYZ  CD1C  CD1E  S100B  CTLA4  HLA-DQA1  CSF2RB  CCL22  CD86  CCR7  LTA  IL2RG  NOD2  CD1B  TNFRSF9    Evidence-Based Complementary and Alternative Medicine further analyzed the correlations between genes in the model and tumor-infiltrating immune cells (Figure 7(g)), and we found that most of the genes were significantly correlated with the follicular helper T cells, activated NK cells, M1 macrophages, resting dendritic cells, and resting mast cells.

Discussion
Extrathyroidal extension of thyroid carcinoma has been demonstrated to show an adverse impact on the prognosis of patients [25]. In this study, we first investigated the genes and pathways that dysregulated in invasive thyroid carcinoma based on the data in the TCGA database (thyroid carcinoma samples with extrathyroid carcinoma present extension status ≥ T3). A large-scale of gene expression pattern changed in invasive thyroid carcinoma compared to control. Additionally, various tumorigenesis-associated pathways were significantly activated, including the p53 signaling pathway, cell adhesion molecules, cytokine-cytokine receptor interaction, and transcriptional misregulation in cancer. erefore, we suggested that these dysregulated genes might mediate the development and progression of invasive thyroid carcinoma by affecting tumorigenesis-associated pathways. Notably, the thyroid hormone synthesis pathway was found to be inhibited in invasive thyroid carcinoma samples. yroid stimulating hormone (TSH) is the most important hormone regulating thyroid gland function. It activates the cAMP pathway and regulates hormone synthesis and proliferation of thyroid follicular cells by binding TSH receptor on the membrane of thyroid follicular cells, thus affecting the onset or progression of follicular cells-originated thyroid cancer [26,27]. Based on the immune enrichment score, all the invasive tumor samples were divided into immune-activated and immunosuppressive groups. Immune-activated samples showed lower tumor purity, but it showed higher immune/ stromal scores and higher expression of T cell stimulators, T cell inhibitors, and histocompatibility complexes.
Moreover, immune-activated samples showed higher infiltration abundance of CD4+ T cells and CD8+ T cells. Increased studies have indicated that patients with high immune infiltration (hot tumors) were more likely to benefit from immunotherapy, while immunosuppressive (cold tumor) tumors were prone to be resistant to immune checkpoint blocker [28,29]. Immune infiltration status and levels in TME showed significant correlations with the prognosis of tumor patients [30]. Tumor infiltrating lymphocytes (TIL) is a heterogeneous lymphocyte group with CD8+ T cells as the main effector cells, exert antitumor roles in TME [31,32]. e loss of TILs function in TME mainly responsible for the tumor progression and failure of cellular immunotherapy [33,34]. CD4+ T cells have the dual function both in killing tumor cells directly and in maintaining and promoting CD8+ T cells survival [35]. Considering their importance in antitumor immunity, we concluded that immune-activated patients might be likely to benefit from immunotherapy. In addition, genes that are implicated in immune-related biological processes (such as T cell activation, leukocyte cell-cell adhesion, lymphocyte proliferation, MHC protein complex, cytokine activity, and immune receptor activity) were significantly upregulated in immune-activated patients, indicating that the immunological grouping is reasonable.
Based on an immune-related gene signature (CCR7, CD1B, CD86, CSF2RB, HCK, HLA-DQA1, LTA, LTB, LYZ, NOD2, TNFRSF9, and TNFSF11), a risk-score model was established to evaluate the immune-activated or immunosuppressive status for invasive thyroid carcinoma. e CCR7 chemokine axis not only mediates the trafficking of effector cells that produce an immune response to growing tumors (helping to combat the spread of cancer) but also controls the migration and metastasis of tumor cells to the lymphatic system (contributing to tumor expansion), and therefore, CCR7 has been considered as a tumor therapeutic target [36]. In PTC, CCR7 has been found to show strong correlations with tumor aggressiveness indicators and tumor size [37]. CD1B and CD86 both function in antitumor   Evidence-Based Complementary and Alternative Medicine immunity by mediating the T cell function and NK cell cytotoxicity [38,39]. CSD2RB has been identified as a recurrent and oncogenic hotspot gain-of-function mutation [40]. HLA-DQA1, a member of MHC Class II molecules, its elevated expression has been identified as an adverse indicator for worse prognosis for esophageal squamous cell carcinoma patients [41]. Loss of NOD2 may promote polarization of tumor-associated macrophages into the M2 phenotype (protumorigenic phenotype) [42]. LTA and LTB encode lymphotoxin alpha and beta proteins, which belong to the tumor necrosis factor (TNF) family together with TNFRSF9 and TNFSF11. Interaction between TNF-TNFR superfamilies can initiate costimulatory signals that play important roles in antitumor immune responses by accelerating the differentiation, survival, and clonal expansion of antigen-primed CD8+ T and CD4+ cells [43]. All of these genes were upregulated in immune-activated samples, suggesting a higher antitumor immunity in these samples. Additionally, several genes were found to be correlated with tumor infiltrating abundance of follicular helper T cells (B cell-help providers) [44], activated NK cells (antitumor cytotoxic lymphocytes) [45], and M1 macrophages (antitumor phenotype) [46], and these cells play important roles in antitumor response. is also confirmed the higher antitumor immunity in immune-activated samples. We further analyzed the drug sensitivity to identify possible drugs for the treatment of invasive thyroid carcinoma. Bortezomib is an inhibitor for proto-oncogene RET has been approved for treatment of medullary thyroid cancer [47]. MG132 is a proteasome inhibitor that can induce tumor cell apoptosis in thyroid carcinoma [48]. AZD8055 inhibits tumor progression in PTC by affecting mTOR activity [49,50]. Staurosporine is a bacterial alkaloid, which can inhibit growth of TPC-derived cell line by the synergy effect with rotenone [51]. Immunosuppressive samples showed higher IC50 for bortezomib, MG.132, staurosporine, and AZD8055, suggesting that patients with immunosuppressive status may be resistant to bortezomib, MG.132, staurosporine, and AZD8055 treatments.
High TMB has been proposed as a useful biomarker for predicting the response to immune checkpoint blockade [52]. However, high TMB has also been reported to be failed to predict the response to immune checkpoint blockade in all solid cancer types [53]. In this study, the TMB showed no significant differences between immune-activated and immunosuppressive samples. Among the top 20 mutation genes, mutation in the BRAF gene accounted for 78% of tumor samples, suggesting that BRAF mutation might be correlated with the development of invasive thyroid carcinoma.
is was consistent with previous studies. BRAF mutation is the most common mutation in PTC, and tumors with BRAF mutation showed significant correlations with advanced T-stage and more frequently central neck dissection [54]. Genes in the model were not found in the top 20

NK cells activated NK cells resting T cells gamma delta T cells regulatory (Tregs) T cells follicular helper T cells CD4 memory activated T cells CD4 memory resting T cells CD4 naive T cells CD8 Plasma cells B cells memory B cells naive
(g) Figure 7: Mutation analysis and drug sensitivity analysis. (a) e drugs with IC50 < 1 for all tumor samples in TCGA cohort; (b) heatmap shows the differences on IC50 value (evaluation for drug sensitivity) for immune-activated and immunosuppressive patients; (c) risk score significantly differed between immune-activated and immunosuppressive patients; (d) tumor mutation burden shows no differences between immune-activated and immunosuppressive patients; (e) statistical result of somatic mutation for all tumor samples in TCGA cohort; (f ) waterfall plot shows the top 20 mutation genes; (g) correlation analysis between tumor-infiltrating immune cells evaluated by CIBERSORTx and genes in model. mutation genes, suggesting that the expression differences of these genes might not be caused by gene mutation.

Conclusion
We identified a 12-gene-based immune signature, which contributes to predicting the immune-activated or immunosuppressive status for patients with invasive thyroid carcinoma. Patients with immune-activated invasive thyroid carcinoma are more likely to benefit from immunotherapy. Immunosuppressive patients might be resistant to bortezomib, MG.132, Staurosporine, and AZD8055 treatments. ese findings will provide clues for improving the clinical management of invasive thyroid carcinoma.
Data Availability e dataset supporting the conclusions of this manuscript is included within the manuscript.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Authors' Contributions
All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Runming Pang and Chunxin Qin. e first draft of the manuscript was written by Runming Pang and all the authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.  Figure S2: GSEA for DEGs between immune-activated and immunosuppressive groups. A, the bubble diagram shows the top 20 activated and suppressed pathways in GSEA analysis; B, the top 10 pathways in GSEA ranked by NES value. Table S1: the 29 immune related gene sets.