Target RNA modification for epigenetic drug repositioning in neuroblastoma: computational omics proximity between repurposing drug and disease

Current RNA-seq sample and CNV sample was suitable for information retrieval in NBL patients

We used a ‘NBL severity signature—target discovery—drug discovery—drug response biomarker’ process to predict repositioning drug (Figure 1A). We collected a list of 20 RNA-MRs regarding differentiation of neural progenitor cells and neurological diseases (Figure 1B), which had functional relationship (Figure 1C). The contributions of 20 RNA-MRs on NBL severity signature was studied as Figure 1D.

Figure 1. Summary for schematic diagram and RNA modification regulators (RNA-MRs) in NBL. (A) The schematic diagram of target mining and drug screening in this study. (B) A schematic view of location and types of predominant RNA modification. W: writer (methyltransferase), R: reader (RNA-binding protein), E: eraser (demethylase). TRDMT1: formerly DNMT2. (C) The protein interaction network among RNA-MRs using STRING analysis. KIAA1429: VIRMA. (D) The detailed schematic diagram for determination of NBL severity signature.

The age, International Neuroblastoma Staging System (INSS) stage, MYCN gene amplification, DNA ploidy, histology, MKI (mitosis karyorrhexis index), and Children Oncology Group (COG) risk score were commonly used as prognostic indexes in clinic for evaluating NBL severity. Three datasets were retrieved from Therapeutically Applicable Research to Generate Effective Treatments (TARGET) program and The Cancer Genome Atlas (TCGA) network: (1) the clinical dataset containing clinicopathological information of 834 NBL patients; (2) the RNA-seq dataset containing transcriptome data and clinicopathological information of 122 NBL patients, and (3) the CNV dataset containing copy number variation data (obtained by whole genome sequencing) and clinicopathological information of 78 NBL patients (Supplementary Figure 1D). To test the reproducibility of patients with RNA-seq/CNV data for overall NBL patients, clinicopathological information and their association with overall survival for 3 datasets is compared. Compared with total patients, RNA-seq dataset contained older patients with higher mortality, COG risks, unfavorable histology, and advanced INSS stage, while CNV dataset contained patients with more significantly amplification of MYCN gene, (Table 1, P<0.05). These indexes showed similar prognostic patterns among three groups (Supplementary Figure 1A–1C). Information retrieval from RNA-seq and CNV data was thus suitable for NBL patients.

Differenced expression of RNA-MRs defined two NBL clusters with different survival and clinicopathology

Twenty RNA-MRs had strong correlation at RNA expression level in NBL patients (Figure 2A). Thus, they shared similar functions in NBL disease course according to the ‘Guilt-by-Association’ principle [18]. The k-means clustering was a most commonly used centroid model in cancer study to detect subtypes, which was useful for prognosis prediction and personalized treatments [19]. As shown in Supplementary Figure 2A, RNA-seq samples could be clustered into 2~5 sub-clusters, with clustering stability increasing from k = 2 to 10. The performance of 2/3/4/5-cluster were quantified using adjusted Rand index (ARI), which evaluated the similarity between the resulting partition and the gold standard partition [20]. Here, two survival state (alive or dead) were assumed to be gold standard partition, and ARI were 0.136/0.076/0.086/0.059 for 2/3/4/5-cluster. Interestingly, the overall clustering performance for 2-cluster was high while the clustering stability was relatively low, compared with 3/4/5-cluster. This showed that a simple classification problem is not necessarily an easy clustering problem. Meanwhile, cluster1 in 2-cluster seemed to have a more significantly poor prognosis (P= 8.825×10^-6vs. 9.028×10^-4/1.485×10^-4/6.387×10^-4 in 3/4/5-cluster), and a clear distinct transcriptional landscape in 2-cluster could be seen using principal components analysis (PCA) (Figure 2B, 2C and Supplementary Figure 2B). Thus, we chose 2-cluster as representative classification for NBL severity. As shown in Figure 2D, cluster1 showed upregulated A-I editing in tRNA modification, m5C, hm5C, and m6A, together with downregulated A-I editing in mRNA modification (P<0.01). And cluster1 was also significantly correlated with higher COG risks, higher MKI, unfavorable histology, MYCN gene amplification and advanced INSS stage (P < 0.05, Table 2).

