Limitations of differential expression for cell-type marker discovery in scRNA-seq atlases

Analysis code and derived results for the manuscript of the same name (Doggett & Zhang, National Library of Medicine, NIH).

Compares Wilcoxon rank-sum differential expression against NS-Forest v4.1 across three annotated single-cell atlases: HLCA Core (lung, 61 clusters), Lake kidney (76 clusters), and Allen Brain MTG (75 clusters) -- 212 clusters in total.

---

Quick start

Reproduce all four manuscript figures from stored results:

export PYTHONPATH=/data/doggettkj/workflow_clean/src
cd /data/doggettkj/workflow_clean/figures
rm -rf __pycache__
/data/doggettkj/envs/nsforest/bin/python -u fig1.py
/data/doggettkj/envs/nsforest/bin/python -u fig2.py
/data/doggettkj/envs/nsforest/bin/python -u fig3.py
/data/doggettkj/envs/nsforest/bin/python -u fig4.py

Output lands in figures/output/ as PDF, PNG, and TIFF, plus a figN_stats.csv per figure containing every number quoted in that figure's caption.

The PYTHONPATH line is required. The figure scripts import the shared modules from src/; without it they fail with ModuleNotFoundError.

The rm -rf __pycache__ line is not optional. Stale bytecode has silently run old function signatures after a source edit more than once. Clear it before every run.

fig4.py additionally requires two compute steps to have been run first -- see the Figures section below.

---

Layout

workflow_clean/
|-- env/            environment record + vendored NS-Forest source
|-- src/            shared modules imported by everything else
|-- compute/        scripts that PRODUCED the results (archival)
|-- results/        derived CSVs -- the inputs to the figures
|-- figures/        scripts that DRAW the figures
|   `-- output/     generated figures and caption statistics
`-- archive_notes/  what was left behind in the pre-cleanup dirs

The dependency direction is one-way: compute/ produced results/, and figures/ turns results/ into the paper.

compute/ is a provenance record. Most of its paths deliberately still point at the original working directories, because that is where those jobs actually ran. Re-running it is neither necessary nor expected; the DE and NS-Forest jobs took hours of Slurm time and their outputs are stored in results/. The two exceptions are compute/ec_umap_coords.py and compute/ec_dotplot_matrix.py, which were written after the cleanup, use clean-tree paths, and do need to be run before Figure 4.

results/ is excluded from version control -- it is roughly 1.8 GB and individual DE tables run 150-220 MB, past GitHub's per-file limit. It is deposited to Zenodo at publication.

---

Environment

Python 3.12.8, at /data/doggettkj/envs/nsforest/bin/python. Package versions in env/requirements.txt (79 packages; key ones are scanpy 1.12.1, numpy, pandas, scipy, matplotlib, anndata).

This is the only interpreter on the machine that can import scanpy, anndata, NS-Forest, or the config modules in src/. The system python3 is useful only for syntax checks.

No environment.yml is included -- conda was not available on the machine where this tree was assembled. requirements.txt plus the Python version is the environment record.

Rebuilding elsewhere

pip install -r env/requirements.txt

Every line resolves from PyPI except NS-Forest, which installs from a pinned public commit (see below). No manual steps are needed.

NS-Forest version

All analyses used NS-Forest 4.1, installed from the main branch of https://github.com/NLM-DIR/NSForest at commit a06d731 (2026-04-21, "update tutorials, add histo plotting").

This is a development state on main, after the v4.1 release tag rather than the tag itself. The package self-reports version 4.1 because the version string had not yet been incremented for the following release. This distinction is recorded because installing the v4.1 tag would give slightly different code from what was actually run.

Differences from the v4.1 tag are confined to the command-line demo script, plotting helpers, and preprocessing utilities not used in this work. The code paths that produce every number in the manuscript -- binary scoring, random forest gene ranking, decision tree evaluation, and On-Target Fraction -- are unchanged between the two.

env/nsforest_vendored/ holds a byte-for-byte snapshot of the installed package, retained as a convenience and as insurance. It is identical to the commit above. To install from it instead of from GitHub:

pip install ./env/nsforest_vendored

Two behaviours of this version are worth knowing before reusing the code:

