Skip to main content

A novel signature to predict thyroid cancer prognosis and immune landscape using immune-related LncRNA pairs

Abstract

Background

Thyroid cancer (TC) is the most common endocrine malignancy worldwide. The incidence of TC is high and increasing worldwide due to continuous improvements in diagnostic technology. Therefore, identifying accurate prognostic predictions to stratify TC patients is important.

Methods

Raw data were downloaded from the TCGA database, and pairwise comparisons were applied to identify differentially expressed immune-related lncRNA (DEirlncRNA) pairs. Then, we used univariate Cox regression analysis and a modified Lasso algorithm on these pairs to construct a risk assessment model for TC. We further used qRT‒PCR analysis to validate the expression levels of irlncRNAs in the model. Next, TC patients were assigned to high- and low-risk groups based on the optimal cutoff score of the model for the 1-year ROC curve. We evaluated the signature in terms of prognostic independence, predictive value, immune cell infiltration, immune status, ICI-related molecules, and small-molecule inhibitor efficacy.

Results

We identified 14 DEirlncRNA pairs as the novel predictive signature. In addition, the qRT‒PCR results were consistent with the bioinformatics results obtained from the TCGA dataset. The high-risk group had a significantly poorer prognosis than the low-risk group. Cox regression analysis revealed that this immune-related signature could predict prognosis independently and reliably for TC. With the CIBERSORT algorithm, we found an association between the signature and immune cell infiltration. Additionally, immune status was significantly higher in low-risk groups. Several immune checkpoint inhibitor (ICI)-related molecules, such as PD-1 and PD-L1, showed a negative correlation with the high-risk group. We further discovered that our new signature was correlated with the clinical response to small-molecule inhibitors, such as sunitinib.

Conclusions

We have constructed a prognostic immune-related lncRNA signature that can predict TC patient survival without considering the technical bias of different platforms, and this signature also sheds light on TC’s overall prognosis and novel clinical treatments, such as ICB therapy and small molecular inhibitors.

Peer Review reports

Introduction

In recent decades, the incidence of thyroid cancer (TC) has been increasing steadily over the past 30 years worldwide [1]. With the development of diagnostic technology, the detection rate of TC in the world has increased year by year, ranking fifth in the incidence of malignant tumors [2], and the growing incidence of TC is raising serious public health issues worldwide. Even though most TCs have a relatively good prognosis, approximately 10% of patients with differentiated thyroid cancer (DTC) may progress to invasive disease, 5% progress to distant metastasis, and approximately 20–30% may relapse [3].

The immune system is now considered to have an important role in the elimination of cancer cells and sheds light on the mechanisms of cancer–immune evasion, contributing to tumor outgrowth [4, 5]. Although immunosurveillance prevents the development of most tumors in normal individuals, cancer cells can still deploy escape strategies such as initiating immune checkpoints [6]. Immune checkpoint inhibitors (ICIs) can restore the normal work of the immune system by suppressing the signal of "rest from work", which is sent by immune checkpoints, and then launch an attack on tumor cells [7]. Immune checkpoint blockade (ICB) immunotherapy has shown durable responses and improved clinical outcomes for patients across most malignancies, including thyroid cancer. Growing evidence has proven that ICI-related genes such as PD-1 and PD-L1 are expressed in ATC and DTC [8]. Pembrolizumab, an anti-PD-1 antibody, may be effective for patients with TC because of the KEYNOTE-158 trial, which shows that approximately 60% of patients in this program achieved disease control [9].

In addition, several studies have demonstrated that tyrosine kinase inhibitors (TKIs) can prevent the proliferation and tumorigenicity of thyroid cancer cells. Treatment with small-molecule tyrosine kinase inhibitors, including gefitinib, pazopanib, lenvatinib, and axitinib, are approved by the FDA to treat thyroid cancer [10]. Some clinical trials were performed to evaluate the efficacy and safety of these inhibitors [11]. Nevertheless, TKIs are promising new agents for the treatment of patients with thyroid cancer.

Long noncoding ribonucleic acids (lncRNAs) belong to the family of noncoding RNAs and are functionally defined as transcripts > 200 nucleotides in length with no protein-coding potential [12]. In the past few years, lncRNAs have attracted much attention due to their previously underappreciated transcriptional regulation, and they can function as both oncogenes and tumor suppressors [13]. It is clear that lncRNAs play crucial roles in the regulation of various biological and pathological behaviors of malignant tumors, especially in tumorigenesis and progression [14, 15].

Moreover, a growing body of studies has shown that the expression of lncRNAs is linked to the immune response and tumor progression. In addition, lncRNAs have an essential role in the development of immune cells and in pathogen response pathways. Notably, individual lncRNAs can act functionally through modular domains and often link protein activity to DNA or RNA targets through interactions with both [16]. Aberrant lncRNA expression has been observed in various cancers [17]. Recently, an increasing number of lncRNAs have been discovered to play an important role in TC tumorigenesis and development. Peng et al. [18] summarized representative lncRNAs in thyroid cancer, such as NEAT1, HOTAIR, PTCSC2, GAS8-AS1, MEG3, BANCR, GAS5, and MALAT1 et al., which are highly linked to the biological behavior of thyroid cancer and show significant value in diagnosis and treatment. Recent study showed that the expression level of lncRNA n384546 was upregulated in TC patients and that its interference inhibited cancer cell proliferation, invasion, and migration [19]. Reliable prognostic models related to tumor immune infiltration may affect the diagnosis, evaluation and treatment decisions of tumors. LncRNAs have recently been used to establish prognostic signatures. Huang et al. [20] constructed an immune-related lncRNA signature to predict the survival outcome of patients with breast cancer. Xu et al. [21] built a signature based on seven immune-related lncRNAs, which showed reliable prognostic value in hepatocellular carcinoma and may predict the outcome of immune checkpoint blockade (ICB) therapy.

In this study, we established a novel method to construct a predictive signature for TC based on immune-related lncRNA pairs, which did not need to consider the technical bias of different platforms. Next, to explore the latent role of this signature, we integrated the immune model with clinical factors of TC patients to build a composite prognostic index, which allowed improved estimation of prognosis. We also explored the correlation between the signature and several aspects, such as immune cell infiltration, immunosuppressive biomarkers, and small-molecule tyrosine kinase inhibitor efficacy. Overall, our new model might provide insight into the prognosis and clinical treatment of TC.

Materials and methods

Data download and preprocessing and differential expression analysis