Figure 2. Differences of gene expression and CNV profile, survival and clinicopathological features in two groups defined by RNA modification regulators (RNA-MRs). (A) Spearman correlation analysis of RNA-MRs. Red: positive relation of gene expression between two genes; blue: negative relation. Color scale: the degree of relation. Cross (×): P>0.05 (no significant difference). (B) Kaplan–Meier overall survival curves for NBL patients in two groups. Ordinal: the percentage of survival; abscissa: survival years. "+" marks in line: censoring samples. (C) Principal component analysis of the transcriptome expression profile in the space of the first two principal components (Comp.). (D) Heatmap comparison of RNA-MRs for two groups. Horizontal color stripe above heatmap: one clinicopathological feature per line. Vertical color stripes: different methylation type regulated by each RNA-MR. Heatmap: expression differences of RNA-MRs (gradient color from green to red in each line showed downregulated to upregulated levels of each gene). *P < 0.05, **P < 0.01, and ***P < 0.001 between Cluster1 and Cluster2. (E) The OncoPrint of CNV pattern in two groups. Blue dot: deletion; green dot: gain; red dot: amplification. The upper barplot: the number of CNVs per patient; the right barplot: the number of genetic mutations per gene (corresponding rate was at left); the bottom color stripe: clinicopathological features for each patient.

These results showed that dysregulated expression of RNA-MRs strongly reflected the changes of total transcriptome, and the disease severity as well. Dysregulated RNA-MRs predicted poor patient prognosis, independently of other prognostic markers such as MYCN amplification status, age at the time of diagnosis, and disease stage.

CNVs of RNA-MRs were associated with their expression

CNV samples belonging to cluster1 was also significantly correlated with higher COG risks, higher MKI, unfavorable histology, and advanced INSS stage (P < 0.0001, Figure 2E and Table 2).

We observed high CNVs of RNA-MRs in genome level of NBL. In detail, the m5C “reader” ALYREF had the highest frequency of CNV events (85%), followed by the m6A “eraser” ALKBH5 (22%), both of which were located in chromosome 17 (Figure 2E and 3A). The amplification of ALYREF was the most frequent alteration in cluster1, and the simultaneous gain of ALYREF and ALKBH5 ranked first in cluster2, implying the importance of CNVs on dysregulated RNA modification network and NBL severity.

Figure 3. Summary of CNV region and correlation between CNV pattern and gene expression. (A) Circos plot illustrating CNV regions for RNA modification regulators (RNA-MRs). Outermost circle: chromosomes; inner circle: the main type and location of CNVs (red: gain, blue: deletion). (B) Correlation between CNV patterns and expression levels. Ordinal: the relative gene expression levels; abscissa: CNV types. * P < 0.05, ** P < 0.01, and *** P < 0.001 among different CNV types.

The effects of CNVs of RNA-MRs on their mRNA expression were further confirmed. Results showed that copy number gain or amplification were related to higher mRNA expression, while deletion resulted in a decline of mRNA expression (Figure 3B). CNVs in drug targets may change their abundance and drug-binding pocket of receptors [21], thus the high frequency of CNVs for RNA-MRs here made inhibitors for their enzymatic activity and antagonists for reader proteins more easily to be resistance.

Differences of functional annotation in two NBL groups

The above findings suggested that the clustering NBL according to gene expression of RNA-MRs was closely correlated to the malignancy of NBL. Then DEGs and functional differences between groups were evaluated. Altogether, 959 and 1154 genes were respectively upregulated and downregulated significantly (log2 (fold change) >1 or <-1, P < 0.05 and the false discovery rate (FDR) <0.05) in cluster1, which showed more severe disease course. In cluster1, DEGs enriching in biological process (BP) of GO (gene ontology) were predominantly associated with malignancy-related biological processes, including cell proliferation, extracellular matrix organization, angiogenesis, and migration (Figure 4A). The DEGs in GO cellular component (CC) were enriched in chromosome part, cytoplasmic membrane and extracellular matrix (Figure 4B). And DEGs in GO molecular function (MF) were associated with chromatin binding, catalytic activity on DNA, and ATPase activity (Figure 4C). Similar changes in corresponding signaling pathways were also observed in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Figure 4D).

