Necroptosis-Related DEGs in OC
A total of 379 patients with OC from the TCGA and 88 normal ovarian tissues from the GTEx were collected, and 139 necroptosis-related genes were downloaded from the KEGG database. The expression of 133 necroptosis-related genes was extracted from the expression matrix of the TCGA and GTEx transcriptome. Through differential gene expression analysis, we identified 120 necroptosis-related DEGs detected in normal ovarian samples and OC samples and visualized them by using heatmaps. As shown in Fig. 2A, some necroptosis genes were up-regulated in the OC cases, including VDAC3, TNF, and STAT1. Other genes were downregulated in the tumor samples, such as JAK1, PYG1, and BCL2.
In addition, the PPI network of 10 hub genes, namely, STAT1, JAK1, TNF, TNFRSF1A, TRADD, FADD, CASP8, TRAF2, RIPK1, and RIPK3, was shown in Fig. 2B. As for the genetic alterations of 10 hub genes, we found that queried genes were altered in 127 (41%) queried samples, mainly including changes in deep deletion or amplification, and missense mutation was the major form of the gene mutations. Finally, 10 genes are shown in Fig. 2C. TNFRSF1 and RIPK1 were the most frequently mutated genes, with the same genetic alteration rate of 10%.
Construction of a Prognostic Risk Model
Given the effect of differential genes on OS of patients with OC, univariate regression survival analysis was used in identifying necroptosis-related genes associated with prognosis. Six genes: IFNB1, JAK1, PYGB, STAT1, STAT4, and ZBP1, were finally screened. Notably, IFNB1, STAT1, STAT4, and ZBP1 belonged to low-risk genes (hazard ratio [HR] < 1, P < 0.05), whereas JAK1 and PYGB belonged to high-risk genes (HR > 1, P < 0.05), which are shown in Fig. 3A (low-risk genes were marked in blue, and high-risk genes were marked in red). According to the expression of the six genes mentioned above, multivariate regression analysis was used in screening candidate prognostic genes. Eventually, three risk genes (i.e., JAK1, PYGB, and STAT1) were used in constructing the model, and the risk scores for each patient were calculated using their gene expression levels and regression coefficients.
First, using the median risk score from the TCGA cohort as a benchmark, we divided the OC cases into high- and low-risk groups. As shown in a risk plot, the high-risk group is shown in yellow, whereas the low-risk group is shown in sky blue (see Fig. 3B). Then, we profiled the survival status of each patient in the TCGA cohort. Patients who died were marked in red, and those who survived were marked in blue. The results showed that the risk of patients increased with the number of deaths (Fig. 3C). Subsequently, PCA and t-SNE analysis (Fig. 3D and 3E) were used in dimension reduction for the visualization of the three genes involved in modeling in a two-dimensional space. Patients from different risk groups were separated in both directions. The genes above can be used in distinguishing between patients in different risk groups.
We drew a KM curve to identify differences in prognosis between low- and high-risk groups. Obviously, the prognostic distinctions were evident. The high-risk group had a poor prognosis, whereas the low-risk group had a favorable prognosis in the TCGA patient cohort (Fig. 4A). Furthermore, our study focused on measuring the area under the curve (AUC) with ROC to validate the accuracy of the model prediction. Specifically, this model can relatively predict the OS of patients with OC at 1, 3, and 5 years, achieving AUC values of 0.573, 0.614, and 0.656 on the test set from the TCGA cohort (Fig. 4B), showing that this model can accurately predict the 5-year OS of patients with OC to some extent.
In addition to the content mentioned above, independent prognostic analysis was needed in assessing the acceptability of the risk score in the prediction of prognosis independent of any other clinical prognostic features, such as age, tumor grade, and stage. Univariate Cox analysis predicted that in the TCGA cohort, age, stage, residual tumor after tumor reduction surgery and risk score were closely related to the OS of patients with OC (P < 0.05, Fig. 4C). On the other hand, when the interplay among various clinicopathological variables was considered, the multivariate Cox analysis was applied, and a similar conclusion was reached—age, residual tumor and risk score were still closely associated with OS (P < 0.01, Fig. 4D).
Validation of the Prognostic Risk Model
To test the predictive power of gene signature, we tested it using the TCGA cohort as an internal test cohort and validated it on another independent GEO patient cohort. Thus, we extracted genes from the three-gene signature model and constructed a multivariate Cox regression model in the GEO database. Likewise, on the basis of the median risk score from the GEO cohort, the OC cases were classified into high- and low-risk groups. Then, we calculated the survival rates of patients according to KM survival in different risk groups. The high-risk group had significantly worse survival outcomes than the low-risk group (Fig. 4E). Furthermore, the ROCs of patients with OC from the GEO dataset (GSE19829) reconfirmed the accuracy of our prediction model to some extent, and the AUCs were 0.73, 0.66, and 0.597 (Fig. 4F). Finally, in the validation of GSE19829, we found that only the risk score can be validated, which offers prognostic information regarding the survival outcomes of patients with OC according to the univariate Cox analysis, compared with known prognostic clinical factors (Fig. 4G).
Construction of a Nomogram
Nomogram was constructed based on various prognostic indicators and risk scores to predict 1-, 3-, and 5-year survival in OC patients, respectively (Fig. 4H). The calibration curve (Fig. 4I) demonstrated the excellent predictive power of the Nomogram.
Identification of High-Risk or Low-Risk genes
As shown in the heatmap (Fig. 5), the OC cases from the TCGA database were grouped into different risk groups. The high expression of each gene involved in model construction was indicated in blue, whereas low expression was given in red. From left to right, the expression of JAK1 and PYGB increased with the risk of patients, whereas STAT1 expression decreased with increasing risk of patients. These results demonstrated JAK1 and PYGB were high-risk genes and STAT1 belonged to a low-risk gene. Unfortunately, no statistical differences in the clinical traits of OC cases (e.g., grade and stage) were found between the low- and high-risk groups.
Functional Enrichment Analysis of DEGs in High- and Low-Risk Groups
Patients with OC in the TCGA dataset were divided into two risk groups, namely low- and high-risk groups. Subsequently, when FDR < 0.05 and log2FC > 0.585, we screened 105 significant necroptosis-related DEGs between the two groups and performed downstream functional enrichment analysis. As shown in Fig. 6A and 6B, Gene Ontology (GO) analysis of the candidate genes revealed strong enrichment of functions associated with type 1 interferon, response to virus, and defense response to symbiont. Only four enriched terms were identified in the subsequent KEGG pathway analyses, including the interactions of viral proteins with cytokines, hepatitis C, influenza A, and the chemokine signaling pathway (Supplementary Fig. 2A and 2B).
Immune-Related Analysis of DEGs in High- and Low-Risk Groups
Notably, apparent differences in immune cell types or immune functions between different risk groups are shown in Fig. 6C and 6D. For example, aDCs, CD8+T cells, Tfh, APC-co-inhibition, inflammation-promoting, and type-I-IFN-response, were upregulated in the low-risk groups, whereas type-II-IFN-response were significantly upregulated in the high-risk groups (***P < 0.001; **P < 0.01; *P < 0.05).
Overexpression of Each Gene Involved in Model Construction in Tumor Tissues
As shown in the previous section, PYGB, JAK1, and STAT1 played important roles in the construction of our necroptosis model. To further demonstrate the level of PYGB in OC, the HPA database was used in confirming the expression of PYGB. As shown in Fig. 7A, the expression of PYGB was enhanced in the OC tissues compared with those in the normal ovarian tissues. Subsequently, the relationship between PYGB expression levels in OC and the prognosis was analyzed. Among the 373 OC cases, 258 cases were highly expressed with a 5-year survival rate of 29%, and 115 cases showed low expression with a 5-year survival rate of 40%. These results indicated that a high PYGB expression tends to be inversely linked to the prognosis of OC patients (Fig. 7B, P = 0.0084). We successively showed the expression levels of JAK1 and STAT1 in normal and tumor tissues. Both genes showed higher expression levels in OC than in normal tissues (Fig. 7C and 7E). As for the relationship between JAK1 and prognosis, patients with OC and high JAK1 expression had significantly reduced 5-year OS rates relative to patients with low JAK1 expression (Fig. 7D, P = 0.0079). By contrast, the 5-year survival rates of STAT1 in the high and low expression groups were 45% and 24%, respectively, indicating that the overexpression of STAT1 is associated with improved prognosis (Fig. 7F, P = 0.00017). Thus, we speculated that JAK1, similar to PYGB, is an unsuitable prognostic marker, whereas STAT1 expression is a favorable biomarker for OC.
All these data suggested that each gene involved in model construction plays a crucial role in the occurrence and progress of OC.
Tumor Typing and Typing Survival Analysis
Through consensus clustering analysis, we observed that k = 2 is the most suitable option for subtyping patients with OC (Fig. 8A). Thus, when clustering stability was considered, patients with OC were separated into two subtypes, namely clusters 1 and 2. The KM curve indicated striking differences (P < 0.05) between the two clusters. Specifically, the survival rate of the patients in cluster 1 was better than that of the patients in cluster 2 (Fig. 8B).
Drug sensitivity analysis
To screen out potential drugs for the treatment of ovarian cancer, we further analyzed and predicted the differences in the sensitivity of different chemotherapy drugs in high and low-risk groups (Fig. 9). The results showed the half-maximal inhibitory concentration (IC50) values of cisplatin and docetaxel in the high-risk group were higher than those in the low-risk group, suggesting that patients in the low-risk group were more sensitive to these drugs. However, no significant differences between different risk groups in the IC50 values of other drugs were found.