Conversation
debug README.md
- Compare PeakATail with 19 single-cell/spatial APA tools - Document PeakATail's unique novelty: first tool to perform de novo unsupervised cell clustering on peak/PAS counts - Include detailed method-by-method analysis - Add feature comparison matrix - Provide literature evidence supporting novelty claim - List recommended improvements for publication Tools analyzed: scMAPA, SCINPAS, SCAPTURE, scTail, scDaPars, Sierra, scAPAtrap, SCAPE, scLAPA, spvAPA, MAAPER, scraps, scAPA, scDAPA, scPAISO, InPACT, SAPAS, vizAPA, stAPAminer Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
- Remove fixed timelines (to be set by user) - Add extensive database validation section - PolyASite 2.0/3.0 validation protocol - PolyA_DB v4 validation - APAeval benchmark integration - Include precision/recall metrics and targets - Add validation pipeline code examples - Document benchmarking strategy based on 19 tool analysis - Define publication figures and manuscript structure - List success criteria and journal targets Database validation resources: - PolyASite: Consolidated 3' end seq atlas - PolyA_DB v4: Systematic PAS identification - APAeval: Community benchmark platform - GENCODE: Annotation validation Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
- Add get_mapping() to indexing.py as public API for barcode mapping - Fix broken cb_total import: use get_mapping() instead - Fix mutable default arg in make_dataframe() (collist=None pattern) - Fix off-by-one guards: len(columns) >= 2 -> >= 3 (4 locations) - Keep sparse matrix: make_dataframe() returns (sparse_csc, pas_ids, collist) instead of dense DataFrame, preventing 4-8GB memory blowup - Rewrite preprocessing() to use scanpy sparse ops (filter_cells/filter_genes) - Fix variable name confusion: pas_count/cell_count -> peaks_per_cell/cells_per_peak - Lower min_genes default from 200 to 50 (APA has 10-100 PAS per cell) - Add --min-pas-per-cell CLI parameter - Fix mutable defaults in Peak.__init__ (peak.py) - Fix bare except in read.py (except: -> except KeyError:) Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
find_close.py:
- Replace intersect with closest (the main bug: PAS near but not
overlapping genes were silently dropped)
- Implement adaptive distance threshold using 3'UTR lengths per gene:
TIER_1 (within UTR), TIER_2 (within UTR x multiplier),
TIER_3 (within max_distance), INTERGENIC (flagged)
- Fix fragile "." filter to only check gene_id column
- Update column indices for closest output schema
gtftobed.py:
- Fix string-not-tuple: ("gene") -> ("gene",) and ("ENSG") -> ("ENSG",)
- Replace fragile positional GTF attribute parsing (biotype_index=17)
with proper key-value parsing of GTF attribute column
- Extract three_prime_utr features and save UTR lengths per gene
annotate.py:
- DO NOT overwrite countmatrix.index with range(1,N) -- preserve PAS IDs
- Fix join to use actual PAS IDs from sparse matrix
- Fix PAS BED filter to use actual PAS IDs
- Remove debug print statements
- Accept new sparse matrix API from make_dataframe()
gtf_cache.py (new):
- GTF processing cache using mtime+size fingerprinting
- Stores parsed results (gene BED, UTR lengths, features) in gtf_cache/
- process_gtf_cached() entry point for background thread pre-processing
main.py:
- Run GTF pre-processing in background thread parallel to peak_calling
- Pass UTR lengths and CLI params to find_close
- Wire up new make_dataframe/annotate/preprocessing API chain
CLI:
- Add --max-gene-distance, --utr-multiplier, --include-extended params
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
…st peak, CB counted past boundary A1: Use `chro` (old chromosome) instead of `chro1` (new chromosome) in the chromosome-change flush blocks, preventing peaks from being written with the wrong chromosome name. A2: Add final flush block after bamfile.close() to emit the last peak of the BAM file, which was silently dropped at EOF. A3: Guard cb_counting/cb_position_counting with `start1 <= l_end` check so reads past the peak boundary are not counted into the peak. E3: Change lambda_window default from 10000 to 5000 to match CLI default. Also close matrix and bedfile handles after the final flush. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
B1: Extract actual non-zero row indices from COO matrix instead of reconstructing sequential arange(1, N+1). This was the root cause of 3838 lost PAS — sequential IDs did not match actual MatrixMarket row indices, causing downstream annotation joins to fail silently. B2: Reset filtered_cb_list at the start of filter_cb() to prevent module-level mutable list from accumulating across multiple calls. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
C1: Add annotated_matrix path to directory_config (annotated_matrix.mtx) so annotate() writes to a separate file instead of overwriting filterd_matrix (filterdmatrix.mtx). Previously the CB-filtered matrix was silently destroyed during annotation. E2: Wire --min-pas-per-cell CLI argument to filter_config.min_pas_per_cell. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
E1: Remove print("bam") and print("seq") debug statements that were
executed on every import of the config module.
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
F: Create ema/output.py with OutputManager class that manages a hierarchical output directory structure (01_peak_calling through 08_differential plus gtf_cache) with save_stats() and save_run_config(). D1: Capture GTF thread exceptions in gtf_error dict and re-raise as RuntimeError after join(), preventing silent failures when GTF processing crashes in the background thread. G: Wire OutputManager into main.py pipeline to save intermediate stats JSON after each stage (peak_calling, cb_filter, gtf_annotation, pas_gene, annotated, preprocessing, clustering) and run_config.json at pipeline start. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
…idation framework Peak Calling: - 4 pluggable strategies (original, lambda_poisson, lambda_gradient, sierra_iterative) - cb_positions tracking for per-position per-CB read assignment - Peak region partitioning by midpoints, no reads lost - Dynamic streaming threshold via window-based local lambda (deque) - Background estimation: Poisson p-value with floor lambda Pipeline Fixes (18 bugs): - peackcalling: wrong chrom on flush, last peak dropped, CB boundary fix - matrixfilter: PAS IDs from COO not shape, mutable defaults, off-by-one - find_close: intersect→closest, adaptive UTR threshold (TIER 1/2/3) - annotate: PAS ID preservation, separate output path - pop_apa + pdui: sparse PAS ID matrix handling - groupby: pandas 2.x axis=1 fix - gtftobed: proper GTF attribute parsing Clustering: - TF-IDF + LSI (TruncatedSVD) + Leiden (replaces Louvain) - Clustering strategy registry (leiden_tfidf, leiden_libsize, external) - Label comparison (ARI, AMI, silhouette, marker peaks) - Gene expression comparison script with Sankey diagram Differential APA: - Fisher's exact test with BH FDR correction - PDUI calculation (proximal-distal usage index) - Delta-PDUI between clusters - Effect size output (odds ratio, delta-proportion) Validation: - Dual database: PolyASite 2.0 (569K) + PolyA_DB v4 (303K) - Precision-focused benchmarking at multiple distance cutoffs - Parameter sweep framework (13 configs tested) - Output validator (12 format/sanity/statistical checks) Infrastructure: - OutputManager with 8-stage hierarchical directory structure - GTF caching (background thread, mtime+size fingerprint) - Structured logging (PeakCallingLogger) - pyproject.toml (uv compatible) - Technical report (17 pages, 16 figures) Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
- Fix read.py: replace fragile query_name.split() with read.get_tag('RG')
matching PI's fix in main branch — works for any BAM source, not just SRA
- New ema/datasets/manager.py: DatasetManager prepares (dataset_id, bam_path)
list from YAML config; merge_before strategy tags reads with RG and merges
BAMs via pysam before peak calling
- New ema/datasets/pas_merge.py: merge_pas_beds() unifies PAS coordinates
across datasets using bedtools merge; concat_matrices() concatenates sparse
MTX files with PAS ID remapping for merge_after strategy
- YAML schema redesign: datasets list with id/merge_strategy/bams per sample;
backward-compatible with legacy --bamDir flag
- main.py: DatasetManager-driven loop over all (dataset_id, bam) pairs;
single-sample path unchanged, multi-sample routes through pas_merge
- Add PyYAML>=6.0 to requirments.txt
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Four merge strategies, selectable per dataset in YAML:
- before: merge BAMs into one tagged BAM, peak call once (more reads)
- after: peak call each BAM separately, merge PAS coordinates via bedtools
- none: single BAM, no merge
- atlas: snap called PAS to a reference atlas via bedtools closest
Each dataset (id) is a clustering unit. 4 BAMs grouped into 2 datasets =
2 cluster sets, each with its own h5ad. Cross-dataset PAS coordinates are
unified (via bedtools merge -s or atlas snap) so cluster comparison shares
the same feature space.
ema/datasets/pas_merge.py
- merge_pas_beds: strand-aware bedtools merge (-s), encodes
{dataset_id}::{pasnumber} in BED col 4, uses -o collapse to track
all original pasnumbers per merged interval, writes
dataset_id\told_pasnumber\tnew_pas_id mapping.
- concat_matrices: handles both legacy headerless COO (peackcalling.py
output) and MatrixMarket-headered files via leading-character peek.
Uses (dataset_id, old_pasnumber) -> new_pas_id lookup, sums duplicates
via tocsr/tocoo round-trip.
ema/datasets/atlas_snap.py (NEW)
- snap_beds_to_atlas: bedtools closest -s -d, distance threshold filter,
same mapping format as merge_pas_beds so concat_matrices works for both.
ema/main.py
- Per-dataset reset_index() to isolate CB column space
- CB list dump per (dataset_id, idx) via get_mapping()
- Atlas vs merge dispatch on directory_config.atlas
- Per-dataset clustering split: extract_per_dataset_mtx via scipy slice,
filter_cb on the dataset's columns, separate h5ad output per dataset.
- Single-sample path unchanged (regression-safe).
ema/matrixfilter.py
- filter_cb gains optional input_matrix_paths and cb_list params for
multi-sample mode while keeping legacy pos+neg behavior as default.
ema/cli.py / ema/config.py / example.yaml
- --atlas, --atlas-distance args + YAML keys
Smoke-tested merge_pas_beds + concat_matrices and snap_beds_to_atlas +
concat_matrices end-to-end with synthetic data; verified strand-aware
merging, distance filtering, and that legacy headerless MTX from
peackcalling.py is parsed correctly without losing the first row.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
End-to-end test on 14.7M-read BAM (treated as 2 datasets) revealed
several bugs that the synthetic smoke test missed:
read.py — RG tag fallback for non-tagged BAMs
When merge_strategy=none/after/atlas, BAMs aren't pre-tagged with RG,
so all cells fell back to sample_id="default" and collided across
datasets. Added module-level set_default_sample_id() that main.py
sets to the current dataset_id between peak_calling calls; read.py
uses it as the RG fallback.
pas_merge.py — integer PAS IDs in BED col 4
annotate.py joins genes (from BED col 4) against pas_ids (from MTX
rows). MTX uses integer rows; BED was writing "PAS_N" strings —
no match. Now writes plain integer IDs (1, 2, 3...) matching MTX
row numbers. _pas_sort_key updated to handle bare integer IDs.
atlas_snap.py — same integer ID fix
Atlas string IDs (e.g. polyasite_chr1_100_+) are preserved in a
new sidecar atlas_pas_id_lookup.tsv. BED col 4 + mapping use the
sequential integer that matches MTX rows.
main.py — split unified BED for find_close
find_close cats pos+neg BED files. In multi-sample we have only
the unified BED, so split it back by strand into pos/neg files
before the find_close call.
main.py — per-dataset clustering stats
output_mgr.save_stats has no slot for "clustering_{ds_id}" — write
a JSON next to each dataset's clusters.h5ad instead.
leiden_tfidf.py — guard against zero-cell clamp
When n_obs - 1 < 2, TruncatedSVD got n_components=-1 and crashed.
Now raises a clear error explaining the dataset is too small for
LSI (only triggers on tiny chromosome-subset test data).
Verified on full BAM (data/SRR8325947, 14.7M reads):
- 2 datasets (same BAM): each clusters to 12 clusters, 662 cells,
cluster sizes byte-identical between datasets (perfect determinism)
- Single-dataset legacy path: 728 cells, 12 clusters
- Total counts in multi-sample sampleA slice == raw peak calling
output (6,337,529 counts, 2,794,774 nnz — no data loss in
extract_per_dataset_mtx or concat_matrices)
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Detailed plan covering: - Multi-PAS PDUI with isoform-aware aggregation (3 strategies: classic, proportion, shannon) - NB regression for differential APA (pairwise + multi-condition) - Cross-dataset cluster matching (3 strategies: marker_overlap, mnn, jaccard) - Multi-sample wiring of switch test - Atlas e2e validation - BAM tagging speedup via samtools addreplacerg - Peak.pasnumber reset hygiene - Cell-gap diagnostic Strategy + registry pattern across all 4 dimensions to avoid if/else chains. HPC section with performance targets, memory budget, profiling discipline, and HPC skills mapped to each agent. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…gents)
Three new strategy registries (consistent with existing peak-calling /
clustering registries) and supporting isoform-level GTF parsing.
PDUI strategies (ema/quantification/strategies/):
classic — current 2-PAS PDUI ported to strategy interface
proportion — per-PAS proportion vector per gene per cell (sums to 1.0)
shannon — Shannon entropy of PAS usage distribution per gene
All support per_isoform / per_gene aggregation + isoform_collapse modes.
Differential APA strategies (ema/switch_test/strategies/):
fisher — wraps existing Fisher exact (no logic change)
nb_pairwise — NB GLM with cluster as binary predictor
nb_multi — NB GLM omnibus LRT across all clusters
Validated: KS p=0.69 (no inflation under null), 20/20 signal recovery,
Fisher confirmed 6x nominal alpha on identical null data.
Cross-dataset cluster matching (ema/clustering/cross_dataset/):
marker_overlap — top-N marker PAS Jaccard + Hungarian assignment
mnn — mutual NN in shared LSI space (re-embedded)
jaccard — direct cell-set Jaccard
Isoform-aware GTF + mapping:
ema/annotate/gtf2isoform_utr.py — per-transcript 3'UTR parser with
mmap + multiprocessing (chunked by chromosome), mtime+size cache
ema/quantification/pas_to_isoform.py — bedtools intersect + pandas
vectorized rank assignment, strand-aware transcript_pos
Cleanups (Agent D):
ema/datasets/manager.py — _tag_bam_with_rg now uses samtools addreplacerg
(multithreaded C path) with version >= 1.10 check
ema/countmatrix/peak.py — Peak.reset_pasnumber() classmethod for
per-dataset isolation in multi-sample mode
Tests: 105/105 passed across all 5 new test files
test_gtf2isoform_utr.py: 22
test_pas_to_isoform.py: 13
test_pdui_strategies.py: 36
test_nb_regression.py: 14
test_cluster_matching.py: 20
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
5-phase optimization plan targeting 8-10x speedup: - Phase 1: pysam threads=4 + bisect→SortedList (~3x free) - Phase 2: 2-bit CB encoding + batch resolve at peak emit (~2x) - Phase 3: 3-stage pipeline (Reader→Finder→Writer, bounded queues, batch=10000) - Phase 4: tile parallelism with per-worker isolation + main-process merge - Phase 0/5: cProfile + memory_profiler baseline + post-validation HPC skills explicitly mapped per phase: python-performance-optimization, python-parallel-data-streaming, python-resource-management, memory-profiling. Streaming bisect endpoint algorithm and strategy registry untouched. Per-worker isolation pattern (same as concat_matrices) used to manage Peak.pasnumber and BarcodeIndex global state under parallelism. Deferred to feature/peak-calling-perf branch after apa-completeness lands. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
The main `ema` pipeline now stops at clustering (per-dataset) +
cross-dataset matching. Heavy compute (PDUI + differential APA across
all clusters) moved to the existing `ema_switch` separate CLI command,
where the user controls scope:
- Which cluster pairs to test (--cell_combinations)
- Which marker subset to use (--marker-top-n N, default 200)
- Which strategy (--diff-method {fisher,nb_pairwise,nb_multi})
- Which PDUI methods (--pdui-method classic,proportion,shannon)
- Resource budget (--max-jobs, --per-worker-mb)
Standard scRNA workflow: cluster → top-N markers per cluster (Wilcoxon
via sc.tl.rank_genes_groups) → union → run NB/Fisher/PDUI on this
~1-2K PAS subset instead of 20K. Reduces NB pairwise scope from
~6 hours to ~5 minutes with no loss of biological information.
ema/utils/resource_manager.py — psutil-based ResourceManager
- get_n_jobs(per_worker_mb): bounded by free RAM + physical_cpus - reserve
- get_batch_size(per_item_mb): for streaming operations
- Conservative: 70% of free RAM, max physical_cpu - 2 workers, hard cap 16
- Falls back to /proc/meminfo if psutil missing
ema/quantification/marker_selector.py
- select_marker_pas(adata, top_n_per_cluster=200, method='wilcoxon')
- restrict_count_matrix(...) — filter to marker subset
- save_markers(...) — write subset to TSV
ema/switch_test/cli.py — full rewrite
- Reads clustered h5ad from `ema` output
- Marker selection via select_marker_pas
- Wraps strategy calls in joblib.parallel_backend(n_jobs=...) so
PDUI strategies that hardcode n_jobs=-1 honor the resource cap
- Wires get_pdui_strategy + get_diff_strategy registries
ema/main.py — revert PDUI+diff blocks (kept cross-dataset matching)
ema/annotate/gtf2isoform_utr.py — use ResourceManager for chunk pool size
ema/datasets/manager.py — fix utf-8 decode bug in samtools --version parse
requirments.txt — add psutil>=5.9
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…prefix Two follow-on fixes after end-to-end test: 1. ema_switch CLI was only loading pas_isoform_map for per_isoform aggregation, but ALL PDUI strategies (classic, proportion, shannon) need it for gene grouping. PDUI was returning 0 rows. Now loads the map whenever any PDUI method is requested, with auto-default to data/Homo_sapiens.GRCh38.99.gtf and emaout/pasbed.bed paths. 2. cross_dataset/jaccard.py _strip_prefix() only split on '#' (Cell Ranger convention). PeakATail uses '_' as the dataset prefix separator (sampleA_AAAA), so jaccard found 0 matches on identical-input test. Now tries '#' first, falls back to '_'. End-to-end verification on full BAM: - Test 1 (ema main): 3m38s, 12 clusters/dataset, byte-identical to baseline, cross-dataset matching: 12 canonical clusters with conf 1.0 - Test 2 (ema_switch fisher, 66 pairs, top-200 markers): 1m08s, PDUI classic 6620 rows, 56026 sig PAS (q<0.05) — Fisher inflation visible - Test 3 (ema_switch nb_pairwise, 3 pairs, top-200): 17s, PDUI proportion + shannon 415K + 406K rows, 486 sig PAS (NB more conservative) - Test 4 (ema_switch nb_multi omnibus, top-200): 48s — PDUI works, but nb_multi diff returned 0 rows (known issue, separate bug) - Test 5 (cross_dataset marker_overlap): 1.6s, 12 canonical, conf 1.0 - Test 6 (cross_dataset mnn): 1.6s, 12 canonical, conf > 1.0 (known normalization bug, separate) - Test 7 (cross_dataset jaccard, after fix): 12 canonical, conf 1.0 ResourceManager picks safe n_jobs (12 on 14-physical / 6 GB free box; was -1 = oversubscribed before). Known issues for next branch: - nb_multi returns 0 rows (likely min_cells_per_group filter too strict) - mnn match_confidence > 1.0 (not normalized) Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…ript
nb_multi: Per-cluster nonzero filter (>=10 cells with nonzero per cluster)
was too strict for sparse single-cell data — every PAS failed because at
least one of 12 clusters always had <10 nonzero cells. Replaced with
scRNA-standard filter:
- >= min_cells_per_group total nonzero cells (default 10)
- >= 2 clusters with any signal
On real data (1440 marker PAS, 12 clusters) now produces 1412 testable rows
(was 0). All p-values valid in [0,1].
mnn: vote_matrix / cluster_size_a gave confidence > 1 when k>1 MNN per cell.
Replaced with row-normalization (vote_matrix / row_totals) so each cluster's
similarity vector is a proper probability distribution summing to 1.
On real data: confidence now in [0.902, 1.0] (was up to 6.4).
scripts/validate_strategies.py — checks output invariants for all 9 strategies:
PDUI: classic (PDUI in [0,1]), proportion (in [0,1] + per (gene,cell) sums to 1),
shannon (entropy >= 0, normalized in [0,1])
Diff: pvalue/qvalue in [0,1] for fisher/nb_pairwise/nb_multi
Match: confidence in [0,1], canonical_cluster count matches expected for
identical-input regression test
Validation result on emaout/per_dataset/sampleA+B/clusters.h5ad:
29 / 29 checks PASS
- PDUI all 3: non-empty, all values in [0,1], proportion sums to 1.0 (99.9%)
- Diff all 3: valid p/q values, fisher 62.8% sig (inflated as expected),
nb_pairwise 21.2% (more conservative), nb_multi 100% on
marker subset (these are pre-selected differential PAS)
- Match all 3: confidence in [0,1], 12 canonical clusters as expected
for identical-input regression
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
bedtools closest output column count = (input_cols + atlas_cols + 1). We hardcode input as BED6 but atlas is variable: PolyASite 2.0 has 12+ columns (score, motif details, region type, etc.). Hardcoded _IDX_DISTANCE=12 only worked for BED6 atlas — for PolyASite 2.0 it pointed at 'NA' (motif column), not the distance integer, so the distance filter dropped every row and atlas_snap returned 0 hits. Fix: derive _IDX_DISTANCE = len(cols) - 1 per row (distance is always last). Atlas-relative indices (chrom=6, start=7, ..., strand=11) are stable since input is fixed BED6. Verified on chr22 + PolyASite 2.0 atlas (50bp distance): 269+263=532 input PAS per dataset, 147 unique atlas PAS hit, 295 mapping rows. Distance distribution: 9 at 0bp (overlap), more <50bp. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Audit table at top shows 8/8 issues addressed with file references and validation evidence. 6 additional bugs discovered during e2e testing also fixed (nb_multi filter, mnn normalization, jaccard prefix, atlas BED12, ResourceManager added, marker_selector added). Architecture correction documented: PDUI + diff moved from main ema loop to ema_switch separate command per user feedback. Performance baseline section lists wall-clock times for all stages. Profiling baseline (cProfile + memory_profiler) deferred to next branch (peak-calling-perf) where it's the entry condition. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Compares two emaout/ directories cell-by-cell to verify regression: - peakcalling/*.bed: chrom/start/end/strand identity (ignoring pasnumber) - peakcalling/*.mtx: total counts + per-column distribution match - peakcalling/*.cb.tsv: same set of CBs (order-insensitive) - unified/*.bed: same peak set - per_dataset/*/clusters.h5ad: same shape + cluster size distribution Used to validate that peak-calling-perf optimizations don't change biological output. Self-comparison test passes (sanity check). Baseline kept locally in reports/baseline_apa_completeness/ (143MB, not committed — too large for git but gitignored). Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…onfig Add a new optional RunConfig.filenames dict (YAML key "filenames", no CLI flag) that lets users rename specific pipeline output files without editing DirectoryConfig. Changes: - DirectoryConfig gains a filenames: dict field (default empty dict). Every overrideable @Property consults self.filenames.get(key, default) so the fallback is always the original hardcoded name. - RunConfig.filenames (Optional[dict], skip_legacy_bridge) is wired manually in main.py: rc.filenames -> set_directory_config(filenames=). - Overrideable keys: pasbed, posbed, negbed, filtered_cb, annotated_matrix, annotatedbed, pas_geneid, endbed, raw_features, utr_lengths, cluster_labels, clusters_h5ad, preprocessed_h5ad. - RunConfig docstring lists all overrideable keys with defaults. - yaml_loader._LIVE_KEYS picks up "filenames" automatically via yaml_keys_from_schema. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…/<ds>/ stage directories This is a BREAKING change to on-disk layout. Existing run outputs at peakatail_runs/<run>/per_dataset/<ds>/ are not readable by the new defaults — users must either re-run or pass explicit --h5ad / --pasbed paths to switch_* subcommands. Migration: every consumer routes through DirectoryConfig accessors, so flipping the layout was a single dataclass change. Stage subdirectories inside each 0X_<stage>/ are created on-demand by writers. Changes: - DirectoryConfig: per-dataset accessors now point to their respective stage dirs (01_peak_calling/<ds>/, 02_cb_filter/<ds>/, etc.). The per_dataset_dir property is removed. setup() no longer pre-creates a per_dataset/ root dir. - outputs.py: all write_* functions import directory_config and route every ds_root through directory_config.*_for(ds_id) instead of hardcoding output_dir/"per_dataset"/ds_id. - downstream_runner.py: per_dataset_dir parameter renamed to output_dir; worker calls set_directory_config(output_dir=...) at startup so all directory_config accessors work in spawned processes. Scratch files (pre_filter.mtx, filtered_matrix.mtx) live in 07_clustering/<ds>/. - main.py: single-sample _ss_h5ad and multi-sample h5ad_paths now use directory_config.clusters_h5ad_for(). Worker tuple passes output_dir instead of per_dataset_dir. Return dicts use clustering_dir key. - viz/pipeline_hooks.py: render_run_outputs() accepts clustering_dir (canonical) + per_dataset_dir (legacy alias for external callers). Dataset h5ad loaded from clustering_dir/<ds>/clusters.h5ad. - viz/run_report.py: dataset inference handles new 0X_<stage>/<ds>/ layout; dataset count uses 07_clustering/ with per_dataset/ fallback. - docs/tutorials/02-quickstart.md: Step 3 updated to show new layout with per-stage subdirectories. Step 4 paths updated accordingly. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…s that slipped in
Untracked (index-only removal, files remain on disk):
- plans/{apa-completeness,global-pool-and-cache,peak-calling-perf}.md
- TECHNICAL_REPORT_COMPLETE.md, README_TECHNICAL_REPORT.md, REPORT_DELIVERY.txt
- reports/technical_report.md, reports/README.md, reports/PDF_CONVERSION_INSTRUCTIONS.txt
- reports/build_report.sh + build_simple.py + convert_to_pdf.py + generate_*.py + run_build_simple.sh
- reports/baseline_apa_completeness/ (h5ad binary fixtures + JSON snapshots)
- run_build.sh, run_generate_figures.sh (root-level report build wrappers)
New .gitignore sections:
- AI agent artefacts: plans/, *plan*.md, *PLAN*.md, *prompt*.md, *.local.md, notes/, scratch/, analysis/
- Agent delivery notes: explicit entries for the three root-level .md/.txt files
- reports/ scratch: technical_report.md, README.md, PDF_CONVERSION_INSTRUCTIONS.txt, build helpers, baseline_apa_completeness/
- Root helpers: run_build.sh, run_generate_figures.sh
- PR scratch: PR_DESCRIPTION.md
- Dev caches: .pytest_cache/, .mypy_cache/, .ruff_cache/, htmlcov/, coverage.xml, .coverage, .DS_Store, Thumbs.db, *~, *.swp
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Phase 1 made directory_config a frozen dataclass instance, but apply_to_legacy_globals still did plain setattr on it, raising FrozenInstanceError at the first ema run (caught in the v9 layout-check smoke run). Branch on holder name: directory_config goes through the helper (which uses object.__setattr__ internally); variable_config / filter_config / args remain on plain setattr since they're not frozen. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
filter_cb, make_dataframe, preprocessing (in matrixfilter.py), and find_close (in annotate/find_close.py) used `directory_config.X` as function default arguments. Python evaluates these once at module import time — when ema.config.directory_config.output_dir is still the unconfigured "emaout" default. set_directory_config(output_dir=run_dir) runs in main() AFTER these modules are imported, so the cascading property accessors never reach the captured defaults. Result on a fresh run: filter_cb() read posmatrix.mtx/negmatrix.mtx from emaout/, found nothing, kept 1 cell out of ~1051, broke downstream clustering with `n_obs=1, n_vars=0`. Fix: defaults are now None. Each function resolves `directory_config.X` at call time inside its body so set_directory_config mutations cascade correctly. Caught by the v9 layout-check end-to-end smoke run. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
feat: PeakATail product CLI — config schema, gene-track viz, docs site, within-gene Fisher, output layout
uv pip install fails on CI runners because there's no .venv (unlike local dev where one already exists). Install into the system Python provided by actions/setup-python instead, and call mkdocs directly (no `uv run` wrapper needed once deps are on the system path). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
The header.png variant carries UCG/UG codon annotations that read as noise at README hero size. main_logo.png is the cleaner snake + "PeakAtail" wordmark and renders better on GitHub / PyPI. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
… glyph A. Single-sample path now registers a dedicated Rich progress bar for each downstream stage (CB filter, GTF annotation, PAS→gene assignment, annotated matrix, preprocessing, clustering) created right before each stage starts and advanced immediately on completion, so the user always sees the current active stage rather than a stale "Peak calling" bar. B. Sequential (non-tile) peak_calling() now accepts an optional progress_client kwarg; the bar total is set to the number of BAM chromosomes once the file is opened and advanced by 1 on each chromosome boundary (plus once for the final chromosome after the last-chrom flush). The pos-strand call per BAM receives the client; the neg-strand call gets None to avoid double-counting. Tile-mode and benchmark callers are unchanged (None default preserves backward compatibility). C. ThickBarColumn replaces Rich's default BarColumn: uses U+2588 FULL BLOCK for filled and U+2591 LIGHT SHADE for unfilled, giving a visually chunkier bar compared to the thin U+2501 heavy-horizontal-line default. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Previous U+2588 FULL BLOCK read as too chunky. Switch filled char to U+2586 LOWER THREE QUARTERS BLOCK (75 % cell height) — clearly thicker than Rich's original U+2501 heavy horizontal (~50 %) but not full cell. Also reduce the reserve for non-bar columns from 50 to 35 cells so the bar runs longer across the terminal. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
switch diff:
- Cluster-pair testing bar (total=len(pairs), >=5 pairs threshold)
advances per pair in both serial and parallel (pool.imap_unordered) paths
- progress_manager kwarg added to run_diff(); ProgressManager wraps the
call in switch_diff.py
switch length:
- PDUI compute bar (total=n_h5ads, >=5 h5ads threshold)
advances once per h5ad after all methods are done for that file
- progress_client kwarg added to run_length(); ProgressManager wraps the
call in switch_length.py
switch match:
- Marker ranking bar (total=n_datasets, >=5 datasets threshold)
- Cross-dataset matching spinner (total=None, indeterminate)
both inside a ProgressManager context in switch_match.py
switch geneview:
- Gene rendering bar (total=len(genes), >=5 genes threshold)
advances per gene (including skipped genes) in switch_geneview.py
All subcommands respect existing --no-progress flag (via common_options).
No stages added for trivial single-step calls (total=1 rule respected).
ProgressClient public API used throughout (no pm._progress direct access).
Also includes two sequenced fixes to ema run peak-calling bar:
- peackcalling.py: denominate set_total by nonempty refs x2 (pos+neg)
so both strand passes fill the same bar correctly
- main.py: pass progress_client to neg-strand peak_calling invocation
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Three coordinated changes: 1. peackcalling.py — denominate the bar by `bamfile.get_index_statistics()` nonempty refs × 2 (one peak_calling invocation per strand), instead of `len(bamfile.references)` which includes every alt/decoy/HLA contig in GRCh38 (~194 names, most carry zero reads). 2. main.py — both strand passes now share the same progress_client. The previous arrangement passed it only to the pos-strand call, so neg-strand work was invisible and the bar capped at ~50 %. 3. progress.py + main.py — add ProgressManager.finish(stage_id) and call it after both peak_calling returns. peak_calling's internal read filtering (CB tag required, MAPQ thresholds) can drop chroms below the chrom-change trigger so the advance count comes up a few short of the apriori total. finish() forces the bar to 100% so the user doesn't see a stalled bar while the rest of the pipeline runs. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…rack stage Three coupled fixes after the user ran single-pair diff and saw no bars: 1. runner.py — drop the `len(pairs) >= 5` threshold on the cluster-pair bar. Always register it; a 1/1 bar still gives uniform visual feedback that work is happening. 2. pipeline_hooks.py — auto-top-N gene_track rendering now creates its own "Gene tracks (top-N)" bar and advances per gene processed (including skips and exceptions, via try/finally). 3. switch_diff.py — widen the `with ProgressManager(...)` scope so it also covers the call to render_switch_diff_outputs. Previously the context exited right after run_diff returned, so volcano + gene_track rendering happened with no active progress manager. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Production-grade benchmark modules that run every differential APA strategy (fisher / nb_pairwise / nb_multi) and every length strategy (classic / proportion / shannon) against a clustered h5ad, emitting per-strategy summary metrics and concordance statistics. Key design decisions: - h5ad loaded once, reused across all strategies to minimise I/O - nb_multi omnibus handled gracefully (no log2fc / delta_proportion fallback) - Jaccard concordance computed only for pairwise strategies that share the same (c1, c2) namespace - length_compare builds per_gene map from adata.var without GTF; per_isoform path lazy-parses GTF only when requested and caches the isoform UTR map across strategies with the same aggregation level Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
- reports/build_pdf.sh (pandoc wrapper added today)
- reports/subset_heavy_strategies.py + explainer_figures.py (local scripts)
- reports/figures/*.svg (alongside the existing *.png rule)
- reports/{subset_,}diff_compare/ + {subset_,}length_compare/ TSV dirs
- reports/subset_curated.h5ad
- gtf_cache/ at repo root (legacy cache location)
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Adds a two-tier merger that runs after every strategy.find_pas(peak) at the
caller layer in peackcalling.py, so every strategy (current and future) gets
within-peak summit merging for free.
Why:
lambda_gradient and sierra_iterative emit "twin" PAS regions 1 bp apart
inside one peak — 371 same-strand pairs <30 bp on the v9 reference run.
These are partition_peak_region midpoint-split artefacts: the smoothed
gradient finds shallow oscillations that each clear per-summit prominence,
but the valley between them is at the local lambda — they are one peak.
No upstream knob (merge_len, pas_gap) addresses within-peak summit spacing.
Design:
Tier 1 — distance gate. Adjacent regions with gap < min_pas_spacing are
merged unconditionally. Default -1 auto-detects median read length per BAM
(98 bp for the reference SRR8325947), since summits closer than one read
length cannot be resolved by any short-read library.
Tier 2 — strategy-driven valley. Each strategy supplies its own valley
threshold via a new PeakFinderStrategy.valley_threshold(peak, user_setting)
method. lambda_poisson and lambda_gradient override to return their own
compute_lambda(heights) — fully dynamic, in lock-step with the strategy's
own background estimate. original and sierra_iterative use the static
min_pas_prominence (default 5.0). Merge when valley < threshold.
Reads inside the merge gap are absorbed automatically: peak.cb_positions
is dense across the parent peak's full span, and reconstruct_cb_dict
slices by inclusive range — the widened survivor region picks up every
read end position between the two pre-merge boundaries without any
explicit summing logic.
Wiring (all three live paths):
- peackcalling.py (monolithic): 5 find_pas call sites + signature params
+ auto-detect resolution after BAM open + forwarding to run_tiled and
run_pipeline dispatchers.
- peak_pipeline.py (--pipeline): merger inside writer_loop and
_flush_peak; run_pipeline signature + writer subprocess args; sentinel
resolved before spawn.
- tile_runner.py (--tiles): JobSpec carries both params; run_job
handles JobSpec + dict modes; pool-wide sentinel resolved once before
spawning workers.
Configuration:
--min-pas-spacing N Tier 1 distance (bp). -1 = auto-detect, 0 = off.
--min-pas-prominence X Tier 2 static fallback. Negative = Tier 2 off.
Both surface via YAML keys (min_pas_spacing, min_pas_prominence) bridged
through RunConfig -> variable_config (legacy_dataclass_attr pattern).
Defaults give the recommended on-on behaviour; setting both to off-
values reverts to bit-exact pre-merger output for reproducibility.
Effect on v9 lambda_gradient run (SRR8325947):
Total PAS 15,948 -> 15,568 (-380)
Close pairs <30bp 371 -> 0 (-371)
Filtered cells 752 -> 752 (no change)
Leiden clusters 11 -> 11 (no change)
Classic length 1,925 -> 1,898 genes (-27 lost as single-PAS)
Proportion length 5,339 -> 5,407 genes (+68, cleaner gene assignment)
Shannon entropy 5,339 -> 5,407 genes (+68, same reason)
Validation:
- 13 unit tests in tests/test_pas_merge.py: both tiers, lambda override,
distance floor, chain merges, empty/single/unsorted inputs.
- All 4 strategies on full SRR8325947 BAM with merger ON: 0 close
pairs <30bp across the board.
- Pipeline == monolithic equivalence on chr22 lambda_gradient: bit-
identical 252-PAS output.
- 665 of the 668 existing tests still pass; 3 pre-existing failures
confirmed unrelated by stash test.
Docs + report:
- New docs/strategies/pas-merger.md page (added to mkdocs nav).
- docs/cli/run.md updated with --min-pas-spacing and --min-pas-prominence.
- reports/PeakATail_Technical_Report.tex: new Section 5 "Post-Detection
PAS-Summit Merger" + pas_merge_explainer.png + Table 6 stats updated to
the new run (1,898 / 5,407 / 5,407 / 20 DPYD PAS).
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
The proportion_heatmap had no genomic-position axis: rows were just ``pas_id`` integers sorted by cross-cluster variance. Two PAS from the same gene 12 kb apart could appear as adjacent rows purely because they both made the top-50 by variance, giving a false impression of genomic proximity. Two changes: 1. Row order: group by gene_id (preserving variance-rank for the first PAS of each gene), then sort by genomic start within each gene. Same-gene PAS now sit together AND in genomic order. 2. Y-axis labels: replace bare ``pas_id`` with ``pas_id gene_id chrN:start-end (strand)`` so each row tells the reader exactly where on the genome the PAS lives. Thin white separator lines drawn between gene groups for quick visual distinction of same-gene vs different-gene rows. Falls back to legacy label format when position metadata is missing from ``pdui_df`` (so old test fixtures keep working). Updates figure metadata: adds ``n_genes`` field, ``row_order`` and ``y_label_format`` strings so downstream consumers (figure index, report captions) know how rows are organised. Regenerated reports/figures/length_proportion_heatmap.png from the lg_merge_2026-05-15 v9 run; technical report PDF rebuilt. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…tions Two visual misleads in the gene-track view: 1. ``pas_positions`` was the midpoint of each PAS region; the actual start/end span was discarded in ``build_gene_panel`` and never reached the renderer. A 500-bp merged PAS rendered identically to a 50-bp singleton. 2. Bar width in the coverage rows was a fixed ``1% × gene_span`` placeholder (matplotlib) or plotly's auto width — neither reflected the PAS region's true extent, so on long genes PAS appeared to crowd together at midpoints without any cue to their real width. Changes: * ``GenePanel`` gains ``pas_starts`` / ``pas_ends`` (same length as ``pas_positions``). ``build_gene_panel`` now copies these straight from the pasbed coords. Defaults to empty lists so older callers that construct ``GenePanel`` directly keep working. * Matplotlib backend renders each coverage bar with ``align="edge"`` at the actual PAS span, clamped to a minimum visible width of ``0.5% × gene_span`` so 1-bp PAS don't disappear on long genes. * New ``_annotate_pas_positions`` helper draws a translucent red stripe + ``pas_id`` label on the gene-structure track at each PAS's true genomic span — readers can immediately see which exon / 3'UTR the PAS occupies and whether two visually-adjacent bars are genuinely close on the genome. * Plotly backend passes the per-bar widths array to ``go.Bar.width`` so the same width information is preserved in the interactive viz. Regenerated reports/figures/geneview_with_isoforms_DPYD.png from the lg_merge_2026-05-15 v9 run; technical report PDF rebuilt. 665 of 668 existing tests still pass (the 3 pre-existing failures from develop are unrelated). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Two correctness issues spotted on the DPYD geneview:
(1) The x-axis range was derived from pasbed coordinates only
(``coords["start"].min()`` / ``coords["end"].max()`` in
``build_gene_panel``). For DPYD on chr1 the pasbed-detected range was
97,169,051 – 97,882,620 (713 kb) but the GTF gene actually spans
97,077,743 – 97,995,000 (917 kb). That's a ~200 kb (22%) chunk of
exon structure clipped off-screen on every render.
Now the renderer also walks ``panel.isoforms`` and expands x_lo/x_hi
to include the GTF-supplied exon backbone, so the full gene structure
is always visible. PAS positions still come from pasbed (correct);
the GTF only contributes the structural backbone.
(2) The x-axis label said "Genomic position (bp)" but tick values for a
chr-1 gene at ~97 Mb were rendered by matplotlib in scientific
notation (``9.7e+07``) — confusing, and the "(bp)" claim was
technically true but easily mistaken for kb or Mb.
Now a FuncFormatter selects Mb / kb / bp based on the actual gene
span (>= 1 Mb -> Mb with 2 decimals; >= 1 kb -> kb with 1 decimal;
< 1 kb -> bp with commas). Axis label is rewritten to
"Genomic position on chrN (Mb)" so the chromosome and unit are both
explicit.
Regenerated reports/figures/geneview_with_isoforms_DPYD.png from the
lg_merge_2026-05-15 v9 run; technical report PDF rebuilt.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
…ility
Three concrete issues spotted on the CLIC2 (3-PAS) walkthrough:
1. The x-axis showed absolute genomic coordinates ("155,276 kb" on
chrX) which look chromosome-sized even though we're zoomed into a
58 kb gene. Switched to "distance from gene start" with the
absolute anchor in the axis label
(``Distance from gene start (kb) -- anchor chrX:155,276,210``)
so the 0-anchored axis makes the gene-scope visually obvious.
Scale heuristic refined to three tiers:
< 5 kb span -> bp
5 kb - 1 Mb -> kb
>= 1 Mb -> Mb
so a 3 kb peak window shows ``0, 500, 1000 bp`` ticks rather
than ``0.0, 0.5, 1.0 kb``.
2. Title showed only the Ensembl gene ID (``Gene ENSG00000155962``).
When the GTF is supplied, the renderer now pulls
``gene_name`` ("CLIC2") via a new ``load_gene_name_from_gtf``
helper and the title leads with the symbol:
``Gene CLIC2 (ENSG00000155962) -- chrX:155,276,211-155,334,657
(- strand, 58.4 kb) -- 3 PAS``.
``gene_name`` is stored on ``GenePanel`` and exposed in the
figure meta sidecar JSON.
3. Real GTF exons (80-200 bp on a 50-900 kb gene = 0.1-0.4 % of
axis width) rendered as 1-2 pixel slivers, easy to mistake for a
rendering bug. Added a minimum visible exon width clamp at
0.4 % of the view span -- the GTF coordinate is unaltered, only
the rectangle width is widened around its true left edge so tiny
exons stay legible. PAS overlay alpha lowered from 0.18 to 0.10
so the red tint doesn't wash out small exons.
Side benefits:
* x_lo / x_hi now take the UNION of pasbed coords and the GTF exon
backbone, so a gene whose 5' end has no detected PAS still
renders the full structure (DPYD: previously clipped 200 kb /
22 %% of the gene).
* meta sidecar gains ``gene_name``, ``view_start``, ``view_end``,
``view_unit`` so downstream consumers know the rendered window
and unit without re-deriving them.
Regenerated reports/figures/geneview_CLIC2_3pas.png +
geneview_with_isoforms_DPYD.png from the lg_merge_2026-05-15 v9
run; technical report PDF rebuilt.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
feat(merger): strategy-agnostic post-detection PAS-summit merger
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Report Contents
Tools Analyzed
scMAPA, SCINPAS, SCAPTURE, scTail, scDaPars, Sierra, scAPAtrap, SCAPE, scLAPA, spvAPA, MAAPER, scraps, scAPA, scDAPA, scPAISO, InPACT, SAPAS, vizAPA, stAPAminer
Key Finding
Test plan
🤖 Generated with Claude Code