Figure 4. GO, KEGG pathway enrichment analysis and the gene set enrichment analysis (GSEA) of differentially expressed genes (DEGs) in two groups. (A) Biological process of GO analysis; (B) Cellular component of GO analysis; (C) Molecular function of GO analysis; (D) KEGG pathways analysis. Y-axis: name of enrichment item; X-axis: gene count in each item (positive value for cluster1 and negative for cluster2); gradient color of the bar chart: -Log2_Pvalue (values >0 for -Log2_Pvalue in Cluster1, and values <0 for -(-Log2_Pvalue) in Cluster2). (E) GSEA enrichment plots in HALLMARK datasets from left to right: MYC targets, E2F targets, unfolded protein response, and MTORC1 signaling. In each plot, top panels: enrichment score (ES) for each gene; bottom panels: the ranking metrics of each gene; Y-axis: ranking metric values; X-axis: ranks for all genes. NES: normalized ES. Norm P (nominal p value): the statistical significance of the observed ES for this plot. FDR q value: a significantly enriched enrichment plot if q< 0.25.

Furthermore, gene set enrichment analysis (GSEA) revealed that the malignant hallmarks of tumors, including MYC targets (NES=2.13, nominal P < 0.001), E2F targets (NES=2.09, nominal P < 0.001), unfolded protein response (NES=1.92, nominal P = 0.002) and MTORC1 signaling (NES=1.64, nominal P=0.028), were significantly associated with cluster1 (Figure 4E).

All of these findings indicated that the dysregulated expression of RNA-MRs was closely correlated with the malignancy of NBL.

Mapping the NBL signature onto drug-perturbed transcriptome profiles in L1000FWD

To evaluate core prognostic gene set as a reduced representation of DEGs, the top 50 records with opposed expression profiles of DEGs and core prognostic gene set were collected (Supplementary Figure 3). These records had similarity scores ranging from 0.090 to 0.141, and from 0.303 to 0.483 (P<0.0001), respectively (Supplementary Tables 1 and 2). Among them 40 records for 32 drugs were overlapped for two queries, implying the representativeness of core prognostic gene set for DEGs and disease severity (Figure 5B). And 11 overlapped drugs were previously reported in NBL cells or clinic. Targets for top 50records were mainly PI3K, topoisomerase, mTOR, estrogen receptor, and CDK. Besides that, reversed expression of RNA-MRs was queried in L1000FWD. Topoisomerase inhibitors, such as amsacrine and idarubicin, were top records with relatively high similarity scores (Supplementary Table 3).

Figure 5. Detailed chart for therapeutic targets discovery (A) and drug discovery (B). (A) The prognostic potential of differentially expressed genes (DEGs) and RNA modification regulators (RNA-MRs) in cluster1 (severe NBL) were analyzed. Genes correlated to survival formed prognostic DEGs set and prognostic RNA-MRs set, respectively. Then function-related prognostic DEGs were identified from prognostic DEGs as genes positively related to prognostic RNA-MRs. The function-related prognostic DEGs and prognostic RNA-MRs together formed core prognostic gene set. (B) The reversed similarities of gene expression pattern in DEGs and core prognostic gene set were calculated in L1000FWD respectively. The same records in two searches were integrated according to target type. The reversed similarities of core prognostic gene set only was calculated in CMap. Then inhibitors from L1000FWD and CMap were summarized according to target type. PR: progesterone receptor; SFU: sodium fluorescein uptake; MR: Mineralocorticoid receptor; MOP: Mitochondrial oxidative phosphorylation; RNAP: RNA polymerase.

Mapping the NBL signature onto drug-perturbed transcriptome profiles in CMap

The opposed expression profile of core prognostic gene set was queried in CMap. There were 157 compounds with tau> 90 (Supplementary Table 4). Together with Perturbagen Class (PCL) analysis, these compounds were enriched in HDAC inhibitor, PI3K inhibitor, mTOR inhibitor, topoisomerase inhibitor, and CDK inhibitor (Figure 5B and Supplementary Table 5).

Overlapped drugs in CMap and L1000FWD

There were totally 88 potential drugs selected from CMap or L1000FWD for further analysis (Figure 5B). Seven of them were obtained in both CMap and L1000FWD query, namely, palbociclib, wortmannin, FCCP, idarubicin, aminopurvalanol-a, 7b-cis, and WZ-3146 (Figure 6A). They belonged to the following target classes: CDK, PI3K, topoisomerase, EGFR, mTOR, exportin, and mitochondrial oxidative phosphorylation uncoupler (Supplementary Table 6). Palbociclib, wortmannin, FCCP, WZ-3146, and idarubicin were reported to inhibit cell cycle and vascularization in clinic or neuroblastoma cells, which meant our results are convinced [22, 23]. Purvalanol A was used in clinic for NBL [24]. Aminopurvalanol A, a homolog of purvalanol a, reversed core prognostic gene expression, seemed to be a choice for NBL via targeting CDK [25]. Exportin 1 (XPO1) was a nuclear receptor exporter involved in the active transport of transcription factors, tumor suppressor proteins, cell-cycle regulators and RNA molecules [26]. The 7b-cis, inhibitor of XPO1, showed significant therapy potential in both CMap and L1000FWD, was worth for further study [27].

