Establishment and validation of a novel immune-related prognostic signature in malignant pleural mesothelioma.

Background: Immune-related genes (IRGs) play an important role in the tumor immune microenvironment and affect tumor prognosis. This study aimed to establish a prognostic signature for malignant pleural mesothelioma (MPM) patients.
Methods: We obtained the relevant data of MPM patients in The Cancer Genome Atlas (TCGA), and univariate and multivariate Cox regression were used to construct the prediction signature and verify it with the external validation dataset GSE2549. A nomogram was then constructed, and its predictive ability was evaluated and analyzed the level of immune cell infiltration in different groups in the signature.
Results: An IRG-related prognostic signature composed of INHBA, CAT, SORT1, TNFSF13B, and BIRC5 was constructed, with patients divided into high-risk and low-risk groups according to the risk score. The survival time of overall survival (OS), progression-free survival (PFS), disease-free interval (DFI), and relapse-free survival (RFS) in low-risk groups was longer than in high-risk groups. Furthermore, the signature had high predictive performance, and the receiver operating characteristic (ROC) of 1, 2, and 3 years could reach 0.853, 0.881, and 0.914, respectively. The predictive accuracy of the signature was verified by using the independent GSE2549 dataset. The levels of activated CD4 T cells, immature dendritic cells, and type 2 T helper cells were higher in high-risk patients. The gene set enrichment analysis (GSEA) analysis showed that a high concentration and P53 signal pathways were found in high-risk groups. Conclusions: This research developed and verified a new type of immune prognostic signature based on five IRGs, which can predict the prognosis of tumor patients and provide new ideas for individualized treatment. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Malignant pleural mesothelioma (MPM) is a rare malignant thoracic cavity cancer closely related to asbestos exposure and has an exceedingly poor prognosis with a median survival of 9.5 months (1). Cisplatin combined with pemetrexed systemic chemotherapy is used as the first-line treatment for the treatment of advanced MPM, and although this extends the median survival of patients who cannot undergo surgical resection (2,3), the overall 5-year survival rate is lower than 5% (4). While the gold standard for diagnosis of MPM is thoracoscopic biopsy, it is not suitable for some patients due to weakness, comorbidities, or personal choice, and there are currently few methods to evaluate and predict the prognosis of the disease. Therefore, there is considerable value in establishing novel and predictive signatures that can reliably estimate clinical outcomes in monitoring personalized prognosis and guiding clinical management of MPM patients.

As a promising antitumor strategy, immunotherapy has provided new methods for treating solid tumors in the past few years and has obtained long-lasting responses and long-term survival benefits for countless cancer patients (5). Immune checkpoints have shown promise as therapeutic targets for various types of tumors, and related studies have shown that genetic changes in MPM, including point mutations, minute deletions, and copy number changes, may lead to new antigens. In view of the crosstalk between tumor cells and the immune microenvironment is indeed closely related to pathogenesis of MPM (6,7), a better understanding of this immune microenvironment with targeted molecular vulnerabilities is particularly associated with the improvement of treatment approaches and patient outcomes. Due to the related factors and the clonal expansion of tumor-infiltrating T lymphocytes support the concept of MPM as an immunogenic disease (8-10), and related studies highlighted the prognostic role of lymphocytes and the occurrence of immunosuppression in the disease (11-13). Therefore, an immune-related prognostic model for MPM is worth exploring. Although there is a similar study has reported a signature for predicting overall survival (OS) of MPM (14), are there other mRNAs that affect the prognosis of patients? Compared with traditional features (such as stage, T-, N-, M-stage), how accurate is the prognosis prediction signature based on these mRNAs? These are the questions that our research needs to explore.

In our study, based on the immune-related genes (IRGs) of MPM in The Cancer Genome Atlas (TCGA) database, we identified and validated an immune-related 5-gene signature and verified the stability of the model through the external validator GSE2549. In addition, we analyzed the correlation between genes in the model and some clinical features and immune cell infiltration. We used comprehensive analysis to develop a new marker based on cancer immune genomics characteristics to predict MPM prognosis and provide targets for its diagnosis and prognosis. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-527/rc).

Methods

Data collection