• The add_fraction wrapper does not work. Call nsforest.nsforesting.calculate_fraction.markers_onTarget directly; its signature is (adata, cluster_header, markers_dict).
• DecisionTree calls os.makedirs(output_folder) unconditionally, even when save=False. Pass a non-empty output_folder or it raises.

One further caution for anyone comparing On-Target Fraction across atlases: markers_onTarget normalizes by however many clusters are present in the object handed to it. HLCA has 61, kidney 76, brain 75, so OTF values are comparable within an atlas but not between atlases. Every cross-organ comparison in this manuscript is of within-atlas method contrasts, never of raw OTF values across atlases.

Method citation: Liu A, Peng B, Pankajam A, Duong TE, Pryhuber G, Scheuermann RH, Zhang Y. Discovery of optimal cell type classification marker genes from single cell RNA sequencing data. BMC Methods 1, 15 (2024). doi:10.1186/s44330-024-00015-2

---

Raw data (not included)

The three atlases are published datasets and are NOT redistributed here. The config files in src/ point at local copies:

• HLCA Core: /data/doggettkj/HLCA_Core/HLCA_core.h5ad 584,944 cells x 27,402 genes, 61 clusters (ann_finest_level). var_names are Ensembl IDs. Pre-normalized; .X is log-normalized. Sikkema et al., Nat Med 2023. Accession: TODO
• Kidney: /data/doggettkj/Kidney/kidney_healthy.h5ad 107,701 cells, 76 clusters. The annotation column is subclass.l3 on disk but subclass_l3 in the derived CSVs; resolve_ann_col() in the config modules handles the difference. Pre-normalized. Lake et al., Nature 2023. Accession: TODO
• Brain MTG: /data/doggettkj/Brain_MTG/brain_mtg_normalized.h5ad 15,603 cells x 50,281 genes, 75 clusters (cluster_id). SMART-seq. Derived locally: the Allen release shipped counts-only in .raw.X with .X empty, so it was normalized here to target_sum=1e4 + log1p to match the convention of the other two atlases, and var_names were set from the gene symbols in adata.var. Hodge et al., Nature 2019 (AIT15.3). Accession: TODO

Most figures never open an h5ad. The exceptions are Figure 2 panel A (cluster size counts) and the two Figure 4 compute steps. Figures 1 and 3 run entirely from the CSVs in results/.

---

Figures

Run every script from figures/, with src/ on the path and the bytecode cache cleared:

export PYTHONPATH=/data/doggettkj/workflow_clean/src
cd /data/doggettkj/workflow_clean/figures
rm -rf __pycache__
/data/doggettkj/envs/nsforest/bin/python -u figN.py

Main figures:

• fig1.py -- Figure 1. The size confound. 1A cluster-size distributions (all three atlases), 1B lung-only DE score and NS-Forest F-beta against log size, 1C counts of practically unrankable genes against size. Reads NS-Forest and DE CSVs plus fig1c_unrankable.csv. Runtime: seconds.
• fig2.py -- Figure 2. Causality and the repair. 2A partial correlations, 2B downsampling dose-response, 2C simulation, 2D the rank-based effect size A. Panel code lives in fig2_panels.py (draw_2a, draw_2b, draw_2c), fig2b_downsampling.py, and fig2d_table.py (draw_2d_table). Reads CSVs plus all three h5ad files. Runtime: several minutes.
• fig3.py -- Figure 3. Specificity. 3A paired violins, 3B grouped four-metric bars. Panel code in fig3_panels.py, which exposes draw_3a(ax) and draw_3b(axes). Reads nsf_native_DE_vs_NSF_Wilcoxon.csv. Runtime: seconds.
• fig4.py -- Figure 4. The endothelial case study. 4A EC arterial UMAP triptych, 4B both genes across all eight endothelial clusters, 4C DE rank trajectory across four background compositions. Reads results/ec_arterial/. REQUIRES the two compute steps below, or it exits with a missing-file error. Runtime: seconds.

Supplementary figures:

• figS1.py -- kidney and brain replication of Figure 1B.
• figS2.py -- kidney and brain replication of Figures 2A and 2B.

S3 and S4 are not yet scripted in this tree.

Checks:

• verify_deauc_ceiling.py -- ceiling-compression counts behind Figure 2D. Reads de_auc_matchedk.csv. Runtime: seconds.