Figure 6. Structures for different kinds of inhibitors. (A) The structure of 7 overlapped drugs in CMap and L1000FWD. (B) Core structures for HDAC inhibitors. (C) Core structures for topoisomerase inhibitors. (D) Core structures for kinase inhibitors.

The MOA and structure analysis for potential inhibitor classes

Among 88 potential predicted drugs, HDAC inhibitors were the largest class. There was a well-admitted pharmacophore model illustrating the structure of HDAC inhibitors (Figure 6B), which consisted of a zinc binding group (ZBG), a linker chain mimicking the lysine side chain, a connecting unit, and a terminal functional “cap” group interacting with the external surface [28]. The ZBGs of most HDAC inhibitors predicted here were hydroxamic acids.

DNA topoisomerase inhibitors were one of the largest predicted classes, which adjusted DNA’s topological structure by cutting, transferring and re-connecting with DNA chains [29]. Here inhibitors for both Topoisomerase I and Topoisomerase II were obtained, which were mainly camptothecin derivatives and anthracyclines (Figure 6C), respectively.

Protein kinases (PKs) were important for signal transduction from extracellular to intracellular, signaling cascade and finally cell response, such as cell proliferation, differentiation, and inflammatory [30]. The preserved hinge region of ATP binding site [31] in all kinases led to the structural similarity of kinase inhibitors (Figure 6D). Inhibitors usually contained a hinge-binding motif, such as a morpholinyl or purinyl substituent, to occupy the adenine binding site of ATP. Hydrophobic pockets were commonly occupied by tolyl group. Allosteric pockets were commonly occupied by terminal piperazinyl-phenyl group, or long-chain such as trifluoromethyl group, which had a hydrophilic group exposing to the solvent region [32]. For kinase consists of “specificity” pocket, such as P-loop in JAK, the terminal group extended into a cleft underneath the P-loop in the N-lobe [33].

Two kinds of PK inhibitors, serine-threonine kinase inhibitors and tyrosine kinase inhibitors, were mainly predicted here. There were four kinds of serine-threonine kinase inhibitors. The mTOR inhibitors were mainly macrolides, morpholino-pyrimidine-phenyl derivatives, and quinoline-imidazoquinoline-phenyl derivatives. The PI3K inhibitors were commonly dual inhibitors and shared similar structures with mTOR inhibitors. The CDKs inhibitors were mainly phenylamino-purine derivatives. Aurora kinase inhibitors were commonly phenylamino-pyrimidine-phenyl derivatives. There were one kind of non-receptor tyrosine kinase (JAK) inhibitor and three kinds of receptor tyrosine kinase inhibitors (FLT3, EGFR and IGF-1) predicted here. JAK inhibitors were triterpenoid and indolone derivatives. FLT3 inhibitors were pyrimidine and staurosporine derivatives. EGFR inhibitors were pyrimidine, quinoline, and quinazoline derivatives. IGF-1 inhibitors were benzimidazole and pyrrolotriazine derivatives (Supplementary Table 6).

Proteomic profile of drugs targeting different targets in CMap

Fifteen compounds with tau> 95 could be find in P100 and GCP Proteomics Connectivity Hubs, and their genomic and proteomic connectivity were shown as similarity matrix (Figure 7A–7C). There were strong expression connections between CDK inhibitors, MTOR inhibitors and topoisomerase inhibitors, and negative connections between them and HDAC inhibitors in genomic and proteomic perturbed profiles.

Figure 7. The relationships among drugs, core prognostic genes, and targets. (A–C) The similarity matrix of genomic, proteomic connectivity for compounds with tau >95 in CMap (A), GCP (B) and P100 (C) database. The map was sorted with the strongest perturbagen connections at the top in red. (D) Schematic representation of the association between RNA modification regulators (RNA-MRs), function-related prognostic DEGs and target of compounds with tau>95. Pink squares: RNA-MRs; green cycles: targets; blue diamond: prognostic DEGs; green diamond: prognostic genes as targets. (E) Schematic representation of targets with a direction relationship with RNA-MRs, together with corresponding inhibitors.

