RNA-seq, clinical data, transcriptome data, SNV, CNV, and methylation data on patients with thyroid cancer were downloaded from The Cancer Genome Atlas on 23 April 2022. For RNA-seq data, samples without clinical follow-up information, survival time, or status were removed; Ensembl was converted into gene symbol, and the median expression of multiple GeneSymbols was used. The pathological types of thyroid cancer are shown in Table 1.

Data on a total of 15 pathways (immune, stromal, DNA damage repair, and oncogenic) and their corresponding genes were obtained from a previous study (Li and Wang, 2021).

ssGSEA analysis was used to calculate the scores of the 15 pathways, EMT pathways (HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION), and cytolytic activity (Rooney et al., 2015) using the R package “GSVA”.

Molecular subtyping was performed separately for the TCGA dataset samples via ConsensusClusterPlus 1.52.0 using the scores for the 15 pathways (Wilkerson and Hayes, 2010). A total of 500 bootstraps were completed with “pam” arithmetic and “pearson” distances. Each bootstrap involved TCGA dataset specimens (≥80%). The cluster number k was set from 2 to 10, and the optimum k was defined as per cumulative distribution function (CDF) and AUC. Differences in survival (KM) curves were analyzed according to the molecular subtypes. Similarly, the distribution differences in clinical characteristics were compared, and chi-square tests were conducted. p < 0.05 was defined as statistically significant. Principal component analysis (PCA) was also performed to test the rationality of the molecular subtype distributions.

The proportions of 22 tumor-infiltrating immune cells (TIICs) were calculated using the CIBERSORT algorithm in all malignant tumor samples. ImmuneScore, StromalScore, ESTIMATEScore, and TumorPurity were determined using the ESTIMATE algorithm. ssGSEA identified scores for 28 kinds of immune cells.

For 172 tumor driver genes (159 of which had copy data) (Gao et al., 2013), we used GISTIC2 to analyze the changes in copy number. Those with a ratio >0.2 were considered Gains, while those with a ratio <0.2 were considered Losses; and the rest were considered to be Diploid. SNV was determined using maftools. Methylation of 450K in seven EMT genes and two mismatch repair genes was determined using the KNN algorithm in the impute package.

