The landscape of tumors-infiltrate immune cells in papillary thyroid carcinoma and its prognostic value

Data sources and processing

We downloaded clinical and RNA sequencing (RNA-seq) data (Fragments Per Kilobase Million (FPKM) values) about PTC from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). From these data, we selected 481 patients with an effective DFS (Kanazawa & Kammori, 2019) for our study. DFS refers to the time from randomization when the disease recurred or when the patient died as a result of disease progression. The microarray datasets (GSE35570, GSE33630, and GSE60542) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). We used the “ComBat” algorithm (Mansbach et al., 2020) to reduce the batch effects from non-biological technical biases between GSE35570, GSE33630, and GSE60542. We combined these three datasets and referred to them as the GEO combined sets.

Consistent clustering of tumor-infiltrating immune cells

Using previous pan-cancer analyses (Charoentong et al., 2017), we obtained a gene set of 28 tumor-infiltrating immune cells. We used a Single Sample Gene Set Enrichment Analysis (ssGSEA) algorithm (Yi et al., 2020) and the ‘GSVA’ package (Hänzelmann, Castelo & Guinney, 2013) in R to obtain the enrichment score of each patient. Each sample’s immune and stromal levels (immune score, stromal score, ESTIMATE score, and tumor purity) were obtained using an estimate algorithm (Ke et al., 2020). We performed 50 iterations of consistent clustering with a resample rate of 80% for each patient based on the immune cell infiltration (ICI) pattern using “ConsensusClusterPlus” package (Wilkerson & Hayes, 2010). We used the “pheatmap” package to visualize the clustering results. We designated the three clusters, ICIcluster 1, ICIcluster 2, and ICIcluster 3, according to their ICI levels (from low to high).

ICI cluster differentially expressed gene (DEG) analysis

Considering that some of the samples had a gene expression level of 0 which would lead to bias in the results, we only selected genes with mean gene expression levels greater than or equal to 0.1 for use in the subsequent study. DEGs were screened and identified from the three ICI clusters using |Log2 fold change (LogFC) > 1.8| and a false discovery rate (FDR) < 0.01,which was implemented by the Wilcox test. Univariate COX analysis (Tibshirani, 2009) was used to determine the prognostic DEGs with p < 0.01.

Identifying gene clustering and ICIscore generation

Using the consistent clustering method, we identified three gene clusters. The genes that positively and negatively correlated with ICIcluster were designated as the upregulated and downregulated genes, respectively. We used a principal-component analysis (PCA) algorithm (Ringnér, 2008; David & Jacobs, 2014; Konishi et al., 2019) to extract the first principal components of the upregulated and downregulated genes. Next, we defined the ICIscore of each patient:

$\mathrm{I}\mathrm{C}\mathrm{I}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum \mathrm{P}\mathrm{C}{\mathrm{A}}_{\mathrm{u}}-\sum \mathrm{P}\mathrm{C}{\mathrm{A}}_{\mathrm{d}}$

where u refers to the upregulated genes and d refers to the downregulated genes. We selected the best cutoff using the “survminer” package. According to the cutoff, we divided the patients into two groups (ICIscore low and ICIscore high).

Constructing the prognostic signature and Riskscore generation

We randomly divided TCGA patients into two cohorts. Cox proportional hazard models are the most widely-used approach for modeling time to event data. They use the principle of recursive elimination to compute the instantaneous rate of an event occurrence (Asano, Hirakawa & Hamada, 2014). We constructed a 5-mRNA prognostic signature using the COX-PH algorithm. Next, we defined the Riskscore (Lee, 2008) of each patient:

$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum _{i=1}^{n}Expi\times Coef$

where n refers to the number of genes in the signature, Expi refers to the level of gene expression in the signature, and Coef refers to the estimated regression coefficient value from the COX-PH algorithm. According to their Riskscore median, the patients were divided into two groups (Riskscore low and Riskscore high). The signature was also fitted in the testing set and the total set.

Statistical analysis

All statistical analyses were conducted using R (3.6.1) software. We used the Mann–Whitney U test to compare two groups, and the Kruskal–Wallis test to compare more than two groups. The Kaplan-Meier (K-M) survival curves (You et al., 2018) were used to generate different groups of DFS curves. The receiver operating characteristic curves (ROC) (Blangero et al., 2020) was used to identify the Riskscore prognosis. Two-tailed P values less than 0.05 were considered statistically significant.

The landscape of tumor-infiltrating immune cells

