Identification of necroptosis-related DEGs
A total of 50 differentially expressed genes between 510 tumors and 58 adjacent nontumor tissues were identified. Out of 67 necroptosis-related genes, 21 were up-regulated and 29 were downregulated in tumors. A heatmap was drawn for showing the 50 genes expression levels (Figure 1A). A PPI network was constructed to show the relationships among all necroptosis-related genes (Figure 1B). GO and KEGG analysis were used to investigate the potential roles of these DEGs, where GO enrichment analysis revealed that these necroptosis-related DEGs are mostly engaged in programmed necrotic cell death, transcription factor activity regulation, and cytokine response. According to KEGG analysis, these genes were primarily involved in cancer and apoptotic pathways (Figure 1C).
Tumor classification based on DEGs
To investigate the association between necroptosis-related DEGs and PTC subtypes, consensus clustering analysis on the 50 DEGs using the “ConsensusClusterPlus” program was performed. The letter “k” denoted the number of clusters. When k was adjusted to 2, intragroup correlations were highest while the intergroup correlations were lowest, showing that the 498 PTC patients were classified into two clusters (Figure 2A), C1 (N = 361) and C2 (N = 137). Additionally, the principal component analysis demonstrated that these patients can be classified into two distinct groups (Figure 2B). However, survival analysis revealed no statistically significant difference in OS between the two clusters (Figure 2C). Following that, a heatmap was created to depict the association between clusters, gene expression profiles, and clinical characteristics (Figure 2D), however, it was discovered that there are minimal clinical distinctions between the two clusters.
Development of a prognostic gene model
The models were constructed using the train set. To begin, a univariate Cox regression analysis was used to find survival-related genes. The results indicated that 16 genes in the training set satisfying a threshold with a p < 0.2 were associated with survival (Figure 2E). More precisely, 12 genes (MAPK8, IPMK, CYLD, AXL, ID1, CDKN2A, BCL2, ATRX, SPATA2, KLF9, BNIP3, and BCL2L11) were detrimental with HR > 1, whereas four genes (FADD, FAS, TNFRSF1B, and MPG) were protective with HR < 1. The overfitting was reduced by using LASSO Cox regression on the optimum λ value and identified 13 genes (FADD, FAS, IPMK, TNFRSF1B, CYLD, AXL, ID1, CDKN2A, MPG, SPATA2, KLF9, BNIP3, and BCL2L11) that were substantially linked with prognosis (Figure 2F & G). Further study of these genes using multivariate Cox regression analysis resulted in the construction of a prognostic signature consisting of seven necroptosis-related genes. For each sample, a risk score was calculated using the following formula: Risk score = FAS*(-1.1721) + IPMK*(2.5404) +TNFRSF1B*(0.9750) + AXL*(1.5358) + CDKN2A*(1.6495) + SPATA2*(4.1484) + KLF9*(2.0030). The patients were classified in the training set into high- and low-risk groups based on their risk score median value and then ranked and assessed their risk score distributions (Figure 3A). Additionally, the survival status of each sample in the training set was visualized using dot plots and discovered that all deceased patients were concentrated in the high-risk group (Figure 3B). A heatmap was created to depict the differential expression of the seven prognostic genes between the two groups (Figure 3C). Survival analysis revealed that patients with high-risk scores had a significantly worse outcome than those with low-risk scores (p = 0.002) (Figure 3J). ROC analysis was utilized to determine the sensitivity and specificity of created prognostic signature, and the areas under the curve (AUC) were found to be 0.834 at one year, 0.938 at three years, and 0.960 at five years respectively (Figure 3M). The PCA revealed that patients with varying risks were dispersed in two directions (Figure 3P).
Following that, the predictive ability of prognostic signature in the test and total sets was validated. The test and total sets were separated into high- and low-risk groups, respectively, based on the median risk score in the training set. The risk score distribution, survival status, and expression of seven necroptosis-related genes are provided below for each patient in the test set and total set (Figure 3D,E,F; Figure 3G, H, I). In both the test and total sets, Kaplan-Meier survival curves indicated that high-risk patients had a shorter OS than low-risk patients (Figure 3K; 3L). In the test set, the one-year AUC was found to be 0.703, the three-year AUC was 0.656, and the five-year AUC was 0.897 in the test set respectively (Figure 3N). In the total set, the AUC for the one-year was found to be 0.795, the AUC for the three years was 0.822, and the AUC for the five years was 0.902. (Figure 3O). PCA verified the separation of patients into two clusters according to their risk status (Figure 3Q; 3R).
Independent prognostic value of the risk signature
The univariate Cox regression analysis revealed that the factors affecting the prognosis were age, clinical stage, and risk grouping (Figure 4A). After adjusting for other confounding factors, it was discovered that age and risk grouping were independent predictors of OS in patients with PTC by multivariate Cox regression analysis (Figure 4B). Furthermore, the ROC curves of risk grouping and other clinicopathological variables showed that risk grouping had a higher prediction level, and its prediction ability increased with time (Figure 4C-E). Furthermore, survival analysis revealed that individuals with high-risk ratings had lower overall survival in all clinical subgroups (Figure S2A-L).
Functional analyses based on the risk signature
A heatmap was created to depict the link between predictive gene expression and clinical features (Figure 5A). Then, using FDR < 0.05 and |log2FC| ≥ 0.585, 268 differentially expressed genes were identified in high- and low-risk groups. Based on DEGs, GO and KEGG analysis was performed to further elucidate the putative biological activities and pathways associated with the risk signature. According to findings, these DEGs were primarily enriched in cell adhesion, immune-related activities, and several infections or immune-related pathways (Figure 5B, C). The GSEA was used to examine the transcript messages of PTC patients who were grouped by risk score into high- and low-risk categories. According to findings, most NRGs prognostic signatures regulated immunological and tumor-related pathways. In the low-risk group, biological pathways such as interferon-alpha and gamma response, allograft rejection, P53 pathway, inflammatory response, and IL6-JAK-STAT3 signaling were enriched. (Figure 5D).
Comparison of the immune activity among subgroups
The ssGSEA was used to compare the enrichment scores of 16 different types of immune cells and the activity of 13 immune-related functions in low- and high-risk groups. The findings revealed that, except for CD8+ T cells, which did not differ substantially between the two groups, other immune cells infiltrate at a higher rate in the low-risk subgroup (Figure 6A). Furthermore, the great majority of immune-related pathways were activated more in the low-risk group than in the high-risk group (Figure 6B). Each sample's immune-score, stromal-score, and ESTIMATES-score were calculated to determine the estimated amount of stromal and immune cells in the tumor sample. It was observed that the low-risk group contained more immune cells than the high-risk group, but the difference in stromal cell composition between the two groups was not statistically significant. Furthermore, ESTIMATES-score was greater in the low-risk group, indicating reduced tumor purity (Figure 6C). The potential differences in immune checkpoint expression between the two groups were investigated, and the results revealed that most immune checkpoints were more expressed in the low-risk group (Figure 6D). Similar outcomes were reached when analyzing the immunological state of GEO cohort. The low-risk group had more immune infiltration and immune checkpoint expression than the high-risk group (Figure 7).
The expression levels of seven prognostic genes
Following that, bioinformatics analysis was performed of mRNA expression data from GEO database and qRT-PCR experiments to validate mRNA expression profiles using 48 paired tumor and normal tissues from PTC patients at China Medical University, in addition to TCGA database, to validate the expression profile of seven prognostic genes in both tumor and normal tissues. The qRT-PCR results revealed that AXL (p = 0.0308), FAS (p = 0.0134), TNFRSF1B (p = 0.0250), and CDKN2A (p = 0.0224) mRNA expression were considerably greater in tumor tissues. In contrast, tumor samples had lower levels of IPMK (p = 0.0010), SPATA2 (p = 0.0071), and KLF9 (p = 0.0002). The expression of these genes was consistent with TCGA and GEO datasets, which adds to the authenticity and confidence of our developed signature (Figure 8).