Prognostic signature detects homologous recombination deficient in glioblastoma

Introduction

Glioblastoma (GBM) is regarded as a catastrophic brain illness due to its low incidence, high fatality rate, fast tumor development, and high genetic instability. Despite multimodality treatment involving temozolomide (TMZ)-based chemotherapy and surgical resection, GBM continues to exhibit very high recurrence and mortality rates, so the therapy choices for GBM remain restricted and exceedingly difficult (1). With the growing understanding of the molecular pathophysiology and biology of GBM over the past several years, tremendous progress has been achieved in improving clinical outcomes. However, the overall survival rate of GBM is substantially below our expectations, and the 5-year survival rate is below 10% (2).

Not only are there more genetic mutations in tumors with significant genomic instability, but the expression of drug-resistant proteins in tumor cells is also increased (3), ultimately encouraging the development of GBM. More than one-third of GBM patients have O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation (4), which results in (5) structural alterations that inhibit transcription factor binding, leading to gene silence and loss of DNA repair activity. Therefore, MGMT methylation is not only a crucial pathogenic component but also a therapeutic target and prognosis indicator (6,7). In a manner analogous to MGMT methylation, homologous recombination deficiency (HRD) can likewise be a significant cause of genomic instability. In recent years, HRD has been identified in several forms of cancer, serving as a powerful prognostic biomarker. Molecular characterization of HRD in ovarian malignancies predicts more effective poly adenosine diphosphate (ADP)-ribose polymerase (PARP) inhibitor therapy and longer life (8), but identification of HRD in adrenal cortical carcinoma predicts a poorer clinical prognosis in comparison to non-HRD patients (9). The genomic alterations produced by HRD, such as identifying genetic mutations, insertion/deletion patterns, chromosomal structural abnormalities, and gene copy number variations, are specific, quantifiable, and permanent alterations that serve as the theoretical foundation for the current clinical detection of HRD (10). HRD as a functional deficiency in homologous recombinant DNA repair may result in a “genomic scar” consisting of loss of heterozygosity (LOH), telomere allele imbalance (TAI), and large-scale state transitions (LST). LOH, TAI, and LST all have distinct definitions and can partially characterize the HRD state of cells (11). While each score has clinical significance on its own, the HRD score, which is derived from the three markers (LOH + TAI + LST), is a more accurate predictor of HRD than any of the individual scores (12,13). In other malignancies, such as ovarian cancer, the association between HRD and tumors has been extensively explored. However, the influence of HRD on prognosis and the tumor immune microenvironment (TIME) in GBM remains unknown.

In recent years, processing single nucleotide polymorphism (SNP) or microarray data by computer technology in cancer diagnosis or prognosis is one of the most important applications. Proper selection of appropriate genes can reduce the difficulty of data analysis and improve the accuracy of analysis results. For example, a social network analysis-based gene selection approach (14) or a graph theoretic-based gene selection (15) can maximize the correlation and minimize the redundancy of selected genes for cancer diagnosis. In this study, we assessed the prognostic role of HRD by analyzing the HRD-related gene expression signatures through machine learning in GBM.

We noticed in this study that the HRD score had a significant predictive influence on GBM. When the HRD score was 17, individuals with GBM had a greater chance of survival. The data of differentially expressed genes (DEGs) related to HRD score were utilized to develop a 14-gene risk scoring model, which might aid in stratifying patients with various survival outcomes. The TIME analysis revealed that the signature of T cells was more prevalent in patients with a favorable prognosis who had a low-risk score. It is envisaged that the outcomes of this study would bring fresh insights into the individualized therapy of GBM. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2077/rc).

Methods

Data collection and processing

