Construction of the Prognostic IRGPs model. TCGA's transcriptome data was used as a training cohort. From the list of immune-related genes(IRGs) obtained by IMMPORT, the genes in the transcriptome data were searched in turn, and a total of 1647 IRGs were found. These genes were also identified in the GSE39582. To ensure relatively high variation in the genes of the two platforms, we retained 325 IRGs with a median absolute deviation(MAD)>0.5. Then, 12,275 immune-related gene pairs (IRGPs) was built based on these 325 IRGs. After univariate Cox regression analysis of these IRGPs in the training group, there are still 28 potential prognostic IRGPs. By using Lasso Cox proportional hazard regression to define the model on the training set, 17 IRGPs were retained to form the final prognostic risk model. These IRGPs composed of 26 unique IRGs, most of them are antibiotics, cytokine receptors and cytokine related molecules(Table 1). Then, the risk score for each patient in the TCGA dataset was calculated based on the model. Finally, we used a time-dependent ROC curve analysis to classify patients into high- or low-immune risk groups. The optimal cut-off value for the risk score was set to -0.576(Figure 1). This value succeeded to stratify the patients in the training cohort into high- and low-risk groups. In other words, the overall survival(OS) of the low-risk group was significantly higher than that of the high-risk group(Figure 2A). We further performed univariate and multivariate Cox proportional hazards analyses to test whether the IRGPs model predicted survival independently of other prognostic factors in the TCGA cohort. Among these analyses, the risk score of the model can be used as an independent prognostic factor(Figure 3A, Figure 3B).
Table 1. List of immune-related genes for constructing prognostic models.
IRG1
|
Category
|
IRG2
|
Category
|
Coefficient
|
CXCL14
|
Cytokines
|
BST2
|
Antimicrobials
|
-0.32
|
RBP7
|
Antimicrobials
|
PTGS2
|
Antimicrobials
|
0.26
|
RBP7
|
Antimicrobials
|
ARG2
|
Antimicrobials
|
0.23
|
APOD
|
Antimicrobials
|
IL17RB
|
Cytokine_Receptors
|
0.05
|
C5AR1
|
Chemokine_Receptors
|
NR3C2
|
Cytokine_Receptors
|
0.18
|
IL10RA
|
Cytokine_Receptors
|
TNFRSF11A
|
Cytokine_Receptors
|
0.12
|
STC2
|
Cytokines
|
HNF4G
|
Cytokine_Receptors
|
0.29
|
RBP1
|
Antimicrobials
|
STC2
|
Cytokines
|
-0.63
|
GNAI1
|
Antimicrobials
|
GRP
|
Cytokines
|
-0.24
|
CCL4
|
Antimicrobials
|
INHBB
|
Cytokines
|
-0.19
|
ABCC4
|
Antimicrobials
|
GRP
|
Cytokines
|
-0.28
|
ARG2
|
Antimicrobials
|
GRP
|
Cytokines
|
-0.37
|
CCR7
|
Antimicrobials
|
INHBB
|
Cytokines
|
-0.31
|
CD86
|
Antimicrobials
|
IL7
|
Cytokines
|
0.28
|
INHBB
|
Cytokines
|
PDGFC
|
Cytokines
|
0.50
|
TNFRSF11A
|
Cytokine_Receptors
|
LCK
|
NaturalKiller_Cell_Cytotoxicity
|
-0.34
|
RORC
|
Cytokine_Receptors
|
PRKCQ
|
TCRsignalingPathway
|
-0.36
|
[i]Abbreviation: IRG, immune-related gene.
Verify the feasibility of IRGPs model for predicting survival. To determine whether the model had consistent prognostic value in different risk groups, we applied the model to GSE39582 as external validation. The patients in the verification cohort are divided into two parts according to the risk score. The OS of subgroups with low risk group increased significantly(Figure 2B). After performing univariate and multivariate Cox proportional hazards analyses in the validation group, we found that the results were similar to the training group, the high risk score of this model suggests a poor prognostic factor(Figure 3C, Figure 3D).
Immune cells infiltration in different risk groups. Previous studies have revealed that tumor infiltrating immune cells is related to prognosis[32]. To determine the infiltration of specific tumor immune cell subsets, we used CIBERSORT to estimate the relative proportion of 22 different immune cells per patient in different risk groups. The radar chart depicts a comparative summary of various immune cells in these two risk groups(Figure 4). In this result, we can find that Dendritic cells activated, Dendritic cells resting, Eosinophils, Macrophage M0, Monocytes, T cells CD4 memory resting and T cells regulatory(Tregs) were enriched in different risk groups. Among them, T cells regulatory (Tregs) and Macrophage M0 were significantly and highly expressed in the high-risk group, and the rest were highly expressed in the low-risk group(Figure 5).
Functional evaluation of the IRGPs model. In order to investigate the expression signatures of genetic perturbations that were significantly altered by the IRGPs model. GSEA was performed between the high-risk and low-risk groups in the TCGA cohort. We use a bubble chart (Figure 6) to show the overall results. In Figure 6, we found that genes in high-risk populations are enriched in stem cells and various advanced tumors. The top five genetic perturbations in the high-risk group were: stem cell up, breast cancer ductal invasive up, multicancer invasiveness signature, gastric cancer advanced vs early up and mammary stem cell up(Figure 7).