Characterization of immune and hypoxia genes
To discover the hub genes which could regulate both immunity and hypoxia process, we screened 31 overlapped genes by intersection of IRGs and HRGs lists (Fig. 1A). Then, we performed function analysis on these 31 genes and found that they were enriched in response to hypoxia, leukocyte migration and regulation of angiogenesis (Fig. 1B). Meanwhile, we created a PPI network to better clarify the interaction of 31 genes at protein level (Fig. 1C).
Consensus cluster analysis
A total of 31 hub genes were incorporated into cluster analysis. The results indicated that CDF value growth was flat when k = 2 and Delta area increased insignificantly at k > 3 (Fig. 2A,B). The fractal matrix showed the favorable intergroup difference and intragroup association, suggesting these pivot genes could categorize all CC samples into two subtypes (Cluster 1 (n=130) and Cluster 2 (n=174)). Therefore, the clustering stability was best for k = 2 (Fig. 2C). Survival analysis illustrated the significant difference in patient outcome between two clusters (Fig. 2D). PCA analysis uncovered the favorable distinction between the two clusters (Fig. 2E). Furthermore, 251 DEGs were collected from differential analysis between two clusters.
Development of a risk classifier
In the training set, we first determined 24 survival-associated indicators based on above 251 DEGs via univariate analysis (Fig. 3A). Then, the candidate genes were enrolled into LASSO regression to remove the over fitting genes (Fig. 3B,C). Finally, multivariate analysis was employed and seven hub genes were selected to develop an IHBRC (Table 1): risk score = (0.0131 × CCL20) + (0.0638 × CXCL2) + (0.2812 × ITGA5) + (0.0340 × PLOD2) + (0.0697 × PTGS2) + (0.0374 × TGFBI) + (0.1113 × VEGFA). In addition, Fig. 3D demonstrated the prognostic power of seven hub predictors.
As suggested by Fig. 4A, high-risk group presented a dismal prognosis benefit in the training set. The AUC values of 1-, 3-, and 5-year survival were 0.845, 0.699, and 0.654, respectively (Fig. 4B). We measured the survival outcome of patients in both groups and found that patients' outcomes were dismal as the risk score elevated (Fig. 4C). Meanwhile, we confirmed the performance of IHBRC in the validation and the entire cohorts using the same analysis described above and obtained the same results for the trend (Fig. 4D-I).
Table 1: Multivariate analysis of the seven model genes in cervical cancer.
Gene
|
Coefficient
|
P-value
|
CCL20
|
0.0131
|
0.007
|
CXCL2
|
0.0638
|
0.001
|
ITGA5
|
0.2812
|
0.001
|
PLOD2
|
0.0340
|
0.001
|
PTGS2
|
0.0697
|
0.008
|
TGFBI
|
0.0374
|
0.001
|
VEGFA
|
0.1113
|
0.001
|
Independent prognostic analysis
To examine the independent value of IHBRC in terms of survival of CC cases, univariate and multivariate analyses were employed. In the training set, univariate analysis demonstrated that low risk score was remarkably correlated with favorable prognosis (Fig. 5A). Furthermore, multivariate analysis still revealed that low risk score was independently associated with favorable outcome of CC patients (Fig. 5B), which could serve as an independent prognostic factor for glioma. These were confirmed by the test and the entire sets (Fig. 5C-F).
GSEA enrichment analysis
To explore the distinction in molecular pathways between the two groups, we applied GSEA based on hallmarks gene sets. The results disclosed that hallmarks including angiogenesis, hypoxia, IL6/JAK/STAT3 signaling, MTORC1 signaling, P53 signaling, TGF β signaling were markedly enriched in high-risk group (Fig. 6).
Immune infiltration analysis
In order to mirror the immune status of two groups, we estimated enrichment value of different immunocytes. Fig. 7A illustrated the relationship between seven model biomarkers and immunocytes. As shown in Fig. 7B, risk score was negatively correlated with the infiltration level of memory B cells, naïve B cells, resting dendritic cells, macrophages M1 and CD8 T cells, while neutrophils were activated in IHBRC-high cohort.
Analysis of immunotherapy and chemotherapy response
Waterfall diagrams indicated the mutational differences in the 20 genes between the two groups. We observed that the IHBRC-high cohort had a higher PIK3CA mutation rate than the IHBRC-low group (31 vs. 20%, Fig. 8A,B). Given the importance of TMB in evaluating immunotherapy response for patients with CC, we observed IHBRC-high group had lower TMB value (Fig. 8C). In addition, high risk score was correlated with a lower IC50 of docetaxel, doxorubicin and gemcitabine (p < 0.05), suggesting that the IHBRC served as a favorable indicator for chemosensitivity (Fig. 8D-F).
Construction of IHBRC-related regulatory network
The reciprocal regulation of mRNA and miRNA is closely bound up with tumor development. Based on the starbase online tool, we identified the target miRNAs of seven model genes with high relevance scores (Fig. 9). Moreover, miRNA set enrichment analysis was performed to explore the function of the target miRNAs by TAM 2.0 tool. The results showed these miRNAs were mainly involved in cell aging, apoptosis, immune response, inflammation and regulation of Stem Cell (Supplementary Table 1).