Migrate mutation data to SQLite (gambl_mutations.db), fix SLMS-3/Publication dedup bug, auto-download via GitHub Releases - #128
Open
rdmorin wants to merge 60 commits into
Open
Migrate mutation data to SQLite (gambl_mutations.db), fix SLMS-3/Publication dedup bug, auto-download via GitHub Releases#128rdmorin wants to merge 60 commits into
rdmorin wants to merge 60 commits into
Conversation
Add full table/column/join-key reference to ?gambl_reference_db and ?gambl_mutations_db (rendered as markdown tables), plus a short "Bundled databases" pointer in the README, so the schema isn't only in the build scripts. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
pull_data(), the publication-samples blocks, and the trios blocks pulled every mutation for a sample set via get_ssm_by_samples() only to filter down to all_lymphoma_genes afterward. Replace with get_ssm_by_regions() against a shared lymphoma_genes_bed_grch37/hg38 (built once via gene_to_region()), so tabix -R restricts to the gene loci up front instead of discarding everything outside them after the fact. cell_lines_data and coding_maf were left on get_ssm_by_samples since they don't filter by gene. Also wrap every get_ssm_by_samples()/get_ssm_by_regions() call site (15 total, including the pre-existing aSHM pulls) in a time_it() helper that logs wall-clock time, to compare against get_ssm_by_samples(subset_from_merge = TRUE/FALSE) on real GSC data. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The Manta SV block (feeding sample_data$grch37/hg38$bedpe) filtered against GAMBLR.data::sample_data$meta$sample_id -- the currently INSTALLED package's bundled metadata, not the metadata this run is actually assembling. That made bedpe row counts depend on whatever GAMBLR.data happened to be installed at build time rather than purely the input data (e.g. 941 vs 2352 rows for the same underlying SV calls, traced to the installed sample set shifting between builds -- confirmed via nrow(GAMBLR.data::sample_data$meta) vs nrow(bundled_meta) diverging). Move the block to run after sample_data$meta is finalized and filter against that local object instead, so bedpe scope is self-consistent with the rest of the bundle for a given run. maf/ashm/seg were never affected -- only this SV block referenced the installed package. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Compares bedpe between the bundled sample_data.rda and a fresh gambl_mutations.db build by sample set, not just row totals: unique sample overlap, and per-sample row-count delta for samples present in both. Used to check whether a bedpe count difference (e.g. after the SV-scoping fix in assemble_bundled_data.R) comes from more samples having SVs at all, or the same samples yielding more SVs than before. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Rename compare_bedpe.R -> compare_bundle_changes.R and extend it to cover
all three data types (SNV: maf+ashm, CNV: seg, SV: bedpe), each compared at
the level of individual events (sample + genome_build + coordinates), not
just row totals:
- per-sample_id counts of gained (new) and lost (missing) events
- one example gained-event row per affected sample, written to
{snv,cnv,sv}_gained_examples.log for quick manual inspection
- full per-sample gained/lost counts written to
{snv,cnv,sv}_persample_counts.tsv
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
write_examples() was already generic (works on any event data frame, not
just "gained"); call it on cmp$lost too, per type, writing
{snv,cnv,sv}_lost_examples.log alongside the existing gained-examples logs
and per-sample counts. One representative lost row per affected sample, so
e.g. samples that lost all their SNVs can be inspected directly instead of
just counted.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
get_gambl_metadata()'s coverage-based QC filter (min_corrected_cov, default 15x) can silently drop samples that were already part of a released bundle out of a fresh rebuild -- e.g. 01-16433_tumorA/B, FFPE genomes with coverage under the bar, taking all their SNVs with them in compare_bundle_changes.R's diff. Shadow get_gambl_metadata() near the top of the script with min_corrected_cov=0 baked in, so every one of the ~13 call sites picks this up in one place: assembly should not lose samples already in the released dataset due to a QC threshold applied after the fact. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
get_ssm_by_regions() with lymphoma_genes_bed restricted by genomic
coordinates alone, which don't always agree with VEP's Hugo_Symbol
assignment at gene boundaries: mutations in overlapping/neighbouring genes
and non-coding loci (AC/AL/AF-prefixed lncRNA transcripts, etc.) were
leaking in as false "gains", confirmed via snv_gained_examples.log showing
136 genes that were never in all_lymphoma_genes.
Fix has two parts, applied to all 7 gene-restricted call sites (pull_data,
both publication-samples blocks, all 4 trios blocks):
- filter(Hugo_Symbol %in% all_lymphoma_genes) restored after the region
pull, on the now-small result (cheap, unlike filtering pre-restriction)
- gene_to_region() calls building lymphoma_genes_bed_grch37/hg38 now pass
pad_length = GENE_PAD_BP (2000bp) so real target-gene mutations just
outside the exact gene span aren't excluded by tabix before the filter
even sees them
aSHM-region call sites are unaffected -- they are legitimately region-
defined regardless of Hugo_Symbol and don't get this filter.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
GENE_PAD_BP=2000 was still narrower than VEP's default 5kb upstream/downstream annotation window, so real Hugo_Symbol-labeled mutations just outside the old pad (e.g. an ID3 variant 4.2kb upstream of its TSS) were dropped by the tabix -R pull before the post-filter ever saw them. 10kb gives a 2x margin; the existing Hugo_Symbol filter keeps over-padding safe. compare_bundle_changes.R now also reports n_retained per sample so a sample's total-loss vs partial-loss status is visible directly in snv/cnv/sv_persample_counts.tsv, instead of needing to be inferred.
Companion to GAMBLR.utils's gene_to_region()/expand_gene_aliases() fix. all_lymphoma_genes feeds 7 filter(Hugo_Symbol %in% all_lymphoma_genes) post-filters in this script; without expansion those still drop rows annotated under a gene's other name (e.g. old/new HGNC histone names) even though the region pull now finds the right coordinates.
maf and ashm are assembled by different code paths in assemble_bundled_data.R (gene-region-restricted pull vs aSHM-region-restricted pull), but the SNV gained/lost report merged both into one undifferentiated bucket -- there was no way to tell which code path produced a given lost/gained example without manual digging. .source (not part of the comparison key) makes this visible directly in snv_persample_counts.tsv / snv_*_examples.log.
Confirmed against every DB-querying function in GAMBLR.open (get_ssm_by_region(s), get_ssm_by_samples, get_coding_ssm, get_all_coding_ssm) via the shared GAMBLR.data::get_ssm_from_db() helper: gene-restricted queries resolve the gene to a region first (same logic used to populate these tables) and filter on (genome_build, Chromosome, Start_Position), never Hugo_Symbol directly. The index was pure overhead -- ~55MB on disk for a column no code path queries.
…MS-3 pulls Root cause of a real, confirmed bug: assemble_bundled_data.R pulled SLMS-3 mutation calls via separate, cohort-specific code blocks, each independently deciding which samples to pull for. A sample belonging to more than one cohort block (e.g. the Dreval FL x Hilton trios overlap) got its mutations pulled and bound multiple times -- verified directly against the previously-shipped gambl_mutations.db: 1,797 exact duplicate maf rows (same pattern in ashm), plus true duplicate rows in sample_meta for the same 5 samples. There was no deduplication safety net anywhere in the pipeline. This surfaced while replacing the Arthur cohort's raw, completely unfiltered strelka flat-file dump (2.83M rows / 65,121 distinct genes, missing the lymphoma-gene-panel filter every other block applies) with a proper SLMS-3 pull. Root cause of both problems is the same: Study/cohort membership was used to gate whether mutations get pulled, instead of being tracked as independent sample metadata. - assemble_bundled_data.R: restructured into phases -- metadata (all cohorts, including Arthur and Hilton, previously added so late that Arthur got zero aSHM coverage and Hilton needed its own dedicated aSHM block to compensate) assembled once; one consolidated SLMS-3 pull for the full deduplicated sample universe instead of ~5 separate per-cohort pulls; Publication-pipeline pulls (genuinely distinct per-paper data, not a redundant recall) kept separate. Arthur's sample_id/Tumor_Sample_Barcode are no longer overwritten to its paper's own patient-level ID (the paper's "Case ID" is now captured as study_id instead); the all-capture cohorts' patient_id was being silently dropped before reaching sample_data$meta, fixed in the same pass since the exact code region was already being touched. - write_mutations_db.R: writes the new sample_study table (sample_id, study, study_id, reference_PMID -- a many-to-many bridge table, one row per (sample_id, study) pair, so multi-study samples need no schema change); adds a write-time defensive dedup on maf/ashm keyed on the same natural mutation-call key GAMBLR.open::get_ssm_by_region() already uses at read time, as a permanent safety net against any future re-introduction of overlapping cohort blocks; drops sample_meta.reference_PMID (single-valued, already demonstrably wrong for multi-study samples) now that sample_study is the source of truth. - get_ssm_from_db.R: this_study now resolves via a join against sample_study instead of a Study column on maf/ashm (which no longer exists), reusing the existing sample_ids/Tumor_Sample_Barcode filtering machinery rather than adding new plumbing. Verified against a synthetic DB: this_study alone, this_study+tool_name, this_study+sample_ids intersection, and an unmatched study (returns zero rows, not an error) all behave correctly. - gambl_mutations_db.R: schema docs updated for the new table and the Study column removal. - test_gambl_db.R: checks updated for the schema change; also fixes a pre-existing bug (asserted idx_maf_gene, which has never existed in write_mutations_db.R -- this check failed regardless of this refactor). Both the write-time dedup and the sample_study many-to-many mechanics were verified directly against synthetic data before this commit (see session record) since the full pipeline requires GSC access to run.
Reddy's raw "Sample ID" column is read as numeric by read_excel(), while every other cohort's study_id is character -- bind_rows() can't reconcile <chr> with <dbl>, confirmed by a real run on the GSC: "Can't combine ..1\$study_id <character> and ..4\$study_id <double>". Coerced study_id to character at every cohort's sample_study frame construction (not just Reddy's), since any of these xlsx-sourced columns could hit the same type-inference issue depending on how the source spreadsheet happens to type that column.
Tracks test_gambl_db.R results and the reasoning behind each of the 4 expected regression-check differences vs the old bundle, for use writing the eventual PR description. Has a running checklist to append further verification (compare_bundle_changes.R, Dreval/Hilton dedup spot-check, Arthur gene-panel restriction) as it's done.
…l/Hilton)
patient_id is patient-level, not sample-level -- ambiguous for any
cohort where a patient can have more than one sample. Confirmed as a
real bug for Hilton specifically (a trios study): LY_RELY_116_tumorA
and LY_RELY_116_tumorB both collapsed to study_id="LY_RELY_116",
making the two sample_study rows indistinguishable. Thomas BL/DLBCL
and Dreval don't currently hit this in practice but share the same
underlying imprecision, fixed for consistency.
All four now use sample_id, already sourced from each xlsx's own
"Genome sample id"/DNAseq_sample_id column -- the closest thing to a
study-native sample-level identifier available for these cohorts.
Arthur (structurally protected by its !grepl("tumor", sample_id)
filter) and Reddy (study_id sourced from the paper's own distinct
"Sample ID" column, already 1:1 with sample_id) were left unchanged.
compare_bundle_changes.R showed patients losing 100% of their SNVs
with zero retention (e.g. patient 08-15460: 19,897 lost, 0 retained).
Root-caused to two compounding bugs in the Arthur block:
1. arthur_meta matched Case ID (from DLBCL_Arthur.xlsx) against
GAMBL's patient_id via `patient_id %in% arthur_case_ids$"Case ID"`.
Same bug class as the study_id chr/dbl mismatch fixed earlier --
read_xlsx() can silently type a leading-zero ID like "08-15460" as
numeric, and %in% fails to match across types without error. Fixed
with explicit character coercion + inner_join() instead of a %in%
filter.
2. `! grepl("tumor", sample_id)` excluded every sample for any patient
with more than one tumor biopsy (..._tumorA/..._tumorB -- the same
multi-sample-per-patient pattern Hilton has). For a patient whose
only GAMBL samples are tumor-suffixed, this matched nothing,
dropping their SLMS-3 coverage entirely. No other cohort in this
script excludes samples this way; removed.
…code
Root-caused the 100%-loss patients (e.g. 08-15460: 19,897 lost, 0
retained) to the now-deleted Arthur raw flat file, which used bare
patient-style Tumor_Sample_Barcode values independent of arthur_meta
entirely -- confirmed by reading the old file directly
(grep("15460", ...) on Tumor_Sample_Barcode found "08-15460" in the
raw file itself). Since that file is already deleted, this specific
loss was never a bug -- it's the expected disappearance of a
non-standard ID convention, correctly superseded by GAMBL's real
sample IDs (confirmed present and correctly sized).
Generalized into a defensive fix for the remaining Publication-pipeline
pulls, which read a paper's own file directly and were trusting its
Tumor_Sample_Barcode values without validation. relabel_to_sample_id()
joins on sample_study$study_id and rewrites Tumor_Sample_Barcode to
the real sample_id wherever they differ -- a no-op for Thomas/Dreval
(study_id already mirrors sample_id there), an actual fix for Reddy's
original-variants file (study_id is the paper's own distinct raw
"Sample ID"). Verified against synthetic data.
Pipeline=="strelka" was Arthur's old, now-deleted unfiltered raw-file dump -- confirmed twice (08-15460, DO52686) that these rows' Tumor_Sample_Barcode values are the raw file's own non-standard IDs, not GAMBL's real sample_id, so every one is expected to show up as "lost" regardless of anything else. Filtered out of the old baseline before the comparison so this already-understood, intentional disappearance stops drowning out genuinely new findings on every run.
compare_bundle_changes.R confirmed the scope reduction flagged during planning: folding cell lines into the panel-restricted consolidated SLMS-3 pull dropped their genome-wide coverage down to ~150-237 genes, showing up as 43k-65k "lost" rows per cell line vs the old bundle. Cell lines (DOHH-2, SU-DHL-10, OCI-Ly10, OCI-Ly3, SU-DHL-4) are now excluded from all_slms3_meta and get their own dedicated, unrestricted get_ssm_by_samples() pull instead, both builds -- matching their original pre-refactor behavior, bound into slms3_grch37/slms3_hg38 alongside everyone else's panel-restricted data.
… calls GAMBLR_examples_output.log surfaced a real bug: get_coding_ssm() for 183 Dreval FL patients returned zero rows, when it used to return real coding-classified data. Traced via elimination (ID matching, t_alt_count type, min_read_support all confirmed fine) to coding_only=TRUE alone zeroing out 28,867 real SLMS-3/grch37 rows -- every one of them was non-coding, which is biologically implausible across 183 real patients and the whole lymphoma gene panel. Root cause: unix_group (and likely other native GAMBL metadata columns) is NA for every cohort built by hand in this script (BL_Thomas, DLBCL_Thomas, FL_Dreval, DLBCL_cell_lines -- 878 samples) rather than pulled live from get_gambl_metadata(), and get_ssm_by_regions() apparently needs it to correctly resolve which merged flat file to read for a sample. This was specific to the two things newly introduced this session -- Phase 4's consolidated SLMS-3 pull and Phase 5's consolidated aSHM pull -- both of which passed sample_data$meta-derived subsets directly to get_ssm_by_regions(). Every other GAMBLR.results call site in this script already re-fetched fresh, complete metadata via get_gambl_metadata() filtered by sample_id first, which is why this only surfaced here. Fixed both: sample_data$meta is now used only to determine which sample_ids belong in each pull; the actual these_samples_metadata passed to get_ssm_by_regions() is always a fresh get_gambl_metadata() call filtered to that same sample_id set.
…, not study GAMBLR.data::sample_data$meta (the installed, stale snapshot this block intentionally reads) has FL_Dreval under the study column, not cohort -- filtering this_study_samples on cohort silently collapsed it to just the 5 cell-line sample_ids, so coding_maf (meant to enrich Dreval's Publication rows with RefSeq/Protein_position) was built from cell-line mutations instead and never matched any Dreval row. Same cohort/study mismatch fixed for the BL_Thomas/DLBCL_Thomas/DLBCL_cell_lines lookups above it. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Prints per-study row/sample counts (and, at the two live-metadata refetch points, exactly which sample_ids get silently dropped) at every major transition: per-cohort metadata assembly, the cohort->study rename/rejoin block, sample_study assembly, the Phase 3 proteinpainter-enrichment joins, the Phase 4 all_slms3_meta refetch and pull outputs, the Phase 5 aSHM pull, and the final per-study coding-classified maf counts. Goal: see exactly where a cohort's count changes during a build instead of only being able to diagnose it after the fact via a separate GSC session against the finished DB. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ode column get_gambl_metadata() already returns a Tumor_Sample_Barcode column alongside sample_id, so renaming sample_id -> Tumor_Sample_Barcode for the all_slms3_meta and ashm_pull_meta diagnostics collided and crashed the build (rename(): "Names must be unique"). Drop the existing column first. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Two DB performance issues, likely contributing to the lag observed while rerunning GAMBLR.open's tests: 1. idx_maf_pipe/idx_vp_pipe were plain b-tree indexes on the raw Pipeline column, but every actual query filters on tolower(Pipeline)/lower(Pipeline) for case-insensitive matching -- SQLite can't use a plain index to satisfy a function-wrapped predicate, so these indexes were never actually usable and every tool_name-filtered query fell back to a full table scan. Fixed at the source instead of with expression indexes: Pipeline is now normalized to lowercase once, at write time, in write_mutations_db() (both maf/ashm and variant_pipeline), so reads can use a plain equality against a plain index. get_ssm_from_db() still accepts any-case tool_name from callers (lowercased before comparing), so this is invisible to existing callers. 2. The new variant_pipeline semi-join (tool_name resolution) had no index covering its join key on either side, meaning that join -- now on the path of every default-tool_name query -- had to scan maf/ashm in full. Added idx_maf_variant_key, idx_ashm_variant_key, idx_vp_variant_key on (Tumor_Sample_Barcode, Chromosome, Start_Position, End_Position, Tumor_Seq_Allele2). Verified with a synthetic build: Pipeline values come out lowercase in both tables, the LOWER()-expression indexes are gone (plain indexes only), and tool_name resolution (both exact-case and caller-supplied mixed-case) still works correctly through the variant_pipeline join. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
get_ssm_by_regions() commonly passes one region per gene in a panel (150-250+ for the full lymphoma gene list). The previous implementation looped and issued one full query per region, each one re-deriving the tool_name/coding_only/sample_ids filters and re-running the variant_pipeline semi-join from scratch -- 150-250x more query overhead than necessary for a single call, and on the hot path of every panel-restricted SLMS-3 pull. Regions are now combined into a single OR'd filter condition and applied in one query per table (maf/ashm), matching the same logical semantics (each region an independent, OR-ed match) but with one round trip instead of N. This also fixes a latent duplicate-row bug: a variant landing inside two overlapping regions (e.g. neighbouring genes' padded windows) was previously returned once per matching region and bind_rows()'d into duplicate rows; a single query's OR condition doesn't double-count a row just because more than one disjunct matches it. Verified with synthetic data: exact expected rows returned across overlapping regions with no duplicates, a non-matching region correctly returns zero rows, and omitting regions entirely is unchanged. Adds rlang to Imports (used directly now via rlang::expr(); previously only a transitive dependency via dplyr). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
bedpe previously had only one index (tumour_sample_id alone). Every get_sv_from_db()/get_manta_sv() call also filters genome_build, and its optional region filter checks both breakpoint ends (CHROM_A/START_A OR CHROM_B/START_B) -- none of which had anything to use, meaning a full table scan on every call regardless of how selective sample_ids/region actually were. Extended idx_bedpe_sample to (tumour_sample_id, genome_build) -- keeping the more selective sample filter leading -- and added genome_build-led indexes on each breakpoint end so SQLite's OR-optimization can use a different index per side of a region search. Verified with synthetic data: all 3 indexes present, sample_ids and region filters both still return correct results. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…tion_id Per review feedback: instead of variant_pipeline repeating the full 5-column natural key (Tumor_Sample_Barcode, Chromosome, Start_Position, End_Position, Tumor_Seq_Allele2, genome_build) on every row, maf/ashm now get a surrogate integer mutation_id assigned at write time, right after deduplication (unique within an elem, across both genome builds). variant_pipeline shrinks to just (mutation_id, elem, Pipeline) -- a plain foreign key instead of 8 columns. The natural-key match doesn't disappear -- write_mutations_db() still has to correlate the pre-dedup, Pipeline-tagged rows against the now-ID-assigned deduped rows, via the same 5-column key, to figure out which mutation_id each pipeline claim belongs to. But that now happens once, at write time (a one-time cost during the already-slow batch build), instead of inside get_ssm_from_db()'s semi-join on every single read query. The read-time join is now a single indexed integer column instead of a 5-column composite (two of them text), which SQLite can resolve far more cheaply, and sidesteps SQLite's lack of native tuple-IN matching entirely rather than working around it. Replaced idx_maf_variant_key/idx_ashm_variant_key/idx_vp_variant_key/ idx_vp_sample with idx_maf_mutation_id/idx_ashm_mutation_id/ idx_vp_mutation_id; idx_vp_pipe drops genome_build (redundant now -- mutation_id already implies build via the maf/ashm row it references). Verified with synthetic data: mutation_id stays unique and correctly isolated across genome builds even when two builds share identical coordinate values; variant_pipeline has exactly the 3 expected columns; a dual-pipeline variant is still found via tool_name resolution despite its surviving maf row being tagged the other pipeline; mutation_id never leaks into caller-facing output; backward-compat fallback (DB built before variant_pipeline existed) still behaves correctly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Records the query-layer indexing work (Pipeline case mismatch, missing variant_pipeline/bedpe indexes, region-query consolidation) and the GAMBLR.open-side calc_mutation_frequency_bin_regions fix, with before/after timing from tools/fast_example_timer.R (~242s -> ~130s total, 142s -> 30.2s for the single largest item). Also updates the earlier variant_pipeline "still to do" note now that it's confirmed fixed on the GSC. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The laptop-vs-GSC comparison (same DB, same code, local disk instead of GSC storage) showed assign_cn_to_ssm dropping from 23.7s to under ~1.1s -- too large a swing for a purely CPU-bound cool_overlaps() join, meaning the GSC-side cost was actually dominated by disk/network I/O for its two underlying queries, not the join itself as originally assessed. calc_mutation_frequency_bin_regions shows the same pattern to a lesser degree. Updates the practical conclusion accordingly: further code-level query optimization has diminishing returns here, the larger remaining lever is GSC storage infrastructure, not this refactor. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Removes the last remaining dependency on GAMBLR.data::sample_data (the installed, bundled object) from assemble_bundled_data.R. By Phase 3 (where this block runs), the local sample_data\$meta already has everything needed -- every cohort, and the study column (renamed from cohort) -- so there's no longer a reason to reach into whatever happens to be installed, which could easily be out of sync with the current run (the exact class of staleness bug that caused the cohort/study mismatch fixed earlier this session). This also removes the last blocker to dropping sample_data.rda from the package entirely: nothing in this script depends on GAMBLR.data::sample_data existing anymore. Since this is a data-raw/ script (sourced directly, not loaded via library()), the fix takes effect on the very next run with no reinstall required, and is independent of whether/when sample_data.rda itself actually gets removed from the package. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
gambl_mutations_db() previously just errored if no local copy was found in
any of its resolved paths, requiring the file to be manually placed there
(scp/rsync -- how it's been distributed all session). Now falls back to
downloading it automatically via the new download_gambl_mutations_db(),
using the ropensci piggyback package to fetch it from this repo's GitHub
Release assets.
Downloads are tagged per package version by default ("data-v<version>"),
so a given install always fetches data built for its own schema rather
than whatever the newest release happens to be -- directly addresses the
"is my DB in sync with my installed code" confusion that caused real
friction earlier this session (e.g. the package-reinstall-vs-DB-rebuild
timing mismatches while chasing the variant_pipeline fix). Auto-download
can be disabled per-call (auto_download = FALSE) or globally
(options(GAMBLR.data.auto_download = FALSE)) for offline/CI environments,
falling back to the original explicit error.
piggyback is a Suggests, not a hard dependency -- most calls to
gambl_mutations_db() never need to download anything, and the function
gives a clear, actionable error if it's missing when actually needed.
Added data-raw/release_mutations_db.R, the corresponding admin-side upload
script (not run as part of this change -- requires a GitHub PAT with write
access to this repo; creates the release/tag if needed and uploads via
piggyback::pb_upload()). Bump DESCRIPTION's Version and run this manually
after each fresh gambl_mutations.db build to keep code and data in
lockstep.
Verified with synthetic tests: skip-if-exists short-circuit, clear error
when piggyback is absent, auto_download=FALSE preserves the exact original
error with no piggyback attempt, the global option override works, and an
explicit db_path still works unaffected.
Also regenerates man/NAMESPACE for several roxygen edits made earlier this
session (Pipeline lowercase normalization, region-query consolidation,
variant_pipeline/mutation_id schema docs) that were pending a
devtools::document() run.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The release tag is derived from DESCRIPTION's Version on disk, but every
OTHER install resolves the same tag from
utils::packageVersion("GAMBLR.data") -- i.e. whatever's actually committed
and installed elsewhere. Running this with a locally-bumped but
uncommitted/unpushed DESCRIPTION would create a release that no install
(including a future reinstall of this exact checkout) could ever find,
since they'd all keep resolving the old version. Now refuses to proceed
if DESCRIPTION has any uncommitted (staged or unstaged) changes, and prints
the commit it's releasing from as a reminder to make sure it's been pushed.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
… steps Documents the version-bump -> commit/push -> release_mutations_db.R sequence discussed in chat, plus prerequisites (piggyback, GitHub PAT) and a troubleshooting section. Deliberately not in README.md -- this is a maintainer-only workflow, not something package users need, but it's co-located with the build/release scripts in data-raw/ so it's easy to find for anyone actually working on the build process. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…download Gates the data-v1.4 GitHub Release tag used by gambl_mutations_db()'s auto-download path.
pb_upload() re-checks the release list to confirm the tag exists, which can momentarily miss a just-created release. piggyback's own retry fallback for this only kicks in when interactive(), which Rscript never is, so add a backoff retry here instead. Confirmed against a real "Release not found" failure while cutting data-v1.4.
devtools::install() was compiling RSQLite from CRAN source using this env's toolchain (gcc/g++ 11.2.0 vs libstdc++-ng runtime 13.2.0), which fails with "'timespec_get' has not been declared in '::'" in RSQLite's C++ SQLite wrapper. Not caused by this branch -- the same workflow passed on an unrelated branch two weeks ago with this same env file. Pinning a conda-forge prebuilt binary (r43-compatible, matching r-base=4.3.3) sidesteps the broken compile path entirely.
Conda-forge's linux-64 builds of r-rsqlite 2.4.2/2.4.3 all require libgcc/libstdcxx >=14, which conflicts with this env's pinned libgcc-ng=13.2.0 (confirmed via LibMambaUnsatisfiableError on the previous push). 2.4.1 only requires >=13, satisfied by the existing pin.
2.4.x builds all depend on the newer libgcc/libstdcxx metapackages (>=13 or >=14), which this env's libgcc-ng=13.2.0 pin doesn't satisfy (confirmed via two failed solves). 2.3.7's r43h0d4f4ea_0 build predates that migration and depends on libgcc-ng/libstdcxx-ng >=12 directly -- the same package names already pinned here, satisfied by 13.2.0. Pinning the full build string since 2.3.7 also has a newer-style build the solver could otherwise pick instead.
…etadata/gambl_metadata assemble_bundled_data.R's Phase 2 rejoin against live GAMBLR.results::get_gambl_metadata() only pulled sample_id/seq_type/cohort, so sample_meta (and, downstream, the now-retired sample_metadata bundled object) was missing ~40 columns the live function actually provides (COO_consensus, DHITsig_consensus, EBV_status_inf, ffpe_or_frozen, time_point, lymphgen_*, normal_sample_id, pairing_status, Tumor_Sample_Barcode, sex, bam_available, genome_build, etc.). Now keeps only genuinely paper-supplement-specific columns from each cohort block (sample_id, study, genetic_subgroup, reference_PMID) and left_joins the full, unrestricted live metadata by (sample_id, seq_type) instead of a narrow select(). Retires two bundled objects entirely, now that GAMBLR.open::get_gambl_metadata() (a separate commit) queries sample_meta from gambl_mutations.db directly: - sample_metadata: had exactly one real consumer anywhere in the GAMBLR.data/.open/.results/.utils/.viz ecosystem (verified precisely, not by substring grep). Drops its build step, data/sample_metadata.rda, R/sample_metadata.R, data-raw/make_sample_metadata.R (now fully dead), and its regression check in data-raw/test_gambl_db.R. - gambl_metadata: a ~3-year-stale (2023-08-28) object with zero real consumers and zero build step already -- fully orphaned dead weight. Drops data/gambl_metadata.rda and its doc entries in R/data.R / R/data-vars.R. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The previous commit widened sample_meta's rejoin against live GAMBLR.results::get_gambl_metadata() to fix missing columns, but captured everything unfiltered -- including internal-only columns like biopsy_id that must never end up in a publicly-distributed database. This matters beyond GAMBLR.open's R-level filtering: anyone can download gambl_mutations.db and query sample_meta directly with SQL, so the restriction has to happen at the point this table is built, not just in one consumer function. Adds data-raw/public_sample_meta_columns.txt: a plain-text whitelist (one column per line, # comments) that a developer can extend for future builds without touching code. assemble_bundled_data.R applies it right after the live metadata join, via select(..., all_of(cols)) rather than any_of() -- a typo'd/renamed column name fails the build loudly instead of silently vanishing or silently including something unintended. sample_id/seq_type/study/reference_PMID stay hardcoded (required for the pipeline's own internal joins/filters), plus genetic_subgroup (paper- supplement-specific, not something the live-metadata join would restore if dropped before it). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Debugging: a full rebuild is producing an unfiltered sample_meta (90 columns) despite the whitelist restriction in the previous commit -- verified correct in isolation (both the file read and the select() logic, tested live against real data on the GSC), and the source on disk matches what's committed byte-for-byte. Something between the Phase 2 fix and the actual dbWriteTable() call is not behaving as the code reads. Adds three [DIAG-COLS] checkpoints printing ncol()/colnames(): right after sample_data$meta is assigned in Phase 2, right before write_mutations_db() is called, and right before the sample_meta dbWriteTable() call inside it. Whichever pair first disagrees pinpoints exactly where the columns come back. Temporary -- remove once the root cause is found. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
They did their job: confirmed sample_meta genuinely has the correct, whitelist-restricted 34-column set at every stage through to the actual dbWriteTable() call. The "90 columns" symptom that prompted them was a stale/wrong file at the destination path during copy-to-laptop, not a bug in this script -- verified end to end with a fresh filename and a raw SQLite read bypassing GAMBLR.open entirely. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Contributor
Author
|
@Kdreval I'm still working on this so no need for review yet but it's looking promising for a lightweight and scalable replacement for our data bundling strategy |
Contributor
Ok, sounds good! Let me know when it’s ready for review! |
get_ssm_by_regions() now stamps maf_seq_type (GAMBLR.results), but the build script was dropping it via select(all_of(all_cols)), and write_mutations_db()'s dedup key and get_ssm_from_db() had no way to use it -- so a sample_id with both genome and capture data was still indistinguishable once it reached the shipped DB, even though the underlying data now carries the distinction. - assemble_bundled_data.R: keep maf_seq_type through all_cols; stamp it explicitly for the cell-line pulls (get_ssm_by_samples() doesn't set it, but cell_lines_data$meta is confirmed genome-only). Publication-pipeline rows are left NA -- no per-row seq_type provenance exists for them, and guessing would be worse than NA. - write_mutations_db.R: add maf_seq_type to the maf/ashm dedup key (so a genome/capture collision at the same position+allele survives as two rows instead of one) and index it. - get_ssm_from_db(): new seq_type filter param, no-ops safely against a gambl_mutations.db built before this column existed. Verified against synthetic sample_data (dual-seq_type collision, NA publication row, pre-existing-schema DB) -- see conversation for the check script. Not yet verified against a real GAMBLR.results rebuild.
…_to_keep selected_columns (used only for the proteinpainter-compatibility enrichment join sources around line 971, not the final maf table) derived its extra columns from maf_columns_to_keep, so it silently picked up maf_seq_type too once that was added there. Two of its three consumers never carry that column (a raw external flat-file read, and a get_ssm_by_samples() pull -- that function doesn't stamp maf_seq_type), so select(all_of(selected_columns)) hard-errored on the very first one. Give selected_columns its own explicit column list instead of inheriting from a list that now serves a different purpose.
maf and ashm dedup independently at write time, so nothing prevented the same real mutation (same sample/position/allele) landing in both tables when a gene is both aSHM-target and coding-panel (BCL2, BCL6, MYC, PIM1, etc.). Three GAMBLR.open accessors (get_ssm_by_samples(), get_ssm_by_regions(), get_ssm_by_region()) hardcode include_ashm=TRUE with no way for a caller to disable it, so this was reaching users on every call, not just as an opt-in edge case. Verified against the real (non-synthetic) rebuilt gambl_mutations.db: 59,457 of 85,884 ashm rows (69%) were exact duplicates of a maf row before this fix; after it, the combined maf+ashm result adds exactly the 26,427 genuinely aSHM-only mutations with zero remaining duplicate natural keys.
maf and ashm had identical schemas, independent mutation_id sequences,
and independent write-time dedup passes, despite nothing anywhere
querying ashm in isolation -- the only reference was
get_ssm_from_db()'s include_ashm toggle, which always bind_rows()'d it
onto maf. Confirmed against the real rebuilt DB: 59,457 of 85,884 ashm
rows (69%) exactly duplicated a maf row (genes that are both aSHM
targets and coding-panel genes: BCL2, BCL6, MYC, PIM1...). A same-
session query-time dplyr::distinct() patch fixed the symptom; this
removes the cause.
Considered and rejected keeping a filterable ashm-provenance flag: the
three GAMBLR.open accessors that ever set include_ashm=TRUE do so
unconditionally (their own regions/sample filters do the real
restricting), and the only place the FALSE default was ever exercised
is get_coding_ssm()/get_all_coding_ssm(), which hardcode
coding_only=TRUE -- checked the Variant_Classification breakdown of
the 26,427 aSHM-only rows: only 991 are coding-classified, the rest
already excluded by coding_only regardless of any aSHM flag. So plain
unconditional pooling isn't just simpler, it's more correct --
get_coding_ssm() gains those 991 previously-invisible real coding
mutations for free.
- assemble_bundled_data.R: aSHM pull rows fold into the same
bind_rows() that builds sample_data$grch37/hg38$maf; no more
separate $ashm slot.
- write_mutations_db.R: one dedup pass, one elem ("maf") instead of
two. variant_pipeline drops its now-always-constant elem column.
- get_ssm_from_db(): include_ashm parameter removed entirely, along
with the gather("ashm")/bind_rows()/cross-table distinct() added
earlier this session (structurally unnecessary now). Backward-compat
guard added for variant_pipeline's elem column in case of an
old-schema DB.
- gambl_mutations_db(): schema docs updated to match.
Bumped Version to 1.5 (schema-breaking, feeds the data-v<Version>
release-tag convention). Verified with a synthetic write_mutations_db()
run (overlapping coding-panel + aSHM row collapses to one; genuinely
aSHM-only row survives; no ashm table, no elem column, no ashm indexes,
get_ssm_from_db() has no include_ashm formal). Not yet verified against
a real GAMBLR.results rebuild.
Both scripts hardcoded assumptions from when maf and ashm were separate tables. Left unfixed, test_gambl_db.R would report false failures (ashm table/index checks) and compare_bundle_changes.R would hard-error outright (SELECT * FROM ashm, now nonexistent). - test_gambl_db.R: drop ashm from table/row-count/index/build_info checks; assert it does NOT exist instead. variant_pipeline's column check drops elem. The legacy sample_data.rda regression check can no longer expect exact equality for maf (the DB's maf now also absorbs what used to be ashm, deduplicated against it) -- replaced with a range check: old maf <= new maf <= old maf + old ashm. - compare_bundle_changes.R: SNV comparison reads maf alone on the new side instead of maf UNION ashm; .source is uniformly "maf" there since the new schema has no retained per-row provenance to trace gained/lost examples back to (the old side, from the .rda, is unaffected and still distinguishes them). Ran both against the real rebuilt gambl_mutations.db: test_gambl_db.R passes every check except the optional legacy-.rda comparisons, whose baseline (data/sample_data.rda, last updated 2025-06-30) is confirmed stale independent of this change -- seg and bedpe, untouched by any of today's work, mismatch it by similarly large margins. compare_bundle_changes.R runs clean; per-event diff against that same stale baseline shows 299 SNV / 0 CNV / 14 SV losses against thousands of gains out of ~889K total maf rows, consistent with ~13 months of accumulated cohort/sample growth rather than a regression.
data-v1.4 was only published a few days ago and this PR hasn't landed on master yet, so there's no meaningful population of installs depending on that exact release asset staying as-is. Reusing 1.4 avoids the awkwardness of a 1.5 bump for what's still all one unmerged, in-flight PR.
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
2 participants
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
.rdasample data togambl_mutations.db(SQLite), with tablesmaf/seg/bedpe/sample_meta/sample_study/variant_pipelinesample_id/Tumor_Sample_Barcodeoverwrite, duplicate SLMS-3 pulls for samples in overlapping cohorts, single-valuedStudycolumn unable to represent multi-cohort samples (replaced by asample_studyjoin table)tool_name="slms-3"queries. Fixed via avariant_pipelinejoin table keyed on a surrogatemutation_idPipelinefilters, bedpe position indexes) -- cuts the example-test suite from ~242s to ~130s (GSC)gambl_mutations.dbis now downloaded on demand via GitHub Releases (piggyback) instead of bundled in the package, with a documented maintainer release process (data-raw/RELEASING.md)Update: seq_type provenance + maf/ashm merge
Investigating whether a sample_id with both genome and capture data could
have its mutations queried separately surfaced two more real bugs, both
fixed here:
maf_seq_typecolumn added tomaf.GAMBLR.results::get_ssm_by_regions()silently merged a sample's genome- and capture-derived mutations with no
way to tell them apart afterward (fixed upstream in GAMBLR.results too).
get_ssm_from_db()gained aseq_typefilter param to query eitherindependently; the write-time dedup key now includes it so a genome call
and a capture call at the same position no longer collapse into one row.
mafandashmmerged into one table. They had identical schemas,independent
mutation_idsequences, and independent write-time dedup,despite nothing anywhere querying
ashmin isolation. Confirmed againsta real rebuild: 69% of
ashmrows exactly duplicated amafrow (genesthat are both aSHM targets and coding-panel genes -- BCL2, BCL6, MYC,
PIM1...). Three GAMBLR.open accessors hardcoded
include_ashm = TRUEwith no way to disable it, so every call to any of them was silently
returning duplicate mutations. Rather than add a filterable
aSHM-provenance flag, checked where the distinction would ever matter and
found nowhere: pooling unconditionally and dropping
get_ssm_from_db()'sinclude_ashmparameter entirely is both simpler and more correct --get_coding_ssm()now also returns ~991 real coding mutations thatpreviously only existed in
ashm, which it never queried.variant_pipelinedrops its now-constantelemcolumn (only ever neededto disambiguate
maf's andashm's independent ID spaces).data-v1.4was only published a few days ago andthis PR hasn't landed on master yet, so there's no install depending on
that exact release asset staying as-is. The
data-v1.4release will beupdated in place with the new merged-schema DB rather than cutting a
1.5.
Test plan
data-raw/test_gambl_db.Rpasses against the rebuiltgambl_mutations.db(updated for the maf/ashm merge; all checks pass except the optional legacysample_data.rdacomparison, whose baseline is confirmed stale -- last updated 2025-06-30, predates this whole migration)data-raw/compare_bundle_changes.Rruns clean (also updated for the merge); per-event diff against that same stale baseline shows far more gains than losses, consistent with ~13 months of accumulated cohort growth rather than a regressionFL_Dreval/DLBCL_cell_linesSLMS-3-tagged coding variants are no longer near-zerodata-v1.4GitHub Release originally cut and verified: fresh install + cleared cache +GAMBLR.open::get_coding_ssm()triggers auto-download and returns correct data with properPipelinetaggingmaf_seq_type: verified against real dual-seq_type samples (e.g. a sample with both genome and capture data shows correctly-separated, non-identical mutation counts per seq_type viaget_ssm_from_db(seq_type = ...))ashmtable,get_ssm_from_db()row count matches the pre-merge deduplicated total exactly, zero duplicate natural keys, spot-checkedBCL6/MYC(the classic overlap genes) viaget_ssm_by_regions()/get_ssm_by_region()with zero duplicatesdata-v1.4release asset not yet updated with the merged-schema DB🤖 Generated with Claude Code