To construct our reference healthy human BM atlas, we used 145 publicly available human BM scRNA-Seq datasets with publicly available sequencing reads from 22 projects ( Figure 1A ). Of the 145 samples, 92 samples have age information, with age ranging from 2 to 84 years and a median of 45 years ( Supplementary Table 1 ). Our collected data also includes 34 foetal BM samples. From samples with gender information, the gender ratio is relatively balanced at 54 female and 45 male samples. 91 samples were sorted using a variety of strategies. All samples were sequenced using the 10x Genomics sequencing platforms.

The raw reads were previously processed as part of the DISCO database (9). We employed our pipeline to remap the reads to a single reference genome, GRCh38 (Ensembl 93), for annotation consistency and to reduce potential technical effects. We then processed the resulting read counts using the Seurat package (10). Standard quality control steps on UMI counts, number of detected genes, and the fraction of mitochondrial and ribosomal genes were applied to filter out low quality cells. To remove the significant batch effects present, we employed our FastIntegration tool developed for atlas-scale integration (11) ( Supplementary Figure 1 ). This was followed by unsupervised clustering and differential gene expression analysis. We then annotated each cluster’s cell type using canonical marker genes ( Supplementary Figure 2 and Supplementary Table 2 ).

Our constructed healthy BM atlas is currently the largest among existing healthy BM atlases (7, 12, 13) in both cell count and number of donor samples. It consists of 673,750 cells with 54 annotated cell types ( Figure 1B ). The cell types present mostly overlap with the previously constructed healthy BM atlas by Hay et al. (7) but are annotated at higher resolution of cell subtypes with known markers. Our annotated cell types can be divided into five major groups: T/NK cells, B cells, monocytes and DCs, progenitor cells, and erythrocytes/megakaryocytes. Other smaller clusters of cells identified are osteoclasts, fibroblast, and mesenchymal stromal cells. For the T/NK cells, we could divide the naïve population into CD4⁺ and CD8⁺ subtypes, the CD8⁺ T cells into GZMK and GZMB subtypes, and NK cells into CD16⁺ and CD56⁺ subtypes. We also identified mucosal-associated invariant T (MAIT), gamma-delta T, and regulatory T (Treg) cells. For B cells, we could identify the subtypes along its developmental path from the common lymphoid progenitors (CLP) to pro-B, pre-B, naïve B, memory B, and plasma cells. The monocyte populations found were the CD14⁺ and CD16⁺ subtypes, while the dendritic cells were divided into classical dendritic cells (cDCs) and plasmacytoid dendritic cells (pDCs). We compared our annotations with available annotations of the contributing samples (GSE185381 (14) and Census of Immune Cells (8)) and we found good concordance among the major cell types ( Supplementary Figure 3 ).

Using the unsorted samples, we computed the distributions of cell type proportions ( Figure 2A ). Mature T cells formed the majority of cells in the bone marrow with naïve CD4⁺ T cells being the most numerous. Naive CD8⁺ T cells expectedly showed a much lower average percentage with this skew being attributed to the high failure rate of CD8⁺ T cells during selection in the thymus (15). Other cell types, namely monocytes, B cells, and NK cells, made up significant fractions. Overall, these proportions are in line with previous studies on the cell types present in the bone marrow (13).

By integrating a large number of cells into a single atlas, there are sufficient cells to populate the secondary developmental pathways. Here our atlas illustrates the different developmental pathways of pDCs (marked by IRF7), namely myeloid and lymphoid origins (16) ( Figure 2B ). The myeloid origin pathway (blue) traces from the DC precursor populations and is marked by the myeloid-associated LYZ. The lymphoid origin pathway (green) traced from the CLP population and is marked by JCHAIN and HIST1H4C. Within the T cell population, we could also distinguish the Treg cell subset ( Figure 2C ). The large cell count enabled us to increase the resolution of clustering and identify rare cell types in the bone marrow, such as osteoclasts ( Figure 2D ) and mesenchymal stromal cells. These cells were not annotated in previous healthy bone marrow atlases. As the integration included foetal data samples, we could identify foetal erythrocytes expressing HBG1 and HBG2 (foetal HBG⁺ erythrocyte).

