Integrated analysis reveals the potential of cluster of differentiation 86 as a key biomarker in high-grade glioma

Data source of HGG patients and expression analysis on CD86

This study first performed the pan-cancer analysis to explore the expression level of CD86 in tumor and normal tissues. Then the expression of CD86 in patients with HGG was explored. The HGG patients were selected from the GBM and LGG datasets in the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A total of 1114 samples of TCGA_GBM and TCGA_LGG were obtained from the TCGA database, and the inclusion criteria of HGG patients were as follows: 1) had the data of tumor grade; 2) WHO grade III and grade IV; 3) had the data of survival time and survival status; 4) survival time≥1 month; 5) had the data of age, gender, IDH mutation status, Chr_1p_19q codeletion, MGMT methylation status, radiotherapy history and chemotherapy history; 6) primary tumor and had gene expression RNA-Seq data. Finally, 298 patients who had the above complete information were enrolled in our further analysis. The baseline data of 298 HGG patients was presented in Table 1.

Then the expression difference of CD86 in normal and HGG groups was analyzed. In addition, the expression pattern of CD86 in subgroup populations of patients with HGG was determined after stratifying patients by clinical characteristics, including different age (<60, ≥60), gender (male, female), IDH mutation status (wildtype, mutant), Chr_1p_19q codeletion (none-codel, codel), MGMT methylation status (unmethylated, methylated), radiotherapy history (yes, no) and chemotherapy history (yes, no). Further, all the HGG patients were divided into 2 groups according to CD86 expression level (high, low), and the association between CD86 expression and the above clinical characteristics was evaluated by Logistic regression analysis.

Prognostic value of CD86

The prognostic value of CD86 in HGG was then explored. Firstly, all HGG patients were divided into CD86 high and low expression groups according to the best cut-off value of CD86 expression level by maxstat R package. Then the Kaplan-Meier analysis by R survival package and log-rank test were performed to explore the association of CD86 expression with the overall survival (OS) time of patients. The receiver operating characteristic curve (ROC) analysis by R package pROC was used to assess the prediction performance of CD86 on OS. In addition, we also verified the prognostic value of CD86 in HGG using the data from the Chinese Glioma Genome Atlas (CGGC) database (http://www.cgga.org.cn/).

To further disclose the prognostic value, we performed a subgroup survival analysis on CD86 regarding different clinical characteristics based on TCGA data. Cox regression analyses by survival R package were then used to explore the association of variables with OS and the independent prognostic value of CD86. We also transformed the CD86 expression value into the ordinal variable and explored the association between CD86 expression groups and OS with or without adjusting confounders by Logistic regression analysis.

Comprehensive nomogram model establishment

Based on the results of multivariable Cox regression analysis, we also established a comprehensive nomogram model containing independent factors (including CD86) and obtained the risk score of each patient by the rms R package. The Kaplan-Meier analysis and log-rank test were used to explore the survival difference between high and low risk score groups. ROC and Decision Curve Analysis (DCA) analyses were performed to evaluate the performance of the nomogram model in HGG.

Pathway enrichment analysis on CD86

Due to the importance of CD86 in HGG, this study subsequently explored the potential regulatory mechanism associated with CD86. The c2.cp.kegg.v7.4.symbols.gmt was obtained from the Molecular Signature Database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). All patients were divided into high and low expression groups according to the best cut-off value of CD86 expression value, and then Gene set enrichment analysis (GSEA) was used to evaluate the significant pathway enriched in 2 groups. Further, a single sample Gene Set Enrichment Analysis (ssGSEA) algorithm was performed to quantify the significant pathway score of each sample in HGG. The correlation between CD86 and pathway score was explored.

Immune infiltration analysis on CD86

This study also explored the correlation of CD86 expression with immune infiltrates in HGG. Through the gene expression profile, we obtained the immune score of 22 immune cells by the CIBERSORT algorithm. The Pearson method was used to assess the correlation of CD86 expression level with the abundance of immune cells. Further, the immune cells with significant correlation with CD86 were then selected for comparing the difference in infiltration level between CD86 high and low expression groups.

