Single-cell I/O: compress intermediates, add tagged and deduplicated BAM outputs

[andrewprzh]

Four independent single-cell I/O changes: two reclaim disk on existing outputs, two add BAM
outputs carrying information the pipeline already computes but only emitted as text.

1. The split-reads FASTA is gzipped

<prefix>.split_reads_<i>.fa.gz unless --no_gzip. Compression happens inside the chunk
workers, and the merge stays a plain byte copy because concatenated gzip members form a valid
stream. Nothing downstream slows down: the only consumer is minimap2, which decompresses
natively.

2. The barcoded read tables are gzipped at the end of the run

They stay plain for the whole run, since split_read_barcode_table and the tagged-BAM build
read them, and are compressed once every sample is done. The per-worker split tables under
aux/ are deliberately left alone: short-lived temporaries whose compression would cost CPU in
every chromosome worker.

3. --large_output tagged_bam (off by default)

<prefix>.tagged.bam — a copy of the input alignments with --barcode_tag / --umi_tag
(CB/UB). Every alignment is kept: primary, secondary, supplementary and unmapped, across all
references including unplaced scaffolds. Purely a side output; nothing downstream reads it.

Skipped with a warning when --barcoded_bam was the input (those alignments already carry the
tags) or in modes with no barcodes.

4. --large_output deduplicated_bam (off by default)

<prefix>.deduplicated.bam — primary alignments of the reads that survived UMI filtering, one
per detected molecule, tagged with barcode, UMI, GX (gene) and TX (transcript). All four
values come from the ReadAssignment already in hand where survivors are chosen, so no extra
pass over the assignments is needed.

Deliberately not wired into fusion detection: UMI filtering keeps one read per assigned
(gene, barcode, UMI) molecule and requires a gene assignment, so the chimeric reads fusion
calling depends on are collapsed or dropped. Deduplication for fusion belongs inside the fusion
algorithm.

Also here

--split_molecules with an already aligned input (--bam) now aborts instead of producing a
FASTA that is never re-aligned. Splitting rewrites the reads so the pieces must be mapped
again, while a BAM says not to map — the requests contradict, so IsoQuant asks the user to
choose rather than guessing. This applies to auto as well as true. Four CI configs that
pair --bam with a splitting-capable mode now pass --split_molecules false; all four hit the
crash before this.

--read_group file_name on a gzipped input changes label: reads.fq.gz grouped as reads.fq
before (one splitext off a two-part extension) and groups as reads now. The split-reads
FASTA needed this to keep a sane group name once it gained .gz, and FileNameGrouper and
StringPoolManager were changed together so they still agree. No CI baseline moves —
STEREO.TOY is the only config pairing file_name with a gzipped input and it asserts on file
existence, not contents.

Review pass

Follow-ups on top of the original branch (163aafde..9dea1a9e):

• The merge could not have worked on a many-contig reference. One fragment per non-empty
  reference plus a single samtools merge call means every fragment is opened at once, and
  samtools gives up just above a thousand handles no matter what ulimit -n says — measured:
  1200 fragments fail at fragment 1019 with RLIMIT_NOFILE at 1048576. GRCh38's full analysis
  set has 3366 contigs, so tagged_bam on it would have crashed at the very end of the run.
  Merging now goes in batches of 500 (_merge_in_rounds), verified on the same 1200 fragments:
  0.6 s, all records present, no intermediates left behind.

• The unmapped-read pass read every BAM twice, end to end, in the parent. Unplaced reads
  sit at the tail of a coordinate-sorted file and fetch(until_eof=True) walked the whole file
  to reach them, once to collect ids and once to copy. unplaced_reads() now uses fetch("*")
  on an indexed BAM, with the linear scan kept as the fallback; a test pins that both paths
  return the same records.