In ageing studies, blood cell population changes are typically analysed using peripheral blood samples. We first examined cell type frequency changes in the bone marrow and compared them against those reported in human blood samples. Due to large variances in cell populations across samples and studies, we only used samples from three studies, namely GSE120221 with 25 samples (13), GSE185381 with 10 samples (14), and the Census of Immune Cells dataset with 8 samples (8). We selected these studies as they had large numbers of unsorted samples with age information. The remaining studies were either composed of sorted samples or had too few donor samples (<5).

We tested the correlations for each annotated cell type using a linear model with the study batch as covariate. Two cell types, effector/memory CD4⁺ T and CD16⁺ NK cells, had statistically significant positive correlations across all samples ( Figure 3A , Supplementary Table 3 ). The effector/memory CD4⁺ T cell’s correlation was statistically significant when combining the regression output from all three studies and even within the GSE120221 and Census study sets, giving us the greatest confidence in this result. Moreover, this correlation has also been found in peripheral blood by Li et al. (17). Increasing proportions of CD16⁺ NK cells have also been reported in the blood of elderly individuals (18).

Among the annotated cell types negatively correlated with age, the naïve CD8⁺ T cell population had the largest average coefficient and was consistent across the three studies. Decline in naïve T cell populations have been reported for peripheral blood in different studies (17, 19, 20) and is well recognised. For the corresponding naïve CD4⁺ T cell population, we obtained a negative correlation with the GSE120221 and GSE185381 sets, but a positive correlation with the Census study set ( Supplementary Figure 4 ). The latter result is unexpected and can be explained by the small number of samples (n=8) coupled with high variability within the data that can mask the expected age-related decline. Another notable cell type negatively correlated with age is the CLP population. As progenitor cell types, changes in these populations have disproportionate impact on mature cell type populations and consequently the immune system’s functions. This decline in the proportion of lymphoid-biased progenitors and consequent myeloid bias in haematopoiesis is widely reported (21). The computed correlations for all cell types are given in Supplementary Table 3 .

We next investigated the cell type-specific genes with expression that correlate with age, using all samples with age metadata ( Figure 3B ). Among the annotated cell types, GMP cells had the largest number of age-associated differentially expressed genes (DEGs) with 663 genes, while IFN-activated T cells had the least with 107 genes. Most of these genes were shared among few cell types, with only 20 genes shared in 10 or more cell types. The most conserved DEG was RPL13A in 17 cell types ( Figure 3C ). RPL13A’s upregulation with age has also been reported for almost all tissue types of mice (22). Other genes of note include IMMP2L and DIP2B which showed downregulation with age in 10 cell types. Switching off IMMP2L signalling has been shown to drive cell senescence (23). Similarly, DIP2B knockout cells expressed senescence markers found in ageing cells (24). We also found positive age correlations for chemokines CCL4, CCL4L2, and CCL5 in tissue-resident NK cells, CCL4L2 in CD16⁺ NK cells, CCL4 in CD56⁺ NK cells, and CCL3 in GZMK⁺ CD8⁺ T cells ( Supplementary Table 4 ).

With ageing, the skew towards myeloid lineage haematopoiesis is well documented (25). Within the MPP population, we investigated transcriptomic changes that can contribute to the decline in differentiation towards the lymphoid lineage. Here we found two notable genes, GADD45B which was upregulated with respect to age, and ZBTB16 that was downregulated ( Figure 3D ). GADD45B has been characterised as a myeloid differentiation gene and plays a role in the response of myeloid cells to stress stimulation (26, 27), while the ZBTB16 gene is a negative regulator of myeloid cell development (28).

We followed up with pathway analysis of the age correlated genes using EnrichR (29). Here we show the enriched Reactome pathways ( Figure 3E ). The different cell types show a wide range of dysregulated pathways. We note that the CD14⁺ MHCII^low monocytes, late haemoglobin⁺ erythrocytes, memory B cells, and pre-B cells similarly show dysregulated translation processes. Other B cell subtypes, namely pro-B cells and naïve B cells, show dysregulated mRNA processing pathways, and the naïve B cell subset also show dysregulated cellular stress pathways. For the interferon-activated T cells, the dysregulated pathways are primarily centred around metabolism and interferon signalling. The metabolic changes in T cells due to ageing is linked to cell senescence and reduced functionality, though the exact mechanisms are under investigation (30).

