Potential Molecular Mechanism of TNF Superfamily-Related Genes in Glioblastoma Multiforme Based on Transcriptome and Epigenome

Objective: This study aimed to investigate the molecular mechanism of tumor necrosis factor (TNF) superfamily-related genes and potential therapeutic drugs for glioblastoma multiforme (GBM) patients based on transcriptome and epigenome.

Methods: Gene expression data, corresponding clinical data, and methylation data of GBM samples and normal samples in the TCGA-GBM and GTEx datasets were downloaded. The TNF-related genes were obtained, respectively, from two groups in the TCGA dataset. Then, the TNF-related differentially expressed genes (DEGs) were investigated between two groups, followed by enrichment analysis. Moreover, TNF superfamily-related gene expression and upstream methylation regulation were investigated to explore candidate genes and the prognostic model. Finally, the protein expression level of candidate genes was performed, followed by drug prediction analysis.

Results: A total of 41 DEGs including 4 ligands, 18 receptors, and 19 downstream signaling molecules were revealed between two groups. These DEGs were mainly enriched in pathways like TNF signaling and functions like response to TNF. A total of 5 methylation site-regulated prognosis-related genes including TNF Receptor Superfamily Member (TNFRSF) 12A, TNFRSF11B, and CD40 were explored. The prognosis model constructed by 5 genes showed a well-prediction effect on the current dataset and verification dataset. Finally, drug prediction analysis showed that zoledronic acid (ZA)-TNFRSF11B was the unique drug–gene relation in both two databases.

Conclusion: Methylation-driven gene TNFRSF12A might participate in the development of GBM via response to the TNF biological process and TNF signaling pathway and significantly associated with prognosis. ZA that targets TNFRSF11B expression might be a potential effective drug for clinical treatment of GBM.

Introduction

Glioblastoma multiforme (GBM) is the most aggressive cancer that represents 15% of all brain tumors (1). The most common length of survival following diagnosis is 12 to 15 months, with fewer than 3% to 5% of people surviving longer than 5 years (2). Typically, surgery after chemotherapy and radiation therapy are commonly used for the treatment of GBM (3). However, the cancer usually recurs due to poor effect of existing drugs or treatment strategies on the diffusive, infiltrative, and metastatic of GBM (4).

Further understanding of the molecular mechanisms involved in the development of GBM may contribute to the development of new therapies and strategies (5). Recent observations after immunotherapies with cytokines suggest an immunological and even clinical response from immunotherapies (6). Actually, members of the tumor necrosis factor (TNF) superfamily (TNFSF) and TNF receptor superfamily (TNFRSF) have crucial roles in both innate and adaptive immunity (7). A previous study shows that tumor necrosis factor (TNF) and the associated receptor superfamily play important roles in the development of GBM (8). Some TNFs such as TNF-α are upregulated in GBM cells, which further play an important role in GBM progression (9). TNFRSF19 is upregulated in advanced glial tumors and promotes glioblastoma cell invasion (10). Furthermore, DNA methylation plays an important role in gene expression regulation during the development of tumor (11). Abnormal epigenetic modification can lead to tumors, genetic disorders, inflammation, aging, and neuropsychiatric abnormalities (12). A previous study shows that epigenetic therapy with inhibitors of histone methylation suppresses DNA damage signaling and increases glioma cell radiosensitivity (13). Methylation profiling can be used to identify different groups of GBM according to their tumorigenesis (14). However, due to the lack of integrated analysis of epigenomic and transcriptome data, the specific role of DNA methylation sites and TNF-related gene expression changes in GBM progress, as well as potential effective drugs associated with these genes, are still unclear.

In this study, the gene expression data, clinical information, and methylation data of GBM tumor samples and normal tissue samples in TCGA-GBM dataset and GTEx dataset were downloaded. The differentially expressed genes (DEGs) were explored between tumor and normal samples, followed by function and pathway enrichment. Then, the expression of TNF superfamily-related genes and upstream methylation regulation mechanism was investigated to explore candidate genes and prognostic models. Finally, the protein expression level of candidate genes was performed, followed by drug prediction analysis. This study hoped to investigate the molecular mechanism of TNF superfamily-related genes and potential therapeutic drugs for GBM patients.

Materials and Methods

Data Acquisition

