Single-Cell RNA-seq Analysis of Human Bone Marrow

An end-to-end scRNA-seq analysis in R/Seurat of human bone marrow and CD34+ hematopoietic progenitor cells: quality control, doublet removal, batch correction, clustering, cell-type annotation, differential expression, pathway enrichment, trajectory inference, and cell-cell communication analysis.

View the full rendered analysis notebook →

Overview

The dataset comprises four 10x Genomics scRNA-seq samples of human bone marrow from Granja et al. 2019: two whole Bone Marrow Mononuclear Cell (BMMC) samples and two CD34+-enriched hematopoietic progenitor samples, spanning three donors.

Sample	Cell type	Donor	Cells (raw)
BMMC_D1T1	Bone Marrow Mononuclear Cells	D1	6,270
BMMC_D1T2	Bone Marrow Mononuclear Cells	D1	6,332
CD34_D2T1	CD34+ Enriched Bone Marrow	D2	2,424
CD34_D3T1	CD34+ Enriched Bone Marrow	D3	5,752

Each sample was profiled across 20,287 genes. After QC filtering and doublet removal, 18,693 cells remained across all four samples.

Pipeline

1. QC & filtering — per-sample thresholds on UMI count, gene count, and mitochondrial percentage, set using median ± k·MAD rather than fixed cutoffs, to adapt to each sample's own distribution.
2. Doublet removal (DoubletFinder) — pK swept per sample via its BCmetric, using the PCA elbow point of each sample (11–19 PCs) as dims.
3. Normalization & feature selection — log-normalization and highly variable gene selection (Seurat defaults), run independently per sample before merging.
4. Batch correction — two parallel merges were built and compared: a naive merge with no correction, and Seurat's anchor-based integration (FindIntegrationAnchors / IntegrateData) across all four samples. Sample-of-origin structure that dominated the naive-merge UMAP was resolved by integration, motivating batch correction as necessary here.
5. Dimensionality reduction & clustering — PCA (30 PCs, chosen from the elbow plot) → UMAP → graph-based (Louvain) clustering.
6. Cell-type annotation — automatic labels via SingleR against the HumanPrimaryCellAtlasData reference (celldex), cross-checked against a manual, marker-driven annotation (module scores over canonical marker sets, e.g. CD3D for T cells, CD19/CD20 for B cells, CD34/CD38 for HSCs). The two approaches were compared with a confusion matrix.
7. Differential expression — pairwise cell-type comparisons (e.g. B vs T cells, T cells vs Monocytes) and BMMC-vs-CD34+ comparisons, visualized as volcano plots.
8. Pathway enrichment — GO Biological Process enrichment (enrichR) on the BMMC-vs-CD34+ DE gene set.
9. Trajectory analysis (Monocle3) — pseudotime ordering of Common Lymphoid Progenitors (CLP), with manually- and automatically-selected root nodes compared.
10. Cell-cell communication (CellChat) — signaling pathway inference within BMMC and within CD34+ samples, restricted to cell types shared by both groups.

Results

16 cell types were resolved across the dataset — B cells, Plasma cells, Pre-B, T cells (CD4+/CD8+), NK cells, CD14+/CD16+ Monocytes, cDC, pDC, Basophils, Erythrocytes, and the stem/progenitor compartment (HSC, LMPP, CLP, GMP/Neutrophils).

Quality control Per-cell UMI counts across the four samples before filtering — the basis for the median/MAD-derived QC thresholds.	Batch correction UMAP of all four samples after Seurat's anchor-based integration, colored by unsupervised cluster (23 clusters at this resolution, later collapsed to 16 cell types).
Cell-type annotation Final manual cell-type annotation, based on differential expression plus a panel of canonical lineage marker genes.	Annotation agreement Confusion matrix comparing the manual, marker-based annotation against SingleR's automatic labels — strong agreement on well-separated types (T cells, CLP, CD16+ Monocytes), more disagreement in the progenitor compartment where marker signal is weaker.
Sample composition Cell-type proportions per sample — BMMC samples are dominated by mature lymphoid/myeloid cells (T, B, CLP, Monocytes), while CD34+ samples are enriched for stem/progenitor populations (HSC, GMP/Neutrophils, LMPP), consistent with CD34+ selection.	Differential expression Volcano plot for B cells vs. T cells — thousands of genes pass the significance threshold given the strong transcriptional separation between lymphoid lineages.
Trajectory analysis Monocle3 pseudotime trajectory restricted to Common Lymphoid Progenitors, rooted at the least-differentiated end of the population.	Cell-cell communication CellChat circle plots showing inferred outgoing signaling strength from selected progenitor and myeloid cell types toward the rest of the cell-type repertoire.

What the data shows

The annotated populations are consistent with active hematopoiesis. Common Lymphoid Progenitors are the single largest population (17.4%, 3,249 cells), followed by Granulocyte-Monocyte Progenitors (12.9%, 2,419 cells). The progenitor compartment (CLP, GMP, HSC, LMPP) together with a full complement of mature lymphoid and myeloid cells is what one expects from bone marrow rather than peripheral blood.

