The potential role of RNA sequencing in diagnosing unexplained insensitivity to conventional chemotherapy in pediatric patients with B-cell acute lymphoblastic leukemia

Pediatric B-cell acute lymphoblastic leukemia (B-ALL) is a highly heterogeneous disease. According to large-scale RNA sequencing (RNA-seq) data, B-ALL patients can be divided into more than 10 subgroups. However, many genomic defects associated with resistance mechanisms have not yet been identified. As an individual clinical tool for molecular diagnostic risk classification, RNA-seq and gene expression pattern-based therapy could be potential upcoming strategies. In this study, we retrospectively analyzed the RNA-seq gene expression profiles of 45 children whose molecular diagnostic classifications were inconsistent with the response to chemotherapy. The relationship between the transcriptome and chemotherapy response was analyzed. Fusion gene identification was conducted for the included patients who did not have known high-risk associated fusion genes or gene mutations. The most frequently detected fusion gene pair in the high-risk group was the DHRSX duplication, which is a novel finding. Fusions involving ABL1, LMNB2, NFATC1, PAX5, and TTYH3 at onset were more frequently detected in the high-risk group, while fusions involving LFNG, TTYH3, and NFATC1 were frequently detected in the relapse group. According to the pathways involved, the underlying drug resistance mechanism is related to DNA methylation, autophagy, and protein metabolism. Overall, the implementation of an RNA-seq diagnostic system will identify activated markers associated with chemotherapy response, and guide future treatment adjustments. Supplementary Information The online version contains supplementary material available at 10.1186/s12920-024-01892-w.


Background
Acute lymphoblastic leukemia (ALL) is the most common malignant disease in children.A high remission rate has been achieved by multiple-drug chemotherapy as an introduction remission therapy for ALL [1].However, approximately 10% of ALL patients fail to achieve remission, indicating a poor response to chemotherapy and a high risk of relapse and mortality [2].ALL is also a highly heterogeneous disorder.In a previous large-scale RNA sequencing (RNA-seq) study, according to the transcriptional landscape of B-ALL, B-ALL patients were divided into more than 10 subgroups with different fusion genes and known gene mutations, which can be used to predict patient prognosis and help risk stratification [3].The underlying genomic defects and mysterious mechanisms of drug resistance are still unclear.Beyond cytogenetics and fluorescence in situ hybridization, gene expression profiling has revealed new cytogenetic subgroups that display specific gene expression patterns.RNA-seq is believed to be implemented within an individual clinical service to enhance the current molecular diagnostic risk classification of leukemia [4][5][6][7].Gene expression patternbased therapy could be a potential upcoming strategy.
Unfortunately, even though we stratified ALL patients according to genomic information, we did not adjust the therapy greatly.Concerning those without the included fusion gene or known gene mutation, physicians cannot stratify ALL patients without access to the results of the morphology of the bone marrow or minimal residual disease (MRD) at the end of induction therapy [2,8,9].
In recent years, small-molecule targeted drugs have been proven effective in clinical trials, regulating a variety of signaling pathways by targeting signaling molecules that play a regulatory role in cell proliferation, differentiation, and apoptosis [10][11][12][13][14][15][16][17][18][19].For those with drugresistant ALL, further exploration of the gene expression patterns of leukemic cells and translation of gene-gene interactions would be meaningful for gene expression pattern-based drug combination therapy development.Genome-wide association studies/transcriptome-wide association studies (GWASs/TWAS) of pharmacogenomic/pharmacotranscriptomic markers can identify relevant markers regardless of whether their function was previously known [20], but these studies have low statistical power due to the number of independent tests performed.Fusion gene detection and transcriptome analysis are approaches for identifying candidate pharmacogenetic/pharmacotranscript markers [21].
In this study, by retrospectively analyzing the clinical data on induction chemotherapy in the SCCLG-2016-ALL collaborative group (ChiCTR2000030357) and the relationship between clinical data and gene expression profiles, we evaluated the possibility of developing a novel pattern of personalized combination therapy in a group of children with unexplained chemotherapy-poor ALL patients.