• Compression level, and it now covers the whole project. Both new writers used Python's
  gzip default of 9. Measured on real ONT outputs, the two kinds of data want different
  levels, so gzip_level_for(name) picks one from the file's extension:
  tables at 6 (a real barcode table runs 16.0 MB/s at level 9 against 39.2 MB/s at 6 for
  3.0% less output; going below 6 was measured and rejected — level 1 is 2.4x faster again but
  8.5% larger, levels 4-5 save almost nothing) and FASTA/FASTQ at 4 (nucleotide data sits
  near gzip's entropy floor, so level 6 grinds at 12.5 MB/s against 65.3 MB/s at 4 for 5.7%
  more output).

  Every gzip writer in the project now goes through open_text_write / gzip_file_in_place
  rather than calling gzip.open itself, so both levels apply everywhere: the allinfo writers,
  the four in convert_grouped_counts, TextFileAssignmentPrinter (read_info /
  read_assignments / corrected_bed / read2transcripts) and the two scripts/ converters.

  One defect fell out of checking the output size rather than the code: the split FASTA is
  compressed in the chunk workers and those temps were named subreads.gz with no .fa, so
  they took the table level while the single-threaded path, which passes the real output
  name, took the sequence level — the level differed by thread count. numbered_chunk_name
  now keeps the whole extension chain last (subreads.fa.gz → subreads_3.fa.gz). Verified
  end to end: 5.28 MB → 5.58 MB, exactly the level-4 number, with identical read-id sets.

• The tagged BAM was rebuilt on every --resume. It is a full copy of the input and it is
  written before read collection, so a crash in the long stages meant re-copying the whole BAM
  on each resume attempt. It now has a resume marker in aux/ (guarded by an existence check
  on the BAM), matching how the other expensive stages behave.

• The reference list is derived once and passed into write_tagged_bam rather than
  recomputed. It has to be the same list that drove the barcode-table split or a fragment finds
  no table and its reads come out untagged; the invariant is now explicit in the signature.
  The split dict also builds its paths with sample.get_barcodes_split_file(), the same helper
  the workers read them back with — the two disagreed for a reference name containing /.

• Barcode calling output order depended on which chunk finished first.
  run_chunks_in_parallel waits on FIRST_COMPLETED and the merging caller appended chunk
  temp names in that order, so the row order of barcoded_reads.tsv and the record order of
  the split FASTA varied run to run — two runs of identical code on identical input differed
  at line 200003, a chunk boundary. Row sets were always identical and every consumer keys by
  read id, so nothing was wrong, but the outputs were not reproducible. handle_result now
  receives (result, chunk_index); the merging caller keys a dict on the index and sorts at
  merge time, and the counting caller ignores it (its selector already sorts by
  (count, barcode), a total order). Scheduling is untouched — waiting for chunks in order
  would idle the pool behind one slow chunk. Two runs now give a byte-identical barcode table
  and a byte-identical decompressed FASTA; the compressed FASTA differs in exactly 10 bytes,
  all gzip header mtime fields, so compare these outputs decompressed.

Also fixed: docs/barcode_calling.md documents the standalone isoquant_detect_barcodes.py,
but named the pipeline's *.split_reads_<i>.fa.gz; the standalone tool writes
<prefix>.split_reads.fasta.gz. Its --no_gzip was undocumented.

Verification

Unit: 1524 tests pass (console tests need minimap2 on PATH), including test_bam_utils.py (35)
and test_file_compression.py (18).

End-to-end on the full SC.Mouse.10x.allinfo CI dataset (2M reads, 3660591 records), which is
where three tagged_bam defects surfaced that a single-chromosome subset could not:

• tagged.bam — 3660591 records, equal to the input; 1655717 secondary, 4874 supplementary and
  3577 unmapped all preserved; CB/UB match the barcode table on all 1794571 barcoded reads
  with zero differences.
• deduplicated.bam — 903570 records, all primary, exactly the read set of
  UMI_filtered.ED3.allinfo; CB/UB/GX match allinfo on every row; records byte-identical to
  the input apart from the tags, SA tags preserved.
• Same checks pass on SC.Mouse.10x.barcoded_bam.allinfo (898640 records), which exercises the
  dedup BAM on the tags-read-from-BAM path.
• Both are pure side outputs: every other file is byte-identical to a run without the flags.
• allinfo baselines pass on every run.

Re-verified after the review pass on chr19 + all 15 unplaced scaffolds + the unmapped tail
(120595 records): tagged BAM 120595 records with secondary/supplementary/unmapped counts
matching the input exactly, CB/UB correct on all 53760 barcoded alignments (0 missing, 0
mismatched, 0 spurious) including 3193 tagged unmapped reads; deduplicated BAM 25056 records,
all primary, read set identical to allinfo with 0 GX and 0 CB mismatches. Resume skips the
tagged BAM with its mtime unmoved. The standalone splitter still emits a valid multi-member
out.split_reads.fasta.gz.

Note that CI produces these BAMs but does not assert on them — the baselines only cover
allinfo, so a crash or knock-on regression is caught but a wrong tag value would not be.

🤖 Generated with Claude Code