Figure 2 is the slow one

Panel A opens each h5ad and bootstraps twelve partial correlations. Expect a long silent pause during the HLCA load; python -u is what makes progress visible.

Figure 4 compute steps

Both read the HLCA h5ad and write into results/ec_arterial/. Run them once, before fig4.py:

export PYTHONPATH=/data/doggettkj/workflow_clean/src
cd /data/doggettkj/workflow_clean
rm -rf compute/__pycache__
/data/doggettkj/envs/nsforest/bin/python -u compute/ec_umap_coords.py
/data/doggettkj/envs/nsforest/bin/python -u compute/ec_dotplot_matrix.py

They produce umap_positive.csv, umap_background.csv, and dotplot_matrix.csv. Both are slow -- these are the last full h5ad reads in the pipeline.

Style

All visual conventions -- widths, colours, fonts, axis helpers -- live in src/figstyle_frontiers.py, which every figure script calls first via fs.apply(). It is the single source of truth; do not set styling in individual scripts. After editing it, run fig1.py as a canary.

Figures carry no axes titles and no suptitles by design: captions carry every claim, and captions fall outside the journal's word limit.

Note on Table 1

An earlier draft carried a Word table of the rank-based effect size results. It is now drawn as Figure 2D by fig2d_table.py, which recomputes its values from de_auc_matchedk.csv rather than transcribing them, so the panel cannot drift from source. There is no separate Table 1 in the current manuscript.

---

Simulation

compute/simulation/sim_tiers.py regenerates results/simulation/sim_tiers_de_vs_nsf.csv, the input to Figure 2C.

cd /data/doggettkj/workflow_clean/compute/simulation
/data/doggettkj/envs/nsforest/bin/python -u sim_tiers.py

42 synthetic clusters (3 marker-quality tiers x 14 log-spaced sizes from 30 to 30,000 cells) over 15 seeds. Only the off-target expression mean differs between tiers; in-target mean is 6.0 throughout. This is a vacuum experiment, not a biologically calibrated one -- that is deliberate, and avoids circularity.

Expected: realized OTF 1.000 / 0.528 / 0.140 (high/med/low), overall r(DE, log size) = 0.899.

Reproducibility note. Re-running reproduces the stored CSV to within ~1e-4 absolute in DE_score and ~1e-7 in NSF_OTF, but not bitwise. The expression matrix is float32, and summing hundreds of thousands of values gives slightly different totals depending on how the numerical library schedules the work. Cluster labels, sizes, and all reported aggregate statistics are identical. The stored CSV is the one the published figure was drawn from.

---

---

Background-composition ladder (Figure 4C)

compute/generalization/ holds the code behind Figure 4C's rank trajectory. ladder.py recomputes DE for a fixed target cluster against progressively wider background groups, walking up the HLCA class hierarchy: endothelium alone, then adding stroma, immune, and the remaining atlas. run_generalization.py drives it across clusters; submit_generalization.sbatch is the Slurm wrapper.

These are provenance, like the rest of compute/. Budget roughly 21 minutes per rung for a 7,000-cell target -- each rung is a full DE recomputation against a different background, not a re-scoring of stored results.

Their outputs for the EC arterial cluster are already in results/ec_arterial/: gen_EC_arterial.json gives the rung structure and NS-Forest's pick at each rung, and de_fullrank/ holds the full DE ranking within each of the four backgrounds. fig4.py reads both.

The same experiment was run across twenty HLCA clusters. NS-Forest's marker was more on-target than DE's in twenty of twenty. That result is not reported in the manuscript and its outputs are not included here; the EC arterial case shown in Figure 4 is the one the paper uses.

Relocating this tree

All paths are hardcoded absolute -- deliberately, so any line can be read without resolving a variable. To move the tree, edit:

• src/config_hlca.py, src/config_kidney.py, src/config_brain.py (lines 6, 7, 9: H5AD_PATH, RESULTS_DIR, FIGURES_DIR)
• the path constants at the top of each script in figures/, and in compute/ec_umap_coords.py and compute/ec_dotplot_matrix.py

H5AD_PATH is the exception worth leaving alone: it points at the published third-party atlases, not at anything this tree owns.

Find them all with:

grep -rn "/data/doggettkj" src/ figures/ compute/

---