Hub genes and their key effects on prognosis of Burkitt lymphoma

Abstract

BACKGROUND

Burkitt lymphoma (BL) is an exceptionally aggressive malignant neoplasm that arises from either the germinal center or post-germinal center B cells. Patients with BL often present with rapid tumor growth and require high-intensity multi-drug therapy combined with adequate intrathecal chemotherapy prophylaxis, however, a standard treatment program for BL has not yet been established. It is important to identify biomarkers for predicting the prognosis of BLs and discriminating patients who might benefit from the therapy. Microarray data and sequencing information from public databases could offer opportunities for the discovery of new diagnostic or therapeutic targets.

AIM

To identify hub genes and perform gene ontology (GO) and survival analysis in BL.

METHODS

Gene expression profiles and clinical traits of BL patients were collected from the Gene Expression Omnibus database. Weighted gene co-expression network analysis (WGCNA) was applied to construct gene co-expression modules, and the cytoHubba tool was used to find the hub genes. Then, the hub genes were analyzed using GO and Kyoto Encyclopedia of Genes and Genomes analysis. Additionally, a Protein-Protein Interaction network and a Genetic Interaction network were constructed. Prognostic candidate genes were identified through overall survival analysis. Finally, a nomogram was established to assess the predictive value of hub genes, and drug-gene interactions were also constructed.

RESULTS

In this study, we obtained 8 modules through WGCNA analysis, and there was a significant correlation between the yellow module and age. Then we identified 10 hub genes (SRC, TLR4, CD40, STAT3, SELL, CXCL10, IL2RA, IL10RA, CCR7 and FCGR2B) by cytoHubba tool. Within these hubs, two genes were found to be associated with OS (CXCL10, P = 0.029 and IL2RA, P = 0.0066) by survival analysis. Additionally, we combined these two hub genes and age to build a nomogram. Moreover, the drugs related to IL2RA and CXCL10 might have a potential therapeutic role in relapsed and refractory BL.

CONCLUSION

From WGCNA and survival analysis, we identified CXCL10 and IL2RA that might be prognostic markers for BL.

Key Words: Burkitt lymphoma; Weighted gene co-expression network analysis; Microarray data; Functional enrichment analysis; Prognosis; Therapeutic target

Core Tip: This study represents the pioneering investigation of gene expression in Burkitt lymphoma (BL) using weighted gene co-expression network analysis, coupled with functional enrichment analysis. In this study, we have successfully identified and validated 10 hub genes. Survival analysis has demonstrated that the overexpression of CXCL10 and IL2RA in BL may serve as robust prognostic indicators. Furthermore, an integrated mRNA signature and age nomogram potentially provide valuable prognostic insights for patients with BLs.

INTRODUCTION

Burkitt lymphoma (BL) is a highly aggressive B-cell non-Hodgkin's lymphoma characterized by the t (8; 14) chromosomal translocation involving the MYC oncogene and the immunoglobulin heavy chain gene (IGH)[1]. Three distinct clinical subtypes of BL have been identified: Namely endemic (African), sporadic (non-endemic), and immunodeficiency-associated. Notably, chronic Epstein-Barr virus infection plays a pivotal role in the pathogenesis of BL, particularly in the endemic subtype[2]. Endemic BL is primarily found in countries located near the equator in Africa. The estimated annual incidence of endemic BL is 3-6 cases per 100000 children in African countries[3], which is approximately 50 times higher than that in the United States[4]. Sporadic BL predominantly occurs in the United States and Western Europe. The annual incidence of BL in the United States is approximately 3 cases per 1 million individuals, while in Europe it stands at around 2.2 cases per 1 million people[5]. Immunodeficiency-associated BL primarily affects individuals with HIV infection, typically those with relatively high CD4 counts and no opportunistic infections[6].

