A comprehensive landscape of transcription profiles and data resources for human leukemia

Leukemia is the main type of hematologic cancer that arises from malignant proliferation and differentiation disorders of blood cells, resulting in abnormal hematopoietic function.^1,2 Generally, leukemia can be classified based on the differentiation stage and the origin of malignancies²; for example, acute lymphoblastic leukemia (ALL), acute myeloid leukemia (AML), chronic lymphocytic leukemia (CLL), and B-cell ALL (B-ALL). Different leukemia subtypes show distinct clinical characteristics and transcriptomic profiles,^3,4 indicating diverse mechanisms underlying their leukemogenesis at the transcription level. Thus, an integrated analysis of leukemia expression data can provide comprehensive insights for elucidating the common and specific regulatory mechanisms among different leukemia subtypes.

Previous studies on leukemia have demonstrated that abnormal gene expression, alternative splicing (AS), and gene fusion play critical roles in the initiation and progression of leukemia.^5,6 For example, the BCR-ABL1 fusion event encodes a continuously activated ABL kinase, which leads to disruptions of key cellular processes and leukemogenesis.⁵ Classical AS of AML1 could produce proteins with different domain structures that regulate cell growth and/or differentiation in leukemogenesis,⁷ and the circular RNA (circRNA) circMYBL2 (the back-splicing [BS] of MYBL2) could regulate FLT3 translation by recruiting PTBP1 to promote AML progression.⁸ Besides, the aberrant expression of long noncoding RNAs (lncRNAs) has demonstrated strong associations with recurrent mutations, clinical features, and outcomes in AML.⁹ Thus, systematic investigation of the expression profiles and alterations of RNA molecules in leukemia samples could improve the understanding of the dysregulated regulation underlying leukemogenesis and facilitate the characterization of the transcriptional features and biomarkers of different leukemia subtypes.

However, most previous leukemia studies focused on 1 specific leukemia subtype or a limited number of subtypes, resulting in massive independent datasets with batch effects, which formed barriers in terms of integrated investigation and systematic comparison of various leukemia subtypes. In addition, further use of these data is severely hampered by a lack of integrative tools for visualization and analysis. Although a few of resources were developed to address this problem, such as LeGenD and Leukemia Gene Atlas,¹⁰ the chief limitations of these resources lie in the small sample size and out-of-date information. For example, the data of Leukemia Gene Atlas were based on gene expression from microarrays from before 2012, and this platform is currently unavailable because of maintenance, whereas the LeGenD database only has deposited information on 70 genes. Therefore, a comprehensive leukemia transcriptome data resource, which contains gene expression, AS events, and fusion profiles of various types of RNAs across leukemia subtypes, is indispensable and an urgent need of leukemia researchers.