Cancer is generally recognised as an ageing-related disease. While ageing increases the risk of cancer, it has also been reported that cancer can alter the expression trends of ageing-related genes (31). We theorised that diseases including cancer can modify the gene expression of diseased cells to appear younger or older than their chronological age. Thus, we investigated how diseased cells’ age predicted by their transcriptome differ from their chronological age. We first trained a cell age predictor using the healthy atlas. The predictor was constructed using elastic net regression and the age-related DEGs identified for each cell type. For each sample, we computed the age of each cell type present, and the median predicted age was used as the overall predicted age. The feature gene set were selected via regularisation in the elastic net regression, and we employed 10-fold cross validation to check the hyperparameters ( Supplementary Figure 5 ). The resulting predictor was able to predict the ages of healthy individuals with a correlation coefficient of 0.92 and an error of 6.3 years ( Figure 4A ). We further tested our predictor using the corresponding cell types from healthy blood samples, obtaining good predictions with CD14⁺ monocytes and Treg cells ( Figure 4B ). List of blood samples are presented in Supplementary Table 5 .

We then used the age predictor to predict the age of patients diagnosed with acute myeloid leukaemia (AML) and multiple myeloma (MM). The predicted ages of the cancer samples showed the trends of older patients being predicted to have a younger age while younger patients were predicted to have an older age ( Figure 4C and Supplementary Table 6 ). For the AML samples, the predicted ages fall in the range of 32 to 51 years old with an average of 41.8. For the MM samples, the predicted ages fall in the range of 49 to 60 years old with an average of 53. This suggests that cancer alters gene expression patterns to portray a cancer-specific apparent age that may be different across different cancer types. In the case of AML, this also affects the paediatric disease subtypes. Future investigation into this apparent cell age phenomenon will encompass other blood cancers.

To demonstrate our healthy BM atlas’ applicability as a reference for identifying diseased cells, we mapped samples of AML and MM onto our atlas. By integrating healthy and diseased samples together, abnormal cells can be identified. We first mapped 10 MM samples from GSE189460 onto our healthy atlas and performed label transfer to annotate the cells ( Figure 5A ). We also computed a cell type prediction confidence score based on each cell’s distance to its neighbours in the reference atlas (Methods). Lower prediction confidence scores denote mapped cells that were phenotypically different from cells found in the healthy atlas, which in this case implied diseased cells. Among the mapped cells, the lowest scores were found among the predicted plasma and erythroblast cells ( Figure 5B ). The predicted plasma cell types correspond to the malignant plasma cells that accumulate in the bone marrow while abnormalities in the erythroblast compartment correspond to the disrupted erythropoiesis process which gives rise to the common MM symptom of anaemia (32).

We compared the malignant plasma cells to their healthy counterparts and obtained 80 DEGs (p< 0.05 and |logFC| > 0.5) ( Figure 5C and Supplementary Table 7 ). The upregulated genes include the previously reported JUND (33), FOSB, PDIA2 (34), and CCL3 (35). Elevated levels of CCL3 in the bone marrow is of particular interest as it has been proposed to suppress erythropoiesis and cause anaemia in MM patients (36). In the list of downregulated DEGs, low CD27 expression in malignant plasma cells has been reported to be correlated with poor prognosis (37) while the tumour suppressor gene WWOX has been found downregulated due to translocations or deletions (38). In the enrichment analysis, dysregulated pathways include TP53 and cytokine signalling ( Figure 5D ). The dysregulation of cytokine signalling is an expected characteristic of malignant plasma cells due to the high levels of cytokines that also serve as chemoattractants that attracts them to the bone marrow (39).

