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

.claude/CLAUDE.md

@@ -345,6 +345,8 @@ See `.claude/BARCODE_CALLING.md` for barcode calling architecture (pipeline inte
 
 See `.claude/CELL_BARCODE_SELECTION.md` for cell barcode detection (branch `badger`, counting and selection ported from https://github.com/algbio/Badger) — the answer to stock 10x whitelists, where `min_score == BARCODE_LEN_10X` makes per-read matching degenerate into exact matching. `--n_cells` decides what `--barcode_whitelist` means: unset it is the cell list (one pass, matched as today); set to a number or `auto` it is a pool, so a first pass extracts barcode windows verbatim (`TenXBarcodeDetector(whitelist_matching=False)`), the counts pick the cell barcodes, and a second ordinary pass matches reads against those. `--barcode_whitelist auto` detects without a pool at all. Measured on three real 10x datasets (ONT StereoQ, ONT cDNA R10.4, concat/split): identical to supplying the true cell list, +13 points of recall over taking the stock list at face value; `-b auto` is worse (R10.4 −5 precision) because counts alone cannot tell a cell from a recurring extraction artifact. Architecture: `CellBarcodeSelector` + `select_cell_barcodes` (`isoquant_lib/barcode_calling/cell_selection.py`), `detect_cell_barcodes` in `isoquant.py`, writing `aux/<prefix>.raw_barcodes_<i>.tsv` → `<prefix>.cell_barcodes.tsv` → `<prefix>.barcoded_reads_<i>.tsv`. `--barcode_correction` is a hidden override. Badger's edit-distance graph correction was implemented, measured against the existing SSW matcher (86.55 vs 87.35 recall given the same barcodes) and dropped: SSW slides to find a shifted window, a fixed-window edit distance cannot. Three gotchas documented there, all found by measurement: raw-mode strand selection must be structural (R1-to-polyT span) — worth 20 points of recall; `score_diff` never rejected anything because the runner-up was untracked and same-offset ties invisible; and `estimate_cell_number` reported a single cell on a flat count distribution.
 
+See `.claude/SC_IO_OUTPUTS.md` for the single-cell I/O outputs (branch `sc_outputs`) — the split-reads FASTA and the barcoded read tables are now gzipped unless `--no_gzip` (compression happens in the chunk workers and at end of run respectively, never on a path anything waits on), plus two optional BAM outputs: `--large_output tagged_bam` (`<prefix>.tagged.bam`, every input alignment with CB/UB, a pure side output) and `--large_output deduplicated_bam` (`<prefix>.deduplicated.bam`, primary-only UMI-survivor subset with CB/UB/GX/TX; deliberately **not** wired into fusion detection — UMI filtering removes the chimeric/inconsistent reads fusion calling depends on, so dedup for fusion belongs inside the fusion algorithm). Shared helpers in `isoquant_lib/utils/bam_utils.py`. Also records the `--split_molecules auto` regression fixed there: splitting with aligned (`--bam`) input produced a FASTA that was never re-aligned.
+
 See `.claude/ANALYSIS_OPTION.md` for the `--analysis` interface — the single option (values `quantification`/`quant`, `transcript_discovery`/`td`, `exon_quantification`/`ex_quant`, `fusion`) that selects pipeline stages, with context-aware defaults. It supersedes the now-hidden/deprecated `--count_exons`, `--count_intron_retentions`, `--fusion`, `--no_model_construction`, which still work. Resolution lives in `resolve_analyses()` (`isoquant.py`), producing internal flags `run_quantification` / `count_exons` / `count_intron_retentions` / `fusion` / `no_model_construction` / `predict_terminal_sites`.
 
 See `.claude/JOINT_EXON_COUNTS.md` for joint exon counts — region-based exon quantification that groups overlapping annotated exons into regions and emits N+1 features per region (N inclusion variants + 1 region-level exclusion). Runs alongside the classic `ExonCounter` when exon quantification is enabled (`--analysis exon_quantification`, aka the deprecated `--count_exons`).

docs/barcode_calling.md

@@ -96,6 +96,11 @@ matches the reads against that much shorter list. See
 
 Supported for `tenX_v3`, `tenX_v2` and `visium_5prime`.
 
+`--no_gzip`
+
+Do not compress the split reads FASTA. It is gzipped by default; aligners read it compressed,
+so there is normally no reason to turn this off.
+
 `--threads` or `-t`
 
 Number of threads for parallel processing (default: 16).
@@ -174,7 +179,11 @@ of them exactly.
 
 Output contains one row per detected molecule. Read IDs include segment coordinates:
 `{original_read_id}_{start}_{end}_{strand}`.
-An additional split FASTA file (`*.split_reads.fasta`) is produced with the extracted cDNA segments.
+An additional split FASTA file (`<prefix>.split_reads.fasta.gz`, numbered
+`<prefix>_<i>.split_reads.fasta.gz` for several inputs) is produced with the extracted cDNA
+segments, gzipped unless `--no_gzip` is set. Aligners read it compressed, so nothing downstream
+is slowed down by this. Inside the IsoQuant pipeline the same file is named
+`<prefix>.split_reads_<i>.fa.gz`.
 
 **Curio** (`curio`):
 

docs/cmd.md

@@ -179,6 +179,10 @@ original file name, and barcode property (e.g. cell type).
   * `corrected_bed` - BED file with corrected read exon coordinates (`*.corrected_reads.bed.gz`);
   * `read2transcripts` - reads assigned to discovered transcript models, in the unified read_info format (`*.transcript_model_reads.tsv.gz`);
   * `allinfo` - old format for UMI filtered reads for single-cell/spatial modes (`*.allinfo`);
+  * `tagged_bam` - single-cell/spatial modes only: a copy of the input alignments with the
+    detected barcode and UMI added as tags (`*.tagged.bam`), see below;
+  * `deduplicated_bam` - single-cell/spatial modes only: the UMI-deduplicated alignments with
+    barcode, UMI, gene and transcript tags (`*.deduplicated.bam`), see below;
   * `none` - do not generate any large output files (not compatible with other values).
 
 Example usage:
@@ -191,6 +195,38 @@ isoquant.py --large_output read_info read_assignments corrected_bed read2transcr
 isoquant.py --large_output none ...
 ```
 
+### Tagged and deduplicated BAM files
+
+Both are indexed BAM files written for single-cell and spatial modes only, and both are off by
+default. Alignment records are copied unchanged apart from the tags, so anything that reads the
+original BAM works on these too.
+
+`tagged_bam` keeps **every** alignment in the input - primary, secondary, supplementary and
+unmapped, across all references including unplaced scaffolds that IsoQuant does not analyse -
+and adds:
+
+  * `--barcode_tag` (`CB` by default) - the detected cell barcode;
+  * `--umi_tag` (`UB` by default) - the detected UMI.
+
+Unmapped reads are tagged like any other: the barcode is called from the read sequence and does
+not depend on the read having aligned. Reads with no barcode are kept without tags. The file is
+redundant when `--barcoded_bam` was used as input, since those alignments already carry the
+tags, and IsoQuant warns and skips.
+
+`deduplicated_bam` keeps only the **primary** alignments of the reads that survived UMI
+filtering - one read per detected molecule - and additionally tags:
+
+  * `GX` - the gene the read was assigned to;
+  * `TX` - the transcript the read was assigned to.
+
+A tag is omitted rather than given a placeholder value when the corresponding value is unknown,
+so novel and ambiguous reads have no `TX` tag.
+
+Note that this file is a **deduplicated** view, not a general-purpose replacement for the input
+alignments: it holds one read per detected molecule, keeps only reads assigned to a gene, and
+drops secondary and supplementary records. Analyses that depend on chimeric or otherwise
+inconsistent reads - fusion detection in particular - must keep reading the original BAM.
+
 The `read_info.tsv` format can be also converted to old formats using the conversion script:
 ```bash
 python -m isoquant_lib.scripts.convert_read_info --read_info SAMPLE.read_info.tsv.gz --format read_assignments --output SAMPLE.read_assignments.tsv
@@ -294,6 +330,11 @@ enabling even on libraries without concatenated molecules: on non-concatenated 1
 recovers about one extra point of recall at unchanged precision, at roughly twice the barcode
 calling runtime.
 
+Splitting requires raw reads. The split molecules have to be aligned, so requesting it
+alongside an already aligned input (`--bam`) is contradictory and IsoQuant stops with an error
+- under `auto` as well as `true`. Supply the reads as FASTQ/FASTA to have them split and
+mapped, or pass `--split_molecules false` to use the alignments as given.
+
 The superseded mode names `tenX_v3_split`, `tenX_v2_split` and `stereoseq_nosplit` still work
 and are translated to the corresponding `--mode` plus `--split_molecules` combination.
 
@@ -508,6 +549,12 @@ We recommend _not_ to modify these options unless you are clearly aware of their
 
 `--no_gzip`
     Do not compress large output files.
+    Compressed outputs use gzip level 6 for tables and level 4 for FASTA, which is where the
+    time/size trade-off sits for each kind of data.
+    This also covers the single-cell outputs: the split-reads FASTA (compressed as it is
+    written, in the barcode-calling workers) and the barcoded read tables (compressed once the
+    run finishes, so they stay readable while the pipeline needs them). Neither slows any
+    subsequent step down.
 
 `--no_gtf_check`
     Do not perform input GTF checks.

docs/formats.md

@@ -168,6 +168,28 @@ Each extra round only writes UMI deduplicated reads in allinfo format to
 statistics to `SAMPLE_ID.UMI_filtered.barcode_barcode_col{C}.ED{N}.stats.tsv`
 (`{C}` = 0-based spot-column index).
 
+### Tagged and deduplicated BAM
+
+Two optional BAM outputs for single-cell and spatial modes, enabled with
+`--large_output tagged_bam` and `--large_output deduplicated_bam` respectively. Both are
+indexed, and alignment records are copied from the input unchanged apart from the added tags.
+
+`SAMPLE_ID.tagged.bam` holds every input alignment - primary, secondary, supplementary and
+unmapped, across all references including unplaced scaffolds - so its record count equals the
+input's. `SAMPLE_ID.deduplicated.bam` holds only the primary alignments of the reads that
+survived UMI deduplication - the same read set as `SAMPLE_ID.UMI_filtered.ED{N}.allinfo`.
+
+| Tag | Present in | Value |
+|-----|------------|-------|
+| `CB` (or `--barcode_tag`) | both | cell barcode |
+| `UB` (or `--umi_tag`) | both | UMI |
+| `GX` | deduplicated only | assigned gene ID |
+| `TX` | deduplicated only | assigned reference transcript ID |
+
+A tag is **omitted** when the value is unknown rather than being set to a placeholder, so a
+read with no barcode carries no `CB`, and a novel or ambiguous read carries no `TX` (`allinfo`
+writes `None` in the corresponding column).
+
 
 ## Quantification formats
 

docs/output.md

@@ -185,8 +185,17 @@ By default, in single-cell and spatial modes IsoQuant only performs quantificati
 UMI-filtered reads will be saved to the same [read_info format](#read-assignments). 
 All counts formats will also be identical (see above).
 
-If IsoQuant detects the barcodes, barcoded reads will be saved in [TSV format](barcode_calling.md#output).
-If barcode calling also splits the reads into individual cDNAs, a FASTA file with cDNAs will be produced.
+If IsoQuant detects the barcodes, barcoded reads will be saved in [TSV format](barcode_calling.md#output),
+gzipped once the run finishes unless `--no_gzip` is set.
+If barcode calling also splits the reads into individual cDNAs, a FASTA file with cDNAs will be
+produced, also gzipped unless `--no_gzip` is set.
+
+Two optional BAM outputs carry the same information on the alignments themselves
+(see [`--large_output`](cmd.md#--large_output) for details):
+
+* `SAMPLE_ID.tagged.bam` - all input alignments with barcode and UMI tags (only with `--large_output tagged_bam`);
+* `SAMPLE_ID.deduplicated.bam` - primary alignments of UMI-deduplicated reads, with barcode, UMI,
+  gene and transcript tags (only with `--large_output deduplicated_bam`).
 
 Note that transcript discovery is performed only in `bulk` mode by default.
 Single-cell and spatial modes require UMI deduplication.

docs/single_cell.md

@@ -347,10 +347,17 @@ Supported for `tenX_v3`, `tenX_v2`, `stereoseq` and `visium_5prime`. `--split_mo
 with any other mode is an error rather than a silent no-op, so a request that cannot be
 honoured never passes unnoticed.
 
-When splitting, IsoQuant writes an additional output file (`*.split_reads.fasta`) containing
+When splitting, IsoQuant writes an additional output file (`*.split_reads_<i>.fa.gz`) containing
 the extracted cDNA segments, and uses it in place of the original reads for alignment. Each
 segment is named with the original read ID plus coordinates and strand:
-`{read_id}_{start}_{end}_{strand}`.
+`{read_id}_{start}_{end}_{strand}`. The file is gzipped unless `--no_gzip` is set; minimap2
+reads it compressed, so alignment is not slowed down.
+
+Splitting rewrites the reads, so the pieces have to be aligned afresh. Supplying an aligned
+BAM (`--bam`) says the opposite - that no mapping should happen - so the two requests
+contradict each other and IsoQuant aborts rather than guessing. This applies to `auto` as well
+as `true`: pass the raw reads if you want the molecules split and re-aligned, or
+`--split_molecules false` to analyse the alignments as they are.
 
 Splitting is worth leaving on even for libraries you do not expect to be concatenated: measured
 on non-concatenated 10x data it recovers about one extra point of recall at unchanged precision

isoquant.py

@@ -77,7 +77,8 @@
 logger = logging.getLogger('IsoQuant')
 
 # Large output file types for --large_output option
-LARGE_OUTPUT_TYPES = ["read_info", "read_assignments", "corrected_bed", "read2transcripts", "allinfo", "none"]
+LARGE_OUTPUT_TYPES = ["read_info", "read_assignments", "corrected_bed", "read2transcripts", "allinfo",
+                      "tagged_bam", "deduplicated_bam", "none"]
 
 
 def bool_str(s):
@@ -579,11 +580,24 @@ def check_and_load_args(args, parser):
             if val not in LARGE_OUTPUT_TYPES:
                 logger.error("Invalid --large_output value: %s. Valid values: %s" % (val, ", ".join(LARGE_OUTPUT_TYPES)))
                 sys.exit(IsoQuantExitCode.INVALID_PARAMETER)
+        _warn_about_unusable_bam_outputs(args)
 
     save_params(args)
     return args
 
 
+def _warn_about_unusable_bam_outputs(args):
+    """Say up front when a requested BAM output cannot be produced, rather than mid-run."""
+    if not args.mode.needs_pcr_deduplication():
+        for output_type in ("tagged_bam", "deduplicated_bam"):
+            if output_type in args.large_output:
+                logger.warning("--large_output %s has no effect in %s mode, which has no barcodes "
+                               "or UMIs; it will be skipped" % (output_type, args.mode.name))
+    elif "tagged_bam" in args.large_output and getattr(args, 'barcoded_bam', False):
+        logger.warning("--large_output tagged_bam is redundant with --barcoded_bam: those "
+                       "alignments already carry the tags; it will be skipped")
+
+
 def load_previous_run(args):
     logger.info("Loading parameters from the previous run")
     logger.error("Only --output/--threads/--debug/--high_memory are compatible with --resume option")
@@ -752,6 +766,24 @@ def _dedup_read_group_specs(args):
         args.read_group = updated_specs
 
 
+def _reject_splitting_aligned_input(args):
+    """Refuse to split molecules when the reads are already aligned.
+
+    Splitting rewrites each read into its constituent cDNAs, so the pieces have to be aligned
+    afresh -- but supplying a BAM says "do not map". The two requests contradict each other,
+    and either way of guessing silently gives the user something they did not ask for, so make
+    them choose.
+    """
+    if not args.split_molecules or args.input_data.input_type.needs_mapping():
+        return
+    logger.critical("Reads cannot be split into separate molecules when they are already "
+                    "aligned (%s input): the split molecules would have to be mapped again. "
+                    "Provide the raw reads instead, or use --split_molecules %s to analyse "
+                    "the alignments as they are."
+                    % (args.input_data.input_type.name, SPLIT_MOLECULES_FALSE))
+    sys.exit(IsoQuantExitCode.INCOMPATIBLE_OPTIONS)
+
+
 def check_input_params(args):
     if not _validate_data_type_and_input(args):
         return False
@@ -762,6 +794,7 @@ def check_input_params(args):
         resolve_deprecated_mode(args)
         args.mode = IsoQuantMode[args.mode]
     resolve_split_molecules(args)
+    _reject_splitting_aligned_input(args)
 
     # translate --analysis (and the deprecated stage flags) into internal booleans
     resolve_analyses(args)

isoquant_detect_barcodes.py

@@ -83,6 +83,8 @@ def add_hidden_option(*args, **kwargs):  # show command only with --full-help
     add_hidden_option("--n_cells_interval", type=int, default=25)
     add_hidden_option("--barcode_correction", type=str, choices=[e.name for e in BarcodeCorrectionMethod],
                       default=BarcodeCorrectionMethod.auto.name)
+    parser.add_argument("--no_gzip", help="do not gzip the split reads FASTA",
+                        dest="gzipped", action='store_false', default=True)
     add_hidden_option('--debug', action='store_true', default=False, help='debug log output.')
 
     args = parser.parse_args(sys_argv)
@@ -103,10 +105,11 @@ def check_args(args):
             args.output_tsv = [args.output + "_%d.barcoded_reads.tsv" % i for i in range(num_files)]
 
     if args.out_fasta is None and args.split_molecules:
+        suffix = ".split_reads.fasta.gz" if args.gzipped else ".split_reads.fasta"
         if num_files == 1:
-            args.out_fasta = [args.output + ".split_reads.fasta"]
+            args.out_fasta = [args.output + suffix]
         else:
-            args.out_fasta = [args.output + "_%d.split_reads.fasta" % i for i in range(num_files)]
+            args.out_fasta = [args.output + "_%d%s" % (i, suffix) for i in range(num_files)]
 
 
 def run_barcode_calling(args):

isoquant_lib/assignment/assignment_io.py

@@ -6,7 +6,6 @@
 ############################################################################
 
 import logging
-import gzip
 
 from isoquant_lib.common import (
     CANONICAL_FWD_SITES,
@@ -17,6 +16,7 @@
     sum_intervals_from_point,
     sum_intervals_to_point
 )
+from isoquant_lib.utils.file_utils import open_text_write
 from isoquant_lib.utils.serialization import (
     write_short_int,
     read_short_int,
@@ -72,9 +72,7 @@ def __init__(self, output_file_name, params, assignment_checker=PrintAllFunctor(
         self.gzipped = gzipped
         if gzipped:
             self.output_file_name += ".gz"
-            self.output_file = gzip.open(self.output_file_name, "wt")
-        else:
-            self.output_file = open(self.output_file_name, "w")
+        self.output_file = open_text_write(self.output_file_name)
 
     def __del__(self):
         self.output_file.close()

isoquant_lib/assignment/assignment_loader.py

@@ -216,7 +216,8 @@ def prepare_read_filter(chr_id, saves_prefix, use_filtered_reads):
         return None
     filtered_reads = set()
     for line in open(filtered_reads_file_name(saves_prefix, chr_id), "r"):
-        filtered_reads.add(line.rstrip())
+        # the file carries tag columns after the read id when deduplicated_bam is requested
+        filtered_reads.add(line.rstrip().split("\t")[0])
     return filtered_reads
 
 

isoquant_lib/assignment/read_groups.py

@@ -12,6 +12,7 @@
 
 import pysam
 
+from isoquant_lib.utils.file_utils import strip_compression_suffix
 from isoquant_lib.utils.error_codes import IsoQuantExitCode
 from isoquant_lib.utils.table_splitter import split_read_table_parallel
 
@@ -86,7 +87,7 @@ def __init__(self, args, sample):
         self.readable_names_dict = {}
         for sample in args.input_data.samples:
             for lib in sample.file_list:
-                readable_name = os.path.splitext(os.path.basename(lib[0]))[0]
+                readable_name = os.path.splitext(os.path.basename(strip_compression_suffix(lib[0])))[0]
                 for f in lib:
                     self.readable_names_dict[f] = readable_name