A total of 96 TNF superfamily-related genes including 18 TNF superfamily (TNFSF), 29 TNF receptor superfamily (TNFRSF), and 49 related downstream signal genes were enrolled for the current analysis based on literature review (15). The RNA-seq data including corresponding methylation and clinical phenotype data of GBM samples in TCGA and GTEx were downloaded from the University of California Santa Cruz (UCSC, http://xena.ucsc.edu/) Genome Browser database (16). A total of 166 GBM samples obtained from the TCGA-GBM dataset were enrolled as the tumor group. Meanwhile, 5 normal paracancerous tissue samples from the TCGA-GBM dataset and 105 normal brain cortex samples from the GTEx dataset were combined as the normal group. Moreover, after being annotated by hg38 gene annotation information (gencode.v23.annotation.gene.probemap) based on the GENCODE database (17) and having filtered low-level expression genes (not expressed in half of all samples), a total of 76 TNF superfamily-related genes (Supplementary Table 1) including 7 TNFSF, 21 TNFRSF, and 48 downstream signal molecules were extracted for subsequent analysis.

The methylation data of GBM patients was obtained by the correspondence between the chip methylation spectrum site and the symbol annotation conversion file downloaded from UCSC-Xena. Furthermore, clinical data including age, race, gender, radiotherapy and chemotherapy information, new tumor information, OS status, and OS time of each patient in the downloaded data were obtained using TCGAbiolinks in R (18).

Differentially Expressed Analysis

The linear regression and empirical Bayesian methods (19, 20) in limma package of R (21) were used to explore DEGs between the tumor group and normal group based on the TNF-related 76 gene expression matrix. The Benjamin & Hochberg method was used for multiple-test correction. The adjP < 0.05 and | log2 fold change (FC)| > 1 were selected as the thresholds for DEG screening. Then, the volcano plots and clustering heat map were constructed using ggplot2 (version: 3.2.1) (22) and using pheatmap (version: 1.0.12) (23), respectively.

Enrichment Analysis of DEGs

Gene ontology-biological process (GO-BP), GO-cellular component (GO-CC), and GO-molecular function (GO-MF) (24), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (25) enrichment analyses of TNF-related DEGs were performed using the Metascape software (parameter: min overlap = 3; P-value cutoff = 0.01; min enrichment = 1.5) (26). P < 0.01 was considered as the cutoff value of significant enrichment results. Moreover, clustering analysis was conducted according to the similarity of genes enriched in each term (similarity of > 0.3). The most statistically significant term (P-value minimum) in each cluster was selected to define this cluster. Then, the top 20 clusters based on the P-value were visualized by a histogram. Furthermore, to further explore the relationship between terms, the interaction network diagram of terms in the top 20 clusters was constructed (inclusion criteria: terms with the most significant P in each cluster; <15 terms in each cluster; no more than 250 terms in total; similarity > 0.3). Finally, the network was constructed by Cytoscape software (version: 3.4.0) (27).

Correlation Between Methylation Level and DEG Expression

Based on the correspondence between methylation sites and genes, all methylation sites corresponding to the differential genes were extracted, and the Pearson correlation coefficient (r) between each site and its corresponding gene expression level was calculated and tested for significance. Finally, the P < 0.05 and r < −0.4 were selected as cutoff values to screen the methylation-related genes.

Prognosis Analysis Based on DEGs and Methylation Sites

The expression value of DEGs in each sample and the associated patient clinical survival information were used for the DEG prognosis analysis. Univariate Cox regression analysis was used to analyze the associations between DEGs and prognosis, and the DEGs with P < 0.05 were considered as the prognosis-related genes. Meanwhile, the methylation value of methylation sites that correspond to prognosis-related genes and the associated patient clinical survival information were used for the methylation site prognosis analysis. The univariate Cox regression analysis was used to analyze the relationship between each methylation site and prognosis. P < 0.05 was considered as the cutoff value for candidate prognostic methylation sites.

Prognostic Model Construction and Verification

The prognostic gene that significantly negatively correlated with the prognostic methylation site was considered as the methylation site-regulated prognosis-related gene. Then, these genes were screened by multivariate Cox expression regression. Based on the prognostic correlation coefficient β and the combination of the expression values of selected genes, the risk score calculation model was defined as

$\begin{array}{l}Riskscore={\beta }_{gene1}\times exp{r}_{gene1}+{\beta }_{gene2}\times exp{r}_{gene2}\\ +\dots +{\beta }_{geneN}\times exp{r}_{geneN}\end{array}$

The corresponding risk score of each sample was calculated, and the samples were divided into high-risk group or low-risk group based on the median risk score. To reveal the relationship between high/low-risk group and prognosis, the Kaplan–Meier (KM) (28) survival curve and heat map were used to assess the survival time distribution and gene expression value of the two groups. To validate the risk model, the expression profile data of WHO IV grade samples, including GBM, rGBM, and sGBM (DataSet ID: mRNAseq_325) (29, 30), were downloaded from the CGGA database (http://www.cgga.org.cn/download.jsp). Clinical information such as gender, age, chemoradiation information, OS status, and OS time in these data were further enrolled in this study.

Independence Analysis of Prognostic Models

To investigate whether the prognostic model could be independent of other clinical variables (including age, gender, etc.), univariate Cox regression analysis was performed based on independent variables including high/low-risk groups, age, and gender. Then, the factors with P < 0.05 were enrolled for the multivariate Cox regression analysis. All investigation was performed based on TCGA and CGGA datasets, followed by visualization with forest plots.

Protein Expression Level Verification

The Human Protein Atlas (HPA) is a database used to study protein expression in different human tissues and organs from RNA and protein levels by transcriptomics and proteomic techniques (31). In order to verify the difference in protein level of the key candidate genes, the HPA database was used to reveal the protein immunohistochemical level of key genes in cortex of normal people and GBM patients.

Drug–Gene Interaction Prediction

The drugs targeted by diseases-related genes were screened using the Drug–Gene Interaction database (DGIdb, version: 2.0) (32). Based on the drug–target gene relations, the drug–target gene interaction network was constructed by using online database STITCH (parameters: species = homo; medium confidence score = 0.4) (http://stitch.embl.de/) (33).

Results

DEGs Between Normal Group and Tumor Group

A total of 41 TNF-related DEGs including 4 TNFSFs, 18 TNFRSFs, and 19 downstream signal molecules were identified between the tumor and normal groups. The volcano plot showed that the upregulated genes and downregulated genes were significantly separated (Figure 1A). The heat map showed that the samples could be obviously distinguished by DEGs (Figure 1B).

FIGURE 1

Figure 1. The volcano plots and heat map for DEGs between the tumor sample and normal sample. (A) The volcano plots of DEGs; the X-axis represents the value of log2 fold change, while the Y-axis represents the value of –log10; red triangles represent the upregulated genes, blue squares represent the downregulated genes, and black nodes represent genes with no significant difference. (B) The heat map for DEGs; colors from blue to yellow indicated low to high representation values. The colored blocks at the top represent samples, of which brilliant blue represents tumor samples and pea green represents normal samples; the colored blocks at the left represent DEGs.

Significant GO Function and KEGG Pathways Enriched by DEGs

The obtained DEGs were significantly enriched in 103 GO-BP, 1 GO-CC, 19 GO-MF, and 32 KEGG pathways in the current functional enrichment analysis. These GO terms and KEGG pathways were clustered into different categories based on the similarity cluster analysis. The top 20 cluster is shown in Figure 2A. The results showed that these DEGs were mainly enriched in GO functions like response to tumor necrosis factor (GO: 0034612), death receptor activity (GO: 0005035), and tumor necrosis factor receptor superfamily binding (GO: 0032813). Meanwhile, these DEGs were mainly enriched in KEGG pathways including the TNF signaling pathway (hsa04668), apoptosis (hsa04210), and NF-kappa B signaling pathway (sha04064) (Figure 2A). Moreover, the investigation of the interaction among terms in each cluster is shown in Figure 2B. Each node represents a term, and the nodes with the same color represent the terms in the same cluster. As expected, the terms with more similarity were always clustered in a functional module, while the terms in different clusters showed less interactions.

FIGURE 2

Figure 2. GO/KEGG pathway enrichment cluster interaction analysis of the differentially expressed genes. (A) The X-axis represents the gene ratio (–log10); the Y-axis represents the different items of functions or pathways. (B) The interactive network among terms; different node colors indicated different clusters, and lines indicated gene similarities among terms.

The Methylation Site–DEG Interaction and Prognostic Gene Investigation

Based on the corresponding relationship between methylation sites and DEGs, all methylation sites corresponding to DEGs were extracted. A total of 504 methylation sites were obtained. The Pearson correlation coefficient (r) between each methylation site and the expression level of its corresponding gene was calculated to screen the methylation-related genes. A total of 39 methylation sites corresponding to 16 DEGs were obtained with the cutoff value of P < 0.05 and r < −0.4.

Furthermore, a total of 14 prognosis-related genes were obtained. The associations between methylation sites that correspond to prognosis-related genes and prognosis were calculated, and 74 methylation sites that correspond to 25 prognosis-related genes were obtained to be associated with prognosis (Supplementary Table 2).

Prognosis Model Constructed by Methylation-Driven Genes

A total of 5 methylation site-regulated prognostic genes including CD40, lymphotoxin beta receptor (LTBR), TNF receptor superfamily member (TNFRSF) 10C, TNFRSF 11B, and TNFRSF12A (Supplementary Figure 1) were revealed. Multivariate Cox regression was performed on these 5 candidate genes, followed by the risk model construction. The results showed that the survival time of the high-risk group was shorter than that of the low-risk group (Figure 3A). With the increase of the risk score, the expression level of these 5 candidate genes was relatively higher, and the survival rate of the high-risk group was significantly lower than that of the low-risk group (Figure 3B). The heat map of these 5 candidate genes in each sample is shown in Figure 3C. Moreover, the GBM samples in the CGGA database were used to evaluate the above risk model. The results showed that the prognosis effect of the risk model on the CGGA database were the same with that on the TCGA database (Supplementary Figure 2), which further indicated that the prognosis model constructed by these 5 candidate genes was effective.

FIGURE 3

Figure 3. Prognostic verification analysis for the current prognostic model based on tumor samples in the TCGA database. (A) Survival analysis for the high-risk group and low-risk group; survival time of the high-risk group was shorter than that of the low-risk group. The X-axis represents the overall survival time (month), while the Y-axis represents the survival rate (percent survival). P < 0.05 was considered to be significantly different. (B) The risk score and follow-up in the high-risk group and low-risk group. (C) The heat map for methylation site-regulated prognosis-related genes including CD40, LTBR, TNFRSF10C, TNFRSF11B, and TNFRSF12A.

Independence Analysis of the Prognosis Model

In order to investigate whether the prognosis model could be independent of other clinical variables, the univariate and multivariate Cox regression analyses based on the TCGA dataset and CGGA dataset were performed. In the TCGA cohort, on univariate Cox regression analysis of clinical valuables and risk score, the results showed that age, relapse or metastasis, drug therapy, radiation therapy, and risk score had a statistically significant impact (P < 0.05). These valuables were further included in multivariate analysis; the results showed that relapse or metastasis (HR 0.49, 95% CI 0.316–0.762, P = 0.002), radiation therapy (HR 0.276, 95% CI 0.127–0.599, P = 0.001), and risk score (HR 1.659, 95% CI 1.080–2.550, P = 0.021) were variables that independently affect survival (Table 1, Figure 4A).

TABLE 1

Table 1. The univariate and multivariate Cox regression analysis results for the TCGA dataset.

FIGURE 4

Figure 4. Forest map of regression analysis based on the TCGA and CGGA datasets. (A) The forest plot for multivariate Cox regression of the TCGA dataset on factors including age, relapse or metastasis, drug therapy, radiation therapy, and risk score. (B) The forest plot for multivariate Cox regression of the CGGA dataset on factors including radiation therapy, chemotherapy, and risk score.

In the CGGA cohort, on univariate Cox regression analysis of clinical valuables and risk score, the results showed that radiation therapy, chemotherapy, and risk score had a statistically significant impact (P < 0.05). These valuables were further included in multivariate analysis; the results showed that radiation therapy (HR 0.580, 95% CI 0.361–0.930, P = 0.024), chemotherapy (HR 0.434, 95% CI 0.288–0.654, P < 0.001), and risk score (HR 1.604, 95% CI 1.102–2.335, P=0.014) were variables that independently affect survival (Table 2, Figure 4B).

TABLE 2

Table 2. The univariate and multivariate Cox regression analysis results for the CGGA dataset.

Moreover, the univariate and multivariate Cox regression analyses based on the TCGA dataset showed that radiation_therapy, Chemo_status, and RiskGroup were potential clinical variables that independently affect survival (Table 2). The forest plot for multivariate Cox regression is shown in Figure 4B.

Protein-Level Verification of Genes in the Prognosis Model

The protein-level verification based on the HPA database showed that the upregulation and downregulation of proteins between normal and tumor samples were consistent with the expression of TNFRSF12A (Figure 5).

FIGURE 5

Figure 5. Immunohistochemical staining results of TNFRSF12A based on the HPA database. Immunohistochemical staining of TNFRSF12A in normal tissue [male, age 45; cerebral cortex (T-X2020), NOS (M-00100)] and in tumor tissue [male, age 47; brain (T-X2000) glioma, malignant, high grade (M-938033)].

Drug Prediction Analysis

A total of 13 gene–drug interaction, including three genes (TNFRSF12A, CD40, and TNFRSF11B) and 13 drug molecules (enavatuzumab, dacetuzumab, pg-102 inhibitor, teneliximab, tetanus toxoid, hydroquinone, streptozotocin, fludarabine, lucatumumab, zoledronic acid, risedronic acid, epinephrine, and testosterone), were explored based on the online database DGIdb (Figure 6A, Supplementary Table 3). Furthermore, the online database STITCH was further used to verify the relationship between drugs and corresponding proteins in genes (Figure 6B, Supplementary Table 4). The result showed that zoledronic acid–TNFRSF11B was the common interaction in both DGIdb database and STITCH database.

FIGURE 6

Figure 6. Drug- gene interaction network. (A) The drug–gene interaction network from the DGIdb database. Red ellipses represent genes, and rectangles represent drugs; the lines represent drug–gene interactions, of which colored lines represent known interactions while gray lines represent unknown interactions. (B) The drug–gene–gene interaction network from the STITCH database. The block represents the drug, while the circle represents gene. Red lines represent the drug–drug interaction, and green lines represent drug–gene interactions, and gray lines represent gene–gene interaction.

Discussion

GBM is the most common malignant brain tumor in adults (34). Although TNF and methylation sites are proved to be associated with the progression of disease (35), the detailed mechanisms are still unclear. In this study, the bioinformatics analysis showed that there were 41 DEGs between two groups. These DEGs were mainly enriched in GO functions, like response to tumor necrosis factor (GO: 0034612), and in KEGG pathways, like the hsa04668~TNF signaling pathway. Moreover, a total of 5 methylation site-regulated prognosis-related genes including TNFRSF12A, CD40, TNFRSF11B, TNFRSF10C, and LTBR were explored. The prognostic model constructed by the five genes was highly correlated with prognosis both in the TCGA cohort and in the CGGA cohort, and the higher risk score indicated lower survival. Finally, zoledronic acid–TNFRSF11B was revealed as key point drug–gene interaction by drug prediction analysis.

It is widely known that the response to TNF signaling plays an important role during peripheral organ inflammation in the brain (36). This signaling is widely proved to participate in the development of various diseases like ovarian cancer and lung cancer (37, 38). The response to TNF signaling mediates primary resistance to epidermal growth factor receptor inhibition in GBM (39). A previous study shows that TNF-α induces angiogenic factor upregulation in malignant glioma cells which play a role in RNA stabilization during GBM (40). Actually, the biological function of response to TNF is commonly realized by certain genes. TNFRSF12A is the sole signaling receptor for the proinflammatory cytokine TWEAK (TNFSF12) (41). Via TNF signaling such as TWEAK–TNFRSF12A, TNFRSF12A can regulate cellular activities including proliferation, migration, differentiation, apoptosis, angiogenesis, tissue remodeling, and inflammation (41). It has been proved that TNFRSF12A, which plays a role in tumor growth and metastasis, is highly expressed in solid tumor types (42). Yang et al. showed that the upregulation of TNFRSF12A contributed to poor prognosis in cancer (43). A previous study shows that differential expression of TNFRSF12A and DNA methylation contributes to the development of brain diseases such as epilepsy (44). Wang et al. proved that there was a close relationship between TNFRSF12A methylation and carcinoma prognosis (45). Based on the TCGA RNA-sequencing and methylation data, a previous study indicates that the methylation and expression levels of TNFRSF12A is significantly associated with prognosis of hepatocellular carcinoma, which can be used as a prognostic risk model (46). In the current study, TNFRSF12A was one of the five methylation site-regulated prognosis-related genes. Meanwhile, the enrichment analysis showed that TNFRSF12A was one of the genes that assembled in response to TNF function. Importantly, the protein-level verification analysis showed that the upregulation and downregulation of proteins between normal and tumor samples were consistent with the expression of TNFRSF12A. Thus, we speculated that methylation-driven gene TNFRSF12A might take part in the progression of GBM through response to the tumor necrosis factor biological process and TNF signaling pathway and significantly associated with the prognosis of GBM.

Zoledronic acid (ZA) is a potent inhibitor currently used in the clinical treatment of cancer (47). A previous study shows that ZA enhances T-lymphocyte antitumor response to human GBM cell lines (48). The antitumor effect of ZA combined with temozolomide can be used to against human GBM cell DNA methyltransferase (49). Actually, the drug effect of ZA is realized via stimulating the expression of certain genes in GBM cells (50). Karabulut et al. indicated that an induction in mRNA levels of TNFRSF family genes was observed in tumor cells under ZA treatment (51). Genetically achieved disturbances to the expression levels of TNFSF11B can modulate the effects of ZA on growing mouse skeletons (52). TNFRSF11B is a potential plasma biomarker for the clinical diagnosis of various cancers (53, 54). TNFRSF11B is proved to be differentially expressed in many immune cells in the brain (55). GBM is resistant to TNF-receptor family gene-induced apoptosis (56). A previous study indicates that the expression of TNF-receptor family genes including CD70 and TNFRSF11B was associated with the progression of GBM (57). In this study, the drug prediction analysis ZA–TNFRSF11B was revealed as the unique drug–gene interaction in both databases. Meanwhile, TNFRSF11B was revealed as one of the five GBM candidate genes in the methylation site-regulated prognosis-related gene analysis. Thus, we speculated that ZA targeting TNFRSF11B expression might be a potentially effective drug for GBM clinical treatment. However, there were some limitations in this study. Firstly, we preliminarily constructed a five-gene prognosis model. The prognostic performance of the five-gene prognosis model should be further confirmed by clinical samples. Secondly, we identified several potential drug targets. The drug–gene interactions should be validated by experiments, and the clinical value should be further investigated.

In conclusion, the methylation-driven gene TNFRSF12A might take part in the progression of GBM through response to the tumor necrosis factor biological process and TNF signaling pathway and significantly be associated with the prognosis of GBM. Moreover, ZA targeting TNFRSF11B expression might be a potentially effective drug for GBM clinical treatment.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

W-cL conceived and designed this research. HX carried out the plan and wrote this paper. CY, J-jL, and Z-yL gave advice and carried out the data analysis. All authors read and approved the final manuscript.

Funding

This work was supported by grants from the Guidance Plan of the Natural Science Foundation of Liaoning Province (Grant Nos. 201602773 and 2019-ZD-0340).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank all members in our lab for the excellent technical help.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2021.576382/full#supplementary-material

Supplementary Figure 1. The Venn diagram for integrate analysis based on prognostic genes, prognostic methylation, and genes negatively related to methylation.

Supplementary Figure 2. Prognostic verification analysis for current prognostic model based on tumor samples in CGGA database. (A) Survival analysis for high-risk group and low-risk group; survival time of high-risk group was shorter than that of low-risk group; The X-axis represented the overall survival time (month), while the Y-axis represented the survival rate (percent survival). P < 0.05 was considered to be significant different. (B) The risk score and follow-up in high-risk group and low-risk group. (C) The heatmap for methylation site regulated prognosis-related genes including CD40, LTBR, TNFRSF10C, TNFRSF11B, and TNFRSF12A.

Supplementary Table 1. The detailed information of the 76 TNF superfamily-related genes.

Supplementary Table 2. The 74 methylation sites were correlated with prognosis.

Supplementary Table 3. The obtained drug-gene interactions from DGIdb database.

Supplementary Table 4. Validation of interactions among the obtained drugs and genes in STITCH database.

References

1. Ostrom QT, Gittleman H, Farah P, Ondracek A, Chen Y, Wolinsky Y, et al. CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United states in 2006-2010. Neuro Oncol. (2012) 14 (Suppl. 5):v1. doi: 10.1093/neuonc/not151

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Young RM, Jamshidi A, Davis G, Sherman JH. Current trends in the surgical management and treatment of adult glioblastoma. Ann Transl Med. (2015) 3:121. doi: 10.3978/j.issn.2305-5839.2015.05.10

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Nam JY, de Groot JF. Treatment of glioblastoma. J Oncol Pract. (2017) 13:629. doi: 10.1200/JOP.2017.025536

CrossRef Full Text | Google Scholar

4. Gallego O. Nonsurgical treatment of recurrent glioblastoma. Curr Oncol. (2015) 22:e273. doi: 10.3747/co.22.2436

CrossRef Full Text | Google Scholar

5. Mao H, LeBrun DG, Yang J, Zhu VF, Li M. Deregulated signaling pathways in glioblastoma multiforme: molecular mechanisms and therapeutic targets. Cancer Invest. (2012) 30:48–56. doi: 10.3109/07357907.2011.630050

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Xu X, Stockhammer F, Schmitt M. Cellular-based immunotherapies for patients with glioblastoma multiforme. Clin Dev Immunol. (2012) 2012:764213. doi: 10.1155/2012/764213

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Collette Y, Gilles A, Pontarotti P, Olive D. A co-evolution perspective of the TNFSF and TNFRSF families in the immune system. Trends Immunol. (2003) 24:387–94. doi: 10.1016/S1471-4906(03)00166-2

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Badr CE, Wurdinger T, Nilsson J, Niers JM, Whalen M, Degterev A, et al. Lanatoside C sensitizes glioblastoma cells to tumor necrosis factor–related apoptosis-inducing ligand and induces an alternative cell death pathway. Neuro Oncol. (2011) 13:1213–24. doi: 10.1093/neuonc/nor067

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kore RA, Abraham EC. Inflammatory cytokines, interleukin-1 beta and tumor necrosis factor-alpha, upregulated in glioblastoma multiforme, raise the levels of CRYAB in exosomes secreted by U373 glioma cells. Biochem Biophys Res Commun. (2014) 453:326–31. doi: 10.1016/j.bbrc.2014.09.068

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Paulino VM, Yang Z, Kloss J, Ennis MJ, Armstrong BA, Loftus JC, et al. TROY (TNFRSF19) is overexpressed in advanced glial tumors and promotes glioblastoma cell invasion via Pyk2-Rac1 signaling. Mol Cancer Res. (2010) 8:1558–67. doi: 10.1158/1541-7786.MCR-10-0334

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Gkountela S, Castro-Giner F, Szczerba BM, Vetter M, Landin J, Scherrer R, et al. Circulating tumor cell clustering shapes DNA methylation to enable metastasis seeding. Cell. (2019) 176:98–112.e14. doi: 10.1016/j.cell.2018.11.046

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ahuja N, Sharma AR, Baylin SB. Epigenetic therapeutics: a new weapon in the war against cancer. Ann Rev Med. (2016) 67:73–89. doi: 10.1146/annurev-med-111314-035900

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gursoy-Yuzugullu O, Carman C, Serafim RB, Myronakis M, Valente V, Price BD. Epigenetic therapy with inhibitors of histone methylation suppresses DNA damage signaling and increases glioma cell radiosensitivity. Oncotarget. (2017) 8:24518–32. doi: 10.18632/oncotarget.15543

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Laffaire J, Everhard S, Idbaih A, Crinière E, Marie Y, de Reyniès A, et al. Methylation profiling identifies 2 groups of gliomas according to their tumorigenesis. Neuro Oncol. (2010) 13:84–98. doi: 10.1093/neuonc/noq110

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Richard AC, Peters JE, Lee JC, Vahedi G, Schäffer AA, Siegel RM, et al. Targeted genomic analysis reveals widespread autoimmune disease association with regulatory variants in the TNF superfamily cytokine signalling network. Genome Med. (2016) 8:76. doi: 10.1186/s13073-016-0329-5

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Goldman M, Craft B, Brooks A, Zhu J, Haussler D. The UCSC xena platform for cancer genomics data visualization and interpretation. BioRxiv. (2018) 326470 doi: 10.1101/326470

CrossRef Full Text | Google Scholar

17. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, et al. GENCODE: the reference human genome annotation for the ENCODE Project. Genome Res. (2012) 22:1760–74. doi: 10.1101/gr.135350.111

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. (2016) 44:e71. doi: 10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Quadrianto N, Buntine WL. Linear regression. In: Sammut C, Webb GI, editors. Encyclopedia of Machine Learning. Boston, MA: Springer (2010). p. 603–6.

Google Scholar

20. Koop G. Bayesian methods for empirical macroeconomics with big data. Rev. Econ. Anal. (2017) 9:33–56.

Google Scholar

21. Smyth GK. Limma: linear models for microarray data. In: Gentleman R, Carey VJ, Huber W, Irizarry RA, Dudoit S, editors. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. New York, NY: Springer (2005). p. 397–420.

Google Scholar

22. Wickham H, Wickham MH. The ggplot Package (2007).

Google Scholar

23. Kolde R, Kolde MR. Package ‘pheatmap' R Package. (2015).

Google Scholar

24. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bland JM, Altman DG. Survival probabilities (the Kaplan-Meier method). BMJ. (1998) 317:1572. doi: 10.1136/bmj.317.7172.1572

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bao Z-S, Chen H-M, Yang M-Y, Zhang C-B, Yu K, Ye W-L, et al. RNA-seq of 272 gliomas revealed a novel, recurrent PTPRZ1-MET fusion transcript in secondary glioblastomas. Genome Res. (2014) 24:1765–73. doi: 10.1101/gr.165126.113

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhao Z, Meng F, Wang W, Wang Z, Zhang C, Jiang T. Comprehensive RNA-seq transcriptomic profiling in the malignant progression of gliomas. Sci data. (2017) 4:170024. doi: 10.1038/sdata.2017.24

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Tissue-based map of the human proteome. Science. (2015) 347:1260419. doi: 10.1126/science.1260419

CrossRef Full Text | Google Scholar

32. Cotto KC, Wagner AH, Feng YY, Kiwala S, Coffman AC, Spies G, et al. DGIdb 3.0: a redesign and expansion of the drug–gene interaction database. Nucleic acids Res. (2017) 46:D1068–73. doi: 10.1093/nar/gkx1143

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kuhn M, von Mering C, Campillos M, Jensen LJ, Bork P. STITCH: interaction networks of chemicals and proteins. Nucleic Acids Res. (2007) 36:D684–8. doi: 10.1093/nar/gkm795

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ni Y, Zhang F, An M, Yin W, Gao Y. Early candidate biomarkers found from urine of glioblastoma multiforme rat before changes in MRI. bioRxiv. (2018) 1–6. doi: 10.1101/117333

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kan S, Chai S, Chen W, Yu B. DNA methylation profiling identifies potentially significant epigenetically-regulated genes in glioblastoma multiforme. Oncol Lett. (2019) 18:1679–88. doi: 10.3892/ol.2019.10512

PubMed Abstract | CrossRef Full Text | Google Scholar

36. D'Mello C, Le T, Swain MG. Cerebral microglia recruit monocytes into the brain in response to tumor necrosis factor signaling during peripheral organ inflammation. J Neurosci. (2009) 29:2089–102. doi: 10.1523/JNEUROSCI.3567-08.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Li BX, Wang HB, Qiu MZ, Luo QY, Yi HJ, Yan XL, et al. Novel smac mimetic APG-1387 elicits ovarian cancer cell killing through TNF-alpha, ripoptosome and autophagy mediated cell death pathway. J Exp Clin Cancer Res. (2018) 37:53. doi: 10.1186/s13046-018-0703-9

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Gong K, Guo G, Gerber DE, Gao B, Peyton M, Huang C, et al. TNF-driven adaptive response mediates resistance to EGFR inhibition in lung cancer. J Clin Invest. (2018) 128:2500–18. doi: 10.1172/JCI96148

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Guo G, Gong K, Ali S, Ali N, Shallwani S, Hatanpaa KJ, et al. A TNF–JNK–Axl–ERK signaling axis mediates primary resistance to EGFR inhibition in glioblastoma. Nat Neurosci. (2017) 20:1074. doi: 10.1038/nn.4584

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Nabors LB, Suswam E, Huang Y, Yang X, Johnson MJ, King PH. Tumor necrosis factor α induces angiogenic factor up-regulation in malignant glioma cells a role for RNA stabilization and HuR. Cancer Res. (2003) 63:4181–7. doi: 10.1002/cncr.22622

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Bhattacharjee M, Raju R, Radhakrishnan A, Nanjappa V, Muthusamy B, Singh K, et al. A bioinformatics resource for TWEAK-Fn14 signaling pathway. J Signal Transduct. (2012) 2012:376470. doi: 10.1155/2012/376470

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Cheng E, Whitsett TG, Tran NL, Winkles JA. The TWEAK receptor Fn14 is an Src-inducible protein and a positive regulator of Src-driven cell invasion. Mol Cancer Res. (2015) 13:575–83. doi: 10.1158/1541-7786.MCR-14-0411

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yang J, Min KW, Kim DH, Son BK, Moon KM, Wi YC, et al. High TNFRSF12A level associated with MMP-9 overexpression is linked to poor prognosis in breast cancer: gene set enrichment analysis and validation in large-scale cohorts. PloS ONE. (2018) 13:e0202113. doi: 10.1371/journal.pone.0202113

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Liu X, Ou S, Xu T, Liu S, Yuan J, Huang H, et al. New differentially expressed genes and differential DNA methylation underlying refractory epilepsy. Oncotarget. (2016) 7:87402–16. doi: 10.18632/oncotarget.13642

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang Y, Zhang S, Xie X, Chen Z, Wu L, Yu Z, et al. Association of TNFRSF12A methylation with prognosis in hepatocellular carcinoma with history of alcohol consumption. Front Genet. (2019) 10:1299. doi: 10.3389/fgene.2019.01299

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Li J, Chen N, Gong X. Prognostic implications of aberrantly expressed methylation-driven genes in hepatocellular carcinoma: a study based on the cancer genome atlas. Mol Med Rep. (2019) 20:5304–14. doi: 10.3892/mmr.2019.10771

CrossRef Full Text | Google Scholar

47. Pietrovito L, Comito G, Parri M, Giannoni E, Chiarugi P, Taddei ML. Zoledronic acid inhibits the RhoA-mediated amoeboid motility of prostate cancer cells. Curr Cancer Drug Targets. (2019) 19:807–16. doi: 10.2174/1568009619666190115142858

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Di Mascolo D, Varesano S, Benelli R, Mollica H, Salis A, Zocchi MR, et al. Nanoformulated zoledronic acid boosts the Vδ2 T cell immunotherapeutic potential in colorectal cancer. Cancers. (2020) 12:104. doi: 10.3390/cancers12010104

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Fukai J, Koizumi F, Nakao N. Enhanced anti-tumor effect of zoledronic acid combined with temozolomide against human malignant glioma cell expressing O6-methylguanine DNA methyltransferase. PloS ONE. (2014) 9:e104538. doi: 10.1371/journal.pone.0104538

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Avci CB, Kurt CC, Tepedelen BE, Ozalp O, Goker B, Mutlu Z, et al. Zoledronic acid induces apoptosis via stimulating the expressions of ERN1, TLR2, and IRF5 genes in glioma cells. Tumor Biol. (2016) 37:6673–9. doi: 10.1007/s13277-015-4519-3

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Karabulut B, Karaca B, Atmaca H, Kisim A, Uzunoglu S, Sezgin C, et al. Regulation of apoptosis-related molecules by synergistic combination of all-trans retinoic acid and zoledronic acid in hormone-refractory prostate cancer cell lines. Mol Biol Rep. (2011) 38:249–59. doi: 10.1007/s11033-010-0102-6

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Vargas-Franco JW, Castaneda B, Gama A, Mueller CG, Heymann D, Rédini F, et al. Genetically-achieved disturbances to the expression levels of TNFSF11 receptors modulate the effects of zoledronic acid on growing mouse skeletons. Biochem Pharmacol. (2019) 168:133–48. doi: 10.1016/j.bcp.2019.06.027

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wang X, Liu Y, Shao D, Qian Z, Dong Z, Sun Y, et al. Recurrent amplification of MYC and TNFRSF11B in 8q24 is associated with poor survival in patients with gastric cancer. Gastric Cancer. (2016) 19:116–27. doi: 10.1007/s10120-015-0467-2

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Luo P, Lu G, Fan LL, Zhong X, Yang H, Xie R, et al. Dysregulation of TMPRSS3 and TNFRSF11B correlates with tumorigenesis and poor prognosis in patients with breast cancer. Oncol Rep. (2017) 37:2057–62. doi: 10.3892/or.2017.5449

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Çogaş P. Evaluation of TNFRSF11B gene polymorphism in patients with acute stroke. J Clin Exp Invest. (2016) 7:184–9. doi: 10.5799/ahinjs.01.2016.02.0594

CrossRef Full Text | Google Scholar

56. Dixit D, Sharma V, Ghosh S, Mehta V, Sen E. Inhibition of casein kinase-2 induces p53-dependent cell cycle arrest and sensitizes glioblastoma cells to tumor necrosis factor (TNF α)-induced apoptosis through SIRT1 inhibition. Cell Death Dis. (2012) 3:e271. doi: 10.1038/cddis.2012.10

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Kenig S, FrangeŽ R, Pucer A, Lah T. Inhibition of cathepsin L lowers the apoptotic threshold of glioblastoma cells by up-regulating p53 and transcription of caspases 3 and 7. Apoptosis. (2011) 16:671–82. doi: 10.1007/s10495-011-0600-6

PubMed Abstract | CrossRef Full Text | Google Scholar