Identified RP2 as a prognostic biomarker for glioma, facilitating glioma pathogenesis mainly via regulating tumor immunity

Data sourcing and processing

Glioma mRNA expression data and clinical samples was collected from the TCGA Database (https://portal.gdc.cancer.gov/) and GTEx Database (https://commonfund.nih.gov) [8]. This study comprised 1158 normal samples and 706 GBMLGG samples for the gene expression profile. Clinical data from 599 GBM patients and 516 LGG patients were collected.

Cell culture

Glioma U251 cell line was purchased from the Chinese Academy of Science (Shanghai, China). U251 cells were cultured in RPMI-1640 medium, and other cell lines were cultured in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin. All cells were incubated in a 37°C incubator with 5% CO2.

Weighted gene co-expression network analysis

WGCNA is an R package that can be used to find clusters of highly correlated genes [9], and the optimal soft threshold β was determined in order to achieve the condition of a scale-free network. In addition, modules were identified using a dynamic tree cutting method. Correlations between Eigengenes and clinical traits were analyzed to determine the modules that were significantly correlated with clinical traits. Gene connectivity was measured by the absolute value of Pearson correlation. Genes with high within-module connectivity were identified as the central genes of the module.

TIMER database analysis

The TIMER (https://cistrome.shinyapps.io/timer/) is a collection of tools for systematically analyzing immune infiltrations in various cancer types [10]. The “Gene” and “scna” modules were used to predict the correlation among the infiltration of immune cells in each tumor sample and RP2 mRNA expression and copy number variation data in GBMLGG. And analyse the relationship between immune cell gene markers and other important contents with the help of the “Correlation” module. We also investigated the relationship between RP2 in pan-cancer and different immune cell infiltration.

SMART

SMART (Shiny Methylation Analysis Resource Tool) (http://www.bioinfo-zs.com/smartapp), an easy-to-use website based on data from TCGA, focuses on the correlation analysis of DNA methylation [11]. We used it to achieve CpG visualization of RP2 on the chromosomes.

MethSurv

MethSurv (https://biit.cs.ut.ee/methsurv/) is a web application that uses DNA methylation data to do variable survival analysis. The prognostic value of hypomethylation at some CpG sites in LGG and GBM were explored.

LinkedOmics

LinkedOmics (http://www.linkedomics.org/login.php), an online analysis site on the basis of TCGA tumor samples, designed to analyze multidimensional data of 32 kinds of tumors [12]. Using the “LinkFinder module”, the co-expressed genes linked to RP2 in the TCGA-GBMLGG database were visualized with volcano plots and heat maps. In addition, the effect of RP2 mRNA high expression on prognosis of glioma patients was analyzed in the same module. What’s more, we made Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis in “LinkInterpreter module”, to seek functional enrichment of co-expressed genes mentioned above. Pathways with FDR <0.05 were considered as meaningful.

GlioVis

The GlioVis (http://gliovis.bioinfo.cnio.es/) is an online platform used by us for data visualization and analysis to explore a wide range of information on the clinicopathological features and prognosis of gliomas [13].

Gene set enrichment analysis

GESA is a method for interpreting gene expression data, embodied in a software package together with a database of 1,325 gene sets. Relevant RNA-seq data of GBMLGG were downloaded from Genomic Data Commons (https://portal.gdc.cancer.gov/) [14]. The samples were divided into two groups according to the expression of RP2, and using this method we analyzed the enrichment of gene set in which RP2 is located. The parameters were established such as: gene set database: h. All. V7.4 Symbols. gmt (Hallmarks); number of permutations: 1,000. P value < 0.05 and FDR < 0.05 were considered as meaningful.

PPI network construction

STRING (https://string-db.org/) is usually adopted for exploring physical interactions or functional associations between proteins. With the help of it, we built a Protein-Protein Interaction Network (PPI) with the top 500 genes most closely related to RP2 from volcano plots. We studied the connection among these genes. The parameter of medium confidence was set at 0.9. Then, the top 500 genes were evaluated by Cytoscape 3.9.1, and a functional cluster analysis were performed with MCODE. The selection criteria are as follows: Max depth = 100, node score cutoff = 0.2, K-core = 2.

GeneMANIA analysis

GeneMANIA (http://www.genemania.org) [15], is a truly powerful online resource for gene function and lists analyses. We used it to draw an interactive functional network of RP2, mainly to find proteins with strong physical interaction with RP2. In the network, we used lines in various thickness and colors to show the functional relationship and correlation strength between the two connected ends.

cBioPortal analysis

The cBioPortal for Cancer Genomics (https://www.cbioportal.org/) is a resource for discovering, visualizing, and analyzing multidimensional genomics data [16]. “Mutations” module was utilized to found out RP2 and its physical interaction protein in glioma.

Protein structure and docking analysis

RCSB Protein Data Bank (https://www.rcsb.org/) enables breakthroughs in exploration of 3D protein structure. The structure of Retinitis Pigmentosa 2 (RP2) was obtained from PDB ID: 2BX6; the structure of ARL3 was obtained from PDB ID: 4GOJ; the structure of WDR83 was obtained from PDB ID: Q9BRX9. The binding patterns between RP2 and ARL3 and between RP2 and WDR83 were studied by docking through the Z DOCK web server (https://zlab.umassmed.edu/zdock/index.shtml) and visualized, analyzed and plotted using PYMOL.

TISIDB

TISIDB (http://cis.Hku.hk/TISIDB/) offers a wealth of information about tumor immunity and makes it possible to do thorough research on tumor-immune interactions [17]. There, we looked for a relationship between RP2 and six aspects of the immune system in gliomas (lymphocytes, immunomodulators, chemokines, etc.).

GSCA analysis

GSCA (http://bioinfo.life.hust.edu.cn/GSCA) is an integrated platform for genomic, pharmacogenomics and immunogenomic genomic gene set cancer analysis [18]. Combining clinical information and small molecule drugs, users can mine candidate biomarkers and valuable drugs for better experimental design and further clinical trials. We used it to uncover the drug sensitivity of ARL3, SLC44A2, WDR83, RP2 and OSTF1.

CTD

The Comparative Toxicogenomics Database (CTD, http://ctdbase.org/) is a digital resource that facilitates the study of novel links in the molecular mechanisms by which chemical substances affect health outcomes. We have used this database to query and visualize the interacting drugs or small molecules of RP2.

Quantitative RT-PCR analysis

Total RNAs were extracted using the QIAGEN RNeasy mini kit, and reverse transcription reactions were performed using the ABI Taqman Reverse Transcription Reagents. After mixing the generated cDNA templates with primers/probes and ABI Taqman Fast Universal PCR Master Mix, reactions were performed with the ABI-7900 Fast Real-time PCR system and SYBR green qPCR Mastermix from Agilent Technologies Stratagene. This technique was adopted in our study to examine the relative mRNA expression of TGF-β and IL-10.

Statistical analysis

R software (version 4.2.1) was used to conduct all statistical analysis in this research. The different expression of RP2 was detected by rank-sum test, meanwhile, the “limma” and “beeswarm” packages were driven. Multivariate Cox regression analysis screened factors significantly related to prognosis (p < 0.05) (Cox model uses the “survival” and “survminer” packages of “R”). SurvivalROC (version 1.0.3) was used to produce the ROC curve for evaluating the ability of RP2 expression level to predict 1-year, 3-year, or 5-year survival. We used the “Survival” and “SurvMiner” packages to produce the survival scene of RP2 in pan cancer. The R software package MAfTools (version 2.8.05) was used to calculate the TMB (Tumor mutation burden) of each Tumor by TMB function. The correlation between RP2 expression and immune checkpoint-related genes using R software (“pheatmap” packages). The particular status of the tumor mutation load was shown by log transformation after we combined the TMB and gene expression data of samples.

Data availability statement

Publicly available datasets were analyzed in this study. The data are accessible in TCGA and GTEx databases. Further inquiries can be directed to the corresponding author.

Weighted gene coexpression network analysis was used to identify the key gene in glioma

To identify pivotal genes in glioma, we first discovered that 8762 genes were upregulated, while 6964 were downregulated in glioma tissue (Figure 1A). Following that, β = 12 was selected as the soft threshold for implementing a scale-free network (Figure 1B). Using the dynamic tree cutting package, 7 modules were defined (Figure 1C). It should be noted that the correlation heatmap for module traits revealed that the turquoise module had the highest association with glioma (Figure 1D). The gene distribution results in the turquoise module displayed that glioma and Module membership (MM) were highly associated, suggested that genes in this module were highly significantly correlated with glioma (Figure 1E). Furthermore, Venn diagram showing 10 intersecting genes, which were most associated with glioma, screened by the WGCNA method combined with Multivariate Analysis, ROC and survival curves (Figure 1F). What’s more, RP2 expression was significantly higher in GBMLGG, founded in the TIMER database (Figure 1G). Finally, by inquiring background information, we selected RP2, which has not been reported in glioma and has the best results, as our research target.

Figure 1. Identification of the key gene modules in WGCNA. (A) Volcano map showed differentially expressed genes. (B) Determination of the soft-thresholding power. (C) Dendrogram of differentially expressed genes clustered based on a dissimilarity measure (1-TOM). (D) The correlation of gene modules with clinical traits. (E) Gene correlation scatter plot of the turquoise module. (F) The Venn diagram showed WGCNA combined genetic features to screen out 10 genes. (G) The expression of RP2 in 34 kinds of cancers.

Association of RP2 expression and clinicopathological characteristics in glioma

In order to determine the relevance of RP2 expression to clinicopathological characteristics, we used R software to analyse. We noticed that RP2 mRNA expression was upregulated in patients above 60, with higher tumour grade, unmutated group of IDH and in the non-coding 1P/19Q (Figure 2A–2D). Futhermore, we discovered that the mRNA expression of RP2 was most pronounced in glioblastoma among all histological types (Figure 2E). These outcomes conferred that the RP2 expression was intimately associated with the clinicopathological characteristics.

Figure 2. The box plot showed the association of RP2 expression with clinicopathological characteristics. (A) Age, (B) grade, (C) IDH status, (D) 1p/19q codeletion, (E) histological type.

Over-expression of RP2 was associated with awful prognosis of glioma patients

To further find the underlying mechanism of RP2 over-expression in glioma patients, we analyzed the effects of RP2 over-expression in TCGA databases on overall survival (OS), disease-free survival (DSS), and progression-free survival (PFS) of glioma patients and plotted survival curves. The findings showed that patients with high RP2 expression had proportionately lower OS, DSS, and PFS times than the low expression group, demonstrating that high RP2 expression was correlated with worse prognosis (Figure 3A–3C). Overall survival was also examined using LinkedOmics to further demonstrate the link between RP2 expression and glioma patient prognosis (Supplementary Figure 1A), and the high and low RP2 expression groups were differentiated by the GlioVis in another way, and both results indicated that the high RP2 expression group had a worse survival expectation than the low RP2 expression group (Supplementary Figure 1B–1D). Besides, in order to forecast the overall survival of glioma patients in the TCGA cohort, we developed a nomogram. The WHO grade for the malignancy, IDH status, 1p/19q codeletion, and RP2 were included in the nomogram as prognostic variables. (Figure 3D). And the calibration curve displayed that the nomogram was credible in predicting possibility of 1-, 3-, 5-years overall survival in glioma (Figure 3E). In conclusion, our data suggested that a worse prognosis for glioma patients is frequently linked to increased RP2 expression.

Figure 3. Relationship between RP2 and prognosis of glioma patients. Glioma patients with lower expression level of RP2 had favorable (A) OS (HR = 2.37, p < 0.001), (B) DSS (HR = 2.54, p < 0.001), (C) PFS (HR = 1.98, p < 0.001). (D) Nomograms and (E) calibration curves.

Correlation between RP2 hypomethylation and prognosis in LGG and GBM

Although we were aware that RP2 was expressed significantly more frequently in gliomas and that this may result in a bad prognosis [19], we didn’t know what was causing it. DNA methylation affects glioma development by altering the expression of key genes [20]. Therefore, we will investigate RP2 related DNA methylation. In LGG, we found a negative correlation of methylation values of cg24511534, cg04586456 and other 7 methylation probes and RP2 expression levels (Figure 4A). In GBM, there was a negative association between the methylation values of the methylation probes cg00433220, cg00347850, and other 10 methylation probes and the expression levels of RP2 (Figure 4B). We further explored how the hypomethylation of these sites would affect the prognosis of glioma with MethSurv, the result of which showed that patients with lower RP2 methylation had worse OS than those with higher RP2 methylation (P < 0.05) (Figure 4C, 4D). Besides, we also visualized the distribution landscape of methylation probes for RP2 on the chromosome, with cg04586456, cg25366157 and other 8 methylation probes on island, cg22522912 and cg14091713 on N_ Shore, and cg00433220 on S_Shore (Figure 4E, 4F).

Figure 4. Correlation between RP2 hypomethylation, and prognosis in GBM and LGG. (A, B) The visualization between the methylation level and the RP2 expression. (C, D) The Kaplan–Meier survival of the promoter methylation of RP2. (E, F) The distribution landscape of methylation probes.

In summary, there was a negative association between the expression of RP2 and its methylation level. It can be concluded that High RP2 expression in gliomas is caused by diminished RP2 methylation, which ultimately results in a bad prognosis for glioma patients.

RP2-associated functional enrichment pathway

Since RP2 presents a clear value in terms of patient prognosis, we produced a landscape of RP2-related genes to further investigate its functional regulatory pathways. The results came out that among the 20119 correlated genes identified, 7081 genes were positively associated (red dots) and 6743 genes were negatively associated (dark green dots) with RP2 expression (FDR <0.01) (Figure 5A). Therefore, we individually took out the top 50 positively and negatively linked genes to make a correlation heat map (Figure 5B, 5C). Gene ontology (GO) showed that RP2 may be involved in processes related to T cell activation, cell cycle, G1/S phase transition, and interferon production (Figure 5D). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis displayed that RP2 was positively associated with processes such as antigen processing presentation and cell cycle (Figure 5E). Additionally, GSEA displayed that RP2 was closely related to protein export, proteasome, and RNA degradation, and so on (Figure 6A–6E). These findings suggested that RP2 may be engaged in processes related to cell cycle and organismal immunity. To further understand the mechanism of RP2 regulatory function, we selected genes with RP2 correlation located in the top 500 to produce clusters of RP2-related protein interactions. The clusters of proteins most associated with RP2 are shown by Figure 6F (yellow). We selected them to visualize them individually (Figure 6G), and the results showed that RP2 had protein interactions with ITGA4 and CDC42, among others, and verified that these proteins show correlation with immune response and metastatic spread of tumor cells, confirming our speculation about the function of RP2. We speculated that RP2 is associated with the cell cycle and immune infiltration.

Figure 5. Potential biological processes that RP2 involved in. (A) Differential gene expression in GBMLGG. (B, C) The heat map respectively displayed the top 50 positively correlated and negatively correlated genes of RP2 in GBMLGG. (D, E) GO and KEGG enrichment analysis in GBMLGG of RP2.

Figure 6. RP2-associated (A–E) pathways and (F) protein interaction networks as well as (G) the most correlated clusters.

Structural and physical interactions of proteins

Proteins often need to dock several tertiary structures to each other to form quaternary structures for biological functions and metabolic reactions [21]. To investigate the potential mechanisms of protein physical interactions, we first selected RP2 to construct a network of proteins with which it can interact (Supplementary Figure 2A). The results showed that RP2 and proteins such as ARL3 and WDR83 can form strong physical interactions. The secondary structures of RP2, ARL3, and WDR83, as well as their possible chemical modifications, such as phosphorylation and ubiquitination were visualized (Supplementary Figure 2B–2D). Besides, the peptide chain encoded by RP2 coiled and folded with ARL3 and WDR83 form a quaternary structure as shown in Supplementary Figure 2E, 2F. We speculate that these two structures may be a form of RP2 involved in the regulation of cell cycle and immunity.

RP2 expression is correlated with immune infiltration in glioma tissues

To investigate the part played by RP2 in influencing tumor immune microenvironment during tumor development, the Immune score (P = 2.5e-18, r = 0.33) and Estimate score (P = 3.5E-19, r = 0.34) of glioma samples have been analyzed. Both showed a clear correlation with RP2 expression (Figure 7A, 7B). Infiltration of six different types of immune cells were shown to be correlated with RP2 expression in glioma. The findings displayed that RP2 expression was definitely correlated with tumor purity in LGG and GBM, and showed an up-regulated relationship. In addition, we analyzed correlations between RP2 expression and multiple gene markers of immune cells in GBM and LGG (Tables 1, 2). As a result, we found RP2 was positively associated with the infiltration of Neutrophils, CD4+T cells, B cells, Dendritic cells, and Macrophages in glioma (Figure 7C, 7D). Since copy number variation (CNV) affects the degree of immune cell infiltration in gliomas [22], we looked into how copy number variation affects immune cell infiltration in LGG and GBM. Our results showed that high amplification notably promoted CD4+ T cell infiltration in LGG, however, in GBM, deep deletion, arm-level deletion and arm-level gain all reduced CD8+ T cell, neutrophil and dendritic cell infiltration to varying degrees (Figure 7E, 7F). In conclusion, our results showed an association between RP2 and tumor cell infiltration in glioma.

Figure 7. RP2 is associated with tumor immune microenvironment in LGG and GBM. (A, B) Correlation between RP2 and immune score and estimated score. (C, D) Correlation between tumor infiltrating immune cells and RP2 expression. (E, F) Correlation between CNV of RP2 and the degree of immune cell infiltration in LGG and GBM.

Infiltration of some immune subtypes under the condition of high RP2 expression

With the purpose of further exploring which immune cells the high expression of RP2 had the highest correlation with and the expression of RP2 in the TME, we performed the following analysis. The TIMER2.0 database was utilized to display the landscapes in which RP2 linked with diverse immune cell infiltrations in malignancies done on several quantitative immune infiltration platforms. The findings reveal that RP2 is associated with the immune infiltration of a wide spectrum of infiltrating cells, including CD4+ T cells, CD8+ T cells, Myeloid dendritic cells, Monocyte, Macrophage, and Neutrophil. Similarly, we found that in LGG and GBM, RP2 is negatively correlated with CD4+ T cells and positively correlated with Monocytes and Macrophages when highly expressed (Figure 8A). After that, to further investigate RP2 expression in the tumor microenvironment of gliomas, we carried out single-cell level analysis. We analysed seven single cell sequencing datasets, including glioma_GSE103224, through the TISCH website, and showed that RP2 is highly expressed in immune cells, particularly in M1 macrophages, M2 Macrophages, and Monocyte cells (Figure 8B, 8C). In conclusion, we discovered that immunosuppressive cell infiltration was positively connected with high RP2 expression in gliomas and was negatively correlated with immune helper cell infiltration. Furthermore, RP2 was significantly and relatively highly expressed in immunosuppressive cells, like M2 Macrophages, in TME.

Figure 8. The correlations of RP2 expression and the infiltration levels of immune cells. (A) The correlations of RP2 expression and the infiltration levels of CD4+ T cells, CD8+ T cells, Myeloid dendritic cells, Monocyte, Macrophage and Neutrophil in cancers. Positive correlation in red and negative correlation in purple. (B, C) Single cell sequencing showed the expression of RP2 in different types of cells.

The dynamics of macrophage cells during glioma progression

To describe the landscape of various cellular presence in the tumor microenvironment, we selected single-gene transcriptomes for clustering analysis, and after quality control to filter out irrelevant cells, Linear dimensionality reduction were used and visualization by t-distributed stochastic neighbor embedding (t-SNE) method. Cells from tumor and normal tissues were divided into 15 clusters (Figure 9A, 9B), and then clustered into six main categories according to the genetic characteristics of different subclusters: astrocytes, macrophages, monocytes, natural killer cells, T cells, and B cells, with monocytes and macrophages accounting for a large proportion of the tumor microenvironment. Meanwhile, we found that RP2 was mainly expressed in monocytes and macrophages (Figure 9C–9F).

Figure 9. Macrophage dynamics in glioma progression. (A, B) Quality control removal of low-quality cells (C) tSNE cell cluster analysis differentiating cell types in the tumor microenvironment. (D) Tumor cells were divided into 13 clusters showing RP2 expression sites in the tumor microenvironment. (E) Cell reclustering was performed by UAMP method, and macrophage types were analyzed. (F) expression of RP2 was shown in the reclustering cells. (G) Evolutionary tracks of all kinds of cells after dimension reduction (H) expression of RP2 changes with pseudo time (I) differentiation tracks of macrophages in glioma, with cluster color code.

We analyzed macrophages by downscaling, unsupervised clustering and trajectory analysis, and found that they could be further identified as Intermediate monocytes and Myeloid dendritic cells (Figure 9G), and the expression of RP2 did not show significant differences among these two types of cells (Figure 9H). The dynamic cellular states in the tumor microenvironment were visualized by our Monocle algorithm, and pseudo-temporal trajectory analysis showed that Intermediate monocytes tend to be depleted in a trajectory with tumor progression, while RP2 would show low expression in the early stages of cell development and higher expression over time, while expression was gradually downregulated after the node (Figure 9I). The above results suggested that the degree of RP2 gene expression influences the progression of the glioma tumor microenvironment and it might play an important part in the evolution of tumor microenvironment in glioma patients.

Chemokines associated with glioma

Macrophages can secrete many cytokines, which play a key role in the tumor microenvironment and development of glioma. To investigate the key effect made by various cytokines produced by macrophages in TME, we visualized the contribution of different types of it in heat map form, and the results showed that: immune inhibitor factors such as TGFBR1 and IL10; MHC molecules such as B2M and HLA-A; immunostimulatory factors such as CD28 and CD48; lymphokines such as Act CD8 and Act CD4 had a positive correlation with RP2 in both GBM and LGG (Figure 10A–10D), Noteworthily, it is conformed that Macrophages may promote the progression of glioma through secrete cytokines like IL10 and TGF-β. We analyzed the correlation between RP2 and these two (Figure 10E, 10F), and found that they were positively correlated. In addition, we also verified them with experiments (Figure 10G, 10H). In conclusion, we speculated that macrophages may cause the poor prognosis of glioma by secreting these cytokines or influenced by them.

Figure 10. The expression of various cytokines when RP2 is highly expressed in LGG and GBM. (A) Immunoinhibitor. (B) MHC molecule. (C) Immunostimulator. (D) Lymphocyte. (E, F) RP2 was associated with the expression of TGF-β and IL-10. (G, H) RP2 was positively correlated with TGF-β and IL-10 mRNA expression.

Effect of over-expression RP2 on immunosuppression and immunotherapy of glioma

In order to further investigate the impact of RP2 upon this immune escape microenvironment and immunotherapy in glioma, we also analyzed the correlation between RP2 and immune checkpoint and the associated TMB, and TIDE. The findings demonstrated a favorable correlation between the expression of most cell checkpoints and the expression of RP2 in GBMLGG, demonstrating that RP2 overexpression might disturb the normal cell cycle and promote the proliferation of glioma (Figure 11A). We also analyzed the relationship between the expression of some classic immune checkpoints and RP2 in glioma, and found that the expression of immune checkpoints increased when RP2 was high, indicating enhanced immune escape (Supplementary Figure 3). The expression of RP2 in GBMLGG was positively related with TMB, suggesting that high expression of RP2 may promote the development of cancer by affecting gene mutations in glioma, meanwhile, a low mutation rate is not conducive to immunotherapy (Figure 11B). Finally, TIDE score was analyzed and it was found that the treatment of the high RP2 group was not sensitive to immune checkpoint blockade (ICB) and the therapeutic effect was poor (Figure 11C). In conclusion, we known that high RP2 expression was conducive to immune escape of glioma, but not conducive to immunotherapy.

Figure 11. RP2 and immunotherapy. (A) Correlation between RP2 and various cell checkpoints. (B) TMB was low when RP2 was highly expressed. (C) TIDE score was high when RP2 expression was high which led to a poor immunotherapy effect.

Gene-drug interactions and drug sensitivity analysis

To explore the drug sensitivity of RP2 and genes highly associated with it, we conducted drug sensitivity analysis using the GDSC database, which showed that RP2 exhibited tolerance to 19 drugs or small molecules when highly expressed. And it showed sensitivity to Trametinib (Figure 12A). To explore the interaction between the pivotal gene and the available therapeutic agents for cancer, we used the CTD database to produce the interaction network between RP2 and drugs, and the results showed that cobaltous chloride and jinfukang were able to inhibit the gene expression of RP2, while six drugs or small molecules such as Methotrexate and Cisplatin were able to promote the expression of RP2 (Figure 12B). We compared three databases, cgp2016, CTD and GDSC, and concluded that Methotrexate was present in all three databases to regulate RP2 at the same time (Figure 12C). Patients are not sensitive to this drug when RP2 is highly expressed and it has been shown to increase RP2 expression, so this drug would not be the preferred option in treatment.

Figure 12. Gene-drug interaction and drug sensitivity analysis. (A) Relationship between RP2 expression and drug sensitivity in the GSCA database. (B) CTD showed the network of action between RP2 and drugs or small molecules. (C) Drug susceptibility results from three databases were displayed.

The relationship between RP2 expression and m6A modification

In several cancer types, the effect of N6-methyladenosine (m6A) RNA modification has been proven [23]. By examining the TCGA-GBMLGG dataset, we explored the relationship among RP2 and the expression of 20 m6A-related genes in GBMLGG. RP2 expression was shown to be strongly and positively linked with 20 regulatory factors except FTO and IGF2BP2 (Supplementary Figure 4A, P < 0.05). We sought to find out if 20 distinct variables expressed differentially in glioma cells with high and low RP2 expression. The findings revealed that 20 variables expressed themselves more strongly in the “High RP2” group as compared to the “Low RP2” group (Supplementary Figure 4B). Meanwhile, based on RP2 expression levels, we separate samples into groups with high and low expression. M6A regulator networks depict a comprehensive picture of the correlation between the expression of m6A regulatory factors and the effect of regulatory factors on the prognosis of glioma. We found significant correlations in expression and prognosis not only between m6A regulatory factors of the same functional class, but also between writers, erasers and readers (Supplementary Figure 4C, 4D, Table 3). According to these findings, RP2 may be correlated with m6A modification in glioma, and the combined effect of METTL14, VIRMA, ALKBH5, YTHDC2 and METTL3 may eventually affect the progression and poor prognosis of glioma.