Transcriptional output tracks cell function. Erythrocytes carry the highest median UMI count (4,358) and the most genes detected per cell (1,905), consistent with the heavy globin transcription of erythropoiesis. CLP cells distribute across 19 separate clusters — more than any other type — indicating substantial heterogeneity in maturation state, whereas CD16+ Monocytes occupy only 9, a comparatively homogeneous population.

CD34+ selection is visible throughout. BMMC samples are dominated by mature lymphoid and myeloid cells, while CD34+ samples are enriched for HSC, GMP/Neutrophil and LMPP progenitors (see the proportions plot above). The same signal appears in the enrichment analysis: the top GO Biological Process term separating the two groups is neutrophil degranulation (GO:0043312) — a mature granulocyte effector program that CD34+-enriched samples largely lack.

Pseudotime recovers B-lymphopoiesis. Rooting the CLP trajectory on an IL7R-high cluster produces an ordering along which CD19 and MS4A1 rise together, from marker-negative progenitors through intermediate states to mature B cells. Because those two markers were not used to build the trajectory, their monotonic increase along pseudotime is independent support that the ordering reflects real differentiation rather than an artifact of the embedding.

Cell-cell signalling separates the two compartments sharply. Restricting to the 16 cell types present in both groups, BMMC supports roughly three times the inferred interactions of CD34+ and four times the total interaction strength:

	BMMC	CD34+
Cells	11,343	7,350
Significant pathways	23	12
Inferred interactions	2,061	622
Total interaction strength	111.7	28.5

Ten pathways are shared. Thirteen are BMMC-specific — including MHC-I, CD45, LCK, CD22, CD23 and ICAM, essentially the lymphocyte antigen-receptor and adhesion machinery of a mature immune compartment. Only two are CD34+-specific, NEGR and CDH, both adhesion families, consistent with progenitors engaging a stromal niche rather than each other.

MIF is also the pathway visualised in detail in the notebook, chosen precisely because it is significant in both groups — a comparison the analysis could not make while it was plotting WNT, which appears prominently in the CellChat reference database but is not inferred as significant in either group here.

The per-pathway strengths sharpen the contrast. MHC-II signalling is roughly nine-fold stronger in BMMC (13.1 vs 1.5), as expected when professional antigen-presenting cells are depleted by CD34+ selection. The one pathway that runs against the trend is midkine (MK), which is stronger in CD34+ than in BMMC (5.2 vs 2.3) — notable because midkine is a growth factor associated with progenitor proliferation, so its relative prominence in the progenitor-enriched samples is biologically coherent rather than incidental.

Repository structure

├── analysis/
│   ├── scrna_seq_analysis.Rmd   # full annotated analysis source
│   ├── report.html              # rendered notebook — start here
│   └── sessionInfo.txt          # exact R + package versions used
├── data/                        # input Seurat objects + celldex reference (Git LFS)
├── envs/
│   ├── environment.yml          # loose conda spec
│   ├── environment.lock.yaml    # pinned conda packages (linux-64)
│   └── r-packages-verified.csv  # manifest of a known-working environment
└── figures/                     # figures used in this README

Reproducing

git clone <this-repo>
cd scrna-seq-bone-marrow-analysis
git lfs pull                      # fetch the data/*.rds inputs

conda env create -f envs/environment.lock.yaml -n single-cell
conda activate single-cell

The conda spec pins R 4.0.5, Seurat 4.0.1 and Matrix 1.4-1, but leaves several transitive dependencies unpinned — and their current versions no longer work on R 4.0.5. Three groups need explicit pins:

# spatstat satellites must match spatstat.core 2.4-4 (2022 generation);
# newer geom moved `integral` to spatstat.univar, which breaks Seurat's load.
remotes::install_version("spatstat.utils",  "2.3-1", dependencies = FALSE)
remotes::install_version("spatstat.data",   "2.2-0", dependencies = FALSE)
remotes::install_version("spatstat.sparse", "2.1-1", dependencies = FALSE)
remotes::install_version("spatstat.geom",   "2.4-0", dependencies = FALSE)
remotes::install_version("spatstat.random", "2.2-0", dependencies = FALSE)

# statnet family: later releases use R 4.1 lambda syntax that R 4.0.5 cannot parse.
remotes::install_version("statnet.common", "4.5.0",  dependencies = FALSE)
remotes::install_version("network",        "1.17.1", dependencies = FALSE)
remotes::install_version("sna",            "2.6",    dependencies = FALSE)
remotes::install_version("ggnetwork",      "0.5.10", dependencies = FALSE)

# CellChat is not on conda-forge/bioconda.
install.packages("cleanrmd")
remotes::install_github("sqjin/CellChat", upgrade = "never", dependencies = FALSE)

The complete 384-package manifest of a verified-working environment is in envs/r-packages-verified.csv.

Then knit analysis/scrna_seq_analysis.Rmd in RStudio, or from the command line:

Rscript -e 'rmarkdown::render("analysis/scrna_seq_analysis.Rmd")'

Tech stack

Seurat · DoubletFinder · SingleR / celldex · Monocle3 · CellChat · enrichR · SingleCellExperiment · ggplot2

Background

This analysis was originally developed as a project for the Single Cell Bioinformatics course at Saarland University; it has since been reorganized and documented for standalone presentation.

License

Code in this repository is released under the MIT License. The dataset belongs to its original authors (Granja et al. 2019); see data/README.md for provenance.