The detailed flow chart of this study is shown in Fig. 1. Using the gene set associated with 28 types of infiltrating immune cells (Table S1) (Charoentong et al., 2017) and the ssGSEA algorithm, we obtained the enrichment score of the tumor-infiltrating immune cells in each PTC sample. Using consistent clustering, we divided the PTC patients into three ICI subtypes (Fig. 2A and Fig. S1). Based on the correlation analysis, we generated a correlation matrix of the tumor-infiltrating immune cells (Fig. 2B). Almost all of the items showed a positive correlation. The most significant correlation was shown between Myeloid-derived suppressor cells (MDSC) and other cells, especially MDSCs and effector memory CD8 T cells (Khaled, Ammori & Elkord, 2013; Lechner et al., 2011). WE used the Kruskal–Wallis test to compare the enrichment score across the three ICI subtypes (Fig. 2C) and we found that there were significant differences across the three ICI groups, indicating that our clustering was successful. Additionally, we performed pairwise combinations of the three ICI subtypes, and the Mann–Whitney U test to compare the expression levels of two immune checkpoint-related genes (PD-L1 and CTLA-4) in every combination (Figs. 2D and 2E). To explore the biological processes of the different ICIclusters, we used GSEA enrichment analysis (Subramanian et al., 2005) on the different ICIclusters. The ICIcluster 2 pathway was significantly enriched compared to ICIcluster 1 (Fig. 3A). Similarly, Fig. 3B shows that the ICIcluster 3 pathway was significantly enriched compared to ICIcluster 2. Detailed enrichment results can be seen in Table S2.

The same approach was applied to the GEO combined sets. We used consistency clustering to divide patients into three groups (Fig. S2A). The correlation analysis of the 28 types of tumor-infiltrating immune cells is shown in Fig. S2B. Significant differences were shown across the three ICIclusters (Fig. S2C), which was consistent with the performance of TCGA set. In the analysis of the immune checkpoint-related genes (PD-L1 and CTLA-4), there were also significant differences across the ICIclusters (Figs. S2D–S2F). GSEA enrichment analysis was also performed for the different ICIclusters, and the results can be found in Figs. S3A and S3B. Detailed enrichment results can be seen in Table S3. In the most significant Top 10 enrichment pathway, the GEO combined set and TCGA set were consistent with one another. This suggests that our ICIclusters represented the tumor-infiltrating immune cell landscape of PTC patients and will lay the foundation for subsequent analysis.

Stratified DFS analysis of the three ICI subtypes

The three ICI subtypes showed different DFS conditions compared to the total set (Fig. 3C). ICIcluster 3 had the most significant immune cell infiltration and was associated with the best prognosis. It is notable that ICIcluster 2 showed moderate immune cell infiltration but was associated with the worst prognosis. Although ICIcluster 1 had the lowest immune cell infiltration, it was associated with a moderate prognosis. Such results were also reflected in the Age <45 and T3 + T4 groups (Fig. 3D). This indicates that high or low immune function has the same influence on tumor prognosis, and will ultimately inhibit tumor development.

Identifying DEGs associated with the three ICI subtypes

Since TCGA set had complete clinical information and prognostic indicators, we conducted subsequent studies using this dataset. In order to make the results more accurate, we selected genes with average gene expression levels greater than or equal to 0.1 in order to determine the underlying molecular differences across the different ICI subtypes. We used the Wilcox-test to identify the DEGs associated with the three ICI subtypes. With |LogFC > 1.8| and FDR < 0.01, we obtained 982 DEGs: 752 upregulated genes and 230 downregulated genes (Fig. S4, Table S4). Upregulated genes were those that were highly expressed in at least two ICI subtypes with higher immune cell infiltration. The rest were labeled downregulated genes. To determine the underlying biological functions of the DEGs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (Xing et al., 2020; Wang et al., 2020) of the upregulated genes and downregulated genes, respectively. The most significant enrichment results are summarized in Figs. 4A–4D, and detailed content is provided in Table S5. Using P < 0.01 as the cutoff, we performed univariate COX analysis to identify the DEGs associated with prognosis. Finally, we obtained 15 DEGs: 10 upregulated and five downregulated genes (Table 1).

Identifying immune gene subtypes

Using our consistent cluster analysis of the 15 genes listed above, we obtained three gene clusters (Fig. 5A). In order to explore the abundance of immune cell infiltration across the three gene subtypes, we performed Kruskal–Wallis tests on the three subtypes. We found that Genecluster 1 and Genecluster 3 had more tumor-infiltrating immune cells and that Genecluster 3 was more significant than Genecluster 1 (Fig. 5B). This difference at the immune level was also shown in the expression levels of the two immune checkpoint-related genes (PD-L1 and CTLA-4; Figs. 5C–5D). Furthermore, we also explored the prognostic differences across the three gene subtypes. Genecluster 1 and 2 were associated with a favorable prognosis while Genecluster 3 was associated with the worst prognosis (Fig. 6A). It is worth noting that in addition to the worst prognosis, Genecluster 3 also had the most immune infiltrating cells. The immune system’s specific mechanism involved with tumor occurrence and development is worthy of further investigation, as its activation may have both favorable and unfavorable results. According to clinical stratification studies, similar results can also be seen in patients who were female, age <45, age >45, N1, Stage I + Stage II, and T1 + T2 (Figs. 6B–6G).

Generation and validation of the ICIscore groups

