Bioinformatics in Cancer and Immune MicroenvironmentView this Special Issue
A Novel Thrombosis-Related Signature for Predicting Survival and Drug Compounds in Glioblastoma
Glioblastoma is the most common primary tumor in the central nervous system, and thrombosis-associated genes are related to its occurrence and progression. Univariate Cox and LASSO regression analysis were utilized to develop a new prognostic signature based on thrombosis-associated genes. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and HALLMARK were used for functional annotation of risk signature. ESTIMATE, MCP-counter, xCell, and TIMER algorithms were used to quantify immune infiltration in the tumor microenvironment. Genomics of Drug Sensitivity in Cancer (GDSC) was used for selecting potential drug compounds. Risk signature based on thrombosis-associated genes shows moderate performance in prognosis prediction. The functional annotation of the risk signature indicates that the signaling pathways related to the cell cycle, apoptosis, tumorigenesis, and immune suppression are rich in the high-risk group. Somatic mutation analysis shows that tumor-suppressive gene TP53 and oncogene PTEN have higher expression in low-risk and high-risk groups, respectively. Potential drug compounds are explored in risk score groups and show higher AUC values in the low-risk score group. A nomogram with valuable prognostic factors exhibits high sensitivity in predicting the survival outcome of GBM patients. Our research screens out multiple thromboses-associated genes with remarkable clinical significance in GBM and further develops a meaningful prognostic risk signature predicting drug sensitivity and survival outcome.
Glioma is the most common primary malignancy in the central nervous system, accounting for approximately 80% of primary malignant brain tumors in adults . According to the 2016 World Health Organization (WHO) category of the Central Nervous System Tumors, glioma can be classified into astrocytoma, oligodendroglioma, oligodendrocyte, and glioblastoma (GBM), including WHO grades I-IV based on its malignancy degree . The overall survival of GBM patients ranges from 1 to 15 years, with great individual differences . At present, the treatment of GBM is mainly based on maximum resection, adjuvant radiotherapy, and TMZ chemotherapy [4, 5]. Due to the highly invasive nature of GBM, it is difficult to completely remove the tumor by neurosurgery . The residual lesions are prone to resistance to radiotherapy and chemotherapy, thus leading to recurrence and malignant progression of GBM. With the development and application of sequencing technology and molecular diagnostic technology, a more objective and accurate tumor classification system has been established clinically at present [6, 7]. However, currently, known molecular markers can only partially explain the prognosis of GBM patients.
Venous thromboembolism (VTE), which comprises deep vein thrombosis (DVT) and pulmonary embolism (PE), occurs extensively in various cancers, including GBM. The incidence of VTE in GBM is high, and complications are easy to occur. VTE has been found in 17% of patients with GBM, and the cumulative probability of VTE was 14.3% at 6 months and 16.8% after 12 months . In addition to VTE, the intratumoral thrombosis in GBMs, which may cause ischemia or necrosis due to vascular obstruction, is also an essential factor affecting the prognosis of GBM patients. Risk factors for GBM-associated thrombosis include general clinical characteristics, such as advanced age, grade, tumor size, reduced activity, and thrombosis-associated biomarkers or genes. Recent studies have shown that the expressions of many thrombosis-associated proteins in GBMs, such as tissue factor (TF), human epidermal growth factor (VEGF), and D-dimer, are correlated with pathological tumor grade and poor prognosis [8, 9]. TF is a transmembrane glycoprotein that plays a vital role in thrombosis, and its expression is significantly upregulated in GBM [10, 11]. In addition, PTEN loss and hypoxia can further upregulate TF expression, thereby inducing angiogenesis and local necrosis of glioblastomas . Previous studies indicated that thrombosis-associated genes mediated GBM metastasis, thrombosis, and other biological processes, suggesting that thrombosis-associated genes changes may play an essential role in diagnosis and prognostic prediction of glioblastomas .
This study aims to identify prognostic thrombosis-associated gene signatures by integrating the transcriptional data of GBM in The Cancer Genome Atlas (TCGA) and the Chinese GBM Genome Atlas (CGGA) databases. We constructed an independent prognostic risk score model based on the mRNA expression of 13 thrombosis-associated genes that were phenotypically significantly associated with immune invasion and somatic mutation in GBM using univariate Cox and LASSO regression analysis. We believe that the findings may provide a theoretical basis for the molecular diagnosis and individualized treatment for GBMs.
2.1. Datasets Source
Gene expression profiles and corresponding clinical information were downloaded from the TCGA database (https://xena.ucsc.edu) and the CGGA database (https://www.cgga.org.cn). Thrombosis-associated genes were obtained from the literature, and the gene list is shown in Table S1. GBM samples with overall survival (OS) of less than 30 days in TCGA and CGGA databases were excluded. The TCGA GBM dataset is randomly divided into the training set (train) and the validation set (test). The total TCGA GBM dataset (sum) and CGGA GBM cohort are the two validation sets in this study.
2.2. Constitution and Validation of the Prognostic Risk Score Model
First, a univariate Cox proportional hazard regression analysis was performed using the “survival” package in R to identify the thrombosis-associated genes involved in prognosis. Genes with value <0.05 were considered to have significant prognostic potential. The least absolute shrinkage and selection operator (LASSO) regression was performed to further screen the genes with independent prognostic value . Based on the highest λ value selected through 1,000 cross-validations in the LASSO method, a set of prognosis genes and their LASSO coefficients (β) were obtained . Then, a gene-based survival risk assessment model was established with LASSO coefficients: , where N = 13, is the expression value of every 13 thrombosis associated genes, and is the corresponding LASSO regression coeﬃcient. Patients in the TCGA and CGGA datasets were divided into high- and low-risk groups according to the median risk score.
2.3. Functional Enrichment Analysis
To evaluate the biological pathway associated with the 13 thrombosis-associated genes, the “GSVA” package was used for gene set variation analysis (GSVA), including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and HALLMARK.
2.4. Immunological Function Analysis
According to the expression characteristics of thrombosis genes, cell components or cellular immune responses in the high-risk and low-risk groups were evaluated by the ESTIMATE , MCP-counter , xCell , and TIMER  algorithms. Heatmaps displayed the immune responses generated under different algorithms.
2.5. Somatic Mutation Analysis
Somatic mutations, including the SNVs, SNPs, and INDELs, were detected using the TCGA mutation data of both the high-risk group (n = 146) and low-risk group (n = 147) with the R package “maftools” . Fisher's exact test was used to determine the pattern of differential mutation. Genes with a value <0.05 were defined as differentially mutated genes. The cooccurrence and mutually exclusive mutations were identified using the CoMEt algorithm . Somatic mutation visualizations were generated using the R package “maftools.”
2.6. Prediction of Chemotherapy Response
According to the public pharmacogenomic databases, PRISM Repurposing dataset (PRISM, https://depmap.org/portal/prism/) and Cancer Therapeutics Response Portal (CTRP, https://portals.broadinstitute.org/ctrp), and according to a previous study, drug sensitivity (IC50) values were predicted by the R package “pRRophetic” .
2.7. Prognostic Model Based on Clinical Features and Risk Score
Univariate Cox and Multivariate Cox proportional hazard regression analysis was performed to identify independent prognostic risk factors, including the risk score and clinical characteristics (age, gender, subtype, IDH mutation, chemotherapy, and radiotherapy) using the “survival” package. The independent prognostic factors were then used to construct a nomogram chart and a calibration curve to evaluate and compare the predicted and actual OS probabilities for GBM patients at 1, 2, and 3 years. The nomogram chart and the calibration curves were constructed using the R package “RMS.”
2.8. Statistical Analysis
All statistical analyses were performed using SPSS 22.0 or R software. The normality of variables was tested using the Shapiro-Wilk normality test. For normally distributed variables, significant quantitative differences between and among groups were determined by a two‐tailed t-test or one‐way ANOVA, respectively. For nonnormally distributed variables, significant quantitative differences between and among groups were determined by a Wilcoxon test or a Kruskal-Wallis test, respectively. The chi-square test was used to analyze the correlation of the classified data. Kaplan-Meier survival curve and log-rank tests were used to detect the prognostic difference in different groups. The R package “survivalROC” was used to plot time-dependent receiver operating characteristic (ROC) curves and calculate the area under the curve (AUC) . values <0.05 were considered statistically significant.
3.1. Identification and Verification of Prognostic Risk Score Model Based on Thrombosis-Associated Gene Signature in GBM
The overall study design is shown in Figure 1. In total, 140 thrombosis-associated genes were identified from the literature, of which 124 genes expressed in both TCGA and CGGA datasets were chosen for subsequent analyses. To identify prognostic thrombosis-associate genes in GBM, LASSO regression analysis was performed using TCGA training, and 13 genes were screened out of 124 genes (Figures 2(a) and 2(b)). The LASSO regression coefficients of 13 thrombosis-associated genes are shown in Table S2. Then, a risk score model was established based on thrombosis-associated gene expression with LASSO coefficients: risk score = (0.1235∗ANXA2 mRNA expression + 0.1175∗C5 mRNA expression +0.0195∗CD59 mRNA expression + 0.0309∗CFH mRNA expression + 0.2343∗CR1 mRNA expression – 0.3472∗F13B mRNA expression + 0.0079∗FAP mRNA expression + 0.2369∗KLKB1 mRNA expression + 0.0247∗LBH mRNA expression + 0.0602∗PDGFA mRNA expression + 0.0104∗PLAT mRNA expression + 0.0419∗SERPING1 mRNA expression – 0.1703∗MASP1 mRNA expression). According to the median risk score, GBM patients in the TCGA training set were divided into high-risk and low-risk groups. The risk score distribution and survival status of patients are shown in Figure 2(c). Kaplan-Meier survival curves indicated that high-risk patients’ OS was significantly worse than that of low-risk patients (Figure 2(d)). A time-dependent ROC curve analysis evaluated the risk score model's predictive accuracy. The result showed that the AUC was 0.752, which proved that the risk score model had good accuracy and predictive ability within the TCGA training cohort (Figure 2(e)). Subsequently, the TCGA test set and TCGA sum set were used to verify the predictive prognosis ability of the risk score model. The risk score distribution and survival status of patients in the TCGA test set and TCGA sum set are shown in Figures 2(f) and 2(i). We also found that the OS of GBM patients in the high-risk group was significantly shorter than that in the low-risk group in the TCGA test set and TCGA sum set (Figures 2(g) and 2(j)). The AUC of the ROC curve in the TCGA test set and sum set were 0.658 and 0.704, respectively (Figures 2(h) and 2(k)).
Further verification in the CGGA dataset showed similar tendencies that the GBM patients in the high-risk group have a poorer prognosis compared to those in the low-risk group (Figure 2(m)). The AUC of the ROC curve was 0.714 (Figure 2(n)). The risk score distribution and survival status of patients in the CGGA dataset are shown in Figure 2(l). These findings revealed that the risk score model based on thrombosis-associated signature has an excellent survival predicting power of patients’ OS.
Additionally, we detected the expression of 13 thrombosis-associated genes in the TCGA and CGGA datasets. Heatmaps of 13 thrombosis-associated genes are shown in Figure S1, in which F13B and MASP1 were regarded as protective genes for expression level decreases with risk score gradually increases, and others were risky genes for growing together. To assess the association between 13 thrombosis-associated genes’ expression and patients’ OS, patients were divided into high-expression and low-expression groups based on the median of gene expression, and Kaplan-Meier curves were performed. As shown in Figure S2, the increased expression of MASP1 and F13B indicated an excellent survival in GBM patients in the TCGA dataset. In contrast, the high expressions of ANXA2, C5, CD59, CFH, CR1, FAP, KLKB1, LBH, and PDGFA predicted a poor prognosis in GBM patients. The association between IDH mutation status and the 13 thrombosis-associated genes was explored in the TCGA and CGGA datasets. ANXA2, C5, CFH, FAP, PDGFA, PLAT, and SERPING1 had higher expressions in the IDH WT tumors in the TCGA and CGGA datasets, which predicted worse survival (Figures S3A and S3B).
3.2. Functional Enrichment Analyses for Thrombosis-Associated Gene Signature
To explore the potential biological function of a prognostic thrombosis-associated gene signature in GBM, GSVA was performed. As shown in Figure 3(a), GO analysis results showed that the thrombosis-associated gene signature was mainly enriched in the positive regulation of Kappa B kinase NF-κB signaling, the regulation of extrinsic apoptotic signaling pathway, the regulation of low-density lipoprotein particle receptor binding, the regulation of response to cytokine stimulus, and so forth. KEGG analysis results indicated that the thrombosis-associated gene signature was primarily concentrated in the lysosome, the apoptosis, the glycosaminoglycan degradation, the focal adhesion, the Toll-like receptor signaling pathway, and so forth (Figure 3(a)). In addition, the HALLMARK analysis results displayed that the thrombosis-associated gene signature was also principally enriched in the apoptosis pathway, the P53 pathway, the TNFα signaling via NF-κB, the epithelial-mesenchymal transition pathway, and so forth (Figure 3(a)).
In addition, to investigate the correlation between the expression of thrombosis-associated gene signature and the known signature, GBM patients in the TCGA dataset were divided into high-risk and low-risk groups based on the median risk score. We found significant differences in APM expression, cell cycle regulation, EMT2, FGFR3 related signature, histones, immune checkpoint, nucleotide excision repair, and Pan_F_TBRs signature between the low-risk group and the high-risk group (Figure 3(b)). Differential expression analysis of m6A-related genes between the low-risk group and the high-risk group showed that the expressions of CBLL1, ELAVL1, LRPPRC, RBM15B, WTAP, YTHDC2, YTHDF1, YTHDF2, YTHDF3, and ZC3H13 in the high-risk group were significantly higher than those in the low-risk group in the TCGA GBM dataset (Figure 3(c)). Notably, the expression of WTAP in the high-risk group was also significantly higher than that in the low-risk group in the CGGA dataset (Figure S4).
3.3. Immunological Function Analysis for Prognostic Thrombosis-Associated Gene Signature
To verify the relationship between prognostic thrombosis-associate gene signature and immune responses, we analyzed the correlation between risk score and immune factors in the TCGA GBM dataset based on ESTIMATE, MCP-counter, TIMER, and xCell. As shown in Figure 4, the risk score showed a negative correlation with Tumor Purity and a positive correlation with Stromal Score, Immune Score, and ESTIMATE Score. Moreover, correlation analysis between risk score and immune cell subpopulations based on the MCP-counter revealed that risk score was negatively correlated with cytotoxic lymphocytes and NK cells but positively associated with monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts (Figure 4). Additionally, correlation analysis between risk score and immune cell subpopulations based on TIMER (Figure 4) and xCell (Figure 4) also indicated that risk score has a correlation with B cells, CD8+ cells, CD4+ cells, neutrophil, macrophage, and dendritic cells (DC). Similar results were obtained in the CGGA dataset (Figure S5).
3.4. Somatic Mutation Analysis for Prognostic Thrombosis-Associated Gene Signature
To identify differences in somatic mutation between the high-risk and low-risk groups, we analyzed the somatic mutation in the TCGA GBM dataset. In the low-risk group, 16 genes were mutated in more than 10% of the samples, while only 12 genes met the criteria in the high-risk group, of which 9 genes overlapped (Figure 5(a)). The top 50 with the highest mutation frequency in the low-risk and high-risk groups are shown in Figure 5(a). Interestingly, TP53 , TTN , PTEN , and EGFR  occupied the top four positions in both low-risk and high-risk groups, and they are interacting with each other to regulate various biological processes related to GBM, suggesting that they may be involved in tumor deterioration. In addition, currently recognized GBM prognostic genes (IDH1, TP53, ATRX, NUP16, TIAM2, NEK10, and ABCA1) showed significant differences in mutation frequency between the high-risk and low-risk groups (Figure 5(b)). Next, we studied cooccurrence and exclusive mutations in the 25 most common mutated genes using the CoMEt algorithm. Except for the prevailing mutually exclusive mutation landscape, there were three unique gene pairs. Two genes exhibited cooccurrence, including HYDIN and FLG, AHNAK and SYNE1, and NF1 and LRP2 (Figure 5(c)), suggesting their probable complementary effect in the same pathway. More interestingly, some genes had differential mutation frequencies between the two groups.
3.5. Identification of Potential Therapeutic Agents for High-Risk Score Patients
To determine potential therapeutic agents for GBM patients with a high-risk score, PRISM and CTRP-derived drug response data were analyzed. We first screened candidate drugs that responded differently between the high-risk and low-risk groups to identify compounds with higher drug sensitivity in GBM patients with a high-risk score. The candidate drug meets two criteria: (1) Drugs respond differently between high-risk and low-risk groups to identify highly drug-sensitive compounds in high-risk GBM patients. (2) Spearman correlation analysis between the estimated AUC values of a candidate drug and risk score was used to select compounds with negative correlation coefficient (Spearman’s r < −0.25 for CTRR or Spearman’s r < −0.35 for PRISM). These analyses yielded nine PRISM-derived compounds (AGM-232, AS-703026, AZD8330, Cobimetinib, Dabrafenib, GDC-0152, Napabucasin, Narasin, and TAK-733) (Figures 6(e) and 6(f)) and three CTRP-derived compounds (Birinapant, Dasatinib, and RITA) (Figures 6(g) and 6(h)). All these compounds had lower estimated AUC values in a high-risk group and a negative correlation with risk score (Figures 6(e)–6(h)).
3.6. Clinical and Molecular Features of Low-Risk and High-Risk GBM Patients
We performed risk stratification analysis on various clinicopathological features (age, IDH mutation, radiotherapy, and chemotherapy) with different groups to evaluate the prognostic value of the risk score model. Older age is a known risk factor for malignant glioblastomas. As shown in Figure 6(a), both GBM patients below 65 years of age and GBM patients aged 65 years or older with a high-risk score had a worse prognosis than those with a low-risk score in the TCGA GBM cohort. A similar result was obtained in the CGGA dataset (Figure S6A). IDH mutation is currently recognized as a molecular marker for prognostic prediction of GBM patients [11, 13, 20]. To determine whether the risk score can be used as a prognostic indicator independent of IDH mutation, we performed a stratified analysis of patients in the high-risk and low-risk groups based on IDH mutation status. In the TCGA GBM and CGGA datasets, IDH wild-type patients with high-risk scores had a worse prognosis than those with a low-risk score, and the same trend was also found in IDH mutation patients (Figures 6(b) and S6B). To investigate the risk score based on thrombosis-associated signature concerning prognosis among GBM patients with radiotherapy and chemotherapy, Figures 6(c) and 6(d) show that patients in the high-risk group without radiotherapy or chemotherapy indicated the worst prognosis. In contrast, patients in the low-risk group with radiotherapy or chemotherapy showed the best forecast. Similar results were obtained in the CGGA dataset (Figures S6C and S6D). The above results proved that the risk score model based on thrombosis-associated genes could serve as an independent prognostic factor for GBM patients.
3.7. Construction and Evaluation of the Clinical-Featured Risk Model of Glioblastoma
Then, we built a nomogram to predict the 1-year, 2-year, and 3-year overall survival with the risk score based on thrombosis-associated signature and clinical factors (age, gender, IDH mutation status, radiotherapy, and chemotherapy) (Figure 7(a)). The calibration chart showed that the nomogram performed well at predicting the 1-year, 2-year, and 3-year OS for the TCGA GBM and CGGA cohorts, and the predicted 2-year OS was approximated to the actual 2-year OS (Figures 7(b) and 7(e)). Kaplan-Meier curves showed a significant difference in OS between the high-risk and low-risk groups based on the clinical-featured risk model (Figures 7(c) and 7(f)). Moreover, ROC curve analysis was performed to estimate the predictive accuracy of the clinical-featured risk model. As shown in Figures 7(d) and 7(g), the AUC of ROC curves were 0.765 and 0.785 in the TCGA and CGGA cohorts, respectively, suggesting that the clinical-featured risk model had better accuracy and stability than the risk score model based on thrombosis-associated gene signature. In conclusion, the above results showed that the clinical-featured risk model had strong prognostic capabilities for GBM patients.
Although significant advances have been made in surgical treatment, radiotherapy, and chemotherapy of GBM, the prognosis of GBM patients is still poor. Recent studies based on genetic and epigenetic biomarkers have improved the accuracy of prognosis prediction of GBM patients and optimized the treatment strategy of GBM [1, 28–33]. For example, IDH mutation, 1p/19q deletion, TP53, and ATRX were widely identified as the prognostic markers for GBM. Currently, thrombosis-associated genes have also been found to play vital roles in tumor progression, metastasis, tumor immunity, and therapeutic resistance. This study aimed to develop a prognostic model for GBM based on thrombosis-associated genes and ultimately help physicians estimate patient outcomes and design appropriate treatment strategies.
This study obtained a novel gene signature composed of thirteen thrombosis-associated genes by analyzing the TCGA GBM and the CGGA datasets. Based on the LASSO regression coefficients of 13 thrombosis-associated genes, we constructed a risk score model that can be used as a prognostic indicator for GBM patients in the TCGA and CGGA datasets. The signature can stratify patients according to their risk score. As results showed, high-risk patients have worsened survival in these datasets. We next comprehensively annotated the risk score’s biological and immune-related functions. Several potential drug compounds related to the risk score were screened out. We finally determined the risk score as an independent prognostic factor in the nomogram incorporating other elements (age, gender, IDH mutation, radiotherapy, and chemotherapy). Among the 13 identified genes, CD59 has long been recognized as the complement membrane regulator of malignant GBM . CFH, a circular RNA complement factor, has promoted GBM progression .
Similarly, hypoxia-induced LBH, a highly conserved transcription cofactor, participates in embryonic development and promotes tumor progression of GBM , and hypoxia upregulates the expression of PLAT in GBM cells . Besides, FAP, overexpressed in several kinds of brain cancers, has been reported to be an excellent target for immunotherapy . MASP-1, a significant mediator in the lectin complement signaling pathway of the innate immune response, has been overexpressed in GBM cell lines . Moreover, a recent study revealed that PDGFA/PDGFRα-regulated GOLM1 critically enhances the tumor progression ability through activating AKT signaling in GBM .
Among the expression level and clinical correlation of thirteen candidate genes in patterns, SERPING1and ANAX2 are of particular interest. SERPING1 (also known as C1-INH) is a plasma protein-encoding C1 inhibitor that plays a vital role in component regulation, coagulation, and fibrinolysis. SERPING1 deficiency has been linked to the development of hereditary angioedema, sepsis, and pancreatic cancer. In line with previous studies, our study also found that the SERPING1 gene is highly expressed in GBM and suggests a poor prognosis in patients with GBM. However, the exact role of SERPING1 in GBMs requires further investigation. ANXA2 is a member of the annexin family encoding a calcium-dependent phospholipid-binding protein that regulates cellular growth and signals transduction pathways. ANXA2 overexpression has been connected with various cancer growths, invasion, metastasis, and drug resistance. Previous evidence indicated that miR155HG was overexpressed and acts as a tumor suppressor to upregulate ANXA2 forming a loop that could promote the GBM’s malignant progression. In this study, we found that the overexpression of ANXA2 predicted a worse prognosis in GBM.
Given the vital roles of the 13 thrombosis-associated genes in GBM, a tumorigenic functional annotation is likely to be observed in the corresponding thrombosis-associated gene signature. As expected, the prognostic signature was enriched in Kappa B kinase NF-κB signaling, the regulation of extrinsic apoptotic signaling pathway, the apoptosis, the focal adhesion, and the Toll-like receptor signaling pathway, the P53 pathway, and the TNFα signaling via NF-κB, all of which were important signaling pathways regulating tumor proliferation and tumor progression. Notably, the prognostic signature was also associated with tumor metastasis. A close connection between m6A methylation regulators and the prognostic signature was observed, and the high score group had higher expression of writer WTAP. This also indicated the hazardous role of the prognostic signature in GBM patients.
Tumor microenvironment (TME) has been proposed to mediate the progression of GBMs critically [41–43]. The various immune infiltrating cells within the TME play an essential role in the crosstalk between GBM cells and TME [41–45]. In this study, multiple immune regulatory cells, including M2 macrophages and fibroblasts, were more expressed in the high-score group. The consistent findings proved that the thrombosis-associated gene signature was involved in an immunosuppressive microenvironment. Since GBMs are highly heterogeneous among individuals, we analyzed the characteristics of somatic mutations in the high-risk and low-risk groups. Our results proved that the tumor-suppressive gene, TP53, was mutated more frequently in the low-risk score group, while the oncogenic gene, PTEN, was mutated more commonly in the high-risk score group. This was by the previous finding, which further proved the prognostic value of the thrombosis-associated gene signature.
The current treatment modality for gliomas, especially GBM, is limited and dismal. To a large extent, novel therapeutic options would help promote the clinical management of GBM patients. In the analysis of potential drug compounds in risk score groups, all identified drugs from PRISM and CTRP have lower AUC values in the high-risk score group, suggesting that GBM patients in the high-risk group have increased sensitivity to the compounds-based therapy. PRISM-derived and CTRP-derived compounds were highly correlated with the thrombosis-associated signature. The above findings indicated that these compounds could be potentially applied in treating GBM patients. These identified compounds could be explored in a clinical trial for more robust verification of their therapeutic efficacy. Finally, a nomogram incorporating risk score, age, gender, IDH mutation, radiotherapy, and chemotherapy was constructed, which showed high sensitivity in predicting survival outcomes of GBM patients.
Although the new signature established by this research provides new biomarkers for thrombosis prevention in GBM patients and gives the foundation for developing anticoagulation therapy, there are still many deficiencies. For example, our risk score model is based on retrospective data and needs future research. Furthermore, apart from the two genes in our study, single genes’ potential function and mechanism need to be further explored.
In sum, an independent prognostic risk score model based on the mRNA expression of 13 thrombosis-associated genes that were phenotypically significantly associated with immune invasion and somatic mutation in GBM was constructed using LASSO regression. We believe that the findings may provide a theoretical basis for the molecular diagnosis and individualized treatment for glioblastomas.
Conflicts of Interest
The authors declare no conflicts of interest.
Wen-Jing Zeng, Zhi-Cheng Gong, and Quan Cheng conceived and designed the study. Wen-Jing Zeng, Hao Zhang, He Li, and Yu-Fang Cao wrote the manuscript. Wantao Wu, Zaoqu Liu, Peng Luo, and Jian Zhang edited the language. Quan Cheng and Hao Zhang crafted the figures and tables. The authors confirm that relevant guidelines and regulations were followed in all methods.
This study was supported by the National Natural Science Foundation of China (no. 81903663), Hunan Provincial Natural Science Foundation of China (nos. 2020JJ5944 and 2018JJ3838), Hunan Provincial Health Committee Foundation of China (202204044869), and Xiangya Hospital Central South University Postdoctoral Foundation.
Figure S1: heatmap depicting the expression difference of thirteen prognostic thrombosis-associated genes in TCGA training set, TCGA test set, TCGA sum set, and CGGA sum set. The column represents each patient and the row represents 13 genes’ expression level. Figure S2: the thirteen prognostic thrombosis-associated genes have independent prognostic value in TCGA dataset and CGGA dataset. ((A)–(M)) Kaplan-Meier curves analysis correlation of OS among different risk groups and 13 genes. Figure S3: expression difference of IDH mutation status in the thirteen prognostic thrombosis-associated genes in (A) TCGA and (B) CGGA, respectively. Figure S4: functional enrichment analyses of thrombosis-associated gene signature in gliomas. (A) GO, KEGG, and HALLMARK analyses for thrombosis-associated prognostic genes in CGGA. (B) The expression of known signature in high-risk and low-risk groups in CGGA. (C) The expression of m6A-related genes in high-risk and low-risk groups in CGGA. GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; m6A, N6-methyladenosine. Figure S5: immunological function analysis of thrombosis-associated gene signature in gliomas. Heatmap shows the differential cellular immune responses between high-risk and low-risk groups analyzed by (A) ESTIMATE, MCP-counter, and TIMER algorithms and (B) xCell algorithm in CGGA. Figure S6: survival analysis of high-risk and low-risk groups with different clinicopathological factors and the prediction of chemotherapy response. (A) Kaplan-Meier survival curves of patients in the high-risk and low-risk groups with aged 65 years or older and those below 65 years of age in the CGGA glioma cohort. (B) Kaplan-Meier survival curves of patients in the high-risk and low-risk group with IDH mutation or wild-type IDH in CGGA glioma cohort. ((C) and (D)) Kaplan-Meier survival of patients in the high-risk and low-risk groups receiving radiotherapy and chemotherapy in the CGGA glioma cohort. Table S1: thrombosis-associated genes list. Table S2: Lasso regression coefficients of 13 thrombosis-associated genes. (Supplementary Materials)
Y. Zhang, H. Li, W. Zhang, Y. Che, W. Bai, and G. Huang, “LASSO based CoxPH model identifies an 11lncRNA signature for prognosis prediction in gastric cancer,” Molecular Medicine Reports, vol. 18, pp. 5579–5593, 2018.View at: Google Scholar
H. Zhang, Z. Chen, Z. Wang et al., “Correlation between APOBEC3B expression and clinical characterization in lower-grade gliomas,” Frontiers in Oncology, vol. 11, Article ID 625838, 2021.View at: Google Scholar
A. Mäenpää, S. Junnikkala, J. Hakulinen, T. Timonen, and S. Meri, “Expression of complement membrane regulators membrane cofactor protein (CD46), decay accelerating factor (CD55), and protectin (CD59) in human malignant gliomas,” American Journal of Pathology, vol. 148, pp. 1139–1152, 1996.View at: Google Scholar
O. H. Minchenko, A. P. Kharkova, K. I. Kubaichuk, D. O. Minchenko, N. A. Hlushchak, and O. V. Kovalevska, “Effect of hypoxia on the expression of CCN2, PLAU, PLAUR, SLURP1, PLAT and ITGB1 genes in ERN1 knockdown U87 glioma cells,” Ukrainian Biochemical Journal, vol. 86, no. 4, pp. 79–89, 2014.View at: Publisher Site | Google Scholar
L. M. Ebert, W. Yu, T. Gargett et al., “Endothelial, pericyte and tumor cell expression in glioblastoma identifies fibroblast activation protein (FAP) as an excellent target for immunotherapy,” Clinical & Translational Immunology, vol. 9, Article ID e1191, 2020.View at: Publisher Site | Google Scholar
V. Pagliara, M. Parafati, A. Adornetto et al., “Dibutyryl cAMP- or Interleukin-6-induced astrocytic differentiation enhances mannose binding lectin (MBL)-associated serine protease (MASP)-1/3 expression in C6 glioma cells,” Archives of Biochemistry and Biophysics, vol. 653, pp. 39–49, 2018.View at: Publisher Site | Google Scholar
N. Zhang, H. Zhang, Z. Wang et al., “Immune infiltrating cells-derived risk signature based on large-scale analysis defines immune landscape and predicts immunotherapy responses in glioma tumor microenvironment,” Frontiers in Immunology, vol. 12, Article ID 691811, 2021.View at: Publisher Site | Google Scholar