The original whole exome sequencing (WES) sequencing data, gene expression profile data and clinical data were downloaded from the publicly available the The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov), the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn), and Gene Expression Omnibus (GEO) database (16,17). Transcriptome expression array data were downloaded from UCSC Xena database (https://xenabrowser.net/datapages/?hub=https://tcga.xenahubs.net:443). In total, 567 and 657 GBM samples with complete clinical data from the TCGA and CGGA databases were obtained. The clinicopathological statistics of the samples from the TCGA are shown in table available at https://cdn.amegroups.cn/static/public/tcr-23-2077-1.xlsx. The clinicopathological statistics of the samples from the CGGA are shown in table available at https://cdn.amegroups.cn/static/public/tcr-23-2077-2.xlsx. We used VarScan2 to make somatic mutation calls, including SNPs and insertion-deletions. We annotated the mutation site and extracted the exon variation by using analysis of variance (ANNOVAR), including single nucleotide variation (SNV), non-frameshift insertion, non-frameshift deletion, frameshift insertion and frameshift deletion. The fraction of genome altered (FGA) score was defined as percentage of copy number altered chromosome regions out of measured regions. Microsatellite instability (MSI) sensor score was defined as MSI detection using paired tumor-normal sequence data. FGA score and MSI sensor score were obtained from the cBioPortal database (http://www.cbioportal.org/study?id = rca_tcga_pan_can_atlas_2018). The receiver operating characteristic (ROC) curve was generated by R software package time ROC. The in vitro experiment data of cell lines were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC2) database (www.cancerRxgene.org) and the Project Score database (https://score.depmap.sanger.ac.uk/). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

HRD score analysis

LOH was defined as the number of counts of LOH greater than 15 MB and less than the entire chromosome length (18). LST were defined as a chromosomal break site between two adjacent regions (the length of both regions is equal to or equal to 10 MB, and the distance between regions was less than 3 MB). The total number of tumor genome cut-off points can be used to describe the genomic instability (19). TAI was defined as the number of chromosomal segment with an allele imbalance extending to one of the sub telomeres but not exceeding the centromere and greater than 11 MB (20). The comprehensive scores of TAI, LST, and LOH were defined as the HRD score. The HRD score of each patient from the TCGA is shown in table available at https://cdn.amegroups.cn/static/public/tcr-23-2077-1.xlsx.

Difference analysis and enrichment analysis

We analyzed the differences of different molecular subtypes through R-package Limma (21). Fold change >1 and adj. P<0.05 were set as the cutoffs to screen for DEGs. In addition, in order to analyze the enrichment of different molecular subtypes in different pathways, the cp.kegg.v7.0.symbols.gmt gene set (22) was used as the reference gene set for gene set enrichment analysis (GSEA). The pathways with P<0.05 and false discovery rate (FDR) <0.25 threshold were considered as significantly enriched.

Construction and evaluation of a prognostic scoring system

There were 491 samples with microarray-seq data samples from the TCGA dataset and 656 from the CGGA dataset grouped into a discovery cohort and a verification cohort respectively. The R package “glmnet” was used to perform least absolute shrinkage and selection operator (LASSO) Cox regression model analysis (23). In LASSO model construction, 10-fold cross-validation was used to find the optimal value of the penalty parameter λ, and then the prognostic genes with regression coefficients were selected based on the optimal λ value. The risk score of each patient was calculated via the non-zero coefficient in Cox regression analysis.

GSEA

GSEA was performed using GSEA software from the Broad Institute (MIT, Cambridge, MA, USA) to identify differential signaling pathways in different groups (24). The normalized enrichment score was calculated for each gene set. GSEA results with a nominal P<0.05 were considered significant.

Tumor microenvironment (TME) cell infiltration

The proportion of immune cells was quantified using the CIBERSORTx algorithm. For CIBERSORTx, the normalized gene expression data were uploaded to the web portal by using LM22 signatures and 1,000 permutations (25).

Statistical analysis

The package “survival” was used to calculate and plot Kaplan-Meier survival curves. Mutect2 is a somatic variant caller that uses local assembly and realignment to detect SNVs and indels (26). Wilcoxon test was used in statistical analysis (27).

Results

Genomic scar-based HRD scores and the mutation landscape in GBM

To investigate the potential relevance of HRD in GBM, our discovery cohort consisted of 567 GBM patients from the TCGA database. We computed HRD scores using WES data and utilized ROC curves to evaluate ideal HRD score cutoffs about GBM survival impact. When sensitivity and specificity were given equal weight, an HRD score of 17 was proposed as the threshold. In this discovery cohort, there was a statistically significant difference in survival probability between groups when utilizing an HRD score of 17 (log-rank test, P=0.002). Notably, HRD ratings were better predictive of survival outcomes with extended follow-up durations (Figure 1A).

Figure 1 Clinical significance of HRD in the TCGA-GBM cohort. (A) Left, Kaplan-Meier survival analysis revealed differences in OS between HRD and HRP patients in the TCGA discovery cohort. Right, 1-, 3-, and 5-year OS curves for patients in the TCGA discovery cohort. (B) Summary of the most prevalent genomic alterations in 567 GBMs. The mutational matrix shows Missense_Mutations (green), Nonsense_Mutations (red), Frame_shift_Deletions (blue), Splice_site mutations (orange), Frame_shift_Insertions (purple), In_Frame_Dels (yellow), In_Frame_Ins (brownness), Translation_Strat_Site mutations (light brown), Nonstop_Mutations (light blue). (C) Histogram of gene mutation frequency in the two subgroups (Wilcoxon signed-rank test, *, P<0.05; ***, P<0.001). (D) Violin plot of fraction genome altered in the HRD group and the HRP group (Wilcoxon signed-rank test, ****, P<0.0001). (E) Violin plot of MSI in the HRD group and the HRP group (Wilcoxon signed-rank test, **, P<0.01). HRD, homologous recombination deficiency; HRP, homologous recombination proficient; AUC, area under the curve; HR, hazard ratio; CI, confidence interval; FGA, fraction of genome altered; MSI, microsatellite instability; TCGA-GBM, The Cancer Genome Atlas-Glioblastoma; OS, overall survival.

In the study population, an integrated analysis of WES data was done. As shown in Figure 1B, similar to the findings of other cohort studies, the five most frequently mutated genes, PTEN, TP53, TTN, EGFR, and MUC16, all exhibited mutation rates of more than 15%. The majority of these changes were missense mutations, followed by nonsense mutations. The other mutation types (frameshift mutations, in-frame insertions and deletions, splice sites, and promoter alterations) occurred less often than 10% of the time. PTEN (Chi-squared test, P<0.001) and EGFR (Chi-squared test, P<0.001) showed considerably greater mutation frequency in the homologous recombination proficient (HRP) subgroup, but isocitrate dehydrogenase 1 (IDH1) (P<0.001), ATRX (Chi-squared test, P<0.001) and TP53 (Chi-squared test, P<0.05) mutations were significantly more prevalent in the HRD subgroup (Figure 1C). Next, we evaluated the additional genetic differences between the two categories, such as FGA and MSI. The hallmarks of genomic instability were significantly higher in the HRD group than in the HRP group (Wilcoxon signed-rank test, P<0.001, Figure 1D,1E).

Discovering the HRD transcriptome signature in GBMs and constructing a genetic prognostic model

To comprehend HRD-specific transcriptome markers in GBMs, we evaluated GBM gene expression profile data (microarray-seq data for 491 cases were available in the TCGA cohort). The limma analysis of gene expression profiling data indicated that 846 genes were differently expressed between HRD and HRP patients (|logFC| >0.58; FDR 0.05). Specifically, 300 DEGs were up-regulated and 546 were down-regulated among the 846 DEGs (Figure 2A).

Figure 2 Analyze the properties of the transcriptome and develop a genetic prognostic model. (A) DEGs between HRD and HRP subtypes, the red triangles represent differentially up-regulated genes and the green triangles represent differentially down-regulated genes. (B) The trajectory of each independent variable: the horizontal axis represents the log value of the independent variable lambda, and the vertical axis represents the coefficient of the independent variable. (C) The confidence interval under each lambda. (D) Kaplan-Meier survival analysis revealed differences in prognosis among the two subgroups. (E) In the TCGA discovery cohort, the risk scores of the cases were arranged in order from low to high, the existential state of the cases and the expression heat maps of the 14 genes. HR, hazard ratio; CI, confidence interval; DEGs, differentially expressed genes; HRD, homologous recombination deficiency; HRP, homologous recombination proficient; TCGA, The Cancer Genome Atlas.

To enable further validation, we selected 30 DEGs based on the requirements |logFC| >1.5 and FDR <0.05 to develop a predictive model. We utilized the R software package glmnet for LASSO Cox regression analysis and incorporated survival time, survival status, and gene expression data to develop a model for predicting the prognosis of GBM. LASSO regression was utilized to address the multicollinearity issue encountered during regression analysis and to minimize the number of genes in the prognostic model. Along with the progressive development of lambda, the coefficients of the independent variable trended toward 0 and increased gradually (Figure 2B). Next, the best model was constructed using a 10-fold cross-test and a confidence interval under each lambda (Figure 2C). Finally, a prognostic model was developed when lambda =0.03 and 14 genes (CHL1, CCL2, PDPN, POSTN, CHI3L1, IL8, FABP5, IGFBP3, LGALS3, EMP3, TAGLN, DCX, GPR37, C21orf62) were selected, the expression levels of the remaining 13 genes being significantly decreased in the HRD subgroup except for DCX. The following is the formula:

RiskScore=0.00840437618881702∗IGFBP3+0.0174609305158104∗LGALS3+0.0218743131736606∗POSTN+0.0867281386871314∗CHL1−0.0530308228032681∗TAGLN+0.0133252395475708∗DCX+0.0137115579495011∗PDPN+0.0105218535965084∗IL8−0.0105297118325562∗FABP5+0.0281053616847262∗CCL2+0.0951958995347424∗EMP3+0.0220799485270586∗CHI3L1−0.0333495503006732∗C21orf62−0.0362471779642503∗GPR37[1]

Patients in the TCGA cohort may be separated into two categories with significantly different survival outcomes based on the prognostic model (Figure 2D). It was not difficult to determine that the proportion of patients with a favorable prognosis according to the prognostic model was much greater than that according to the HRD score, 24.8% vs. 11.5%. In other words, the 14-gene prognostic model outperformed the HRD score alone in prognostic stratification, as it can identify not only individuals with an HRD score of 17 but also those whose gene transcriptome patterns indicated HRD status.

The samples were then inserted into the risk score calculation method. Patients with a higher survival rate had a lower risk score than those with a higher risk score. Moreover, the expression trend of these 14 genes was considerably lower in the subgroup with the lowest risk score, except DCX (Figure 2E).

Analyzing the genomic and TIME characteristics of patients with GBM based on the prognostic models

By studying the peculiarities of each patient’s TIME, we endeavored to determine the potential causes driving the varying prognoses. The median HRD score of patients in the low-risk category was considerably greater than that of individuals in the high-risk cohort (Figure 3A). The GSEA revealed significant enrichment of DNA damage repair (DDR) signaling pathways in the low-risk score subgroup, including the homologous recombination (HR) pathway and the non-homologous end joining (NHEJ) pathway (Figure 3B), which infers that the 14-gene prognostic model can accurately assess the HRD status of patients.

Figure 3 The prognostic categorization of glioblastoma patients based on the prognostic model. (A) The violin diagram of the two subgroups’ HRD scores. (B) GSEA analysis of the most enriched pathways in the subgroup with the low-risk score relative to the subgroup with the high-risk score. (C) Immune cell infiltration signal enrichment violin diagram for the two TCGA subgroups in the tumor microenvironment (Wilcoxon signed-rank test, · and - P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001). (D) The heat map of the TCGA cohort’s immunity infiltrates cell signal enrichment. KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer; HRD, homologous recombination deficiency; GSEA, gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

HRD has been demonstrated to change TIME and enhance the clinical outcome of breast cancer patients treated with immune checkpoint inhibitors (28). HRD induces the collapse of the replication fork to activate the cGAS-STING pathway associated with tumor immune infiltration, which promotes the production of IF-β to initiate anti-tumor innate immune responses (29). To investigate further the influence of HRD on TIME in GBM, two subgroups of tumor-infiltrating immune cells were evaluated using gene expression deconvolution methods (CIBERSORTx) (Figure 3C). We discovered that the signature of T cells was considerably enhanced in the subgroup with the lowest risk score (Figure 3D), suggesting that there may be a greater invasion of γδT cells.

Previous in vitro study showed that γδT cell-mediated tumor death is sensitive in GBM and that the killing effect can be augmented by zoledronate to cure and prolong the lives of patients (30). The profile of γδT cells was considerably enriched in the subgroup with a low-risk score and a favorable prognosis, consistent with earlier findings.

The prognostic model was validated in the CGGA cohort

We evaluated the 14-gene predictive model using GBM patients from the CGGA cohort to establish the model’s robustness and broad applicability. The prognostic model was also able to split the CGGA cohort into two subgroups, with 361 cases allocated to the low-risk-score subgroup and 295 to the high-risk-score subgroup. The findings of the survival analysis are depicted in Figure 4A; the low-risk-score subgroup had better outcomes than the high-risk-score subgroup, consistent with the TCGA cohort results. The GESA results indicated that the HR pathway and NHEJ pathway were highly active in the low-risk-score group. These pathways are essential for coordinating efficient cellular genomic instability repair responses (Figure 4B).

Figure 4 The prognostic model was validated in the CGGA cohort. (A) A Kaplan-Meier survival analysis for the two groups. (B) GSEA analysis of the major enriched pathways of the low-risk subgroup in the CGGA cohort compared to the high-risk subgroup. (C) Immunity cell infiltration signal enrichment violin diagram for the two CGGA subgroups in the tumor microenvironment (Wilcoxon signed-rank test, · and - P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001). HR, hazard ratio; CI, confidence interval; KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer; CGGA, Chinese Glioma Genome Atlas; GSEA, gene set enrichment analysis.

Analysis using CIBERSORTx revealed that the signature of T cells was similarly considerably enriched in the subgroup with the lowest risk score (Figure 4C). Combining the findings of CIBERSORTx analysis and survival analysis revealed that one of the explanations for the better prognosis indicated by our prognostic model may be the enrichment of T lymphocytes, which have a powerful tumor-killing impact.

The prognostic model may serve as an independent predictor of the prognosis for GBM patients

In GBM tumor genomes, abnormal gene expression and mutations are prevalent, and some of them have critical prognostic effects and are employed therapeutically as prognostic markers (31). It is known that GBM patients with the IDH1/2 mutation or the 1p19q co-deletion have a better prognosis than those with the wild-type (32-34). The MGMT methylation status is an additional molecular biomarker that has a significant effect on predicting patient survival and treatment response (35). However, some alterations would make GBM more aggressive and result in a worse prognosis (36). Even though the fact that present prognostic and predictive indicators have a role in directing the clinical care of patients, there is still a need to investigate new effective prognostic biomarkers.

To test the independence and clinical utility of the 14-gene prognostic model, univariate and multivariate Cox regression analyses were performed on clinical factors, as shown in Table 1. These analyses indicated that the 14-gene marker model was substantially linked with survival rate. It was hypothesized that our 14-gene prognostic model may predict the survival and prognosis of GBM patients independently. The predictive performance of the 14-gene model was superior to those of prognostic parameters such as age [univariable odds ratio (OR) =2.19, 95% confidence interval (CI): 1.51–3.19, P<0.001] and IDH mutation (univariable OR =0.24, 95% CI: 0.17–0.34, P<0.001). Primary versus Recurrence, World Health Organization (WHO) stage II against III versus IV, MGMT promoter methylation yes versus no, and 1p19q deletion yes versus no (Figure 5A). However, there was no significant difference in overall survival between males and females, suggesting that gender is not a prognostic factor (Figure S1).

Figure 5 HRD can serve as a biomarker for PARP inhibitors in GBM. (A) Forest plot of clinical variables and the prognostic of GBM. (B) Loss of fitness score histograms of the two groups when PARP1 or PARP2 was eliminated. (C) Histograms of IC₅₀ values for olaparib and niraparib in the two groups. OR, odds ratio; CI, confidence interval; IDH, isocitrate dehydrogenase; Rx, primary or recurrent tumor; MGMT, O6-methylguanine-DNA methyltransferase; IC₅₀, half maximal inhibitory concentration; HRD, homologous recombination deficiency; GBM, glioblastoma.

Predicting the therapy with PARP inhibitors

Given that HRD status might determine whether individuals with breast and ovarian cancer would benefit from PARP inhibitor treatment, we investigated the responsiveness of GBM tumor cell lines to PARP inhibitors using in vitro experiment data (37-40). First, knocking out the PARP1 gene resulted in a considerably greater loss of fitness score in the low-risk group than in the high-risk group (Wilcoxon signed-rank test, P=0.02). When the PARP2 gene was silenced, the two groups exhibited a similar pattern, but without statistical significance (Figure 5B).

Additionally, to determine the independent cytotoxic impact of PARP inhibitors in GBM, the half maximal inhibitory concentration (IC₅₀) values of olaparib and niraparib in GBM cell lines were analyzed. Both olaparib and niraparib showed considerably lower IC₅₀s in the group with a low-risk score compared to the group with a high-risk score (Figure 5C). These findings imply a possible therapeutic advantage of PARP inhibitors for GBM patients with a low-risk score.

Discussion

HRD is a common genetic alteration that is commonly seen in malignancies. HRD research contributes to the advancement of our understanding of cancers. Genomic instability is linked to the incidence of GBM, although no research has elucidated the impact and importance of HRD in GBM. GBM is a tumor without BRCA1/2 mutations but with additional genomic instability characteristics. In addition, no gene panel gold standard for HRD detection has been provided in GBM. The HRD score is a mutational signature-based method for predicting HRD. It identifies the effects of HRD rather than its fundamental cause and can predict HRD in GBM more accurately. Our study is the first to investigate the predictive significance of HRD in GBM.

In our discovery group, 11.5% of GBM patients were HRD-positive (HRD score 17), and these patients exhibited considerably enriched DDR signaling pathways and a higher likelihood of survival. TMZ is the most often utilized chemotherapeutic agent for GBM patients. TMZ is an oral alkylation drug that covalently links a methyl group to a guanine base in tumor DNA, resulting in DNA mismatches during cell replication. The mismatch repair system might result in the death and regression of tumor cells (41). HRD can possibly generate a high level of genomic instability and boost the activity of TMZ. On the other hand, this also suggests a better outcome for the HRD group. Previous research has demonstrated that platinum-based chemotherapy (another alkylation agent) or PARP inhibitors, which break DNA, can exert maximal synthetic lethality in HRD patients, therefore destroying breast cancer tumors (42). In the treatment of GBM patients, platinum-based medicines have a tumultuous history, with early benefits but dose-limiting toxicity. Currently, platinum medication therapy for GBM is showing signs of development and effectiveness (43). The relationship between HRD-positive GBM patients, particularly those without MGMT methylation, and the efficacy of platinum treatment require additional investigation.

IDH1, TP53, and ATRX were shown to have considerably greater mutation frequencies in HRD subgroups as compared to HRP subgroups when comparing gene mutation profiles. IDH1 is involved in DNA demethylation, histone demethylation, and an additional essential epigenetic regulatory function. Over 80% of low-grade gliomas (LGGs) and subsequent GBMs include mutations in IDH1/2. Mutations in IDH1 lead to dysregulation of histone and DNA methylation, as well as the potential activation of endogenous retroviruses (ERVs) (44), which may result from IDH1 mutations, resulting in genomic instability. ATRX is a chromatin regulator that operates as a member of Switch/sucrose non-fermentable (SWI/SNF) helicase family (45). ATRX is necessary for homologous recombination, facilitates activities such as DNA repair synthesis and sister chromatid exchange, and its deficiency promotes genomic instability (46). The TP53 gain-of-function mutation promotes inflammation in GBM, as evidenced by the upregulation of C-C motif chemokine ligand 2 (CCL2) and tumor necrosis factor alpha (TNFA) via nuclear factor kappa B (NFB) signaling, resulting in an increase in microglia and monocyte-derived immune cell infiltration (47). Although TP53 mutations have been reported as a kind of poor prognostic factor for GBM (48), TP53 mutations coincide closely with ATRX and IDH1 mutations (49), therefore it is logical that TP53 mutations are highly concentrated in the HRD subgroup with a better prognosis. It is consistent with the finding that PTEN and EGFR mutations are indicators of a poor prognosis since the mutation frequency of PTEN and EGFR genes were significantly greater in the HRP subgroup with a bad prognosis. The deletions or mutations of PTEN are associated with treatment resistance in GBMs with a poor prognosis (50), and EGFR mutations are more aggressive (51).

We performed a LASSO Cox regression analysis of GBM cases according to the 14-gene prognostic model constructed by HRD-related DEGs. These 14 genes were CHL1, CCL2, PDPN, POSTN, CHI3L1, IL8, FABP5, IGFBP3, LGALS3, EMP3, TAGLN, DCX, GPR37, C21orf62; and the expression levels of all the genes were significantly decreased in the low-risk-score subgroup, except for DCX. Transcriptome analysis of long-term surviving GBM patients showed that reduced messenger RNA (mRNA) expression levels of EMP3, LGALS3, and IGFBP3 genes might predict poorer survival outcomes (52). High expression of IGFBP3 messenger RNA (mRNA) is positively correlated with tumor grade and predicted shorter overall survival. Decreased or lost expression of IGFBP3 can inhibit cell proliferation and induce G2/M arrest of cell cycle and apoptosis of glioma cells, as well as inducing accumulation of DNA damage (53). A study showed that the high expression of LGALS3 can promote the treatment resistance of GBM, and is related to tumor risk and prognosis (54). LGALS3 may play a role through cell death inhibition pathways that may involve the Bcl-2 family and promote chemoradiotherapy resistance by preferentially activating DNA damage checkpoint responses (55). EMP3 not only plays a role in promoting tumor invasion mediated by the PI3K/Akt pathway (56) but also plays an immunosuppressive role with the high expression of EMP3 accompanied by the recruitment of a large number of M2 macrophages and the reduction of T cell infiltration (57). The high expression and nuclear accumulation of DCX have improved the ability of invasive gliomas in vitro and in vivo (58), while, the specific mechanism by how DCX regulates tumor cell survival in GBM has not been investigated. So far, the study on the role of CHL1 in GBM has suggested that high expression of CHL1 could promote the occurrence of malignant tumors (59). CCL2 can stimulate GBM cancer cell migration and invasion by increasing angiogenesis (60). Furthermore, CCL2-CCR2 signaling activated by CCL2 overexpression would lead to evasion of immune-mediated tumor cell killing (61). PDPN is overexpressed in several solid tumors associated with tumor malignant progression, epithelial-mesenchymal transition, and metastasis. PDPN is primarily expressed in the mesenchymal type of GBM, and positively correlated with tumor malignancy (62). Silencing POSTN can visibly reduce tumor-associated macrophages, thereby inhibiting tumor growth in xenografts (63). The overexpression of IL8 can induce endothelial cell proliferation and tubular cell formation, promoting angiogenesis progression in GBM (64). The functions of these genes in GBM that we summarized also proved the rationality of the genes selected to construct the prognostic model indirectly. The relationship between CHI3L1, FABP5, C21orf62, TAGLN, and GPR37 genes and GBM has not been described in the literature.

Through our 14-gene prognostic model, both the TCGA cohort and the CGGA cohort can be divided into two subgroups in GBM, the low-risk-score group with better survival outcomes. It was possible because the low-risk-score group had a higher HRD score, and was more associated with many oncological features, such as significantly enriched DDR pathway and γδT cell signal. γδT cells are immune cells that both kill tumor cells and tumor stem cells and recognize antigens (65). γδT cells can directly kill tumor cells by secreting cytokines (IFN-γ, TNF-α), NK cell receptors on the cell surface, and antibody-dependent cytotoxicity (ADCC) effect (66). Alternatively, B/DC/αβT/NK cells can be activated in various ways, such as acting as antigen-presenting cells to activate αβT cells, or inducing NK-mediated antitumor cytotoxicity by the 4-1BB costimulation pathway, to achieve indirect tumor killing (67). Enrichment of γδT cells is associated with a good prognosis in some cancers. In lung adenocarcinoma (68), particularly in long-surviving patients, elevated γδT cells in tumors are mainly IL-17A-releasing γδT17 cells. In tumor xenografts of mice, γδT cells completely eradicated preexisting tumors with an initial size of approximately 5 mm. It may be partly due to the interaction of NKG2D receptors on γδ-T cells with their tumor-expressed ligands, thereby overcoming the inhibitory signals of MHC I molecules.

GBM is a highly invasive illness with a poor prognosis and high tumor heterogeneity, hence it is very important to undertake precise prognostic stratification on patients. This is the first study to investigate the prognostic significance of HRD and HRD-related gene expression in GBM. According to our research, HRD-related traits may function as a novel biomarker. In addition, a 14-gene prognosis classification model based on HRD-related DEGs might be employed as a prognostic factor in its own right. Innovatively, we have developed a GBM prognostic model based on genomic instability that is highly stable, has broad application, and makes it easier to identify patients with favorable prognoses. It outlines potential prerequisites for translating this model into clinically applicable applications.

It is believed that the outcomes of this study may give crucial insights and lead to a more effective therapeutic therapy and prognosis for GBM. As our study was done in a cohort of a public database, we did not gather thorough clinical data regarding the clinical care of patients and therapy alternatives. The prognostic mechanism of genomic instability in GBM requires more research. In this review, the analysis and discussion of the HRD-related genome, transcriptome, and TIME would give a new technique for identifying therapy and prognosis-related biomarkers.

Conclusions

Genomic instability is a classical characteristic of tumors, and the HRD score is an important method to evaluate genomic instability. We found for the first time that the HRD score has an important prognostic role in GBM when the HRD score was greater than 17, the survival results of patients were significantly better than those with an HRD score less than 17. Our research provides a new perspective for understanding GBM. Next, by the HRD-related transcriptome characteristics analysis, the prognosis model based on DEGs constructed by the machine learning method can predict patient survival more accurately than the HRD score. Finally, the exploration results of the TIME suggest that the enrichment of γδ T cell characteristic signals may be the specific mechanism of different prognoses in GBM.

Acknowledgments

Funding: This work was supported by grants from the Project of Guangzhou Science and Innovation Commission (No. 201707010380); Guangdong Province Medical Science and Technology Research Fund (No. A2024087); Guangdong Province Medical Science and Technology Research Fund (No. B202411); Guangzhou Health Science and Technology Project Fund (No. 20241A011093).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2077/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2077/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2077/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Ylanan AMD, Pascual JSG, Cruz-Lim EMD, et al. Intraoperative radiotherapy for glioblastoma: A systematic review of techniques and outcomes. J Clin Neurosci 2021;93:36-41. [Crossref] [PubMed]
2. Wen PY, Weller M, Lee EQ, et al. Glioblastoma in adults: a Society for Neuro-Oncology (SNO) and European Society of Neuro-Oncology (EANO) consensus review on current management and future directions. Neuro Oncol 2020;22:1073-113. [Crossref] [PubMed]
3. Wang Y, Xu H, Liu T, et al. Temporal DNA-PK activation drives genomic instability and therapy resistance in glioma stem cells. JCI Insight 2018;3:e98096. [Crossref] [PubMed]
4. Wu S, Li X, Gao F, et al. PARP-mediated PARylation of MGMT is critical to promote repair of temozolomide-induced O6-methylguanine DNA damage in glioblastoma. Neuro Oncol 2021;23:920-31. [Crossref] [PubMed]
5. Ibrahim Al-Obaide MA, Arutla V, Bacolod MD, et al. Genomic Space of MGMT in Human Glioma Revisited: Novel Motifs, Regulatory RNAs, NRF1, 2, and CTCF Involvement in Gene Expression. Int J Mol Sci 2021;22:2492. [Crossref] [PubMed]
6. Rao AM, Quddusi A, Shamim MS. The significance of MGMT methylation in Glioblastoma Multiforme prognosis. J Pak Med Assoc 2018;68:1137-9. [PubMed]
7. Omabe K, Uduituma S, Igwe D, et al. Deeper Insight in Metastatic Cancer Progression; Epithelial-to- Mesenchymal Transition and Genomic Instability: Implications on Treatment Resistance. Curr Mol Med 2021;21:860-71. [Crossref] [PubMed]
8. Vanacker H, Harter P, Labidi-Galy SI, et al. PARP-inhibitors in epithelial ovarian cancer: Actual positioning and future expectations. Cancer Treat Rev 2021;99:102255. [Crossref] [PubMed]
9. Nguyen L, W M, Martens J, Van Hoeck A, et al. Pan-cancer landscape of homologous recombination deficiency. Nat Commun 2020;11:5584. [Crossref] [PubMed]
10. Knijnenburg TA, Wang L, Zimmermann MT, et al. Genomic and Molecular Landscape of DNA Damage Repair Deficiency across The Cancer Genome Atlas. Cell Rep 2018;23:239-254.e6. [Crossref] [PubMed]
11. Stover EH, Fuh K, Konstantinopoulos PA, et al. Clinical assays for assessment of homologous recombination DNA repair deficiency. Gynecol Oncol 2020;159:887-98. [Crossref] [PubMed]
12. Sztupinszki Z, Diossy M, Krzystanek M, et al. Migrating the SNP array-based homologous recombination deficiency measures to next generation sequencing data of breast cancer. NPJ Breast Cancer 2018;4:16. [Crossref] [PubMed]
13. Takaya H, Nakai H, Takamatsu S, et al. Homologous recombination deficiency status-based classification of high-grade serous ovarian carcinoma. Sci Rep 2020;10:2757. [Crossref] [PubMed]
14. Rostami M, Forouzandeh S, Berahmand K, et al. Gene selection for microarray data classification via multi-objective graph theoretic-based method. Artif Intell Med 2022;123:102228. [Crossref] [PubMed]
15. Marquard AM, Eklund AC, Joshi T, et al. Pan-cancer analysis of genomic scar signatures associated with homologous recombination deficiency suggests novel indications for existing cancer drugs. Biomark Res 2015;3:9. [Crossref] [PubMed]
16. Zhao Z, Zhang KN, Wang Q, et al. Chinese Glioma Genome Atlas (CGGA): A Comprehensive Resource with Functional Genomic Data from Chinese Glioma Patients. Genomics Proteomics Bioinformatics 2021;19:1-12. [Crossref] [PubMed]
17. Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn) 2015;19:A68-A77. [Crossref] [PubMed]
18. Abkevich V, Timms KM, Hennessy BT, et al. Patterns of genomic loss of heterozygosity predict homologous recombination repair defects in epithelial ovarian cancer. Br J Cancer 2012;107:1776-82. [Crossref] [PubMed]
19. Manié E, Popova T, Battistella A, et al. Genomic hallmarks of homologous recombination deficiency in invasive breast carcinomas. Int J Cancer 2016;138:891-900. [Crossref] [PubMed]
20. Birkbak NJ, Wang ZC, Kim JY, et al. Telomeric allelic imbalance indicates defective DNA repair and sensitivity to DNA-damaging agents. Cancer Discov 2012;2:366-75. [Crossref] [PubMed]
21. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
22. Kanehisa M, Furumichi M, Tanabe M, et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017;45:D353-61. [Crossref] [PubMed]
23. Hastie T, Qian J, Tay K. An Introduction to glmnet. CRAN R Repositary 2021. Available online: https://cran.r-project.org/web/packages/glmnet/vignettes/glmnet.pdf
24. Subramanian A, Kuehn H, Gould J, et al. GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics 2007;23:3251-3. [Crossref] [PubMed]
25. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
26. BenjaminDSatoTCibulskisKCalling somatic SNVs and indels with Mutect2.BioRxiv 2019:861054. 10.1101/861054
27. Shi Z, Qu X, Guo C, et al. Identification of clinical trait-related small RNA biomarkers with weighted gene co-expression network analysis for personalized medicine in endocervical adenocarcinoma. Aging (Albany NY) 2021;13:22361-74. [Crossref] [PubMed]
28. Yang C, Zhang Z, Tang X, et al. Pan-cancer analysis reveals homologous recombination deficiency score as a predictive marker for immunotherapy responders. Hum Cell 2022;35:199-213. [Crossref] [PubMed]
29. Li W, Lu L, Lu J, et al. cGAS-STING-mediated DNA sensing maintains CD8(+) T cell stemness and promotes antitumor T cell therapy. Sci Transl Med 2020;12:eaay9013. [Crossref] [PubMed]
30. Nakazawa T, Nakamura M, Park YS, et al. Cytotoxic human peripheral blood-derived γδT cells kill glioblastoma cell lines: implications for cell-based immunotherapy for patients with glioblastoma. J Neurooncol 2014;116:31-9. [Crossref] [PubMed]
31. Zheng S, Alfaro-Munoz K, Wei W, et al. Prospective Clinical Sequencing of Adult Glioma. Mol Cancer Ther 2019;18:991-1000. [Crossref] [PubMed]
32. Bhavya B, Anand CR, Madhusoodanan UK, et al. To be Wild or Mutant: Role of Isocitrate Dehydrogenase 1 (IDH1) and 2-Hydroxy Glutarate (2-HG) in Gliomagenesis and Treatment Outcome in Glioma. Cell Mol Neurobiol 2020;40:53-63. [Crossref] [PubMed]
33. Wang LB, Karpova A, Gritsenko MA, et al. Proteogenomic and metabolomic characterization of human glioblastoma. Cancer Cell 2021;39:509-528.e20. [Crossref] [PubMed]
34. Louis DN, Perry A, Reifenberger G, et al. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol 2016;131:803-20. [Crossref] [PubMed]
35. Mansouri A, Hachem LD, Mansouri S, et al. MGMT promoter methylation status testing to guide therapy for glioblastoma: refining the approach based on emerging evidence and current challenges. Neuro Oncol 2019;21:167-78. [Crossref] [PubMed]
36. Fang R, Chen X, Zhang S, et al. EGFR/SRC/ERK-stabilized YTHDF2 promotes cholesterol dysregulation and invasive growth of glioblastoma. Nat Commun 2021;12:177. [Crossref] [PubMed]
37. Dréan A, Lord CJ, Ashworth A. PARP inhibitor combination therapy. Crit Rev Oncol Hematol 2016;108:73-85. [Crossref] [PubMed]
38. Lord CJ, Ashworth A. Targeted therapy for cancer using PARP inhibitors. Curr Opin Pharmacol 2008;8:363-9. [Crossref] [PubMed]
39. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
40. Dwane L, Behan FM, Gonçalves E, et al. Project Score database: a resource for investigating cancer cell dependencies and prioritizing therapeutic targets. Nucleic Acids Res 2021;49:D1365-72. [Crossref] [PubMed]
41. Lee SY. Temozolomide resistance in glioblastoma multiforme. Genes Dis 2016;3:198-210. [Crossref] [PubMed]
42. Telli ML, Timms KM, Reid J, et al. Homologous Recombination Deficiency (HRD) Score Predicts Response to Platinum-Containing Neoadjuvant Chemotherapy in Patients with Triple-Negative Breast Cancer. Clin Cancer Res 2016;22:3764-73. [Crossref] [PubMed]
43. Robert M, Patsouris A, Frenel JS, et al. Emerging PARP inhibitors for treating breast cancer. Expert Opin Emerg Drugs 2018;23:211-21. [Crossref] [PubMed]
44. Turcan S, Makarov V, Taranda J, et al. Mutant-IDH1-dependent chromatin state reprogramming, reversibility, and persistence. Nat Genet 2018;50:62-72. [Crossref] [PubMed]
45. Argentaro A, Yang JC, Chapman L, et al. Structural consequences of disease-causing mutations in the ATRX-DNMT3-DNMT3L (ADD) domain of the chromatin-associated protein ATRX. Proc Natl Acad Sci U S A 2007;104:11939-44. [Crossref] [PubMed]
46. Juhász S, Elbakry A, Mathes A, et al. ATRX Promotes DNA Repair Synthesis and Sister Chromatid Exchange during Homologous Recombination. Mol Cell 2018;71:11-24.e7. [Crossref] [PubMed]
47. Ham SW, Jeon HY, Jin X, et al. TP53 gain-of-function mutation promotes inflammation in glioblastoma. Cell Death Differ 2019;26:409-25. [Crossref] [PubMed]
48. Zhang Y, Dube C, Gibert M Jr, et al. The p53 Pathway in Glioblastoma. Cancers (Basel) 2018;10:297. [Crossref] [PubMed]
49. Liu XY, Gerges N, Korshunov A, et al. Frequent ATRX mutations and loss of expression in adult diffuse astrocytic tumors carrying IDH1/IDH2 and TP53 mutations. Acta Neuropathol 2012;124:615-25. [Crossref] [PubMed]
50. Zhang P, Meng X, Liu L, et al. Identification of the Prognostic Signatures of Glioma With Different PTEN Status. Front Oncol 2021;11:633357. [Crossref] [PubMed]
51. Eskilsson E, Røsland GV, Solecki G, et al. EGFR heterogeneity and implications for therapeutic intervention in glioblastoma. Neuro Oncol 2018;20:743-52. [Crossref] [PubMed]
52. Gerber NK, Goenka A, Turcan S, et al. Transcriptional diversity of long-term glioblastoma survivors. Neuro Oncol 2014;16:1186-95. [Crossref] [PubMed]
53. Chen CH, Chen PY, Lin YY, et al. Suppression of tumor growth via IGFBP3 depletion as a potential treatment in glioma. J Neurosurg 2019;132:168-79. [Crossref] [PubMed]
54. Hu WM, Yang YZ, Zhang TZ, et al. LGALS3 Is a Poor Prognostic Factor in Diffusely Infiltrating Gliomas and Is Closely Correlated With CD163+ Tumor-Associated Macrophages. Front Med (Lausanne) 2020;7:182. [Crossref] [PubMed]
55. Wang H, Song X, Huang Q, et al. LGALS3 Promotes Treatment Resistance in Glioblastoma and Is Associated with Tumor Risk and Prognosis. Cancer Epidemiol Biomarkers Prev 2019;28:760-9. [Crossref] [PubMed]
56. Martija AA, Pusch S. The Multifunctional Role of EMP3 in the Regulation of Membrane Receptors Associated with IDH-Wild-Type Glioblastoma. Int J Mol Sci 2021;22:5261. [Crossref] [PubMed]
57. Chen Q, Jin J, Huang X, et al. EMP3 mediates glioblastoma-associated macrophage infiltration to drive T cell exclusion. J Exp Clin Cancer Res 2021;40:160. [Crossref] [PubMed]
58. Ayanlaja AA, Ji G, Wang J, et al. Doublecortin undergo nucleocytoplasmic transport via the RanGTPase signaling to promote glioma progression. Cell Commun Signal 2020;18:24. [Crossref] [PubMed]
59. Yang Z, Xie Q, Hu CL, et al. CHL1 Is Expressed and Functions as a Malignancy Promoter in Glioma Cells. Front Mol Neurosci 2017;10:324. [Crossref] [PubMed]
60. Lim SY, Yuzhalin AE, Gordon-Weeks AN, et al. Targeting the CCL2-CCR2 signaling axis in cancer metastasis. Oncotarget 2016;7:28697-710. [Crossref] [PubMed]
61. Chang AL, Miska J, Wainwright DA, et al. CCL2 Produced by the Glioma Microenvironment Is Essential for the Recruitment of Regulatory T Cells and Myeloid-Derived Suppressor Cells. Cancer Res 2016;76:5671-82. [Crossref] [PubMed]
62. Shiina S, Ohno M, Ohka F, et al. CAR T Cells Targeting Podoplanin Reduce Orthotopic Glioblastomas in Mouse Brains. Cancer Immunol Res 2016;4:259-68. [Crossref] [PubMed]
63. Zhou W, Ke SQ, Huang Z, et al. Periostin secreted by glioblastoma stem cells recruits M2 tumour-associated macrophages and promotes malignant growth. Nat Cell Biol 2015;17:170-82. [Crossref] [PubMed]
64. Conroy S, Kruyt FAE, Wagemakers M, et al. IL-8 associates with a pro-angiogenic and mesenchymal subtype in glioblastoma. Oncotarget 2018;9:15721-31. [Crossref] [PubMed]
65. Capsomidis A, Benthall G, Van Acker HH, et al. Chimeric Antigen Receptor-Engineered Human Gamma Delta T Cells: Enhanced Cytotoxicity with Retention of Cross Presentation. Mol Ther 2018;26:354-65. [Crossref] [PubMed]
66. Kanazawa T, Hiramatsu Y, Iwata S, et al. Anti-CCR4 monoclonal antibody mogamulizumab for the treatment of EBV-associated T- and NK-cell lymphoproliferative diseases. Clin Cancer Res 2014;20:5075-84. [Crossref] [PubMed]
67. Maniar A, Zhang X, Lin W, et al. Human gammadelta T lymphocytes induce robust NK cell-mediated antitumor cytotoxicity through CD137 engagement. Blood 2010;116:1726-33. [Crossref] [PubMed]
68. Cairo C, Surendran N, Harris KM, et al. Vγ2Vδ2 T cell Costimulation Increases NK cell Killing of Monocyte-derived Dendritic Cells. Immunology 2014; Epub ahead of print. [Crossref] [PubMed]