Two takes on the road into acute myeloid leukemia: a single-cell analysis of preleukemic mouse blood progenitors, and a bulk RNA-seq differential-expression pipeline built as a workflow-engineering exercise. The first asks which cells change before leukemia; the second asks which genes separate AML from healthy blood in bulk cohorts.

Single-cell RNA-seq: preleukemic populations

A reanalysis of 38 mouse bone-marrow HSPC samples across eight preleukemic mutation models (Calr, Dnmt3a, Ezh2, Flt3-ITD, Idh1, Jak2, Npm1c, Utx) from Isobe et al., with a trajectory arm, a TCGA-LAML survival arm, and an R Shiny dashboard. Each stage is a Quarto document, and the rendered HTML files are the analysis record: they explain each method choice in place.

The first version of this analysis used one fixed QC cutoff for all 38 samples, anchor-based integration to an arbitrary reference sample, and cell-level tests for composition and differential expression. This version is a statistical revision. The question is the same; how it is answered changed:

  • Cell calling with emptyDrops (FDR ≤ 0.001) instead of a fixed 200-gene floor, so each barcode is tested against the ambient profile.
  • Per-sample adaptive QC (3 MADs on the log scale) plus doublet removal with scDblFinder. Libraries differ, and a fixed cutoff conflates quality with biology.
  • Harmony integration over 50 PCs with sample as the batch, so no sample is picked as an arbitrary reference and condition is not treated as nuisance.
  • SingleR annotation against the Dahlin 2018 mouse HSPC atlas, with pruned scores and marker verification, so each label carries a confidence.
  • Composition tested with propeller and differential expression with pseudobulk edgeR, both on ~ model + condition. The mouse, not the cell, is the replicate.
  • Trajectory by diffusion pseudotime and CellRank fate probabilities from an Hlf-high HSC root. There are no spliced counts, so no velocity; that limitation is stated rather than worked around.
  • Survival by age-adjusted Cox proportional hazards on the continuous signature score, with Kaplan-Meier curves for display. A median split discards information, and age is the dominant confounder.

What it found. 276,294 barcodes were called as cells and 230,684 remained after QC and doublet removal. Harmony mixed the 38 samples well (per-cluster sample-mixing entropy 0.95 to 0.97, where 1 is even mixing). Only 0.24% of cells were low-confidence under SingleR, and the label distribution matches an LK sort (4.6% HSCs). Then most of the results the original pipeline reported as significant went away. A chi-square on pooled cells gives p < 10⁻¹⁵ for composition, but propeller at the sample level finds no cell type at FDR 0.05. Pseudobulk differential expression finds no shared mutant-versus-WT genes in HSCs and at most 13 in any cell type. A program shared across eight different mutations is not detectable at n = 38, and per-model effects cannot be tested with 2 to 3 mice per arm. In TCGA-LAML (n = 151), neither the paper’s PLPS nor Stem11 signature is associated with overall survival after age adjustment (hazard ratio per SD 0.97 for both; p = 0.75 and 0.81), and the unadjusted result is null too.

Once the animal or the patient is the unit of inference and the covariates are included, most of the original findings do not hold. That is the correct result, not a disappointing one.

Platforms & Tools: R, Python, Quarto, R Shiny, Seurat, DropletUtils, scDblFinder, Harmony, SingleR, propeller / limma, edgeR, Scanpy, CellRank, lifelines, Conda, shinyapps.io

Source data drawn from Isobe et al., Cell Genomics (2023) (GEO GSE227026), and TCGA-LAML clinical data from NCI via cBioPortal. The Quarto stages, helper scripts, and the dashboard live in bioinformatics-public/preleukemia_analysis.

The deployed R Shiny dashboard for the revised analysis. Tabs cover per-sample QC, a UMAP of the integrated HSPC atlas colored by cell type, condition, model, or pseudotime, propeller composition tests, a pseudobulk edgeR volcano per cell type, and TCGA-LAML Kaplan–Meier curves with the age-adjusted Cox results. The app opens in its own window at naraenp2.shinyapps.io/dashboard.

Bulk RNA-seq differential expression (Nextflow)

A compact Nextflow DSL2 pipeline for bulk RNA-seq differential expression, AML vs. healthy, run on real public RNA-seq cohorts. It complements the single-cell work: where that analysis follows preleukemic cells in mice, this one checks that the canonical AML markers separate AML from healthy blood in bulk human cohorts. AML samples come from TCGA-LAML and healthy controls from GTEx whole blood, both pulled from the recount3 project, which re-aligns and re-quantifies TCGA and GTEx through one uniform Monorail / STAR / GENCODE v26 pipeline so the gene-level counts are directly comparable across the two sources.

It’s a workflow-engineering exercise: a small, readable pipeline (channels, processes, publishDir, profile-driven config) on top of a transparent, dependency-light biology layer: library-size CPM normalization, a per-gene Welch t-test on log2-CPM, and a hand-rolled Benjamini-Hochberg FDR. The interactive volcano below labels the canonical AML markers (FLT3, KIT, MEIS1, HOXA9, MPO, CD34, …), which sit cleanly above the significance line.

Four stages:

  1. LOAD_COUNTS: join the TCGA-LAML + GTEx gene sums on Ensembl ID, map to HGNC symbols via GENCODE v26, subsample to balanced groups, and filter low-expression genes on pooled expression so the filter stays independent of the group contrast.
  2. NORMALIZE_COUNTS: library-size CPM, then log2(CPM + 1).
  3. RUN_DE: per-gene Welch t-test with BH-adjusted p-values.
  4. MAKE_VOLCANO: an interactive Plotly volcano.

The real-data inputs (~130 MB from recount3 + the GENCODE annotation) are fetched once with a small fetch_real_data.sh helper, and the whole thing runs in seconds on a laptop. Pinned conda env, project-relative paths, and fast data-free unit tests for the DE math.

Comparator caveat: GTEx has no bone-marrow tissue, so whole peripheral blood is the closest large healthy comparator. Cohort is therefore confounded with disease, tissue, and collection protocol at once. The AML markers still recover cleanly, but progenitor genes can read as “up in AML” simply because mature blood has no progenitor populations, so the direction of a fold change should be read with that in mind. Swapping in a healthy bone-marrow cohort is the natural next step.

Platforms & Tools: Nextflow DSL2, Python (numpy / pandas / scipy / plotly), recount3, GENCODE v26, conda, pytest

The pipeline source and the main.nf workflow live in bioinformatics-public/aml_rnaseq_nf; see docs/REPORT.md for a full run report: dataset provenance, the embedded volcano, and a runtime profile.

Interactive volcano of the differential-expression results, with the canonical AML markers labeled above the significance line.