The source of radiomics data

Due to the potential of CD86 as an important biomarker in HGG, we innovatively explored the possibility of radiomics characteristics of MRI images for predicting CD86 expression without invasiveness. A total of 458 patients with MRI images in TCGA_GBM and TCGA_LGG datasets were obtained from the Cancer Imaging Archive (TCIA) (https://public.cancerimagingarchive.net/nbia-search/) database. The inclusion criteria of images were as follows: 1) enhanced MRI CE-T1WI images; 2) preoperative images; 3) had sufficient quality; 4) had corresponding TCGA data. Finally, 89 HGG patients with complete TCGA data and MRI images were selected for further radiomics analysis.

Region of interest segmentation and features extraction of MRI images

Two experienced radiologists, blinded to the clinical data of the patients, performed region of interest (ROI) segmentation using 3D Slicer software. Then the open-source Python package PyRadiomics was used to extract the radiomic features. The extracted features included (1) shape-based features, such as three-dimensional size and shape of the ROI; (2) first-order features, describing the distribution of the voxel intensities within the ROI using commonly used and basic metrics; and (3) texture features, including the gray-level co-occurrence matrix (GLCM), gray level dependence matrix (GLDM), gray-level run-length matrix (GLRLM), gray-level size-zone matrix (GLSZM), neighboring gray-tone difference matrix (NGTDM). In total, 107 radiomics features were extracted from each HGG image. The intraclass correlation efficient (ICC) was performed to evaluate the consistency of the (ROI) segmentation by two physicians. Finally, 84 radiomics features with ICC≥0.8 were selected for further analysis.

Radiomics model establishment and performance validation

Before model establishment, 89 HGG patients with complete TCGA data and MRI images were assigned to the training set and validation set according to the 6:4 ratio, and their baseline distribution between the training set and validation set had no difference. Subsequently, a radiomics model was constructed in the training set. To determine the more important features among 84 radiomics features, we used the mRMR (Maximum relevance, minimum redundancy) algorithm by mRMRe R package to filter radiomics features among them. After mRMR analysis, 20 radiomics features were identified. Subsequently, we performed RFE (Recursive feature elimination) algorithm on 20 features by caret package to obtain the optimal feature sets. RFE analysis can rank all the radiomics features and gradually remove the less important factors. After RFE analysis, the optimal features were determined and used for the further prediction model establishment.

Further, we used the optimal features to construct the radiomics model for predicting CD86 expression in the training set. Support Vector Machine (SVM) and Logistic regression (LR) were the common algorithms to establish the machine learning model in recent studies [14–16]. Therefore, this study only used these 2 methods to develop the prediction model based on the radiomics features. SVM and LR algorithms were performed with caret and glm R packages, respectively, and the variables’ importance in the two models was assessed. Further, the radiomics score (RS) integrating the optimal features in 2 models was calculated. Calibration curve, ROC, precision-recall (PR) curve, and DCA were performed to evaluate the performance of 2 radiomics models in the training set. Further, ROC, PR, and DCA were also performed in the validation set to assess the stability of our established models. Further, the RS between CD86 high and low expression groups was also compared based on the whole population.

Comprehensive model establishment containing radiomics score and clinical characteristics

This study further explored the clinical value of the radiomic model in HGG. Based on the RS from the SVM method, the Kaplan-Meier method and log-rank test were performed to explore the survival difference between RS high and low groups. The multivariable Cox regression analysis was conducted to identify the independent factors associated with the OS among RS and clinical characteristics. Then the significant independent factors were enrolled into the establishment of a comprehensive nomogram model. ROC and Calibration curves were performed to evaluate the performance of the comprehensive model for predicting the prognosis of HGG patients.

Statistical analysis

The data used in this study were analyzed with SPSS software and R packages. The t-test and χ² test were used to compare the quantitative and qualitative data between the 2 groups, respectively. The association of CD86 with clinical characteristics and death risk was analyzed with Logistic regression. The Pearson method was used to explore the correlation between CD86 expression and the infiltration level of immune cells. The survival difference of the Kaplan-Meier curve was compared with the log-rank test. Univariable Cox regression analysis was used to explore the correlation of variables with OS time, and multivariable Cox regression analysis was used to determine the independent factors. P value<0.05 was considered as a statistically significant difference.

Availability of data and materials

The dataset used and/or analyzed during the current study is available from the corresponding author on reasonable request.

The expression of CD86

This study first explored the expression of CD86 in human cancers. The pan-cancer analysis showed that CD86 was aberrantly expressed in most cancers compared with that in normal tissues (Figure 1A). Regarding HGG, CD86 was significantly upregulated in tumor tissues (Figure 1B, P<0.001). We also compared the expression difference of CD86 among HGG patients based on different clinical characteristics, finding that the expression of CD86 showed no difference in different genders, with or without radiotherapy and chemotherapy groups (all P>0.05). Higher expression was significantly observed in patients with age≥60 years, non-codel of Chr_1p_19q, unmethylated MGMT promoter, wildtype IDH, and grade IV (all P<0.05). Logistic regression analysis further indicated the correlation of CD86 with these characteristics (Table 2).

The prognostic value of CD86

This study further assessed the potential clinical value of CD86 from the aspect of prognostic impact. The survival analysis showed that high expression of CD86 was significantly related to the poor OS time of HGG patients (Figure 2A, HR=2.96, P<0.001). The ROC analysis showed the prediction ability of CD86 on the 1, 3, and 5-year overall survival. We also used the data of CGGA to verify the clinical value of CD86 in HGG, and the results were consistent with that of TCGA (Figure 2B). These results confirmed the clinical value of CD86 in HGG.

We also assessed the prognostic value of CD86 on the overall survival in different subgroups based on TCGA data. Survival analysis further showed that CD86 high expression was significantly associated with poor prognosis in all subgroups (Figure 3, all HR>1, all P<0.05). These results indicated that CD86 can favorably distinguish the patient’s prognosis.

To further verify the prognostic value of CD86, we performed a Cox regression analysis. The univariate Cox regression analysis (Table 3) showed that 7 variables including CD86 (HR=1.604, P<0.001) were significantly related to the OS of HGG patients. We then enrolled the variables with P<0.05 into the multivariable Cox regression analysis, finding that CD86 was independently related to the prognosis of HGG patients (HR=1.213, P=0.035). We then translated the CD86 expression value into a qualitative variable and explored the correlation of the CD86 group with patients’ death risk, finding that the death risk of HGG patients in the crude model increased with the increase of CD86 expression (Table 4, P for trend<0.001). After adjusting the confounders of age and gender, the correlation was still observed (P for trend=0.002). Especially when the expression of CD86 was larger than 3.745, HGG patients had twice as death risk as these patients with CD86 expression less than 2.678 (OR=2.657, P=0.002).

Comprehensive nomogram model establishment containing CD86 and clinical characteristics

Based on the multivariable Cox regression analysis, we established a comprehensive nomogram model containing the independent factors associated with the OS of patients. The nomogram analysis showed that the comprehensive model had 0.1-0.95, 0.05-0.95, and 0.05-0.9 probabilities for predicting the 1, 3, and 5-year survival, respectively (Figure 4A). Overall C-index value of nomogram model was 0.819 (95%CI: 0.782-0.855, P<0.001). The risk score of the nomogram model was significantly related to the OS, and patients with higher risk scores had poorer prognoses (Figure 4B, HR=10.82, P<0.001). The ROC and DCA analyses showed that risk score had favorable prediction performance (Figure 4C, 4D).

Potential mechanism exploration of CD86

The above analyses have disclosed the potential of CD86 as an important biomarker in HGG, and we then predicted the possible regulatory mechanism associated with CD86. The GSEA enrichment analysis indicated that CD86 was significantly related to the pathway of natural killer (NK) cell-mediated cytotoxicity (Figure 5A, P<0.001), and this way was enriched in the high expression class. We then performed the single sample GSEA and obtained the score of the natural killer cell-mediated cytotoxicity pathway in each sample, finding that the pathway score was positively associated with CD86 expression in HGG patients (Figure 5B, P<0.05). The enrichment analysis showed that CD86 might be related to immunity.

To explore the indicative roles of CD86 on the tumor microenvironment (TME), CIBERSORT was adopted to evaluate the relative proportion of 22 kinds of immune cells. The analysis indicated that CD86 expression showed a negative correlation with the infiltration level of Plasma_cells, T_cells_CD8, T_cells_CD4_naïve, T_cells_follicular_helper, NK_cells_activated, and Dendritic_cells_resting, but positive correlation with Monocytes infiltration (Figure 6A, all P<0.05). Further comparison analysis (Figure 6B) found that the CD86 high expression group had higher infiltration of Monocytes (P<0.001), but lower infiltration of the antitumor immune cell of NK_cells_activated (P=0.02). It followed that CD86 might be significantly involved in the tumor immune microenvironment in HGG.

Radiomics features extraction and prediction model establishment

The above results have disclosed the potential of CD86 as an immune and prognostic biomarker in HGG. Detection and prediction of CD86 can help us to monitor disease progression and evaluate the immune response. However, the detection of CD86 was an invasive. Therefore, the innovative point of this study is the construction of a noninvasive radiomics model for predicting CD86 expression. We used 2 common modeling approaches to evaluate the application of the model with radiomics features.

A total of 89 HGG patients with MRI images and 107 radiomics features were enrolled in the analyses. The patients were first divided into a training set and a validation set according to the 6:4 ratio. Then we developed a model based on the radiomics features for predicting the CD86 expression in the training set. To filter effective features among all the radiomics features, the mRMR algorithm was initially performed and 20 radiomics features were identified. Following by RFE algorithm, 5 most important features were finally determined, including original ngtdm Contrast, original glcm Maximum Probability, original glszm GrayLevelNonUniformity, original gldm DependenceNonUniformity, original shape Maximum2DDiameterColumn (Figure 7A). Further, we established 2 radiomics models with 5 features for predicting CD86 expression by common SVM and LR methods, respectively. The Calibration curves showed that the predicted probabilities of 2 models have favorable consistency with actual probabilities (Figure 7B, 7D). The importance of 5 variables in 2 models was assessed and presented in Figure 7C, 7E.

Further, this study evaluated the prediction performance of 2 models in the training set. ROC analysis (Figure 8A, 8B) indicated that the SVM model has the same prediction ability as the LR model both in training (AUC: 0.769 vs 0.767) and validation sets (AUC: 0.707 vs 0.703). PR curve showed that AUC among the 2 models was similar (0.788 vs 0.797). DCA curve showed that both 2 models achieved more clinical net benefit. We also performed the ROC, PR, and DCA to evaluate the model’s stability for predicting the CD86 expression in the validation set (Figure 8A, 8B), finding that the performance of established models by SVM and LR in the validation set was similar with that in the training set. These results indicated the superiority and stable performance of the radiomics model established with 5 features to predict the CD86 expression. Further, we calculated the radiomics score (RS) according to 5 features in each patient, finding that RS in the CD86 high expression group was significantly higher than that in the low expression group in the 2 models (Figure 8C, all P<0.001).

To further explore the clinical value of the radiomics model in HGG, we analyzed the association of RS with prognosis in HGG. Due to the performance similarity of SVM and LR, we just only evaluate the clinical value of SVM model. Kapan-Meier analysis indicated that a higher RS was related to the shorter OS time (Figure 9A, P<0.001). We then integrated the RS and significant independent prognostic factors, and constructed a comprehensive nomogram model (Figure 9B), finding that the comprehensive model has good predicted probability (Figure 9C). ROC analysis indicated that the nomogram model can favorably predict the 1, 3, and 5-year survival of HGG patients (Figure 9D).