Integrating RNA-seq and scRNA-seq to explore the

IntroductionAcute lymphoblastic leukemia (ALL) relates to the situation that immature lymphoid progenitor cells are expanded clonally and lymphocyte sources become abnormal (1). The American Cancer Society (ACS) reported 5,290 new adult and pediatric cases of ALL that resulted in 1,580 deaths in 2021 (2). Despite the obviously declined ALL incidence with time, ALL remains the primary tumor type for children (3). ALL is composed of the B-cell and T-cell lineages, with the former taking up about 85% of pediatric ALL (4). Chemotherapy and the Hematopoietic Stem Cell Transplantation (HSCT) technique have developed rapidly, as a result, pediatric ALL presents the cure rate of nearly 90% (5). Its overall survival has progressed remarkably in the long term, while 15–20% of patients still suffer ALL relapse, that largely explains the increase in mortality for ALL patients (6). Therefore, the identification of potential biomarkers of relapse in B-ALL patients is crucial to improve the prognosis.Nicotinamide adenine dinucleotide (NAD +) acts as a pivotal coenzyme in the redox reaction, and also crucially constitutes the energy metabolism. In traditional concept, nicotinamide nucleotide metabolism is very static, which mainly puts emphasis on the way different forms of NAD (oxidized or reduced) interconvert with the nicotinamide adenine dinucleotide phosphate(NADP) (7, 8). Nevertheless, according to researches in the past thirty years, NAD possesses complex and dynamic biological progresses, such as metabolism, transport, and function. NAD is capable of being converted into NADP, nicotinic acid adenine dinucleotide phosphate (NAADP), and cyclic ADP-ribose (cADPR), thereby remarkably impacting the energy transduction and the cell signaling. Besides, the degradation products of NAD, like nicotinamide and positive methyl nicotinamide, are also pivotal modulators of the epigenetics, energy metabolism, as well as disease states (9, 10). The pathway metabolites NAD can be the substrates for various enzymes like PARPs, which can impact the cellular homeostasis from various aspects (11, 12). What’s more, KPT-9274 acts as a different nicotinamide phosphoribosyl transferase (NAMPT) inhibitor, and NAMPT limits the NAD rescue biosynthesis rate, and crucially affects the energy metabolism. KPT-9274 could inhibit NAMPT for NAD + depletion, thereby inhibiting B-ALL cell growth (13). Despite this, researches fail to well explain the impact exerted by NAD + metabolism-related genes on the recurrence of B-ALL.Taken together, this study explored the biological significance of NAD + metabolism in the initial diagnosis and recurrence of B-ALL based on the differential NAD + metabolism-related genes in the initial diagnosis and recurrence of B-ALL patients, searched potential biomarkers of B-ALL, revealed the molecular mechanism at the single-cell level, and deepened the understanding of the pathogenesis of B-ALL.Materials and methodsData sourcesThe workflow chart of this study was shown in Supplementary Figure S1. Transcriptome data of B-ALL patients were selected from GSE3912, which contained Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) gene-expression matrix of bone marrow samples in 32 first diagnosis children patients and 54 relapses samples, and was sourced via GEO database. GSE130116, single-cell RNA-sequencing libraries of bone marrow from 7 pairs of pediatric B-ALL patients with initial diagnosis or relapse with gene expression matrix as FPKM value, was also acquired from GEO database.In addition, according to the research of Li et al., NAD+ metabolization-related genes (NMRGs) were selected from KEGG database (Pathway: hsa00760) and Reactome Database (R-HSA-196807) (14). After combination and duplication eradication, a total number of 51 NMRGs were obtained. Moreover, complete clinical and expression information of 9 pediatric B-ALL patients was acquired from the TARGET database. Clinical and Demographic characteristics of B-ALL patients from GSE3912 dataset, GSE130116 dataset, and TARGET database were exhibited in Supplementary Table S1.Screening of differentially expressed NMRGsBased on 51 obtained NMRGs, the expression levels of NMRGs were found in the GSE3912 dataset, and the differentially expressed NMRGs (DE-NMRGs) between the initial diagnosis and relapse samples were compared using R-package limma (Version 3.48.3), with |log2FC| > 0.5, p Value < 0.05 as screening criteria (15). Next, the screened DE-NMRGs were performed functional enrichment analyses by ClusterProfiler (version 4.0.2). The enriched results satisfied p < 0.05 and count ≥ 1 were regarded as significantly enriched.Identification of biomarkers and correlation with clinical characteristicsR package Boruta (version 7.0.0) was firstly used to determine the importance of each DE-NMRGs through random forest (RF) algorithm and select DE-NMRGs with the confirmed importance (16). Then, the repeated cross validation was applied to the selected DE-NMRGs to further screen the most important genes which were considered as biomarkers. In order to investigate the ability of biomarkers to extinguish initial diagnosis from relapse samples, logistic regression fitting was employed on the biomarkers, and receiver operating characteristic (ROC) curves for logistic regression model was plotted using R package pROC (Version 1.18.0) (17). The ggplot2 (version 3.3.5) was employed to visualize the expressions of biomarkers in initial diagnosis and relapse samples (18).Furthermore, to investigate the difference in biomarker expressions between different clinicopathological characteristics of B-ALL patients, biomarkers expression information and clinical characteristics (CNS Status, Ethnicity, Gender, Race) of the 9 samples from TARGET database were extracted and compared the expression levels of biomarkers in different clinical traits by ggplot2 (Version 3.3.5).Tumor microenvironment analysisTME cells were important components in tumor tissues, and an increasing number of evidence had proven their clinicopathological significance in prognosis and treatment of cancers. Cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) was first employed to compute the proportions of 22 types of immune cells in the samples from initially diagnosed and relapse groups, and the results were visualized by tidyverse (version 1.3.1) (19). Then, the percentage of each immune cell was compared between the first diagnosis and relapse group and plotted by vioplot (version 0.3.7) (https://github.com/TomKellyGenetics/vioplot) to distinguish the differentially expressed immune cells (DEIs) (20). Furthermore, corrplot (version 0.92) (https://github.com/taiyun/corrplot) was applied to detect the correlations between the 22 types of immune cells and the biomarkers, to find the immune cells significantly correlated with biomarkers (p < 0.05) (21).In addition, the overlap analysis was conducted between the DEIs and the immune cells correlated with biomarkers, and the intersection was regarded as key DEIs. Subsequently, based on the expression data of the biomarkers, the correlations between biomarkers and corresponding gene markers of key DEIs were determined by Pearson correlation analysis.Cluster and pseudotime analysisSeurat (version 4.1) in R package was employed to conduct scRNA-seq quality control on the 7 pairs of B-ALL samples in GSE130116, with nFeature_RNA > 100, percent.mt < 5, and nCount_RNA > 3 as screening criteria (22). Then, the filtered genes were further conducted variance analysis, and the top 2,000 genes whose expression varied significantly among cells were identified and used for subsequent cell type identification. The principal components for highly variable genes were calculated. Based on the PCA result, principal components with p values less than 0.05 were used to identify clusters using the UMAP2 algorithm. Finally, ‘FindAllMarkers’ function was applied to search cluster markers with min.pct=0.2 and only.pos=TRUE. Cluster cell types were assigned according to cluster markers and cluster labels from the R package ‘SingleR’ (version 1.6.1) (23). In addition, the pseudotime analysis was employed by R package ‘monocle 2’ (version 2.20.0), and the cells were visualized in trajectory (24).Ligand-receptor analysisThe cell communication analysis was conducted through CellPhone DB database. Firstly, the numbers of ligand-receptor interaction and polymer between cell types were counted, and the
Read More