Pr90 testland#133
Merged
Merged
Conversation
Implements STARsolo Phases 14.1–14.4 (the 10x Chromium Gene-count MVP) plus
the CellRanger 4.x/5.x-matching flag set, all ported faithfully from STAR's
source and verified byte-identical against real STARsolo.
New `src/solo/`:
- mod.rs — SoloType/params plumbing, barcode-read input (SoloReadReader),
SoloContext, per-read processing, CellRanger4 adapter clip
- whitelist.rs — 2-bit barcode packing + sorted-whitelist load + read-stage
CB matching (exact/1MM/1MM_multi) + UMI checks (STAR
SoloReadBarcode_getCBandUMI.cpp)
- gene.rs — per-read gene assignment (--soloStrand), reuses
quant::overlapping_genes
- count.rs — UMI dedup (Exact/NoDedup/1MM_All/1MM_Directional/1MM_CR),
MultiGeneUMI(_CR) filtering, 1MM_multi_Nbase_pseudocounts CB
posterior, raw matrix.mtx/barcodes.tsv/features.tsv writer
Driver: new align_reads_solo loop in lib.rs (reads cDNA + barcode in lockstep,
aligns cDNA, quantifies per cell); solo params + validation in params/mod.rs.
CellRanger-matching flags (--clipAdapterType CellRanger4, --outFilterScoreMin,
--soloCBmatchWLtype 1MM_multi_Nbase_pseudocounts, --soloUMIfiltering
MultiGeneUMI_CR, --soloUMIdedup 1MM_CR) produce a matrix byte-identical to real
STARsolo (3/3 deterministic) — verified via a Linux-container differential
harness (test/solo_cellranger_diff.py + Dockerfile.solodiff + solo_diff_docker.sh),
since STAR 2.7.11b reads 0 reads on Apple-Silicon macOS.
Also adds a CellRanger-vs-STARsolo-vs-rustar runtime/stats benchmark scaffold
(test/solo_bench.py + Dockerfile.bench).
Tests: 479 lib + 11 integration (incl. test_starsolo_cellranger_style_matrix),
0 clippy warnings. Docs in docs-old/phase14_starsolo.md + ROADMAP.md.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Performance + scalability fixes for STARsolo matrix output, validated by a CellRanger-vs-STARsolo-vs-rustar benchmark on a real 5' mouse 10x dataset. - Buffered I/O: write_barcodes / write_matrix_mtx / write_features wrote to a raw std::fs::File (one syscall per line). For the full 3.69M-barcode barcodes.tsv that was ~3.7M syscalls, dominating runtime (esp. over virtiofs). Wrap in BufWriter + add CbWhitelist::unpack_barcode_into (no per-line String alloc). On the mouse bench, raw-matrix write dropped 1306s -> 3s, output byte-identical. - build_matrix Step 1 (per-cell): sort the flat record list by cell barcode and process one cell's contiguous slice at a time, so peak memory is a single cell's umi->gene maps rather than a global cell->umi->gene nested map over all records (the previous version's tens-of-GB allocator blowup). Mirrors STAR SoloFeature_collapseUMIall.cpp. Output identical. - Benchmark harness test/solo_bench.py (+ Dockerfile.bench): runs all three tools under /usr/bin/time -v and reports runtime/peak-RSS/matrix stats; results + method in docs-old/phase14_benchmark.md. 479 lib + 11 integration tests, 0 clippy warnings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
build_matrix previously materialized a global cell -> (gene -> count) map for the whole raw matrix before writing. Replace it with stream_matrix: process one cell at a time (records already sorted by cb) and write each cell's entries straight to a temporary MatrixMarket body, counting nnz on the fly, then emit matrix.mtx as the header (rows cols nnz) followed by the body (the BySJout temp-file pattern). The global output map is never built, so matrix-output memory is bounded by a single cell regardless of how many cells the raw whitelist matrix spans. Verified byte-identical to the prior output on the real 10M-read mouse benchmark (matrix.mtx / barcodes.tsv / features.tsv all identical, 2,650,591 entries). Note: this bounds matrix-output memory for large raw matrices but does not cut the current ~37 GB peak, which is dominated by the loaded index (~27 GB SA) plus alignment buffers — mmap'ing the index is the lever for that. 479 lib + 11 integration tests, 0 clippy warnings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Loading the suffix array and SAindex with std::fs::read materialized them as
anonymous Vecs (~21 GB + ~1.8 GB for mouse). Anonymous memory is not
reclaimable, so under pressure the OS swaps it — which froze the host (Jetsam)
on a real solo run.
Back PackedArray with PackedBytes { Owned(Vec) | Mapped(Arc<Mmap>) }: build still
uses Owned (write() requires it); load now memory-maps the SA and the SAindex
packed-data region (read normally past the small header) with MADV_RANDOM. The
read path is unchanged (operates on &[u8] either way). Arc<Mmap> keeps the
two-pass GenomeIndex clone cheap.
mmap pages are file-backed: demand-loaded and reclaimable (dropped, never
swapped) under pressure, so a multi-GB SA no longer pins anonymous RAM and the
run degrades to paging instead of OOM-killing.
Validated on the real 10M-read mouse benchmark (GRCm39, 21 GB SA): matrix.mtx /
barcodes.tsv / features.tsv byte-identical to the in-RAM result, and the run
COMPLETES under a 24 GB container cap (peak RSS 23.4 GB, OOMKilled=false) where
the previous ~37 GB anonymous footprint would have been OOM-killed. (Genome's
forward+RC buffer, ~5.6 GB, is still an owned Vec — Phase 2.)
479 lib + 11 integration tests, 0 clippy warnings.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…se 2)
load_genome built a 2*n_genome Vec (~5.6 GB for mouse) holding the forward
strand plus a precomputed reverse complement. Replace Genome.sequence's type
with GenomeSeq { Owned(Vec) | Mapped { fwd: Arc<Mmap>, n_genome } }:
- Build stays Owned (the full forward+RC buffer is needed for SA construction
and the on-disk write, and only Owned supports slicing/mutation).
- Load memory-maps the on-disk Genome file, which holds ONLY the forward strand
(n_genome bytes). The reverse-complement half is computed on access in
GenomeSeq::base(i) = complement(fwd[2n-1-i]) for i >= n. So the ~2.8 GB RC
buffer is never materialized, and the forward bytes are reclaimable
file-backed pages instead of an anonymous Vec.
Hot alignment access was already single-byte (seed extension, SAindex lookup),
so it routes through base()/get(); build-time slicing routes through as_slice()
(always Owned). PartialEq/Debug are hand-implemented (memmap2::Mmap has neither).
Combined with the SA + SAindex mmap (commit 4480806), all three large index
components are now reclaimable file-backed memory; the anonymous floor of a solo
run is just the alignment working buffers.
Validated on the real 10M-read mouse benchmark (GRCm39, 2.79 Gb genome):
matrix.mtx / barcodes.tsv / features.tsv byte-identical to the precomputed-RC
result — i.e. RC-on-access is correct for every reverse-strand alignment —
with OOMKilled=false.
479 lib + 11 integration tests, 0 clippy warnings.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
STARsolo's `--soloFeatures GeneFull` (CellRanger's include-introns default):
a read counts toward a gene if it overlaps the gene's full body (exons +
introns), so purely intronic reads are counted — unlike exonic `Gene`.
- GeneAnnotation gains chr_gene_body (one [min exon start, max exon end) span
per gene) + overlapping_genes_full(); overlapping_genes() refactored to a
shared interval-overlap helper.
- gene.rs: SoloFeature{Gene,GeneFull}; assign_gene_se() takes the feature and
picks exon vs gene-body overlap.
- SoloContext now holds a list of features, each with its own SoloRecorder;
process_read does the shared CB-match + UMI check once and assigns a gene per
feature. write_gene_matrix writes one Solo.out/<feature>/raw/ per feature, so
`--soloFeatures Gene GeneFull` produces both matrices in one pass.
- params: validate soloFeatures (only Gene/GeneFull supported).
480 lib + 11 integration tests (incl. genefull_counts_intronic_read), 0 clippy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Faithful port of STAR's SoloFeature_emptyDrops_CR.cpp (CellRanger's EmptyDrops
variant) as a standalone `emptydrops` binary that post-processes any raw count
matrix (MatrixMarket genes×cells + barcodes/features tsv, plain or .gz):
1. guaranteed cells from the CellRanger-2.2 knee (nExpectedCells/maxPct/ratio);
2. ambient RNA profile from barcode ranks [indMin,indMax) with a Good-Turing
P0 unseen-mass correction (approximating STAR's SGT smoothing);
3. candidate barcodes (rank >= nSimple, total >= max(umiMin, frac*median));
4. per-candidate multinomial log-likelihood under the ambient profile;
5. Monte-Carlo p-values — simN ambient-drawn barcodes, running log-prob per
count reused across candidates; p = (1+#sim<obs)/(1+simN);
6. Benjamini-Hochberg FDR; cells with padj <= FDR added to the guaranteed set.
Defaults mirror `--soloCellFilter EmptyDrops_CR 3000 0.99 10 45000 90000 500
0.01 20000` with seed 19760110 and simN 10000. Writes barcodes.tsv + cells.txt
+ emptydrops.json. Lets us call cells on rustar/STARsolo/CellRanger raw matrices
with one algorithm for an apples-to-apples filtered comparison.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
- solo_genefull_compare.py: GeneFull raw-count parity (total UMI / genes /
per-cell correlation) across rustar/STARsolo/CellRanger + EmptyDrops cell-set
overlap (Jaccard) vs CellRanger's filtered barcodes.
- solo_compare_h5ad.py: optional --{rustar,starsolo,cellranger}-cells to filter
each raw matrix by an EmptyDrops cells.txt instead of the CR2.2 knee, so the
written h5ad reflect EmptyDrops cell calls for the CellRanger comparison.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
solo_genefull_h5_compare.py loads matrices one at a time (memory-careful) and reports the intron effect (Gene vs GeneFull total UMI), raw-count parity vs CellRanger, EmptyDrops cell-set agreement (same-algorithm + vs CR native filtered), per-cell UMI correlation on shared cells, and writes the EmptyDrops-filtered h5ad for both tools. Result on mouse 5k-PBMC (10M subsample): GeneFull +18% UMI over Gene, within 2.8% of CellRanger; rustar-ED captures 3857/3858 CR native cells; per-cell UMI r=0.9995; same-algorithm Jaccard 0.9749. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
CellRanger appends '-1' to barcodes; STARsolo/rustar do not. Normalize on load so the rustar/STARsolo vs CellRanger per-cell correlations match cells instead of reporting shared=0. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…fers
Two behavior-preserving speedups (count matrices verified byte-identical on the
mouse 5k-PBMC 10M-read set):
1. `--outSAMtype None` (previously errored "not yet implemented") now routes to
NullWriter, and the SE/PE/solo loops skip building SAM RecordBufs entirely
via params.emits_alignments(). Quant/solo runs that only need the count
matrix no longer pay SAM text formatting + the multi-GB Aligned.out.sam
write. ~12-15% faster natively; far more in a container where the SAM write
crosses virtiofs.
2. Gene assignment (assign_gene_se) ran once per feature per read — twice/read
for `Gene GeneFull` — allocating fresh Vecs each call. It now reuses two
thread-local scratch buffers, and GeneAnnotation gains
overlapping_genes{,_full}_into(transcript, &mut buf) so the overlap query
writes into a caller-provided buffer instead of allocating.
480 lib + 11 integration tests pass, 0 clippy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
`memmap2::Advice` / `Mmap::advise` (madvise) are Unix-only, so the SA + SAindex mmap load failed to build on windows-x86_64. Wrap the best-effort `advise(Advice::Random)` in a cfg(unix) `advise_random` helper with a no-op cfg(not(unix)) stub. No behavior change on Unix. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Each `Solo.out/<feature>/Summary.csv` now reports the sequencing/mapping funnel and per-cell statistics, matching STARsolo's format: - Number of Reads, Reads With Valid Barcodes, Sequencing Saturation - Reads Mapped to Genome (Unique+Multiple / Unique), Reads Mapped to <feature> - genome -> Exonic / Intronic / Intergenic confident-mapping split (computed from the Gene vs GeneFull read counts: exonic=Gene, intronic=GeneFull-Gene, intergenic=genome_unique-GeneFull) - Estimated Number of Cells (CR2.2 knee), reads/UMIs in cells, fraction in cells, mean/median reads/UMI/genes per cell, total genes detected Mechanism: SoloContext gains per-feature `feature_reads` atomics (incremented in process_read on a unique gene assignment); stream_matrix returns per-cell (reads, UMIs, genes) + genes-detected; write_gene_matrix takes AlignmentStats and writes the summary. Matrix output is unchanged (byte-identical verified). test/solo_summary_compare.py cross-tabulates rustar / STARsolo Summary.csv and CellRanger metrics_summary.csv. On mouse 5k-PBMC (10M): saturation 14.5/14.5/14.6%, genome 92.0/92.0/91.9%, median genes/cell 601/601/599, median UMI/cell 915/916/908; rustar exonic/intronic/intergenic 49.5/8.8/16.9% vs CellRanger 53.5/11.4/15.9%. 482 lib tests (+median/knee helpers), 0 clippy. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…sense) The Summary.csv mapping funnel previously derived exonic/intronic from the unique-gene-assignment counts (Gene, GeneFull), which dropped multi-gene (ambiguous) reads into the intergenic residual and folded antisense reads in too. Replace it with CellRanger's method: classify each uniquely-mapped read by position (exonic if it overlaps any exon, else intronic if in a gene body, else intergenic) independent of strand, and report antisense (maps to a gene body on the opposite strand) as a separate orientation metric. `classify_read` does this in a single pass, sharing the two overlap queries with the per-feature gene assignment, so process_read does no more gene-model work than before (the two assign_gene_se calls are replaced by one classify_read). Region/antisense are counted over uniquely-mapped reads regardless of barcode (CellRanger's "confidently mapped" = MAPQ 255); matrix output is unchanged (byte-identical verified). Result vs CellRanger (mouse 5k-PBMC, 10M): intronic 10.8% vs 11.4%, intergenic 16.1% vs 15.9%, antisense 5.9% vs 6.0% — all within ~0.5pp (were 8.8/16.9/—). The residual exonic gap (48.4 vs 53.5%) is the genome-mapping denominator (STAR unique 75.2% vs CellRanger confident 80.8%) + CellRanger's >=50%-exon / splice-compatibility rule, not the binning method. Runtime unchanged (~78s/10M), 483 lib tests, 0 clippy. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Phase 14.6 — write Solo.out/<feature>/filtered/ (the called cells), not just
raw/. `--soloCellFilter`:
- CellRanger2.2 <nExp> <maxPct> <ratio> (default) — CR2.2 knee
- TopCells <N> — top-N barcodes by UMI
- None — skip filtered/
- EmptyDrops_CR — knee-guaranteed cells (Monte-Carlo rescue stays in the
standalone `emptydrops` binary; logged)
Mechanism: stream_matrix is split into build_matrix_body (per-cell dedup → shared
plain temp body + per-cell CellStat{cb,reads,umis,genes}) and finalize_matrix
(header + verbatim body for raw, or cb-remapped/renumbered body for filtered).
Raw and filtered share the one streaming pass; the filtered matrix re-reads the
temp body filtering by called cb. Raw output byte-identical (verified).
`--soloOutGzip yes` gzips matrix.mtx/barcodes.tsv/features.tsv (raw + filtered)
and appends .gz (CellRanger-style); default `no` keeps STARsolo's plain files so
the byte-for-byte STARsolo comparison still holds. write_file() finishes the gzip
stream explicitly.
Verified on mouse 5k-PBMC (10M): GeneFull filtered 3821 cells vs STARsolo 3820
(Jaccard 0.9997, 0 STARsolo-only); filtered UMI sum == Summary "UMIs in Cells";
all 6 .gz files pass gunzip -t; runtime unchanged (~85s/10M). 484 lib + 12
integration tests, 0 clippy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
`--soloFeatures SJ` writes Solo.out/SJ/raw/{features.tsv,barcodes.tsv,matrix.mtx}
where rows are the SJ.out.tab junctions (STARsolo symlinks features → SJ.out.tab;
rustar writes the identical 9-column lines in the same chr/start/end-sorted
order). Per uniquely-mapped read the crossed junctions (absolute intron coords
from extract_junction_keys) are recorded with the resolved CB+UMI, then collapsed
per (cell, junction) by --soloUMIdedup; junctions filtered out of SJ.out.tab are
dropped.
Plumbing: SoloContext gains sj_enabled + sj_records; process_read takes the read's
junctions and emits SjCountRecord; the solo loop extracts junctions for unique
spliced reads. SpliceJunctionStats gains sj_feature_order() (row order) and a
shared write_sj_lines() (used by both SJ.out.tab and the SJ features.tsv). The
post-run solo output moved into run_single_pass (write_solo_output) where sj_stats
is live. Respects --soloOutGzip.
Verified on mouse 5k-PBMC (10M): SJ features.tsv byte-identical to SJ.out.tab
(100,970 junctions), matrix 100,970 × 3,686,400, 1.41M entries / 1.87M molecules,
all rows/cols in range; runtime unchanged (~82s). Integration test
test_starsolo_sj_feature plants a spliced read across the GT-AG intron and asserts
the one-junction matrix. 501 tests, 0 clippy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Reads mapping to multiple genes (gene-ambiguous), previously dropped, are now
distributed across their gene set and written as real-valued
UniqueAndMult-<method>.mtx alongside the unique matrix.mtx in raw/.
- classify_read gains want_multi → returns the sense gene set for ambiguous reads
(Gene + GeneFull); process_read records a MultiGeneRecord (resolved CB) per
ambiguous read into recorder.multi_gene.
- build_multi_matrices re-reads the raw matrix body (per-cell unique counts u_g,
cb-ascending) and merges each cell with its ambiguous molecules (deduped by UMI,
gene set = union), then per method:
Uniform — 1/N to each gene in the set
PropUnique — proportional to unique counts (uniform fallback if none)
Rescue — proportional to unique + a uniform multi prior
EM — iterate theta_g = u_g + (multi ∝ theta) to convergence
Each ambiguous molecule contributes total mass 1.0, so every method conserves
the grand total (only the spread differs). Cells with only multi reads are
skipped (no unique row).
Verified on mouse 5k-PBMC (10M, Gene): all four matrices sum to 4,527,229
(= 4,219,608 unique + 307,621 multi molecules); Uniform/Rescue spread wider (2.80M
nnz) than PropUnique/EM (2.77M); values real (e.g. 0.5, 1.33333). Runtime
unchanged (~79s, multi is a second pass over the raw body). Unit test
distribute_multi_methods + integration test test_starsolo_multimappers (overlapping
genes G1/G3). 503 tests, 0 clippy.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
solo_sj_multi_compare.py compares the SJ matrix (junctions coord-matched, since each tool has its own SJ.out.tab; STARsolo symlinks features→SJ.out.tab) and the UniqueAndMult-*.mtx matrices (genes/barcodes align by GTF/whitelist order). Result on mouse 5k-PBMC (10M, CellRanger params): multiMappers correlate r=0.98-0.998 vs STARsolo; SJ junction sets overlap 99.9% with per-junction r=0.93. Totals run ~4% (Gene) / ~10% (SJ) higher — the known "rustar maps slightly more reads" gap, amplified for splice-spanning reads — not a distribution-logic difference. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Plate-well full-length protocols (Smart-seq2): no cell barcodes or UMIs in the reads. A --readFilesManifest TSV (read1<TAB>read2<TAB>cellID) lists one library per cell; reads are aligned, gene-assigned (Gene/exon), and counted per cell with no UMI dedup (each uniquely-gene-assigned read = a count). - src/solo/smartseq.rs: manifest parser (single-end MVP; read2 must be '-'), SmartSeqCounts (per-cell gene→count), and the genes×cells matrix writer (barcodes.tsv = manifest cell IDs; respects --soloOutGzip). - run_smartseq in lib.rs builds the gene model from --sjdbGTFfile, iterates each cell's reads (rayon per batch), and writes Solo.out/Gene/raw/. Dispatched ahead of the droplet solo_ctx (which needs a whitelist). - params: --readFilesManifest; SmartSeq exempt from the --readFilesIn / whitelist-correction / two-read-file rules; requires the manifest. Verified on mouse cDNA reads (2 synthetic cells → per-cell read counts) and an integration test test_starsolo_smartseq (Exon1 reads → G1 counts 5/3). 504 tests, 0 clippy. PE SmartSeq is a follow-up. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Chemistries whose cell barcode is split into several fixed-position segments (sci-RNA-seq, SPLiT-seq, …), each with its own whitelist. - SoloBarcodeLayout is now an enum: Simple (CB_UMI_Simple, unchanged) and Complex, which parses --soloCBposition/--soloUMIposition (startAnchor_startDist_endAnchor_ endDist; read-start anchoring) and assembles the CB by concatenating the segment slices from the barcode read. - CbWhitelist::load_complex builds the combined whitelist as the cartesian product of the per-segment whitelists (concatenated, packed). Matching the assembled CB against this is equivalent to STARsolo's per-segment matching for both Exact and 1MM (a 1MM in the concatenation is a 1MM in exactly one segment), so the rest of the pipeline (correction, UMI dedup, matrix) is reused unchanged. - params: --soloCBposition / --soloUMIposition; CB_UMI_Complex requires one position + one whitelist per segment. Adapter-anchored / variable positions are a follow-up (read-start only for now). Unit tests (parse_position, complex_layout_assembles_segments) + integration test test_starsolo_cb_umi_complex (2 segments × 2 whitelists → 4-cell product, CB AAGG matched, 1 molecule). 507 tests, 0 clippy. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
--soloCellFilter EmptyDrops_CR now writes filtered/ with the full rescue (knee guaranteed cells + cells whose profile is significantly non-ambient), instead of just the knee. emptydrops_called re-reads the raw matrix body for the ambient (rank [indMin,indMax)) and candidate (rank >= nSimple, total >= minUMI) cell profiles, then runs the same multinomial Monte-Carlo + Benjamini-Hochberg as the standalone `emptydrops` binary (seed 19760110, Good-Turing P0 ambient). Params: EmptyDrops_CR nExpected maxPct maxMinRatio indMin indMax umiMin umiMinFracMedian candMaxN FDR [simN]. Verified on mouse 5k-PBMC (10M, GeneFull, umiMin=100): 3821 knee + 629 rescued = 4450 filtered cells — identical to the standalone binary's result on the same matrix. With STAR's default umiMin=500 there are no sub-knee candidates so it reduces to the knee. 507 tests, 0 clippy. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
--readFilesManifest read2 may now be a mate-2 file (not just '-'). For PE cells, run_smartseq reads mate pairs in lockstep (PairedFastqReader), aligns each with align_paired_read, and counts the fragment once toward the gene from the union of both mates' overlaps (both-mapped → both transcripts; half-mapped → the mapped mate). SE manifests are unchanged (read counts). Verified end-to-end on real Smart-seq2 mouse data (GEO GSE228456 / SRP429940, 3 sorted single monocytes, 1M reads each): PE fragment counts run slightly below SE read counts (stricter proper-pairing, one count/fragment) and detect more genes (read2 covers extra regions) — as expected. Integration test test_starsolo_smartseq_paired (mate1 Exon1 + mate2 rc(Exon2) → proper FR pair on G1, 4 fragments). 508 tests, 0 clippy. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
--soloFeatures Velocyto writes Solo.out/Velocyto/raw/{spliced,unspliced,ambiguous}
.mtx (genes × cells), the scVelo/dynamo-ingestible input for RNA velocity. Unlike
classic Velocyto's spliced/unspliced heuristic, this uses the Sullivan et al. (NAR
2025) mature/nascent/ambiguous classification, which rustar computes exactly from
the alignment (no pseudoalignment D-list tricks needed — alignment already places
reads by coordinate):
- spliced = the read splices (a junction in the CIGAR) → processed mRNA;
- unspliced = no junction but an aligned block leaves the exons into an intron;
- ambiguous = no junction, block wholly within an exon (origin indistinguishable
— kept separate rather than folded into spliced, the paper's key point).
Gene is assigned by gene-body overlap (so intronic/nascent reads count). Per
(cell,gene) each UMI is resolved to one category (priority unspliced > spliced >
ambiguous) then UMI-deduped per category.
Mechanism: GeneAnnotation gains per-gene merged exon intervals + block_is_exonic;
gene.rs velocyto_category(); SoloContext.velocyto_records; process_read records on
gene_full assignment; build_velocyto_matrices writes the 3 matrices. Respects
--soloOutGzip.
Verified on mouse 5k-PBMC (10M): spliced 1.77M / unspliced 1.29M / ambiguous 1.84M
molecules, all in range (high ambiguous = short reads within exons, as the paper
expects). Integration test test_starsolo_velocyto (one read per category → one
molecule each). 509 tests, 0 clippy. Closes the last major STARsolo feature gap.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Swap the deflate backend from flate2/miniz_oxide to libdeflate: - solo matrix output (--soloOutGzip): write_file now collects the body and gzip-compresses it in one libdeflater call (level 6) instead of streaming flate2. libdeflate has no streaming API, but matrix/barcode bodies are tens of MB and fit in memory. - BGZF/BAM: add noodles-bgzf with its `libdeflate` feature; Cargo feature unification routes all BGZF deflate (BAM output, stdout) through libdeflate — no code change in io/bam.rs. Measured on real rustar data (mouse 5k-PBMC GeneFull matrix + full whitelist): - matrix.mtx (49MB) L6: 2.9x faster (1031→361ms), -3% size - barcodes.tsv (63MB) L6: 3.2x faster (2060→647ms), +6% size - BGZF blocks L1 (BAM default): 1.4x faster AND -21% size (libdeflate's L1 beats miniz_oxide's both ways); L6: 3.4x faster, same size Output is standard gzip/BGZF (zcat/samtools read it unchanged); 509 tests pass incl. BAM round-trip through libdeflate BGZF. (Evaluated fastrgz first — it's a FASTQ-only compressor with a non-transparent FASTR transform, slower decode, and no library API, so it doesn't fit rustar's matrix/BGZF gzip; libdeflate does.) Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ops MC - align_reads_solo: background producer thread overlaps the single-threaded gzip inflate with the parallel alignment (was a per-batch serial section stalling every rayon worker). Warm-cache A/B on mouse 10M / 4 features / 10 threads: alignment phase ~68s -> ~54s (~20% faster). - count.rs build_matrix_body & build_velocyto_matrices: par_sort the record vec and run per-cell UMI dedup/format in parallel, emitting pre-formatted per-cell byte buffers sequentially in CB order (output byte-identical). par_sort also in build_sj_matrix / build_multi_matrices. - emptydrops_called: per-simulation splitmix-seeded RNG + rayon parallel Monte-Carlo. Deterministic and thread-count-independent (verified run A == run B). Per-sim seeding draws differ from the old single stream, so the rescued set shifts by Monte-Carlo noise (Gene 408->407, GeneFull 629->628; 1 boundary cell each at FDR=0.01); knee cells and raw matrices unchanged. Raw/filtered matrices byte-identical for #1-2; only EmptyDrops rescue differs. 509 tests passing, 0 clippy warnings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Velocyto ambiguous handling (rustar extension beyond STARsolo): - --soloVelocytoAmbiguous yes (default) keeps the STARsolo 3-matrix output (spliced/unspliced/ambiguous). `no` folds exon-only molecules into spliced and writes only spliced/unspliced, no ambiguous.mtx (an exon-only read is most likely mature mRNA; cf. He, Soneson & Patro 2023). unspliced is untouched; verified spliced gains exactly the folded ambiguous mass. Fix (introduced by the pipelined-decode change in the prior commit): the background FASTQ-decode producer thread fed a bounded depth-2 channel, but on an early consumer exit (--readMapNumber, or a `?` error) the consumer stopped draining while the producer blocked on a full channel forever, deadlocking the scope join. Full runs never hit it (consumer reads to EOF). The consumer now runs in an IIFE and drops the receiver on every exit path, waking a blocked producer (its send returns Err) so the scope joins cleanly. Tests: test_starsolo_velocyto_fold_ambiguous (fold output) and test_starsolo_readmapnumber_terminates (30k-read early-exit guard, sized to fill the channel). 511 tests passing, 0 clippy warnings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…its own file Per PR #90 review (Psy-Fer): for the faithful-port / drop-in release, the STARsolo Summary.csv must not be altered. The CellRanger-style mapping funnel (Reads Mapped Confidently to Exonic/Intronic/Intergenic Regions + Reads Mapped Antisense to Gene) is moved out of Summary.csv into a separate additional file, Solo.out/<Feature>/CellRanger.summary.csv, written only when Gene + GeneFull run. Adding a new file doesn't change the existing STARsolo outputs, so drop-in compatibility holds; output-changing features come in later PRs. - write_summary no longer takes/writes the region funnel - new write_cellranger_summary writes the funnel to CellRanger.summary.csv - test_starsolo_summary_split: Summary.csv has no funnel rows; the funnel rows live in CellRanger.summary.csv 512 tests passing, 0 clippy warnings. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Profiling (warm align, mouse 10M, AVX-512 EC2) showed cluster_seeds as the #1 self-time hotspot (~19%), with hashbrown reserve_rehash + SipHash hash_one in its call tree. The win_bin (strand,bin)->window and diag_ranges (diag,mate) maps have integer keys but used the default SipHash hasher. Switch both to rustc_hash FxHashMap and pre-size to kill rehashing. Output is bitwise-identical (512 tests pass); pure speedup. (x86-64-v4/AVX-512 build was ruled out — no effect.) Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Profiling (warm align, mouse 10M, AVX-512 EC2) showed clear_page_erms (kernel page zeroing) at ~12% of align time plus ~10% in mimalloc. mimalloc's default ~10ms purge delay returns freed pages to the OS; under the per-read alloc/free churn they are immediately re-faulted and re-zeroed. Set mi_option_purge_delay = -1 (never purge, retain freed pages) via libmimalloc-sys FFI at startup — mi_option_set overrides the cached default even though the allocator initialized before main(). Local mouse-10M warm align: sys time ~22s -> ~11s, wall ~8% faster, Max RSS unchanged (the 22 GB mmap'd index dominates; retained heap is negligible). Output unchanged; 512 tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Targets the residual hashbrown reserve_rehash (~2.5%) remaining after P1/P2: the per-read seed dedup HashSet in find_seeds used default SipHash, unsized. Switch to a pre-sized FxHashSet. Also pre-size the cluster_seeds `windows` Vec (<= one per anchor) to avoid growth reallocs. Output unchanged; 512 tests pass. Deeper per-read scratch-buffer reuse (thread-local) was evaluated and deferred: the remaining cluster_seeds cost is logic, not allocation, and a thread-local rewrite of the core aligner risks faithfulness bugs for ~5%. The next larger contained win is 2-bit genome for compare_seq_to_genome (~12%). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
P2 set purge_delay=-1 (never purge) for speed, but that makes mimalloc retain a resident free-page reserve that inflates Max RSS (reusable scratch, not data — the real heap is ~3 GB next to a reclaimable ~25 GB mmap'd index). Read an optional RUSTAR_PURGE_DELAY_MS env override (default -1) so users can cap that reserve at a small wall-time cost. No behavior change by default. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The 42 GB Max RSS on the Linux mouse-solo benchmark broke down (measured via /proc and macOS vmmap) as ~24 GB reclaimable mmap'd index + ~17 GB anonymous heap. But RssAnon was 17 GB from the first ~5 s and flat thereafter, while the actual solo records are ~170 MB and macOS shows only ~2.7 GB resident anon for the same run. The 17 GB was mimalloc's eager *arena* commit: on Linux (an overcommit OS) the default (arena_eager_commit=2) commits a large arena pool up front that is never touched. Set mi_option_arena_eager_commit=0 so arenas commit on demand. Measured (mouse 10M, 4 features + EmptyDrops): RssAnon 17 GB -> 2 GB, Max RSS 41 -> 26 GB (= ~24 GB reclaimable index + ~2 GB live heap), wall +~3% (the on-demand commit faults). Other commit knobs (eager_commit, eager_commit_delay) had no effect — it is specifically arena eager-commit. See test/aws/mem_tune.sh. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the O(genes-before-position) linear scan in GeneAnnotation::overlapping_in_into with a max-end segment-tree traversal (build_maxend_seg + explicit-stack DFS that prunes any subtree whose max exon/gene-body end is <= the query start). This was the #1 hotspot in a STARsolo-parity solo run (perf: 18.5% of samples). On the 1M-read PBMC benchmark it drops to ~2.5% and cuts solo wall time ~14% (122s -> 104.5s), narrowing the gap to STARsolo from ~1.40x to ~1.20x. Gene + GeneFull count matrices are byte-identical before/after (verified on 1M real reads); a new fuzz test cross-checks the segment-tree query against brute force over 6000 randomized differential cases. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add STAR-faithful sparse-SA support: --genomeSAsparseD=D keeps every Dth suffix over STAR's doubled fwd+RC coordinate space. STAR stores the mirrored position w = 2*nGenome-1-ii, so the keep-rule inverts the mirror stride; the resulting D=2 SA is byte-identical to STAR's D=2 SA for the yeast genome. On the mm10 PBMC index D=2 cuts peak RSS ~31% (25.7 GB -> 17.7 GB) with output-identical alignments. Also add cache-line prefetch hooks used by the seed search: PackedArray::prefetch / SuffixArray::prefetch, and MADV_RANDOM (advise_random) on the genome mmap so the kernel stops read-ahead on the random SA-driven access pattern. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
New src/align/simd_scan.rs: find_stop() scans a read/genome span for the first mismatch using arch-specific SIMD (SSE2 _mm_cmpeq_epi8 / NEON vceqq_u8, scalar fallback), fuzz-tested to 20k cases and verified on x86_64 + aarch64. seed.rs: SA prefetch in the max_mappable_length binary-search loop and a SIMD forward-only fast path in compare_seq_to_genome. stitch.rs: expand_hit returns a fixed [Option<...>; 2] (no heap Vec), FxHash for post-recursion dedup, and the per-window by_len / diag_ranges scratch buffers are hoisted to per-read scope and cleared+reused. cpu.rs: prefetch_read() (x86 _mm_prefetch, aarch64 prfm). genome/mod.rs: #[inline] on get_base. All output-identical. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ered FASTQ Restructure align_reads_single_end / align_reads_paired_end into 3-stage scoped pipelines: a producer thread decodes FASTQ, rayon aligns, and a dedicated writer thread drains a bounded sync_channel(2) via a new `trait AlignmentWriter: Send`. This overlaps decode + align + write instead of serializing them. Switch flate2 to the pure-Rust zlib-rs backend (~2-3x faster inflate/deflate on the FASTQ decode + BGZF paths, no C toolchain) and give the FASTQ reader a 512 KiB BufReader. solo/count.rs threads the pipelined reader through the solo path. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add a PGO pipeline for the 5 natively-built Linux release targets. src/bin/gen_pgo_fixture.rs synthesizes a small genome + reads fixture; scripts/pgo-build.sh runs the instrumented (profile-generate) build, aligns the fixture to collect a profile, then rebuilds the shipped binary with profile-use. release.yml gains a `pgo: true` matrix flag, the llvm-tools component, and a "Build (PGO)" step. ~10-20% faster alignment, output-identical. Restricted to Linux because PGO must run the instrumented binary natively (the macOS targets are cross-compiled on Apple Silicon runners). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add test/aws/: scripts to provision an EC2 instance (setup_aws.sh, launch.sh, bench-user-policy.json), stage genomes/reads onto NVMe instance-store (stage.sh, prep_human.sh, prep_10x.sh), and run the benchmark matrix -- old-vs-new (bench_oldnew.sh), the 3-way rustar/STARsolo/CellRanger comparison (bench_3way.sh), human WGS (bench_human_nvme.sh), memory-trajectory (mem_traj.sh, mem_tune.sh), and perf profiling (profile.sh, profile_solo.sh, profile2.sh). Plus scrape_results.py and the rendered benchmark SVGs. Also gitignore .DS_Store. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…RUSTSEC advisories Clippy (x86_64): use pointer::cast instead of `as` in prefetch_read. Audit: bump anyhow 1.0.103 (RUSTSEC-2026-0190), memmap2 0.9.11 (RUSTSEC-2026-0186), crossbeam-epoch 0.9.20 (RUSTSEC-2026-0204). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
README: add a Performance & Benchmarks section (native three-way vs STARsolo 2.7.11b and CellRanger 10.0.0 on 5k_Mouse_PBMCs_5p_gem-x_GEX), a STARsolo usage example, a Single-cell feature entry (Gene/GeneFull/SJ/ Velocyto, CB_UMI_Simple/Complex/SmartSeq, CB correction, UMI dedup, multi-mapper resolution, EmptyDrops_CR/CellRanger2.2 cell calling) and a sparse-SA entry; refresh the status line and test count (516) and replace the "no STARsolo" limitation. ROADMAP: mark Phase 14 substantially complete (Gene matrix byte-identical to STARsolo). Update the pipeline diagram, phase summary row, and the 14.5-14.11 sub-phase table (all done except 14.7 CB/UB SAM tags). Add a 2026-07 entry covering the feature set, the segment-tree gene-overlap optimization, sparse SA, and the native three-way benchmark that supersedes the old Docker-emulation numbers. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Split Performance & Benchmarks into two subsections. The existing mouse 5k-PBMC three-way (vs STARsolo/CellRanger) becomes "Mouse". Add "Human — GRCh38, memory-bound": a full snRNA-seq sample (SRR13278444, ~155M reads) whose 32 GB dense index does not stay resident in 48 GB RAM, so the run is page-fault-bound (9.5 h production run, 558M faults). --genomeSAsparseD 2 shrinks the index to 16 GB so it fits: on a 500k-read subset, 54 s vs 165-460 s dense with 16x fewer page faults, sensitivity-equal (77.6% vs 77.0% unique). Documents that on a RAM-constrained box index size, not raw aligner speed, is the dominant cost for a human genome. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
byte_char_slices: [b'A',b'C',b'G',b'T'] -> *b"ACGT" (gen_pgo_fixture x3, alignment_features), [b'S',b'M']/[b'L',b'B'] -> *b"SM"/*b"LB" (sam test). question_mark: if-let-else return None -> ? in transcriptome projection. manual_assert_eq: assert!(x == 0) -> assert_eq! in suffix_array test. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…rsolo Integration branch combining two independent streams on top of the CI/clippy fix (f7b33d1): 1. Cross-batch align pipeline — run_batch_pipeline<In,Out> driver replacing the triplicated per-batch inflight-queue plumbing in the SE/PE/solo loops (rayon::spawn + ordered forward/drain), so a slow read never idles the pool at a batch barrier. Byte-identical output. 2. --soloBarcodeMate 1 (5' cellgeni PE-solo) — barcode read from mate 1, both mates aligned as a pair, unioned genes at mate-1 strand (process_read_pe). align_reads_solo_pe is written on the same run_batch_pipeline driver as the SE-solo path (not the original 2-stage form), so both solo paths share the pipeline. Per-mate clip5p/clip3p, soloBarcodeMate/soloBarcodeReadLength/ genomeLoad params, SoloPairedReader. Not included (per scope decision): the stitch.rs recursion stall guard. Verified: 494 tests + integration, 0 clippy, fmt clean. Byte-identical to the serial baseline on SE (3178), PE (100010), and solo (all Solo.out + SAM), and byte-identical to real Linux STARsolo on both the 5' PE and 3' SE differentials. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Merging #90 into main, with some fixes and tests against STARsolo