Patient characteristics
Forty-five children with common-B ALL were treated.Bone marrow samples from patients before chemotherapy and when their disease met the criteria for remission according to the SCCLG-2016-ALL protocol were used as test samples and control samples, respectively.The risk classification of this protocol is provided in the supplementary information (SI 1).All patients underwent molecular diagnostic analyses, including karyotype and FISH analysis for ETV6::RUNX1, BCR::ABL1, and KMT2A rearrangements.Our 45 samples represented a selected subset deliberately biased toward samples without a positive molecular diagnostic classification, which explained the poor response to chemotherapy.The clinical and cytogenetic characteristics of the cohort are shown in Table 1 and are inconsistent with the established features of childhood ALL [3,4].
The prednisone response was evaluated on Day 8 after induction by peripheral blood leukemic cell counts.MRD was assessed on Day 15 and Day 33 of induction by flow cytometry.The prednisone response test and MRD status are shown in Table 1.The rates of positive MRD on Day 15 and Day 33 were higher than those in all ALL patients managed in our center, reflecting the selection of nonstandard patients.The included cohort comprised 62.2% non-high-risk and 37.8% high-risk B-ALL patients.None of the included datasets were analyzed in previous publications.
Bone marrow/peripheral blood mononuclear cells were extracted, and total RNA was extracted for quality determination.After qualification, cDNA was synthesized by hybridization, and the original data were obtained by RNA sequencing.The original data were converted by statistical software for statistical analysis, and the differential expression profiles of mRNAs between the B-ALL patients and the bone marrow of patients in remission were obtained.Based on the NCBI Ref Seq, UCSC, RNAdb, and other database resources, bioinformatics analysis was conducted on the chip results, such as Gene Ontology and Pathways, and the gene network map of mRNA was established by calculating the Pearson correlation coefficient.

Definitions
In this study, risk classification was based on the treatment response to corticosteroids according to the prednisone test and MRD level on Day 15 and Day 33.After seven days of prednisone, prednisone poor response (PPR) was defined as ≥ 1 × 10 9 /L blasts in the peripheral blood.High risk (HR) (n = 17) was defined as poor  In the analysis of RNA-seq transcriptome data, we defined upregulated genes as genes whose expression increased twofold and whose Benjamini-Hochberg adjusted p values were greater than those of the control, while downregulated genes were defined as genes whose expression decreased twofold and whose Benjamini-Hochberg adjusted p values were greater than those of the control.

RNA-seq
Patient blood and bone marrow samples were obtained from the Cancer Tissue Bank of Sun Yat-sen Memorial Hospital.The details of sample handling, RNA extraction, library preparation, and sequencing parameters are provided in the supplementary information (SI 2).