We used the PCA algorithm to calculate the ICIscore for each PTC sample. The ICIscore consisted of two parts: (1) the PCAu calculated from upregulated genes and (2) the PCAd calculated from downregulated genes. The PCAu and PCAd of each sample were the first principal components of the relevant gene. Finally, using the ICIscore formula, we obtained the ICIscore of each sample. With the “survminer” package, we calculated the optimal ICIscore cutoff as 2.17 (Fig. 7A) and divided the patients into two groups (ICIscore high and ICIscore low). Since the ICIscore correlated with immune status and the immune score reflected the level of immune cell infiltration, it was reasonable to believe that there was a correlation between these values. In order to explore this relationship, we drew a scatter plot (Fig. 7B) that showed a highly positive correlation between the two, indicating that ICIscores could reflect the level of immune cell infiltration to some extent. Although immune scores do not have prognostic value, we found that ICIscores did correlate with prognosis in a follow-up study. The distribution of all patients is shown in Fig. 7C. Patients with high ICIscores were mainly distributed in Genecluster 3 with poor prognosis predicted. To determine whether different ICIscore groups had different prognoses, we plotted the K-M survival curve. We found that the ICIscore high group showed significantly poorer prognoses (Fig. 7D). To evaluate whether ICIscores reflected differences in immune levels across the groups, we compared the tumor-infiltrating immune cells between the two groups using a Mann–Whitney U test (Fig. 7E). Additionally, we selected PDCD1, CD274, CTLA-4, LAG3, and IDO1 as the immune checkpoint-related genes and IL1A, IL1B, IL6, IL8, IFNG, and TNF as the immune active factors (Ayers et al., 2017; Deng et al., 2012). The Mann–Whitney U test showed that most immune checkpoint-related genes, except PDCD1 and TNF, and immune active factors were highly expressed in the ICIscore high group (Fig. 7F). Based on the GSEA results, we also found that the “allograft rejection” and “graft-versus-host disease” pathways were significantly enriched in the ICIscore high group, while the “maturity onset diabetes of the young” and “olfactory transduction” pathways were significantly enriched in the ICIscore low group (Figs. 7G–7H, Table S6).

Generation and validation of the Riskscore group

To further explore the potential prognostic value of these genes, we randomly divided TCGA PTC patients into two groups at a 1:1 ratio. We used the COX-PH algorithm to construct a 5-mRNA prognostic signature (Fig. 8A, Table 2). We specified the Riskscore median cutoff of the training set to fit the testing set and the total set. Riskscore = MMP9 × 0.019 + AKR1C1 × −2.358 + PLA2G2E × 0.230 + CARTPT × 0.017 + SLC5A1 × 1.914. Previous literature suggested that these five genes were closely related to cancer prognosis. In addition, it has been reported that the elevated MMP-9 level in PTC patients is defined as a malignant factor of PTC (Zarkesh et al., 2018) and may be involved in the mechanism of ROCK1 in the occurrence and development of PTC (Luo et al., 2017). AKR1C1, PLA2G2E, CARTPT, and SLC5A1 have been shown to be abnormally expressed in cancer and influence prognosis (Zeng et al., 2017; Sato et al., 2014; Burgos, Iresjö & Smedh, 2016; Gao et al., 2019; Mojica, Luna-Vital & Gonzalez de Mejia, 2018; Lei et al., 2016).

The distribution of all patients is shown in Fig. 8B. Patients with high Riskscores had poor prognosis. Before verifying the signature’s prognostic effect, we compared the immune cell infiltration between the two Riskscore groups (Fig. 8C). There was a significant difference in the level of immune cell infiltration between the two groups, and the Riskscore high group showed more immune cell infiltration. Additionally, most immune checkpoint-related genes and immune active factors, except IL1A and IL6, were highly expressed in the Riskscore high group (Fig. 8D). We also performed GSEA enrichment analysis on the two groups. The “allograft rejection” and “asthma” pathways were significantly enriched in the Riskscore high group, while the “synthesis and degradation of ketone bodies” and “ on-homologous and end-joining” pathways were significantly enriched in the Riskscore low group (Figs. 8E and 8F). More details are shown in Table S7.

By plotting K-M survival curves for the training, testing, and the total sets, we found that the Riskscore high group had a significantly poor prognosis (Figs. 9A–9C). In order to reveal the prognostic value of the signature, we verified the prognostic Riskscore value by drawing ROC curves for the different time periods (1-year, 3-years, and 5-years; Figs. 9D–9F). The AUC of the training, testing, and total set were mostly greater than 0.7, indicating that our signature had good prognostic value.

ICIscore and Riskscore clinical correlation analysis

While exploring potential independent prognostic factors, we did not include pM stages in this study due to the presence of many missing values in pM stages. Using univariate and multivariate COX analysis, we found that ICIscores and Riskscores could be used as independent prognostic factors (Figs. 10A–10D). Additionally, pT, pN, and pStage may also be potential prognostic factors according to univariate COX analysis with P value < 0.05. The nomograms showed the different clinical traits of patients without scores,with ICIscore or with Riskscore (Figs. 10E–10G). This provides a clinical reference for these patients’ prognosis.