The relationship between targets, function-related prognostic genes and RNA modification regulators

The drug-perturbed profiles at gene and protein levels were clustered into two groups by their connections, which implied there were not only one potential role for RNA-MRs. To confirm this hypothesis, we construct network between drugs, targets, prognostic genes and RNA-MRs (Figure 7D, 7E). Figure 7D showed these genes were clustered into two groups. One cluster was mainly regulated by TET2 and enriched in JAK and PI3K signaling pathways. The other cluster was mainly regulated by NSUN2/ALYREF/ADAR regulators, and enriched in cell cycle and DNA replication process. In Figure 7E, sub-cluster at left was genes enriched in CDK, topoisomerase, MTOR related signaling. They had connections with ALYREF, NSUN2, and ADARB1 via PLK1, BIRC5, CDK1, AURKA. These genes were both drug targets and prognostic genes, which had closely relationship with RNA-MRs. The other sub-cluster at right showed genes enriched in HDAC, JAK, MEK related signaling. They had connections with TET2 via JAK2 and FLT3. Thus, there were two enriched functions of core prognostic genes, and combined use of different inhibitors may be better to reverse expression profile and block disease course.

Experimental validation of the predicted drugs in neuroblastoma SH-SY5Y cells

To experimentally validate the antitumor activities of predicted drugs, as well as their possible regulatory effects on RNA-MRs, we tested 5 of the 7 overlapped drugs in CMap and L1000FWD. Treatments with drugs for 24h inhibited cellular viability of SH-SY5Y cells in a dose-dependent manner (Figure 8A–8E). The IC50 were 1.26×10^-4, 2.65×10^-7, 6.39×10^-6, 9.13×10^-7, 1.92×10^-6 mol/L for purvalanol A, idarubicin, WZ3146, AZD-8055, TG-101348, respectively. These drugs induced decreased expressions of ALYREF, ADAT1, ADAT3, and NSUN2. They also induced increased expressions of ALKBH5 and TET2 (Figure 8F). The validated 5 drugs had no previous report on regulation of RNA-MRs. Here the method was a relatively high accuracy model to predict inhibitors for RNA-MRs.

Figure 8. The cellular viability and dysregulation of RNA modification regulators (RNA-MRs) in human neuroblastoma cell line SH-SY5Y cells. (A–E) Concentration-dependent inhibition of cellular viability by purvalanol A, idarubicin, WZ3146, AZD-8055, TG-101348 (10^-8~10^-5 mol/L) at 24 h, measured by CCK8 assay. (F) The expression levels of ALYREF, ADAT1, ADAT3, NSUN2, ALKBH5, TET2 detected by quantitative real-time PCR. *p < 0.05, **p < 0.01 versus control group.

Evaluation of the role of RNA-MRs in NBL disease course