Patients with BL frequently exhibit rapid tumor growth, spontaneous tumor lysis, and elevated levels of serum lactate dehydrogenase. Currently, patients with BL necessitate a high-intensity multi-drug regimen in conjunction with adequate intrathecal central nervous system prophylaxis. However, the absence of an established standard treatment protocol for BL persists[7]. BL is an aggressive lymphoma, which can potentially be cured; however, patients with refractory and relapsed disease have an extremely poor prognosis[8]. Therefore, it is important to identify robust biomarkers for predicting the prognosis of BLs and discriminating patients who might benefit from therapy. The development of BL depends on the constitutive expression of the MYC gene located on chromosome 8q24, which encodes the transcription factor protein MYC[9]. MYC orchestrates the expression of target genes, regulating a variety of cellular processes, including cell growth, division, apoptosis, metabolism, adhesion, and motility[10]. MYC gene rearrangements are seen in the vast majority of BLs, and factors other than MYC translocation need to be present in the process of BL. However, it is not clear why and how B cells develop genetic alterations that result in increased MYC expression and ultimately lead to BL.

The Gene Expression Omnibus (GEO) is an international public repository constructed and maintained by the National Center for Biotechnology Information[11]. At the time of writing, the GEO database hosts more than 194000 public series. Weighted gene co-expression network analysis (WGCNA) is a widely used systematic biological method for generating gene co-expression networks[12,13]. In this study, WGCNA was first used to analyze genes of BL samples mined from the GEO database. Subsequently, we identified these hub genes and conducted a functional enrichment analysis. Additionally, a survival analysis was conducted to identify an mRNA signature that exhibits a significant association with prognosis. Finally, a prognostic nomogram was established based on the combination of gene signature and clinical characteristics.

MATERIALS AND METHODS

Data collection and preprocessing

The raw gene expressions and the corresponding clinical follow-up data of GSE4475 and GSE69051 were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/. Accessed Jan 20, 2023)[14], and the two datasets were built based on the GPL96 platform (HG-U133A) and GPL14951 platform (Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip) respectively. Analysis was performed on the raw gene expression data of the BL datasets and the corresponding clinical follow-up obtained from GSE4475, which included a total of 36 BL samples. The survival data of the hub genes was verified by downloading another dataset, GSE69051, which included 77 BL samples. The mRNA sequencing data annotation information was used to match the probe with the corresponding gene and transform the gene name into a gene symbol. Probes that corresponded to more than one gene were excluded from the dataset.

Co-expression network construction

WGCNA converts gene expression data into co-expression modules, establishing relationships between genes and focusing on gene modules rather than individual genes[15]. Besides, WGCNA can identify the gene modules related to clinical traits and has been widely used in cancer research. In this study, the top 5000 most variable genes were used to construct a co-expression network by using the package of WGCNA in R[13]. The power value was filtered out during the module construction process using the WGCNA algorithm. The mean connectivity and scale independence of network modules were analyzed using the gradient test under different power values, which ranged from 1 to 20. When the degree of independence was 0.85, the appropriate power value was determined. Then, the soft threshold test was performed. In this study, the soft threshold β was 12, and the network type was “signed”. The WGCNA algorithm further identified co-expression modules under these conditions. The minimum size of the gene group was set at 100 to ensure the reliability of the results for this module. Then, the correlation between the characteristics of the module-trait association module and clinical traits was visually expressed. The relationship between the expression profile and traits was analyzed to make a scatter plot between gene significance and module membership.

Hub genes identification and functional analysis

