Integrated bioinformatics analysis reveals dynamic candidate genes and signaling pathways involved in the progression and prognosis of diffuse large B-cell lymphoma

Data collection

We evaluated and downloaded the raw gene expression profiles from the National Center for Biotechnology Information-Gene Expression Omnibus (NCBI-GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The series was based on GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. GSE12453 (Brune et al., 2008) was used to identify differentially expressed and co-expressed genes. It was also used to predict tumor-infiltrating immune cells (TIICs). It contained 11 DLBCL cases and 25 normal controls. The controls were non-neoplastic B lymphocytes isolated from healthy donors’ blood or tonsils from routine tonsillectomy patients. Similar expression profiles on the selected GPL570 platform contained normal but reactive controls or non-human controls, such as cell lines or limited controls, and were excluded. We used the GSE31312 (Visco et al., 2012) with 498 DLBCL cases for survival analysis and identifying immune cells infiltrating the TME.

Study design and data pre-processing

The GSE12453 CEL file was pre-processed using the Affymetrix package (Gautier et al., 2004) in R software (https://www.r-project.org/) version 4.0.2 (Team, 2019). The procedures included background correction, log₂ transformation, followed by quantile normalization. We performed a standard quality assessment, including scaling factors and NUSE plots, and hierarchical clustering to identify outliers. The study’s design flowchart is shown in Fig S1.

Screening of differentially expressed genes (DEGs)

The DEGs between DLBCL cases and non-cancer controls were screened with R package limma (https://www.bioconductor.org/), Release 3.11 (Ritchie et al., 2015). The cutoff criteria were p < 0.05, |log₂fold change (FC)| > 1.0. The pheatmap package was used to generate hierarchical clustering, and ggplot2 (Wickham, 2016) was used to show the volcano plot in R. We used the resulting data output tables that included gene ID, log₂FC, unadjusted and adjusted p-values in subsequent analyses.

Functional and pathway enrichment analysis

We investigated the functional roles and pathway signaling relevance of the DEGs. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the R clusterprofiler (https://bioconductor.org/), Release 3.12 (Yu et al., 2012). A p-value < 0.05 was considered significantly enriched. The GO categories included biological process (BP), molecular function (MF), and cellular component (CC).

Gene set enrichment analysis (GSEA)

We used the normalized expression dataset of GSE12453 for the GSEA (www.gsea-msigdb.org/gsea/index.jsp), version 4.1.0. We followed the recommended protocol. This included “gct” and “cls” file formats for the expression dataset and phenotype labels, respectively (Subramanian et al., 2005). Significant gene sets had a false discovery rate (FDR) <0.25 and a nominal p-value < 0.05.

Protein-protein interaction (PPI) network of DEGs and hub modules selection

We applied the STRING database (https://string-db.org), version 11.0 (Szklarczyk et al., 2019), to assess and map the identified DEGs into a human PPI network. Then, we uploaded the resulting data output into Cytoscape software (https://cytoscape.org/), version 3.8.2 (Shannon et al., 2003), and used MCODE with default parameters (Bader & Hogue, 2003) and the Cytohubba (Chin et al., 2014) plug-ins to select significant modules and top-ranked genes. Nodes and edges represented genes and their interactions.

Predicting tumor-infiltrating immune cells (TIICs)

We applied the Cibersort algorithm (Newman et al., 2019) in R to predict TIICs using the normalized GSE12453 dataset, according to the CIBERSORT instructions. We used data with a p-value < 0.05 for further analysis.

Co-expression network construction and identification of modules related to DLBCL

We conducted gene co-expression analysis on the processed GSE12453 data using the weighted correlation network analysis (WGCNA) (Langfelder & Horvath, 2008; Langfelder & Horvath, 2012) in R (Team, 2019). We followed the standard protocols, including quality control procedures. We performed a power β transformation on the computed Pearson correlation matrix to ensure a scale-free topology. The minimum number of module genes was set at 30. The WGCNA R package then generated a co-expression network from the resulting adjacency matrix. We applied the dynamic Tree Cut package (Team, 2019) to create the co-expression modules from the color-coded hierarchical clustering dendrogram. We assessed clinical module-trait relationships with Pearson’s correlation. Gene significance (GS) and module membership (MM) were also analyzed for their correlation in modules. Statistically significant modules were defined as p < 0.05.

Identifying differentially co-expressed hub module genes

The hub genes were identified from the intersection between the DEGs and genes significant in the WGCNA modules. We then analyzed the hub genes with the STRING and Cytoscape databases. We further conducted functional and pathway analysis on the hub genes to determine the relevant genes that impact DLBCL.

Validation and survival analysis of core hub genes

We validated the expression levels of the hub genes in the GEPIA2 (gepia2.cancer-pku.cn/#index) database (Tang et al., 2019). We retrieved the pre-processed quantile normalized series matrix file (GSE31312, Affymetrix HG-U133 Plus 2.0 GeneChips) containing 498 de-novo adults DLBCL from the NCBI-GEO database. We used it for the prognostic value analysis (overall survival (OS) and progression-free survival (PFS)). We plotted Kaplan–Meier survival curves with ggplot2 in R. The genes with p-value < 0.05 were selected as core hub genes. Also, receiver operating characteristic (ROC) curves were plotted with the pROC R package (Team, 2019) to validate the diagnostic value of the core genes. Then, the human protein atlas (proteinatlas.org/) was used to confirm their protein expressions in some selected lymph node samples (Uhlen et al., 2017).

We used the ssGSEA (single-sample Gene Set Enrichment Analysis) package (Subramanian et al., 2005) in R on the GSE31312 dataset to investigate immune-associated core hub genes in the TME.

In addition, we also employed the GEPIA2 and Oncomine (Rhodes et al., 2004) databases to determine the expression levels and relevance of the core hub genes in other tumors.

Differentially expressed genes (DEGs)

According to standard protocols, one poor quality control sample (GSM312887) was excluded from the GSE12453 dataset after pre-processing. We identified 1,260 DEGs between the 11 DLBCL and 24 non-tumor controls. They comprised 1,014 up-regulated and 246 down-regulated genes (p < 0.05 and |log₂FC| > 1). The gene list, Affymetrix probe ID, and log FC are shown in (Excel S1). The expression patterns of the DEGs are shown in Fig. 1. Figure 1A shows the volcano plot of all expressed genes. The DLBCL cases showed a distinctive gene expression profiling. As shown in Fig. 1B, the heatmap of the top 169 DEGs with |log₂FC| > 2 suggested that the identified DEGs expression levels could differentiate DLBCL from non-tumor samples. We utilized the CytoHubba application in Cytoscape, employing five calculation methods: the Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), Degree, Edge Percolation Component (EPC), and EcCentricity to rank the top 250 DEGs. The genes from the five methods were intersected using the Venn diagram software (http://bioinformatics.psb.ugent.be/webtools/Venn/). Most of the intersecting genes (Fig. 1C) were associated with significant and valuable pathways such as the proteasome, spliceosome, and viral protein interaction with cytokine and cytokine receptor (Fig. S2). Intersecting genes are common genes with a high degree of interconnection and are more likely to represent key candidate genes with important biological regulatory functions.

Functional and pathway analysis

We applied the enrichGO or enrichKEGG function of the clusterProfiler package (Yu et al., 2012) in R to investigate the biological functions of all the DEGs (p < 0.05). The GO and KEGG pathway analysis results are shown in Fig. 2. In GO analysis, the DEGs were primarily enriched in ATP synthesis coupled electron transport and mitochondria ATP synthesis coupled electron transport for biological processes (BP). Their cellular components (CC) were mainly related to the mitochondria protein complex, mitochondria inner membrane, etc. Their molecular functions (MF) consisted of structural constituents of ribosomes and NADH dehydrogenase activity (Fig. 2A). The KEGG pathway analysis showed significant enrichment in ribosome and oxidative phosphorylation (Fig. 2B). Strikingly, the DEGs involved in the KEGG pathway and GO enrichment analysis were mainly upregulated genes, with few down-regulated genes.

Further analysis showed that the down-regulated genes were enriched in B cell activation (adjust p-value = 0.003) for GO-BP; our data showed no other significant GO and KEGG enrichment for the down-regulated genes. The GO and KEGG enrichment analysis for the upregulated genes was similar to the analysis we performed earlier. The GO and KEGG analysis classifies genes into functional categories to help understand their functions and regulatory pathways. Hence, with additional investigations, the genes in these pathways might throw more light on the pathogenesis of DLBCL.

Gene Set Enrichment Analysis (GSEA) of DLBCL expression dataset (GSE12453)

We performed a GSEA analysis to compare the DLBCL and non-tumor controls’ expression profiles to understand better the biological functions of the relevant genes discovered. We analyzed all the qualified genes in the GSE12453 expression dataset. The KEGG output from the GSEA was similar to our previous pathway analysis and confirmed our earlier results. The Hallmark gene sets showed immune and metabolic-related signaling predominance, including estrogen response late, epithelial-mesenchymal transition, and UV response up (Table S1). Interestingly, genes defining late response to estrogen, epithelial-mesenchymal transition (such as in fibrosis and metastasis), and genes upregulated in response to ultraviolet (UV) radiation have not been fully elucidated in DLBCL. However, evidence suggests these genes have essential roles in oncogenesis. These findings provide evidence that can drive future research with therapeutic implications.

Predicting tumor-infiltrating immune cells (TIICs)

Our GO, KEGG, and GSEA analyses showed that the DEGs were enriched in some immune-related biological functions. Immune cell infiltration into tumors plays an essential role in tumorigenesis and metastasis. So, we applied the CIBERSORT algorithm to the GSE12453 dataset to predict immune cells infiltrating the TME. Among the immune subsets analyzed, activated memory CD4⁺ T cells, CD8⁺ T cells, and M0, M1, and M2 macrophages were the most represented cell fractions within the DLBCL microenvironment (Figs. 3A, 3B). The correlations among the TIICs ranged from high to negligible (Fig. S3). Regulatory T cells (Tregs) showed a moderate negative correlation with activated memory CD4⁺ T cells. However, M0 macrophages had a high positive correlation with activated memory CD4⁺ T cells and M1 macrophages. The varied infiltrating immune cell types could be reflective of the complexity and the unusual behavior of DLBCL. Uncommitted macrophages (M0) can polarize into the M1 (considered tumoricidal) and M2 (pro-tumorigenic) to show paradox effects on tumor prognosis (Dancsok et al., 2020). Memory B cells and activated dendritic cells were the most represented fractions in the non-tumor controls (Fig. 3A). The findings suggest that TIICs may be closely associated with clinical outcomes. Future studies, including their correlation to DLBCL disease stages, will be meaningful, particularly for immunotherapy.

Protein-protein interaction (PPI) network and hub modules establishment

We constructed the PPI network of all DEGs using the STRING database’s multiple proteins function to determine genes likely to perform biological functions together. The highest confidence of 0.9 was set with unconnected nodes taken out. As shown in Fig. 4A, it yielded 609 nodes and 8,583 edges. These genes were highly inter-connected than expected (PPI enrichment p-value < 1.0e−16). We observed six important clusters with a k-score > 10 when we analyzed the tab-separated values (tsv) file with Cytoscape’s MCODE. The largest cluster (#1) had the highest score of 33.17 (Fig. 4B); it mainly comprised Ribosome genes (FDR = 3.49e−83) in KEGG analysis. Cluster 2 (Fig. 4C) was enriched in genes associated with oxidative phosphorylation (FDR, 7.23e−56) and Parkinson’s disease (FDR, 2.62e−55). Cluster three (Fig. 4D) genes were mainly involved in chemokine signaling (FDR, 1.06e−22). We identified 109 highly connected DEGs (Excel S1) from these top three clusters. These DEGs could play essential roles in DLBCL, so they were selected for further hub gene screening.

In this study, the down-regulated DEGs were not part of the constructed PPI, so we determined their relevance in DLBCL. We ranked the top 60 DEGs (30 up-and 30 down-regulated) by log₂FC and analyzed them in the STRING database (Fig. S4A). The only down-regulated DEG in the network built was identified upregulated when verified in the GEPIA2 database. To further investigate the down-regulated DEGs’ functional relatedness, a PPI was constructed for all the down-regulated DEGs (Fig. S4B). We found that the down-regulated genes were enriched in B cell activation for GO-BP, shown in (Fig. S4C). Finally, we used the degree method of the CytoHubba application to predict the top 250 important genes (Fig. 4E). Genes with a high degree of centrality are vital since they have many direct interacting gene partners. If confirmed, these critical findings could improve the general understanding and the potential causes of variation in the clinical prognosis of DLBCL.

Weighted gene co-expression network (WGCNA) analysis

To identify co-expression modules that could share similar biological functions or regulatory mechanisms with clinical relevance to DLBCL, we applied the WGCNA package (Langfelder & Horvath, 2008; Langfelder & Horvath, 2012) in R (Team, 2019). The GSE12453 dataset we processed was used. We carried out quality control procedures, including inspecting good genes and sample hierarchical clustering to detect potential outliers, but no obvious outliers were found (Fig. 5A). The 35 samples yielded two main clusters. We applied the WGCNA on the top 25% of the 21,654 expressed genes ranked by the largest variance. To satisfy a scale-free network topology, we choose the soft-threshold power β of eight with R2 = 0.86 (Figs. 5B, 5C). Hierarchical clustering and the dynamic tree-cutting yielded 18 modules of co-expressed genes (Fig. 5D). Finally, we visualized the top 1,000 significantly expressed genes with a heatmap (Fig. 5E); they represent interesting genes for further analysis.

To investigate the molecular mechanisms of the traits, we correlated each Module Eigengenes (ME) to disease status (DLBCL and non-tumor controls). The results are shown in Fig. 6. The ME turquoise and green (Figs. 6A, 6B) containing 292 and 72 genes, respectively, strongly correlated with DLBCL. ME dark magenta with five genes had the strongest negative correlation. The cut-off was set at gene significance (GS) value >0.8, and absolute Module Membership (MM) value >0.7. Besides, the GS versus MM plots for these three modules were highly correlated (Fig. 6C), reflecting their high association with DLBCL. We selected these three clinically significant modules with the 369 high connectivity genes (gene list shown in Excel S1) for further analysis. The highly connected genes are often the most important (central) elements of the respective modules and tend to play key roles in the biological processes. The ME genes are listed in Excel S1. Altogether, these co-expressed genes might provide new clues to understand the biology of DLBCL in the future.

PPI and functional enrichment analysis of the WGCNA relevant modules

The 369 high connectivity genes from the three relevant modules were filtered in the STRING database followed by the Cytoscape; the network yielded a PPI with 195 nodes and 1,579 edges. The PPI and gene list are detailed in Fig. S5 and Excel S1. These 195 genes were considered functionally important. As presented in (Fig. 7), functional annotation revealed that these genes were involved in viral transcription and viral gene expression in the BP category. In KEGG analysis, the genes were primarily enriched in ribosome, coronavirus disease-COVID-19, and oxidative phosphorylation. The identified pathways were roughly consistent with that of the DEGs. These processes and signaling pathways are usually disrupted in cancer and could provide an insight into the pathogenesis of DLBCL.

Identification of hub genes and pathways

Eventually, we identified 46 important differentially co-expressed genes by the Venn diagram software, as shown in (Fig. 8A). These 46 genes were common between the DEGs and WGCNA hub module genes and were regarded as hub genes. We re-analyzed the 46 genes with the STRING and Cytoscape databases, and the PPI is shown in (Fig. 8B). Their GO and KEGG enrichments in the R software (Table 1) were similar to the other analyses. The KEGG common genes in the ribosome, COVID-19, and oxidative phosphorylation pathways are shown in Table 2. The above pathway genes play essential roles in metabolic reprogramming and tumor-promoting inflammation of cancer and warrant further studies.

Validation of expressions and prognostic analysis of core hub genes

We applied GEPIA2 to validate the reliability and authenticity of the 46 hub genes in the cancer genome atlas (TCGA) dataset. We identified 44 prognostic genes with higher expression (consistent with that in the GSE12453 dataset) in DLBCL tissues than the non-tumor control tissues (p < 0.01) (Fig. 9) (15 genes shown) and (Excel S1). Kaplan–Meier survival analysis on the GSE31312 showed 15 of these genes (p < 0.05) correlated with patient outcomes (Fig. 10 & Table S2). Except for RPL11, the patients with high expressions had significantly shorter 5-year OS and PFS, suggesting these genes are potential oncogenes and have a role in DLBCL development and/or progression. The 15 genes (Fig. 11A) were considered as the core hub genes. Moreover, ROC curve analysis for their diagnostic potentials obtained AUCs ranging from 0.992 to 1.00, indicating optimal performance to accurately differentiate DLBCL from non-tumor control cases (Fig. S6). Also, immunohistochemistry data from the human protein atlas (HPA) database demonstrated the protein expressions of some of the genes in some lymph node samples with cytoplasmic/membranous localization (Fig. S7, four genes shown); data were retrieved from https://www.proteinatlas. The genes included RPS21, MRPS28, RPL31, and RPL30. As expected, they would be involved in metabolic pathways such as glycolysis and processes including signal transduction and cell division.

The core hub genes’ functional annotation was mainly associated with Ribosome and Coronavirus disease-COVID-19 (Figs. 11B, 11C; Table 3). To assess the tumorigenic potentials of the COVID-19 genes regarding immune cells infiltrating the TME, the ssGSEA analysis was applied (GSE31312). As shown in Fig. S8, seven out of the nine COVID-19 pathway genes negatively correlated with mast cells, five with immature dendritic cell (iDC), and three genes negatively correlated with plasmacytoid DC. RPL30, RPL31, RPS25, and FAU positively correlated with tumor-infiltrating lymphocytes (TIL) and macrophages. RPL30, RPL31 correlated with Tregs. These infiltrating immune cells may be involved in regulating tumor proliferation, dormancy, and drug resistance.

Finally, we determined whether the core genes were upregulated in other tumors. The GEPIA2 database revealed that all the 15 core genes were upregulated in thymoma (THYM), and 11 genes were upregulated in testicular germ cell tumors (TGCT). Notably, RPL30 and FAU genes were consistently upregulated in all six different cancer types identified (Table S3). In the Oncomine database, some of the core genes were upregulated in various lymphoma datasets and other cancers, including Sarcoma (Fig. S9). The results suggest that the upregulation of these 15 hub genes may not be limited to DLBCL.