The TIDE (Jiang et al., 2018; Fu et al., 2020) algorithm (http://tide.dfci.harvard.edu) was used to evaluate the exclusion of CTL and dysfunction of tumor infiltration cytotoxic T lymphocytes (CTL) by immunosuppressive factors.

The sensitivity to traditional medicines (IC50 values) was predicted using pRRophetic (Geeleher et al., 2014).

In the TCGA dataset, thyroid cancer samples were randomly grouped into the training and test datasets in a 1:1 ratio. In the TCGA dataset, we identified pathway genes and pathways with Pearson correlations below the threshold |R| > 0.4, p < 0.05 to obtain related genes. In the training dataset, univariate Cox analysis was performed to screen genes related to prognosis. LASSO Cox regression in the glmnet package in R language and stepAIC in the MASS package were performed to select the best prognostic genes. A penalty coefficient λ of the optimal value and genes for the model development were determined through 10-fold cross-validation for a total of 1000 iterations. The risk scores for each were calculated using the following formula: I M s c o r e = ∑ β i × E x p i ,

where βi refers to the Cox hazard ratio coefficient of mRNA and Expri is the expression level of a gene. Samples in the training dataset were assigned into two groups of high-risk and low-risk based on the optimal segmentation point cutoff, which was determined using the survminer package. Simultaneously, the effectiveness and robustness of the prognostic risk model were validated in test and entire TCGA datasets. Survival differences among the risk groups were evaluated using Kaplan–Meier (KM) curves combined with log-rank tests. The performance of IMscore in pan-cancer, immunotherapy datasets (IMvigor210 and GSE91061) was also evaluated.

Sangerbox assisted with this article (Shen et al., 2022).

The software packages used in this study were implemented in R software (version 4.2.2; https://www.r-project.org/). A p-value < 0.05 was considered statistically significant.

Based on scores in 15 pathways, three molecular subtypes (Immune-Enrich (E), Stromal-Enrich (E), and Immune-Deprived (D)) in thyroid cancer were identified by ConsensusClusterPlus for k = 3 (Figure 1A). The PCA results showed that the three molecular subtypes had distinct boundaries, indicating the rationality of the subtype classification (Figure 1B). Samples in Immune-D showed better OS, while those in Immune-D showed better progression-free survival (PFS) (Supplementary Figure S1A). The distribution of clinical features of the three molecular subtypes indicated the significance of the T and N stages (Supplementary Figure S1B).

The results of the ESTIMATE analysis showed higher and lower ImmuneScore, StromalScore, and ESTIMATEScore in Immune-E and Immune-D, respectively (Figures 2A–C). TumorPurity was lower in Immune-E (Figure 2D). A higher EMT score was observed in Stromal-E (Figure 2E). The cytolytic activity score was increased in Immune-E (Figure 2F). In total, 28 kinds of immune cells showed higher scores in Immune-E (Figure 2G), while 18 of 22 immune cells also had higher scores in Immune-E (Figure 2H).

Furthermore, the expression levels of PDCD1, CTLA4, LAG3, and CD274(PD-L1) were upregulated in Immune-E (Figures 3A–D). The expression analysis of MHC genes showed increases in 21 genes in Immune-E (Figure 3E). These results indicated higher immune infiltration in Immune-E.

Next, we analyzed the gene mutations among the molecular subtypes. The results demonstrated that 60 genes among 172 tumor-driving genes showed varying degrees of mutation in the three molecular subtypes (Figure 4A). The TMB was higher in Immune-D compared to Stromal-E (Figure 4B). Tumor driver gene mutations and wild-type samples used for KM analysis showed better survival outcomes in samples with CSMD1 and ERBB3 wildtype compared to samples with CSMD1 and ERBB3 mutations (Figure 4C). CNV analysis of 159 genes showed copy number amplification and deletion in 22 genes in the three molecular subtypes (Figure 5A). Expression analysis of the corresponding genes in CNV groups of DOLPP1, PLEKHA6, PTEN, and MNDA demonstrated that the four genes had higher expression levels in the Gain group and low expression in the Loss group (Figure 5B).

A total of seven EMT genes and two mismatch repair genes were used to calculate the 450K beta values. The beta values of ZEB1, TW1ST1, CDH2, CDH1, and MLH1 differed among the three molecular subtypes (Figure 6A). Pearson correlation analysis of gene expressions and beta values showed that ZEB1, VIM, CDH2, CDH1, and CLDN1 expressions were negatively correlated with beta value (Figure 6B). The beta value of the cg probe site in CDH1 was higher in Immune-E (Figure 6C). Similarly, the beta value of the cg probe site was negatively correlated with CDH1 expression (Figure 6D).

We used TIDE (http://tide.dfci.harvard.edu/) software to evaluate the potential clinical effect of immunotherapy according to the molecular subtypes. The TIDE score was lower for Immune-E, indicating that Immune-E may be more suitable for immunotherapy. Moreover, 47% of samples in Immune-E showed immunotherapy response, a proportion higher than those in Stromal-E and Immuno-D (Figure 7A). The IC₅₀ values for cisplatin, erlotinib, sunitinib, paclitaxel, saracatinib, and dasatinib were lower in Immune-E, suggesting that Immune-E is more sensitive to those chemotherapeutic drugs (Figure 7B).

In the TCGA dataset, Pearson correlation analysis between genes in pathways and pathways identified 1784 genes with |R|>0.4 and p < 0.05. Then, in the TCGA training dataset, seven prognosis genes (p < 0.05) for thyroid cancer were screened from 1784 genes using univariate Cox analysis. Finally, four genes were used to construct a prognostic model (IMScore = 0.732*HSPA6+0.917*FLNC 1.083*CLDN2 0.966*E2F1) through lasso analysis and the stepAIC method.

In the TCGA training, testing, and entire TCGA datasets, samples were classified into high and low IMScore groups using the cutoff. KM curve analysis showed that patients in the low group had longer survival times. Moreover, in terms of PFS, DFI, and DSS, the low group showed better PFI (p = 0.04) and DSS (p < 0.0001) (Figure 8A). The IMScores were higher in the Immune-D and Stromal-E subtypes (Figure 8B).

Among 32 cancer types in the TCGA dataset, high IMScore survival times were shorter than low IMScore survival times except for TGCT and UCS (Supplementary Figure S2). We validated the prediction effect of IMScore in the immunotherapy datasets IMvigor210 and GSE91061. In the IMvigor210 dataset, samples with a low IMScore had better survival outcomes, and the 0.5-, 1-, and 1.5-year AUCs were 0.58, 0.64, and 0.65, respectively (Figure 9A). The samples with low TIDE had better survival outcomes, and the 0.5-, 1-, and 1.5-year AUCs were 0.54, 0.57, and 0.57, respectively (Figure 9B). Samples with low PD-L1 also had better survival outcomes, and the 0.5-, 1-, and 1.5-year AUCs were 0.6, 0.6, and 0.59, respectively (Figure 9C). The prediction of the response to treatment showed AUCs of TIDE, PD-L1, and IMScore of 0.58, 0.57, and 0.67, respectively (Figure 9D). In the GSE91061 dataset, samples with a low IMScore had better survival outcomes, and the 0.5-, 1-, and 1.5-year AUCs were 0.59, 0.75, and 0.75, respectively (Figure 9E). The samples with low TIDE had better survival outcomes, and the 1-, 2-, 2.5-year AUCs were 0.61, 0.58, and 0.59, respectively (Figure 9F). The survival outcomes did not differ significantly between the low- and high- PD-L1 groups, and the 0.5-, 1-, 1.5-year AUCs were 0.54, 0.57, and 0.57, respectively (Figure 9G). The prediction of the response to treatment showed AUCs of TIDE, PD-L1, and IMScore of 0.58, 0.55, and 0.61, respectively (Figure 9H). The results of the aforementioned analyses demonstrated the better prediction effect of the IMScore compared to TIDE.

First, univariate analysis showed that age, gender, TNM stage (p < 0.001), stage, and IMScore were significantly associated with a shorter OS in patients with thyroid cancer (Figure 10A). Then, we established a nomogram model that included the important predictors in the Cox analysis to predict the prognosis of thyroid cancer (Figure 10B). The calibration curve showed good concordance between the predicted and observed values of 1-, 3-, and 5-year OS (Figure 10C). The decision curve showed that the nomogram had the best prediction performance for the prognosis of thyroid cancer (Figure 10D).