The gene expression data and clinical data of TC patients were downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/) to perform a comprehensive analysis. The RNA-seq profiles of 558 cases comprised 500 thyroid carcinoma samples and 58 normal thyroid samples. We also included clinical data from 504 TC patients with an overall survival time of 0 days. Samples from the TCGA database were divided randomly into a training set (n = 252) and an internal test set (n = 252) at a ratio of 1:1. Next, the data were annotated to distinguish the mRNAs and lncRNAs for further analysis. We retrieved 2483 immune-related genes (IRGs) from the ImmPort database (https://immport.niaid.nih.gov) and performed Pearson correlation to identify the association between immune-modulating genes and all lncRNAs. LncRNAs with a correlation coefficient |R|≥ 0.4 and P ≤ 0.001 were considered related to immune genes and were used for further analysis. Screening of the DEirlncRNAs was based on logFC > 1 and FDR < 0.05 using the limma R package [22].

Pairing immune-related lncRNAs

We performed cyclical pairwise comparisons between the DEirlncRNA expression values to generate a score for each irlncRNA pair and established a 0-or-1 matrix. If the expression level of lncRNA1 was higher than that of lncRNA2, this pair was assigned a score of 1; otherwise, the score was 0. To prevent biases and unrepeatability, we validated the expression quantities of the lncRNA pairs. After removing lncRNA pairs scoring 0 or 1 in more than 80% and less than 80% of the total pairs, the remaining pairs were considered valid matches to build the prognostic signature.

Construction of a novel signature based on the DEirlncRNA pairs

Univariate Cox regression analysis (P < 0.001) was first implemented to identify the prognostic irlncRNA pairs with TC by using R (survival package). Subsequently, least absolute shrink age and selection operator (LASSO) Cox regression was conducted by utilizing the Glmnet R package to reduce the number of pairs [23]. Lasso regression was carried out to acquire a well-balanced prognostic model by running 1000 cycles. Then, 14 pairs of immune-related lncRNAs were ultimately identified as our prognostic model, and the risk score of each TC patient was calculated based on the following formula: RiskScore = . Here, n was the number of selected lncRNA pairs, βi was the coefficient of pair i, and Pi was the expression of lncRNA pair i. To validate the sensitivity and accuracy of the model, time-dependent ROC curve analysis was performed by utilizing the SurvivalROC package in R. The AUC was calculated for ROC curves, and sensitivity and specificity were calculated to assess the performance of the risk score. The optimal cutoff score based on the 1-year ROC curve was identified to separate the patients into low-risk and high-risk groups.

Reverse-transcription quantitative PCR

To validate the bioinformatics analysis results, 14 single lncRNAs of 14 DEirlncRNA pairs were selected to perform quantitative real-time PCR (qRT‒PCR) analyses. Twelve matched tumor and peritumor samples of PTC were collected from the First Affiliated Hospital of China Medical University. Total RNA was extracted from tissue samples using TRIzol (Invitrogen, United States), and then RNA was reverse transcribed into cDNA with the QuantiTect Reverse Transcription Kit (Takara, Shiga, Japan). Gene expression was analyzed by qRT‒PCR, which was performed with SYBR Premix Ex TaqII (Takara) and a LightCycler 480 system (Roche, Indianapolis, IN, USA). GAPDH was used for data normalization, and these data were analyzed by 2−ΔΔCT. The primer sequences are listed in Additional file 1: Table S1.

Validation of the risk assessment model

We estimated the prognostic capability of the risk score, sex, age and tumor stage (TNM stage and clinical stage) for overall survival in terms of time-dependent AUC values. The survival difference between subgroups was evaluated by Kaplan‒Meier survival analysis. The survival curve was visualized using the survival and survminer packages in R. The specific risk score value of each sample based on the signature was also used for visualization in R. We further performed principal component analysis (PCA) to assess the accuracy of the classification according to different risk scores. To predict the reliability and stability of the constructed model, Pearson’s χ2 test was used to assess the association between the signature and clinical features. The risk score differences between groups of conventional clinical features were compared by the Wilcoxon signed rank test. The analysis results were visualized by a strip chart and box diagrams. Moreover, we used univariate and multivariate Cox proportional hazards analyses of the risk score and other clinical characteristics to examine whether the model could be used as an independent variable. Then, the prognostic model was further validated in the internal validation set, including the training set and test set.

Correlation between the risk model and immune cell infiltration

The CIBERSORT algorithm was performed to estimate the tumor-infiltrating immune cell profiles of the samples in the TC dataset, followed by quality filtering, and only 509 tumor samples with P < 0.05 were selected for the following analysis. To further explore whether the risk model was related to immune cell infiltration in TC, other methods, such as TIMER, XCELL, QUANTISEQ, MCP-counter, EPIC, and CIBERSORT-ABS, were also used. We applied the Wilcoxon signed-rank test to compare the tumor-infiltrating immune cells among different risk groups. A box chart was used to visualize the results. To assess the relationship between the model and the infiltrated immune cells, Spearman correlation analysis was performed. The results of the Spearman rank correlation coefficients are shown in the lollipop-style mutation diagram. We visualized the results using the ggplot2 package in R.

Correlation between the risk model and the molecules related to ICIs

An ESTIMATE algorithm [24] was used to estimate the proportion of stromal cells (stromal score) and immune cells (immune score) that infiltrated the tumor tissue between low-risk and high-risk cases. We also compared the expression of the major markers of cytolytic activity (GZMA,PRF1) and the HLA gene between the high- and low-risk groups. The limma package was applied to explore the expression of ICI-related molecules among the different risk groups, and the results are shown in violin plots developed using the ggpubr package in R.

Analyses of the model ability in clinical treatment

Small-molecule inhibitors such as gefitinib, pazopanib, lenvatinib, axitinib, AMG-706 (also known as motesanib), tipifarnib, sunitinib, and sorafenib, which are now approved by the FDA, have shown great curative effects in thyroid cancer. To evaluate the ability of the signature to treat thyroid cancer, we calculated the IC50 values of these drugs in the THCA dataset. We performed the Wilcoxon signed-rank test to analyze the differences among groups, and the results are displayed in box plots utilizing the pRRophetic, ggpubr, and ggplot2 R packages.

Results

Identification of differentially expressed irlncRNAs (DEirlncRNAs)

This study was carried out based on the flowchart shown in Fig. 1a. First, RNA sequencing (RNA-seq) transcriptome data of TC were downloaded from the TCGA database. A total of 505 available TC patients were included in our study, and the baseline characteristics of all the patients are shown in Table 1. A total of 1148 irlncRNAs were identified by performing a coexpression analysis (shown in Additional file 2: Table S2). Next, we set the thresholds as the log fold change (FC) > 1.0 or < -1.0 and false discovery rate (FDR) < 0.05 and identified 200 DEirlncRNAs between tumor and normal thyroid tissues (Fig. 1b), including 117 upregulated and 83 downregulated genes (Fig. 1c).

Fig. 1
figure 1

Identification of immune-related lncRNAs in TC. (a) Analysis workflow of this study. (b) The heatmap shows the expression level of differentially expressed immune-related lncRNAs (DEirlncRNAs) in the TCGA dataset (blue, lower expression; red, higher expression). (c) The volcano plot represents the result of DEirlncRNA analysis of the THCA dataset (green: downregulated genes; black: no differentially expressed genes; red: upregulated genes)

Table 1 Clinical characteristics of the thyroid cancer patients

Construction of a novel prognostic signature based on DEirlncRNA pairs for TC

First, a total of 200 DEirlncRNAs were screened by pairwise comparison, and 13,349 DE-irlncRNA pairs were constructed. With single factor Cox regression analysis, prognosis-related DE-irlncRNA pairs were further selected. Then, lasso-penalized Cox analysis was conducted to narrow the number of lncRNA pairs to 14 over 1000 repetitions (Fig. 2a, b). Next, we constructed a risk assessment model that contained these 14 DEirlncRNA pairs using a risk score method and calculated the risk score of each sample. We performed receiver operating characteristic (ROC) analysis at 1 year (area under the curve (AUC) = 0.973) and identified the optimal cutoff point at 0.868 to divide 497 patients into high- and low-immune risk groups (Fig. 2c, d). We further investigated the prognostic value of the signature by plotting the distributions of the risk score and survival time (Fig. 2e, f). Low-risk TC patients exhibited a superior clinical outcome compared with high-risk TC patients. Then, 3-, 5-, and 10-year ROC curves were plotted and showed that the model had the ability to predict the survival outcome of TC patients with high accuracy and sensitivity (Fig. 2g).

Fig. 2
figure 2

Construction of a risk assessment model and confirmation of the signature. (a, b) Results of the lasso regression. (c) The cutoff value of the one-year ROC curve was used to separate the patients. (d) Fourteen DEirlncRNA pairs were identified to construct the signature. (e) Distribution of lncRNA model risk scores. (f) Survival status of TC patients in the subgroups. (g) ROC analysis of the risk scores for overall prognosis prediction. The 3-(red), 5-(green), and 10-year (blue) ROC curves of the model suggested that all AUC values were over 0.95. (h) PCA plot of the lncRNA model. (i) Kaplan‒Meier curve presenting survival in the high-risk and low-risk sets

Expression levels of fourteen DEirlncRNAs

We further validated the expression of 14 single lncRNAs of DEirlncRNA pairs in 12 matched tumor and peritumor samples using qRT‒PCR analysis. Compared with the peritumor controls, the expression levels of LINC00900, DCST1-AS1, HAGLROS, LINC02560, AL158206.1, SMIM25, AC090673.1, AC012038.2, and LINC01614 was significantly increased (P < 0.05), whereas ELN-AS1, AC093585.1, AC005237.1, and LBX2-AS1 were significantly decreased in PTC tissues (P < 0.05) (Additional file 3: Fig. S1), which was consistent with the bioinformatics results obtained from the TCGA dataset.

Clinical assessment and evaluation of the signature

Based on the validated optimal cutoff point, 45 TC patients with higher risk scores were included in the high-risk group, and 452 TC patients with lower risk scores were included in the low-risk group. The PCA results confirmed the reliability of this kind of classification (Fig. 2h). Kaplan‒Meier survival analysis showed that low-risk TC patients exhibited a significantly better prognosis than high-risk TC patients (P < 0.001, Fig. 2i). Furthermore, we assigned THCA patients into different subgroups according to clinical characteristics, including sex (male and female), N stage (N0 and N1), stage (stage I + II and stage III + IV) and T stage (T1 + T2 and T3 + T4). KM analysis also showed that high-risk patients had a worse overall survival than low-risk patients, which indicated that the signature prediction was excellent (P < 0.001, Fig. 3).

Fig. 3
figure 3

The prognostic value of the signature. THCA patients were assigned to different subgroups, including sex (male and female), N stage (N0 and N1), stage (stage I + II and stage III + IV) and T stage (T1 + T2 and T3 + T4)

Signature as an independent prognostic predictor

A series of chi-square tests were performed to explore the relationship between the risk of TC and common clinical features, including age, sex, TNM stage, and clinical stage. The strip chart (Fig. 4a) and scatter diagrams showed that age, N stage, clinical stage, and survival status (Additional file 4: Fig. S2) were significantly related to the risk of TC. Furthermore, we compared the AUCs of the 1-year ROC curves of the risk score to those of other clinical features and found that the risk score achieved a higher AUC value among these factors (Fig. 4b). Next, univariate and multivariate Cox regression analyses were performed to validate prognostic factors associated with TC. Univariate analysis showed that age (P < 0.001, hazard ratio (HR) = 1.153, 95% confidence interval (CI) [1.095 ~ 1.213]), stage (P < 0.001, HR = 2.929, 95% CI [1.896–4.527]), T stage (P = 0.005, HR = 2.384, 95% CI [1.291–4.405]), and risk score (P < 0.001, HR = 4.714, 95% CI [3.015–7.371]) were significantly associated with prognosis, whereas only age (P = 0.040, HR = 1.061, 95% CI [1.003–1.122]) and risk score (P = 0.005, HR = 3.409, 95% CI [1.750 − 6.641]) independently predicted the clinical outcome of TC patients by multivariate Cox regression analysis (Fig. 4c, d, Table 2).

Fig. 4
figure 4

The clinical value of the signature. (a) Strip chart, which was labeled as follows: < 0.001 = ***, < 0.01 = **, and < 0.05 = *. (b) Comparison of the 1-year ROC curves of the risk score with those of other clinical features showed the superiority of the risk score. (c, d) Univariate and multivariate Cox regression analyses were applied to identify prognostic factors associated with TC

Table 2 Independence of the signature for predicting the clinical outcomes of TC

Internal validation of the signature

We further validated the prognostic model in the training (Additional file 5: Fig. S3) and testing sets (Additional file 6: Fig. S4). The results from two Kaplan‒Meier survival analyses indicated that the risk score was closely linked to OS (P < 0.001) in the training and testing sets. Moreover, the PCA revealed a significant disparity in patients between the high- and low-risk groups. We next performed a time-dependent ROC analysis in the training cohort, and the results showed an area under the curve (AUC) of 0.929, 0.956, and 0.951 at 1 year, 3 years, and 5 years, respectively. In the testing cohort, the predictive efficacy (AUC values of 0.976, 0.942, and 0.981, respectively) was validated similarly. Univariate and multivariate Cox regression analyses showed that the risk score could independently predict the clinical outcome of TC patients in both sets. This result confirmed the overall accuracy and validity of the prognostic signature.

Correlation of the risk assessment model with immune cell infiltration in TC

The TME is essential for cancer prognosis and is composed of tumor cells, stromal cells, and infiltrating immune cells. We next investigated whether the signature was associated with the TME in TC patients. The proportion of tumor-infiltrating immune subsets was determined based on the CIBERSORT algorithm, and 21 kinds of immune cell profiles were detected in the TC samples (Fig. 5a, b). We found that the risk score had a significant relationship with the infiltration levels of immune cells. Correlation analysis was conducted using the Spearman correlation test, and the results were displayed in a lollipop diagram (Fig. 6a). The data are listed in Additional file 7: Fig. S5. The high-risk samples showed a negative association with tumor-infiltrating immune cells, such as CD8 + T cells, CD4 + T cells, neutrophils, B cells, NK cells and M1 macrophages, compared with the low-risk samples, whereas they were positively associated with Tregs, myeloid dendritic cells, monocytes, and cancer-associated fibroblasts, as analyzed by the Wilcoxon signed-rank test (Fig. 6b, Additional file 8: Table S3).

Fig. 5
figure 5

TIC profile in tumor samples and correlation analysis. (a) Proportions of the 21 kinds of TICs explored by the CIBERSORT algorithm. (b) Correlations between 21 immune cell components

Fig. 6
figure 6

Immune infiltration status in different risk groups. a Correlation between tumor immune infiltration and the immune-related lncRNA signature. b The P value comparing the risk score and tumor-infiltrating immune cells

Assessment of immune status in the THCA subtypes

We further evaluated the transcript levels of granzyme A (GZMA) and perforin (PRF1) in different risk groups. We found increased expression of GZMA in the low-risk group compared with the high-risk group (P < 0.01, Fig. 7a), whereas PRF1 (P > 0.05, Fig. 7b) showed no significant difference. The levels of human leucocyte antigen (HLA)-related gene expression were further studied, as shown in Fig. 7c. The results suggested that the low-risk group had higher expression levels of HLA genes than the high-risk group. Through the ESTIMATE algorithm, we explored the different characteristics of the immune microenvironment in the two risk groups. Compared with the high-risk group, the ESTIMATE score, immune score, and stromal score were all higher in the low-risk group (Fig. 7d). The results indicated that this phenomenon might be associated with the immunosuppressive microenvironment.

Fig. 7
figure 7

Differential immune features in the high-risk and low-risk subgroups. (a, b) The gene expression of GZMA and PRF1 in high- and low-risk groups. (c) The expression of HLA-related genes in high- and low-risk groups. (d) Stromal Score, Immune Score and ESTIMATE Score in the two risk subgroups (< 0.001 = ***, < 0.01 = **, and < 0.05 = *)

Correlation of the signature with genes related to ICIs in TC

Immune checkpoint blockade has revolutionized cancer treatment, and ICB therapy has already been applied in thyroid cancer patients. Therefore, we employed seven key ICI-related genes, PDCD1 (also called PD-1), CD274 (also called PD-L1), PDCD1LG2 (also called PD-L2), CTLA-4, LAG3, CD74, and IDO1, to investigate whether these ICI-related molecules were related to our new signature. The results indicated that the high-risk score group showed a negative correlation with the expression of PD-1 (P < 0.01, Fig. 8a), PD-L1 (P < 0.001, Fig. 8b), LAG3 (P < 0.01, Fig. 8c), CTLA-4 (P < 0.001, Fig. 8e), PD-L2 (P < 0.05, Fig. 8f), and CD74 (P < 0.001, Fig. 8g), whereas IDO1 (P > 0.05, Fig. 8d) showed no significant difference.

Fig. 8
figure 8

Correlation between ICI target genes and the immune-related lncRNA signature. (ag) The high-risk score group showed a negative correlation with the expression of PD-1 (a), PD-L1 (b), LAG3 (c), CTLA-4 (e), PD-L2 (f), and CD74 (g), whereas IDO1 (d) showed no significant difference. (< 0.001 = ***, < 0.01 = **, and < 0.05 = *)

The signature as an indicator in TC treatment

As a personalized medical treatment, small-molecule inhibitors have been widely applied in the clinical treatment of thyroid cancer. We next investigated the relationship between the signature and the half maximal inhibitory concentration (IC50) of some common inhibitors in the THCA dataset. As shown in Fig. 9, the high-risk samples were positively related to the IC50 of gefitinib (P < 0.01, Fig. 9a), sunitinib (P < 0.001, Fig. 9b), and tipifarnib (P < 0.001, Fig. 9c) but negatively related to axitinib (P < 0.001, Fig. 9e), AMG-706 (P < 0.001, Fig. 9f), and pazopanib (P < 0.01, Fig. 9g). In addition, lenvatinib and sorafenib showed no relationship with the signature (P > 0.05, Fig. 9d, h).

Fig. 9
figure 9

Correlation between common TKIs and the immune-related lncRNA signature. (ah) The high-risk samples were positively related to the IC50 of gefitinib (a), sunitinib (b), and tipifarnib (c) but negatively related to axitinib (e), AMG-706 (f), and pazopanib (g). Lenvatinib (d) and sorafenib (h) had no relationship with the signature. (< 0.001 = ***, < 0.01 = **, and < 0.05 = *)

Discussion

Thyroid cancer is the most common type of cancer of the endocrine system, and its incidence has increased almost threefold over the past decades [25]. Traditional clinical characteristics, such as TNM stage, can be used to predict the severity related to TC, but it is difficult to accurately estimate the risk of recurrence [26]. Thus, it is imperative to establish powerful tools that can be effectively applied to aid in the diagnosis, prognosis, and treatment of patients with TC.

Accumulating evidence reveals that since lncRNAs have high tissue and cancer specificity, they might play an active role in cancer initiation, development and progression. An increasing number of studies have shown that lncRNAs promote tumor cell proliferation, invasion, metastasis, and angiogenesis and can serve as an excellent tool to modulate therapeutic decisions in cancer. Increasing evidence suggests that lncRNAs are involved in TC tumorigenesis and progression as important regulatory factors [27]; thus, lncRNAs have attracted much attention as potential targets in the diagnosis and treatment evaluation of TC. For example, Liu et al. found that MALAT1 may have an oncogenic function in PTC and may thus be a potential diagnostic marker for PTC [28]. In our research, a few of the DEirlncRNAs in the model have already been revealed to play roles in various cancers, such as HAGLROS, LBX2-AS1, LINC00900, LINC01614, AC090673.1, and especially TC, while others were found to be related to TC for the first time. For instance, Li et al. demonstrated that LBX2-AS1 activated FSTL3 by binding to the TF RARα to hasten the proliferation, migration, and invasion of thyroid cancer [29]. HAGLROS was included in a prognostic model for TC, which was significantly correlated with TC recurrence [30]. Li et al. [31]constructed a prognostic model composed of seven lncRNAs, including LINC01614, AC090673.1, and LINC00900, which could serve as potential biomarkers for THCA prognosis. Nevertheless, it is necessary to validate whether this immune-related lncRNA model could be a helpful predictive indicator in TC. Many researchers are currently focused on constructing signatures with both coding genes and noncoding RNAs, which can assess the survival status of patients with malignant carcinoma [32, 33]. Unlike most traditional risk models, our newly constructed signature involved two-lncRNA pairwise comparisons and relative ranking based on gene expression entirely from the same TC patient. Although from different sequencing platforms, our prognostic model does not require gene expression data normalization. Previous studies have supported the effectiveness of this method [34, 35].

In the current study, we established an immune-related lncRNA model and evaluated its prognostic value as well as its correlation with immune cell infiltration, ICI-related genes, and TKIs in TC. First, we performed a differential coexpression analysis to identify DEirlncRNAs based on data from TCGA. LncRNA pairs were systematically identified through pairwise comparisons in the same sample without the need for data normalization. In addition, univariate analysis with Lasso regression analysis was performed on the pairs to validate the most suitable variables. Fourteen significant DEirlncRNA pairs with maximum prognostic values were determined with multiple repeats and random stimulation. Next, we used these pairs to develop the predictive risk score model. Then, we calculated not only the 3-, 5-, and 10-year AUC values of the ROC curve but also identified the optimal cutoff point of the 1-year ROC curve to separate TC patients into high- and low-risk groups. Furthermore, Kaplan‒Meier curves, time-dependent ROC curves, and Cox proportional hazards regression analysis showed that the model had independent predictive value for TC prognosis. Finally, we evaluated the relationship between this novel model and tumor-infiltrating immune cells, ICI-related molecules, and small-molecule inhibitor validity.

Immune infiltrates in the tumor microenvironment (TME) play a vital role in tumor development and progression and affect the clinical outcomes of cancer patients [36]. Dysfunction of the immune status in the TME contributes to the development and progression of cancer, and this is the basis of many immunotherapy studies. Moreover, tumor immunotherapy is now considered to have an important role in the elimination of cancer cells and sheds light on the mechanisms of cancer immune evasion, contributing to tumor outgrowth. Recent studies have suggested that lncRNAs play a central role in innate and adaptative cancer immunity regulation [37]. Immune-related lncRNA pairs as signatures are better at predicting prognosis than single lncRNAs. Therefore, it is necessary to explore more immune-related lncRNAs in tumors for future clinical practice. In this research, we carried out pairwise comparisons of a given set of immune-related lncRNAs and expression values. Thus, our prognostic signature could help address batch effects between different platforms and overcome the reprocessing and normalization of data.

Immune cell infiltration reflects the TME and reportedly impacts the outcome of TC progression. Immune-related lncRNAs are correlated with the development of TC. To explore the relationship between the prognostic model and immune-infiltrating cells, we applied seven commonly accepted methods, including TIMER [38, 39], CIBERSORT [40], XCELL [41], QUANTISEQ [42], MCP-counter, [43] EPIC [44], and CIBERSORT-ABS [45]. By integrating analyses, we found that the levels of Tregs, myeloid dendritic cells, monocytes and cancer-associated fibroblasts in the high-risk group were higher than those in the low-risk group, while the levels of neutrophils, M1 macrophages, NK cells, CD8 + , CD4 + T cells, and B cells were significantly negatively correlated with the risk of signature. These results were consistent with the findings of some previous experimental studies, which aimed to determine the correlation between each cell type and the aggressiveness of TC. For example, it has been reported that neutrophils play an antitumor role and can be beneficial to the prognosis of TC [46], which is consistent with the findings of our analysis. Our results also revealed that the abundance of Tregs was more associated with the high-risk group, which is similar to the findings of previous literature. In those studies, the levels of Tregs were higher in PTC than in multinodular goiter patients, and Tregs were consistently present in extraglandular invasion and lymph node metastasis [47, 48]. Monocytes have been observed to promote the occurrence and development of tumors, and their high density is closely related to thyroid tumor invasion and reduced survival, which was also confirmed in our study [49]. Dendritic cells (DCs) are the sentinel antigen-presenting cells (APCs) of the immune system [50], which suggests that the high level of DCs observed in the high-risk group could favor antigen presentation to T cells. Existing evidence shows that CD8 + , CD4 + T cells, and B cells play a protective role in the PTC tumor microenvironment (TME) [6]. In a wide immunohistochemical characterization of the immune network in patients with chronic lymphocytic thyroiditis concurrent with DTC, high CD8 + T lymphocyte infiltration was associated with improved disease-free survival. In addition, CD4 + T cells often differentiate into helper T cells, producing cytokines, and cooperating with CD8 + T cells to kill cancer cells [48]. A study indicated that B cells could promote carcinogenesis by inducing immunosuppression [51]. Natural killer (NK) cells are a family of innate immune cells, which play a central role in antiviral immunity and tumor immunosurveillance. A recent experimental study found that IL-12 immunotherapy could inhibit tumor growth and prolong survival by reactivating both CD8 + T and NK cells [52]. These observations can be further explored for a holistic understanding of the nuances of immune cell infiltration in the TC microenvironment.

The granzyme-perforin pathway is a primary method by which cytotoxic lymphocytes destroy cancer cells [53]. GZMA is a member of the granzyme family and is mainly secreted by cytotoxic cells, such as cytotoxic T cells and natural killer (NK) cells. PRF1 form membrane pores releasing granzymes leading to the cytolysis of the target cells [54]. We found that GZMA was highly expressed in the low-risk group, which is consistent with the higher NK cell percentage in the low-risk group (Additional file 7: Fig. S5). The ESTIMATE algorithm was used to explore the underlying mechanisms that cause PTC immune differences, and we found that the low-risk subset had higher stromal and immune scores. This finding is consistent with previous studies concluding that patients with high TME scores exhibit a stronger antitumor immune response, re likely to benefit more from immunotherapy, and survive longer [55]. Pioneering investigations have revealed that immunotherapy targeting immune checkpoints and human leukocyte antigen (HLA) provide great hope for the clinical treatment of human cancers [56]. The results showed that the expression levels of HLA-related genes were significantly higher in the low-risk group. Subsequently, we studied the interaction between the risk score and immune-related functions. The above findings further explained the reasons for the tumor-promoting status in the low-risk group and revealed the immunosuppressive microenvironment present in this group.

Recently, ICB immunotherapy has been viewed as a promising cancer therapeutic modality for malignant tumors. The identification of PD-L1 as an immunostat blockade has led to the development of several cancer immunotherapies. For RAI-refractory PTC patients, recent evidence has shown that overexpression of PD-L1 together with lymphocyte infiltration into the tumor TME is significantly associated with the effectiveness of ICB [57, 58]. Liotti et al. [59] observed a positive PD-1 effect on TC cell proliferation and migration through SHP-2 action on the Ras pathway. In this study, we found that patients in the low-risk subgroup had higher expression of common immune checkpoint molecules such as PD-1, PD-L1, LAG3, CTLA-4, PD-L2, and CD74, which are commonly expressed in human cancer. This phenomenon could indicate that low-risk patients might have a better, more beneficial response from ICB immunotherapy. Tyrosine kinase inhibitors (TKIs) are an innovative personalized strategy that target pro-oncogenic kinases, including EGFR, MET, PDGFR, VEGFR-1, VEGFR-2, RAF, FGFR, and RET. Our signature showed that the risk score was related to some of these inhibitors, such as gefitinib, sunitinib, and tipifarnib, indicating that this new model might be a novel method for assessing the efficacy of systemic therapy based on a genetic understanding of TC. In addition, some clinical trials have applied PD-1 and PD-L1 inhibitors in combination with TKIs, RAIs, or chemotherapy to manage and defeat deadly TC. Our signature may provide new insight to predict which patients are more suitable for these treatments, either alone or in combination.

To the best of our knowledge, a prognostic model based on irlncRNA pairs in TC has not been reported to date. Our predictive model is based on a 0-or-1 matrix and could be applied in an individualized manner while eliminating batch bias. In addition, our signature first combined DElncRNA pairs with ICB and TKI efficacy for analysis. Various additional methods were used to support the prognostic value and feasibility of this new model.

However, this lncRNA-based prognostic signature had several limitations. First, the establishment and validation of the model was based only on the TCGA database, which might lead to selection bias. To verify the predictive values of the risk assessment model, a larger external dataset of TC should be analyzed, preferably validated using GEO datasets. Unfortunately, no survival information of THCA could be obtained from the GEO cohort. Second, the research on the relationship between DElncRNA pairs and ICB and TKI efficacy was based on the inference of several algorithms and has not been experimentally verified. In addition, the calculation formulas of this prognostic signature may be too complex for clinical application. However, more experimental studies are still needed to validate our observations.

Conclusion

In conclusion, we constructed an immune-related lncRNA pair model without considering the technical bias of different platforms. As such, our novel model could be used without the need to eliminate batch effects and could serve as an independent single-sample estimate of the survival risk subgroup of TC patients. The model may provide new possibilities for translation to clinical practice for TC patients and help in distinguishing those who could benefit from ICB immunotherapy and TKI therapy.

Availability of data and materials

The generated and analyzed datasets of the current research are available in TCGA database (http://cancergenome.nih.gov/), and the immune-related genes were obtained from the ImmPort database (https://immport.niaid.nih.gov).

Abbreviations

TC:

Thyroid cancer

DEirlncRNA:

Differentially expressed immune-related lncRNA

ICIs:

Immune checkpoint inhibitors

DTC:

Differentiated thyroid cancer

ICB:

Immune checkpoint blockade

TKIs:

Tyrosine kinase inhibitors

lncRNAs:

Long noncoding ribonucleic acids

TCGA:

The cancer genome Atlas

LASSO:

Least absolute shrink age and selection operator

qRT-PCR:

Quantitative real-time PCR

PCA:

Principal component analysis

ROC:

Receiver operating characteristic

TME:

Tumor microenvironment

References

  1. Vaccarella S, Dal Maso L, Laversanne M, Bray F, Plummer M, Franceschi S. The impact of diagnostic changes on the rise in thyroid cancer incidence: a population-based study in selected high-resource countries. Thyroid. 2015;25(10):1127–36.

    PubMed  Article  Google Scholar 

  2. Rahib L, Smith B, Aizenberg R, Rosenzweig A, Fleshman J, Matrisian L. Projecting cancer incidence and deaths to 2030: the unexpected burden of thyroid, liver, and pancreas cancers in the United States. Can Res. 2014;74(11):2913–21.

    CAS  Article  Google Scholar 

  3. Chowdhury S, Veyhl J, Jessa F, Polyakova O, Alenzi A, MacMillan C, Ralhan R, Walfish P. Programmed death-ligand 1 overexpression is a prognostic marker for aggressive papillary thyroid cancer and its variants. Oncotarget. 2016;7(22):32318–28.

    PubMed  PubMed Central  Article  Google Scholar 

  4. Schreiber R, Old L, Smyth M. Cancer immunoediting: integrating immunity’s roles in cancer suppression and promotion. Science (New York, NY). 2011;331(6024):1565–70.

    CAS  Article  Google Scholar 

  5. Grivennikov S, Greten F, Karin M. Immunity, inflammation, and cancer. Cell. 2010;140(6):883–99.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  6. Galdiero MR, Varricchi G, Marone G. The immune network in thyroid cancer. Oncoimmunology. 2016;5(6): e1168556.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  7. Chen DS, Mellman I. Oncology meets immunology: the cancer-immunity cycle. Immunity. 2013;39(1):1–10.

    PubMed  Article  CAS  Google Scholar 

  8. Bastman J, Serracino H, Zhu Y, Koenig M, Mateescu V, Sams S, Davies K, Raeburn C, McIntyre R, Haugen B, et al. Tumor-infiltrating T cells and the PD-1 checkpoint pathway in advanced differentiated and anaplastic thyroid cancer. J Clin Endocrinol Metab. 2016;101(7):2863–73.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. Mehnert J, Varga A, Brose M, Aggarwal R, Lin C, Prawira A, de Braud F, Tamura K, Doi T, Piha-Paul S, et al. Safety and antitumor activity of the anti-PD-1 antibody pembrolizumab in patients with advanced, PD-L1-positive papillary or follicular thyroid cancer. BMC Cancer. 2019;19(1):196.

    PubMed  PubMed Central  Article  Google Scholar 

  10. Cabanillas M, Habra M. Lenvatinib: role in thyroid cancer and other solid tumors. Cancer Treat Rev. 2016;42:47–55.

    CAS  PubMed  Article  Google Scholar 

  11. Brose MS, Nutting CM, Jarzab B, Elisei R, Siena S, Bastholt L, de la Fouchardiere C, Pacini F, Paschke R, Shong YK, et al. Sorafenib in radioactive iodine-refractory, locally advanced or metastatic differentiated thyroid cancer: a randomised, double-blind, phase 3 trial. Lancet. 2014;384(9940):319–28.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, Barrette TR, Prensner JR, Evans JR, Zhao S, et al. The landscape of long noncoding RNAs in the human transcriptome. Nat Genet. 2015;47(3):199–208.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  13. Sahu A, Singhal U, Chinnaiyan A. Long noncoding RNAs in cancer: from function to translation. Trends Cancer. 2015;1(2):93–109.

    PubMed  PubMed Central  Article  Google Scholar 

  14. Fang Y, Fullwood M. Roles, functions, and mechanisms of long non-coding RNAs in cancer. Genom Proteom Bioinform. 2016;14(1):42–54.

    Article  Google Scholar 

  15. Mercer T, Dinger M, Mattick J. Long non-coding RNAs: insights into functions. Nat Rev Genet. 2009;10(3):155–9.

    CAS  PubMed  Article  Google Scholar 

  16. Chen Y, Satpathy A, Chang H. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol. 2017;18(9):962–72.

    CAS  PubMed  Article  Google Scholar 

  17. Bhan A, Soleimani M, Mandal SS. Long noncoding RNA and cancer: a new paradigm. Cancer Res. 2017;77(15):3965–81.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. Peng X, Zhang K, Ma L, Xu J, Chang W. The role of long non-coding RNAs in thyroid cancer. Front Oncol. 2020;10:941.

    PubMed  PubMed Central  Article  Google Scholar 

  19. Feng J, Zhou Q, Yi H, Ma S, Li D, Xu Y, Wang J, Yin S. A novel lncRNA n384546 promotes thyroid papillary cancer progression and metastasis by acting as a competing endogenous RNA of miR-145-5p to regulate AKT3. Cell Death Dis. 2019;10(6):433.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  20. Huang Z, Xiao C, Zhang F, Zhou Z, Yu L, Ye C, Huang W, Li N. A novel framework to predict breast cancer prognosis using immune-associated LncRNAs. Front Genet. 2020;11: 634195.

    CAS  PubMed  Article  Google Scholar 

  21. Xu Q, Wang Y, Huang W. Identification of immune-related lncRNA signature for predicting immune checkpoint blockade and prognosis in hepatocellular carcinoma. Int Immunopharmacol. 2021;92: 107333.

    CAS  PubMed  Article  Google Scholar 

  22. Ritchie M, Phipson B, Wu D, Hu Y, Law C, Shi W, Smyth G. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucl Acids Res. 2015;43(7): e47.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  23. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22.

    PubMed  PubMed Central  Article  Google Scholar 

  24. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    PubMed  Article  CAS  Google Scholar 

  25. Fagin JA, Wells SA Jr. Biologic and clinical perspectives on thyroid cancer. N Engl J Med. 2016;375(23):2307.

    PubMed  Article  Google Scholar 

  26. Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, Pacini F, Randolph GW, Sawka AM, Schlumberger M, et al. 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;26(1):1–133.

    PubMed  PubMed Central  Article  Google Scholar 

  27. Kentwell J, Gundara J, Sidhu S. Noncoding RNAs in endocrine malignancy. Oncologist. 2014;19(5):483–91.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. Liu J, Dong H, Yang Y, Qian Y, Liu J, Li Z, Guan H, Chen Z, Li C, Zhang K, et al. Upregulation of long noncoding RNA MALAT1 in papillary thyroid cancer and its diagnostic value. Future Oncol. 2018;14(29):3015–22.

    CAS  PubMed  Article  Google Scholar 

  29. Li J, Shen J, Qin L, Lu D, Ding E. LBX2-AS1 activates FSTL3 by binding to transcription factor RARα to foster proliferation, migration, and invasion of thyroid cancer. Front Genet. 2021;12: 765033.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. Guo K, Chen L, Wang Y, Qian K, Zheng X, Sun W, Sun T, Wu Y, Wang Z. Long noncoding RNA RP11–547D24.1 regulates proliferation and migration in papillary thyroid carcinoma: Identification and validation of a novel long noncoding RNA through integrated analysis of TCGA database. Cancer Med. 2019;8(6):3105–19.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. Li Z, Wang H, Deng X, Zhang J, Wang L, Tang W, Youa W, Nian W, Construction and validation of an immune-related lncRNA prognosis model for thyroid cancer. Comb Chem High Throughput Screen 2022;25. https://doi.org/10.2174/1386207325666220223123917.

  32. Zhu X, Tian X, Yu C, Shen C, Yan T, Hong J, Wang Z, Fang JY, Chen H. A long non-coding RNA signature to improve prognosis prediction of gastric cancer. Mol Cancer. 2016;15(1):60.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  33. Gu J, Zhang X, Miao R, Xiang X, Fu Y, Zhang J, Liu C, Qu K. Six-long non-coding RNA signature predicts recurrence-free survival in hepatocellular carcinoma. World J Gastroenterol. 2019;25(2):220–32.

    PubMed  PubMed Central  Article  Google Scholar 

  34. Peng PL, Zhou XY, Yi GD, Chen PF, Wang F, Dong WG. Identification of a novel gene pairs signature in the prognosis of gastric cancer. Cancer Med. 2018;7(2):344–50.

    CAS  PubMed  Article  Google Scholar 

  35. Li B, Cui Y, Diehn M, Li R. Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol. 2017;3(11):1529–37.

    PubMed  PubMed Central  Article  Google Scholar 

  36. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol. 2020;17(8):807–21.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  37. Denaro N, Merlano M, Lo Nigro C. Long noncoding RNAs as regulators of cancer immunity. Mol Oncol. 2019;13(1):61–73.

    CAS  PubMed  Article  Google Scholar 

  38. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucl Acids Res. 2020;48(W1):W509–14.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  41. Aran D. Cell-type enrichment analysis of bulk transcriptomes using xCell. Methods Mol Biol. 2020;2120:263–76.

    CAS  PubMed  Article  Google Scholar 

  42. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, Krogsdam A, Loncova Z, Posch W, Wilflingseder D, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 2019;11(1):34.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  43. Dienstmann R, Villacampa G, Sveen A, Mason MJ, Niedzwiecki D, Nesbakken A, Moreno V, Warren RS, Lothe RA, Guinney J. Relative contribution of clinicopathological variables, genomic markers, transcriptomic subtyping and microenvironment features for outcome prediction in stage II/III colorectal cancer. Ann Oncol. 2019;30(10):1622–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. Racle J, de Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. 2017. https://doi.org/10.7554/eLife.26476.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Tamminga M, Hiltermann TJN, Schuuring E, Timens W, Fehrmann RS, Groen HJ. Immune microenvironment composition in non-small cell lung cancer and its association with survival. Clin Transl Immunol. 2020;9(6): e1142.

    CAS  Article  Google Scholar 

  46. Mao M, Huang R, Zheng J, Liang H, Huang W, Liu J, Li J. OGDHL closely associates with tumor microenvironment and can serve as a prognostic biomarker for papillary thyroid cancer. Cancer Med. 2021;10(2):728–36.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. Liu Y, Yun X, Gao M, Yu Y, Li X. Analysis of regulatory T cells frequency in peripheral blood and tumor tissues in papillary thyroid carcinoma with and without Hashimoto’s thyroiditis. Clin Transl Oncol. 2015;17(4):274–80.

    PubMed  Article  CAS  Google Scholar 

  48. French JD, Weber ZJ, Fretwell DL, Said S, Klopper JP, Haugen BR. Tumor-associated lymphocytes and increased FoxP3+ regulatory T cell frequency correlate with more aggressive papillary thyroid cancer. J Clin Endocrinol Metab. 2010;95(5):2325–33.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  49. Zhou C, Duan D, Liu S. Predictive value of a prognostic model based on lymphocyte-to-monocyte ratio before radioiodine therapy for recurrence of papillary thyroid carcinoma. Technol Cancer Res Treat. 2021;20:15330338211027910.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. Dudek AM, Martin S, Garg AD, Agostinis P. Immature, semi-mature, and fully mature dendritic cells: toward a DC-cancer cells interface that augments anticancer immunity. Front Immunol. 2013;4:438.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  51. Wennhold K, Shimabukuro-Vornhagen A, von Bergwelt-Baildon M. B cell-based cancer immunotherapy. Transfus Med Hemother. 2019;46(1):36–46.

    PubMed  PubMed Central  Article  Google Scholar 

  52. Parhar RS, Zou M, Al-Mohanna FA, Baitei EY, Assiri AM, Meyer BF, Shi Y. IL-12 immunotherapy of Braf(V600E)-induced papillary thyroid cancer in a mouse model. Lab Invest. 2016;96(1):89–97.

    CAS  PubMed  Article  Google Scholar 

  53. Voskoboinik I, Whisstock JC, Trapani JA. Perforin and granzymes: function, dysfunction and human pathology. Nat Rev Immunol. 2015;15(6):388–400.

    CAS  PubMed  Article  Google Scholar 

  54. Voskoboinik I, Dunstone MA, Baran K, Whisstock JC, Trapani JA. Perforin: structure, function, and role in human immunopathology. Immunol Rev. 2010;235(1):35–54.

    CAS  PubMed  Article  Google Scholar 

  55. Zeng D, Zhou R, Yu Y, Luo Y, Zhang J, Sun H, Bin J, Liao Y, Rao J, Zhang Y, et al. Gene expression profiles for a prognostic immunoscore in gastric cancer. Br J Surg. 2018;105(10):1338–48.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  56. Olson B, Li Y, Lin Y, Liu ET, Patnaik A. Mouse models for cancer immunotherapy research. Cancer Discov. 2018;8(11):1358–65.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  57. Cunha L, Marcello M, Morari E, Nonogaki S, Conte F, Gerhard R, Soares F, Vassallo J, Ward L. Differentiated thyroid carcinomas may elude the immune system by B7H1 upregulation. Endocr Relat Cancer. 2013;20(1):103–10.

    CAS  PubMed  Article  Google Scholar 

  58. Shi R, Qu N, Luo T, Xiang J, Liao T, Sun G, Wang Y, Wang Y, Huang C, Ji Q. Programmed death-ligand 1 expression in papillary thyroid cancer and its correlation with clinicopathologic factors and recurrence. Thyroid: Off J Am Thyroid Assoc. 2017;27(4):537–45.

    CAS  Article  Google Scholar 

  59. Liotti F, Kumar N, Prevete N, Marotta M, Sorriento D, Ieranò C, Ronchi A, Marino FZ, Moretti S, Colella R, et al. PD-1 blockade delays tumor growth by inhibiting an intrinsic SHP2/Ras/MAPK signalling in thyroid cancer cells. J Exp Clin Cancer Res. 2021;40(1):22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We acknowledge the TCGA database for providing their platforms and the contributors for uploading their meaningful datasets. The authors hereby confirm that neither the manuscript nor any part of it has been published or is being considered for publication elsewhere. We acknowledge that all authors participated sufficiently in the work and take public responsibility for its content.

Funding

This research was funded by The National Natural Science Foundation of China, Grant Number 81800700) and The Clinical Research Fund of Chinese Medical Association (Grant No. 15010010589).

Author information

Authors and Affiliations

Authors

Contributions

BS is the first author of this study. WT is the corresponding author supervising this work. BS designed the study and wrote the first manuscript; BG performed the statistical analysis; WT reviewed the manuscript; LT, FZ and ZL wrote sections of the manuscript. TL performed the project administration. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Weiping Teng.

Ethics declarations

Ethics approval and consent to participate

Informed consent was obtained from all individual participants included in the study. All the individual participants consented to the publication of medical data (including figures from diagnostic imaging results and from histological examination results). The research was approved by the institutional review board at Medical Science Research Ethics Committees of the First Affiliated Hospital of China Medical University (protocol number 275 from 3th of July 2020) and the 1964 Helsinki declaration and its later amendments or comparable ethical standards.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

A total of 1148 immune-related lncRNAs were validated by co-expression analysis (submitted as a separate EXCEL file).

Additional file 2: Table S2.

Premier Sequences for qRT‒PCR Analysis.

Additional file 3: Figure S1.

The expression levels of 14 single lncRNAs of DEirlncRNA pairs quantified using qRT-PCR analysis in 12 paired thyroid cancer tissues and no-tumorous samples.

Additional file 4: Figure S2.

Boxplots indicated that (a) age, (b) survival status, (c) N stage, and (d) clinical stage were related to the risk score.

Additional file 5: Figure S3.

Confirmation of the Signature in the Training Set.

Additional file 6: Figure S4.

Confirmation of the Signature in the Testing Set.

Additional file 7: Figure S5.

The Specific Correlations between Tumor Infiltrating Immune Cells and the Signature.

Additional file 8: Table S3.

The P value of Comparing Risk Sore and Tumor Infiltrating Immune Cells.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Song, B., Tian, L., Zhang, F. et al. A novel signature to predict thyroid cancer prognosis and immune landscape using immune-related LncRNA pairs. BMC Med Genomics 15, 183 (2022). https://doi.org/10.1186/s12920-022-01332-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-022-01332-7

Keywords

  • TCGA
  • lncRNA pairs
  • Immune infiltration
  • Thyroid carcinoma
  • Immune checkpoint inhibitors