We obtained the sequencing data and corresponding clinical information of MPM in TCGA database (https://cancergenome.nih.gov/) and acquired the microarray dataset GSE2549 (15), GSE51024 (16) from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). As all data are publicly available, ethics committee approval was not required. After converting the ID in TCGA into gene symbols according to the GENCODE file, a list of IRGs was obtained from the Immunology Database and Analysis Portal (ImmPort) (17). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data analyzing and identification of IRGs

Using the “limma” (18) package in R to obtain differential expression genes (DEGs) from GSE51024, the cut-off criteria for screening DEGs were |log2(fold change)| ≥1 and P<0.05. Differential expression IRGs (DEIRGs) were obtained through the intersection of IRGs and DEGs from GSE51204.

Establishment of DEIRGs prognostic signature

We established the DEIRGs signature for predicting the prognosis of MPM patients. Firstly, the “survival” R package was used to perform univariate Cox analysis to screen out genes that were significantly relevant to the OS (P<0.05), and multivariate Cox regression analysis was then performed to obtain the coefficient. According to the multivariate Cox regression analysis results, a prognostic signature was established, and the risk score calculated based as follows = regression coefficients1 × expression level of mRNA1 + regression coefficients2 × expression level of mRNA2 +… + regression coefficientsn × expression level of mRNAn.

Validation of the prognostic signature

To verify the prediction accuracy of the signature, the risk score of MPM patients was obtained from the prognostic signature, and patients were divided into high-risk groups and low-risk groups according to the median. The receiver operating characteristic (ROC) curve and area under the ROC curve (AUC value) were performed to evaluate the predictive value of the prognostic signature, and the external validation dataset GSE2549 was used to validate it.

Validation of the independence of the prognostic signature and nomogram analysis

To detect whether the risk score of patients evaluated by the prognostic signature of IRGs was independent of the clinical characteristics as a prognostic factor of MPM, the risk score and clinical characteristics (age, gender, T stage, N stage, and clinical stage) were combined. The “pheatmap” R package was used to analyze the relationship between the molecular characteristics of MPM patients and other clinical variables in the high-risk and low-risk groups, and the “rms” (19) R package based on the multivariate Cox regression analysis was used to construct the nomogram based on independent prognostic factors, including risk score and M stage. The performance of the nomogram was evaluated by using the ROC curve, and the predictive accuracy of the signature we constructed was evaluated by comparing it with the signature built by previous research (20,21).

Analysis of signature and clinicopathological characteristics

To further explore the correlation between the different groups and clinicopathological characteristics, patients were grouped according to their clinical information and the survival status of the high-risk and low-risk patients was compared.

Immune cells infiltration analysis

To explore the relationship between genes and immune cell infiltration levels in the constructed signature, we used the single-sample gene set enrichment analysis (ssGSEA) algorithm in the “GSVA” package (22) to quantify the relative abundance of 10 immune cell types (23) infiltration in MPM. Pearson correlation analysis was used to detect the correlation between the gene expression level and the infiltration of each immune cell.

GSEA enrichment analysis

Pathway differences between the high-risk and low-risk groups were investigated by GSEA enrichment analysis (24,25). mRNA expression microarray in the high- and low-risk groups in MPM patients was taken as the expression dataset, and c2.cp.kegg.v7.4.symbols.gmt was selected as the enrichment analysis gene set to run GSEA software. Displacement test times were set as 1,000, and a gene set with false discovery rate (FDR) <0.25 was considered significant enrichment.

Statistical analyses

All statistical analyses and related visualization were conducted to determine independent prognostic factors using the R package (R software version 3.6.3) and Perl (5.30.1) software.

Results

Identification of IRGs

MPM data in TCGA has only tumor patients and lacks normal tissues for comparison. Therefore, GSE51024 was analyzed by using the “limma” package in R, and by setting the criterion of |log2(fold change)| >1 and P<0.05 and comparing the tumor tissue with normal tissue, a total of 1,212 DEGs were obtained. We then intersected the 1,811 IRGs downloaded from the Immport database with the 1,212 DEGs to obtain 101 differentially expressed IRGs (DEIRGs) (Figure S1, ). Relevant information of patients is shown in .

Table 1

Interaction of DEGs from GSE51024 and immune-related genes

Differential expression immune-related genes

EDN1, IGF1, CXCL3, FGR, CXCL13, A2M, TEK, BIRC5, C3, SORT1, KL, THRB, EDNRA, PPBP, OXTR, FGFR3, IL6R, ADRB1, CTSB, SEMA3G, PDGFD, SPP1, SFTPD, PTGS2, AGER, IL20RA, IL18, BMP2, TNFRSF19, FGFR4, AREG, CXCL17, NR3C2, CRABP2, IL17RD, SCG2, CXCL2, WFDC2, FABP6, TGFBR3, CAT, LYZ, MARCO, CLEC11A, S100A12, BMP5, AGTR2, PROCR, PTGER3, CXCL14, PLAU, S100P, CLDN4, PLA2G2A, PLTP, HBEGF, TNFSF10, SLC40A1, IL1RL1, SDC4, S100A14, CSF3, EDNRB, CCL20, ARRB1, SCGB3A1, ADRB2, MUC4, S100A8, FGF14, ANGPT1, DUOX1, PDGFRL, RORA, TNFSF4, INHBA, LIFR, DMBT1, S1PR1, FGF9, SDC1, CDH1, TNFSF13B, AZGP1, CX3CL1, ICAM1, CALCRL, FGF18, LTF, CXCL1, DKK1, IL7R, VCAM1, ROBO2, IL33, PGC, CXCL5, FABP4, KITLG, VIPR1, SEMA3C

DEGs, differential expression genes.

Table 2

Clinical features of the MPM patients involved in this study

Variables	TCGA	GSE2549
Number of patients	84	39
Age (years), median [range]	64 [28–82]	NA
Gender
Female	15	NA
Male	69	NA
Status
Alive	12	8
Dead	72	31
Clinical stage
1	10	NA
2	15	NA
3	43	NA
4	16	NA
T stage
T1	14	NA
T2	24	NA
T3	31	NA
T4	13	NA
Unknow	2	NA
N stage
N0	43	NA
N1	10	NA
N2	25	NA
N3	3	NA
Unknown	3	NA
M stage
M0	77	NA
M1	7	NA

MPM, malignant pleural mesothelioma; TCGA, The Cancer Genome Atlas; NA, not available.

Development and validation of DEIRGs signature

The prognostic value of DEIRGs in MPM was discussed according to their prominent role in the tumor microenvironment. We performed univariable Cox regression to analyze the 101 DEIRGs and survival data from 82 patients, and based on the univariable Cox regression, obtained 21 prognosis-related DEIRGs. Subsequently, based on multivariate Cox regression analysis, we finally obtained five DEIRGs; BIRC5, SORT1, CAT, INHBA, and TNFSF13B, to construct a prognostic signature. Based on the regression coefficients, risk score = (0.0889 × expression of BIRC5) + (0.0523 × expression of SORT1) + (−0.0358 × expression of CAT) + (0.0446 × expression of INHBA) + (−0.0222 × expression of TNFSF13B). According to the median value of the risk score, all patients were grouped into the high-risk and low-risk groups. By drawing the survival curve, the OS of the high-risk group was significantly lower than the low-risk group (). According to Liu’s research (26), we found that in multiple clinical evaluation indicators such as relapse-free survival (RFS), progression-free survival (PFS), and disease-free interval (DFI), patients with high-risk groups had shorter survival time compared to patients with low-risk groups (). The ROC curve used to evaluate the predictive value of the signature showed good accuracy with 1 year (AUC =0.853), 2 years (AUC =0.881), and 3 years (AUC =0.914) (). The heatmap shows and compares the expression characteristics of the five genes signature between high-risk groups and low-risk groups (), and the death sample survival time decreased with an increase in the risk score are illustrated in . The results show the constructed signature in the independent validation set GSE2549 and demonstrate patients in the high-risk group had shorter OS compared with low-risk groups (). The GSE2549 dataset was also used to evaluate the performance of the DEIRGs signature (), and each patient’s risk score was calculated with a signature.

Figure 1

Construction of DEIRGs signature and verification. (A-D) Kaplan-Meier curves for patients in the high-risk group and low-risk group in the TCGA cohort: (A) Kaplan-Meier curves for the OS of patients; (B) Kaplan-Meier curves for the RFS of patients; (C) Kaplan-Meier curves for the PFS of patients; (D) Kaplan-Meier curves for the DFI of patients. (E) ROC curve for the prediction of MPM survival for 1, 2, and 3 years. (F) Heatmap showing the prognostic DEIRGs signature in the TCGA cohort. (G) Overall survival time of MPM patients in TCGA cohort. (H) The risk score distribution in TCGA cohort. (I) Kaplan-Meier curves for the OS of patients in the high-risk and low-risk groups in GSE51024. (J) Heatmap showing the prognostic DEIRGs signature in the GSE51024. (K) Overall survival time of MPM patients in GSE51024. (L) The risk score distribution in GSE51024. DEIRGs, differential expression immune-related genes; TCGA, The Cancer Genome Atlas; OS, overall survival; RFS, relapse-free survival; PFS, progression-free survival; DFI, disease-free interval; ROC, receiver operating characteristic; MPM, malignant pleural mesothelioma.

Independent prognostic factor analysis and construction of nomogram

To evaluate the role of the DEIRGs signature in some common clinicopathological features such as age, sex, T stage, N stage, M stage, and tumor stage, univariate and multivariate Cox regression analyses were performed, and the results confirm that the risk score was an independent prognostic indicator of the TCGA cohort (). Taking into account the significance of the clinical utility we constructed, we drew the decision curve analysis (DCA) curve, and as shown in , the signature combined with the risk was more beneficial than other factors (age, gender, stage). Further, based on multivariate Cox regression analysis results, we established a novel nomogram incorporating two independent factors to predict the 1-, 2-, and 3-year survival rates of MPM patients (), and the AUC values for OS were 0.852, 0.875, and 0.878, respectively (). In comparing the prediction performance of multiple signatures, these results show the prediction performance of our signature is superior to those previously published ().

Figure 2

Evaluation of the signature. (A,B) Univariate and multivariate Cox regression analyses regarding OS in MPM patients. (C) DCA for assessing the clinical utility of the nomogram. (D) Nomogram for the prediction of OS at 1, 2, and 3 years. (E-G) Time-dependent ROCs for 1-, 2-, and 3-year OS of the nomogram. OS, overall survival; MPM, malignant pleural mesothelioma; DCA, Decision Curve Analysis; ROC, receiver operating characteristic; AUC, area under the ROC curve.

Figure 3

Comparing the prediction performance of multiple signatures. (A-C) Time-dependent ROC analysis of the signature for predicting the OS of patients. (D-F) Clinical-related ROC analysis of the signature. (A,D) The signature in this study. (B,E) The signature of Bai’s research (20). (C,F) The signature of Zhou’s research (21). ROC, receiver operating characteristic; OS, overall survival; AUC, area under the ROC curve.

Subgroup analysis of DEIRGs signature for survival prediction

Subgroup analysis was performed to assess the prognostic value of the DEIRGs signature within the same clinicopathological risk factors by assigning patients into different subgroups, including a younger group (age ≤65 years), an elder group (age >65 years), a male group, an advanced-stage group (stage III & IV), a group with patients of T1 & T2, a group with patients of T3 & T4, and a group with patients of M0. Subgroup analysis demonstrated the DEIRGs signature could still distinguish patients into different survival groups within all subgroups with statistically significant prognostic ability ().

Figure 4

Subgroup analysis for the DEIRGs signature. Patients were stratified into seven subgroups for survival analysis based on age ≤65 years (A), age >65 years (B), patients with M0 (C), gender (male) (D), stage III and IV (E), T1 and T2 (F), T3 and T4 (G). Group 1: low-risk group; Group 2: high-risk group. DEIRGs, differential expression immune-related genes.

Signature associated with clinical characteristics

The heatmap revealed the levels of genes in the high-risk group and low-risk groups in the TCGA MPM cohort were different. However, the association between the high-risk and low-risk groups and T, N, M, stage, and other clinicopathological characteristics were not statistically significant (). Through survival analysis of the genes in the signature, we found that the high expression of BIRC5, SORT1, and INHBA was related to poor prognosis in MPM, while the high expression of CAT and TNFSF13B was closely related to a better OS (). In addition, compared with patients with lymph node metastasis (Nx), the expression level of CAT in patients without lymph node metastasis (N0) was lower, and the expression level of SORT1 in patients without distant metastasis was higher than in patients with distant metastasis. BIRC5, INHBA, and risk score had higher expression levels in patients who died than in those who survived ().

Figure 5

Clinical evaluation by the risk assessment signature. (A) Heatmaps of various clinicopathological characteristics in high- and low-risk groups. (B-F) Kaplan-Meier survival curves of five genes in the signature: (B) BIRC5; (C) SORT1; (D) CAT; (E) INHBA; (F) TNFSF13B. (G-K) Gene expression in different groups: (G) CAT with N; (H) SORT1 with M; (I) BIRC5 with survival state; (J) INHBA with survival state; (K) risk score with survival state.

Immune infiltration analysis

By analyzing the data of immune infiltration results downloaded from TIMER, we found that, according to the difference of high- and low-risk, the different types of immune cells calculated by algorithms in different databases (TIMER, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC) were also different (). To further explore the correlation between the DEIRG signature and the level of immune infiltration, ssGSEA was used to calculate ten immune cell infiltration scores. The high-risk group was positively correlated with the scores of three immune cells, including activated CD4 T cells, immature dendritic cells, and type 2 T helper cells, and negatively associated with natural killer cell (P<0.05, ). We further explored the correlation between the expression level of genes in the signature and the level of immune cell infiltration, and the results showed that the gene in the signature had a significant correlation with natural killer cells and immature dendritic cells except for CAT. In addition, the correlation between BIRC5 and type 2 T helper cells was as high as 0.93, and there was a negative correlation between CAT and activated CD4 T cell and type 2 T helper cells ().

Figure 6

Estimation of tumor-infiltrating cells by signature. (A) Heatmap for immune responses based on CIBERSORT, ESTIMATE, MCPcounter, ssGSEA, and TIMER algorithms among high- and low-risk groups. (B) ssGSEA scores in the high- and low-risk patients in the TCGA cohort. (C-G). Correlation between genes and immune infiltrating cells: (C) BIRC5; (D) INHBA; (E) CAT; (F) SORT1; (G) TNFSF13B. *, P<0.05; ***, P<0.001. ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; ns, no significance.

GSEA analysis between high-risk and low-risk groups

We also performed GSEA to explore the biological signal pathway between the high-risk and low-risk groups. The results showed that the significantly enriched pathways in the high-risk group involved cell cycle, DNA replication, extracellular matrix (ECM) receptor interaction, focal adhesion, homogeneous recombination, the P53 signaling pathway, regulation of actin cytoskeleton, and the transforming growth factor-β (TGF-β) signaling pathway (FDR <0.01) ().

Figure 7

GSEA analysis in the TCGA based on the high-risk and low-risk groups. GSEA, gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

Discussion

MPM is a very rare malignant tumor, for which current treatments of surgery, chemotherapy, and radiotherapy provide limited clinical efficacy and a poor prognosis. As exploring molecular biomarkers is of great significance to the early diagnosis, prognosis prediction, and treatment strategies of MPM, we constructed a prognostic signature in this study consisting of five immune genes (INHBA, CAT, SORT1, TNFSF13B, and BIRC5).

Based on the relevant data of MPM in TCGA, we used univariate and multivariate cox regression analysis to build a prognostic signature by five genes, including INHBA, CAT, SORT1, TNFSF13B, and BIRC5. In this signature, the ROC values of 1, 2, and 3 years were 0.853, 0.881, and 0.914, respectively, showing good accuracy. In comparing the prediction performance of multiple signatures, these results show the prediction performance of our signature is superior to those previously published. Furthermore, low-risk patients had better OS than those in the high-risk group, and compared with high-risk patients, low-risk patients also showed the advantage of longer survival in RFS, PFS, and DFI. We also used an independent GEO dataset to verify the predictability of the signature, and a nomogram was drawn to provide a personalized estimate of the potential survival rate and contribute to the personalized management of MPM. Tumor-infiltrating cells have been proven to have vital functions in cancer, and we analyzed the immune infiltration levels of patients in different risk groups. This showed they were in activated CD4 T cells, immature dendritic cells, and type 2 T helper cells, and the fraction of the high-risk group was higher than that of the low-risk group. In natural killer T cells, the fraction of the low-risk group was higher than that of the high-risk group, while GSEA revealed that patients in high-risk and low-risk groups had different levels of enrichment pathways. Compared with the prognostic characteristics of MPM previously constructed, the signature established in this study had higher predictive performance.

Furthermore, in this prognostic signature, five genes (INHBA, CAT, SORT1, TNFSF13B, and BIRC5) were selected as crucial genes. SORT1 is a protein of the sortilin family, and according to existing reports, has different roles in different types of cancer. Targeted SORT1 has been reported to inhibit tumor metastasis and promotes tumor apoptosis in breast (27) and ovarian (28) cancer, while high expression of SORT1 may promote tumor progression in gastric cancer (29) and lung cancer (30). TNFSF13B is a member of the TNF ligands family and has been shown to play an important role in the proliferation and differentiation of B cells (31). While previous research on the TNFSF13B mainly focused on the immune system and hematological tumors, one study showed that as a B cell-activating cytokine, TNFSF13B could regulate T cell function by increasing T cell activation and TH1 polarization, enhancing the expression of pro-inflammatory leukocyte trafficking chemokine CCR6 and promoting memory phenotype, thereby enhancing antitumor immunity. Furthermore, TNFSF13B also has an apparent immune regulation function, promoting CD4+Foxp3+ spleen and Treg in the tumor microenvironment (32). In this study, TNFSF13B was used as a protective factor in MPM, and in the TISIDB database (33), we found that belimumab is an intravenous immunosuppressant for adjuvant treatment of systemic lupus erythematosus by acting on the target of TNFSF13B, which also provides new ideas for treating the disease. BIRC5, also known as survivin, is a negative regulatory protein that prevents apoptosis (34) and can regulate cancer development by inhibiting cell apoptosis and inducing cell proliferation (35). The abnormal amplification of BIRC5 has been observed in a variety of malignant tumors, and its overexpression is closely related to tumor occurrence and development, radiotherapy and chemotherapy resistance, and the poor prognosis of cancer patients (36-38). Recently, BIRC5 has attracted much attention as a molecular target for cancer therapy. CAT (catalase) is a key enzyme for H2O2 metabolism, and its expression and localization changes significantly in tumors (39). It can also ameliorate this hypoxic environment and reverse the malignant phenotype of tumor cells (40). INHBA is a member of the TGF-β superfamily of the gene encoding protein, and increasing reports have demonstrated it is a crucial factor in some solid malignant tumors, including prostate cancer (41), esophagus squamous cell carcinoma (42), and bladder cancer (43). Chen’s study found that knocking down the expression of INHBA could inhibit gastric cancer progression through the TGF-β signaling pathway (44), and another study found that interfering with the expression of INHBA could significantly reduce the degree of cell proliferation, invasion, and migration of colon cancer (45). In this study, we found a high expression of INHBA was related to the poor prognosis of MPM patients, which raises the question of whether reducing its expression could prolong their survival.

The GSEA showed patients in the high-risk group concentrated on the P53 signaling pathway, which is probably the most commonly inactivated protein in cancer (46,47), and TP53 mutations are associated with poor OS (48). Related study has shown that miR-320a is significantly inhibited in MPM samples, P53 can target PD-L1 by mediating miR-320a (49), and PD-L1 is highly expressed in mesothelioma and tumor stroma. P53 is associated with poor survival (50-52), providing a new direction for treating MPM. In addition, microRNA-215-5p can be used to suppress the progression of mesothelioma through the MDM2-p53-signal axis (53).

Inevitably, this study has some limitations. Firstly, the prognostic model proposed was established and verified using data from public databases and requires verification in clinical trials. Secondly, there remains a need for in vivo and in vitro experiments to further understand the biological function of the role of the DEIRG signature in MPM.

Conclusions

We constructed a prognostic signature and nomogram of MPM patients which can predict the prognosis of individual cases. Moreover, exploring the relationship between high and low-risk group immune cell infiltration and clinicopathologic features provides a new direction for guiding individualized and accurate treatment.

The article’s supplementary files as