The mRNA expression data were calculated from HTSeq-FPKM release of RNA-seq data. The log scale was applied for analyzing. The interaction and correlation among RNA-MRs were firstly analyzed by STRING (http://www.string-db.org/) and spearman correlation analysis, respectively. Consensus clustering was a popular algorithm for unsupervised classing in cancer research [47]. RNA-seq sample was then clustered into two groups (cluster1 and cluster2) by analyzing expression of RNA-MRs with “ConsensusClusterPlus” (K-means, 50 iterations, resample rate of 80%, and Pearson correlation) for R v3.6.0.

The Kaplan–Meier method with a two-sided log-rank test was used to compare the overall survival of the patients in two groups. The transcriptome differences between two groups were analyzed using PCA analysis (prcomp function in R). The DEGs between groups were calculated using limma in R. The consensus expression of RNA-MRs and the distribution of clinicopathological features in two groups were shown by pheatmap in R.

Evaluation of CNVs for RNA-MRs and its relationship with their expression profile

When genes are drug targets, their CNVs affected gene expression abundance and protein structure, thus altered drug-binding pocket of receptors, therapeutic efficacy and safety of drugs [21]. Here CNVs were identified using segmentation analysis and GISTIC algorithm. The CNV patterns for each patient, as well as corresponding clinicopathological features were shown by ComplexHeatmap in R. The chromosome structure, location of CNV were represented using RCircos in R. To evaluate effects of CNVs on expressions of RNA-MRs, the expression levels were compared using One-way ANOVA or t-tests among NBLs with different CNV patterns.

Functional differences in two NBL clusters defined by RNA-MRs

To functionally annotate DEGs, GO for BP, MF, and CC, together with KEGG pathway analysis were performed using the clusterProfiler, org.Hs.eg.db genome-wide annotation, and the enrichplot package in R. The top ten items were visualized by ggplot2 package.

To investigate the functional difference of whole transcriptome between clusters, GSEA (Broad Institute) was performed with hallmark gene set (v7.0). Gene sets with nominal P-value <0.05, and FDR <0.25 were considered to be significantly enriched.

Evaluation of prognostic RNA-MRs

The prognostic values of RNA-MRs were evaluated by univariate Cox regression and LASSO analyses (Figure 9). ADAR, ADARB1, ADAT1, ADAT3, NSUN2, ALYREF, TET2, and ALKBH5 were associated with survival significantly (P<0.05). Because L1000 high-throughput gene expression assay directly measured 978 genes and calculated related 12 000 genes instead of all transcriptome, there were only 10 RNA-MRs containing in CMap and L1000FWD. And 3 of them were prognostic regulators (ADAR, ADARB1, ADAT1), most of which were tumor suppressors and unsuitable for inhibitor design. Therefore, query by only RNA-MRs could not get convinced information retrieval of drugs targeting RNA-MRs, and this query result was used just as a reference.

Figure 9. Prognostic RNA modification regulators (RNA-MRs) and their prognostic values. Prognostic RNA-MRs had a P<0.05 in univariate Cox regression analyses (A) and LASSO analysis (B, C). The riskScore calculated using prognostic RNA-MRs generated a convinced ROC curve (D, E).

Preparation a core prognostic gene set for drug-perturbed expression profile analysis

Because the degree of co-expression was extremely high if genes were related to the same function [48], we then calculated co-expression correlation coefficient to obtain downstream or upstream related genes for prognostic RNA-MRs. Firstly, the association between overall survival and DEGs was evaluated by univariate Cox regression analyses, and 515 DEGs (P < 0.05) from total 2113 DEGs were considered as prognostic genes. Secondly, correlation between prognostic genes and prognostic RNA-MRs was evaluated by spearman correlation analysis (correlation coefficient >0.5 and p<0.001), and 213 genes from 515 DEGs were considered as function-related prognostic genes. Finally, core prognostic gene set consisted of 213 function-related prognostic genes and 8 prognostic RNA-MRs (Figure 5A).

Transcriptome proximity analysis in L1000 FWD

To evaluate the importance of RNA-MRs and their downstream or upstream related genes, information retrieval of core prognostic gene set was analyzed as a reasonable reduced representation of total DEGs. The reversed expression profile of total 2113 DEGs and core prognostic gene set were used separately to query L1000 FWD (http://amp.pharm.mssm.edu/L1000FWD). The similarity scores of predicted drugs were given to quantify the similarity (i.e. overlap) between the input pair of up/down gene sets and the signature up/down genes in L1000FWD. Fisher exact test and subsequent Benjamini-Hochberg procedure were used to assess the statistical significance (P value). Then Z-score (z) was calculated for quantifying the deviation of the observed rank from the expected rank. The combined score (c) was calculated as “c=z·log10(P)”. Drugs with higher similarity scores (ranging from 0 to 1, P<0.05) were considered mimicking reversed expression profile of total DEGs or core prognostic gene set, and as potential therapeutically treatments.

Transcriptome proximity analysis in CMap

The CMap (https://clue.io) was a large-scale compendium of functional perturbations (compounds, shRNAs, cDNAs, and biologics) coupled to an information-rich omics profile, which supplied a comprehensive analysis for perturbational class as well. However, there was a technical limit for input gene set (150 up/down-regulated genes). As a complementary of query in L1000FWD, the core prognostic gene set were used to query CMap.

The similarity of the query to each CMap signature was computed and yielded a rank-ordered list of the signatures. Three connectivity score metrics constituted a statistical framework to provide a holistic quantification of the similarity: (1) a nominal P value using the Kolmogorov-Smirnov enrichment statistic; (2) a FDR accounting for multiple hypothesis testing; (3) a connectivity score (tau) comparing reversed core prognostic gene set to perturbagen set (named “Touchstone”, TS) via weighted Kolmogorov-Smirnov enrichment statistical analysis. A tau>90 was used as convinced strong scores. To assess the universal targets of drugs, PCL describing compounds with the same mechanism of action (MOA) or biological functions was used to cluster perturbagens.

Evaluation of merging drugs in CMap and L1000FWD

Firstly, structure and MOAs of drugs retrieved from both CMap and L1000FWD were analyzed. Then, structure and MOAs of drugs were analyzed when they shared same targets: 1) targets retrieved from both CMap and L1000FWD, or 2) drugs with similarity>0.3 in L1000FWD, or 3) PCLs with tau>99 or drug count>5 in CMap. A literature research for these drugs and targets were used to query their use in clinic or experiment for NBL.