Fusion detection
RNA sequence fastq.gzfiles were aligned to the human genome (hg38) using STAR aligner (version 2.7.7a) in 2-pass mode with a parsed version of the comprehensive GENCODE 38 annotation.36.The parameter details are provided in the supplemental materials.Detection of the fusion genes was performed using Arriba v1.1.0(https:// github.com/suhrig/arriba/)with aligned bam files as the input data with default parameters.The visualization of fusion genes was performed using circlize v0.4.13 (https://github.com/jokergoo/circlize).

Gene expression analysis
RSEM (v1.3.3) was used to calculate the expression levels (reads count, TPM[transcripts per million], FPKM[fragments per kilobase of transcript per million fragments mapped]) of all the genes in each sample and to generate an expression matrix with the bam files generated in the Fusion detectionsection as the input data.

Gene expression classifier
Classification based on gene expression profiles (TPM values) for each sample was performed using Consensus-ClusterPlus [22].The 5,000 genes with the largest absolute deviations among all samples were selected.The expression levels were normalized to the median.The parameter details of ConsensusClusterPlus are provided in the supplemental materials.

Differential expression analysis
Differential expression analysis was performed by DESeq2 v1.34.0 with a read count matrix as the input data.The significantly differentially expressed genes were defined as genes with adjusted p values less than 0.05 and twofold or greater differences in expression according to the DESeq2 results.Then, all significantly differentially expressed genes were annotated, and pathway analysis was performed with ClusterProfiler (https://github.com/YuLab-SMU/clusterProfiler) or Metascape (https:// Metascape.org/).

Statistical analyses
Comparisons of categorical variables were ascertained by Pearson's χ2 test or Fisher's exact test.Two-sided p values are reported.Analyses were performed with R (v4.0.2).

Clinical characteristics and cytogenetic features of the cohort
The clinical characteristics and cytogenetic features of the cohort are shown in Table 1.In addition to the above information, information on variations in fusion genes, gene mutations, and karyotypes was also collected.Gene , and BTG1 (n = 1) have been reported.Regarding the Ras signaling pathway, FLT3 (n = 5), KRAS (n = 12), and NRAS (n = 5) were commonly detected alone or together in 14 patients.Hyperdiploidy of chromosomes + 4, +10, and/or + 17 was detected in 9 patients.ZNF384-TCF3 fusion (n = 1), EP300-ZNF384 (n = 1), and ZNF384-SYNRG (n = 1) were positive in 3 patients at the first diagnosis.Since most included patients did not have known high-risk associated fusion genes or gene mutations, fusion gene identification was conducted.

Comparison of fusion genes between the high-risk group and the non-high-risk group
According to the initial chemotherapy response defined by the prednisone test and MRD, but not by fusion genes or gene mutations, we classified the cohort into a highrisk group and a non-high-risk group.To determine the differences in biological characteristics between the highrisk group and the non-high-risk group, the differential expression profiles were compared.Clinical characteristics and cytogenetic feature comparisons between the high-risk group and the non-high-risk group are shown in Table 2.
To identify potential high-risk-related fusion genes, the genes that were differentially expressed between the highrisk group and the non-high-risk group were further analyzed (Fig. 3).Fusions involving ABL1, LMNB2, NFATC1, PAX5, and TTYH3 at onset were more frequently detected in the high-risk group, while fusions involving LFNG, TTYH3, and NFATC1 were frequently detected in the relapse group.In addition, there are several newly discovered fusion genes with high frequencys, such as DMD::AL591378, MRPL13::ZNF740, and DHRSX duplication (Fig. 4).DMD::AL591378 fusion was frequently found in onset samples and disappeared at remission (Fig. 4A).The MRPL13::ZNF740 fusion was frequently found in the non-high-risk group at onset, but it did not disappear after treatment (Fig. 4B).DHRSX-DHRSX presented a high frequency in the high-risk group (Fig. 4C).DHRSX-DHRSX is AL683807.2(442,713) and DHRSX (8017)-DHRSX.During analysis, when fusion sites are located between genes, upstream or downstream genes are selected according to the location of fusion sites in the transcript to represent the fusions.In the output result, for DHRSX-DHRSX, AL683807.2(442,713) is the only form, while DHRSX (8017)-DHRSX does not exist.According to the results, at the locus 8017 upstream of DHRSX, there is a fusion within the DHRSX gene, which is a duplication.

Elevated gene expression and low gene expression in the cohort
Compared to the gene expression profile before remission, 3824 genes were significantly upregulated and 5853 genes were significantly downregulated at remission (Fig. 5A).The differentially expressed genes were annotated functionally.The most highly upregulated pathways are illustrated in Fig. 5B and are mainly involved in the activation of neutrophils in the immune response, hemostasis, and stress response.The most downregulated pathways are illustrated in Fig. 5C and are involved mainly in herpes simplex virus infection and DNA methylation.DNA methylation played a very special role in this cohort.Moreover, the genes related to homophilic cell adhesion via plasma membrane adhesion molecules, development growth, regulation of cell differentiation, and B-cell proliferation were involved in leukemia progression in the cohort.By comparing the overlap between these DEGs and fusion genes, we identified several shared genes (Fig. 5D and E).However, the overlaps were not significant.

New classifier for identifying subgroups in B-ALL patients
There were 102 upregulated genes and 56 downregulated genes in the high-risk group (Fig. 5F).The pathways involved are illustrated in Fig. 5G-H.Pathways related to exocytosis, cell chemotaxis, inflammatory response, and cell proliferation and survival were ultimately upregulated.These differences are likely attributed to the leukemia stem cells and relapse clones retained in the high-risk group.The pathways related to the regulation of ion transmembrane transport, potassium ion transmembrane transport, negative regulation of cellular amide  metabolic processes, and protein-DNA complex assembly were suppressed, which likely suppressed protein metabolism.By comparing the overlap between these DEGs and fusion genes, we identified several shared genes (Fig. 5I and J).However, the overlaps were not significant.RAS signaling pathway gene mutations are controversial factors regarding patient prognosis.Thus far, we do not know the necessity of targeted therapy for ALL patients with RAS signaling pathway gene mutations, nor do we know the targets that should be focused on for ALL patients with RAS signaling pathway gene mutations.Therefore, we investigated the general characteristics of the RAS signaling pathway gene mutation subgroup.Thirteen patients carrying RAS signaling pathway genes were identified in the cohort, 4 of which presented high-risk characteristics, while 9 presented of nonhigh-risk characteristics.There were 193 upregulated genes and 108 downregulated genes in the RAS signaling pathway gene mutation subgroup (Fig. 6A).The pathways involved are illustrated in Fig. 6B-C.By comparing the overlap between these DEGs and fusion genes, we identified several shared genes (Fig. 6D and E).However, the overlaps were not significant.
For comparison between the RAS signaling pathway genes mutated high-risk and RAS pathway genes mutated non-high-risk, there were 193 upregulated genes and 108 downregulated genes in the RAS signaling pathway gene mutation subgroup (Fig. 6F).The pathways involved are illustrated in Fig. 6G-H.By comparing the overlap between these DEGs and fusion genes, we identified several shared genes (Fig. 6I and J).However, the overlaps were not significant.Interestingly, the fusion genes MRPL13::ZNF740 (n = 4), DHRSX duplication (n = 4), CYR5A::DIPK1C (n = 4) and CXorf21::MRPS16A Fig. 3 Comparison of differentially expressed fusion genes between the high-risk group and non-high-risk group.The table on the right shows detailed information on 22 genes shared among the high-risk group, onset group and relapse group, but not among the genes in the medium-risk group.The number in the table indicates the number of patients with the corresponding fusion genes

Characteristics
Non-high risk (n =  (n = 3) were dominant in RAS signaling pathway-mutated B-ALL (Fig. 6K).For B-ALL without known molecular markers, molecular subtyping according to expression profiles would help identify new subgroups (Fig. 7).The AUCs of the consistency indices between K = 3 and K = 10 were not significantly different (Fig. 7A-C).Therefore, they were subtyped into three groups: Group 1, Group 2 and Group 3 (Fig. 7D).Differentially expressed genes were compared among these subgroups.A pairwise comparison of  the three subgroups is illustrated in Tables 3 and Fig. 8.A comparison of Group 1 and Group 2, revealed and the DEGs and pathways involved in these genes affected RNA polymerase I promoter opening, eukaryotic translation elongation, and the electron transport chain of the OXPHOS system in mitochondria (Fig. 8A).A comparison of Group 1 and Group 3 revealed that the DEGs and involved pathways affected eukaryotic translation elongation and the electron transport chain of the OXPHOS system in mitochondria (Fig. 8B).The expression profiles of Group 1 and Group 3 were very similar, as shown by the transcriptome data.Group 3 should be regarded as deriving from Group 1.A comparison of Group 2 and Group 3 revealed that the DEGs and pathways involved affect eukaryotic translation elongation and that HDACs deacetylate histones (Fig. 8C).Gene Cluster 1, involving pathways of SRP-dependent cotranslational protein targeting membrane and signal sequence recognition, was upregulated in Group 1. Gene Cluster 2, involving pathways of mRNA 5'-splice site recognition pathways, NABA proteoglycans, and nucleosome assembly, was upregulated in Group 3. The expression of genes in Cluster 3, involving antigen receptor-mediated signaling pathway, mitochondrial transmembrane transport, eukaryotic translation elongation, and mitochondrial protein pathways, was upregulated in Group 2 (Fig. 8D).
To determine whether these pathways are related to the therapeutic response, we further validated the model by comparing MRD and prednisone responses between clusters and subgroups (Fig. 8E).Group 1 is more likely to be high-risk, while Group 3 is less likely to be highrisk.Although these two groups showed similar downregulated transcriptomes, the differential expression of gene Cluster 3 was likely the most important factor related to the different treatment responses.The upregulated pathways in the high-risk group were related to the circulatory system process, positive regulation of cell motility, signaling by GPCRs, tube morphogenesis, NABA microsome association, regulation of cell adhesion, and response to the hormones.The downregulated pathways in the high-risk group were related to the immune response-regulating signaling pathway.To eliminate the bias caused by mutations of P16 deletion, KRAS mutations, and NRAS mutations, which were frequent in the cohort, the distribution differences were compared between subgroups, and the results presented no significant difference (p > 0.05).

Discussion
This is an exploratory study, with the main finding that RNA sequencing can reveal previously unreported fusion genes in some pediatric B-ALL patients whose chemotherapy response is inconsistent with the risk stratification corresponding to routine molecular biological features.Abnormal regulation of previously neglected signaling pathways can also be found.These new abnormalities may be related to pathogenesis and treatment response.Further comparison of groups of patients with different treatment responses may reveal new classifications of expression profiles.A novel approach for targeted treatment of B-ALL based on the transcriptome is discussed.
In a previous study, according to the transcriptional landscape, B-ALL patients were divided into subgroups with different fusion genes and known gene mutations, which included MEF2D fusions, TCF3-PBX1, ETV6::RUNX1 or ETV6::RUNX1-like, DUX4 fusions, ZNF384 fusions, BCR::ABL1 or BCR::ABL1like, high hyperdiploidy, KMT2A fusions, PAX5 and CRLF2 fusions, PAX5 (p.P80R) mutations, IKZF1 (p.N159Y) mutations, ZEB2 (p.H1038R)/IGH::CEBPE, TCF3/4::HLF, and NUTM1 fusions [3].Gene mutations among signaling molecules, epigenetic factors, and transcription factor genes are enriched in different subgroups, which means that, except for subgroup identifiers, enriched mutations in different functional genes are attributed to the gene expression files together [23].However, regarding a random individual patient, it does not always copy the genomic and epigenetic characteristics of certain subgroups, and there would be confounding variations [24][25][26][27][28].It is still necessary to discover the underlying genomic defects and mysterious mechanisms in individuals for treatment response [29].By using an intertomics-based approach, transcriptome-based subcluster identification can be achieved directly at the clinical level, which may be critical for predicting potential drug targets for individual subclusters [5].Therefore, RNA-seq can be implemented in the diagnostic workflow of ALL and enhances the individual molecular diagnostic risk classification of ALL [4,30].
In this study, we explored the gene expression pattern in a group of childhood B-ALL patients whose molecular biological characteristics and chemotherapy response did not completely match the above subgroup identifiers.RNA-seq-based subgroups were used stratify patients into low-risk, intermediate-risk, and high-risk groups according to survival [3].BCR::ABL1-positive B-ALL is usually regarded as a "high-risk" reference [31][32][33], while  ETV6::RUNX1-positive B-ALL is a "low-risk" reference [34].Subgroups with TCF3::PBX1, ETV6::RUNX1-like, ZNF384 fusions, DUX4 fusions [34], or high hyperdiploidy are low-risk according to the BCR::ABL1-positive subgroup and ETV6::RUNX1-positive subgroup.PAX5 and CRLF2 fusions and hyperdiploidy (≤ 50 chromosomes) are classified into the intermediate-risk group due to an inferior 5-year overall survival compared with that of ETV6::RUNX1-positive/ETV6::RUNX1-like.MEF2D fusions and KMT2A fusions tend to be associated with high risk [3].Usually, the low-risk group responds well to chemotherapy.It does not need highly intensified treatment to induce or maintain long-term remission [35,36].
In contrast, the high-risk group is often insensitive to chemotherapy and needs intensified treatment with chemotherapy, immunotherapy, or small molecular targeted medicine [37].We did not stratify B-ALL patients simply according to the above criteria.Post-chemotherapy treatment response evaluation through prednisone response and/or MRD corresponds well to the long-term prognosis.Therefore, in this study, we retrospectively evaluated the included patients in terms of prednisone response and MRD status, ignoring the molecular biological characteristics indentified by the enrollment test, to analyze the transcriptome characteristics related to treatment response-.
Gene expression profiling has revealed new cytogenetic subgroups that display certain specific gene expression patterns.In this study, DNA methylation played a very special role in this cohort.Moreover, the relationships between genes involved in homophilic cell adhesion via plasma membrane adhesion molecules, development growth, regulation of cell differentiation, B-cell proliferation, and leukemia progression in the cohort are not well understood and should be further evaluated for treatment designation.There is a potential role for DNA methylation-targeting therapy in this cohort.Furthermore, pathways related to exocytosis, cell chemotaxis, inflammatory response, and cell proliferation and survival were upregulated, which was probably attributed to leukemia stem cell and relapse clone retention in the high-risk group.The suppressed pathways related to the regulation of ion transmembrane transport, potassium ion transmembrane transport, negative regulation of cellular amide metabolic processes, and protein-DNA complex assembly likely suppress protein metabolism and provide self-protection.
RAS signaling pathway mutational status of NRAS, KRAS, and PTPN11 genes is associated with genetic/ cytogenetic features in children with B-precursor acute lymphoblastic leukemia.Previous reports have shown that 70% of ALL patients with damaging germline ETV6 variants exhibit hyperdiploid karyotypes with characteristic recurrent mutations in NRAS, KRAS, and PTPN11 [38].The RAS signaling pathway is related to the pathogenesis, prognosis, and relapse process of B-ALL.However, in the analysis of the subgroup expression profile of the cohort, we found that the distribution of RAS signaling pathway mutations among subgroups was not significantly different, which indicated that RAS signaling pathway mutations did not play key roles in controlling the variations in the expression profile.This finding does not support the idea that the mutation status of the RAS signaling pathway may be involved in therapy designation in the future.
Fusion gene detection of candidate pharmacogenes/ pharmacotranscripts cannot revealed new genes or transcripts [21].The transcriptomes of fusion genes may be more informative of disease development, cell proliferation, differentiation, and apoptosis.In this study, we identified new fusion genes without identifying their biological functions.Seven patients in the highrisk group carried DHRSX duplications.According to a previous report, using overexpression and knockdown analyses, Zhang et al. showed that DHRSX promoted starvation-induced autophagy in HeLa and U2OS cells [39].The promotion of autophagy by DHRSX involves the downregulation of AKT/mTOR phosphorylation and the upregulation of beclin-1, which are highly related to the antiapoptotic function and chemotherapy resistance of leukemia cells [40][41][42].This result provides new insight for further exploration of a new mechanism for drug resistance in B-ALL.Whether DHRSX duplication enhances the function of starvation-induced autophagy in B-ALL remains to be determined in future studies.
Although ALL patients are stratified according to genomic information, we cannot adjust the therapy without MRD evaluation at the end of induction therapy [2,8,9].Treatment of childhood refractory and relapsed acute lymphoblastic leukemia (R/R ALL) patients with refractory and poor prognostic genes cannot be resolved by chemotherapy alone [43].In recent years, smallmolecule targeted drugs have been proven effective in clinical trials [44][45][46][47][48]. Small molecule targeted drugs regulate a variety of signaling pathways by targeting signaling molecules that play a regulatory role in cell proliferation, differentiation, and apoptosis [10][11][12][13][14][15][16][17][18][19].For example, targeting BCR::ABL1 and BCR::ABL1-like [29,46,[49][50][51] significantly affects patient prognosis.Therefore, RNA-seq-based therapy adjustments may be more beneficial in the ETV6::RUNX1-positive subgroup.The "low-risk" ETV6::RUNX1-positive subgroup was not always low-risk, similar to the patients in this study.It has been reported that the intensity of chemotherapy [52] in the ETV6::RUNX1-positive subgroup should not decrease for those who have positive MRD postinduction chemotherapy.Discontinuation of L-asparaginase and poor response to prednisolone is associated with poor outcomes in ETV6::RUNX1-positive pediatric B-ALL patients, indicating drug insensitivity [53].Autophagy inhibition [54], spleen tyrosine kinase (SYK) [55], and IGF2BP1 [56] are potential future therapeutic targets for ETV6::RUNX1-driven B-cell precursor acute lymphoblastic leukemia due to their roles in ETV6::RUNX1 cell survival and prognosis.Through RNA-seq, we can identify upregulated and activated pathways at diagnosis and design combination chemotherapeutics to sensitize resistant primary cells to conventional drugs.There were 4 ETV6::RUNX1-positive high-risk patients for whom further analysis of the mechanism underlying the poor response to therapy was needed.For those with drugresistant ALL, further exploration of the gene expression patterns of leukemic cells and translation of gene-gene interactions would be meaningful for gene expression pattern-based drug combination therapy development.

Conclusions
Although large cohort studies have systematically classified ALL, some patients often cannot use these classification models to explain their illness and treatment response.When their insensitivity to treatment cannot be explained by the existing classification system, we can analyze their expression profiles through RNA-seq by implementing a diagnostic system to identify activated markers related to leukemogenesis and drug resistance.These findings will guide the selection of targeted drugs and the design of treatment plans in the future.The implementation of an RNA-seq diagnosis system and a new concept for treatment deserve a prospective study.

Fig. 1
Fig. 1 Overview of gene fusions in 45 pediatric B-ALL patients in the cohort.A. Illustration of the comparison between fusion break points at onset and at the end of remission induction.B-D.The locus distribution of different mutation types on different chromosomes before and after remission

Fig. 2
Fig. 2 Overview of gene fusions in 45 pediatric B-ALL patients in the cohort.(A) Comparison of break points in patients before and after remission.(B) Fusion gene pairs in B-ALL detected in the onset samples.(C) Fusion gene pairs in non-high-risk B-ALL patients detected in the onset samples.(D) Fusion gene pairs in high-risk B-ALL patients detected in the onset samples.(E) Fusion gene pairs in B-ALL detected in remission samples

Fig. 4
Fig. 4 Newly discovered fusion genes with high frequency.(A)DMD::AL591378 fusion in onset samples, but not in remission samples.(B)MRPL13::ZNF740 fusion in onset samples and remission samples.(C)DHRSX::DHRSX in the high-risk group

Fig. 5 Fig. 6 Fig. 7
Fig. 5 Genetic expression alterations in 45 pediatric B-ALL patients in the cohort.(A) Volcano plot comparing the gene expression profiles before and after remission.(B) The most highly upregulated pathways.(C) The most downregulated pathways.(D) Venn diagram showing the overlap between fusion genes and upregulated genes after remission.(E) Venn diagram showing the overlap between fusion genes and downregulated genes after remission.(F) Volcano plot comparing the gene expression profiles in the high-risk group before and after remission.(G) The most highly upregulated pathways in the high-risk group.(H) The most downregulated pathways in the high-risk group.(I) Venn diagram showing the overlap between fusion and upregulated genes in the high-risk group after remission.(J) Venn diagram showing the overlap between the fusion genes and downregulated genes in the high-risk group after remission

Fig. 8
Fig. 8 Expression profiles in molecular subgroups.A-C.Pairwise comparison of expression profiles between the two subgroups and the involved pathways.D. Comparison of expression profiles among the three subgroups and the involved pathways.E. Validation of the molecular subtyping model by comparing MRD and prednisone response between clusters and subgroups

Table 1
Baseline Characteristics of Study Participants

Table 1 (
continued) variations, including exon deletions or single nucleotide variations, in CDKN2A

Table 2
Baseline Characteristics of Study Participants in the analysis of fusion genes

Table 3
Comparison for gene expression across groups of different mRNA expression profiles Note: In parentheses are the number of annotated genes