bedfile plots with gtars - #127
Conversation
Release 0.14.10
Release 0.14.11
Release 0.14.12
Release 0.15.0
There was a problem hiding this comment.
🟡 Changes recommended
There are confirmed correctness/documentation issues in the new aggregation output (partition scaling and per-chrom n handling) that should be addressed before merge.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Review details
- Files reviewed: 5/5 changed files
- Comments generated: 7
- Review effort level: Lite
| for row in rows: | ||
| entry = result.setdefault(row.chrom, {"mean": [], "sd": [], "n": int(row.n)}) | ||
| while len(entry["mean"]) <= row.bin_idx: | ||
| entry["mean"].append(0.0) | ||
| entry["sd"].append(0.0) | ||
| entry["mean"][row.bin_idx] = float(row.mean) | ||
| entry["sd"][row.bin_idx] = float(row.sd) |
There was a problem hiding this comment.
n counts how many member files contributed a value at a given bin. The arrays are dense from index 0 — jsonb_array_elements_text ... WITH ORDINALITY yields bin_idx = ordinality - 1 — so every file with data for a chromosome appears in bin 0. That makes bin 0 the bin with the largest count, and that count is the number of files with data for that chromosome.
The query orders by bin_idx, so the first row for each chromosome is bin 0. The stored n is that count, not an arbitrary bin.
Agreed it is implicit, and it would break if the ORDER BY were ever dropped. Leaving it as-is for now since nothing reads the field — bedbase-ui passes the top-level n_files into the region-distribution plot.
| # Fallback: wrap old scalar columns. | ||
| return BedSetDistributions( | ||
| n_files=0, | ||
| scalar_summaries=_old_stats_to_scalar_summaries( | ||
| bedset_object.bedset_means, | ||
| bedset_object.bedset_standard_deviation, | ||
| ), | ||
| ) |
There was a problem hiding this comment.
Fixed in 92fe156. get_distributions now passes bedset_object.bedfile_count through as n_files and into each scalar's n.
The fallback also omitted histogram entirely while the aggregation path always emits it, so a consumer could only tell the two shapes apart by probing for absence — bedbase-ui destructured it and threw. It now emits histogram: None, and the UI guards on it.
| def get_batch( | ||
| self, | ||
| identifiers: list, | ||
| full: bool = False, | ||
| distributions: bool = False, | ||
| ) -> BedBatchResult: |
There was a problem hiding this comment.
Agreed, and not addressed here. The blocker is the fixtures rather than the test itself: tests/utils.py inserts a single bed with no distributions blob, so n=1 collapses the statistics — stddev is null, and col_min == col_max sends _scalar_histogram down its single-bin early return — and the defer(BedStats.distributions) path has nothing to defer.
A meaningful test needs a second bed with different values and a distributions blob on at least one. That is a bigger change than belongs in this PR; worth doing separately.
| distributions: bool = False, | ||
| ) -> BedBatchResult: | ||
| """ | ||
| Get multiple bed file records by identifiers in a single DB round-trip. |
There was a problem hiding this comment.
Taken as suggested in 5304f56.
On the count, for the record: SelectInLoader chunks parent keys at 500 (sqlalchemy/orm/strategies.py:2967), so get_batch costs three statements up to 500 identifiers and five at the MAX_BATCH_SIZE ceiling of 1000. The new wording avoids committing to a number, which is the right call.
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
5edc86f removed the `agg_columns.extend([` wrapper along with the `* 100` multipliers, leaving an orphaned list literal. aggregation.py has not parsed since, and both bedsets.py and bedfiles.py import it, so `import bbconf` fails outright — conftest.py cannot load and no tests are collected. Restore the wrapper. The multipliers stay off, but for a different reason than the autofix assumed: the bed_stats.*_percentage columns hold a fraction, not a percentage. Both producers divide by the region count — regionstat.R:212 stores Freq/length(query), gtars_backend.py:227 stores count/total — and bedbase-ui renders the per-file values as value * 100. Since the aggregation no longer rescales, `mean_pct`/`sd_pct` would have been misnamed, so they become `mean`/`sd`. That also matches the shape already used by `scalar_summaries`, leaving every value in BedSetDistributions on the same scale as the column it aggregates. bedbase-ui scales to percent at the plot layer, where the local comparison path already does. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
The scalar was mapped to bed_stats.tssdist, a leftover TSS column that nothing writes: BedStatsModel has never declared the field, so the value bedboss computed was discarded by extra="ignore" before insert. count() returned 0, `if not n: continue` fired, and the key was silently absent from every scalar_summaries payload since the aggregation was written. Rather than adding a column, a model field and a migration to wire it up, drop it. Remove the aggregation entry, the matching key in the _old_stats_to_scalar_summaries fallback, and the BedSetDistributions docstring line that pointed readers at it as the replacement for the neighbor_distances KDE. That line now states the real reason the KDE is not aggregated: each file's curve is fit over its own log10 range, so bin i means a different bp value per file. bed_stats.tssdist is left in place; removing it needs its own migration. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
Four fixes, none behavioural beyond the values reported: region_distribution set each chromosome's "n" from whichever (chrom, bin) row arrived first and never updated it, so a chromosome whose member files have ragged arrays reported the count for an arbitrary bin. Take the max across bins, which is the number of files contributing anything to that chromosome. get_distributions' pre-aggregation fallback reported n_files=0 and n=0 for bedsets that do have members. The bedset row carries bedfile_count, so use it. The fallback also omitted "histogram" entirely while the aggregation path always emits it — two shapes for one field, which a consumer can only distinguish by probing for absence. Emit histogram: None instead; there are no per-file values to bin, and the key now matches. The BedSetDistributions docstring described tss_histogram as summed per-bin counts; _aggregate_tss_histogram computes AVG and STDDEV. get_batch's docstring claimed a single DB round-trip; the two selectinload options make it three. It is batched, not single. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
Same defect Copilot flagged in region_distribution, one function over and unreported: "n" was read from rows[0], so it described whichever bin sorted first rather than the collection. Member files emitting different-length counts arrays give bins different contributing-file counts, and the reported value was then arbitrary. Take the max, matching what _aggregate_region_distribution now does. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
Both queries order by bin_idx, and every file with a non-empty array for a chromosome contributes bin 0, so the first row already carried the number of contributing files — max() computed the same value. The change was robustness against a future reordering, not a fix, and nothing reads either field: bedset-plots.ts passes the top-level n_files into regionDistributionSlot. Reverts 3e608c9 and the aggregation.py hunk of 92fe156. The fallback n_files, histogram shape, and the two docstring corrections from 92fe156 stand. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
"with batched DB queries" is tighter than "without a per-record query" and matches the review comment's own phrasing. Neither commits to a query count, which is the right call: selectinload chunks parent keys at 500, so get_batch costs three statements up to 500 identifiers and five at the MAX_BATCH_SIZE ceiling of 1000. The `or 0` on bedfile_count was dead: the column is nullable=False in the initial migration, so the guard implied a NULL that cannot occur. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PRGhWjGiUXFoygd35QsKQ7
Changes:
_aggregate_region_distribution() (aggregation.py:245) takes all member files of a bedset and, per (chromosome, bin index), computes the mean, sd, and n across files:
Concerns:
The scalar stats are one number per file (a plain AVG over N rows).
region_distributionis a nested array per file, so aggregating means materializing files × chromosomes × bins rows in SQL before the GROUP BY — that product is what can blow up for 10K files, depending on how many bins gtars emits per file.TODO:
__version__.pyfile