diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 05f83dfd..b15ea681 100644 --- a/.claude/CLAUDE.md +++ b/.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/.raw_barcodes_.tsv` → `.cell_barcodes.tsv` → `.barcoded_reads_.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` (`.tagged.bam`, every input alignment with CB/UB, a pure side output) and `--large_output deduplicated_bam` (`.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`). diff --git a/.claude/SC_IO_OUTPUTS.md b/.claude/SC_IO_OUTPUTS.md new file mode 100644 index 00000000..71c3f175 --- /dev/null +++ b/.claude/SC_IO_OUTPUTS.md @@ -0,0 +1,347 @@ +# Single-cell I/O outputs + +Branch `sc_outputs` (based on `badger`, PR #420). Four independent changes to what the +single-cell pipeline writes to disk: two reclaim space on existing outputs, two add BAM +outputs carrying information the pipeline already computes but previously only emitted as text. + +## 1. The split-reads FASTA is gzipped + +`.split_reads_.fa.gz` unless `--no_gzip` (dest `args.gzipped`, default `True`). +The name is built in `call_barcodes` (`isoquant.py`), which picks the suffix from `args.gzipped`. + +Compression happens **inside the chunk workers**, not serially in the parent: +`detect_barcodes.py` writes each per-chunk temp through `open_text_write`, and the merge that +follows is a plain byte copy. That works because concatenated gzip members form a valid gzip +stream — the merge does not decompress and recompress. `numbered_chunk_name()` keeps `.gz` +last so the temps are recognised as compressed. + +Nothing downstream slows down: the only consumer is minimap2, which gets a path on the command +line and decompresses natively. No Python code opens this file. + +Two things the extension change touched, both fixed: + +- `string_pools.py` and `read_groups.py` derive file-name read groups with a single + `os.path.splitext`, so `.fa.gz` would have left a trailing `.fa` in the group name. Both now + call `strip_compression_suffix()` first. `test_file_compression.py` pins this. +- `read_mapper.py` derives the BAM name and already handled `.fa.gz` correctly (verified, not + changed). + +## 2. The barcoded read tables are gzipped at the end of the run + +`.barcoded_reads_.tsv` (`sample.barcodes_tsv`) stays **plain for the whole run** — +`split_read_barcode_table` reads it, and so does the tagged-BAM build. `compress_barcode_tables` +runs in `process_all_samples` after every sample is done and before `clean_up()`. + +The per-worker split tables under `aux/` (`sample.barcodes_split_reads + "_"`) are +deliberately **left uncompressed**: short-lived temporaries read back by Python, so compressing +them would cost CPU in every chromosome worker for a transient saving. + +Resume hazard, handled: `call_barcodes` short-circuits on the `barcodes_done` markers and +re-populates `sample.barcoded_reads` with the *uncompressed* names, which no longer exist after +compression. Readers go through `resolve_optionally_gzipped()`, which returns whichever of +`` / `.gz` exists. + +Helpers live in `isoquant_lib/utils/file_utils.py`: `open_text_write`, `open_text_read`, +`resolve_optionally_gzipped`, `gzip_file_in_place`, `strip_compression_suffix`. + +Both writers pass `compresslevel=GZIP_LEVEL` (6). Python's `gzip` defaults to 9. Measured on a +real 25 MB barcode table from 411K ONT reads (UUID read ids, so realistic entropy — an earlier +synthetic table was near-incompressible and badly overstated the gap): + +| level | barcode table | vs L6 size | split FASTA | vs L6 size | +|-------|---------------|-----------|-------------|-----------| +| 1 | 94.6 MB/s | +8.5% | 92.5 MB/s | +13.8% | +| 4 | 55.8 MB/s | +2.0% | 65.3 MB/s | +5.7% | +| 5 | 44.2 MB/s | +0.4% | 32.5 MB/s | +3.4% | +| 6 | 39.2 MB/s | — | 12.5 MB/s | — | +| 9 | 16.0 MB/s | −3.0% | 2.7 MB/s | −3.9% | + +So the two kinds of output take different levels, picked by `gzip_level_for(name)` from the +file's own extension: + +- **tables** (TSV, BED, MTX, allinfo) — `GZIP_LEVEL = 6`. Going below is not worth it: + `compress_barcode_tables` is the one serial compressor, but at 39 MB/s a ~100 GB table + (roughly a billion reads) costs ~45 min against ~105 min at level 9, and 6→1 would save + another ~25 min for 8.5% more disk while 4-5 save almost nothing. +- **sequences** (FASTA/FASTQ) — `GZIP_LEVEL_SEQUENCES = 4`. Nucleotide data sits near gzip's + entropy floor so the high levels grind for nothing: 5.2x the throughput for 5.7% more output. + +Every gzip *writer* in the project now goes through `open_text_write` / `gzip_file_in_place`, +so both levels apply everywhere: the allinfo writers in `dataset_processor`, the four in +`convert_grouped_counts`, `TextFileAssignmentPrinter` (read_info / read_assignments / +corrected_bed / read2transcripts), and the two `scripts/` converters. None of them call +`gzip.open` any more. + +### Barcode calling output order + +`run_chunks_in_parallel` waits on `FIRST_COMPLETED` and handed results to `handle_result` in +**completion order**, and `_process_single_file_in_parallel` appended the chunk temp names to a +list in that order — so which chunk finished first decided the row order of +`barcoded_reads.tsv` and the record order of the split FASTA. Two runs of identical code on +identical input differed (measured: at line 200003, a chunk boundary). Row sets were always the +same and every consumer keys by read id, so nothing was wrong; the outputs simply were not +reproducible. + +`handle_result` now takes `(result, chunk_index)`. The merging caller keys a dict on the index +and sorts once at merge time; the counting caller ignores it, being order-insensitive already +(`CellBarcodeSelector.sorted_barcodes` sorts by `(count, barcode)`, a total order). Scheduling +is untouched — waiting for chunks *in order* instead would idle the pool behind one slow chunk, +which is the property `run_chunks_in_parallel` exists to provide. + +Verified: two runs now produce a byte-identical barcode table and a byte-identical decompressed +FASTA, with the same row and read-id sets as before the change. The compressed FASTA still +differs in exactly 10 bytes, all of them gzip header `mtime` fields — the deflate streams are +identical. Zeroing that field was tried and **dropped**: `gzip.open()` cannot set it, so the +writers would have to build `gzip.GzipFile` objects by hand, which is more machinery than +byte-comparable `.gz` outputs are worth. Compare them decompressed. + +Only barcode calling had this pattern; every other parallel stage uses `proc.map`, which +preserves input order. + +**Gotcha, caught only by measuring the output size.** The split FASTA is compressed in the +chunk workers, and those temps were called `subreads.gz` — no `.fa`, so the inference gave them +the *table* level while the single-threaded path, which passes the real output name, gave them +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`) and the temp is named for +what it holds. Verified end to end: the FASTA went 5.28 MB → 5.58 MB, exactly the level-4 +number, with identical read-id sets. + +The table has to stay plain during the run for a second reason worth recording: +`split_read_table_parallel` streams it line by line in *every* worker, so a gzipped table would +be decompressed once per thread rather than read once. + +## 3. `--large_output tagged_bam` + +`.tagged.bam` — a copy of the input alignments with `--barcode_tag` / `--umi_tag` +(CB/UB) added. Off by default. + +Named `tagged_bam` rather than `barcoded_bam` because `--barcoded_bam` is already an *input* +flag meaning the opposite direction. + +**Purely a side output.** The barcode-table split happens regardless, and nothing downstream +reads the result — building this BAM needs those split tables in the first place, so there was +never a split to save by reusing its tags. + +Built in `DatasetProcessor.write_tagged_bam`, right after the split-table block in +`process_sample` while the tables still exist. One fragment per chromosome via +`map_over_chromosomes(write_tagged_bam_in_parallel, ...)`, then merged and indexed. + +Guarded by its own resume marker, `tagged_bam_lock_filename(sample)` (in `aux/`, so it outlives +`clean_up`) plus an existence check on the BAM itself. This output is a full copy of the input, +the most expensive thing on the branch, and it is written *before* read collection — without +the marker every `--resume` after a crash in the long stages re-copied the whole BAM. The +marker follows the `barcodes_done` precedent and is deliberately not deleted at the end of +`process_sample`. + +The reference list is computed **once** in `process_sample` and passed into `write_tagged_bam`. +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; computing it twice left that invariant implicit. + +Every alignment is kept — primary, secondary, supplementary — plus a separate +`write_unmapped_bam` pass, because `fetch(chr)` never returns unmapped reads. + +The chromosome list comes from `references_with_alignments(bam_files)`, **not** from +`get_chr_list()`. This was a bug found only on the full CI dataset: IsoQuant analyses the +22 assembled mouse chromosomes, so iterating the analysed list silently dropped 1177 records +sitting on unplaced scaffolds (`GL456382.1`, `JH584299.1`, …) — in a file documented as a copy +of the input. The chr19 subset used during development contained no scaffolds and could not +surface it. Empty references are filtered out via the BAM index so a fragmented assembly does +not spawn a task per contig. + +Restoring the records is only half of it: they came back *untagged*, because the +per-chromosome barcode split is driven by `split_barcodes_dict`, which was keyed on the +analysed chromosomes too. Barcode calling runs over the whole input before any chromosome +filtering, so the tags exist in `.barcoded_reads_.tsv` — they simply never reached a +split file the tagging worker could read. `process_sample` now widens that dict to +`references_with_alignments(...)` when `tagged_bam` is enabled (and only then, since it costs +an extra pass over the barcode table). + +### Unmapped reads were never tagged + +A separate bug, on the same output, found the same way. `write_unmapped_bam` copied unplaced +records verbatim with **no tags at all**, so on the CI dataset 3193 of the 3577 unmapped reads +lost barcodes they genuinely had. A barcode is called from the read sequence and does not +depend on the read aligning anywhere, so those tags are meaningful — arguably more so, since +recovering unmapped reads per cell is a reason to want this file. + +They belong to no chromosome and so appear in no split table. `write_tagged_bam` therefore +collects the unmapped read ids first (`collect_unmapped_read_ids`, a few thousand) and scans +the whole-sample barcode table for just those (`load_barcode_umi_tags(..., read_ids=...)`), +which keeps memory at the size of the unmapped set rather than the 1.8M-row table. + +Worth recording how this was nearly missed: after the scaffold fix the tag-mismatch count +stayed at exactly 3193, and the first reading was that the scaffolds explained it — the +arithmetic `1794571 - 1791378 = 3193` matched. It was a coincidence of two unrelated numbers. +Checking which read ids were actually untagged showed all 3193 were unmapped, and that the +scaffold records had contributed no new barcoded read ids at all. + +With `--barcoded_bam` as input there are no split tables and the input already carries the +tags, so IsoQuant warns and skips; likewise in modes with no barcodes at all. + +## 4. `--large_output deduplicated_bam` + +`.deduplicated.bam` — **primary alignments only**, restricted to the reads that survived +UMI filtering, tagged with barcode, UMI, `GX` (gene) and `TX` (transcript). Off by default. + +All four values are in hand exactly where survivors are chosen, so no extra pass over the +assignments is needed: `UMIFilter._process_chunk` already holds the `ReadAssignment`. +`UMIFilter._survivor_record` writes them as extra tab-separated columns on the existing +per-chromosome survivors file (`/.save_filtered_`) when `output_read_tags` +is set — which `parallel_workers.filter_umis_in_parallel` derives from +`large_output_enabled(args, "deduplicated_bam")`. + +Its one other consumer, `prepare_read_filter` (`assignment_loader.py`), now takes +`line.split("\t")[0]`, which makes it tolerant of both the bare-id and the tagged format. + +Two ordering constraints, both respected in `process_sample`: + +- built **after** `filter_umis` and **before** `clean_up`, which deletes `out_raw_file + "_*"` + including the survivors files; +- only the first edit distance writes those files, and the `barcode2barcode` rounds never do, + so the subset is defined by the primary dedup round. + +### Deliberately *not* wired into fusion detection + +An earlier version of this branch auto-enabled `deduplicated_bam` for fusion runs and fed it to +`FusionDetector` in place of the original BAM, reasoning that PCR duplicates would otherwise +inflate breakpoint support. **That was wrong and has been reverted**: fusion evidence lives in precisely the reads UMI filtering removes. The filter keeps +one read per assigned (gene, barcode, UMI) molecule and requires a gene assignment, so chimeric +reads spanning two genes — the inconsistent reads fusion calling is built on — are collapsed or +dropped outright. + +(The secondary/supplementary dimension was in fact fine: `fusion_detector.py:470` skips both and +takes breakpoints from the SA tag on the primary record. That is not what makes the subset +unsuitable — the read-level filtering is.) + +Deduplicating for fusion, if wanted at all, belongs **inside** the fusion algorithm where it can +see the chimeric reads before they are filtered. Fusion detection therefore reads the original +BAMs, exactly as it did before this branch — `get_bam_files_from_samples` is untouched. + +## Shared code + +`isoquant_lib/utils/bam_utils.py`: + +- `index_bam` — BAI with a CSI fallback for references too long for BAI +- `merge_bam_files` — merge per-chromosome fragments and index; single fragment is moved, not + merged; `None` and missing fragments are skipped (chromosomes with nothing to write). + Merging goes through `_merge_in_rounds` in batches of `BAM_MERGE_BATCH` (500): samtools opens + every input at once and gives up just above a thousand handles regardless of `ulimit -n` + (measured: 1200 fragments fail at fragment 1019 with `RLIMIT_NOFILE` at 1048576). One + fragment per non-empty reference means GRCh38's full analysis set, at 3366 contigs, would + have crashed the merge. A lone leftover batch is carried into the next round rather than + copied; the intermediates are deleted once the final merge succeeds. +- `unplaced_reads` — the unmapped records with no coordinates. `fetch("*")` seeks straight to + them on an indexed BAM, with a `fetch(until_eof=True)` filter as the fallback. The two + callers (`collect_unmapped_read_ids` and `write_unmapped_bam`) previously each scanned the + entire BAM linearly just to reach the tail, both serially in the parent +- `write_tagged_chromosome_bam` — the one copy loop, parameterised by `primary_only` and + `keep_untagged`; feature 3 uses `(False, True)`, feature 4 uses `(True, False)` +- `references_with_alignments` — every reference carrying reads, empty ones dropped via the + index; what the tagged BAM iterates instead of the analysed chromosome list +- `write_unmapped_bam` +- `load_survivor_tags` / `load_barcode_umi_tags` — the two tag sources +- `PLACEHOLDERS = {"*", ".", "None"}` — a tag is **omitted** rather than carrying a + placeholder. `"None"` is in there because IsoQuant stores the literal string `"None"` as + `assigned_transcript` for novel and ambiguous reads (visible in `allinfo` too), and + `TX:Z:None` would be a trap for anything reading the tag. + +`DatasetProcessor.map_over_chromosomes(worker, sample, *extra)` factors out the +`ProcessPoolExecutor` boilerplate the per-chromosome stages all repeat. + +## Regression fixed along the way + +`--split_molecules auto` (default since the `badger` branch) made every splitting-capable mode +split, including when input was an **aligned BAM**. Splitting rewrites the reads, so the pieces +need re-aligning, but `--bam` input skips the mapping stage — `call_barcodes` replaced +`sample.file_list` with the FASTA and `get_chromosome_ids` then tried to open it as a BAM +(`ValueError: file has no sequences defined`). This broke `-m tenX_v3 --bam ...`, which the +`SC.Mouse.10x.allinfo` CI config uses. + +Fixed by `_reject_splitting_aligned_input` in `isoquant.py`, run right after +`resolve_split_molecules`: splitting with aligned input aborts, under `auto` as well as `true`. +The check lives in the pipeline rather than in `options.resolve_split_molecules` because that +module is barcode-calling helpers, shared with the standalone `isoquant_detect_barcodes.py`, +which has no `input_data` at all. + +Aborting on `auto` rather than quietly not splitting is deliberate: passing a BAM says "do not +map", asking for splitting says "rewrite the reads", and either way of guessing hands the user +something they did not ask for. Four CI configs (`SC.Mouse.10x.allinfo`, +`SC.Mouse.10x.barcoded_bam.allinfo`, `GROUP12/13.SC.SIRVs.R10`) pair `--bam` with a splitting +mode and now carry an explicit `--split_molecules false`; all four were already broken by the +crash before this. + +## Verification performed + +Unit: `isoquant_tests/test_bam_utils.py` (35 tests), `isoquant_tests/test_file_compression.py` +(18 tests). The merge test monkeypatches `BAM_MERGE_BATCH` to 3 and runs 11/12/13 fragments, +straddling the batch boundary where a lone leftover has to be carried forward; a separate test +pins that `unplaced_reads` returns the same records with and without an index. + +End-to-end on chr19 of `Mouse.10x.5k.ONT_cDNA.R10.4.no_trunc.bam` (115841 records, 55392 +primary), `-m tenX_v3` with the 5K whitelist — and then on the **full** CI dataset via +`SC.Mouse.10x.allinfo`, which is what caught the scaffold bug above: + +- `deduplicated_bam`: 25056 records, 0 secondary/supplementary; read-id set is **exactly** the + 25056 rows of `allinfo`; CB/UB/GX all match `allinfo` on every row; records byte-identical to + the input apart from tags; SA tags preserved. +- `tagged_bam`: 115841 records — identical to the input count, secondary and supplementary + included; CB/UB match the barcode table exactly on all 49756 barcoded reads. +- Both are **pure side outputs**: every other file in the run is byte-identical to a run + without the flag, apart from two `.gz` files whose embedded mtime differs (decompressed + content identical). +- Fusion: verified untouched — `--analysis fusion` neither enables nor reads either BAM, and + the fusion code path is byte-identical to the branch base. + +On the full CI dataset (2M reads, 3660591 records) via `SC.Mouse.10x.allinfo`: + +- `deduplicated_bam`: 903570 records, all primary, exactly matching that run's own + `Total reads saved`. +- `tagged_bam`, after all three fixes: 3660591 records (exact), 1655717 secondary and 4874 + supplementary preserved, 3577 unmapped preserved of which 3193 barcoded, and CB/UB matching + the barcode table on **all 1794571** barcoded reads with 0 differences. +- Tag *values* were never wrong at any stage — 0 conflicting and 0 extra tags throughout; every + defect was a read or a tag going missing, never a wrong one. +- The `allinfo` baselines pass within tolerance on every run. + +### CI coverage + +`SC.Mouse.10x.allinfo` requests `tagged_bam deduplicated_bam`; `SC.Mouse.10x.barcoded_bam.allinfo` +requests `deduplicated_bam` only — `tagged_bam` is deliberately omitted there because +`--barcoded_bam` skips the split-table block it reads from, so it would warn and skip. That +second config is worth having because it exercises the dedup BAM on the tags-read-from-BAM +path, where the tags never pass through a barcode table at all. + +Note that CI *produces* these BAMs but does not assert anything about them — the baselines only +cover `allinfo`. A crash or a knock-on regression would be caught; a wrong tag value would not. + +### File-name read groups + +`strip_compression_suffix` was added for the `.fa.gz` split reads, but it changes +`--read_group file_name` for **any** gzipped input: `reads.fq.gz` used to group as `reads.fq` +(one `splitext` off a two-part extension) and now groups as `reads`. That is the intended +name, and `FileNameGrouper` and `StringPoolManager` were changed together so they still agree, +but it is a visible change to group labels for existing bulk runs. No CI baseline moves: +`STEREO.TOY` is the only config pairing `file_name` with a gzipped input and it is +`run_type: void` (checks that the grouped count files exist, not their contents). + +### Resume + +`--resume` restores `large_output` from the pickled `.params` (the resume parser rejects it on +the command line), so the survivors-file format is always consistent within a run and +`load_survivor_tags` can never meet a format the run did not write. + +Resuming an **interrupted** run works, and the tagged BAM is now skipped rather than rebuilt +(verified: the skip is logged and the file's mtime does not move). + +Resuming a **completed** run does not, and did not +before this branch either: `clean_up()` deletes `out_raw_file + "_*"`, which includes the +survivors files, and `prepare_read_filter` opens them without an existence guard. Verified on +the branch base (b59bfbcc), where the same scenario fails even earlier. `clean_up` and that +guard are untouched here. + +The mechanism, for whoever fixes it: `clean_up` removes the *global* UMI lock but not the +per-edit-distance one (`umi_filtered_lock_file_name`, in `aux/`), so a resumed `filter_umis` +returns early without rewriting the survivors files it just deleted. `write_deduplicated_bam` +then finds nothing and warns "No reads survived UMI filtering", which is a misdiagnosis — but +the run dies seconds later in `prepare_read_filter` for the same underlying reason, so nothing +was added here to paper over it. diff --git a/docs/barcode_calling.md b/docs/barcode_calling.md index a80472fe..56a0bdee 100644 --- a/docs/barcode_calling.md +++ b/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 (`.split_reads.fasta.gz`, numbered +`_.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 +`.split_reads_.fa.gz`. **Curio** (`curio`): diff --git a/docs/cmd.md b/docs/cmd.md index b401e9ca..8c1b995f 100644 --- a/docs/cmd.md +++ b/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. diff --git a/docs/formats.md b/docs/formats.md index 7df8243f..ea12e03b 100644 --- a/docs/formats.md +++ b/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 diff --git a/docs/output.md b/docs/output.md index 53867396..1174d46b 100644 --- a/docs/output.md +++ b/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. diff --git a/docs/single_cell.md b/docs/single_cell.md index def2711f..295eed1d 100644 --- a/docs/single_cell.md +++ b/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_.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 diff --git a/isoquant.py b/isoquant.py index 1e82cefc..d14fb0d1 100755 --- a/isoquant.py +++ b/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) diff --git a/isoquant_detect_barcodes.py b/isoquant_detect_barcodes.py index bb3b30fd..712d9386 100755 --- a/isoquant_detect_barcodes.py +++ b/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): diff --git a/isoquant_lib/assignment/assignment_io.py b/isoquant_lib/assignment/assignment_io.py index 4361631a..caf4d6f2 100644 --- a/isoquant_lib/assignment/assignment_io.py +++ b/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() diff --git a/isoquant_lib/assignment/assignment_loader.py b/isoquant_lib/assignment/assignment_loader.py index 5e8dc115..516b9b30 100644 --- a/isoquant_lib/assignment/assignment_loader.py +++ b/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 diff --git a/isoquant_lib/assignment/read_groups.py b/isoquant_lib/assignment/read_groups.py index 16aae8fb..5ae94b15 100644 --- a/isoquant_lib/assignment/read_groups.py +++ b/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 diff --git a/isoquant_lib/barcode_calling/detect_barcodes.py b/isoquant_lib/barcode_calling/detect_barcodes.py index bdee20e6..4779056e 100644 --- a/isoquant_lib/barcode_calling/detect_barcodes.py +++ b/isoquant_lib/barcode_calling/detect_barcodes.py @@ -28,6 +28,8 @@ from ..modes import IsoQuantMode from isoquant_lib.utils.error_codes import IsoQuantExitCode from ..common import setup_worker_logging, _get_log_params +from isoquant_lib.utils.file_utils import (GZIP_SUFFIX, open_text_write, + strip_compression_suffix) from .common import reverese_complement, load_barcodes from .cell_selection import NOSEQ, CellBarcodeSelector, select_cell_barcodes from . import ( @@ -156,7 +158,7 @@ def __init__(self, output_file_name, barcode_detector, header=False, output_sequ self.output_sequences_file = None self.process_function = self._process_read_split if split_reads else self._process_read_normal if self.output_sequences: - self.output_sequences_file = open(self.output_sequences, "w") + self.output_sequences_file = open_text_write(self.output_sequences) if header: self.output_file.write(barcode_detector.header() + "\n") self.read_stat = ReadStats() @@ -336,10 +338,22 @@ def setup_detector_worker(log_file, log_level, barcode_detector): _WORKER_DETECTOR["detector"] = barcode_detector +def numbered_chunk_name(file_name, num): + """Append a chunk index, keeping the extensions last. + + The suffixes have to survive: a chunk of a FASTA is still a FASTA, and that is what decides + its compression level. + """ + base = strip_compression_suffix(file_name) + compression = file_name[len(base):] + root, extension = os.path.splitext(base) + return "%s_%d%s%s" % (root, num, extension, compression) + + def process_chunk(read_chunk, output_file, num, out_fasta=None, split_reads=False, barcode_detector=None): output_file += "_" + str(num) if out_fasta: - out_fasta += "_" + str(num) + out_fasta = numbered_chunk_name(out_fasta, num) counter = 0 if barcode_detector is None: @@ -391,8 +405,10 @@ def count_barcodes_in_reads(args, barcode_length): def submit(pool, chunk, num): return pool.submit(count_chunk, chunk, split_reads, barcode_length) - def handle_result(result): + def handle_result(result, chunk_index): counts, malformed, reads = result + # order-insensitive: the counts are summed and CellBarcodeSelector sorts them by + # (count, barcode), so the chunk index is of no use here selector.merge_counts(counts, malformed) return reads @@ -509,9 +525,13 @@ def open_read_chunks(input_file): def run_chunks_in_parallel(read_chunk_gen, args, barcode_detector, submit, handle_result): """Feed read chunks to a worker pool, keeping args.threads tasks in flight. + Results are handed to handle_result in completion order, so a consumer whose output depends + on the order gets the chunk index and has to reorder by it -- see _process_single_file_in_parallel. + Waiting for the chunks in order instead would idle the pool behind a single slow one. + Args: submit: (pool, chunk, chunk_index) -> Future - handle_result: (result) -> number of reads processed + handle_result: (result, chunk_index) -> number of reads processed """ # Clean up parent memory before spawning workers gc.collect() @@ -523,10 +543,10 @@ def run_chunks_in_parallel(read_chunk_gen, args, barcode_detector, submit, handl mp_context=mp_context, initializer=setup_detector_worker, initargs=(log_file, log_level, barcode_detector)) as proc: - future_results = [] + future_results = {} # future -> the index of the chunk it is processing chunk_counter = 0 for chunk in read_chunk_gen: - future_results.append(submit(proc, chunk, chunk_counter)) + future_results[submit(proc, chunk, chunk_counter)] = chunk_counter chunk_counter += 1 if chunk_counter >= args.threads: break @@ -539,12 +559,11 @@ def run_chunks_in_parallel(read_chunk_gen, args, barcode_detector, submit, handl for c in completed_features: if c.exception() is not None: raise c.exception() - read_counter += handle_result(c.result()) + read_counter += handle_result(c.result(), future_results.pop(c)) sys.stdout.write("Processed %d reads\r" % read_counter) - future_results.remove(c) if reads_left: try: - future_results.append(submit(proc, next(read_chunk_gen), chunk_counter)) + future_results[submit(proc, next(read_chunk_gen), chunk_counter)] = chunk_counter chunk_counter += 1 except StopIteration: reads_left = False @@ -566,28 +585,41 @@ def _process_single_file_in_parallel(input_file, output_tsv, out_fasta, args, ba os.makedirs(tmp_dir) tmp_barcode_file = os.path.join(tmp_dir, "bc") - tmp_fasta_file = os.path.join(tmp_dir, "subreads") if out_fasta else None - output_files = [] + tmp_fasta_file = None + if out_fasta: + # compress the per-chunk temps exactly when the final file is compressed: the work + # then happens in the workers, and merging stays a byte concat + # named for what it holds, so the chunk inherits the FASTA compression level + tmp_fasta_file = os.path.join(tmp_dir, "subreads.fa") + if out_fasta.endswith(GZIP_SUFFIX): + tmp_fasta_file += GZIP_SUFFIX + chunk_outputs = {} def submit(pool, chunk, num): return pool.submit(process_chunk, chunk, tmp_barcode_file, num, tmp_fasta_file, split_reads) - def handle_result(result): + def handle_result(result, chunk_index): tmp_out_file, tmp_out_fasta, read_count = result - output_files.append((tmp_out_file, tmp_out_fasta)) + chunk_outputs[chunk_index] = (tmp_out_file, tmp_out_fasta) return read_count run_chunks_in_parallel(read_chunk_gen, args, barcode_detector, submit, handle_result) + # by chunk index, not by completion: the merge below defines the row order of the barcode + # table and the record order of the FASTA, and both have to be the same on every run + output_files = [chunk_outputs[i] for i in sorted(chunk_outputs)] with open(output_tsv, "w") as final_output_tsv: - final_output_fasta = open(out_fasta, "w") if out_fasta else None + # binary: concatenating gzip members byte-wise gives a valid multi-member stream, + # so per-chunk compression in the workers merges without re-compressing here + final_output_fasta = open(out_fasta, "wb") if out_fasta else None header = barcode_detector.header() final_output_tsv.write(header + "\n") stat_dict = defaultdict(int) for tmp_file, tmp_fasta in output_files: shutil.copyfileobj(open(tmp_file, "r"), final_output_tsv) if tmp_fasta and final_output_fasta: - shutil.copyfileobj(open(tmp_fasta, "r"), final_output_fasta) + with open(tmp_fasta, "rb") as tmp_fasta_handle: + shutil.copyfileobj(tmp_fasta_handle, final_output_fasta) for line in open(stats_file_name(tmp_file), "r"): v = line.strip().split("\t") if len(v) != 2: diff --git a/isoquant_lib/barcode_calling/pipeline.py b/isoquant_lib/barcode_calling/pipeline.py index 121067c8..2af84118 100644 --- a/isoquant_lib/barcode_calling/pipeline.py +++ b/isoquant_lib/barcode_calling/pipeline.py @@ -118,7 +118,10 @@ def call_barcodes(args): output_fasta_list = None new_reads = [] if args.split_molecules: - output_fasta_list = [sample.split_reads_fasta + "_%d.fa" % i for i in range(len(input_files))] + # minimap2 reads gzipped FASTA natively, so compressing costs nothing downstream + fasta_suffix = ".fa.gz" if args.gzipped else ".fa" + output_fasta_list = [sample.split_reads_fasta + "_%d%s" % (i, fasta_suffix) + for i in range(len(input_files))] new_reads = [[fasta] for fasta in output_fasta_list] # Check if all files were already processed during resume diff --git a/isoquant_lib/barcode_calling/umi_filtering.py b/isoquant_lib/barcode_calling/umi_filtering.py index 433823ef..7c78960a 100644 --- a/isoquant_lib/barcode_calling/umi_filtering.py +++ b/isoquant_lib/barcode_calling/umi_filtering.py @@ -144,7 +144,7 @@ class UMIFilter: def __init__(self, umi_length: int = 0, edit_distance: int = 3, disregard_length_diff: bool = True, only_unique_assignments: bool = False, only_spliced_reads: bool = False, - barcode_remap: Optional[Dict[str, str]] = None): + barcode_remap: Optional[Dict[str, str]] = None, output_read_tags: bool = False): """ Initialize UMI filter. @@ -155,6 +155,8 @@ def __init__(self, umi_length: int = 0, edit_distance: int = 3, disregard_length only_unique_assignments: Only process uniquely assigned reads only_spliced_reads: Only process spliced reads barcode_remap: Optional mapping from barcode to spot ID for spot-level dedup + output_read_tags: also record barcode/UMI/gene/transcript for each surviving read, + so a deduplicated BAM can be tagged without re-reading the assignments """ self.umi_length = umi_length self.max_edit_distance = edit_distance @@ -162,6 +164,7 @@ def __init__(self, umi_length: int = 0, edit_distance: int = 3, disregard_length self.only_unique_assignments = only_unique_assignments self.only_spliced_reads = only_spliced_reads self.barcode_remap = barcode_remap + self.output_read_tags = output_read_tags self.selected_reads: Set[str] = set() self.stats: Dict[str, int] = defaultdict(int) @@ -368,6 +371,24 @@ def _process_gene(self, gene_dict: Dict[str, List[ReadAssignment]]): for r in self._process_duplicates(gene_dict[barcode]): yield r + def _survivor_record(self, read_assignment) -> str: + """One line of the survivors file: the read id, plus tags when they are wanted. + + Everything comes from the assignment already in hand, so no second pass is needed. + """ + if not self.output_read_tags: + return read_assignment.read_id + if read_assignment.isoform_matches: + gene_id = read_assignment.isoform_matches[0].assigned_gene or "." + transcript_id = read_assignment.isoform_matches[0].assigned_transcript or "." + else: + gene_id, transcript_id = ".", "." + return "\t".join([read_assignment.read_id, + read_assignment.barcode or "*", + read_assignment.umi or "*", + gene_id, + transcript_id]) + def _process_chunk(self, gene_barcode_dict: Dict[str, Dict[str, List[ReadAssignment]]], allinfo_outf, read_ids_outf=None) -> Tuple[int, int]: """ @@ -400,7 +421,7 @@ def _process_chunk(self, gene_barcode_dict: Dict[str, Dict[str, List[ReadAssignm spliced_count += 1 if read_ids_outf: - read_ids_outf.write(read_assignment.read_id + "\n") + read_ids_outf.write(self._survivor_record(read_assignment) + "\n") self.selected_reads.add(read_assignment.read_id) diff --git a/isoquant_lib/dataset_processor.py b/isoquant_lib/dataset_processor.py index de868833..88c971aa 100644 --- a/isoquant_lib/dataset_processor.py +++ b/isoquant_lib/dataset_processor.py @@ -7,7 +7,6 @@ import gc import glob -import gzip import itertools import logging import multiprocessing @@ -34,7 +33,11 @@ write_string, ) from isoquant_lib.utils.stats import EnumStats -from isoquant_lib.utils.file_utils import merge_files, merge_counts +from isoquant_lib.utils.file_utils import (merge_files, merge_counts, gzip_file_in_place, + open_text_write, resolve_optionally_gzipped) +from isoquant_lib.utils.bam_utils import (PLACEHOLDERS, collect_unmapped_read_ids, + load_barcode_umi_tags, merge_bam_files, + references_with_alignments, write_unmapped_bam) from .alignment.alignment_processor import AlignmentType from .assignment.read_groups import prepare_read_groups, get_grouping_strategy_names from .assignment.assignment_io import IOSupport, ReadInfoPrinter, VoidPrinter @@ -52,6 +55,8 @@ umi_filtered_global_lock_file_name, umi_filtered_lock_file_name, umi_output_prefix, + tagged_bam_fragment_name, + tagged_bam_lock_filename, ) from isoquant_lib.model_construction.transcript_printer import GFFPrinter, VoidTranscriptPrinter from .barcode_calling.umi_filtering import create_transcript_info_dict @@ -62,6 +67,8 @@ collect_reads_in_parallel, construct_models_in_parallel, filter_umis_in_parallel, + write_deduplicated_bam_in_parallel, + write_tagged_bam_in_parallel, ) logger = logging.getLogger('IsoQuant') @@ -134,6 +141,8 @@ def process_all_samples(self, input_data): logger.info("Secondary alignments will%s be used" % ("" if self.args.use_secondary else " not")) for sample in input_data.samples: self.process_sample(sample) + for sample in input_data.samples: + self.compress_barcode_tables(sample) self.clean_up() logger.info("Processed " + proper_plural_form("experiment", len(self.input_data.samples))) @@ -177,8 +186,15 @@ def process_sample(self, sample): if self.args.barcoded_reads: sample.barcoded_reads = self.args.barcoded_reads - for chr_id in self.get_chr_list(): - split_barcodes_dict[chr_id] = sample.barcodes_split_reads + "_" + chr_id + # a tagged BAM copies every reference, including the unplaced scaffolds IsoQuant + # does not analyse, so those reads need a barcode table of their own to be tagged. + # The same list drives the split and the copy, so every fragment finds its table. + tagged_bam_references = None + if large_output_enabled(self.args, "tagged_bam"): + tagged_bam_references = references_with_alignments([f[0] for f in sample.file_list]) + split_chr_ids = self.get_chr_list() if tagged_bam_references is None else tagged_bam_references + for chr_id in split_chr_ids: + split_barcodes_dict[chr_id] = sample.get_barcodes_split_file(chr_id) barcode_split_done = split_barcodes_lock_filename(sample) if self.args.resume and os.path.exists(barcode_split_done): logger.info("Barcode table was split during the previous run, existing files will be used") @@ -188,6 +204,19 @@ def process_sample(self, sample): self.split_read_barcode_table(sample, split_barcodes_dict) open(barcode_split_done, "w").close() + # nothing to tag with otherwise; the user was warned about that at startup + if tagged_bam_references is not None: + tagged_bam_done = tagged_bam_lock_filename(sample) + if (self.args.resume and os.path.exists(tagged_bam_done) + and os.path.exists(sample.out_tagged_bam)): + logger.info("Tagged BAM was written during the previous run, keeping %s" + % sample.out_tagged_bam) + else: + if os.path.exists(tagged_bam_done): + os.remove(tagged_bam_done) + self.write_tagged_bam(sample, tagged_bam_references) + open(tagged_bam_done, "w").close() + if self.args.read_assignments: saves_file = self.args.read_assignments[0] logger.info('Using read assignments from {}*'.format(saves_file)) @@ -201,6 +230,9 @@ def process_sample(self, sample): if self.args.mode.needs_pcr_deduplication(): self.filter_umis(sample) + # the survivors files live under out_raw_file and are deleted by clean_up + if large_output_enabled(self.args, "deduplicated_bam"): + self.write_deduplicated_bam(sample) total_assignments, polya_found, self.all_read_groups = self.load_read_info(saves_file) @@ -603,9 +635,7 @@ def filter_umis(self, sample): allinfo_fname = output_prefix + ".allinfo" if self.args.gzipped: allinfo_fname += ".gz" - allinfo_outf = gzip.open(allinfo_fname, "wt") - else: - allinfo_outf = open(allinfo_fname, "w") + allinfo_outf = open_text_write(allinfo_fname) for all_info_file_name, stats_output_file_name, umi_filter_done in results: if save_allinfo: @@ -693,9 +723,7 @@ def filter_umis(self, sample): allinfo_fname = output_prefix + ".allinfo" if self.args.gzipped: allinfo_fname += ".gz" - allinfo_outf = gzip.open(allinfo_fname, "wt") - else: - allinfo_outf = open(allinfo_fname, "w") + allinfo_outf = open_text_write(allinfo_fname) for all_info_file_name, stats_output_file_name, umi_filter_done in results: if save_allinfo: @@ -727,11 +755,94 @@ def filter_umis(self, sample): open(umi_filtering_done, "w").close() + def map_over_chromosomes(self, worker, sample, *extra_args, chr_ids=None): + """Run worker(sample, chr_id, *extra_args) for every chromosome, in parallel. + + Defaults to the chromosomes this run analysed; pass chr_ids to cover others. + """ + gen = (worker, itertools.repeat(sample), + self.get_chr_list() if chr_ids is None else chr_ids, + *(itertools.repeat(a) for a in extra_args)) + if self.args.threads > 1: + gc.collect() + mp_context = multiprocessing.get_context('fork') + log_file, log_level = _get_log_params() + with ProcessPoolExecutor(max_workers=self.args.threads, mp_context=mp_context, + initializer=setup_worker_logging, + initargs=(log_file, log_level)) as proc: + return list(proc.map(*gen, chunksize=1)) + return list(map(*gen)) + + def write_tagged_bam(self, sample, references): + """A copy of the input BAM(s) with barcode and UMI tags, keeping every alignment. + + Purely a side output: the barcode table split it reads from happens regardless, and + nothing downstream looks at the result. references is the list the barcode table was + split over, so each fragment has a table to read its tags from. + """ + logger.info("Writing tagged BAM") + bam_files = [f[0] for f in sample.file_list] + fragments = self.map_over_chromosomes(write_tagged_bam_in_parallel, sample, self.args, + chr_ids=references) + + # fetch() by chromosome never returns unmapped reads, so they need a pass of their own. + # They belong to no chromosome and so appear in no split table, but a barcode is called + # from the read sequence and does not need an alignment -- look theirs up directly. + unmapped_fragment = tagged_bam_fragment_name(sample.out_raw_file, "unmapped") + unmapped_ids = collect_unmapped_read_ids(bam_files) + unmapped_tags = {} + for table in sample.barcoded_reads: + unmapped_tags.update(load_barcode_umi_tags(resolve_optionally_gzipped(table), + read_ids=unmapped_ids)) + if write_unmapped_bam(bam_files, unmapped_fragment, unmapped_tags, + self.args.barcode_tag, self.args.umi_tag): + # a table row exists for every read, but its barcode may be the uncalled placeholder + barcoded = sum(1 for tags in unmapped_tags.values() if tags[0] not in PLACEHOLDERS) + logger.info("Copied %d unmapped reads, %d of them barcoded" + % (len(unmapped_ids), barcoded)) + fragments.append(unmapped_fragment) + else: + os.remove(unmapped_fragment) + + merged = merge_bam_files(sample.out_tagged_bam, fragments, self.args.threads) + if merged is None: + logger.warning("No alignments to tag, %s was not written" % sample.out_tagged_bam) + else: + logger.info("Tagged BAM saved to %s" % merged) + return merged + + def write_deduplicated_bam(self, sample): + """A primary-only BAM of the reads that survived UMI filtering, carrying their tags.""" + logger.info("Writing deduplicated BAM") + fragments = self.map_over_chromosomes(write_deduplicated_bam_in_parallel, sample, self.args) + merged = merge_bam_files(sample.out_deduplicated_bam, fragments, self.args.threads) + if merged is None: + logger.warning("No reads survived UMI filtering, %s was not written" + % sample.out_deduplicated_bam) + else: + logger.info("Deduplicated BAM saved to %s" % merged) + return merged + + def compress_barcode_tables(self, sample): + """Gzip the barcoded read tables once nothing needs them any more. + + They stay plain for the whole run because split_read_barcode_table reads them, and + only tables this run produced are touched -- files passed via --barcoded_reads belong + to the user and are left alone. + """ + if not self.args.gzipped: + return + produced = sorted(glob.glob(sample.barcodes_tsv + "_*.tsv")) + for table in produced: + gzipped = gzip_file_in_place(table) + logger.info("Compressed barcode table to %s" % gzipped) + def split_read_barcode_table(self, sample, split_barcodes_file_names): logger.info("Splitting read barcode table") # Supports both IsoQuant's 6-column format and third-party 3-column format # Only columns 0, 1, 2 (read_id, barcode, umi) are required and preserved - split_read_table_parallel(sample, sample.barcoded_reads, split_barcodes_file_names, self.args.threads, + barcode_tables = [resolve_optionally_gzipped(f) for f in sample.barcoded_reads] + split_read_table_parallel(sample, barcode_tables, split_barcodes_file_names, self.args.threads, read_column=0, group_columns=(1, 2), delim='\t') logger.info("Read barcode table was split") diff --git a/isoquant_lib/parallel_workers.py b/isoquant_lib/parallel_workers.py index 694aeb0b..086f1e77 100644 --- a/isoquant_lib/parallel_workers.py +++ b/isoquant_lib/parallel_workers.py @@ -27,6 +27,8 @@ read_stat_file_name, transcript_stat_file_name, umi_filtered_lock_file_name, + dedup_bam_fragment_name, + tagged_bam_fragment_name, allinfo_file_name, allinfo_stats_file_name, ) @@ -38,6 +40,8 @@ from .assignment.assignment_aggregator import ReadAssignmentAggregator from isoquant_lib.utils.string_pools import setup_string_pools from .common import large_output_enabled +from isoquant_lib.utils.bam_utils import (load_barcode_umi_tags, load_survivor_tags, + write_tagged_chromosome_bam) logger = logging.getLogger('IsoQuant') @@ -392,7 +396,11 @@ def filter_umis_in_parallel(sample, chr_id, chr_ids, args, edit_distance, output # When using barcode_remap (spot-based dedup), never produce filtered reads if barcode_remap: output_filtered_reads = False - umi_filter = UMIFilter(args.umi_length, edit_distance, barcode_remap=barcode_remap) + # a deduplicated BAM needs the tags of every surviving read, and this is the only place + # where the assignments carrying them are already loaded + output_read_tags = output_filtered_reads and large_output_enabled(args, "deduplicated_bam") + umi_filter = UMIFilter(args.umi_length, edit_distance, barcode_remap=barcode_remap, + output_read_tags=output_read_tags) filtered_reads = sample.get_filtered_reads_file(chr_id) if output_filtered_reads else None umi_filter.process_single_chr(chr_id, sample.out_raw_file, transcript_type_dict, @@ -405,3 +413,36 @@ def filter_umis_in_parallel(sample, chr_id, chr_ids, args, edit_distance, output logger.info("PCR duplicates filtered for chromosome " + chr_id) return all_info_file_name, stats_output_file_name, umi_filtered_done + + +def write_deduplicated_bam_in_parallel(sample, chr_id, args): + """Write one chromosome of the deduplicated BAM: primary alignments of surviving reads. + + Returns the fragment path, or None when nothing survived on this chromosome. + """ + read_tags = load_survivor_tags(sample.get_filtered_reads_file(chr_id)) + if not read_tags: + return None + + bam_files = [f[0] for f in sample.file_list] + fragment = dedup_bam_fragment_name(sample.out_raw_file, chr_id) + written = write_tagged_chromosome_bam(bam_files, fragment, chr_id, read_tags, + barcode_tag=args.barcode_tag, umi_tag=args.umi_tag, + primary_only=True, keep_untagged=False) + logger.info("Deduplicated %d alignments for chromosome %s" % (written, chr_id)) + return fragment + + +def write_tagged_bam_in_parallel(sample, chr_id, args): + """Write one chromosome of the tagged BAM: every alignment, barcoded ones carrying tags. + + A chromosome whose reads were all unbarcoded has no split table; its alignments are still + copied, just without tags, so no alignment is lost for want of a barcode. + """ + read_tags = load_barcode_umi_tags(sample.get_barcodes_split_file(chr_id)) + bam_files = [f[0] for f in sample.file_list] + fragment = tagged_bam_fragment_name(sample.out_raw_file, chr_id) + written = write_tagged_chromosome_bam(bam_files, fragment, chr_id, read_tags, + barcode_tag=args.barcode_tag, umi_tag=args.umi_tag) + logger.info("Tagged %d alignments for chromosome %s" % (written, chr_id)) + return fragment diff --git a/isoquant_lib/quantification/convert_grouped_counts.py b/isoquant_lib/quantification/convert_grouped_counts.py index e423934b..47855294 100644 --- a/isoquant_lib/quantification/convert_grouped_counts.py +++ b/isoquant_lib/quantification/convert_grouped_counts.py @@ -13,9 +13,10 @@ from traceback import print_exc import pandas import logging -import gzip import gffutils +from isoquant_lib.utils.file_utils import open_text_write + try: from isoquant_lib.utils.error_codes import IsoQuantExitCode from isoquant_lib.utils.file_naming import ( @@ -87,7 +88,7 @@ def convert_to_matrix(input_linear_counts, output_file_path, feature_id_to_name_ output_file_path += ".tsv" output_file_path += ".gz" if gzipped else "" - with gzip.open(output_file_path, 'wt') if gzipped else open(output_file_path, 'w') as outfile: + with open_text_write(output_file_path) as outfile: # Write the header with group_ids columns = list(sorted(count_matrix.columns)) if num_groups > GROUP_COUNT_CUTOFF: @@ -137,7 +138,7 @@ def convert_to_mtx(input_linear_counts, output_file_prefix, feature_id_to_name=N gene_name = feature_id_to_name.get(gene_id, gene_id) if feature_id_to_name is not None else gene_id ft_out.write(f"{gene_id}\t{gene_name}\n") - with gzip.open(mtx_file + ".gz", 'wt') if gzipped else open(mtx_file, 'w') as mtx_out: + with open_text_write(mtx_file + ".gz" if gzipped else mtx_file) as mtx_out: # Write the header mtx_out.write("%%MatrixMarket matrix coordinate real general\n") mtx_out.write(f"{len(unique_genes)} {len(unique_groups)} {df.shape[0]}\n") @@ -231,7 +232,7 @@ def convert_profile_to_matrix(input_linear_counts: str, output_file_path: str, logger.warning("You have %d groups in your matrix, conversion might take a lot of time " "and the output file can be very large" % num_groups) - with gzip.open(output_file, 'wt') if gzipped else open(output_file, 'w') as outfile: + with open_text_write(output_file) as outfile: outfile.write(feature_type + '_id\t' + '\t'.join(columns) + '\n') for feature_id in feature_order: incl_row = incl.loc[feature_id] @@ -276,7 +277,7 @@ def convert_profile_to_mtx(input_linear_counts: str, output_file_prefix: str, excl_nonzero = df[df['exclude_counts'] != 0] def _write_mtx(path, frame, value_col): - with gzip.open(path + ".gz", 'wt') if gzipped else open(path, 'w') as out: + with open_text_write(path + ".gz" if gzipped else path) as out: out.write("%%MatrixMarket matrix coordinate real general\n") out.write(f"{len(feature_order)} {len(unique_groups)} {frame.shape[0]}\n") for _, row in frame.iterrows(): diff --git a/isoquant_lib/scripts/convert_read_info.py b/isoquant_lib/scripts/convert_read_info.py index 5d098040..35f81efb 100644 --- a/isoquant_lib/scripts/convert_read_info.py +++ b/isoquant_lib/scripts/convert_read_info.py @@ -25,6 +25,7 @@ from typing import TextIO from isoquant_lib.common import junctions_from_blocks +from isoquant_lib.utils.file_utils import open_text_write # read_info column indices @@ -49,6 +50,8 @@ def _open_file(path: str, mode: str = "r") -> TextIO: + if mode.startswith("w"): + return open_text_write(path) if path.endswith(".gz"): return gzip.open(path, mode + "t") return open(path, mode) diff --git a/isoquant_lib/scripts/exon_splice_site_to_group_lists.py b/isoquant_lib/scripts/exon_splice_site_to_group_lists.py index 732d961f..dcec8ad0 100644 --- a/isoquant_lib/scripts/exon_splice_site_to_group_lists.py +++ b/isoquant_lib/scripts/exon_splice_site_to_group_lists.py @@ -32,8 +32,12 @@ import gzip from typing import TextIO +from isoquant_lib.utils.file_utils import open_text_write + def _open(path: str, mode: str = "r") -> TextIO: + if mode.startswith("w"): + return open_text_write(path) if path.endswith(".gz"): return gzip.open(path, mode + "t") return open(path, mode) diff --git a/isoquant_lib/utils/bam_utils.py b/isoquant_lib/utils/bam_utils.py new file mode 100644 index 00000000..943ac272 --- /dev/null +++ b/isoquant_lib/utils/bam_utils.py @@ -0,0 +1,249 @@ +############################################################################ +# Copyright (c) 2025-2026 University of Helsinki +# # All Rights Reserved +# See file LICENSE for details. +############################################################################ + +""" +Writing BAM outputs: merging per-chromosome fragments, and copying reads with tags. + +The pipeline processes chromosomes in parallel, so any BAM it emits is built as one fragment +per chromosome and merged at the end. Records are copied verbatim apart from the tags being +added, because downstream consumers depend on fields IsoQuant does not itself use -- fusion +detection in particular reads breakpoints from the SA tag and from terminal soft clips. +""" + +import logging +import os +import shutil +from typing import Dict, Iterable, Iterator, List, Optional, Set, Tuple + +import pysam + +from isoquant_lib.utils.file_utils import open_text_read + +logger = logging.getLogger('IsoQuant') + +# Tags for the gene and transcript a read was assigned to, following the 10x convention. +GENE_TAG = "GX" +TRANSCRIPT_TAG = "TX" + +# read_id -> (barcode, umi, gene, transcript); any element may be None +ReadTags = Dict[str, Tuple[Optional[str], Optional[str], Optional[str], Optional[str]]] + +# Stand-ins for "nothing here" in the tables these tags are read from: "*" for a barcode or +# UMI that was not called, "." for a feature this code left unset, and the literal "None" that +# IsoQuant writes for a read with no assigned transcript (novel and ambiguous reads). A tag is +# omitted rather than carrying any of them. +PLACEHOLDERS = frozenset(("*", ".", "None")) + +# How many fragments one samtools merge call may take; see _merge_in_rounds. +BAM_MERGE_BATCH = 500 + + +def index_bam(bam_file: str, threads: int = 1) -> None: + """Index a BAM, falling back to CSI for references too long for BAI.""" + try: + pysam.index('-@', str(threads), bam_file) + except pysam.SamtoolsError as err: + logger.info("Samtools failed to generate a .bai index: %s" % err) + logger.info("Trying to build a CSI index instead") + try: + pysam.index('-@', str(threads), '-c', bam_file) + except pysam.SamtoolsError as csi_err: + logger.error("Failed to create CSI index: %s" % csi_err) + + +def _merge_in_rounds(output_bam: str, fragments: List[str], threads: int) -> List[str]: + """Merge fragments into output_bam, in rounds when there are too many for one call. + + samtools opens every input at once and gives up a little above a thousand handles, which a + reference with many unplaced scaffolds reaches easily. Returns the intermediate files the + rounds produced, for the caller to delete. + """ + intermediates = [] + current = list(fragments) + round_index = 0 + while len(current) > BAM_MERGE_BATCH: + merged = [] + for start in range(0, len(current), BAM_MERGE_BATCH): + batch = current[start:start + BAM_MERGE_BATCH] + if len(batch) == 1: + # a lone leftover: carry it into the next round rather than copying it + merged.append(batch[0]) + continue + part = "%s.merge%d_%d.bam" % (output_bam, round_index, start // BAM_MERGE_BATCH) + pysam.merge("-f", "-@", str(threads), part, *batch) + intermediates.append(part) + merged.append(part) + current = merged + round_index += 1 + pysam.merge("-f", "-@", str(threads), output_bam, *current) + return intermediates + + +def merge_bam_files(output_bam: str, input_bams: Iterable[str], threads: int = 1, + remove_fragments: bool = True) -> Optional[str]: + """Merge per-chromosome BAM fragments into one indexed BAM. + + The fragments share their header (each is written from the same template) and cover + disjoint references, so a plain coordinate merge is correct. + """ + fragments = [bam for bam in input_bams if bam and os.path.exists(bam)] + if not fragments: + logger.warning("No BAM fragments to merge into %s" % output_bam) + return None + + if len(fragments) == 1: + shutil.move(fragments[0], output_bam) + else: + for intermediate in _merge_in_rounds(output_bam, fragments, threads): + os.remove(intermediate) + if remove_fragments: + for fragment in fragments: + os.remove(fragment) + + index_bam(output_bam, threads) + logger.info("Merged %d BAM fragments into %s" % (len(fragments), output_bam)) + return output_bam + + +def references_with_alignments(input_bams: Iterable[str]) -> List[str]: + """Every reference that actually carries alignments, across all input BAMs. + + A tagged BAM is meant to be a copy of the input, so it has to cover references IsoQuant + itself skipped -- unplaced scaffolds hold alignments too. Empty references are dropped so + a fragmented assembly does not spawn thousands of tasks that copy nothing. + """ + references = set() + for bam_file in input_bams: + with pysam.AlignmentFile(bam_file, "rb") as inf: + try: + references.update(stat.contig for stat in inf.get_index_statistics() + if stat.mapped + stat.unmapped > 0) + except ValueError: + # no index: fall back to the header, empty references cost an empty fragment + references.update(inf.references) + return sorted(references) + + +def load_survivor_tags(survivors_file: str) -> ReadTags: + """Read a UMI-filtering survivors table into read_id -> (barcode, umi, gene, transcript). + + The table carries the tag columns only when a tagged BAM was requested; without them the + read ids are still the definition of the surviving subset, so they map to empty tags. + """ + read_tags = {} + if not os.path.exists(survivors_file): + return read_tags + with open(survivors_file) as handle: + for line in handle: + fields = line.rstrip("\n").split("\t") + if len(fields) < 5: + read_tags[fields[0]] = (None, None, None, None) + else: + read_tags[fields[0]] = (fields[1], fields[2], fields[3], fields[4]) + return read_tags + + +def load_barcode_umi_tags(barcode_table: str, read_ids: Optional[Set[str]] = None) -> ReadTags: + """Read a barcode table into read_id -> (barcode, umi, None, None). + + read_ids restricts the result to those reads, so the whole-sample table can be scanned for + a handful of reads without holding millions of entries in memory. + """ + read_tags = {} + if not os.path.exists(barcode_table): + return read_tags + with open_text_read(barcode_table) as handle: + for line in handle: + if line.startswith("#"): + continue + fields = line.rstrip("\n").split("\t") + if len(fields) < 3: + continue + if read_ids is None or fields[0] in read_ids: + read_tags[fields[0]] = (fields[1], fields[2], None, None) + return read_tags + + +def unplaced_reads(inf: pysam.AlignmentFile) -> Iterator[pysam.AlignedSegment]: + """The unmapped records that carry no coordinates, which fetch() by chromosome never yields. + + They sit at the end of a coordinate-sorted file, so an index turns reaching them into a + seek; without one the whole file has to be read. + """ + if inf.has_index(): + return inf.fetch("*") + return (read for read in inf.fetch(until_eof=True) + if read.is_unmapped and read.reference_id == -1) + + +def collect_unmapped_read_ids(input_bams: Iterable[str]) -> Set[str]: + """Ids of the unplaced reads, so their tags can be looked up before they are copied.""" + read_ids = set() + for bam_file in input_bams: + with pysam.AlignmentFile(bam_file, "rb") as inf: + for read in unplaced_reads(inf): + read_ids.add(read.query_name) + return read_ids + + +def _apply_tags(read, tags, barcode_tag: str, umi_tag: str) -> None: + for tag_name, value in zip((barcode_tag, umi_tag, GENE_TAG, TRANSCRIPT_TAG), tags): + if value and value not in PLACEHOLDERS: + read.set_tag(tag_name, value) + + +def write_tagged_chromosome_bam(input_bams: Iterable[str], output_bam: str, chr_id: str, + read_tags: ReadTags, barcode_tag: str = "CB", umi_tag: str = "UB", + primary_only: bool = False, keep_untagged: bool = True) -> int: + """Copy one chromosome's alignments into output_bam, adding tags from read_tags. + + primary_only drops secondary and supplementary records. keep_untagged decides whether a + read absent from read_tags is copied without tags (a tagged copy of the input) or skipped + (a subset restricted to the listed reads). + """ + bam_list = list(input_bams) + written = 0 + with pysam.AlignmentFile(bam_list[0], "rb") as template: + with pysam.AlignmentFile(output_bam, "wb", template=template) as outf: + for bam_file in bam_list: + with pysam.AlignmentFile(bam_file, "rb") as inf: + for read in inf.fetch(chr_id): + if primary_only and (read.is_secondary or read.is_supplementary): + continue + tags = read_tags.get(read.query_name) + if tags is None and not keep_untagged: + continue + if tags is not None: + _apply_tags(read, tags, barcode_tag, umi_tag) + outf.write(read) + written += 1 + return written + + +def write_unmapped_bam(input_bams: Iterable[str], output_bam: str, + read_tags: Optional[ReadTags] = None, + barcode_tag: str = "CB", umi_tag: str = "UB") -> int: + """Copy unmapped reads, which fetch() by chromosome never returns. + + Only truly unplaced records (reference_id == -1) are taken: an unmapped read parked at a + reference position does come back from fetch(), and copying it here too would duplicate it. + + These reads get tagged like any other: a barcode is called from the read sequence, so it + does not depend on the read having aligned anywhere. + """ + bam_list = list(input_bams) + written = 0 + with pysam.AlignmentFile(bam_list[0], "rb") as template: + with pysam.AlignmentFile(output_bam, "wb", template=template) as outf: + for bam_file in bam_list: + with pysam.AlignmentFile(bam_file, "rb") as inf: + for read in unplaced_reads(inf): + tags = read_tags.get(read.query_name) if read_tags else None + if tags is not None: + _apply_tags(read, tags, barcode_tag, umi_tag) + outf.write(read) + written += 1 + return written diff --git a/isoquant_lib/utils/file_naming.py b/isoquant_lib/utils/file_naming.py index 5e05d994..943b8854 100644 --- a/isoquant_lib/utils/file_naming.py +++ b/isoquant_lib/utils/file_naming.py @@ -29,6 +29,12 @@ def split_barcodes_lock_filename(sample): return sample.barcodes_split_reads + "_lock" +def tagged_bam_lock_filename(sample): + # in the aux dir, so it outlives clean_up and a resumed run does not copy the whole + # input BAM again for an output that is already complete + return sample.barcodes_split_reads + "_tagged_bam_done" + + def clean_locks(chr_ids, base_name, fname_function): for chr_id in chr_ids: fname = fname_function(base_name, chr_id) @@ -48,6 +54,15 @@ def filtered_reads_file_name(out_raw_file: str, chr_id: str): return out_raw_file + "_filtered_" + chr_id +def dedup_bam_fragment_name(out_raw_file: str, chr_id: str): + # under out_raw_file so clean_up removes it even if the merge never happens + return out_raw_file + "_dedup_" + convert_chr_id_to_file_name_str(chr_id) + ".bam" + + +def tagged_bam_fragment_name(out_raw_file: str, chr_id: str): + return out_raw_file + "_tagged_" + convert_chr_id_to_file_name_str(chr_id) + ".bam" + + def umi_filtered_reads_file_name(out_umi_filtered_tmp: str, chr_id: str, edit_distance: int): return out_umi_filtered_tmp + ("_%s_ED%d" % (chr_id, edit_distance)) diff --git a/isoquant_lib/utils/file_utils.py b/isoquant_lib/utils/file_utils.py index 0ffe471f..2d8dbe40 100644 --- a/isoquant_lib/utils/file_utils.py +++ b/isoquant_lib/utils/file_utils.py @@ -5,6 +5,7 @@ # See file LICENSE for details. ############################################################################ +import gzip import logging import os import re @@ -17,6 +18,85 @@ logger = logging.getLogger('IsoQuant') +GZIP_SUFFIX = ".gz" + +# Levels for the two kinds of output IsoQuant writes, measured on real ONT data (see below). +# Python's gzip defaults to 9, which is a bad trade for both. +# +# Tables (TSV, BED, MTX, allinfo): 6. On a real barcode table level 9 runs at 16 MB/s against +# 39 MB/s at 6 for 3% less output, and going below 6 saves little -- level 5 is 11% faster for +# 0.4% more, level 1 is 2.4x faster for 8.5% more. +GZIP_LEVEL = 6 +# Sequences (FASTA/FASTQ): 4. Nucleotide data sits near gzip's entropy floor, so the high +# levels grind: 12.5 MB/s at 6 against 65.3 MB/s at 4, for 5.7% more output. These are also +# the largest files IsoQuant writes. +GZIP_LEVEL_SEQUENCES = 4 + +SEQUENCE_SUFFIXES = (".fa", ".fasta", ".fq", ".fastq") + + +def gzip_level_for(file_name): + """Compression level for an output, chosen by the kind of data its name implies.""" + if strip_compression_suffix(file_name).endswith(SEQUENCE_SUFFIXES): + return GZIP_LEVEL_SEQUENCES + return GZIP_LEVEL + + +def open_text_write(file_name, compresslevel=None): + """Open for text writing, compressing when the name says so. + + The level defaults to what the name implies; pass one to override. + """ + if file_name.endswith(GZIP_SUFFIX): + if compresslevel is None: + compresslevel = gzip_level_for(file_name) + return gzip.open(file_name, "wt", compresslevel=compresslevel) + return open(file_name, "w") + + +def open_text_read(file_name): + """Open for text reading, decompressing when the name says so.""" + if file_name.endswith(GZIP_SUFFIX): + return gzip.open(file_name, "rt") + return open(file_name, "r") + + +def resolve_optionally_gzipped(file_name): + """Return the existing path among and .gz. + + Outputs that are compressed once the run finishes are still referred to by their plain + name (in resumed runs, for instance), so readers have to accept either. + """ + if os.path.exists(file_name): + return file_name + gzipped = file_name + GZIP_SUFFIX + if os.path.exists(gzipped): + return gzipped + return file_name + + +def gzip_file_in_place(file_name, keep_original=False): + """Compress a finished output to .gz. Returns the resulting path.""" + if file_name.endswith(GZIP_SUFFIX): + return file_name + if not os.path.exists(file_name): + return file_name + gzipped = file_name + GZIP_SUFFIX + with open(file_name, "rb") as inf, \ + gzip.open(gzipped, "wb", compresslevel=gzip_level_for(file_name)) as outf: + shutil.copyfileobj(inf, outf) + if not keep_original: + os.remove(file_name) + return gzipped + + +def strip_compression_suffix(file_name): + """Drop a trailing compression suffix so extension logic sees the real one.""" + for suffix in (GZIP_SUFFIX, ".gzip", ".bgz"): + if file_name.endswith(suffix): + return file_name[:-len(suffix)] + return file_name + def check_file_exists(file_path: str, description: str): """Check that a file exists, exit with error if not.""" diff --git a/isoquant_lib/utils/input_data_storage.py b/isoquant_lib/utils/input_data_storage.py index 8c66061c..fc751c81 100644 --- a/isoquant_lib/utils/input_data_storage.py +++ b/isoquant_lib/utils/input_data_storage.py @@ -96,6 +96,8 @@ def _init_paths(self): self.out_cell_barcodes_tsv = self._make_path(self.prefix + ".cell_barcodes.tsv") self.out_cell_barcodes_stats = self._make_path(self.prefix + ".cell_barcodes.stats") self.out_umi_filtered = self._make_path(self.prefix + ".UMI_filtered") + self.out_tagged_bam = self._make_path(self.prefix + ".tagged.bam") + self.out_deduplicated_bam = self._make_path(self.prefix + ".deduplicated.bam") self.out_umi_filtered_tmp = self._make_aux_path(self.prefix + ".UMI_filtered") self.out_umi_filtered_done = self._make_aux_path(self.prefix + ".UMI_filtered.done") self.split_reads_fasta = self._make_path(self.prefix + ".split_reads") diff --git a/isoquant_lib/utils/string_pools.py b/isoquant_lib/utils/string_pools.py index cfd64f5b..356c454b 100644 --- a/isoquant_lib/utils/string_pools.py +++ b/isoquant_lib/utils/string_pools.py @@ -37,6 +37,7 @@ import logging from typing import List, Dict +from isoquant_lib.utils.file_utils import strip_compression_suffix from isoquant_lib.assignment.read_groups import AbstractReadGrouper, get_grouping_pool_types from isoquant_lib.assignment.assignment_loader import load_genedb from .serialization import write_int, write_string, read_int, read_string @@ -238,7 +239,7 @@ def build_file_name_pool(self, sample): # file_list is a list of lists (libraries), each library has one or more files for lib in sample.file_list: # Get basename without extension, matching FileNameGrouper logic - 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] file_names.add(readable_name) # Add in sorted order for deterministic IDs diff --git a/isoquant_tests/github/configs/GROUP12.SC.SIRVs.R10.yaml b/isoquant_tests/github/configs/GROUP12.SC.SIRVs.R10.yaml index 33a93ad5..90765064 100644 --- a/isoquant_tests/github/configs/GROUP12.SC.SIRVs.R10.yaml +++ b/isoquant_tests/github/configs/GROUP12.SC.SIRVs.R10.yaml @@ -12,7 +12,7 @@ bam: /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc. /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.group_SIRV6.bam /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.group_SIRV7.bam datatype: ont -isoquant_options: '"-t 8 --complete_genedb --force --mode tenX_v3 --read_group barcode --barcoded_reads +isoquant_options: '"-t 8 --complete_genedb --force --mode tenX_v3 --split_molecules false --read_group barcode --barcoded_reads /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.barcodes.tsv --barcode2spot /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.barcode2spot.tsv "' diff --git a/isoquant_tests/github/configs/GROUP13.SC.SIRVs.R10.yaml b/isoquant_tests/github/configs/GROUP13.SC.SIRVs.R10.yaml index 6d79f316..145c822a 100644 --- a/isoquant_tests/github/configs/GROUP13.SC.SIRVs.R10.yaml +++ b/isoquant_tests/github/configs/GROUP13.SC.SIRVs.R10.yaml @@ -12,7 +12,7 @@ bam: /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc. /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.group_SIRV6.bam /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.group_SIRV7.bam datatype: ont -isoquant_options: '"-t 8 --complete_genedb --force --mode tenX_v3 --read_group barcode --barcoded_reads +isoquant_options: '"-t 8 --complete_genedb --force --mode tenX_v3 --split_molecules false --read_group barcode --barcoded_reads /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.barcodes.3cols.tsv --barcode2spot /abga/work/andreyp/ci_isoquant/data/groups/mix/SRIV.simulated_R10_no_trunc.barcode2spot.multicolumn.tsv:0:1,2,3 "' diff --git a/isoquant_tests/github/configs/SC.Mouse.10x.allinfo.yaml b/isoquant_tests/github/configs/SC.Mouse.10x.allinfo.yaml index 451a699c..a43f0b08 100644 --- a/isoquant_tests/github/configs/SC.Mouse.10x.allinfo.yaml +++ b/isoquant_tests/github/configs/SC.Mouse.10x.allinfo.yaml @@ -5,8 +5,8 @@ genome: /abga/work/andreyp/data/reference/mouse/GRCm39.primary_assembly.genome.f genedb: /abga/work/andreyp/data/reference/mouse/gencode.vM36.basic.annotation.gtf bam: /abga/work/andreyp/ci_isoquant/data/barcodes/Mouse.10x.5k.ONT_cDNA.R10.4.no_trunc.bam datatype: ont -isoquant_options: '"-m tenX_v3 --complete_genedb --barcode_whitelist /abga/work/andreyp/ci_isoquant/data/barcodes/10xMultiome_5K.tsv - --large_output allinfo read2transcripts -t 10"' +isoquant_options: '"-m tenX_v3 --split_molecules false --complete_genedb --barcode_whitelist /abga/work/andreyp/ci_isoquant/data/barcodes/10xMultiome_5K.tsv + --large_output allinfo read2transcripts tagged_bam deduplicated_bam -t 10"' edit_distance: '3' baselines: allinfo: diff --git a/isoquant_tests/github/configs/SC.Mouse.10x.barcoded_bam.allinfo.yaml b/isoquant_tests/github/configs/SC.Mouse.10x.barcoded_bam.allinfo.yaml index 9a8f8d13..3913174a 100644 --- a/isoquant_tests/github/configs/SC.Mouse.10x.barcoded_bam.allinfo.yaml +++ b/isoquant_tests/github/configs/SC.Mouse.10x.barcoded_bam.allinfo.yaml @@ -5,8 +5,8 @@ genome: /abga/work/andreyp/data/reference/mouse/GRCm39.primary_assembly.genome.f genedb: /abga/work/andreyp/data/reference/mouse/gencode.vM36.basic.annotation.gtf bam: /abga/work/andreyp/ci_isoquant/data/barcodes/Mouse.10x.5k.ONT_cDNA.R10.4.no_trunc.tagged.bam datatype: ont -isoquant_options: '"-m tenX_v3 --complete_genedb --barcoded_bam --strip_barcode_suffix - --large_output allinfo read2transcripts -t 10"' +isoquant_options: '"-m tenX_v3 --split_molecules false --complete_genedb --barcoded_bam --strip_barcode_suffix + --large_output allinfo read2transcripts deduplicated_bam -t 10"' edit_distance: '3' baselines: allinfo: diff --git a/isoquant_tests/test_bam_utils.py b/isoquant_tests/test_bam_utils.py new file mode 100644 index 00000000..ca6b50ca --- /dev/null +++ b/isoquant_tests/test_bam_utils.py @@ -0,0 +1,333 @@ +############################################################################ +# Copyright (c) 2025-2026 University of Helsinki +# All Rights Reserved +# See file LICENSE for details. +############################################################################ + +"""BAM outputs: per-chromosome tagging and merging.""" + +import os +import shutil + +import pysam +import pytest + +import isoquant_lib.utils.bam_utils as bam_utils +from isoquant_lib.utils.bam_utils import ( + GENE_TAG, + TRANSCRIPT_TAG, + collect_unmapped_read_ids, + load_barcode_umi_tags, + load_survivor_tags, + merge_bam_files, + references_with_alignments, + write_tagged_chromosome_bam, + write_unmapped_bam, +) + +HEADER = {"HD": {"VN": "1.0"}, + "SQ": [{"SN": "chr1", "LN": 10000}, {"SN": "chr2", "LN": 10000}]} + + +def make_read(name, chr_index, start=100, secondary=False, supplementary=False, unmapped=False): + read = pysam.AlignedSegment() + read.query_name = name + read.query_sequence = "ACGT" * 10 + read.query_qualities = pysam.qualitystring_to_array("I" * 40) + if unmapped: + read.is_unmapped = True + read.reference_id = -1 + read.reference_start = -1 + else: + read.reference_id = chr_index + read.reference_start = start + read.mapping_quality = 60 + read.cigarstring = "40M" + read.is_secondary = secondary + read.is_supplementary = supplementary + return read + + +def write_bam(path, reads): + with pysam.AlignmentFile(path, "wb", header=HEADER) as outf: + for read in reads: + outf.write(read) + pysam.index(path) + return path + + +def read_names(path): + with pysam.AlignmentFile(path, "rb") as inf: + return [read.query_name for read in inf.fetch(until_eof=True)] + + +@pytest.fixture +def input_bam(tmp_path): + """chr1 has a primary, a secondary and a supplementary record for r1, plus r2 and r3.""" + reads = [ + make_read("r1", 0, 100), + make_read("r1", 0, 200, secondary=True), + make_read("r1", 0, 300, supplementary=True), + make_read("r2", 0, 400), + make_read("r3", 0, 500), + make_read("r4", 1, 100), + make_read("unmapped1", 0, unmapped=True), + ] + return write_bam(str(tmp_path / "in.bam"), reads) + + +class TestWriteTaggedChromosomeBam: + def test_tags_are_applied(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + tags = {"r1": ("ACGTACGT", "TTTTTT", "GENE1", "TX1")} + write_tagged_chromosome_bam([input_bam], out, "chr1", tags) + with pysam.AlignmentFile(out, "rb") as inf: + tagged = [r for r in inf.fetch(until_eof=True) if r.query_name == "r1"] + assert tagged, "r1 must be present" + for read in tagged: + assert read.get_tag("CB") == "ACGTACGT" + assert read.get_tag("UB") == "TTTTTT" + assert read.get_tag(GENE_TAG) == "GENE1" + assert read.get_tag(TRANSCRIPT_TAG) == "TX1" + + def test_custom_tag_names(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + write_tagged_chromosome_bam([input_bam], out, "chr1", {"r2": ("AAAA", "CCCC", None, None)}, + barcode_tag="XC", umi_tag="XM") + with pysam.AlignmentFile(out, "rb") as inf: + read = next(r for r in inf.fetch(until_eof=True) if r.query_name == "r2") + assert read.get_tag("XC") == "AAAA" + assert read.get_tag("XM") == "CCCC" + assert not read.has_tag("CB") + + @pytest.mark.parametrize("placeholder", ["*", ".", "None"]) + def test_placeholders_are_not_written(self, input_bam, tmp_path, placeholder): + """'*' means an uncalled barcode, '.' an unset feature, and IsoQuant writes the + literal 'None' for the transcript of a novel or ambiguous read.""" + out = str(tmp_path / "out.bam") + write_tagged_chromosome_bam([input_bam], out, "chr1", + {"r2": (placeholder,) * 4}) + with pysam.AlignmentFile(out, "rb") as inf: + read = next(r for r in inf.fetch(until_eof=True) if r.query_name == "r2") + assert not read.has_tag("CB") + assert not read.has_tag("UB") + assert not read.has_tag(GENE_TAG) + assert not read.has_tag(TRANSCRIPT_TAG) + + def test_keeps_untagged_reads_by_default(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + written = write_tagged_chromosome_bam([input_bam], out, "chr1", {"r1": ("A", "C", "G", "T")}) + assert written == 5 # every chr1 alignment, tagged or not + assert sorted(set(read_names(out))) == ["r1", "r2", "r3"] + + def test_subset_when_untagged_are_dropped(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + written = write_tagged_chromosome_bam([input_bam], out, "chr1", {"r2": ("A", "C", "G", "T")}, + keep_untagged=False) + assert written == 1 + assert read_names(out) == ["r2"] + + def test_primary_only_drops_secondary_and_supplementary(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + written = write_tagged_chromosome_bam([input_bam], out, "chr1", {"r1": ("A", "C", "G", "T")}, + primary_only=True, keep_untagged=False) + assert written == 1 + with pysam.AlignmentFile(out, "rb") as inf: + read = next(inf.fetch(until_eof=True)) + assert not read.is_secondary and not read.is_supplementary + + def test_other_chromosomes_are_untouched(self, input_bam, tmp_path): + out = str(tmp_path / "out.bam") + write_tagged_chromosome_bam([input_bam], out, "chr2", {}) + assert read_names(out) == ["r4"] + + def test_records_are_copied_verbatim(self, input_bam, tmp_path): + """Fusion detection reads SA tags and soft clips, so nothing but tags may change.""" + source = pysam.AlignedSegment() + source.query_name = "sa_read" + source.query_sequence = "ACGT" * 10 + source.query_qualities = pysam.qualitystring_to_array("I" * 40) + source.reference_id = 0 + source.reference_start = 700 + source.mapping_quality = 60 + source.cigarstring = "10S30M" + source.set_tag("SA", "chr2,100,+,30M10S,60,0;") + bam = write_bam(str(tmp_path / "sa.bam"), [source]) + + out = str(tmp_path / "out.bam") + write_tagged_chromosome_bam([bam], out, "chr1", {"sa_read": ("A", "C", "G", "T")}, + primary_only=True, keep_untagged=False) + with pysam.AlignmentFile(out, "rb") as inf: + read = next(inf.fetch(until_eof=True)) + assert read.get_tag("SA") == "chr2,100,+,30M10S,60,0;" + assert read.cigarstring == "10S30M" + assert read.query_sequence == "ACGT" * 10 + + def test_several_input_bams(self, tmp_path): + first = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + second = write_bam(str(tmp_path / "b.bam"), [make_read("r2", 0)]) + out = str(tmp_path / "out.bam") + assert write_tagged_chromosome_bam([first, second], out, "chr1", {}) == 2 + assert sorted(read_names(out)) == ["r1", "r2"] + + +class TestReferencesWithAlignments: + def test_only_references_that_carry_reads(self, input_bam, tmp_path): + """chr2 has one read, so both chromosomes are listed; an empty one would not be.""" + assert references_with_alignments([input_bam]) == ["chr1", "chr2"] + + def test_empty_references_are_dropped(self, tmp_path): + """A fragmented assembly must not spawn a task per empty contig.""" + bam = write_bam(str(tmp_path / "one.bam"), [make_read("r1", 0)]) + assert references_with_alignments([bam]) == ["chr1"] + + def test_union_across_bams(self, tmp_path): + first = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + second = write_bam(str(tmp_path / "b.bam"), [make_read("r4", 1)]) + assert references_with_alignments([first, second]) == ["chr1", "chr2"] + + def test_covers_references_isoquant_would_skip(self, input_bam): + """The regression: scaffolds absent from the analysed chromosome list still have reads.""" + assert "chr2" in references_with_alignments([input_bam]) + + +class TestWriteUnmappedBam: + def test_only_unmapped_reads(self, input_bam, tmp_path): + """fetch() by chromosome never returns these, so they need their own pass.""" + out = str(tmp_path / "unmapped.bam") + assert write_unmapped_bam([input_bam], out) == 1 + assert read_names(out) == ["unmapped1"] + + def test_unmapped_reads_get_tagged(self, input_bam, tmp_path): + """A barcode is called from the sequence, so an unaligned read still has one. + + These reads belong to no chromosome and appear in no split table; copying them + untagged silently dropped the tags of every unmapped read. + """ + out = str(tmp_path / "unmapped.bam") + write_unmapped_bam([input_bam], out, {"unmapped1": ("ACGT", "TTTT", None, None)}) + with pysam.AlignmentFile(out, "rb") as inf: + read = next(inf.fetch(until_eof=True)) + assert read.get_tag("CB") == "ACGT" + assert read.get_tag("UB") == "TTTT" + + def test_untagged_unmapped_reads_are_still_copied(self, input_bam, tmp_path): + out = str(tmp_path / "unmapped.bam") + assert write_unmapped_bam([input_bam], out, {"someone_else": ("A", "C", None, None)}) == 1 + with pysam.AlignmentFile(out, "rb") as inf: + assert not next(inf.fetch(until_eof=True)).has_tag("CB") + + +class TestCollectUnmappedReadIds: + def test_finds_the_unplaced_reads(self, input_bam): + assert collect_unmapped_read_ids([input_bam]) == {"unmapped1"} + + def test_none_when_all_aligned(self, tmp_path): + bam = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + assert collect_unmapped_read_ids([bam]) == set() + + +class TestUnplacedReads: + def test_same_reads_with_and_without_an_index(self, input_bam, tmp_path): + """An index turns reaching the unplaced reads into a seek; the result must not change.""" + unindexed = str(tmp_path / "noindex.bam") + shutil.copyfile(input_bam, unindexed) + with pysam.AlignmentFile(input_bam, "rb") as indexed: + assert indexed.has_index() + with_index = [read.query_name for read in bam_utils.unplaced_reads(indexed)] + with pysam.AlignmentFile(unindexed, "rb") as plain: + assert not plain.has_index() + without_index = [read.query_name for read in bam_utils.unplaced_reads(plain)] + assert with_index == without_index == ["unmapped1"] + + +class TestMergeBamFiles: + def test_merges_and_indexes(self, tmp_path): + first = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + second = write_bam(str(tmp_path / "b.bam"), [make_read("r4", 1)]) + out = str(tmp_path / "merged.bam") + assert merge_bam_files(out, [first, second]) == out + assert sorted(read_names(out)) == ["r1", "r4"] + assert os.path.exists(out + ".bai") or os.path.exists(out + ".csi") + assert not os.path.exists(first) and not os.path.exists(second) + + def test_single_fragment_is_moved(self, tmp_path): + only = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + out = str(tmp_path / "merged.bam") + assert merge_bam_files(out, [only]) == out + assert not os.path.exists(only) + assert read_names(out) == ["r1"] + + def test_missing_and_none_fragments_are_skipped(self, tmp_path): + """Chromosomes with no surviving reads contribute None.""" + only = write_bam(str(tmp_path / "a.bam"), [make_read("r1", 0)]) + out = str(tmp_path / "merged.bam") + assert merge_bam_files(out, [None, only, str(tmp_path / "gone.bam")]) == out + assert read_names(out) == ["r1"] + + def test_nothing_to_merge(self, tmp_path): + out = str(tmp_path / "merged.bam") + assert merge_bam_files(out, [None, None]) is None + assert not os.path.exists(out) + + @pytest.mark.parametrize("count", [11, 12, 13]) + def test_more_fragments_than_one_merge_call_takes(self, tmp_path, monkeypatch, count): + """samtools opens every input at once, so a reference with many contigs needs rounds. + + The counts straddle a multiple of the batch size, which is where a lone leftover + fragment has to be carried into the next round. + """ + monkeypatch.setattr(bam_utils, "BAM_MERGE_BATCH", 3) + fragments = [write_bam(str(tmp_path / ("f%d.bam" % i)), [make_read("r%d" % i, 0, 100 + i)]) + for i in range(count)] + out = str(tmp_path / "merged.bam") + assert merge_bam_files(out, fragments) == out + assert sorted(read_names(out)) == sorted("r%d" % i for i in range(count)) + assert not any(os.path.exists(f) for f in fragments) + # the intermediates of every round are cleaned up as well + assert sorted(os.path.basename(f) for f in tmp_path.glob("merged.bam*")) == \ + [os.path.basename(out), os.path.basename(out) + ".bai"] + + +class TestLoadSurvivorTags: + def test_with_tag_columns(self, tmp_path): + path = tmp_path / "s.tsv" + path.write_text("r1\tACGT\tTTTT\tGENE1\tTX1\nr2\tCCCC\tGGGG\t.\t.\n") + assert load_survivor_tags(str(path)) == { + "r1": ("ACGT", "TTTT", "GENE1", "TX1"), + "r2": ("CCCC", "GGGG", ".", "."), + } + + def test_read_ids_only(self, tmp_path): + """Without deduplicated_bam the table is a bare id list; the ids still define the subset.""" + path = tmp_path / "s.tsv" + path.write_text("r1\nr2\n") + assert load_survivor_tags(str(path)) == {"r1": (None, None, None, None), + "r2": (None, None, None, None)} + + def test_missing_file(self, tmp_path): + assert load_survivor_tags(str(tmp_path / "nope.tsv")) == {} + + +class TestLoadBarcodeUmiTags: + def test_reads_the_split_table(self, tmp_path): + path = tmp_path / "b.tsv" + path.write_text("#read_id\tbarcode\tUMI\nr1\tACGT\tTTTT\nr2\tCCCC\tGGGG\n") + assert load_barcode_umi_tags(str(path)) == {"r1": ("ACGT", "TTTT", None, None), + "r2": ("CCCC", "GGGG", None, None)} + + def test_filtering_by_read_ids(self, tmp_path): + """The whole-sample table is scanned for a handful of unmapped reads.""" + path = tmp_path / "b.tsv" + path.write_text("r1\tACGT\tTTTT\nr2\tCCCC\tGGGG\nr3\tAAAA\tCCCC\n") + assert load_barcode_umi_tags(str(path), read_ids={"r2"}) == { + "r2": ("CCCC", "GGGG", None, None)} + assert load_barcode_umi_tags(str(path), read_ids=set()) == {} + + def test_short_lines_are_skipped(self, tmp_path): + path = tmp_path / "b.tsv" + path.write_text("r1\tACGT\tTTTT\nbroken\n") + assert list(load_barcode_umi_tags(str(path))) == ["r1"] + + def test_missing_file(self, tmp_path): + assert load_barcode_umi_tags(str(tmp_path / "nope.tsv")) == {} diff --git a/isoquant_tests/test_chunk_ordering.py b/isoquant_tests/test_chunk_ordering.py new file mode 100644 index 00000000..c7d65e38 --- /dev/null +++ b/isoquant_tests/test_chunk_ordering.py @@ -0,0 +1,60 @@ +############################################################################ +# Copyright (c) 2025-2026 University of Helsinki +# All Rights Reserved +# See file LICENSE for details. +############################################################################ + +"""Chunk results reach the caller in completion order, carrying the index to reorder by.""" + +import time +from types import SimpleNamespace + +from isoquant_lib.barcode_calling.detect_barcodes import run_chunks_in_parallel + +CHUNKS = 6 +THREADS = 3 + + +def finish_in_reverse(index: int) -> int: + """Make later chunks finish first, so completion order is not submission order.""" + time.sleep(0.05 * (CHUNKS - index)) + return index + + +class TestRunChunksInParallel: + def test_every_chunk_is_handled_once_with_its_index(self): + completion_order = [] + collected = {} + + def submit(pool, chunk, num): + return pool.submit(finish_in_reverse, num) + + def handle_result(result, chunk_index): + completion_order.append(chunk_index) + collected[chunk_index] = result + return 1 + + run_chunks_in_parallel(iter(range(CHUNKS)), SimpleNamespace(threads=THREADS), + None, submit, handle_result) + + assert sorted(collected) == list(range(CHUNKS)) + # the index really identifies its own chunk, so a caller can key on it + assert all(index == result for index, result in collected.items()) + # and it was needed: the chunks did not come back in the order they went out + assert completion_order != sorted(completion_order) + + def test_submission_order_is_recoverable(self): + """What _process_single_file_in_parallel does to keep its merge deterministic.""" + outputs = {} + + def submit(pool, chunk, num): + return pool.submit(finish_in_reverse, num) + + def handle_result(result, chunk_index): + outputs[chunk_index] = "chunk_%d" % result + return 1 + + run_chunks_in_parallel(iter(range(CHUNKS)), SimpleNamespace(threads=THREADS), + None, submit, handle_result) + + assert [outputs[i] for i in sorted(outputs)] == ["chunk_%d" % i for i in range(CHUNKS)] diff --git a/isoquant_tests/test_file_compression.py b/isoquant_tests/test_file_compression.py new file mode 100644 index 00000000..758e030a --- /dev/null +++ b/isoquant_tests/test_file_compression.py @@ -0,0 +1,173 @@ +############################################################################ +# Copyright (c) 2025-2026 University of Helsinki +# All Rights Reserved +# See file LICENSE for details. +############################################################################ + +"""Compression helpers for large single-cell outputs.""" + +import gzip +import os +import random + +import pytest + +from isoquant_lib.barcode_calling.detect_barcodes import numbered_chunk_name +from isoquant_lib.utils.file_utils import ( + gzip_file_in_place, + open_text_read, + GZIP_LEVEL, + GZIP_LEVEL_SEQUENCES, + gzip_level_for, + open_text_write, + resolve_optionally_gzipped, + strip_compression_suffix, +) + + +class TestOpenHelpers: + def test_plain_round_trip(self, tmp_path): + path = str(tmp_path / "x.tsv") + with open_text_write(path) as handle: + handle.write("a\nb\n") + assert not os.path.exists(path + ".gz") + with open_text_read(path) as handle: + assert handle.read() == "a\nb\n" + + def test_gzip_round_trip(self, tmp_path): + path = str(tmp_path / "x.tsv.gz") + with open_text_write(path) as handle: + handle.write("a\nb\n") + # really compressed, not just named .gz + with gzip.open(path, "rt") as handle: + assert handle.read() == "a\nb\n" + with open_text_read(path) as handle: + assert handle.read() == "a\nb\n" + + +class TestStripCompressionSuffix: + @pytest.mark.parametrize("name, expected", [ + ("S.split_reads_0.fa.gz", "S.split_reads_0.fa"), + ("S.split_reads_0.fa", "S.split_reads_0.fa"), + ("reads.fastq.gzip", "reads.fastq"), + ("reads.fastq.bgz", "reads.fastq"), + ("reads.bam", "reads.bam"), + ]) + def test_suffixes(self, name, expected): + assert strip_compression_suffix(name) == expected + + def test_read_group_name_is_unaffected_by_compression(self): + """File-name read groups must not change when an input gets compressed.""" + plain = os.path.splitext(strip_compression_suffix("S.split_reads_0.fa"))[0] + gzipped = os.path.splitext(strip_compression_suffix("S.split_reads_0.fa.gz"))[0] + assert plain == gzipped == "S.split_reads_0" + + +class TestResolveOptionallyGzipped: + def test_prefers_the_plain_file(self, tmp_path): + plain = tmp_path / "t.tsv" + plain.write_text("x") + (tmp_path / "t.tsv.gz").write_bytes(b"") + assert resolve_optionally_gzipped(str(plain)) == str(plain) + + def test_falls_back_to_gz(self, tmp_path): + """A resumed run refers to the plain name after the table was compressed.""" + gzipped = tmp_path / "t.tsv.gz" + gzipped.write_bytes(b"") + assert resolve_optionally_gzipped(str(tmp_path / "t.tsv")) == str(gzipped) + + def test_missing_returns_the_original(self, tmp_path): + missing = str(tmp_path / "nope.tsv") + assert resolve_optionally_gzipped(missing) == missing + + +class TestGzipLevelFor: + """Sequence data sits near gzip's entropy floor, so it gets a lower level than tables do.""" + + @pytest.mark.parametrize("name", ["reads.fa.gz", "reads.fasta.gz", "reads.fq.gz", + "reads.fastq.gz", "S.split_reads_0.fa.gz", + "plain.fasta"]) + def test_sequences(self, name): + assert gzip_level_for(name) == GZIP_LEVEL_SEQUENCES + + @pytest.mark.parametrize("name", ["counts.tsv.gz", "reads.bed.gz", "S.allinfo.gz", + "m.mtx.gz", "S.barcoded_reads_0.tsv"]) + def test_tables(self, name): + assert gzip_level_for(name) == GZIP_LEVEL + + def test_the_level_reaches_gzip(self, tmp_path): + """Sizes, because a GzipFile does not expose the level it was built with.""" + rnd = random.Random(0) + payload = ">r\n" + "".join(rnd.choice("ACGT") for _ in range(200000)) + "\n" + + def written(name, **kwargs): + path = tmp_path / name + with open_text_write(str(path), **kwargs) as handle: + handle.write(payload) + return os.path.getsize(str(path)) + + by_name = written("a.fa.gz") + assert by_name == written("b.fa.gz", compresslevel=GZIP_LEVEL_SEQUENCES) + assert by_name > written("c.fa.gz", compresslevel=9) + # a table of the same bytes takes the table level instead + assert written("d.tsv.gz") == written("e.tsv.gz", compresslevel=GZIP_LEVEL) + + +class TestGzipFileInPlace: + def test_compresses_and_removes_the_original(self, tmp_path): + path = tmp_path / "t.tsv" + path.write_text("read1\tACGT\tTTTT\n") + result = gzip_file_in_place(str(path)) + assert result == str(path) + ".gz" + assert not path.exists() + with gzip.open(result, "rt") as handle: + assert handle.read() == "read1\tACGT\tTTTT\n" + + def test_keep_original(self, tmp_path): + path = tmp_path / "t.tsv" + path.write_text("x\n") + gzip_file_in_place(str(path), keep_original=True) + assert path.exists() + + def test_already_compressed_is_a_no_op(self, tmp_path): + path = tmp_path / "t.tsv.gz" + with gzip.open(path, "wt") as handle: + handle.write("x\n") + assert gzip_file_in_place(str(path)) == str(path) + + def test_missing_file_is_a_no_op(self, tmp_path): + missing = str(tmp_path / "nope.tsv") + assert gzip_file_in_place(missing) == missing + assert not os.path.exists(missing + ".gz") + + +class TestNumberedChunkName: + def test_plain(self): + assert numbered_chunk_name("/tmp/subreads", 3) == "/tmp/subreads_3" + + def test_a_fasta_chunk_is_still_a_fasta(self): + """The chunk inherits the level from its name, so the extension has to survive.""" + assert numbered_chunk_name("subreads.fa.gz", 3) == "subreads_3.fa.gz" + assert gzip_level_for(numbered_chunk_name("subreads.fa.gz", 3)) == GZIP_LEVEL_SEQUENCES + + def test_keeps_the_compression_suffix_last(self): + """Otherwise the per-chunk temp would not be recognised as compressed.""" + assert numbered_chunk_name("/tmp/subreads.gz", 3) == "/tmp/subreads_3.gz" + + +class TestMultiMemberConcatenation: + def test_byte_concatenated_members_read_as_one_stream(self, tmp_path): + """Chunks are compressed in the workers and merged with a byte copy. + + That only works because concatenated gzip members form a valid stream. + """ + merged = tmp_path / "merged.fa.gz" + with open(merged, "wb") as out: + for chunk in range(3): + part = tmp_path / ("part_%d.gz" % chunk) + with open_text_write(str(part)) as handle: + handle.write(">read%d\nACGT\n" % chunk) + with open(part, "rb") as inf: + out.write(inf.read()) + with gzip.open(merged, "rt") as handle: + assert handle.read() == ">read0\nACGT\n>read1\nACGT\n>read2\nACGT\n" diff --git a/isoquant_tests/test_split_molecules.py b/isoquant_tests/test_split_molecules.py index 714880b4..1ae7209c 100644 --- a/isoquant_tests/test_split_molecules.py +++ b/isoquant_tests/test_split_molecules.py @@ -10,6 +10,7 @@ import pytest +import isoquant from isoquant_lib.barcode_calling import options from isoquant_lib.modes import ( DEPRECATED_MODE_ALIASES, @@ -23,12 +24,16 @@ NON_SPLITTING_MODES = ["curio", "visium_hd", "custom_sc"] -def resolve(mode, split_molecules=None): - """Run the same two steps check_input_params does, and return (mode, split flag).""" - args = argparse.Namespace(mode=mode, split_molecules=split_molecules) +def resolve(mode, split_molecules=None, needs_mapping=True): + """Run the same steps check_input_params does, and return (mode, split flag).""" + input_data = argparse.Namespace( + input_type=argparse.Namespace(name="fastq" if needs_mapping else "bam", + needs_mapping=lambda: needs_mapping)) + args = argparse.Namespace(mode=mode, split_molecules=split_molecules, input_data=input_data) options.resolve_deprecated_mode(args) args.mode = IsoQuantMode[args.mode] options.resolve_split_molecules(args) + isoquant._reject_splitting_aligned_input(args) return args.mode, args.split_molecules @@ -93,6 +98,32 @@ def test_alias_defaults_survive_resume(self, alias): assert first is second +class TestAlignedInput: + """Aligned input and molecule splitting are contradictory requests. + + A BAM says "do not map"; splitting rewrites the reads so the pieces must be mapped again. + Guessing either way gives the user something they did not ask for, so both auto and true + abort. Splitting anyway used to replace file_list with a FASTA that was opened as a BAM. + """ + + @pytest.mark.parametrize("mode", SPLITTING_MODES) + @pytest.mark.parametrize("requested", [None, SPLIT_MOLECULES_AUTO, SPLIT_MOLECULES_TRUE]) + def test_splitting_aborts_on_aligned_input(self, mode, requested): + with pytest.raises(SystemExit) as excinfo: + resolve(mode, requested, needs_mapping=False) + assert excinfo.value.code != 0 + + @pytest.mark.parametrize("mode", SPLITTING_MODES) + def test_false_is_accepted_for_aligned_input(self, mode): + """The documented way out: analyse the alignments as they are.""" + assert resolve(mode, SPLIT_MOLECULES_FALSE, needs_mapping=False)[1] is False + + @pytest.mark.parametrize("mode", NON_SPLITTING_MODES) + def test_unsupported_modes_are_unaffected(self, mode): + """Nothing would have been split anyway, so there is no contradiction to report.""" + assert resolve(mode, needs_mapping=False)[1] is False + + class TestDeprecatedModeAliases: @pytest.mark.parametrize("alias, expected_mode, expected_split", [ ("tenX_v3_split", IsoQuantMode.tenX_v3, True),