Thus, in this study, we attempted to address this issue by developing the LeukemiaDB data repository (http://bioinfo.life.hust.edu.cn/LeukemiaDB) to curate a comprehensive transcription landscape of leukemia from large-scale RNA sequencing (RNA-seq) data (>3000 samples) of 14 subtypes and 53 related cell lines. In addition, we systematically investigated the expression characteristics and structure variants of RNA molecules across patients with, and cell lines of, different leukemia subtypes.

RNA-seq data were collected from the NCBI GEO and SRA databases using the keywords "leukemia AND Homo sapiens AND RNA-seq." The detailed metadata of samples were retrieved using crawler modules of Beautiful Soup 4.4.0 and urllib2 in Python. The disease information of samples was manually confirmed based on the metadata and related published papers, and irrelevant samples (eg, nonleukemia or ambiguous disease) were discarded. After these procedures, a total of 3504 samples were collected from public data portals. Transcriptome data were downloaded and unpacked through SRA Toolkit (version 2.9.0-ubuntu64), of which data from the non-Illumina platform were removed. Adapter sequence removal and reads trimming were accomplished via Trimmomatic¹¹ (version 0.32) with default parameters, and FastQC (version 0.11.5) was used for quality control. HISAT2 (version 2.0.4)¹² was used to quickly align clean reads to the human reference genome GRCh38. Samples with an alignment rate of <75% were excluded from our study.

To better estimate the abundance of lncRNAs, we merged the lncRNA annotations files of Ensembl version 84 and NONCODE version 6.0 by GffCompare (version 0.9.8) as previously described.¹³ In brief, class codes from the output of GffCompare were used to explore relationships between lncRNAs of NONCODE and genes of Ensembl: (1) lncRNAs of NONCODE with class codes of “=” and “c” were deemed redundant transcripts from Ensembl, and their NONCODE IDs were replaced by Ensembl ones (their expression will be merged with Ensembl genes); (2) lncRNAs with other types of class codes were considered novel lncRNAs compared with those found in the Ensembl annotation. In addition, we inspected flanking genes around lncRNAs via a sliding window using bedtools (version 2.17). Finally, the merged GTF annotation file was used for abundance estimation of lncRNAs in further analysis.

Raw counts of protein-coding genes and lncRNAs were obtained via featureCounts (version 2.0.3),¹⁴ and the abundance was estimated via StringTie (version 1.3.0) using the GRCh38 and merged GTF annotation file mentioned earlier. Samples with a small amount of detected lncRNA (<1000) and coding gene (<3000) were excluded, with the aimed of ensuring that samples comparable and have enough molecule information for the further analysis. Finally, 3068 samples were included in LeukemiaDB. Moreover, to better integrate the expression data from different projects, batch effect removal (BER) was performed on protein-coding genes or lncRNAs via the R package ComBat function¹⁵ for patient and cell line samples, respectively. Principal component analysis and uniform manifold approximation and projection (UMAP) were used to evaluate and visualize the power of the BER strategy. Principal component analysis was performed via R package FactoMineR and factoextra, whereas UMAP was implemented using R package umap with the “n_neighbors = 15” and “method = naive” parameters. The results of UMAP analysis for samples (the axes information in UAMP) were deposited in the supplemental Material 1.

Specifically expressed genes (SEGs) were detected via the SEGtool¹⁶ using the sensitive mode (detect_mod = 1) for each cell line or leukemia subtype. Here, the average expression matrix of genes in all samples of the given condition (eg, the chronic myeloid leukemia [CML] disease) was used as input for SEG detection. The symbol SEG_high indicated that the gene was specifically expressed in a small number of leukemia subtypes, whereas SEG_low represented that the expression level of the gene was at the extremely low level in the given condition(s) when compared with that of other genes.

To achieve reliable precision and balanced performance, CIRI2 (version 2.0.6),¹⁷ CIRCexplorer2 (version 2.3.6),¹⁸ and circRNA_finder¹⁹ were combined to detect circRNAs. Firstly, samples with a number of BS reads <10 were excluded from circRNAs detection. Furthermore, circRNAs were retained if they fit these criteria: (1) circRNAs with BS reads ≥2; (2) circRNAs detected by at least 2 tools; and (3) circRNAs found in >10% of samples (at least 3 samples) in a certain leukemia subtype or cell line. In addition, circRNAs with ≤2 mismatches across 3 tools were deemed identical. The abundance of circRNAs was estimated by BS reads and normalized using reads-per-million mapped reads. We defined circRNA observed in only 1 leukemia subtype or cell line as a specific circRNA.

Two types of coexpression pairs (messenger RNA [mRNA]-lncRNA and mRNA-circRNA) were identified. The candidate coexpression pairs were defined by a |Spearman correlation coefficient| (|R|) > 0.5 and a P value < .01. The detailed information of 1665 human transcription factors (TFs) and 1025 TF cofactors was obtained from AnimalTFDB 3.0,²⁰ and a number of 2 712 247 TF–target pairs were retrieved from the hTFtarget database.²¹ 

The rMATS (docker version: rMATS-docker 0.1beta) was used to detect 5 types of AS events²²: (1) alternative 5′ donor sites (A5SS), shortened 5′ exons via internal splice sites; (2) alternative 3′ donor sites (A3SS), similar to A5SS but occurring in 3′ exons; (3) skipped exon (SE) means exon(s) excluded in the transcript; (4) retained intron (RI) means intron retained in transcripts; and (5) mutually exclusive exons, 1 exon is included whereas the other excluded compared with that in putative transcripts.²³ We also defined the splicing ratio of a AS event, which is the ratio of samples with this AS event to all samples of a cell line or leukemia subtype. AS events that met the following criteria were retained: (1) the AS events were detected in at least 3 samples and (2) occurred in >40% samples of a leukemia subtype or cell line. The percent splice in value (calculated by rMATS) and a Student t test (adjusted P values using Benjamini Hochberg method) were used to assess the level and significance of AS events. AS events that occurred in only 1 leukemia subtype or cell line were defined as specific.

Fusion genes were detected from paired-end RNA-seq data using EricScript²⁴ with GRCh38 reference, in which, R (version 3.4.3), Samtools (version 0.1.19), BWA (version 0.7.17), and Seqtk (version 1.3-r106) were used. The EricScore estimated using EricScript was identified as the fusion score (a range of ∼0-1), which was positively correlated to the reliability of fusion events. The occurrence rate for a fusion event was calculated as the number of fusion-detected samples divided by the total number of samples in a given leukemia subtype. The reliable fusion genes were retained based on the following criteria: (1) fusion genes with spanning reads ≥2, crossing reads ≥1, fusion score >0.7 and (2) occurred in at least 2% or 15% (the higher one) of samples of a leukemia subtype or cell line. Fusion events that occurred in only 1 leukemia subtype or cell line were considered specific.

The Kaplan-Meier survival analysis was carried out using R packages (survival and survminer), whereas the log-rank test was used to assess the statistical significance of the survival curves in 2 user-defined groups. The relationship between gene expression and drug response was evaluated using Spearman correlation coefficient of the 50% inhibitory concentration data and gene expression profiles, as we have previously used in GSCALite,²⁵ which were retrieved from cancer cell lines from the Genomics of Drug Sensitivity in Cancer (GDSC)²⁶ and the Cancer Therapeutics Response Portal (CTRP) databases.²⁷ Drug-gene pairs with |R| > .2 and a false discovery rate of <0.05 were curated in LeukemiaDB, and a gene with positive (or negative) Spearman correlation coefficient in a given drug-gene pair was defined as drug resistant (or sensitive) in that condition.

LeukemiaDB was constructed with Python Flask-RESTful API frameworks (https://flask-restful.readthedocs.io/), AngularJS structural framework (https://angularjs.org), and Bootstrap (version v3.3.7, a free and open-source CSS/HTML front-end framework, https://getbootstrap.com/). A source-available and cross-platform NoSQL database program MongoDB (version linux-x86_64-3.3.6) (https://www.mongodb.com/) was used for data deposition. The LeukemiaDB provides a stable and user-friendly service with the Apache HTTP Server.

In this study, we depicted the transcriptome profiles of 3036 samples from 187 public data sets, including 14 leukemia subtypes, 53 leukemia cell lines, and 92 normal samples (Figure 1A-B). The workflow of the transcriptome analysis is shown in Figure 1C, and the workflow is described in detail in “Methods.” Previous studies demonstrated that batch effects, which often result in discrepancies in the statistical distributions across data from different data sets, may be the largest barrier to integration analysis for data from different studies.²⁸ Thus, the well-treated (BER) data are the basis of integration analysis, and all deposited data in our LeukemiaDB were processed with BER.

After the BER procedure, we found that most samples were well clustered based on their disease subtypes (supplemental Figures 1A-B, 2A-B, and 3), especially for CML, acute monocytic leukemia, and acute megakaryoblastic leukemia, implying the success of our BER strategy. We investigated the expression profiles of protein-coding genes and lncRNAs for each leukemia subtype and cell line. Here, we detected ∼20 000 protein-coding genes in all samples from patients with leukemia or leukemia cell lines, of which nearly 5% were highly expressed with a transcripts per million (TPM) value of >100, and 25% to 35% were expressed at a low expression level (TPM, 0-1; Figure 2A-B). In addition, ∼61 000 and ∼46 000 lncRNAs were detected in patients with leukemia and leukemia cell lines, respectively, of which ∼95% were expressed with a TPM value of <1 (Figure 2C-D). Different leukemia subtypes or cell lines (excluding atypical CML [aCML] and the ARH77 cell line) exhibited similar expression distribution, in which ∼5% of the protein-coding genes were highly expressed with TPM values of >100 (supplemental Figure 4A-B), and ∼95% lncRNAs were lowly expressed with TPM values <1 (supplemental Figure 4C-D).

In addition, high-quality circRNAs were identified for each leukemia subtype (see “Methods”). Although the numbers of detected circRNAs varied greatly among the 3 tools used (Figure 2E), we found that the distributions of circRNA total BS counts in all samples were similar (supplemental Figure 4E). In total, 128 138 circRNAs were identified by at least 2 tools, of which only 6.9% (n = 8882) were retained via rigorous strategies (Figure 2E). Most (n = 6105) of the 8882 circRNAs were derived from cell lines, of which 2527 circRNAs were shared with patients (Figure 2F). For each cell line and leukemia subtype, the number of detected circRNAs varied greatly, ranging from 0 to 4574 (supplemental Figure 4F) in cell lines and 0 to 2403 in leukemia subtypes (supplemental Figure 4G).

We also investigated the AS transcripts in leukemia samples. Here, we found that the SE event was the most prevalent AS event among 5 AS types (47.9%, n = 13 607), and the RI event was the least prevalent (9.9%, n = 2821; Figure 2G). In addition, we found that 459 genes covered 5 AS types in all leukemia samples (Figure 2G), which may be prone to splicing genes in leukemia. Although the numbers of AS events varied greatly for different leukemia subtypes or cell lines (supplemental Figure 4H-I), the distribution of AS events in different leukemia subtypes or cell lines exhibited a similar tendency to the earlier description, that is, SE events accounted for nearly 50% of the total AS events (supplemental Figure 4J-K). Interestingly, we found that the AS profiles of the RS4;11 cell line differed from those of other cell lines, and the occurrence of the RI event accounted for ∼50% of the total AS events (supplemental Figure 4J).

Fusion genes were also investigated in paired-end sequencing samples from patients with leukemia and leukemia cell lines (Figure 2H-I). Here, we observed that 68 fusion genes were identified in patients with B-ALL, and there were <10 fusion genes in patients with AML and CML (Figure 2H). Notably, far more fusion genes were observed in the HL−60/S4 cell line compared with other cell lines (Figure 2I), which may be because of the characteristics of the samples; the HL−60/S4 samples were taken from a single study that had high sequencing depth, thus, more stable fusion genes were detected.

We integrated the transcriptome profiles of functional RNAs in leukemia subtypes and cell lines as well as the TCGA-AML survival information and drug sensitivity data from GDSC and CTRP to develop the comprehensive LeukemiaDB. The LeukemiaDB data portal provides information about 5 types of molecules, that is, protein-coding genes, lncRNA, circRNA, fusion genes, and AS, by integrating expression abundance, coexpression modules, TF-target regulation, and leukemia-specific molecules in leukemia subtypes or cell lines (Figure 3A-B). In addition, detailed sample information is displayed by LeukemiaDB or linked to relevant websites, and LeukemiaDB also curates the relationships between gene expression and patient survival from TCGA-AML data as well as the associations between gene expression and drug resistance from GDSC and CTRP data. Three main sections are provided for searching, browsing, and data visualization. With the “browse by modules” section, users can click a diagram and access the molecules of interest directly. For example, the protein-coding gene module allows users to quickly search gene expression profiles via a gene symbol or Ensembl ID. The protein-coding gene module also allows users to browse the average expression of a target gene in 14 leukemia types (Figure 3B). With the “browse by subtypes” section, users can easily obtain all transcriptome profiles of a selected leukemia type, including the expression level and regulatory information of protein-coding genes, lncRNAs, and circRNAs as well as AS events, fusion genes, and subtype-specific molecular profiles (Figure 3B). Via this section, LeukemiaDB provides the aforementioned transcriptomic molecular features for a given cell line (Figure 3B).

Further comparative analyses were performed to identify commonalities and differences among different leukemia subtypes and cell lines. Firstly, we detected SEGs in each leukemia subtype and cell line. In a previous study, we demonstrated that SEG_high genes may be crucial molecules in oncogenesis and serve as potential biomarkers for the diagnosis and prognosis of cancers.²⁹ Interestingly, in this study, we observed that the number of SEG_high genes was predominant both in cell lines and patient samples (Figure 4A-B), which may provide important clues for leukemia research. For example, NOTCH3, as a major oncogenic driver in T-cell ALL (T-ALL),³⁰ was identified as a SEG_high gene in multiple T-ALL cell lines (supplemental Figure 5A), and several SEG_high genes were identified as potential markers in CLL,³¹ for example, ABCA6, ARHGAP44, and WNT3 (supplemental Figure 5B). Notably, B-ALL and pre–B-ALL shared 39 of 40 (97.5%) SEG_high genes, which suggests high similarity in the molecular characteristics and leukemogenesis between these 2 leukemia subtypes. In addition, BLACE, which is identified as the SEG_high gene in B-A (supplemental Figure 5C), was reported to be expressed exclusively in B-ALL but rarely expressed in other types of hematologic malignancies.³² Note that additional SEGs were included in LeukemiaDB for further study.

Next, we explored specific circRNAs in leukemia subtypes or cell lines; however, the number of circRNAs varied greatly (Figure 4C-D). Nonetheless, some specific circRNAs have been reported to play important roles in leukemogenesis. For example, several T-ALL–specific circRNAs identified in our study, for example, TASP1- and PSEN1-derived circRNAs, have been reported to be dysregulated in T-ALL.³³ Notably, TASP1-derived circRNA was identified as a specific circRNA in T-ALL cell lines, for example, CUTLL1 and HPB-ALL. In addition, the PVT1-derived circRNA, which was previously reported to affect ALL cell growth and pose a potential therapeutic target,³⁴ was detected in T-ALL.

Regulatory network analysis of TF-target, mRNA-lncRNA, and mRNA-circRNA was performed to better elucidate the regulatory relationship of protein-coding genes, lncRNAs, circRNAs, and TFs. In patient samples, 2 297 453 to 21 289 406 highly correlated mRNA-mRNA coexpression pairs were detected (supplemental Figure 5D), and 3 994 844 to 37 707 656 pairs were found in cell lines (supplemental Figure 5E). In addition, we highlighted TF-target pairs with high correlation (|R| > .5; P value < .05) and kept 19 687 to 167 951 pairs in patients with leukemia and 43 876 to 325 966 pairs in cell lines (Figure 5A-B). The number of mRNA-lncRNA pairs varied from 1 400 200 to 71 435 444 in patients with leukemia (supplemental Figure 5F) and 243 373 to 277 991 449 in cell lines (supplemental Figure 5G). Notably, we identified fewer mRNA-circRNA pairs than other coexpression types, ranging from 0 to 234 in patients and 0 to 216 in cell lines (Figure 5C-D).

Previous studies have demonstrated the critical roles of mRNA-circRNA, mRNA-lncRNA, and TF-target pairs in leukemia progression.^35,36 In this study, we observed that the circRNA derived from DNAH14 was highly correlated to the host gene in T-ALL (R = 0.73; P value = .02), which implies the regulation by the circRNA of its host gene DNAH14 (supplemental Figure 5H). In addition, a previous study showed that the circular isoform of DNAH14 may be involved in the regulation of tumor-associated pathways,³⁷ which may function by regulating the host gene. We also inspected the potential roles of coregulation modules in leukemia, and we found that several HOXA9-related pairs formed a HOXA9 module (Figure 5E), in which the HOXA9/MEIS1 pair was frequently coexpressed in human myeloid leukemia³⁸ and related to leukemogenesis.^39,40 Another study reported that PBX3, which is highly correlated with HOXA9 (R = 0.68; P value < 2.0 × 10¹⁶), is an important cofactor of HOXA9 in AML.⁴¹ In this study, we observed that a crucial driver gene, that is, NOTCH1 of T-ALL, was coexpressed with multiple lncRNAs, for example, NONHSAG114941.1 (R = 0.61; P value = 6.2 × 10⁹) and NONHSAG103232.2 (R = 0.50; P value = 6.6 × 10⁶), which has been reported to promote cell proliferation in T-ALL by interacting with NOTCH1.⁴² Note that the developed LeukemiaDB includes additional modules for the research community.

To better elucidate the characteristics of AS events among patients with leukemia or leukemia cell lines, we first detected the aforementioned AS events, of which SE AS events with a high splicing ratio (>0.9) accounted for >50% (Figure 6A,C). In addition, specific AS events were identified, and the number of specific AS events varied from 7 (myelodysplastic syndrome) to 4551 (AML) in patients (Figure 6B) and from 1 (eg, Jurkat and MOLM-1) to 3293 (HL-60/S4) in cell lines (Figure 6D). The specific AS events may contribute to the progression of specific leukemia subtypes, for example, we found that the SE event of SLC38A9 was AML-specific and occurred in 49% of AML samples, which has been reported as a risk factor in AML.⁴³ 

We also found that the AS of MYB and HNRNPH1 occurred frequently in all leukemia subtypes and several cell lines (Figure 6E; supplemental Figure 5A), which implies that the dysfunction of AS in these genes may play important roles or function as valuable biomarkers in leukemogenesis. MYB is considered a potential therapeutic target in leukemia,⁴⁴ and HNRNPH1 has been reported as a crucial regulator of cellular proliferation and disease progression in CML.⁴⁵ In addition, we identified the top 10 genes with the most AS events across leukemia subtypes. Here, we found that the AS of several genes appeared in most leukemia subtypes, for example, A3SS of MYL6 (supplemental Figure 6B), A5SS of MYB (supplemental Figure 6C), RI of HMBS (supplemental Figure 6D), mutually exclusive exons of GTPBP10 (supplemental Figure 6E), and SE of POLL (supplemental Figure 6F). These top genes may contribute to the progression and treatment of leukemia.

Gene fusion, as a main driver event in leukemia, can serve as an important biomarker in the diagnostic therapy and prognosis of leukemia, for example, BCR-ABL in CML.⁴⁶ LeukemiaDB includes deposited comprehensive fusion events from patients with various leukemia subtypes and corresponding cell lines. Most of the fusion events mentioned in a previous study⁴⁷ can be found in our LeukemiaDB, both the common (BCR-ABL1, PML-RARA) and rare (or defined as novel: EBF1-PDGFRB) fusion events. Further analysis of fusion events demonstrated that 10 of 196 fusion events in all the patients with leukemia were shared among different leukemia subtypes, and 493 of 1694 in all cell lines samples were shared among different cell lines, whereas other fusion events are specific in different leukemia subtypes (Figure 7A-B). The occurrence rates of fusion genes ranged from 0.15 to 0.71 in patients with leukemia (supplemental Figure 7A) and 0.41 to 1.00 among cell lines (supplemental Figure 7B). Notably, we found that several specific fusion genes play important roles in the leukemogenesis of given leukemia subtypes. For example, as a driven event in CML,⁴⁶ BCR-ABL1 was identified as a specific fusion gene in patients with CML (Figure 7C; occurrence rate = 0.71) and CML cell lines (Figure 7D; occurrence rate = 0.65-1.00). In addition, the fusion score of BCR-ABL1 in 10 of 14 patients with CML was >0.74 (Figure 7E), and its fusion score in most CML cell line samples was >0.90 (Figure 7F). The fusion gene KMT2A-AFF1, which was previously reported to be related to poor prognosis of pro–B-ALL⁴⁸ and occurred in subtype 4 or subtype 5 AML, was detected in 22 of 28 (78.6%) of samples of SEM (pro–B-ALL cell line) and 9 of 16 (56.3%) of samples in MV4-11 (subtype 5 AML cell line; Figure 7G). Notably, the average fusion scores of KMT2A-AFF1 were >0.9 in SEM and MV4-11 cell lines (Figure 7H-I).

Leukemia is the main group of hematopoietic malignancies with high histological and functional heterogeneity. Over the past decade, many studies have focused on different leukemia subtypes, resulting in vast amounts of transcriptome data. However, these data are dispersed among thousands of published studies and cannot be shared directly or repurposed in the research community because of technical barriers, such as batch effects and bioinformatic skills. Therefore, in this study, we systematically integrated and analyzed large-scale high-quality RNA-seq data of human leukemia and the information regarding the influence of genes on drug response and survival, and then we developed the comprehensive leukemia database LeukemiaDB, which allows users to investigate the important features of RNAs underlying different leukemia subtypes at the transcriptional level.

A systematic investigation of expression profiles of protein-coding genes and noncoding RNAs and the regulatory mechanisms across leukemia subtypes can help us better understand leukemogenesis and progression, thereby facilitating effective leukemia diagnosis and therapy. In this study, we detected ∼20 000 protein coding genes, 60 000 lncRNAs, 8882 circRNAs, 5 AS event types, and various fusion genes (Figure 2). Our data have demonstrated that different leukemia subtypes or cell lines exhibit both specific and common characteristics. For example, aCML is a rare subtype of myelodysplastic/myeloproliferative neoplasm that shares clinical and laboratory features with CML.⁴⁹ The ARH77 cell line is derived from a patient with plasma cell leukemia; however, it exhibits the characteristics of a healthy B-lymphoblastoid and those of Epstein-Barr virus–transformed B-lymphoblastoid cell lines.^50,51 These particularities may cause aCML, and the ARH77 cell line shows great differences from patients with other leukemias and other leukemia cell lines. Note that the detected genes or AS events are dependent on the sequencing depth of the sample. For example, we observed that more fusion genes and AS events were identified in the HL-60/S4 cell line (Figure 2I; supplemental Figure 4I), which may be caused by samples with the highest sequencing depth. In addition, our results highlight the importance of lncRNAs, circRNAs, and TFs in the gene dysregulation and pathogenesis of leukemia, for example, HOXA9 in AML (Figure 5E) and NONHSAG114941.1, NONHSAG103232.2 (Figure 5F), and circDNAH14 in T-ALL. Massive regulatory pairs, including both reported and newly identified pairs, are provided to the research community in the developed LeukemiaDB data portal.

The structural and splicing alteration of RNAs as oncogenic drivers play important roles in leukemogenesis and progression.^43,46 Our integrated analysis of large-scale RNA-seq data provides comprehensive views for gene fusion and AS events in leukemia subtypes and cell lines. For example, fusion gene BCR-ABL1, as the driven event and therapy target in CML,⁴⁶ was detected in 96% of the CML samples (Figure 5C). In addition, several specific fusion events with high occurrence rates and fusion scores were detected in pro–B-ALL cell lines, for example, KMT2A-AFF1 (Figure 5F). Furthermore, other recenty discovered fusion genes are included in LeukemiaDB, which may provide new insights into the mechanism of leukemogenesis in different leukemia subtypes.

In this study, we applied rigorous strategies to curate and integrate leukemia RNA-seq data from different studies; however, several limitations can be identified and discussed. Firstly, most of the considered data sets were generated using different library protocols and sequencing platforms, thereby representing technological limitations. The severe inhomogeneity of sequencing data can influence the detection sensitivity of RNAs and downstream analysis, for example, inefficient detection of circRNAs on poly-A–enriched RNA-seq data and the influence of read length on fusion gene detection.⁵² Secondly, even for the same leukemia subtype or cell line, several factors, for example, different treatment conditions, genetic characteristics, and unavailable clinical features (representing biological limitations), can affect data integration from different studies. In some cases, batch effects, which may be the largest barrier to integration analysis for data from different studies, often result in discrepancies in the statistical distributions across data from different studies or laboratories. Although we used several well-known algorithms and tools to perform BER and different leukemia subtypes were grouped together in the UMAP plot, several batch effects could be observed as well, such as the CLL and CML (supplemental Figure 3). Moreover, we also observed that several samples were confounded with other subtypes of samples, for example, a few of T-ALL samples in CLL cluster (supplemental Figure 3) and pre–B-ALL samples mixed with B-ALL, which may be attributes of those samples (we rechecked their metadata) and a limitation of available BER methods.⁵³ Finally, the unbalanced sample sizes of different leukemia subtypes may influence the power of statistical analysis. For example, massive circRNAs were detected in K562 (supplemental Figure 4D), which may be caused by the large amount of public RNA-seq data for the K562 line. Besides, although we had carefully controlled batch effects during data integration from different studies, challenges remain in batch effects that require further investigation. In addition, because of the incomplete clinical metadata of public data sets, samples had not been further subdivided by cytogenetics, morphology, and immunophenotype profiles, which is a chief limitation of integration analysis from public data and further direct application of our data on precision medicine in Leukemia researchers. However, further use of these data is severely hampered by a lack of integrative tools for visualization and analysis, and this integrated resource in our study is the optimal solution.

Collectively, we have provided comprehensive and multidimensional transcriptome profiles (protein coding genes, lncRNAs, circRNAs, AS events, and fusion genes) of patients with leukemia and leukemia cell lines. In the future, we will plan to update the LeukemiaDB data portal along with the massive cohorts released, and further analysis will be performed to solve urgent problems in the field of leukemia research.

The authors are grateful to the laboratory members and the providers of public data for their contribution.

This work was supported by the National Key R&D Program of China (2021YFF0703704), the Natural Science Foundation for Distinguished Young Scholars of Hubei Province of China (2020CFA070), the Science, Technology and Innovation Commission of Shenzhen Municipality (JCYJ20210324141814037), the China Postdoctoral Science Foundation (2021M701793, 2021M824702, and 2022M713541), and the Science and Technology Project of Nantong City (MS22022109).

Contribution: A.-Y.G. and Q.Z. designed this project; M.L. and Q.Z. collected data and performed analysis; Y.J.K. worked on the fusion identification for leukemia subtypes; M.L. and Y.M. constructed the database; and A.-Y.G., M.L., and Q.Z. wrote and revised the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Qiong Zhang, Research Center of Clinical Medicine, Affiliated Hospital of Nantong University, Nantong 226001, China; e-mail: transomics@gmail.com; and An-Yuan Guo, College of Life Science and Technology, Huazhong University of Science and Technology, Wuhan 430074, China; e-mail: guoay@hust.edu.cn.