The Protein-Protein Interaction (PPI) network of the interested module was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (available from https://cn.string-db.org/. Accessed 25 Jan 2023)[16] and presented by Cytoscape software. The cytoHubba tool was used to screen the hub genes. Then, the selected hub genes were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using Database for Annotation, Visualization, and Integrated Discovery (DAVID v.6.8: available from https://david.ncifcrf.gov/. Accessed Jan 29, 2023)[17]. The possible functions were analyzed by biological processes, cellular components, and molecular functions, while the potential signal pathways were analyzed using KEGG.

Construction of hub genes PPI and genetic interaction network

The PPI network was used to analyze the hub genes at the protein level, and the STRING database (available from https://string-db.org/. Accessed Jan 31, 2023) was used to check and predict the interaction between proteins[16]. The genetic interaction (GI) network, constructed using gene function prediction, aims to understand the complex interactions between genes. We used the Gene Multiple Association Network Integration Algorithm (GeneMANIA, available from https://genemania.org/. Accessed Jan 31, 2023) to analyze the hub genes[18]. The threshold of a collective score of 0.15 was implemented.

Statistical analysis

Based on the 50^th percentile cut-off value of each hub gene mRNA, patients were divided into high-expression and low-expression groups. Log-rank test and Kaplan-Meier estimation were performed to obtain log-rank P value and evaluate hub genes in overall survival (OS). Cox regression analysis was conducted to examine the association between the risk score and clinical information, and a nomogram was generated. The survival curve and nomogram were carried out by R version 4.2.1. Additionally, P < 0.05 was statistically significant.

Drug-gene interaction

The DGIdb database (available from https://dgidb.org/. Accessed Feb 20, 2023) was used to investigate drug-gene interaction to identify drugs associated with hub genes[19]. The interaction network was visualized via Cytoscape.

RESULTS

Construction and screening of BL co-expression module

In this study, we obtained the BL dataset from GSE4475, resulting in a total of 13514 gene expression values. The clinical features of the BL samples are listed in Table 1. Then, we selected a total of 5000 genes with the highest average expression values for cluster analysis. Firstly, the clustering tree of 36 samples of BL was extracted from GSE4475 (Figure 1). Secondly, we calculated the soft threshold (power value), and when the weight was equal to 12, the independence exceeded 0.85, indicating higher average connectivity (Supplementary Figure 1). By utilizing this power value for hierarchical clustering analysis and combining similar analysis results, a total of 8 modules were identified, including black (1073 genes), blue (967 genes), brown (853 genes), green yellow (140 genes), grey (1019 genes), magenta (219 genes), pink (267 genes) and yellow (462 genes) (Figure 2A). Genes in grey were not included in any module, thus we analyzed the interactive relationships underlying the 7 co-expression modules (Figure 2B). Given the well-established association between age and prognosis in BL patients, we opted to investigate the module that exhibited the strongest correlation with age for subsequent analysis[20,21]. A significant correlation between the yellow module and age was discovered (Figure 3). The correlation between modules and samples is shown in Supplementary Figure 2. Finally, we conducted a scatter diagram of the correlation between the yellow module and age (Figure 4).

Figure 1 Clustering tree of 36 samples of Burkitt lymphoma extracted from GSE4475. Red indicated more gene expression, white less, and grey indicated deletion. CCS: Chromosomal Complexity Score.

Figure 2 Sample clustering to detect outliers and construction of co-expression modules. A: The constructed co-expression modules of Burkitt lymphoma genes by weighted gene co-expression network analysis; B: Interaction analysis between gene co-expression modules. The heatmap showed the Topological Overlap Matrix among genes in the analysis. Different colors on the x-axis and y-axis represented different modules. The intensity of inter-module connections was visually represented by the yellow brightness in the central region, gradually transitioning into deeper shades of orange.

Figure 3 Module-trait association. Correlation thermography between modular feature genes and clinical features of Burkitt lymphoma. Each row corresponded to a module feature, and the column corresponded to a clinical feature. Each cell contained the correlation and the corresponding P value. CCS: Chromosomal Complexity Score; ME: Module membership.

Figure 4 The scatter plot of the correlation for an age-related gene between module membership and gene significance in the yellow module.

Table 1 Clinical features of Burkitt lymphoma patients.

Clinical features	Total (n = 36)
Age, mean (range)	31.0 (2-90)
Gender
Male	24
Female	11
Unknown	1
Stage
I	4
II	10
III	5
IV	6
Unknown	11
Survival status
Alive	20
Dead	7
Unknown	9
Ki 67
≤ 75%	4
75%-90%	9
> 90%	22
Unknown	1
CCS
< 10	29
≥ 10	5
Unknown	4

CCS: Chromosomal Complexity Score.

Hub genes identification

All of the genes from the yellow module were uploaded to the STRING database, and a PPI network was constructed using Cytoscape software (Supplementary Figure 3). And the top 10 hub genes (SRC, TLR4, CD40, STAT3, SELL, CXCL10, IL2RA, IL10RA, CCR7 and FCGR2B) were screened out by cytoHubba tool. GeneMANIA showed the GI network of hub genes interaction at the mRNA expression level (Figure 5A). The STRING database generated the PPI co-expression network by analyzing the hub genes at the protein level (Figure 5B).

Figure 5 Genetic and Protein-Protein interaction network of hub genes. A: GeneMANIA was used to construct a genetic interaction network. The black nodes with a slash represent the query gene, while the other nodes represent the predicted genes. The purple edges indicate co-expression, whereas the blue edges signify co-localization; B: A physically and functionally connected Protein-Protein Interaction network implemented common goals through Search Tool for the Retrieval of Interacting Genes/Proteins, where nodes represented proteins and edges represented pairs of interactions between proteins. Node size and color indicated richness, while edge size and color reflected combined scores.

Functional and pathway enrichment analysis

Enrichment analyses of GO and KEGG were conducted to explore potential pathways of the hub genes. Forty-five GO-enriched terms were shown in Supplementary Table 1. The top 10 GO terms (Figure 6A) included inflammatory response, external side of plasma membrane, plasma membrane, cellular response to lipopolysaccharide, positive regulation of interleukin-12 production, receptor binding, positive regulation of MAP kinase activity, positive regulation of JNK cascade, immune response and positive regulation of humoral immune response. In the KEGG analysis, 14 pathways enriched by genes in the yellow module were shown in Supplementary Table 2, and the top 10 KEGG terms (Figure 6B) included viral protein interaction with cytokine and cytokine receptor, cytokine-cytokine receptor interaction, toxoplasmosis, measles, tuberculosis, chemokine signaling pathway, lipid and atherosclerosis, Toll-like receptor signaling pathway, Hepatitis B and JAK-STAT signaling pathway.

Figure 6 Functional enrichment analysis results of hub genes. A: The top 10 gene ontology terms of hub genes; B: The top 10 Kyoto Encyclopedia of Genes and Genomes pathways of hub genes. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Survival analysis

Additional survival analysis was conducted on the hub genes to evaluate their impact on BL patients' survival. Due to the small sample size of GSE4475, we opted for GSE69051 for survival analysis (Figure 7). Two of the 10 hub genes were significantly associated with OS: CXCL10 (P = 0.029, Figure 7F) and IL2RA (P = 0.0066, Figure 7G).

Figure 7 Kaplan–Meier survival curve. A to J: Kaplan–Meier survival curve of identified hub genes in GSE69051.

Establishment of the nomogram and assessment of predictive value

Based on the hub genes and clinical data of the patients, a nomogram was developed to predict the 1- and 3-year OS of BL patients (Figure 8A). The model had a c-index of 0.84, and the calibration curve demonstrated strong agreement between predicted and observed survival times for both 1- and 3-year OS probabilities in the BL cohort (Figure 8B).

Figure 8 Nomogram and calibration plot for GSE69051 cohort. A: The nomogram was constructed to predicting1, 3-year survival rate of Burkitt lymphoma patients; B: The calibration curves for predicting patient survival at 1 and 3 years in the cohort. OS: Overall survival.

Identification of associated drugs

The drugs related to IL2RA and CXCL10 were identified by the DGIdb database, as these were the only significant results from survival analysis (Figure 9). These results may provide new ideas for the treatment of BL with poor prognosis.

Figure 9 Drugs related to IL2RA and CXCL10.

DISCUSSION

BL, a highly aggressive lymphoma identified and described by Denis Burkitt in the last century, continues to be the most common childhood malignancy in Africa nowadays[22]. A defining feature of BL is the translocation between the c-MYC gene and the IgH gene, which is found in 80% of cases [t (8; 14)], or between c-MYC and either the kappa or lambda light chain gene, which is found in the remaining 20% [t (2; 8) or t (8; 22)][23]. The proliferation rate and apoptosis rate of BL tumor cells are extremely high, indicating that nearly 100% of the cells are positive for Ki-67. Intensive, short-course combination chemotherapy is recommended for most BL patients. DA-EPOCH (dose-adjusted etoposide, doxorubicin, cyclophosphamide, vincristine, prednisone) + rituximab may be an option for patients who cannot tolerate more aggressive regimens[7]. As the standard of treatment for BL has not yet been established, strictly controlled clinical trials are also recommended. The prognosis of BL patients is associated with both clinical and laboratory characteristics[8,24]. The BL international prognostic index can be used to assess the prognosis of adult patients with disseminated or immunodeficiency-related BL, but it is not currently used for stratifying BL treatment[20]. Previous studies have demonstrated associations between MYC rearrangements, TCF3 mutations or ID3 alterations (its negative regulator), TP53 modifications, CCND3 and CDKN2A changes, as well as non-antigen-dependent B cell receptor signaling (tonic B cell receptor signaling) with the development and prognosis of BL; however, a comprehensive investigation into the prognostic significance of molecular events associated with BL is lacking[25].

As a bioinformatics algorithm, WGCNA offers numerous advantages over conventional methods for differential expression analysis. It primarily focuses on elucidating co-expression patterns, facilitating the identification of biologically relevant modules comprising interconnected genes, and enabling the detection of pivotal hub genes[26-28]. So far, gene modules related to several cancers have been analyzed and verified using WGCNA[29,30].

In this study, 8 modules were obtained through WGCNA. As many prior studies have shown a strong correlation between age and the prognosis of BL patients[20,21], we chose the yellow module that had the strongest correlation with age for further analysis. Ten hub genes (SRC, TLR4, CD40, STAT3, SELL, CXCL10, IL2RA, IL10RA, CCR7 and FCGR2B) were identified using cytoHubba. GO and KEGG functional analyses were conducted on hub genes, and the PPI and GI analysis of these hub genes revealed their related biological functions. Based on survival analysis, CXCL10 and IL2RA have been identified as genes that affect survival. Afterward, we used a nomogram to develop a new risk assessment system for BL patients based on the aforementioned genes and age, aiming to aid in identifying high-risk groups for this disease.

CXCL10 is one of the three ligands for CXCR3, which is a chemokine receptor[31]. Various studies have demonstrated that in addition to attracting CD8+ and CD4+ effector T cells to tumor sites and sites of inflammation, CXCL10 also governs the polarization and enhances the biological functionality of these cells. This makes CXCL10 a key chemokine driver and a valid target for the therapy of autoimmune diseases such as Inflammatory Bowel Disease, Multiple Sclerosis, Rheumatoid arthritis, and others. Previous studies have also found that chemokines and their receptors are involved in supporting tumor development and metastatic spread[32-35]. In addition to inducing effector TH1 cells, CXCL10 has recently been proven to be associated with the recruitment of CXCR3+ CD8+ T cells to the tumor site and the induction of Granzyme B production by these cells, thereby enhancing their anti-tumor activities[36]. Barreira da Silva et al[37] used Dipeptidyl peptidase 4 inhibitors to increase the endogenous level of CXCL10, thereby suppressing experimental melanoma. It has also been demonstrated that the combination of CXCL10 gene therapy and radiotherapy improves therapeutic efficacy in cervical cancer using a HeLa cell murine xenograft tumor model[38]. Numerous studies have demonstrated a positive correlation between increased expression of CXCL10 at the tumor site and improved prognosis in various human cancers[39-41]. However, the biological functions of CXCL10 in BL have not been addressed so far. Our study initially discovered that the high expression of CXCL10 appeared to be associated with a better prognosis. In our prognostic model, CXCL10 outperforms age, which is an accepted prognostic factor for BL. Further studies are required to investigate and validate the mechanism of CXCL10 in BL.

IL2RA (CD25) is a low-affinity receptor for its ligand interleukin 2 (IL2), but when combined with IL2RB (CD122) and IL2RG (CD133), it forms the high-affinity IL2 receptor[42]. The binding of IL2 to IL2 receptor activates JAK1 and JAK3, which in turn activate several pathways that regulate cell survival and proliferation, including the PI3K/AKT, RAS/RAF/MEK/ERK, and STAT5 pathways[43]. IL2RA expression is elevated in a variety of cancers, especially hematologic tumors[44-46]. Fell et al[47] studied 69 patients with leukemia, lymphoma, or multiple myeloma and found that the expression of IL2RA in T cells was associated with frailty independent of age. This means that patients with high IL2RA expression showed better tolerance to chemotherapy and thus might have a superior prognosis. However, another study demonstrated that high IL2RA mRNA expression was an independent and adverse prognostic factor in acute myeloid leukemia, specifically stratifying patients into a worse prognosis[48], while reports on IL2RA in chronic myelogenous leukemia (CML) were controversial, they described it as either a promoter or an inhibitor of CML cell proliferation and disease aggressiveness[45,46]. This study demonstrated that BL patients with high expression of IL2RA exhibited a better prognosis. Due to the controversial reports on the function of IL2RA and the lack of research on BL, further studies are required to validate the prognostic value of IL2RA in BL.

As a predictive statistical tool, a nomogram visually displays the significant factors that impact outcomes in multifactor regression analyses and simplifies survival probability prediction through an easy-to-understand graphical representation[49]. The construction of the nomogram model in this study is based on age, IL2RA, and CXCL10. The nomogram effectively visualizes the impact of identified hub genes and facilitates survival prediction, with the multivariate regression analysis serving as the fundamental component of this model. However, the nomogram would benefit from a validation cohort to enhance its current model. Therefore, it is recommended that additional patients with long-term follow-up be included in future studies.

Based on CXCL10 and IL2RA, we have also identified some drugs that may potentially play a therapeutic role in relapsed and refractory BL, which require further research on pharmacology and treatment protocols. There are also some limitations of the present study. Firstly, the sample size may not be sufficient and could result in selection bias. Secondly, the three different clinical types of BL have the same histological features and similar clinical behavior but differ in epidemiology, clinical presentation, and genetic characteristics, which might need to be classified and analyzed separately. What's more, additional genetic and experimental studies are required to explain the mechanism and the function of these hub genes in the carcinogenesis and progression of BLs. Due to the limited experimental conditions, our study exclusively utilized data sourced from publicly available databases. However, further validation is needed in larger samples or more external datasets.

CONCLUSION

In conclusion, this study is the first to investigate gene expression in BL using WGCNA. These findings provide a framework for identifying co-expression gene modules and discriminating key pathways and hub genes in BL. In the present study, we identified and verified 10 hub genes. Survival analysis showed that overexpression of CXCL10 and IL2RA in BL may serve as superior prognostic indicators. Additionally, an integrated mRNA signature and age nomogram potentially offer prognostic value for patients with BLs.

ARTICLE HIGHLIGHTS

Research background

Burkitt lymphoma (BL) is an exceptionally aggressive malignant neoplasm originating from either the germinal center or post-germinal center B cells. However, a standardized treatment regimen for BL has yet to be established. The utilization of microarray data and sequencing information retrieved from public databases presents promising prospects for the identification of novel diagnostic or therapeutic targets.

Research motivation

It is crucial to identify biomarkers that can predict the prognosis of BLs and distinguish patients who would benefit from specific therapies.

Research objectives

The aim of our study was to identify hub genes and conduct gene ontology analysis specifically in BL, as well as perform functional enrichment analysis. Additionally, we performed survival analysis and developed a novel prognostic model incorporating candidate genes along with clinical features.

Research methods

The gene expression profiles and clinical traits of BL patients were obtained from the Gene Expression Omnibus database. Weighted gene co-expression network analysis (WGCNA) was employed to construct gene co-expression modules, while the cytoHubba tool was utilized to identify hub genes. Prognostic candidate genes were identified through overall survival (OS) analysis. A nomogram was developed to evaluate the predictive value of the hub genes.

Research results

In this study, we identified 8 modules through WGCNA analysis and found a significant correlation between the yellow module and age. By using the cytoHubba tool, we identified 10 hub genes (SRC, TLR4, CD40, STAT3, SELL, CXCL10, IL2RA, IL10RA, CCR7, and FCGR2B). Among these hubs, two genes (CXCL10 with P = 0.029 and IL2RA with P = 0.0066) were associated with OS based on our survival analysis.

Research conclusions

This study is the first to investigate gene expression in BL using WGCNA. We have identified and validated 10 hub genes, demonstrating that the overexpression of CXCL10 and IL2RA in BL can serve as robust prognostic indicators. Furthermore, the integration of an mRNA signature with age nomogram holds promising potential for predicting patient outcomes in BLs.

Research perspectives

Further genetic and experimental investigations are imperative to elucidate the underlying mechanism and functional significance of these hub genes in the carcinogenesis and progression of BLs.