A Novel Immune-Related Prognostic Signature for Thyroid Carcinoma.

BACKGROUND: Recent studies have confirmed that immune-associated genes perform a crucial function in recurrence and metastasis of thyroid carcinoma. A reliable immune-related prognostic signature for patients with thyroid cancer is needed. This study constructed a novel immune-related prognostic signature for thyroid cancer and evaluated its prognostic value by bioinformatics analysis.
METHODS: In this study, we anatomized differentially expressed immune-associated genes from The Cancer Genome Atlas database. The samples from The Cancer Genome Atlas database were randomly divided into training set and test set. A novel immune-related prognostic signature for thyroid cancer was developed by least absolute shrinkage and selection operator and Cox regression analysis: Risk score = (0.6846 × expression value of C-X-C motif chemokine ligand 5 [CXCL5]) + (1.1556 × expression value of Azurocidin 1 [AZU1]) + (-0.3156 × expression value of nucleotide binding oligomerization domain containing 1 [NOD1] + (0.0542 × expression value of TNF Receptor Superfamily Member 11b [TNFRSF11B]) + (0.0952 × expression value of VGF nerve growth factor inducible [VGF]). The established prognostic signature was evaluated based on training set and test set by survival curves, receiver-operator characteristic curves, risk score, survival status, heatmap, and independent prognostic analysis. Meanwhile, we appraised the correlation between target immune-associated genes and clinical stage, tumor-infiltrating immune cells respectively.
RESULTS: Five immune-associated genes were used for constructing an immune-related prognostic signature by least absolute shrinkage and selection operator, univariate, and multivariate analysis. Survival curves, receiver-operator characteristic curves, and independent prognostic analysis showed the signature had significant prediction value. Clinical and immune cell correlation analyses indicated that target immune-associated genes may participate in tumor immune infiltration and tumor progression.
CONCLUSIONS: We constructed a novel 5 immune-associated genes signature for predicting the prognosis of patients with thyroid cancer, which may help clinical workers evaluate individualized therapy and prognosis.

Introduction

Thyroid cancer is the most common endocrine tumor[1] that has been increasing fast in the recent years all over the world.[2,3] Thyroid cancer screening and overdiagnosis may be the cause of rising morbidity.[4] Although most patients have a favorable prognosis, there are still few recurrences or deaths that occurred after the treatment of thyroid cancer.[5] There are many significant prognostic factors for patients with thyroid cancer, including age, gender, tumor histology, and presence of lymph node or distant metastases.[6] Moreover, recent studies have demonstrated that biomarkers play a significant role in the occurrence and progression of thyroid cancer.[7]

Now, high-throughput expression data provide significant capabilities to appraise the relationship between genes and clinical features of patients with thyroid cancer. Pak et al established a risk scoring system that will help to choose suitable treatment for papillary thyroid carcinoma (PTC).[8] Kim et al identified 3 immune-related metagene signature can classify thyroid cancer into several molecular subtypes.[9] Brennan et al developed a prognostic signature for intermediate-risk thyroid cancer.[10] Ma et al suggested that metabolic deregulations mediate dedifferentiation of thyroid cancer.[11]

Immune-related biomarkers and tumor-immune microenvironment act a significant role in the occurrence and evolution of thyroid cancer. Immune-related key gene CLDN10 has a strong correlation with lymphatic metastasis of thyroid cancer.[12] Immune-associated genes (IAGs) play an extremely essential role in the classification of papillary thyroid cancer.[13] Distinct immune microenvironments show significant value with pathological invasiveness and risk stratification of thyroid cancer.[14]

In this study, we utilized transcriptome data from The Cancer Genome Atlas (TCGA) and constructed an independent prognostic IAGs signature, which showed favorable value for evaluating the prognosis of patients with thyroid cancer.

Materials and Methods

Data Acquisition

We got the immune gene list from Immunology Database and Analysis Portal (ImmPort; https://immport.niaid.nih.gov).[15] Data of messenger RNA (mRNA) expression and the corresponding clinical data were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/). Tumor Immune Estimation Resource (TIMER) is an online resource for analysis of immunocyte infiltration,[16] which has been applied in previous research.[17] In this study, we evaluated potential association between immune infiltrates and IAGs expression by TIMER.

Differential Expression of IAGs Filtrating

Gene expression data were downloaded in Fragments per Kilobase Million format. The genes were filtrated by the cutoff criteria of P < .05 and |log2 fold change (FC)| > 2 between normal samples and thyroid cancer samples by the R package limma.[18,19]

Functional Enrichment Analysis

Differential expressed IAGs were annotated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis through R package clusterprofiler.

Construction of Immune-Related Prognostic Signature

First, samples were randomly divided into training set and test set by R package createDataPartition. Training set was used for constructing signature, while test set was used to evaluate signature. Then, the samples from training set were used to filtrate the IAGs by univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO). The selected IAGs were analyzed by multivariate stepwise Cox regression analysis. Finally, an immune-related prognostic signature was constructed by the genes and its regression coefficient (β). β is the weighted correlation coefficients in the multiple Cox regression analysis. The risk score = Σ(βi × Expression i) (i = the prognostic IAGs). Furthermore, we calculated the risk score of each sample from training set and test set, and samples were divided into high-risk group and low-risk group by median risk score of training set.

Evaluation of Immune-Related Prognostic Signature

Survival curves of training set and test set were drawn by Kaplan-Meier method and compared by log-rank test. Receiver–operating characteristic curves (ROC) of training set and test were plotted by R package survivalRoc.[20] Meanwhile, risk score, survival status, and heatmap of signature in training set and test set were used to visually evaluate the performance of this signature. Independent prognostic analysis was used to predict whether risk score was an independent prognostic factor in patients with thyroid cancer.

Clinical and Immune Cell Correlation Analysis

Immune cell correlation analysis was conducted by online tool (TIMER) and clinical correlation analysis was conducted by R package Beeswarm, which helped us understand the roles of target IAGs in tumor progression and immunocyte infiltration.

Statistical Analyses

All statistical analyses were performed by version 3.5.3 of R software (https://www.r-project.org/) and IBM SPSS statistic version 21. Quantitative variables were compared using the t test, and categorical variables were compared using the Fisher exact test or Pearson χ2 test. Survival curves were plotted by the Kaplan-Meier method and compared by the log-rank test using the R package survival.

Results

Differential Expressed IAGs

We downloaded mRNA expression data from TCGA data set, which include 509 thyroid cancer samples and 58 normal samples. Then, we combined clinical characteristics with mRNA expression data (Table 1). In the meanwhile, we obtained 2498 IAGs from the ImmPort. Finally, 82 differentially expressed IAGs between thyroid cancer samples and normal samples were selected for further research (P value < .05 and |log2 FC| > 2), which included 63 upregulated and 19 downregulated genes. The expression of 82 differentially expressed IAGs were visualized by Volcano plot and Heatmap plot (Figure 1).

Table 1.

Clinical Characteristics of Patients With Thyroid Cancer.

Clinical characteristic	Variables	Total (500)	Percentages (%)
Age at diagnosis	47.32 ± 15.84
Gender	Female	367	73.4
	Male	133	26.6
Stage	I	280	56
	II	52	10.4
	III	113	22.6
	IV	54	10.8
	Unknow	1	0.2
T	T1	140	28
	T2	166	33.2
	T3	169	33.8
	T4	23	4.6
	Tx	2	0.4
N	N0	226	45.2
	N1	225	45
	Nx	49	9.8
M	M0	280	56
	M1	9	1.8
	Mx	211	42.2

Figure 1.

Volcano plot (A) and Heatmap (B) of 82 differentially expressed immune-associated genes (IAGs) in normal samples and thyroid cancer samples.

Eighty-two differential expressed IAGs were annotated according to GO database. Top 10 results of the GO functional analysis are shown in Figure 2A, which showed that cell chemotaxis, myeloid leukocyte migration, leukocyte chemotaxis, and leukocyte migration were relatively enriched for biological process. As for cellular component, receptor complex was relatively enriched, while for molecular functions were receptor ligand activity, cytokine activity, and cytokine receptor binding. Besides, these genes were also annotated according to KEGG pathways. As shown in Figure 2B, 7 KEGG pathways were most positively related to IAGs, such as cytokine–cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, and interleukin (IL)-17 signaling pathway.

Figure 2.

Top 10 results of the GO functional enrichment analysis of differentially expressed immune-associated genes (IAGs). The node color changes gradually from red to blue in ascending order according to the adjusted P values. The size of the node represents the number of counts (A). Enrichment of 7 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of differential expressed IAGs (B).

First, 82 differentially expressed IAGs from training set were subjected to univariate Cox regression analysis. The results showed that 14 IAGs were closely related to the survival of patients with thyroid cancer (Table 2). Then, LASSO analysis was used to remove highly relevant IAGs, which helped us avoid the overfitting of signature. Finally, 9 IAGs were subjected to multivariate stepwise Cox regression (Figure 3). The results showed that 5 IAGs (CXCL5, AZU1, NOD1, TNFRSF11B, and VGF) were prognostic genes, which included 4 latent dangerous IAGs and 1 latent protective IAG (Figure 4). And immune-related signature was constructed as follows: risk score = (0.6846 × expression value of CXCL5) + (1.1556 × expression value of AZU1) + (−0.3156 × expression value of NOD1) + (0.0542 × expression value of TNFRSF11B) + (0.0952 × expression value of VGF). We used this signature to calculate risk scores for each sample in training set and test set. And samples in training set and test set were divided into high-risk group and low-risk group by median risk score of training set.

Table 2.

Univariate Cox Regression Analysis.

id	HR	HR.95L	HR.95H	P value
CXCL5	1.605863	1.29425	1.992502	1.68E-05
TNFRSF11B	1.051041	1.022537	1.08034	.000387
OGN	1.169625	1.070269	1.278205	.000542
CSF2	1.297897	1.103618	1.526377	.001623
VGF	1.073086	1.023027	1.125594	.003804
RETN	1.103924	1.031776	1.181117	.004143
BMP2	1.061575	1.018776	1.106173	.004429
APOD	1.036986	1.008222	1.066571	.01139
NOD1	0.8621	0.751897	0.988455	.033473
RXRG	0.950694	0.906711	0.99681	.036423
MET	0.979853	0.961044	0.99903	.039581
SDC4	0.996341	0.992836	0.999857	.041424
AZU1	1.856132	1.019263	3.380112	.04314
PLXNC1	0.802616	0.64578	0.997542	.047462

Abbreviations: CXCL5, C-X-C motif chemokine ligand 5; NOD1, nucleotide binding oligomerization domain containing 1; TNFRSF11B, TNF receptor superfamily member 11b; HR, Hazard ratio; OGN, Osteoglycin; CSF2, Colony stimulating factor 2; VGF, VGF nerve growth factor inducible; RETN, Resistin; BMP2, Bone morphogenetic protein 2; APOD, Apolipoprotein D; RXRG, Retinoid X receptor gamma; MET, MET proto-oncogene, receptor tyrosine kinase; SDC4, Syndecan 4; AZU1, Azurocidin 1; PLXNC1, Plexin C1.

Figure 3.

Least absolute shrinkage and selection operator (LASSO) analysis.

Figure 4.

Multivariate analysis between immune-associated genes (IAGs) and overall survival (OS) in patients with thyroid cancer. Four latent dangerous IAGs.

Survival curves of training set and test set showed that patients in high-risk group had worse overall survival compared to those in low-risk group (P value < .05), and ROC curves of training set and test set indicated that this signature had strongly predictive power (Figure 5). In Figure 6, risk score, survival status, and heatmap of signature in training set and test set showed that this signature had significant value for evaluating the prognosis of patients with thyroid cancer. The results of independent prognostic analysis attested that risk score was a significant prognostic factor for thyroid cancer (Figure 7).

Figure 5.

Survival curves of training set (A) and test set (B). Receiver–operator characteristic (ROC) curves of training set (C) and test set (D).

Figure 6.

Risk score, survival status, and heatmap of signature in training set (A) and test set (B).

Figure 7.

Independent prognostic analysis for thyroid cancer samples. Univariate independent prognostic analysis (A) and multivariate independent prognostic analysis (B).

As shown in Figure 8, AZU1 and CXCL5 had a correlation with clinical stage (P value < .05), which indicated that AZU1 and CXCL5 may act a key role in tumor progression. Meanwhile, CXCL5, NOD1, AZU1, and TNFRSF11B showed a strong correlation with immune cells (Figure 9), which indicated that these genes may participate in the process of immune infiltration of thyroid tumor (correlation coefficient > 0.3). C-X-C motif chemokine ligand 5 and NOD1 showed a strong correlation with neutrophil and myeloid dendritic cells. AZU1 was closely related to myeloid dendritic cells, and TNFRSF11B was closely connected with CD8+ T cells.

Figure 8.

Clinical correlation analysis. Clinical correlation analysis of AZU1 and C-X-C motif chemokine ligand 5 (CXCL5) showed a strong correlation with clinical stage.

Figure 9.

Immune cell correlation analysis. C-X-C motif chemokine ligand 5 (CXCL5) and nucleotide binding oligomerization domain containing 1 (NOD1) showed strong correlation with neutrophil and myeloid dendritic cells. AZU1 had a strong correlation with myeloid dendritic cells. TNF receptor superfamily member 11b (TNFRSF11B) had a strong.

Discussion

Papillary thyroid carcinoma is the most common type, accounting for more than 75% of total thyroid cancer, which has the least aggressive ability.[21] Although it has an excellent prognosis, its recurrence and metastasis still happen. The American Joint Committee on Cancer TNM classification fails to adequately predict the prognosis of patients with papillary cancer, and consequently it is urgent to develop a novel prognostic signature. In this study, we anatomized 82 differentially expressed IAGs and explored the functions of these IAGs by the GO and KEGG pathway analysis. After univariate, LASSO, and multivariate analysis, 5 IAGs were finally selected to construct an immune-related prognostic signature. Survival curves, ROC curves, risk score, survival status, heatmap, and independent prognostic analysis demonstrated that this signature had significant value for thyroid cancer. Clinical and immune cell correlation analysis indicated that several IAGs may participate in tumor immune infiltration and tumor progression.

Functional enrichment analysis suggested that these IAGs are involved in significant biological process and signaling pathway as following: cytokine–cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, IL-17 signaling pathway, natural killer cell–mediated cytotoxicity, hematopoietic cell lineage, and pertussis. Some signaling pathways identified in this study have been reported in thyroid cancer in recent years. Feng et al suggest that HOXC10 upregulation contributes to human thyroid cancer and indicates poor survival outcomes through cytokine–cytokine receptor interaction and chemokine signaling pathway.[22] Chemokine signaling pathway is transduced by chemokine receptors expressed in the immune cells. When receptor is activated, the α and β subunits of G protein dissociate to activate diverse downstream pathways, which lead to cellular polarization and actin reorganization.[23,24]

The prognostic signature included 5 IAGs (CXCL5, AZU1, NOD1, TNFRSF11B, VGF). C-X-C motif chemokine ligand 5 encodes a protein that is a member of the CXC subfamily of chemokines, which plays a role in cancer cell proliferation, migration, and invasion.[25] Cui et al show that CXCL5–CXCR2 axis promotes the growth of PTC cells. Mechanistically, the CXCL5 promoted the proliferation of PTC cells, accelerated the G1/S transition, upregulated the expression of a group of S (DNA synthesis) or M (mitosis)-promoting cyclins and cyclin-dependent kinases (CDKs), and downregulated CDK inhibitors in PTC cells.[26] Nucleotide binding oligomerization domain containing 1 encodes a member of NOD-like receptor family of proteins. This protein plays an important role in innate immunity.[27,28] Mey et al suggest that NOD1 being significantly lower in tumor tissue compared to healthy tissues. Further investigation showed that NOD1 inhibits proliferation by directly targeting proto-oncogene SRC and inducing cell cycle arrest at G1 phase.[29] AZU1 encodes multifunctional inflammatory mediators, Azurophil granules, which play an important role in the innate immunity of humans.[30] VGF is specifically expressed in neuroendocrine cells. Recent studies have found that VGF differentially expressed and methylated in normal tissue and several tumor tissue.[31,32] TNFRSF11B is also known as OPG and TR1, encoding a member of the TNF-receptor superfamily. Deligiorgi et al show that OPG involving the mutated KRas proteins and c-Fos overexpression is related to central lymph node metastasis in thyroid cancer.[33]

In conclusion, we identified an immune-related signature based on training set and test set from TCGA database which showed a favorable value for predicting the prognosis of patients with thyroid cancer. When patients after treatment intend to know their survival risk, we could provide individual advice by this signature.

However, there are still some limitations in this study. First, because Gene Expression Omnibus lacks corresponding survival data, we are unable to evaluate the signature again by external validation. Second, some of these IAGs have not previously been identified as biomarkers for thyroid cancer, and further experimental exploration is deserved.

Conclusions

We constructed a novel immune-related prognostic signature for thyroid cancer, which may help clinical workers evaluate the individualized prognosis of patients.