Proteomics proximity analysis in CMap

In CMap, two Proteomics Connectivity Hubs, P100 and GCP were queried to evaluate MOA similarities of potential drugs (tau>95). Drugs with different targets induced different protein change profiles, which represented different signaling pathways being affected. A similarity matrix of proteomic connectivity in available drugs was calculated by Pearson correlation.

Network analysis for inhibitors of downstream or upstream effectors

The network analysis was used to analyze relationship between drug and targets, as well as between drug response biomarkers and prognostic RNA-MRs. All targets of drugs (tau>90 in CMap, or merging drugs in L1000FWD) were searched in STRING, and their basic interaction were downloaded. Interaction between function-related prognostic genes and RNA-MRs, and between RNA-MR themselves, were merged together using spearman correlation scores calculated above (score>0.5). Cytoscape was used to construct the network.

Antitumor bioassay

CCK8 assay

The human neuroblastoma cell line SH-SY5Y was cultured in RMPI 1640 medium containing 15% fetal bovine serum (FBS). Cells were seeded into 96-well plates (10 000 per well), treated for 24 h with different concentration of purvalanol A, idarubicin, WZ3146, AZD-8055, TG-101348 (10^-8~10^-5 mol/L, all purchased from Selleck, Shanghai, China). The cells were incubated at 37 °C for 2 h with CCK-8 Kit (Dojindo, Tokyo, Japan). Medium without cells was incubated with CCK8 as blank well. The optical densities (ODs) were checked at 450 nm using a microplate reader. The inhibition rate was calculated as following: (ODdrug-ODblank)/(ODcontrol-ODblank) × 100%.

Quantitative real-time PCR

The mRNA expression of AlyREF, ADAT1, ADAT3, NSUN2, ALKBH5, TET2 and GAPDH were detected with qRT-PCR. Cells were seeded into 6-well plates (300 000 per well), treated with purvalanol A, idarubicin, WZ3146, AZD-8055, TG-101348 (5×10^-6, 7×10^-8, 5×10^-6, 3×10^-7, and 5×10^-6 mol/L) for 24 h. The total RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and the cDNA was reverse transcribed using ReverTra Ace (Toyobo, Osaka, Japan). The data were detected using SYBR Green PCR Master Mix (Toyobo) by ABI 7500 Fast (Applied Biosystems). Relative quantification of mRNA ratio of the target gene to GAPDH was calculated using the 2^−ΔΔCt method. The primers were in Table 3.

Statistical analysis

The chi-square tests were used to compare the distribution of gender, age, INSS, MYCN, ploidy, histology, grade, MKI, and COG risks.

The association between expression of RNA-MRs and clinicopathological characteristics, as well as CNVs and clinicopathological characteristics, were analyzed with chi-square test or Mann-Whitney U test.

All statistical analyses and figures were conducted using R v3.6.0 (https://www.r-project.org/), SPSS 20.0 (IBM, Chicago, USA) and GraphPad Prism 6.0 (GraphPad Software, La Jolla, CA, USA). All statistical results with a P <0.05 were considered to be significant.

Ethics statement and dataset collection

Three datasets, clinical dataset, RNA-seq dataset, and CNV dataset were retrieved from TARGET and TCGA databases (Supplementary Figure 1D). The data was accessed according to open access guidelines where written informed consents were obtained in accordance with the local institutional review boards. Clinicopathological information was summarized in Table 1. To evaluate the information retrieval of patients with RNA-seq and CNV data for overall NBL patients, univariate and multivariate Cox regression analyses were performed, and the association between clinicopathological factors and overall survival of patients in 3 datasets was determined.