We next mapped 10 AML samples from GSE185381 (14) onto our BM atlas ( Figure 5E ). We performed label transfer to label the AML sample cells and computed the prediction confidence score for each cell. Comparison of the computed scores with the malignant cell labels in the original annotation showed good correlation with the malignant cells have low scores ( Figure 5F ). The malignant cells’ transferred labels are primarily progenitor types and myeloid types (monocytes and dendritic cells), which have been previously reported (40). This presents an approach towards identifying malignant cells as an alternative to employing mutation detection.

We retrieved healthy bone marrow sample datasets from the DISCO database (9). The datasets were previously preprocessed from raw reads and mapped onto the human reference genome, GRCh38 (Ensembl 93). Except for data integration with FastIntegration, we employed the Seurat package (10) for downstream data analyses. For each sample, we filtered the cells based on their unique molecular identifier (UMI) counts, the number of detected genes, and the fraction of mitochondrial and ribosomal genes. As the data were acquired in different experiments, we utilised the distribution of QC metrics from each sample to manually determine the cut-offs applied for filtering. Subsequently, samples with less than 200 cells were removed and the gene expression for each cell was normalised to the total expression using the “NormalizeData” function found in the Seurat package.

We next applied our FastIntegration algorithm to integrate the retained data samples. For each sample, we first identified the top 3000 highly variable features. Thereafter, we merged the lists from all samples and selected the top 3000 most common ones. These highly variable features were used to identify the anchors between samples and the anchors were then used for the subsequent integration and batch correction steps. After integration, the batch-corrected gene expression values were standardised using the “ScaleData” Seurat function and Principal Component Analysis (PCA) was performed. The first 30 PCs were then used for Uniform Manifold Approximation and Projection (UMAP) to enable visualisation. For clustering, we built the KNN graph based on the Euclidean distance in PCA space and applied the Louvain algorithm. The Wilcoxon rank sum tests were used for identifying differentially expressed genes (DEGs) in each cluster, which were subsequently utilised for manual cell type annotation.

We selected three projects (GSE185381, GSE120221, and HCA_HematopoieticImmuneCellAtlas) that had unsorted data and sample sizes larger than 5. For each cell type, we performed linear regression on the proportion of each cell type with age with the following model:

where x_j is the percentage in sample j. A Spearman correlation p value was computed for each project and an overall p value for the linear model.

To identify age-correlated genes, we averaged the batch-corrected gene expression values of each cell type in each sample. We only considered cell types found in at least 20 samples and with at least 20 cells. Only the genes expressed in > 10% cells were used. Subsequently, we correlated gene expression with age using the Pearson correlation and retaining genes with p value less than 0.01 as age-related genes. For gene set enrichment analysis, we used Enrichr (29) with gene sets from the KEGG, GO, Reactome, and TRUST databases. Pathways with an adjusted p value less than 0.01 were selected for visualisation. We also compared our identified age-related genes with the gene list in GenAge's database (42) and genes identified in epigenetic clocks (43).

Based on the identified age-related genes, we trained a model to predict the age of a sample. As the age-related genes we identified were linearly correlated with age, we applied elastic net linear regression from the glmnet package (44) to build an age prediction model for each cell type. The performance of this model was evaluated using 10-fold cross validation ( Supplementary Figure 5 ), and the median predicted age of all cell types was taken as the predicted age for the sample. To estimate the model’s accuracy, we calculated the median absolute difference between the predicted age and chronological age. Prior to predicting the age of new cell samples, we first mapped them to the reference atlas and using only the cells that were confidently mapped for the age prediction.

We downloaded single-cell RNA-Seq data of acute myeloid leukaemia and multiple myeloma samples from the DISCO database. For each disease, we first integrated the data of all samples using FastIntegration. Subsequently, another round of integration was performed to integrate the diseased samples together with the healthy bone marrow atlas with the latter serving as the reference. Finally, we used the integrated data for PCA, UMAP generation, and clustering. We annotated the disease sample cells based on the most prevalent cell type among its 30 nearest neighbours in the reference atlas. We then compute a confidence score by taking the inverse of the distance between each cell and its 30 nearest neighbours, and then normalising it to a range of 0 to 1. The confidence scores were fit into a two-component Gaussian mixture distribution, with cells in the first component being deemed as confidently assigned.