diff --git a/.Rbuildignore b/.Rbuildignore index b7dc945..1203c55 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -2,3 +2,4 @@ ^\.Rproj\.user$ ^LICENSE\.md$ ^data-raw$ +^gambl_mutations\.db$ diff --git a/.gitignore b/.gitignore index cd67eac..cad1feb 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ .Rproj.user +gambl_mutations.db diff --git a/DESCRIPTION b/DESCRIPTION index 286baf9..0f95aa8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GAMBLR.data Title: Collection of Curated Data for Genomic Analysis of Mature B-cell Lymphomas in R -Version: 1.3 +Version: 1.4 Authors@R: c( person(given = "Ryan", @@ -25,15 +25,20 @@ Description: The package contains manually curated data for the genomic Analysis License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 Depends: R (>= 3.5) LazyData: true Imports: + DBI, + dbplyr, dplyr, ggplot2, + rlang, + RSQLite, tidyr LazyDataCompression: xz Suggests: + piggyback, testthat (>= 3.0.0) Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index 6381e6b..03d6819 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,12 @@ # Generated by roxygen2: do not edit by hand +export(download_gambl_mutations_db) +export(gambl_mutations_db) +export(gambl_reference_db) +export(get_cn_segments_from_db) +export(get_genes) +export(get_ssm_from_db) +export(get_sv_from_db) import(dplyr) import(ggplot2) import(tidyr) diff --git a/R/data-vars.R b/R/data-vars.R index 11cf5ab..f2c79e0 100644 --- a/R/data-vars.R +++ b/R/data-vars.R @@ -19,12 +19,13 @@ if (getRversion() >= "2.15.1") { "SCORE", "START_A", "START_B", "Start_Position", "Start_Position.x", "Tumor_Sample_Barcode", "VAF_tumour", "Variant_Classification", "bin", "category", "chrom", "cohort", "colour", "curated", "end", - "ensembl_gene_id", "gambl_metadata", "gene", "genome_build", + "ensembl_gene_id", "gene", "genome_build", "grch37_ashm_regions", "group", "head", "hg38_ashm_regions", - "hot_spot", "hotspot_regions_grch37", "hotspot_regions_hg38", - "is_alias", "log.ratio", "mutated", "mutation_count", "n_mut", - "name", "pair_status", "pathology", "patient_id", "region", - "row_id", "sample_id", "seq_type", "start", "t_alt_count", - "tumour_sample_id", "window_end", "window_start" + "hot_spot", "hotspot_regions_grch37", "hotspot_regions_hg38", + "is_alias", "log.ratio", "maf_seq_type", "mutated", "mutation_count", "n_mut", + "name", "pair_status", "pathology", "patient_id", "region", + "row_id", "sample_id", "seq_type", "start", "t_alt_count", + "tumour_sample_id", "window_end", "window_start", + "study", "Pipeline", "elem" )) } diff --git a/R/data.R b/R/data.R index 02179a8..dad8c53 100644 --- a/R/data.R +++ b/R/data.R @@ -933,45 +933,6 @@ "somatic_hypermutation_locations_GRCh38_v_latest" -#' GAMBL Metadata -#' -#' A data frame with metadata for a collection of GAMBL samples. -#' -#' @format ## `gambl_metadata` -#' A data frame with 4785 rows and 27 columns -#' \describe{ -#' \item{sample_id}{Sample identifier.} -#' \item{patient_id}{Patient identifier.} -#' \item{pathology}{Pathology.} -#' \item{seq_type}{Sample sequencing type.} -#' \item{genome_build}{Genome build the sample coordinates are in reference to.} -#' \item{pairing_status}{Matched or unmatched.} -#' \item{Tumor_Sample_Barcode}{Sample ID in another column, needed for certain functions.} -#' \item{age_group}{Sample age group.} -#' \item{compression}{The compression available for a particular sample.} -#' \item{bam_available}{Boolean.} -#' \item{pathology_rank}{Pathology rank.} -#' \item{cohort}{Sample cohort} -#' \item{COO_consensus}{COO consensus.} -#' \item{DHITsig_consensus}{DHIT signature consensus.} -#' \item{EBV_status_inf}{EBV status.} -#' \item{ffpe_or_frozen}{FFPE or frozen.} -#' \item{fl_grade}{FL grade.} -#' \item{hiv_status}{Sample HIV status.} -#' \item{lymphgen}{Lymphgen.} -#' \item{lymphgen_cnv_noA53}{Lymphgen with CNV no A53.} -#' \item{lymphgen_no_cnv}{Lymphgen no CNV.} -#' \item{lymphgen_with_cnv}{Lymphgen with CNV} -#' \item{lymphgen_wright}{Lymphgen Wright.} -#' \item{molecular_BL}{Molecualr BL.} -#' \item{normal_sample_id}{Normal sample ID} -#' \item{sex}{Female or Male} -#' \item{time_point}{Smaple timepoint.} -#' } -#' @keywords internal -"gambl_metadata" - - #' Mirage Metrics #' #' A data frame with Quality Control metrics for a collection of GAMBL samples. diff --git a/R/data_comp.R b/R/data_comp.R index a9d3d58..ae31fdf 100644 --- a/R/data_comp.R +++ b/R/data_comp.R @@ -3,197 +3,124 @@ #' #' @description Retrieve gene names from bundled lymphoma gene lists. #' -#' @details Complete lists of genes described as significantly mutated in large -#' lymphoma studies have been curated and provided with this package. +#' @details Gene lists are read from the bundled normalized reference database +#' (`inst/extdata/gambl_reference.db`), which is built from the LLMPP curated +#' lists. Confidence tier and data version are columns in that database, so the +#' previous per-entity, per-version `.rda` objects are no longer required. #' #' @param entities Optional vector specifying one or more lymphoma entities #' e.g. MCL, DLBCL, BL. -#' @param curated_only Specify FALSE to retrieve all genes or leave default -#' for the curated subset. +#' @param curated_only When TRUE (default) return the curated set (Tier 1 and +#' Tier 2) and drop retired Tier 3 genes. Set FALSE to include all tiers. +#' Ignored when `tier` is supplied. +#' @param tier Optional integer vector of confidence tiers to keep, e.g. +#' `tier = 1` for high-confidence genes only, or `tier = c(1, 2)`. When +#' supplied it takes precedence over `curated_only`. #' @param gene_format Specify what to return as output. Can be one of: #' * "symbol" (the default): list of gene symbols -#' * "ensembl": list of ENSEMBLE IDs -#' * "data.frame": data frame with column Gene and per-entity gene status -#' @param version Specify which version to return. Currently supported versions -#' are 0.0 (legacy version from original GAMBLR), 0.1, and _latest. The -#' latter will always point to the highest numeric version of the genes. +#' * "ensembl": list of Ensembl IDs +#' * "data.frame": one row per gene with a `_Tier` column per +#' requested entity (the gene's tier there, or NA if absent). Filter +#' these columns to focus on tiers, e.g. +#' `dplyr::filter(df, DLBCL_Tier == 1, FL_Tier == 1) |> dplyr::pull(Gene)`. +#' @param version Data version to return. Defaults to "_latest", the current +#' LLMPP-sourced version in the reference database. Legacy numeric versions +#' (0.0, 0.1) are no longer bundled. #' -#' @return A character vector of gene symbol or Ensembl IDs or a data frame. +#' @return A character vector of gene symbols or Ensembl IDs, or a data frame. #' #' @import dplyr tidyr #' +#' @export +#' #' @examples -#' all_dlbcl_genes <- get_genes(entities = "DLBCL", curated_only = FALSE) +#' # high-confidence (Tier 1) DLBCL genes +#' dlbcl_tier1 <- get_genes(entities = "DLBCL", tier = 1) +#' +#' # genes that are Tier 1 in DLBCL, FL and BL simultaneously +#' gene_df <- get_genes(entities = c("DLBCL", "FL", "BL"), gene_format = "data.frame") +#' shared_tier1 <- dplyr::filter(gene_df, DLBCL_Tier == 1, FL_Tier == 1, BL_Tier == 1) +#' #' all_curated_genes <- get_genes() get_genes <- function( entities = c("DLBCL", "MCL", "BL"), curated_only = TRUE, + tier = NULL, gene_format = "symbol", version = "_latest" ) { - # We need to handle the legacy version (0.0) separately - # because it was not pathology-specific - if (version == "0.0") { + # Gene lists are sourced from the normalized reference database (built from + # the LLMPP curated lists) rather than per-entity, per-version .rda objects. + # "version" and "tier" are columns, so lymphoma_genes_*_v0.1/v0.2/v_latest + # are no longer needed. + if (!version %in% c("_latest", "latest")) { message( paste( - "You have requested the legacy version of the lymphoma genes.", - "It is only provided here for backwards compatibility.", - "We recommend you use the most recent version by keeping the", - "default value of the argument version. Use at your own risk." + "Legacy per-version gene lists (e.g. 0.0, 0.1) are no longer", + "bundled. Data is now sourced from the LLMPP curated lists;", + "returning the current version." ) ) - - legacy_lymphoma_genes <- eval( - parse( - text = paste0( - "GAMBLR.data::", - "lymphoma_genes_lymphoma_genes_v0.0") - ) - ) - - legacy_lymphoma_genes <- legacy_lymphoma_genes %>% - dplyr::select( - Gene, ensembl_gene_id, - intersect( - colnames(.), - entities - ) - ) - if (curated_only) { - # drop any row where all pathologies have FALSE - legacy_lymphoma_genes <- legacy_lymphoma_genes %>% - dplyr::filter( - ! if_all(3:ncol(.), ~ . == "FALSE") - ) - } - - if (gene_format == "symbol") { - return(legacy_lymphoma_genes$Gene) - } else if (gene_format == "ensembl") { - return(legacy_lymphoma_genes$ensembl_gene_id) - } else if (gene_format == "data.frame") { - return(legacy_lymphoma_genes) - } else { - stop( - "You requested output format that is not supported." - ) - } } - #construct file name using entity - entities <- tolower(entities) - - r_objects <- entities %>% - paste0( - "lymphoma_genes_", - ., - "_v", - version - ) - - # check for unsupported gene sets - all_files <- system.file( - "extdata", - package = "GAMBLR.data" - ) %>% - list.files( - recursive = TRUE, - full.names = TRUE - ) - - all_files <- gsub( - ".*extdata/", - "", - all_files - ) - - all_files <- all_files[grepl("lymphoma_genes", all_files)] - - all_files <- all_files[grepl(version, all_files)] - - available_entities <- gsub("(.*/\\s*(.*$))", "\\2", all_files) + # `tier`, when supplied, selects exactly those confidence tiers (e.g. + # tier = 1 for high-confidence genes only) and takes precedence over + # `curated_only`. Otherwise curated_only == TRUE keeps the curated set + # (Tier 1 + Tier 2) and drops retired Tier 3 genes; FALSE returns all tiers. + keep_tiers <- if (!is.null(tier)) { + as.integer(tier) + } else if (curated_only) { + c(1L, 2L) + } else { + c(1L, 2L, 3L) + } - available_entities <- gsub(".tsv", "", available_entities) + con <- gambl_reference_db() + entities_u <- toupper(entities) - missing_sets <- setdiff( - entities, - available_entities - ) + available <- dplyr::tbl(con, "gene_entity") %>% + dplyr::distinct(entity) %>% + dplyr::pull(entity) + missing_sets <- setdiff(entities_u, available) if (length(missing_sets) > 0) { warning( paste( "The gene set for the entity", - missing_sets, + paste(missing_sets, collapse = ", "), "is not available and will not be returned." ) ) - r_objects <- r_objects[ - grepl( - paste( - available_entities, - collapse="|" - ), - r_objects - ) - ] - entities <- entities[entities %in% available_entities] + entities_u <- intersect(entities_u, available) } - # Combine all lists into one df - # Do it in a way that when GAMBLR.data is ot imported, the data objects - # are still available - all_entities_data <- list() - - for (i in seq_along(r_objects)) { - all_entities_data[[i]] <- eval( - parse( - text = paste0( - "GAMBLR.data::", - r_objects[i] - ) - ) - ) - } - - all_entities_data <- all_entities_data %>% - # only select necessary columns - lapply( - ., - `[`, - , - c("ensembl_gene_id", "Gene", "curated") - ) - - names(all_entities_data) <- toupper(entities) - - all_entities_data <- all_entities_data %>% - bind_rows(.id = "entity") - - if (curated_only) { - # drop any row where curated is FALSE - all_entities_data <- all_entities_data %>% - dplyr::filter( - curated == "TRUE" - ) - } + dat <- dplyr::tbl(con, "gene_entity") %>% + dplyr::filter(entity %in% entities_u, tier %in% keep_tiers) %>% + dplyr::select(entity, ensembl_gene_id, Gene = gene, tier) %>% + dplyr::collect() if (gene_format == "symbol") { - return(all_entities_data$Gene %>% unique %>% sort) + return(dat$Gene %>% unique %>% sort) } else if (gene_format == "ensembl") { - return(all_entities_data$ensembl_gene_id %>% unique %>% sort) + return(dat$ensembl_gene_id %>% unique %>% sort) } else if (gene_format == "data.frame") { + # one row per gene; one _Tier column per requested entity + # holding the gene's tier there (NA if absent), so callers can do + # filter(df, DLBCL_Tier == 1, FL_Tier == 1, ...) %>% pull(Gene) return( - all_entities_data %>% - select(-curated) %>% - mutate(is_gene = "TRUE") %>% - tidyr::pivot_wider( - names_from = "entity", - values_from = "is_gene" - ) %>% - replace(is.na(.), "FALSE") + dat %>% + dplyr::distinct(Gene, ensembl_gene_id, entity, tier) %>% + tidyr::pivot_wider( + id_cols = c(Gene, ensembl_gene_id), + names_from = "entity", + values_from = "tier", + names_glue = "{entity}_Tier" + ) %>% + dplyr::select(Gene, ensembl_gene_id, + dplyr::any_of(paste0(entities_u, "_Tier"))) %>% + dplyr::arrange(Gene) ) } else { stop( diff --git a/R/download_gambl_mutations_db.R b/R/download_gambl_mutations_db.R new file mode 100644 index 0000000..c202237 --- /dev/null +++ b/R/download_gambl_mutations_db.R @@ -0,0 +1,69 @@ +#' @title Download the GAMBL mutations database. +#' +#' @description Downloads `gambl_mutations.db` from this package's GitHub +#' Release assets, if not already present. Normally called automatically by +#' [gambl_mutations_db()] the first time it's needed; call it directly to +#' pre-fetch the file (e.g. before a batch job on a machine with no other +#' network access) or to force a re-download. +#' +#' @details Requires the `piggyback` package (Suggests, not a hard +#' dependency of GAMBLR.data -- most calls to [gambl_mutations_db()] never +#' need to download anything). Downloads are tagged per package version by +#' default (`"data-v"`), so a given install always fetches data +#' built for its own schema, not whatever happens to be the newest release +#' -- see `data-raw/release_mutations_db.R` for the corresponding upload +#' step (admin-only, run manually after a fresh build; bump +#' `DESCRIPTION`'s `Version` and cut a matching release whenever the schema +#' changes, so the two stay in lockstep). +#' +#' @param dest_path Where to save the downloaded file. Defaults to the +#' standard cache location, `tools::R_user_dir("GAMBLR.data", "cache")`. +#' @param tag GitHub release tag to download from. Defaults to +#' `paste0("data-v", utils::packageVersion("GAMBLR.data"))`. +#' @param repo GitHub repo in `"owner/repo"` form. +#' @param overwrite Re-download even if `dest_path` already exists. Default +#' `FALSE`. +#' +#' @return `dest_path`, invisibly. +#' +#' @export +download_gambl_mutations_db <- function( + dest_path = file.path(tools::R_user_dir("GAMBLR.data", "cache"), "gambl_mutations.db"), + tag = paste0("data-v", utils::packageVersion("GAMBLR.data")), + repo = "morinlab/GAMBLR.data", + overwrite = FALSE) { + + if (file.exists(dest_path) && !overwrite) { + message("gambl_mutations.db already exists at ", dest_path, + "; set overwrite = TRUE to re-download.") + return(invisible(dest_path)) + } + + if (!requireNamespace("piggyback", quietly = TRUE)) { + stop("Downloading gambl_mutations.db requires the 'piggyback' package. ", + "Install it with install.packages(\"piggyback\"), or point ", + "gambl_mutations_db() at an existing copy via the db_path argument, ", + "options(GAMBLR.data.mutations_db = ...), or the GAMBLR_MUTATIONS_DB ", + "environment variable.", call. = FALSE) + } + + dir.create(dirname(dest_path), recursive = TRUE, showWarnings = FALSE) + message(sprintf( + "Downloading gambl_mutations.db from %s @ %s into %s -- this is a large file (several hundred MB) and may take a while...", + repo, tag, dirname(dest_path) + )) + piggyback::pb_download( + "gambl_mutations.db", + repo = repo, + tag = tag, + dest = dirname(dest_path), + overwrite = TRUE + ) + if (!file.exists(dest_path)) { + stop("Download appeared to succeed but ", dest_path, " still doesn't exist -- ", + "check that the release tag \"", tag, "\" actually has a gambl_mutations.db asset.", + call. = FALSE) + } + message("Done: ", dest_path) + invisible(dest_path) +} diff --git a/R/gambl_mutations_db.R b/R/gambl_mutations_db.R new file mode 100644 index 0000000..2722808 --- /dev/null +++ b/R/gambl_mutations_db.R @@ -0,0 +1,147 @@ +#' @title Connect to the GAMBL mutations database. +#' +#' @description Returns a (cached) read-only DBI connection to `gambl_mutations.db`, +#' the SQLite database holding the large sample-level frames (MAF, seg, +#' bedpe) built from the in-memory `sample_data` object assembled by +#' `data-raw/assemble_bundled_data.R`, via the shared `write_mutations_db()` +#' helper (`data-raw/write_mutations_db.R`). +#' +#' @details Unlike the small bundled [gambl_reference_db()], this file is large +#' (several hundred MB) and is NOT shipped inside the package. It is +#' distributed as a GitHub Release asset (see [download_gambl_mutations_db()]) +#' and cached locally. The path is resolved in this order: +#' \enumerate{ +#' \item the `db_path` argument; +#' \item `getOption("GAMBLR.data.mutations_db")`; +#' \item the `GAMBLR_MUTATIONS_DB` environment variable; +#' \item the user cache dir, `tools::R_user_dir("GAMBLR.data", "cache")`; +#' \item `gambl_mutations.db` in the working directory (development). +#' } +#' If none of these exist and `auto_download` is `TRUE` (the default), the +#' file is downloaded automatically into the user cache dir via +#' [download_gambl_mutations_db()]. Set `auto_download = FALSE`, or +#' `options(GAMBLR.data.auto_download = FALSE)` to disable this globally +#' (e.g. offline/CI environments) and get an explicit error instead. +#' +#' # Table reference +#' +#' Both genome builds are stacked into the same table (not split across +#' separate tables); filter on `genome_build` (`"grch37"` / `"hg38"`) to select +#' one. `maf` holds both coding-panel and aSHM-region mutations, pooled +#' unconditionally at write time (a mutation found by both pulls -- common +#' for genes that are both aSHM targets and coding-panel genes, e.g. BCL2, +#' BCL6, MYC, PIM1 -- collapses to one row, not two; see +#' [get_ssm_from_db()]). +#' +#' | Table | Rows (typical) | Grain | Key columns | +#' | --- | --- | --- | --- | +#' | `maf` | ~3.7M | one somatic mutation call | `mutation_id` (surrogate integer key, unique within `maf` across both genome builds -- assigned at write time, after deduplication; `variant_pipeline` references it as a foreign key), `Hugo_Symbol, Chromosome, Start_Position, End_Position, Tumor_Sample_Barcode, Variant_Classification, HGVSp_Short, t_alt_count, n_alt_count, Pipeline (lowercase, e.g. `"slms-3"`/`"publication"` -- normalized at write time so it can be matched with a plain, indexed equality), genome_build, maf_seq_type` (`"genome"`/`"capture"`, so a Tumor_Sample_Barcode with both isn't indistinguishable once merged; `NA` for publication-pipeline rows, which don't carry per-row seq_type provenance -- see [get_ssm_from_db()]'s `seq_type` parameter) (+ ~35 more MAF-standard columns) | +#' | `seg` | ~126k | one copy-number segment | `ID, chrom, start, end, LOH_flag, log.ratio, CN, genome_build` | +#' | `bedpe` | ~900 | one Manta structural-variant breakpoint pair | `CHROM_A, START_A, END_A, CHROM_B, START_B, END_B, manta_name, SCORE, STRAND_A, STRAND_B, tumour_sample_id, normal_sample_id, VAF_tumour, DP, pair_status, FILTER, genome_build` | +#' | `sample_meta` | ~3.3k | one sample | `patient_id, sample_id, Tumor_Sample_Barcode, seq_type, pathology, cohort, study, ...` (its own `genome_build` column records the sample's native alignment build, unrelated to the per-row build stamp used in the other tables) | +#' | `sample_study` | varies | one (sample, study) membership fact -- many-to-many, a sample belonging to N studies is N rows | `sample_id, study, study_id, reference_PMID`. `study_id` is that study's own identifier for the sample where it differs from GAMBL's `sample_id` (e.g. a paper's own case/patient ID); NA where no study-specific ID has been sourced. | +#' | `variant_pipeline` | varies | one (mutation, Pipeline) fact -- many-to-many, a variant independently called by N pipelines is N rows | `mutation_id, Pipeline`. `mutation_id` is a foreign key to `maf.mutation_id` -- not a natural key, so this table doesn't repeat the variant's own position/sample columns. | +#' | `build_info` | 9 | key/value | provenance (`source`, `built_at`, `builder`) and expected row counts (`n_maf`, `n_seg`, `n_bedpe`, `n_samples`, `n_sample_study`, `n_variant_pipeline`) used by the `data-raw/test_gambl_db.R` regression checks | +#' +#' `maf` does NOT carry a `Study` column. Cohort/study membership is +#' tracked once per sample in `sample_study`, not once per mutation row -- +#' a mutation row's `Study` couldn't represent a sample belonging to more +#' than one study anyway, and a per-row tag caused the same sample's +#' mutations to be pulled and duplicated once per cohort-specific code path +#' that claimed it (see `GAMBLR.data::get_ssm_from_db()`'s `this_study` +#' parameter for how to filter by study post-refactor). +#' +#' Similarly, `maf` only ever carries ONE `Pipeline` value per variant, +#' even though the same real mutation can be independently produced by more +#' than one pipeline (e.g. a cohort's own published/curated maf and our +#' separate SLMS-3 recall both calling the same position) -- when that +#' happens, `write_mutations_db()`'s write-time dedup keeps one row +#' arbitrarily (whichever pipeline's data was assembled first), which would +#' otherwise make the variant invisible to a query for the *other* pipeline +#' even though it's still fully present in the table. `variant_pipeline` +#' records every pipeline that actually produced each variant, independent +#' of which one "won" as `maf`'s single representative row; see +#' `GAMBLR.data::get_ssm_from_db()`'s `tool_name` parameter, which resolves +#' against this table rather than `maf`'s own `Pipeline` column. +#' +#' ## Join keys +#' `maf` and `bedpe` (via `tumour_sample_id`) key to +#' `sample_meta$Tumor_Sample_Barcode` / `sample_meta$sample_id`; `seg$ID` keys +#' to `sample_meta$sample_id`; `sample_study$sample_id` keys to +#' `sample_meta$sample_id`; `variant_pipeline$mutation_id` keys to +#' `maf.mutation_id`. +#' +#' ## Indexes +#' `maf`: `(genome_build, Chromosome, Start_Position)`, +#' `Tumor_Sample_Barcode`, `Pipeline`, `maf_seq_type`, `mutation_id`. No index on +#' `Hugo_Symbol`: gene-restricted queries resolve the gene to a region first +#' (the same logic used to populate this table) and filter on +#' `(genome_build, Chromosome, Start_Position)` instead -- see +#' `GAMBLR.data::get_ssm_from_db()`. `seg`: `(genome_build, ID)`, +#' `(genome_build, chrom, start)`. `bedpe`: `(tumour_sample_id, +#' genome_build)`, `(genome_build, CHROM_A, START_A)`, `(genome_build, +#' CHROM_B, START_B)` (one per breakpoint end, so a region search that OR's +#' both ends can use a different index per side). +#' `sample_meta`: `sample_id`, `Tumor_Sample_Barcode`. `sample_study`: +#' `sample_id`, `study`. `variant_pipeline`: `mutation_id`, `Pipeline`. +#' +#' @param db_path Optional explicit path to the .db file. +#' @param auto_download When no existing copy is found, download it +#' automatically via [download_gambl_mutations_db()]. Default `TRUE`; +#' can also be disabled globally with +#' `options(GAMBLR.data.auto_download = FALSE)`. +#' +#' @return A read-only DBIConnection. +#' +#' @examples +#' \dontrun{ +#' con <- gambl_mutations_db() +#' DBI::dbListTables(con) +#' +#' # coding mutations in TP53, grch37, one pipeline (Pipeline is stored +#' # lowercase -- see the Pipeline column note above) +#' dplyr::tbl(con, "maf") |> +#' dplyr::filter(genome_build == "grch37", Hugo_Symbol == "TP53", +#' Pipeline == "slms-3") |> +#' dplyr::collect() +#' } +#' @export +gambl_mutations_db <- function(db_path = NULL, auto_download = TRUE) { + if (is.null(db_path)) { + opt <- getOption("GAMBLR.data.mutations_db", default = NULL) + env <- Sys.getenv("GAMBLR_MUTATIONS_DB", unset = "") + cache_path <- file.path(tools::R_user_dir("GAMBLR.data", "cache"), "gambl_mutations.db") + candidates <- c( + opt, + if (nzchar(env)) env, + cache_path, + "gambl_mutations.db" + ) + candidates <- candidates[!is.null(candidates)] + hit <- candidates[file.exists(candidates)] + if (length(hit) == 0) { + can_auto_download <- isTRUE(auto_download) && isTRUE(getOption("GAMBLR.data.auto_download", TRUE)) + if (can_auto_download) { + db_path <- download_gambl_mutations_db(dest_path = cache_path) + } else { + stop("gambl_mutations.db not found. Set options(GAMBLR.data.mutations_db = \"/path\"), ", + "the GAMBLR_MUTATIONS_DB env var, place it in ", + tools::R_user_dir("GAMBLR.data", "cache"), + ", or call gambl_mutations_db() with auto_download = TRUE (the default) ", + "to fetch it automatically.", call. = FALSE) + } + } else { + db_path <- hit[1] + } + } + stopifnot(file.exists(db_path)) + + cached <- getOption("gamblr.mutations.con") + if (!is.null(cached) && DBI::dbIsValid(cached) && + identical(normalizePath(cached@dbname), normalizePath(db_path))) { + return(cached) + } + con <- DBI::dbConnect(RSQLite::SQLite(), db_path, flags = RSQLite::SQLITE_RO) + options(gamblr.mutations.con = con) + con +} diff --git a/R/get_cn_segments_from_db.R b/R/get_cn_segments_from_db.R new file mode 100644 index 0000000..3546fc5 --- /dev/null +++ b/R/get_cn_segments_from_db.R @@ -0,0 +1,25 @@ +#' @title Query copy-number segments from the mutations database. +#' +#' @description Data-access helper that retrieves copy-number segments from the +#' `seg` table of `gambl_mutations.db`. Backs the GAMBLR.open `get_cn_segments` +#' and `get_sample_cn_segments` accessors. +#' +#' @param projection Genome build / `genome_build` column value. Default "grch37". +#' @param sample_ids Optional character vector of segment `ID`s (sample ids) to keep. +#' @param con Optional DBI connection (defaults to [gambl_mutations_db()]). +#' +#' @return A data frame of seg rows (the `genome_build` helper column is dropped; +#' chromosome-prefix handling is left to the caller, matching the bundled seg). +#' +#' @export +get_cn_segments_from_db <- function(projection = "grch37", + sample_ids = NULL, + con = NULL) { + if (is.null(con)) con <- gambl_mutations_db() + q <- dplyr::tbl(con, "seg") %>% + dplyr::filter(genome_build == projection) + if (!is.null(sample_ids)) q <- dplyr::filter(q, ID %in% sample_ids) + res <- dplyr::collect(q) + res$genome_build <- NULL + res +} diff --git a/R/get_genes_db.R b/R/get_genes_db.R new file mode 100644 index 0000000..b291f8f --- /dev/null +++ b/R/get_genes_db.R @@ -0,0 +1,54 @@ +#' @title Connect to the bundled GAMBLR reference database. +#' +#' @description Returns a (cached) read-only DBI connection to the normalized +#' SQLite reference database built from the LLMPP curated gene lists. This is the +#' shared backing store used by [get_genes()] (and intended for the Lymphopedia +#' wiki / data API), replacing the versioned `lymphoma_genes_*` .rda objects. +#' +#' @details The database is rebuilt from source by `data-raw/build_reference_db.R` +#' and is small (well under 1 MB), so — unlike [gambl_mutations_db()] — it is +#' bundled inside the package (`inst/extdata/gambl_reference.db`). +#' +#' # Table reference +#' +#' | Table | Grain | Key columns | +#' | --- | --- | --- | +#' | `gene_entity` | one row per gene x lymphoma entity | `gene, ensembl_gene_id, entity, tier, ashm, qc, citekey, pmid, mutation_effect, extra` (JSON of entity-specific fields), `source_version` | +#' | `ashm_regions` | one row per gene x aSHM region x genome build | `gene, chrom, start, end, region, regulatory_comment, build, source_version` | +#' | `build_info` | key/value | provenance: `source` (LLMPP resources/curated), `source_version`, `built_at`, `builder` | +#' +#' `gene_entity$entity` takes values like `"DLBCL"`, `"BL"`, `"FL"`, `"PMBL"`, +#' `"MCL"`, `"MZL"`; `tier` is `1` (high-confidence), `2` (lower-confidence), or +#' `3` (retired). `ashm_regions$build` is `"grch37"` or `"hg38"`. +#' +#' @param db_path Optional explicit path to the .db file. Defaults to the copy +#' installed with the package; falls back to `inst/extdata/gambl_reference.db` +#' during development. +#' +#' @return A DBIConnection. +#' +#' @examples +#' \dontrun{ +#' con <- gambl_reference_db() +#' dplyr::tbl(con, "gene_entity") +#' +#' # high-confidence DLBCL genes +#' dplyr::tbl(con, "gene_entity") |> +#' dplyr::filter(entity == "DLBCL", tier == 1) |> +#' dplyr::collect() +#' } +#' @export +gambl_reference_db <- function(db_path = NULL) { + if (is.null(db_path)) { + db_path <- system.file("extdata", "gambl_reference.db", package = "GAMBLR.data") + if (!nzchar(db_path) || !file.exists(db_path)) { + db_path <- file.path("inst", "extdata", "gambl_reference.db") + } + } + stopifnot(file.exists(db_path)) + cache <- getOption("gamblr.reference.con") + if (!is.null(cache) && DBI::dbIsValid(cache)) return(cache) + con <- DBI::dbConnect(RSQLite::SQLite(), db_path, flags = RSQLite::SQLITE_RO) + options(gamblr.reference.con = con) + con +} diff --git a/R/get_ssm_from_db.R b/R/get_ssm_from_db.R new file mode 100644 index 0000000..0749282 --- /dev/null +++ b/R/get_ssm_from_db.R @@ -0,0 +1,155 @@ +#' @title Query SSMs from the mutations database. +#' +#' @description Shared data-access helper that retrieves simple somatic mutations +#' (SSM, MAF format) from `gambl_mutations.db`, pushing the common filters down to +#' indexed SQL. This is the single query path behind the GAMBLR.open SSM +#' accessors (get_coding_ssm, get_ssm_by_samples, get_ssm_by_region(s), ...). +#' +#' @param projection Genome build / `genome_build` column value. Default "grch37". +#' @param sample_ids Optional character vector of `Tumor_Sample_Barcode` to keep. +#' @param tool_name Pipeline to keep, matched case-insensitively (the value +#' you pass is lowercased before comparing; the stored Pipeline column is +#' itself normalized to lowercase at write time, so this is a plain, +#' indexed equality check rather than a function-wrapped one). Default +#' "slms-3"; set NULL to skip the Pipeline filter. Resolved via the +#' `variant_pipeline` join table (see [gambl_mutations_db()]'s schema +#' docs) rather than `maf`'s own Pipeline column, so a variant +#' independently called by more than one pipeline is still found by a +#' query for either one, even though its single maf row can only +#' carry one Pipeline value. Falls back to filtering the column directly +#' against DBs built before variant_pipeline existed. +#' @param coding_only When TRUE, keep only coding `Variant_Classification`s. +#' @param include_silent When FALSE (and `coding_only`), drop Silent mutations. +#' @param min_read_support Keep only variants with `t_alt_count` >= this value. +#' @param this_study Optional single study to restrict to, matched against +#' `sample_study$study` (not a `maf` column -- see +#' [gambl_mutations_db()]'s schema docs). Resolves to the set of sample_ids +#' belonging to that study, intersected with `sample_ids` if both are +#' supplied; composes with `tool_name`/Pipeline exactly as before (e.g. +#' `this_study="Reddy"` alone returns Reddy's SLMS-3 recall by default, +#' while `this_study="Reddy", tool_name="publication"` returns Reddy's +#' as-published rows only). +#' @param seq_type Optional character vector to restrict rows to one or more +#' `maf_seq_type` values (`"genome"`/`"capture"`), so a sample_id with both +#' can be queried separately instead of getting both back merged together. +#' Default `NULL` (no filter, matching prior behaviour). Rows predating +#' `maf_seq_type` (e.g. publication-pipeline rows, which don't carry +#' per-row seq_type provenance) have `NA` here and are excluded whenever +#' this filter is set, same as any other `NA` under a `%in%` filter. Silently +#' ignored against a `gambl_mutations.db` built before this column existed. +#' @param regions Optional data frame with columns `chrom`, `start`, `end`; +#' rows are OR-ed and applied as a single combined query (not one query per +#' region), so passing hundreds of regions (e.g. one per gene in a panel) +#' doesn't cost hundreds of separate round trips. +#' @param con Optional DBI connection (defaults to [gambl_mutations_db()]). +#' +#' @return A data frame of MAF rows (the `genome_build`/`mutation_id` helper +#' columns are dropped). +#' +#' @export +get_ssm_from_db <- function(projection = "grch37", + sample_ids = NULL, + tool_name = "slms-3", + coding_only = FALSE, + include_silent = TRUE, + min_read_support = 0, + this_study = NULL, + regions = NULL, + seq_type = NULL, + con = NULL) { + if (is.null(con)) con <- gambl_mutations_db() + cc <- if (include_silent) coding_class else coding_class[coding_class != "Silent"] + tn <- if (!is.null(tool_name)) tolower(tool_name) else NULL + + # this_study resolves to sample_ids via a join against sample_study, done + # once here rather than as a per-query filter inside base_query() -- reuses + # the existing sample_ids/Tumor_Sample_Barcode machinery below instead of + # requiring maf to carry its own Study column (it doesn't; see + # gambl_mutations_db()'s schema docs), and keeps this_study (which + # samples) and tool_name/Pipeline (which rows for those samples) composable + # exactly as they were. + if (!is.null(this_study)) { + study_sample_ids <- dplyr::tbl(con, "sample_study") %>% + dplyr::filter(study == this_study) %>% + dplyr::distinct(sample_id) %>% + dplyr::pull(sample_id) + sample_ids <- if (is.null(sample_ids)) study_sample_ids + else intersect(sample_ids, study_sample_ids) + } + + # A variant that's independently produced by more than one pipeline (e.g. a + # cohort's own published maf and our SLMS-3 recall both calling the same + # position) only has room for one Pipeline value on its single maf + # row -- write_mutations_db() picks one arbitrarily when it deduplicates. + # variant_pipeline (mutation_id, Pipeline) records every pipeline + # that actually produced each variant, keyed on the surrogate mutation_id + # surviving rows are assigned at write time (not the 5-column natural key + # -- see write_mutations_db.R for why), so tool_name is resolved against + # it via a single-column semi-join instead of the maf row's own + # (possibly-not-representative) Pipeline value. Falls back to the old + # direct-column filter against DBs built before this table existed. + has_variant_pipeline <- !is.null(tn) && "variant_pipeline" %in% DBI::dbListTables(con) + # elem distinguished maf's and ashm's independent mutation_id spaces before + # the two tables were merged into one -- only present on a DB built before + # that merge, kept here purely for backward compatibility with one. + has_vp_elem <- has_variant_pipeline && "elem" %in% DBI::dbListFields(con, "variant_pipeline") + + # apply the filters that are common to every query against this table, + # regardless of how many regions (if any) are requested + base_query <- function(table_name) { + q <- dplyr::tbl(con, table_name) %>% + dplyr::filter(genome_build == projection) + if (!is.null(tn)) { + if (has_variant_pipeline) { + matching_ids <- dplyr::tbl(con, "variant_pipeline") + matching_ids <- if (has_vp_elem) { + dplyr::filter(matching_ids, elem == table_name, Pipeline == tn) + } else { + dplyr::filter(matching_ids, Pipeline == tn) + } + matching_ids <- dplyr::distinct(matching_ids, mutation_id) + q <- dplyr::semi_join(q, matching_ids, by = "mutation_id") + } else { + q <- dplyr::filter(q, Pipeline == tn) + } + } + if (coding_only) q <- dplyr::filter(q, Variant_Classification %in% cc) + if (min_read_support > 0) q <- dplyr::filter(q, t_alt_count >= min_read_support) + if (!is.null(sample_ids)) q <- dplyr::filter(q, Tumor_Sample_Barcode %in% sample_ids) + if (!is.null(seq_type) && "maf_seq_type" %in% DBI::dbListFields(con, table_name)) { + q <- dplyr::filter(q, maf_seq_type %in% seq_type) + } + q + } + + # All regions are combined into a single OR'd condition and applied in one + # query, not one query per region. get_ssm_by_regions() commonly passes a + # region per gene in a panel (150-250+ for the full lymphoma gene list) -- + # looping here used to mean that many separate round trips, each one + # redoing the tool_name/coding_only/sample_ids filters and re-running the + # variant_pipeline semi-join from scratch for every single region. A single + # combined query also fixes a latent duplicate-row bug the loop had: 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 two identical rows; SQL's OR doesn't double-count a + # row just because more than one disjunct matches it. + gather <- function(table_name) { + q <- base_query(table_name) + if (!is.null(regions) && nrow(regions) > 0) { + region_exprs <- lapply(seq_len(nrow(regions)), function(i) { + rc <- as.character(regions$chrom[i]) + rs <- as.numeric(regions$start[i]) + re <- as.numeric(regions$end[i]) + rlang::expr(Chromosome == !!rc & Start_Position > !!rs & Start_Position < !!re) + }) + combined <- Reduce(function(a, b) rlang::expr(!!a | !!b), region_exprs) + q <- dplyr::filter(q, !!combined) + } + dplyr::collect(q) + } + + res <- gather("maf") + res$genome_build <- NULL + res$mutation_id <- NULL + res +} diff --git a/R/get_sv_from_db.R b/R/get_sv_from_db.R new file mode 100644 index 0000000..2b2544b --- /dev/null +++ b/R/get_sv_from_db.R @@ -0,0 +1,41 @@ +#' @title Query structural variants (SV) from the mutations database. +#' +#' @description Data-access helper that retrieves Manta SV breakpoints (BEDPE +#' format) from the `bedpe` table of `gambl_mutations.db`, pushing the common +#' filters down to SQL. Backs the GAMBLR.open `get_manta_sv` accessor. +#' +#' @param projection Genome build / `genome_build` column value. Default "grch37". +#' @param sample_ids Optional character vector of `tumour_sample_id` to keep. +#' @param min_vaf Minimum `VAF_tumour`. Default 0. +#' @param min_score Minimum Manta `SCORE`. Default 0. +#' @param pass_only When TRUE, keep only `FILTER == "PASS"` breakpoints. +#' @param region Optional list with `chrom`, `start`, `end`; keeps breakpoints +#' whose A- or B-end falls in the range (positions compared numerically). +#' @param con Optional DBI connection (defaults to [gambl_mutations_db()]). +#' +#' @return A data frame of BEDPE rows (the `genome_build` helper column is dropped). +#' +#' @export +get_sv_from_db <- function(projection = "grch37", + sample_ids = NULL, + min_vaf = 0, + min_score = 0, + pass_only = FALSE, + region = NULL, + con = NULL) { + if (is.null(con)) con <- gambl_mutations_db() + q <- dplyr::tbl(con, "bedpe") %>% + dplyr::filter(genome_build == projection, + VAF_tumour >= min_vaf, + SCORE >= min_score) + if (!is.null(sample_ids)) q <- dplyr::filter(q, tumour_sample_id %in% sample_ids) + if (pass_only) q <- dplyr::filter(q, FILTER == "PASS") + if (!is.null(region)) { + rc <- region$chrom; rs <- as.numeric(region$start); re <- as.numeric(region$end) + q <- dplyr::filter(q, (CHROM_A == rc & START_A >= rs & START_A <= re) | + (CHROM_B == rc & START_B >= rs & START_B <= re)) + } + res <- dplyr::collect(q) + res$genome_build <- NULL + res +} diff --git a/README.md b/README.md index d2d051f..7742701 100644 --- a/README.md +++ b/README.md @@ -15,3 +15,17 @@ Lab and without GSC access. For detailed information about how to install, bundled resources, functionality overview, tutorials, frequently asked questions and much more please visit the website [morinlab.github.io/GAMBLR.data/](https://morinlab.github.io/GAMBLR.data/) + +## Bundled databases + +Curated and sample-level data are served from two SQLite databases rather than +per-object `.rda` files. See `?gambl_reference_db` and `?gambl_mutations_db` +for the full table reference (columns, join keys, indexes). + +- **`gambl_reference_db()`** — small, bundled with the package. Curated gene + lists and aSHM regions (from LLMPP); backs [get_genes()]. +- **`gambl_mutations_db()`** — large (several hundred MB), distributed as a + cached release asset rather than shipped in the package. MAF, aSHM-region + mutation, copy-number segment, and structural-variant calls, plus sample + metadata; backs the GAMBLR.open `get_ssm_*`, `get_cn_segments*`, and + `get_manta_sv` accessors. diff --git a/data-raw/PR_verification_notes.md b/data-raw/PR_verification_notes.md new file mode 100644 index 0000000..bde22dd --- /dev/null +++ b/data-raw/PR_verification_notes.md @@ -0,0 +1,499 @@ +# Verification notes: sample↔study join table / deduplicated SLMS-3 pulls + +Working notes for the PR description covering the `assemble_bundled_data.R` / +`write_mutations_db.R` / `get_ssm_from_db.R` refactor (branch `rmorin-dev`, +commits `1a8cfb6`, `017d3db`). Append further findings (e.g. +`compare_bundle_changes.R` output) below as they come in. + +## Background + +Root cause: SLMS-3 mutation calls were pulled via separate, cohort-specific +code blocks, each independently deciding which samples to pull for. A sample +belonging to more than one cohort (e.g. the Dreval FL × Hilton trios overlap) +got its mutations pulled and bound multiple times. Confirmed 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. + +Fix: `Study`/cohort membership moved off `maf`/`ashm` entirely into a new +`sample_study` many-to-many join table; SLMS-3 is now pulled exactly once per +deduplicated sample regardless of cohort membership. 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. + +## `test_gambl_db.R` run — 2026-07-14 + +All structural/integrity checks passed: + +``` +== tables present == 7/7 PASS (maf, ashm, seg, bedpe, sample_meta, sample_study, build_info) +== row counts (all > 0) == 6/6 PASS + maf 543,214 rows + ashm 134,826 rows + seg 128,075 rows + bedpe 2,352 rows + sample_meta 3,346 rows + sample_study 2,880 rows +== both genome builds present == 4/4 PASS +== required columns present == 5/5 PASS +== indexes present == 7/7 PASS (incl. new idx_study_sample, idx_study_study) +== pipelines present == 2/2 PASS -- Pipeline values: Publication, SLMS-3 (no more "strelka") +== spot query: MYC locus, grch37, slms-3 == PASS -- 4,549 variants found +== build_info counts match actual table counts == 4/4 PASS +== sample_metadata.rda matches sample_meta table == PASS (3,346 == 3,346) +``` + +### Regression vs. old `sample_data.rda` (informational, all 4 differences expected) + +| Table | Old (`sample_data.rda`) | New (`gambl_mutations.db`) | Δ | Explanation | +| --- | --- | --- | --- | --- | +| `maf` | 3,663,468 | 543,214 | −85% | Almost entirely Arthur's unfiltered `strelka` dump being removed (2,827,527 of the old rows, 65,121 distinct genes vs. the ~150-237 gene panel) + the Dreval×Hilton dedup. This is the fix working as intended, not a regression. | +| `ashm` | 134,868 | 134,826 | −0.03% | Noise-level. | +| `seg` | 125,982 | 128,075 | +1.7% | Unrelated to this refactor (CNV/seg blocks untouched) -- likely normal cohort growth in GAMBL's live metadata since the old snapshot. | +| `bedpe` | 941 | 2,352 | +150% | Bonus fix, not explicitly planned: Arthur and Hilton were previously added to `sample_data$meta` *after* the Manta SV block ran, so they got zero SV coverage -- same root cause as the aSHM gap that *was* explicitly fixed. Now that metadata assembly happens once, up front, they get SV coverage too. | + +## Still to do + +- [ ] Run `compare_bundle_changes.R` for the granular per-sample gained/lost/retained view (append results below). +- [ ] Confirm Dreval×Hilton overlap samples (`05-32150T`, `08-15460T`, `09-33003T`, `15-13383T`, `17-36275T`) show 2 `sample_study` rows each and roughly halved mutation counts vs. the old bundle. +- [ ] Confirm Arthur's `maf` rows are gene-panel-restricted (~150-237 distinct genes, not 65,121) and now have hg38 + aSHM coverage. +- [x] Functional check: `get_ssm_by_samples()`/`get_ssm_by_patients()` against real Reddy (capture) data -- 20,480 rows returned, confirming the query layer and underlying data are correct end to end (see note below; the initial zero-result was unrelated to this refactor). + +## Non-issue: `get_ssm_by_samples()`/`get_ssm_by_patients()` initially returned 0 rows for Reddy + +Root cause turned out to be unrelated to this refactor: both functions +default `this_seq_type = "genome"` and filter `these_samples_metadata` by +it *in addition to* whatever metadata you pass in, rather than inferring it. +Reddy is a capture cohort (`seq_type == "capture"` for all 999 samples), so +the default silently dropped every row before the query ever reached the +database. Confirmed the data itself was fine throughout: a raw SQL query +against `maf` for Reddy's `Tumor_Sample_Barcode` values returned 47,626 +rows the whole time. Passing `this_seq_type = "capture"` explicitly fixed +both functions (20,480 rows each, matching after `get_ssm_from_db()`'s +Pipeline/genome_build/read-support filters narrow the raw 47,626). + +Separate, pre-existing GAMBLR.open UX gap worth considering independently +of this PR: neither function warns when the `seq_type` filter drops every +row, so this fails silently rather than with a clear message. + +## Fix: Arthur's Case ID -> patient_id matching, and a `!grepl("tumor", ...)` filter that dropped multi-sample patients entirely + +`compare_bundle_changes.R` showed some patients losing 100% of their SNVs +with zero retention (e.g. patient `08-15460`: 19,897 lost, 0 retained -- +found via `egrep "Tumor_Sa|05-32150|08-15460" snv_persample_counts.tsv`). +Confirmed as an Arthur-specific bug, two compounding causes: + +1. `arthur_meta` matched Arthur's own "Case ID" (from `DLBCL_Arthur.xlsx`) + against GAMBL's `patient_id` via `patient_id %in% arthur_case_ids$\`Case ID\``. + Same class of bug as the `study_id` chr/dbl mismatch fixed earlier: + `read_xlsx()` can silently type a leading-zero ID like `08-15460` as + numeric, and a numeric-vs-character `%in%` comparison fails to match + without any error or warning. Fixed by coercing both sides to character + and using an explicit `inner_join()` instead of a `%in%` filter. +2. The block also had `! grepl("tumor", sample_id)`, which excludes every + sample for any patient with more than one tumor biopsy (e.g. + `..._tumorA`/`..._tumorB` -- the same multi-sample-per-patient pattern + Hilton has). For a patient whose *only* GAMBL samples are tumor-suffixed, + this filter matched nothing, dropping their SLMS-3 coverage entirely. No + other cohort in this script excludes samples this way, and no comment + explained the original intent -- removed. + +**Update**: rebuilt with this fix and re-ran -- `08-15460` (bare) still +showed the identical 100%-loss numbers, byte-for-byte. That ruled out +`arthur_meta`/`arthur_case_ids` as the source: those only ever fed +`sample_data$meta`, never `sample_data$maf` directly. Traced the real +source instead by reading the *old raw flat file itself*: +```r +arthur_raw <- readr::read_tsv("inst/extdata/studies/DLBCL_Arthur.maf.gz", + col_select = "Tumor_Sample_Barcode") +grep("15460", unique(arthur_raw$Tumor_Sample_Barcode), value = TRUE) +# [1] "08-15460" +``` +Confirmed: the deleted raw dump used bare patient-style IDs for +`Tumor_Sample_Barcode` (not GAMBL's real per-sample naming), directly and +independently of `arthur_meta`. **Conclusion: this specific finding was +never a bug.** Every "100%-loss" row keyed to `08-15460` (bare) is simply +the expected disappearance of the old file's non-standard ID convention, +now correctly superseded by GAMBL's real sample IDs (`08-15460T`: 1,515 +retained; `08-15460_tumorB`: 665 retained -- both present and correctly +sized in both old and new bundles). + +Generalized this into a defensive fix rather than leaving it Arthur-specific: +added `relabel_to_sample_id(df, study_name)`, applied to all four +Publication-pipeline pulls (Thomas BL, Thomas DLBCL, Dreval, Reddy's +original-variants file). It 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/Hilton (`study_id` already mirrors `sample_id` +there) and an actual fix for Reddy (`study_id` is the paper's own raw +"Sample ID", genuinely distinct from `sample_id`). Verified against +synthetic data: relabels when a study_id match exists, no-ops when +`Tumor_Sample_Barcode` already is the real sample_id, and passes through +unmatched values unchanged rather than dropping them. + +Still need to: rebuild with this fix and re-run `compare_bundle_changes.R` +to confirm `08-15460` (bare) is gone from the loss list entirely, and +check whether the same bare-vs-suffixed pattern explains the +`DO52686`/`07-35482`-style samples from the earlier, larger (3.1M row) +loss count. + +## Fix: consolidated Phase 4/5 pulls passed sample_data$meta directly to get_ssm_by_regions() -- missing coding calls for hand-built cohorts + +`GAMBLR_examples_output.log`'s regenerated run showed a `get_coding_ssm()` +example for 183 Dreval FL patients (sample_id matching `"SP"`) return **zero +rows**, when it used to return real coding-classified data. Traced through +several rounds of elimination (ID matching confirmed fine, `t_alt_count` +confirmed numeric, `min_read_support` confirmed innocent) to: `coding_only += TRUE` alone reduced 28,867 real SLMS-3/grch37 rows for these patients to +zero -- every one of them was non-coding (`Intron`/`5'Flank`/`3'Flank`/ +`5'UTR`/`3'UTR`). For 183 real FL patients to have literally zero coding +mutations across the whole lymphoma gene panel is biologically implausible. + +Root cause: `unix_group` (and likely other native GAMBL metadata columns) +is `NA` for every cohort built by hand in this script rather than pulled +live from `get_gambl_metadata()` -- confirmed directly: +```r +filter(sample_data$meta, is.na(unix_group)) %>% count(study) +# BL_Thomas 234, DLBCL_Arthur 153, DLBCL_Thomas 43, DLBCL_cell_lines 5, FL_Dreval 443 +filter(sample_data$meta, !is.na(unix_group)) %>% count(study) +# DLBCL_Hilton 159, NCI_DLBCL_Golub 124, dlbcl_chapuy 233, dlbcl_reddy 999, dlbcl_schmitz 951 +``` +`get_ssm_by_regions()` apparently relies on columns like `unix_group` to +determine which underlying merged flat file to read for a sample -- with it +missing, it was silently returning an incomplete (non-coding-only) subset +for affected samples instead of erroring. This was specific to Phase 4 (the +consolidated SLMS-3 pull) and Phase 5 (the 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 (the SV/bedpe block, the cell-lines pull, the Publication-pipeline +enrichment steps) already re-fetched fresh, complete metadata via +`get_gambl_metadata()` filtered by `sample_id` first, which is exactly why +this only surfaced in the two things newly introduced this session. + +Fixed both: `sample_data$meta` is now used only to determine which +`sample_id`s 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. General principle going +forward, not just for these two spots: never pass `sample_data$meta` (or +anything derived from it) directly to a GAMBLR.results function. + +**Update**: rebuilt and re-ran `compare_bundle_changes.R` broadly (not just +the one example). Confirmed fixed -- samples that previously showed +catastrophic 0-retained losses now show large, healthy retained counts +(`SP116715`: 2,901 retained; `16-11636T`: 1,624; `14-20552_tumorB`: 670; +`04-24937T`: 990), with only small gains/losses on top (single digits to +low teens). Total SNV loss across the whole bundle dropped from 3.1M rows +to 299, spread across 113 samples (~2.6 rows/sample average) -- a +different order of magnitude, consistent with normal residual noise rather +than a bug. + +## Current state as of this point in the session (superseded -- see below) + +- **SNV**: 348 samples w/ gains, 113 samples w/ losses, 299 total lost rows + -- converged to normal residual noise. The gained-only, 0-retained + samples at the top of `snv_persample_counts.tsv` are the known "left out + of the last release" bonus fix (unrelated, already understood). +- **CNV (seg)**: unchanged from prior runs (5 cell lines w/ gains, 0 w/ + losses) -- expected, nothing in this session touched seg assembly. +- **SV (bedpe)**: unchanged from prior runs (433 w/ gains, 7 w/ losses) -- + expected, consistent with the earlier Arthur/Hilton SV coverage bonus + fixes already being stable. + +`compare_bundle_changes.R`'s row-count-based diffing was genuinely clean at +this point -- but it can't detect a variant surviving with the *wrong* +Pipeline tag (the row is still there, at the same position, just labeled +differently), which is exactly the bug found afterward (see "Fix: +write-time dedup silently relabeled real SLMS-3 calls as 'Publication'" +below). Do not treat a clean `compare_bundle_changes.R` run alone as +sufficient sign-off for Pipeline-level correctness -- cross-check with a +`tool_name`-filtered query (e.g. `get_all_coding_ssm()` or +`get_ssm_from_db(tool_name=...)`) per cohort as well. + +## Fix: cell lines need a separate, genome-wide SNV pull + +After the strelka-exclusion fix, `compare_bundle_changes.R` came back clean +except for the 5 cell lines (`SU-DHL-4`, `OCI-Ly3`, `OCI-Ly10`, `SU-DHL-10`, +`DOHH-2`), each showing ~43k-65k lost rows against a few hundred/~1k +retained. This is the scope-reduction flagged during planning, now +confirmed as something to actually fix rather than accept: the original +pre-refactor cell-line pull used `get_ssm_by_samples()` with no gene-panel +filter (genome-wide), but folding cell lines into the consolidated, +panel-restricted SLMS-3 pull (Phase 4) silently reduced them to panel-only +coverage. + +Fixed by excluding cell lines from `all_slms3_meta` (so they're not also +pulled panel-restricted) and adding back a dedicated, unrestricted +`get_ssm_by_samples()` pull for them specifically, both builds, bound into +`slms3_grch37`/`slms3_hg38` alongside the panel-restricted data for +everyone else. Still tagged `Pipeline = "SLMS-3"` (same underlying pipeline, +just unrestricted scope for these 5 samples). + +**Update**: rebuilt and re-ran -- cell-line losses are gone, confirming the +fix worked. One more thing surfaced in the same run, another bonus rather +than a bug: the top of `snv_persample_counts.tsv` now shows several samples +(e.g. `15-18723_CLC02138`, `15-43891_tumorB`, `07-17613T`) with real gained +counts (64-630) but `n_lost = 0` and `n_retained = 0` -- meaning they had +*zero* rows in the old bundle at all, not a mismatch. Confirmed these were +samples left out of the last release for some reason and are now correctly +included -- consistent with the pattern already seen for Arthur/Hilton +(aSHM, SV coverage) where the old per-cohort-block gating silently excluded +samples that the consolidated pull now correctly picks up. + +## Fix: `study_id` used `patient_id` instead of `sample_id` for Thomas/Dreval/Hilton + +Caught in review (not by the automated checks above): the first pass at +populating `sample_study$study_id` used `patient_id` for Thomas BL, Thomas +DLBCL, Dreval, and Hilton. For Hilton specifically this is a real bug, not +just imprecision -- Hilton is a trios study, so a single patient can have +multiple samples (e.g. `LY_RELY_116_tumorA` and `LY_RELY_116_tumorB`), and +`patient_id` collapsed both to the same `study_id`, making the two +`sample_study` rows ambiguous. Thomas/Dreval don't currently have that +ambiguity in practice, but the same fix applies for consistency: all four +now use `sample_id` (already sourced from each xlsx's own "Genome sample +id"/`DNAseq_sample_id` column, i.e. already the closest thing to a +study-native sample-level identifier) instead of `patient_id`. Arthur and +Reddy were left unchanged -- Arthur's block already excludes +`tumor`-suffixed sample_ids per patient (structurally safe from the same +ambiguity), and Reddy's `study_id` comes from the paper's own distinct +"Sample ID" column, already 1:1 with `sample_id`. + +## Fix: `this_study_samples`/`these_samples`/`these_samples_dlbcl` filtered on `cohort` instead of `study` + +Found via a fresh GSC test (`get_all_coding_ssm()` against the built DB, +joined to metadata): `FL_Dreval` showed only 7 of 441 samples with any +coding-classified variant at all, each with 1-2 rows; `DLBCL_cell_lines` +looked similarly too low. Traced to the "proteinpainter compatibility" +block (Phase 3), which intentionally reads the *installed, stale* +`GAMBLR.data::sample_data$meta` rather than the locally-built one (a +pre-existing, documented, out-of-scope inconsistency) -- but filtered that +stale snapshot on a column called `cohort`, when this script's own +convention (and, empirically, the stale snapshot itself) uses `study`. +`GAMBLR.data::sample_data$meta %>% filter(cohort %in% c("FL_Dreval", +"DLBCL_cell_lines"))` returned only the 5 cell-line ids; switching to +`filter(study %in% ...)` returned the correct 400+ ids. Fixed at all three +call sites (`these_samples`, `these_samples_dlbcl`, `this_study_samples`). + +Mechanistically this bug only affected the RefSeq/Protein_position +enrichment `left_join()` for Dreval/cell-line Publication rows (a +non-match there just leaves those two columns `NA`, it can't drop rows) -- +so on its own it was not expected to explain the *row-count* collapse +`get_all_coding_ssm()` showed. Confirmed by rebuilding: see below. + +## Real-time build diagnostics added to `assemble_bundled_data.R` + +To stop diagnosing via a slow guess-rebuild-test loop against the GSC, added +`diag_summary()` / `diag_summary_maf()` / `diag_missing()` helper calls at +every major checkpoint in the script (every Phase 1 cohort, the +cohort-\>study rename/rejoin block, `sample_study` assembly, the Phase 3 +enrichment joins, the Phase 4 `all_slms3_meta` live-metadata refetch and +pull outputs, the Phase 5 aSHM pull, and the final per-study +coding-classified `maf` counts). Prints per-study row/sample counts and, +at the two live-metadata refetch points, exactly which sample_ids get +silently dropped (if any) -- the log now shows where a cohort's count +changes during the build itself, not just after the fact. + +### First full run with diagnostics -- 2026-07-14 + +The `cohort`/`study` fix produced healthy-looking **aggregate** per-study +coding counts in the in-memory `sample_data` object, before the SQLite +write: + +| Study | grch37 maf, coding-classified only | | hg38 maf, coding-classified only | | +| --- | --- | --- | --- | --- | +| | rows | samples | rows | samples | +| Dreval | 39,765 | 441 | 4,982 | 437 | +| Arthur | 17,855 | 159 | 2,084 | 157 | +| Reddy | 17,805 | 981 | 10,816 | 947 | +| Schmitz | 8,352 | 468 | 8,320 | 468 | +| Thomas | 3,577 | 276 | 46,207 | 277 | +| Chapuy | 3,226 | 230 | 3,189 | 230 | +| Hilton | 3,078 | 158 | 2,607 | 158 | +| NCI_Golub | 2,030 | 114 | 2,017 | 114 | +| cell lines (no `sample_study` row, shown as `NA`) | 1,211 | 5 | 1,214 | 5 | + +Also visible in this run, two small pre-existing gaps (not blocking, not +part of this refactor's scope): +- `Reddy meta -> reddy_meta_gambl join`: 2 of 1001 Reddy sample_ids + (`Reddy_832T`, `Reddy_3813T`) have no matching live cohort. +- `Hilton trios`: 1 of 160 (`14-27873_tumorA`) not found in live metadata. + +**This table turned out to mask a real, separate problem** -- see below. +It aggregates Publication + SLMS-3 pipeline rows together, and Dreval's +Publication-pipeline pull (the paper's own curated maf) dominates the +total, so the healthy-looking aggregate hid the fact that the *SLMS-3* +pipeline's own contribution was nearly empty. + +## Fix: write-time dedup silently relabeled real SLMS-3 calls as "Publication" + +A follow-up GSC test using `get_all_coding_ssm()` (which specifically +restricts to `Pipeline == "slms-3"`, unlike the aggregate table above) +still showed `FL_Dreval` at 7/441 samples with any coding variant, each +with 1-2 rows -- i.e. the `cohort`/`study` fix above did not actually +resolve the originally-reported symptom; it only fixed a separate, +real-but-narrower bug (broken RefSeq/Protein_position enrichment). + +Extensive isolation via direct SQL against the built DB and side-by-side +`GAMBLR.results::get_ssm_by_regions()` calls ruled out, in order: a stale +DB path/cache (confirmed same file, fresh mtime), `tool_name` +case-sensitivity (already handled via `tolower()` on both sides), +`coding_class` casing/definition mismatches (byte-identical, confirmed via +`LENGTH()`/`HEX()` in SQL), an oversized `sample_ids` `IN`-clause hitting a +SQLite variable limit (reproduced with just Dreval's own 441 ids), a +scale/multi-cohort-merge bug in `get_ssm_by_regions()` (incrementally +rebuilt the exact Phase 4 sample scope -- Dreval alone, +Hilton, +Thomas, ++Arthur -- and every combination pulled healthy coding calls interactively), +and the Hugo_Symbol post-filter dropping bystander genes (a single sample's +raw coding calls were all confirmed to be on canonical panel genes and all +survived the filter). + +Root cause, confirmed empirically: `sample_data$grch37$maf <- +bind_rows(grch37_publication_rows, slms3_grch37)` puts Publication rows +first, and `write_mutations_db()`'s dedup +(`distinct(..., .keep_all = TRUE)`, keyed on `Tumor_Sample_Barcode, +Chromosome, Start_Position, End_Position, Tumor_Seq_Allele2`) keeps +whichever row comes first for a given key. Dreval's own published maf +(`fl_data$ssm_to_bundle`) was itself originally SLMS-3-called, so a real +mutation legitimately appears in *both* `grch37_publication_rows` and +`slms3_grch37` with an identical key -- the dedup silently kept the +Publication-tagged copy and dropped the SLMS-3-tagged one every time. The +variant was never lost from the bundle, but it became invisible to +`tool_name = "slms-3"`-style queries. Confirmed directly: of 6,807 real, +undeduplicated SLMS-3 coding calls for Dreval (from a raw +`get_ssm_by_regions()` pull), 6,791 (99.8%) have an exact key match among +the DB's `Pipeline = "Publication"` rows. + +### Fix: `variant_pipeline` join table + +Same architectural pattern as `sample_study`: a single-valued `Pipeline` +column on `maf`/`ashm` can't represent a variant independently called by +more than one pipeline, any more than a single-valued `study` column could +represent multi-cohort sample membership. Added `variant_pipeline` +(`Tumor_Sample_Barcode, Chromosome, Start_Position, End_Position, +Tumor_Seq_Allele2, genome_build, elem, Pipeline` -- one row per +(variant, Pipeline) pair, `elem` distinguishing `maf`/`ashm`), captured in +`write_mutations_db()` from the full pre-dedup data before the existing +dedup collapses each variant to one representative row. `get_ssm_from_db()`'s +`tool_name` filter now resolves via a semi-join against `variant_pipeline` +instead of `maf`/`ashm`'s own `Pipeline` column, with a fallback to the old +direct-column filter for DBs built before this table existed. + +Verified via synthetic data: a variant called by both "Publication" and +"SLMS-3" pipelines (identical key) is correctly deduplicated to one +`maf` row but produces two `variant_pipeline` rows; `get_ssm_from_db(tool_name += "slms-3")` correctly finds it even though its surviving `maf` row is +tagged "Publication"; `tool_name = "publication"` also finds it; and with +`variant_pipeline` dropped (simulating an old DB), the query correctly +falls back to the old (narrower, expected) behavior rather than erroring. + +**Update**: rebuilt and confirmed on the GSC. `get_all_coding_ssm()` and +direct SQL against the fresh DB both show healthy Dreval SLMS-3-tagged +coding counts through the `variant_pipeline` join, resolving the original +symptom. + +## Performance: query-layer indexing and GAMBLR.open query patterns + +While retesting the fix above, `GAMBLR.open`'s `tools/logExampleOutputs.R` +example suite runtime came back noticeably regressed against its +pre-refactor baseline: 51.2s (`b061311`, 2026-06-12, pre-`gambl_mutations.db` +sample_study/variant_pipeline work) vs. 4m37s-11m35s across several +mid-refactor reruns. Investigated and fixed several distinct, compounding +issues: + +- **`Pipeline` case mismatch broke every index on it.** `idx_maf_pipe`/ + `idx_vp_pipe` were plain 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 b-tree index to + satisfy a function-wrapped predicate, so these indexes were never + actually usable; every `tool_name`-filtered query (the default) 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 both `maf`/`ashm` and `variant_pipeline`, so reads use a + plain equality against a plain index. +- **No index at all on the new `variant_pipeline` semi-join's key.** Added + composite indexes on the variant key to both sides of the join. +- **`variant_pipeline`'s 5-column natural key replaced with a surrogate + `mutation_id`** (see "Fix: `variant_pipeline` join table" above) -- a + single indexed integer join instead of a 5-column composite (two of them + text). +- **`bedpe` had only one index** (`tumour_sample_id` alone) despite every + `get_sv_from_db()`/`get_manta_sv()` call also filtering `genome_build` + and its region search checking both breakpoint ends. Added + `(tumour_sample_id, genome_build)` and a `genome_build`-led index per + breakpoint end. +- **`get_ssm_from_db()`'s `regions` parameter looped, one query per + region.** `get_ssm_by_regions()` commonly passes 100+ regions (one per + gene in the aSHM/lymphoma panel) in a single call -- consolidated into + one query with a combined OR'd condition instead of N round trips. Also + fixed a latent duplicate-row bug this had: a variant inside two + overlapping regions was previously returned once per matching region and + `bind_rows()`'d into duplicate rows. +- **`GAMBLR.open::calc_mutation_frequency_bin_regions()` looped over every + region calling a function that queries the database per region** (not a + `gambl_mutations.db`-side issue, but the single largest cost found: 142s + of a ~240s total example-suite run). Its own per-region helper already + supported a pre-fetched `maf_data` short-circuit that skipped the + database entirely; the plural function just never used it. Fixed to + bulk-fetch once via `get_ssm_by_regions()` before the per-region loop. + Also replaced `mclapply()` with `lapply()` here, at the user's request -- + once each iteration is fast in-memory subsetting rather than a database + call, there's little left to parallelize, and `mclapply()`'s + fork-per-worker model both multiplies peak memory (each fork copies the + parent's memory) and forking a live DB connection across processes is a + correctness/contention risk this avoids entirely. + +All fixes verified with synthetic-data/stubbed-out control-flow tests +(see individual commits) rather than against the live GSC, since none of +this repo's automated tests have GSC access. + +### Results + +Measured with `tools/fast_example_timer.R` (GAMBLR.open), a leaner +per-example-timed replacement for `logExampleOutputs.R` written during +this investigation -- `devtools::run_examples()` has `document = TRUE` by +default, meaning every run also pays for a full roxygen regeneration of +every `.Rd` file + `NAMESPACE` plus two full package reloads, unrelated to +anything under test here. + +| Stage | Total suite time | `calc_mutation_frequency_bin_regions` | +| --- | --- | --- | +| Mid-refactor, before indexing fixes | ~242s (4m22s incl. `run_examples()` overhead) | 142s | +| After `calc_mutation_frequency_bin_regions` bulk-fetch + `lapply` fix (indexing fixes already included), on the GSC | ~130s | 30.2s | +| Same DB synced to a local laptop, same code, local disk instead of GSC storage | 33.8s (incl. 5.8s package load) | 12.2s | + +On the GSC, the next-largest remaining items after the fix above were +`assign_cn_to_ssm` (23.7s) and `get_ashm_count_matrix` (19.2s). Investigated +both -- neither has the same "N queries instead of 1" bug already fixed +above; both already make single, consolidated database calls +(`get_cn_segments()` + `get_ssm_by_samples()`, and `get_ssm_by_regions()` +respectively). Initially assessed `assign_cn_to_ssm`'s cost as inherent +R-side computation (its `cool_overlaps()` call is a chromosome-level +many-to-many join between mutations and CN segments, scaling with +mutations x segments per chromosome per sample) rather than something an +index could fix. + +**That assessment was wrong, or at least incomplete** -- the laptop +comparison above shows `assign_cn_to_ssm` dropping out of the slowest-10 +list entirely (under ~1.1s) on the same DB, same code, local disk. A +`cool_overlaps()` join being purely CPU-bound R computation wouldn't +plausibly speed up ~20x just from moving the DB file to local storage; the +dominant cost on the GSC was disk/network I/O for its two underlying +queries, not the join itself. `calc_mutation_frequency_bin_regions` shows +the same pattern to a lesser degree (30.2s -> 12.2s on the identical single +bulk query + windowing code) -- some of its remaining cost is I/O-bound +too, not purely the sliding-window computation. + +Practical implication: further **code**-level query optimization on the +GSC has sharply diminishing returns at this point. The larger remaining +lever is infrastructure -- where `gambl_mutations.db` physically lives on +GSC storage (e.g. local scratch/SSD vs. a networked home-directory mount) +-- not something fixable from this refactor. + +Not a like-for-like comparison to the original 51.2s baseline: that +predates `gambl_mutations.db` entirely (in-memory, lazy-loaded +`sample_data.rda`, not a SQL database), and cohort sizes have grown +materially since then regardless (e.g. `DLBCL_Gascoyne` 21->42, +`DLBCL_GenomeCanada` 59->110). A SQL-backed system pays real per-query +overhead (connection, parse, plan, result marshaling) on every one of the +dozens of separate calls across the example suite that an in-memory R +object never had to pay; indexing closes the "full scan vs. indexed +lookup" gap, not that architectural difference. diff --git a/data-raw/RELEASING.md b/data-raw/RELEASING.md new file mode 100644 index 0000000..d2f586f --- /dev/null +++ b/data-raw/RELEASING.md @@ -0,0 +1,81 @@ +# Releasing a new `gambl_mutations.db` + +This is the maintainer-only process for publishing a freshly-built +`gambl_mutations.db` so that `GAMBLR.data::gambl_mutations_db()` can +auto-download it for new installs. It's not part of the main `README.md` +because it's not something package users ever need to do -- only whoever +is building and shipping a new version of the bundled data. + +## Background + +`gambl_mutations.db` is too large (several hundred MB) to ship inside the +package itself. Instead it's distributed as a GitHub Release asset and +downloaded on demand: + +- `data-raw/assemble_bundled_data.R` (+ `data-raw/write_mutations_db.R`) + builds the file from GSC data. +- `data-raw/release_mutations_db.R` (this doc) uploads it to a GitHub + Release. +- `R/download_gambl_mutations_db.R` downloads it back down for any install + that doesn't already have a copy -- called automatically by + `R/gambl_mutations_db.R` the first time it's needed. + +Releases are tagged `data-v`. This is deliberate: +`download_gambl_mutations_db()` resolves the tag to fetch from +`utils::packageVersion("GAMBLR.data")` on the *installing* machine, so a +given install always gets data built for its own schema -- not whatever +the newest release happens to be. This means **a schema change requires +both a version bump and a new release**, not just a new upload. + +## Steps + +1. **Build a fresh `gambl_mutations.db`** via `assemble_bundled_data.R` (on + the GSC). Run `data-raw/test_gambl_db.R` and `data-raw/compare_bundle_changes.R` + against it first -- don't release something that hasn't passed those. + +2. **Bump `DESCRIPTION`'s `Version`.** Required any time the DB schema + changed since the last release (new/changed tables, columns, or + indexes) -- if in doubt, bump it. The release tag is derived from this. + +3. **Commit and push the version bump.** `release_mutations_db.R` will + refuse to run otherwise. The reason: the tag it creates is only + findable by other installs if the version that generated it is what's + actually committed and installed elsewhere -- an uncommitted or + unpushed bump would create a release nothing can ever resolve to, + including your own next reinstall of this checkout. + +4. **Install `piggyback`** if you haven't already: `install.packages("piggyback")`. + +5. **Make sure you have a GitHub PAT with write access to this repo.** + Either already configured (`gitcreds::gitcreds_set()`) or set via the + `GITHUB_PAT`/`GITHUB_TOKEN` environment variable. This is the step most + likely to trip you up if it's not already set -- `piggyback::pb_upload()` + fails with an auth error without it, not a helpful one. + +6. **Run the release script from the repo root:** + ```bash + Rscript data-raw/release_mutations_db.R /path/to/gambl_mutations.db + ``` + This creates the `data-v` release/tag if it doesn't already + exist and uploads the file to it. There's no dry-run mode -- double + check the path argument before running, since this creates a real, + public release with a several-hundred-MB asset attached. + +7. **Verify the download side works** by clearing your local cache + (`unlink(file.path(tools::R_user_dir("GAMBLR.data", "cache")), recursive = TRUE)`) + and calling `GAMBLR.data::gambl_mutations_db()` fresh -- it should + report downloading, then connect successfully. + +## Troubleshooting + +- **"DESCRIPTION has uncommitted changes"** -- commit and push the version + bump first (step 3). +- **Auth error from `piggyback::pb_upload()`/`pb_new_release()`** -- your + GitHub PAT isn't configured, or doesn't have write access to this repo + (step 5). +- **A user reports a download that 404s or an outdated schema** -- check + that their installed `packageVersion("GAMBLR.data")` actually matches a + tag that exists (`piggyback::pb_releases(repo = "morinlab/GAMBLR.data")`); + usually means either the release for their version was never cut, or + they installed from a commit whose `DESCRIPTION` version doesn't have a + matching release yet. diff --git a/data-raw/assemble_bundled_data.R b/data-raw/assemble_bundled_data.R index 597a7a9..a1fba17 100644 --- a/data-raw/assemble_bundled_data.R +++ b/data-raw/assemble_bundled_data.R @@ -1,10 +1,36 @@ -setwd("~/my_dir/repos/GAMBLR.data/") +# --- Working directory & config ------------------------------------------- +# Run from the GAMBLR.data checkout (holds inst/extdata and data/); override +# with the GAMBLR_DATA_ROOT env var. No hard-coded paths. +PKG_ROOT <- Sys.getenv("GAMBLR_DATA_ROOT", unset = getwd()) +stopifnot(dir.exists(PKG_ROOT)) +setwd(PKG_ROOT) library(readxl) library(GAMBLR) library(parallel) library(tidyverse) +# GAMBLR.results locates GSC data via config.yml (repo_base / project_base). +# check_config_and_value() looks first in the working directory, then falls back +# to the copy shipped in GAMBLR.results/inst/extdata. As a last resort (and to +# cover any direct config::get() calls), point R_CONFIG_FILE at that shipped +# config when nothing else provides one. +if (!file.exists("config.yml") && !nzchar(Sys.getenv("R_CONFIG_FILE"))) { + Sys.setenv(R_CONFIG_FILE = system.file("extdata", "config.yml", + package = "GAMBLR.results")) +} + +# get_gambl_metadata()'s min_corrected_cov QC filter (default 15x) can drop +# samples that were previously part of a released bundle out of a fresh one +# (e.g. 01-16433_tumorA/B: FFPE genomes with coverage under the bar, silently +# excluded, taking all their SNVs with them). Disable it for bundle assembly +# so a rebuild doesn't lose samples that are already in the released dataset. +# Shadows the package function for the rest of this script, so every call +# site (there are ~13) picks this up in one place. +get_gambl_metadata <- function(...) { + GAMBLR.results::get_gambl_metadata(..., min_corrected_cov = 0) +} + # Global variables definition colnames_for_bundled_meta <- c( "patient_id", @@ -35,7 +61,12 @@ pmids <- list( maf_columns_to_keep <- c( "RefSeq", - "Protein_position" + "Protein_position", + # Carries genome/capture provenance through every select(all_of(all_cols))/ + # select(any_of(all_cols)) call below, so a sample_id with both seq_types + # doesn't become indistinguishable once its rows are merged (see + # GAMBLR.results::get_ssm_by_regions(), which now stamps this per-row). + "maf_seq_type" ) all_cols <- c( @@ -43,11 +74,184 @@ all_cols <- c( maf_columns_to_keep ) -# restrict to the most inclusive DLBCL gene list -all_lymphoma_genes <- lymphoma_genes_comprehensive$Gene +# restrict to the most inclusive DLBCL gene list. Expanded to include known +# aliases (e.g. old/new HGNC histone names) so the Hugo_Symbol %in% +# all_lymphoma_genes filters below don't drop rows annotated under a gene's +# other name. +all_lymphoma_genes <- GAMBLR.utils::expand_gene_aliases(lymphoma_genes_comprehensive$Gene) + +# Captured once, up front -- referenced by pull_data()'s diagnostics below. +# GAMBLR.data:::coding_class is a lazy-loaded internal object; reading it on +# every pull_data() call (deep into Phase 4) instead of once here means any +# mid-run change to the installed GAMBLR.data package (e.g. a concurrent +# devtools::install()) can hit a missing/replaced .rdb file well after this +# script started, rather than failing fast at startup if at all. +diag_coding_class <- GAMBLR.data:::coding_class + +# Built once, reused by every get_ssm_by_regions() call below that wants only +# lymphoma-gene mutations. Restricting via tabix -R (region) instead of +# pulling every mutation for a sample set is the performance win; every call +# site still adds `filter(Hugo_Symbol %in% all_lymphoma_genes)` afterward on +# the now-small result, because gene_to_region()'s coordinates and VEP's +# Hugo_Symbol assignment don't always agree at gene boundaries -- without the +# post-filter, mutations from overlapping/neighbouring genes and non-coding +# loci (e.g. AC/AL/AF-prefixed lncRNA transcripts) leak in as false "gains". +# GENE_PAD_BP widens the tabix window itself so real target-gene mutations +# just outside gene_to_region()'s exact span (promoter/UTR/annotation-source +# discrepancies) aren't lost before the Hugo_Symbol filter even sees them; +# the filter makes over-padding cheap (extra I/O, not incorrect inclusion). +# 10kb is a 2x margin over VEP's default 5kb upstream/downstream annotation +# window (the likely source of Hugo_Symbol on these rows) -- e.g. an ID3 +# variant 4.2kb upstream of its TSS (well inside VEP's 5kb default, but +# outside a 2kb pad) was confirmed lost under the old 2000bp value. +GENE_PAD_BP <- 10000 +lymphoma_genes_bed_grch37 <- create_bed_data( + gene_to_region(gene_symbol = all_lymphoma_genes, projection = "grch37", + return_as = "bed", pad_length = GENE_PAD_BP), + genome_build = "grch37" +) +lymphoma_genes_bed_hg38 <- create_bed_data( + gene_to_region(gene_symbol = all_lymphoma_genes, projection = "hg38", + return_as = "bed", pad_length = GENE_PAD_BP), + genome_build = "hg38" +) + +# Wraps an expression, printing its wall-clock time; returns the expression's +# value unchanged. Used below to benchmark get_ssm_by_samples() +# (subset_from_merge TRUE vs FALSE, per-sample loop) against +# get_ssm_by_regions() (tabix -R) at each call site, to decide which +# approach to standardize on for each use case. +time_it <- function(label, expr) { + t0 <- Sys.time() + result <- expr + elapsed <- as.numeric(difftime(Sys.time(), t0, units = "secs")) + message(sprintf("[TIMING] %-55s %8.1fs", label, elapsed)) + result +} +##### Importing SLMS-3 variants +# Reused by the single consolidated SLMS-3 pull in Phase 4 below (and by +# nothing else -- every cohort-specific SLMS-3 pull that used to call this +# once per cohort block has been removed; see Phase 4). +pull_data <- function( + pull_meta, + pull_projection = "grch37" + ){ + lymphoma_genes_bed <- if(pull_projection == "grch37") lymphoma_genes_bed_grch37 else lymphoma_genes_bed_hg38 + + # label identifies this specific call (projection x seq_type) in the log, + # since pull_data() is called 4 times per build (genome/capture x + # grch37/hg38) -- without this, diagnostics from different calls would + # be indistinguishable in the log. + call_label <- sprintf("pull_data(%s, seq_type=%s, n=%d samples)", + pull_projection, + paste(unique(pull_meta$seq_type), collapse = "/"), + n_distinct(pull_meta$sample_id)) + message(sprintf("[DIAG] %s: starting", call_label)) + diag_summary(pull_meta, paste0(call_label, ": input metadata"), study_col = NULL) + + raw <- get_ssm_by_regions( + regions_bed = lymphoma_genes_bed, + these_samples_metadata = pull_meta, + basic_columns = FALSE, + projection = pull_projection + ) + + # Catches exactly the kind of silent, no-error data loss this + # investigation is chasing: if get_ssm_by_regions() only partially read + # an underlying file, readr/vroom's parsing-problem log (normally only + # surfaced as a generic console warning) is captured here per-call + # instead, so a future run can show precisely which rows/columns had + # problems rather than just that "one or more parsing issues" occurred + # somewhere in the whole build. + probs <- tryCatch(readr::problems(raw), error = function(e) NULL) + n_probs <- if (!is.null(probs)) nrow(probs) else NA_integer_ + message(sprintf("[DIAG] %s: raw get_ssm_by_regions() = %d rows, %d distinct Tumor_Sample_Barcode, %s parsing problem(s)", + call_label, nrow(raw), n_distinct(raw$Tumor_Sample_Barcode), + if (is.na(n_probs)) "unknown (problems() not applicable to this object)" else n_probs)) + if (!is.null(probs) && nrow(probs) > 0) { + message(sprintf("[DIAG] %s: parsing problem detail (up to 20 rows):", call_label)) + print(utils::head(probs, 20)) + } + diag_summary_maf(raw, paste0(call_label, ": raw, pre Hugo_Symbol filter, all Variant_Classifications")) + diag_summary_maf(raw %>% filter(Variant_Classification %in% diag_coding_class), + paste0(call_label, ": raw, pre Hugo_Symbol filter, coding-classified only")) + + slms3 <- raw %>% + filter(Hugo_Symbol %in% all_lymphoma_genes) %>% + select( + all_of(all_cols) + ) -# Importing BL data from Thomas et al + diag_summary_maf(slms3, paste0(call_label, ": final, post Hugo_Symbol filter, all Variant_Classifications")) + diag_summary_maf(slms3 %>% filter(Variant_Classification %in% diag_coding_class), + paste0(call_label, ": final, post Hugo_Symbol filter, coding-classified only")) + + return(slms3) +} + + +# --- Diagnostics ------------------------------------------------------------- +# Prints per-study row/sample counts at a labeled checkpoint, so a build's +# log shows exactly where a cohort's count changes (or unexpectedly drops) +# in real time, instead of only being discoverable afterward via a separate +# GSC session against the finished DB. sample_col/study_col differ across +# checkpoints (metadata frames use sample_id + study or cohort; maf-shaped +# frames use Tumor_Sample_Barcode and carry no inline study column at all -- +# use diag_summary_maf for those). +diag_summary <- function(df, label, sample_col = "sample_id", study_col = "study") { + message(sprintf("[DIAG] %s: %d rows, %d distinct %s", label, nrow(df), n_distinct(df[[sample_col]]), sample_col)) + if (!is.null(study_col) && study_col %in% names(df)) { + summary_tbl <- df %>% + group_by(.data[[study_col]]) %>% + summarise(n_rows = n(), n_samples = n_distinct(.data[[sample_col]]), .groups = "drop") %>% + arrange(desc(n_rows)) + print(as.data.frame(summary_tbl)) + } + invisible(df) +} + +# Same idea for maf/ashm-shaped frames: no inline study column, so study is +# looked up via sample_study. A sample in >1 study is counted once per study +# it belongs to (matching how sample_study itself works), so totals across +# studies can exceed nrow(df) -- that's expected, not a bug. +diag_summary_maf <- function(df, label, study_lookup = sample_data$sample_study) { + message(sprintf("[DIAG] %s: %d rows, %d distinct Tumor_Sample_Barcode", label, nrow(df), n_distinct(df$Tumor_Sample_Barcode))) + summary_tbl <- df %>% + left_join(study_lookup %>% select(sample_id, study), by = c("Tumor_Sample_Barcode" = "sample_id"), relationship = "many-to-many") %>% + group_by(study) %>% + summarise(n_rows = n(), n_samples = n_distinct(Tumor_Sample_Barcode), .groups = "drop") %>% + arrange(desc(n_rows)) + print(as.data.frame(summary_tbl)) + invisible(df) +} + +# Prints which of before_ids are missing from after_ids -- for catching +# silent-filter row loss (e.g. a get_gambl_metadata() re-fetch that quietly +# drops sample_ids not present in the live source) at the exact step it +# happens, instead of guessing after the fact. +diag_missing <- function(before_ids, after_ids, label) { + missing <- setdiff(before_ids, after_ids) + message(sprintf("[DIAG] %s: %d of %d input ids missing after this step%s", + label, length(missing), length(before_ids), + if (length(missing) > 0) paste0(" (e.g. ", paste(head(missing, 5), collapse = ", "), ")") else "")) + invisible(missing) +} + + +# ============================================================================ +# Phase 1: metadata (and, where genuinely distinct per-paper data exists, +# Publication-pipeline SSM calls) for every cohort. No SLMS-3 pulling here -- +# that used to happen per-cohort-block, which silently duplicated mutation +# rows whenever a sample belonged to more than one cohort (confirmed: e.g. +# the Dreval FL x Hilton trios overlap). Every cohort below also builds a +# small (sample_id, study, study_id, reference_PMID) frame, assembled into +# the new sample_study table in Phase 2 -- this replaces the old inline +# `Study =` tagging on mutation rows entirely, for every pipeline, so cohort +# membership is tracked once per sample instead of once per mutation row. +# ============================================================================ + +# --- Thomas BL (grch37, genome) -------------------------------------------- # It has more patients and also contains sample ids, not just patient ids bl_data <- list() @@ -88,6 +292,8 @@ bl_data$meta_to_bundle <- read_xlsx( mutate(genetic_subgroup = Subgroup) %>% select(- Subgroup) +diag_summary(bl_data$meta_to_bundle, "bl_data$meta_to_bundle", study_col = "cohort") + bl_data$ssm_to_bundle <- read_xlsx( "inst/extdata/studies/BL_Thomas.xlsx", sheet = 6 @@ -113,8 +319,20 @@ bl_data$cnv_to_bundle <- read_xlsx( filter(ID %in% bl_data$meta_to_bundle$sample_id) %>% mutate(CN = round(2 * 2^log.ratio)) +# study_id mirrors sample_id here (not patient_id): the xlsx's own +# "Genome sample id" column was already adopted directly as sample_id, and +# patient_id is patient-level, not sample-level -- ambiguous for any cohort +# where a patient could have more than one sample (see hilton_study below +# for a case where that ambiguity is real). +thomas_bl_study <- data.frame( + sample_id = bl_data$meta_to_bundle$sample_id, + study = "Thomas", + study_id = as.character(bl_data$meta_to_bundle$sample_id), + reference_PMID = pmids$Thomas_BL +) + -# Importing FL data from Dreval et al +# --- Dreval FL (grch37, genome) --------------------------------------------- fl_data <- list() fl_data$meta <- read_xlsx("inst/extdata/studies/FL_Dreval.xlsx", sheet = 1) @@ -140,6 +358,7 @@ fl_data$meta_to_bundle <- fl_data$meta_to_bundle %>% filter(! sample_id == "NA") %>% arrange(sample_id) +diag_summary(fl_data$meta_to_bundle, "fl_data$meta_to_bundle", study_col = "cohort") fl_data$ssm_to_bundle <- read_xlsx( "inst/extdata/studies/FL_Dreval.xlsx", @@ -158,6 +377,8 @@ fl_data$ssm_to_bundle <- fl_data$ssm_to_bundle %>% select(names(GAMBLR.helpers:::maf_header[1:45])) +diag_summary_maf(fl_data$ssm_to_bundle, "fl_data$ssm_to_bundle (raw, as read from xlsx sheet 2)", study_lookup = data.frame(sample_id = fl_data$meta_to_bundle$sample_id, study = "FL_Dreval")) + fl_data$cnv_to_bundle <- read_xlsx( "inst/extdata/studies/FL_Dreval.xlsx", sheet = 3 @@ -167,9 +388,16 @@ fl_data$cnv_to_bundle <- read_xlsx( ) %>% mutate(CN = round(2 * 2^log.ratio)) +# study_id mirrors sample_id, not patient_id -- see thomas_bl_study above. +dreval_study <- data.frame( + sample_id = fl_data$meta_to_bundle$sample_id, + study = "Dreval", + study_id = as.character(fl_data$meta_to_bundle$sample_id), + reference_PMID = pmids$Dreval_FL +) -# Importing DLBCL data in hg38 from Thomas et al +# --- Thomas DLBCL (hg38, genome) -------------------------------------------- dlbcl_data <- list() dlbcl_data$meta <- read_xlsx("inst/extdata/studies/BL_Thomas.xlsx", sheet = 2) @@ -195,6 +423,8 @@ dlbcl_data$meta_to_bundle <- dlbcl_data$meta_to_bundle %>% filter(! sample_id == "NA") %>% arrange(sample_id) +diag_summary(dlbcl_data$meta_to_bundle, "dlbcl_data$meta_to_bundle", study_col = "cohort") + dlbcl_data$ssm_to_bundle <- read_xlsx( "inst/extdata/studies/BL_Thomas.xlsx", sheet = 6 @@ -219,38 +449,28 @@ dlbcl_data$cnv_to_bundle <- read_xlsx( filter(ID %in% dlbcl_data$meta_to_bundle$sample_id) %>% mutate(CN = round(2 * 2^log.ratio)) -##### Importing SLMS-3 variants -pull_data <- function( - pull_meta, - pull_projection = "grch37" - ){ - slms3 <- get_ssm_by_samples( - these_samples_metadata = pull_meta, - basic_columns = FALSE, - projection = pull_projection - ) %>% - select( - all_of(all_cols) - ) %>% - filter( - Hugo_Symbol %in% all_lymphoma_genes - ) - return(slms3) -} +# study_id mirrors sample_id, not patient_id -- see thomas_bl_study above. +thomas_dlbcl_study <- data.frame( + sample_id = dlbcl_data$meta_to_bundle$sample_id, + study = "Thomas", + study_id = as.character(dlbcl_data$meta_to_bundle$sample_id), + reference_PMID = pmids$Thomas_BL +) -# Importing DLBCL capture data +# --- DLBCL capture cohorts: Reddy, Schmitz, Chapuy, Golub ------------------- reddy_data <- list() schmitz_data <- list() chapuy_data <- list() golub_data <- list() # Importing metadata from Reddy et al and updating IDs to be consistent with GAMBL metadata -reddy_meta <- read_excel( +reddy_meta_full <- read_excel( "inst/extdata/studies/DLBCL_Reddy.xlsx", sheet = 1 ) %>% mutate( + study_id = `Sample ID`, patient_id = paste0( "Reddy_", `Sample ID` @@ -269,7 +489,18 @@ reddy_meta <- read_excel( "UNCLASS", COO_consensus ) - ) %>% + ) + +# Reddy's own paper-specific sample identifier ("Sample ID"), captured +# before it's dropped by the select() below -- see reddy_study. +reddy_study <- data.frame( + sample_id = reddy_meta_full$sample_id, + study = "Reddy", + study_id = as.character(reddy_meta_full$study_id), + reference_PMID = pmids$Reddy_DLBCL +) + +reddy_meta <- reddy_meta_full %>% dplyr::select( sample_id, patient_id, @@ -278,13 +509,17 @@ reddy_meta <- read_excel( COO_consensus ) -setwd("/projects/rmorin/projects/gambl-repos/gambl-kdreval/") +# (config is resolved via config.yml / the GAMBLR.results fallback; no chdir to +# a gambl repo is needed to locate data — repo_base is absolute) +# patient_id is deliberately retained here (unlike the original version of +# this select(), which dropped it) so all-capture samples get a real +# patient_id in sample_data$meta instead of NA. reddy_meta_gambl <- get_gambl_metadata() %>% dplyr::filter(cohort == "dlbcl_reddy") %>% dplyr::select( - sample_id, lymphgen, EBV_status_inf, cohort, pathology, seq_type, - unix_group, genome_build, pairing_status, normal_sample_id + sample_id, patient_id, lymphgen, EBV_status_inf, cohort, pathology, + seq_type, unix_group, genome_build, pairing_status, normal_sample_id ) %>% mutate(reference_PMID = pmids$Reddy_DLBCL) @@ -293,6 +528,9 @@ reddy_data$meta_to_bundle <- left_join( reddy_meta_gambl ) +diag_missing(reddy_meta$sample_id, reddy_data$meta_to_bundle$sample_id[!is.na(reddy_data$meta_to_bundle$cohort)], "Reddy meta -> reddy_meta_gambl join (missing cohort after join)") +diag_summary(reddy_data$meta_to_bundle, "reddy_data$meta_to_bundle", study_col = "cohort") + schmitz_data$meta <- get_gambl_metadata() %>% dplyr::filter(cohort == "dlbcl_schmitz") %>% mutate(reference_PMID = pmids$Schmitz_DLBCL) %>% @@ -327,11 +565,41 @@ all_capture_meta <- bind_rows( reddy_meta_gambl ) -# warning: this is very slow! -all_capture_grch37_ssm_to_bundle <- pull_data(all_capture_meta) -all_capture_hg38_ssm_to_bundle <- pull_data(all_capture_meta, "hg38") +diag_summary(all_capture_meta, "all_capture_meta", study_col = "cohort") + +# No paper-specific ID has been parsed anywhere in this script for these +# three cohorts (unlike Reddy/Thomas/Dreval/Arthur, they're pulled straight +# from GAMBL's own metadata, never cross-referenced against each paper's own +# supplementary sample table) -- study_id is NA until that sourcing is done. +schmitz_study <- data.frame( + sample_id = schmitz_data$meta$sample_id, + study = "Schmitz", + study_id = NA_character_, + reference_PMID = pmids$Schmitz_DLBCL +) +chapuy_study <- data.frame( + sample_id = chapuy_data$meta$sample_id, + study = "Chapuy", + study_id = NA_character_, + reference_PMID = pmids$Chapuy_DLBCL +) +golub_study <- data.frame( + sample_id = golub_data$meta$sample_id, + study = "NCI_Golub", + study_id = NA_character_, + reference_PMID = pmids$Chapuy_other +) -# Importing DLBCL cell lines +# Add data from Reddy paper's own original variant calls (distinct from the +# SLMS-3 recall above -- this is the as-published set, kept as its own +# Publication-pipeline pull, not deduplicated against SLMS-3). +reddy_original_maf <- read_tsv( + "inst/extdata/studies/reddy_original_variants_with_VAF.maf.gz" +) %>% + select(any_of(all_cols)) + + +# --- DLBCL cell lines -------------------------------------------------------- cell_lines_data <- list() cell_lines_data$meta <- get_gambl_metadata(seq_type_filter = "genome") %>% @@ -340,6 +608,9 @@ cell_lines_data$meta <- get_gambl_metadata(seq_type_filter = "genome") %>% )) %>% arrange(sample_id) +message(sprintf("[DIAG] cell_lines_data$meta: %d rows (expect 5): %s", + nrow(cell_lines_data$meta), paste(cell_lines_data$meta$sample_id, collapse = ", "))) + cell_lines_data$meta_to_bundle <- data.frame( cell_lines_data$meta$patient_id, cell_lines_data$meta$sample_id, @@ -357,42 +628,239 @@ cell_lines_data$meta_to_bundle <- data.frame( colnames(cell_lines_data$meta_to_bundle) <- colnames_for_bundled_meta -cell_lines_data$grch37$ssm_to_bundle <- get_ssm_by_samples( +cell_lines_data$grch37$cnv_to_bundle <- get_cn_segments( these_samples_metadata = cell_lines_data$meta, - basic_columns = FALSE -) %>% select(all_of(all_cols)) + projection="grch37" +) %>% + dplyr::select(all_of(c("ID","chrom","start","end","LOH_flag","log.ratio","CN","seg_seq_type"))) -cell_lines_data$hg38$ssm_to_bundle <- get_ssm_by_samples( +cell_lines_data$hg38$cnv_to_bundle <- get_cn_segments( these_samples_metadata = cell_lines_data$meta, - projection = "hg38", - basic_columns = FALSE -) %>% select(all_of(all_cols)) + projection="hg38" +) %>% + dplyr::select(all_of(c("ID","chrom","start","end","LOH_flag","log.ratio","CN","seg_seq_type"))) + +# No sample_study row for cell lines -- no published cohort applies to them +# (matches the previous Study=NA behaviour for this group). + + +# --- Arthur (grch37, genome) ------------------------------------------------- +# Previously bundled from a raw, completely unfiltered flat file +# (inst/extdata/studies/DLBCL_Arthur.maf.gz, tagged Pipeline="strelka") with +# no sample-list gate and no lymphoma-gene-panel restriction at all -- +# 2.83M rows / 65,121 distinct genes, ~18,480 "mutations" per sample vs +# ~137 for a properly gene-panel-restricted SLMS-3 sample. That block is +# removed entirely; Arthur's samples now flow through the same consolidated +# SLMS-3 pull as every other cohort (Phase 4). +# Case ID -> patient_id is done as an explicit join (not a +# `patient_id %in% arthur_case_ids$\`Case ID\`` filter) so a type mismatch +# between the xlsx's Case ID and GAMBL's patient_id (e.g. numeric vs +# character -- the same class of bug already fixed for study_id/Reddy's +# Sample ID elsewhere in this script; read_xlsx() can silently type a +# leading-zero ID like "08-15460" as numeric) can't silently drop matching +# patients from a %in% comparison. transmute() keeps only what's needed for +# the join, since nothing else from this sheet is used downstream. +arthur_case_ids <- read_xlsx( + "inst/extdata/studies/DLBCL_Arthur.xlsx", + sheet = 1 +) %>% filter(`WGS data` == 1) %>% + transmute(patient_id = as.character(`Case ID`)) + +# sample_id/Tumor_Sample_Barcode are GAMBL's own real values throughout -- +# never overwritten to the paper's own patient-level ID (as the removed +# code used to do). Arthur's own "Case ID" is captured separately as +# study_id in arthur_study below, joined to metadata by sample_id like any +# other study-specific identifier. +# +# No `!grepl("tumor", sample_id)` filter here (the previous version of this +# block had one): it excluded every sample for any patient with more than +# one tumor biopsy (e.g. ..._tumorA/..._tumorB -- the same multi-sample- +# per-patient pattern Hilton has), dropping those patients' SLMS-3 coverage +# entirely. No other cohort in this script excludes samples this way. +arthur_meta <- get_gambl_metadata() %>% + mutate(patient_id = as.character(patient_id)) %>% + inner_join(arthur_case_ids, by = "patient_id") %>% + filter(seq_type == "genome") %>% + mutate( + cohort = "DLBCL_Arthur", + reference_PMID = pmids$Arthur_DLBCL + ) + +diag_missing(arthur_case_ids$patient_id, arthur_meta$patient_id, "Arthur Case ID -> patient_id join (WGS-flagged patients with no matching genome sample in get_gambl_metadata())") +diag_summary(arthur_meta, "arthur_meta", study_col = "cohort") + +arthur_study <- data.frame( + sample_id = arthur_meta$sample_id, + study = "Arthur", + study_id = as.character(arthur_meta$patient_id), + reference_PMID = pmids$Arthur_DLBCL +) + + +# --- Hilton trios (genome + capture, both projections) ---------------------- +trios_samples <- read_xlsx( + "inst/extdata/studies/DLBCL_Hilton.xlsx" +) %>% +drop_na(DNAseq_sample_id) + +# (config is resolved via config.yml / the GAMBLR.results fallback; no chdir to +# a gambl repo is needed to locate data — repo_base is absolute) + +trios_meta <- get_gambl_metadata() %>% + filter( + seq_type %in% c("genome", "capture"), + sample_id %in% trios_samples$DNAseq_sample_id + ) %>% + mutate( + cohort = "DLBCL_Hilton", + reference_PMID = pmids$Hilton_DLBCL + ) -cell_lines_data$grch37$cnv_to_bundle <- get_sample_cn_segments( - these_sample_ids = cell_lines_data$meta$sample_id +diag_missing(trios_samples$DNAseq_sample_id, trios_meta$sample_id, "Hilton trios DNAseq_sample_id -> get_gambl_metadata() (samples in the trios sheet not found live)") +diag_summary(trios_meta, "trios_meta", study_col = "cohort") + +# study_id uses sample_id, not patient_id: Hilton is a trios study, so a +# single patient can have multiple samples (e.g. LY_RELY_116_tumorA and +# LY_RELY_116_tumorB) -- patient_id would collapse them to the same +# study_id, making the two rows ambiguous. sample_id already carries the +# distinguishing suffix and matches what the paper itself would call each +# sample. +hilton_study <- data.frame( + sample_id = trios_meta$sample_id, + study = "Hilton", + study_id = as.character(trios_meta$sample_id), + reference_PMID = pmids$Hilton_DLBCL ) -cell_lines_data$hg38$cnv_to_bundle <- get_sample_cn_segments( - these_sample_ids = cell_lines_data$meta$sample_id, - projection = "hg38", - with_chr_prefix = TRUE + +# ============================================================================ +# Phase 2: assemble sample_data$meta and sample_data$sample_study once, now +# that every cohort's metadata (including Arthur and Hilton, previously +# added much later -- after the main aSHM pull had already run without them) +# is available from the start. +# ============================================================================ +sample_data <- list() + +sample_data$meta <- bind_rows( + bl_data$meta_to_bundle, + fl_data$meta_to_bundle, + dlbcl_data$meta_to_bundle, + cell_lines_data$meta_to_bundle ) -cell_lines_data$grch37$sv_to_bundle <- get_manta_sv( - these_samples_metadata = cell_lines_data$meta, +sample_data$meta <- sample_data$meta %>% + select(-COO_consensus, -lymphgen, -EBV_status_inf) %>% + left_join( + ., + get_gambl_metadata() %>% + select( + sample_id, + COO_consensus, + lymphgen, + EBV_status_inf + ) + ) %>% + select(all_of(colnames_for_bundled_meta)) + +sample_data$meta <- bind_rows( + sample_data$meta, + all_capture_meta, + arthur_meta, + trios_meta ) -cell_lines_data$hg38$sv_to_bundle <- get_manta_sv( - these_samples_metadata = cell_lines_data$meta, - projection = "hg38" +diag_summary(sample_data$meta, "sample_data$meta (all cohorts bound, pre metadata-fixing)", study_col = "cohort") + +### begin metadata fixing +# This preserves the original cohort column and ensures there are no duplicates +# in the metadata +fix <- sample_data$meta +fix <- fix %>% rename(study = cohort) + +diag_summary(fix, "fix (post rename cohort->study)", study_col = "study") + +# Keep only sample_id/seq_type/study/reference_PMID/genetic_subgroup from +# the per-block frames -- everything else comes fresh from the join below +# instead of whatever narrow/partial copy a given cohort block happened to +# carry. genetic_subgroup is kept here specifically because it's paper- +# supplement-specific (from the Thomas/Dreval xlsx sheets, not central +# tracking) and would NOT be reintroduced by that join -- live +# get_gambl_metadata() has no such column. +fix <- fix %>% + select(sample_id, seq_type, study, reference_PMID, genetic_subgroup) + +pre_join_ids <- fix$sample_id +fix <- left_join( + fix, + get_gambl_metadata(), + by = c("sample_id", "seq_type") ) -# Adding the manta SVs for published studies +message(sprintf("[DIAG] fix (post live get_gambl_metadata() join): %d rows (was %d rows pre-join -- any increase means the join fanned out)", nrow(fix), length(pre_join_ids))) +diag_summary(fix, "fix (post live get_gambl_metadata() join)", study_col = "study") + +# sample_meta ends up in the PUBLICLY-DISTRIBUTED gambl_mutations.db (a +# GitHub Release asset anyone can download and query directly with SQL) -- +# restrict it to an explicit whitelist rather than everything live +# get_gambl_metadata() happens to return, which includes internal-only +# columns (biopsy_id, data_path, fastq_data_path/fastq_link_name, +# link_name, library_id, res_id, unix_group, ...) that must never end up +# in a public artifact. The whitelist lives in a plain text file, not +# here, so a developer can add a column to future builds by editing that +# file -- no code change needed. all_of() (not any_of()) deliberately +# errors on a typo'd/renamed column name rather than silently dropping it. +public_metadata_cols <- readLines("data-raw/public_sample_meta_columns.txt") +public_metadata_cols <- trimws(public_metadata_cols) +public_metadata_cols <- public_metadata_cols[ + nzchar(public_metadata_cols) & !startsWith(public_metadata_cols, "#") +] +fix <- fix %>% + select(sample_id, seq_type, study, reference_PMID, all_of(public_metadata_cols)) + +diag_summary(fix, "fix (post public-column whitelist)", study_col = "study") + +fix <- fix %>% filter(!is.na(study)) + +diag_missing(pre_join_ids, fix$sample_id, "fix: post-join !is.na(study) filter (sample_ids dropped here)") +diag_summary(fix, "fix (post filter !is.na(study))", study_col = "study") + +fix <- distinct(fix) + +diag_summary(fix, "fix (post distinct -- this becomes sample_data$meta)", study_col = "study") + +sample_data$meta <- fix +### end metadata fixing + +# sample_study: a many-to-many bridge table (one row per (sample_id, study) +# pair, not one row per sample), replacing the old inline `Study =` tagging +# on mutation rows. A sample belonging to N studies is simply N rows here -- +# no schema change needed as multi-study overlap becomes more common. +sample_data$sample_study <- bind_rows( + thomas_bl_study, + thomas_dlbcl_study, + dreval_study, + reddy_study, + schmitz_study, + chapuy_study, + golub_study, + arthur_study, + hilton_study +) %>% distinct() + +diag_summary(sample_data$sample_study, "sample_data$sample_study", study_col = "study") + + +#add SVs +# Adding the manta SVs for published studies. Scoped to sample_data$meta (the +# metadata just finalized above, for this run) rather than +# GAMBLR.data::sample_data$meta (whatever happens to be installed) so the SV +# sample scope is self-consistent with the rest of this bundle instead of +# drifting with installed-package state across runs. full_genome_meta <- get_gambl_metadata(seq_type_filter = "genome") bundled_meta <- full_genome_meta %>% filter( - sample_id %in% GAMBLR.data::sample_data$meta$sample_id + sample_id %in% sample_data$meta$sample_id ) full_sv_to_bundle <- get_manta_sv( @@ -449,92 +917,81 @@ annotated_sv_keep_grch37 <- left_join( filter(!is.na(partner)) %>% select(c(1:16)) - -# Combine everything together -sample_data <- list() - -sample_data$meta <- bind_rows( - bl_data$meta_to_bundle, - fl_data$meta_to_bundle, - dlbcl_data$meta_to_bundle, - cell_lines_data$meta_to_bundle -) - -sample_data$meta <- sample_data$meta %>% - select(-COO_consensus, -lymphgen, -EBV_status_inf) %>% - left_join( - ., - get_gambl_metadata() %>% - select( - sample_id, - COO_consensus, - lymphgen, - EBV_status_inf - ) - ) %>% - select(all_of(colnames_for_bundled_meta)) - -sample_data$meta <- bind_rows( - sample_data$meta, - all_capture_meta -) - -sample_data$hg38$maf <- bind_rows( - bl_data$ssm_to_bundle %>% mutate( - Pipeline = "Publication", - Study = "Thomas" - ), - dlbcl_data$ssm_to_bundle %>% mutate( - Pipeline = "Publication", - Study = "Thomas" - ), - cell_lines_data$hg38$ssm %>% mutate( - Pipeline = "SLMS-3", - Study = NA - ), - all_capture_hg38_ssm_to_bundle %>% mutate( - Pipeline = "SLMS-3", - Study = case_when( - Tumor_Sample_Barcode %in% reddy_data$meta$sample_id ~ "Reddy", - Tumor_Sample_Barcode %in% schmitz_data$meta$sample_id ~ "Schmitz", - Tumor_Sample_Barcode %in% chapuy_data$meta$sample_id ~ "Chapuy", - Tumor_Sample_Barcode %in% golub_data$meta$sample_id ~ "NCI_Golub" - ) - ) -) +sample_data$grch37$bedpe <- annotated_sv_keep_grch37 +sample_data$hg38$bedpe <- annotated_sv_keep sample_data$grch37$seg <- bind_rows( fl_data$cnv_to_bundle, - cell_lines_data$grch37$cnv + cell_lines_data$grch37$cnv_to_bundle ) sample_data$hg38$seg <- bind_rows( bl_data$cnv_to_bundle, dlbcl_data$cnv_to_bundle, - cell_lines_data$hg38$cnv + cell_lines_data$hg38$cnv_to_bundle ) -#add SVs -sample_data$grch37$bedpe <- annotated_sv_keep_grch37 -sample_data$hg38$bedpe <- annotated_sv_keep -# This is needed for the proteinpainter compatibility -GAMBLR.data::sample_data$meta$cohort %>% table +# ============================================================================ +# Phase 3: Publication-pipeline SSM data. Genuinely distinct per-paper data +# (not a redundant SLMS-3 recall), kept as separate, un-deduplicated pulls -- +# tagged Pipeline="Publication" only; cohort membership lives in +# sample_data$sample_study, not an inline Study column. +# ============================================================================ + +# Publication-pipeline data is read directly from each paper's own +# supplementary file, which isn't guaranteed to use GAMBL's own +# Tumor_Sample_Barcode convention -- Arthur's now-removed raw dump was the +# clearest example, using bare patient-style IDs instead of GAMBL's real +# per-sample naming. Defensive relabel: for any row whose +# Tumor_Sample_Barcode matches a study's own identifier +# (sample_study$study_id) rather than GAMBL's real sample_id, replace it +# with the real sample_id. A no-op wherever Tumor_Sample_Barcode already IS +# the real sample_id -- true for Thomas/Dreval/Hilton, where study_id +# mirrors sample_id (see thomas_bl_study above). +relabel_to_sample_id <- function(df, study_name) { + lookup <- sample_data$sample_study %>% + filter(study == study_name, !is.na(study_id)) %>% + select(study_id, .gambl_sample_id = sample_id) + df %>% + left_join(lookup, by = c("Tumor_Sample_Barcode" = "study_id")) %>% + mutate(Tumor_Sample_Barcode = coalesce(.gambl_sample_id, Tumor_Sample_Barcode)) %>% + select(-.gambl_sample_id) +} + +# This is needed for the proteinpainter compatibility. Previously read +# from the *installed* GAMBLR.data::sample_data rather than the +# locally-built sample_data above -- a stale-data risk (the installed copy +# could be out of sync with this run, and is the reason sample_data.rda +# was still shipped in the package at all) fixed by switching to the local +# object, which by this point (Phase 2 already ran) has every cohort this +# block needs. +sample_data$meta$study %>% table selected_columns <- c( "Tumor_Sample_Barcode", "Hugo_Symbol", "NCBI_Build", "Chromosome", "Start_Position", "End_Position", - "Tumor_Seq_Allele2", maf_columns_to_keep + "Tumor_Seq_Allele2", "RefSeq", "Protein_position" + # Deliberately NOT maf_columns_to_keep (which now also carries + # maf_seq_type) -- these are join-source columns for the + # proteinpainter-compatibility enrichment below, not the final maf + # table. Two of the three sources selected against this (the raw + # BL_Thomas manuscript flat file, and the get_ssm_by_samples() FL + # pull below) never carry maf_seq_type at all, so all_of() on the + # full list would hard-error here every time. ) -these_samples <- GAMBLR.data::sample_data$meta %>% - filter(cohort %in% c("BL_Thomas")) %>% +these_samples <- sample_data$meta %>% + filter(study %in% c("BL_Thomas")) %>% pull(sample_id) -these_samples_dlbcl <- GAMBLR.data::sample_data$meta %>% - filter(cohort %in% c("DLBCL_Thomas", "DLBCL_cell_lines")) %>% +these_samples_dlbcl <- sample_data$meta %>% + filter(study %in% c("DLBCL_Thomas", "DLBCL_cell_lines")) %>% pull(sample_id) +message(sprintf("[DIAG] these_samples (BL_Thomas): %d ids", length(these_samples))) +message(sprintf("[DIAG] these_samples_dlbcl (DLBCL_Thomas + DLBCL_cell_lines): %d ids", length(these_samples_dlbcl))) + coding_maf <- read_tsv("/projects/adult_blgsp/results_manuscript/BL.hg38.CDS.maf") %>% # get from flat maf file to show SSM in hg38 coordinates similar to the original manuscript filter(Tumor_Sample_Barcode %in% these_samples & # drop BL58 cell line ! str_detect(Tumor_Sample_Barcode, "^SP|^06")) %>% # drop ICGC and 1 LLMPP case @@ -558,409 +1015,296 @@ coding_maf <- bind_rows( coding_maf_dlbcl ) -dim(GAMBLR.data::sample_data$hg38$maf) +message(sprintf("[DIAG] coding_maf (BL+DLBCL Thomas, hg38 proteinpainter enrichment source): %d rows, %d distinct Tumor_Sample_Barcode", nrow(coding_maf), n_distinct(coding_maf$Tumor_Sample_Barcode))) -sample_data$hg38$maf <- sample_data$hg38$maf %>% +# hg38 Publication rows (Thomas BL + Thomas DLBCL), enriched via the +# proteinpainter-compatibility join above. This used to be applied directly +# to sample_data$hg38$maf once it existed early in the script; now it's +# applied to the Publication-only rows before Phase 6 combines them with the +# consolidated SLMS-3 pull. +hg38_publication_rows <- bind_rows( + bl_data$ssm_to_bundle %>% mutate(Pipeline = "Publication"), + dlbcl_data$ssm_to_bundle %>% mutate(Pipeline = "Publication") +) %>% + relabel_to_sample_id("Thomas") %>% left_join(coding_maf) -this_study_samples <- GAMBLR.data::sample_data$meta %>% - filter(cohort %in% c("FL_Dreval", "DLBCL_cell_lines")) %>% +diag_summary_maf(hg38_publication_rows, "hg38_publication_rows (Thomas BL+DLBCL, final)") + +this_study_samples <- sample_data$meta %>% + filter(study %in% c("FL_Dreval", "DLBCL_cell_lines")) %>% pull(sample_id) +message(sprintf("[DIAG] this_study_samples (FL_Dreval + DLBCL_cell_lines): %d ids", length(this_study_samples))) + # FLs in grch37 -coding_maf <- get_ssm_by_samples( - these_samples_metadata = get_gambl_metadata() %>% - filter(sample_id %in% this_study_samples), - basic_columns = FALSE) %>% - select( - all_of(selected_columns) - ) +coding_maf <- time_it("coding_maf get_ssm_by_samples", { + get_ssm_by_samples( + these_samples_metadata = get_gambl_metadata() %>% + filter(sample_id %in% this_study_samples), + basic_columns = FALSE) %>% + select( + all_of(selected_columns) + ) +}) + +message(sprintf("[DIAG] coding_maf (FL_Dreval + DLBCL_cell_lines, grch37 proteinpainter enrichment source): %d rows, %d distinct Tumor_Sample_Barcode", nrow(coding_maf), n_distinct(coding_maf$Tumor_Sample_Barcode))) + +diag_summary_maf(fl_data$ssm_to_bundle, "fl_data$ssm_to_bundle (pre relabel/enrichment)", study_lookup = data.frame(sample_id = fl_data$meta_to_bundle$sample_id, study = "FL_Dreval")) fl_data$ssm_to_bundle <- fl_data$ssm_to_bundle %>% + relabel_to_sample_id("Dreval") %>% dplyr::left_join( coding_maf ) %>% distinct() - -sample_data$grch37$maf <- bind_rows( - fl_data$ssm_to_bundle %>% mutate( - Pipeline = "Publication", - Study = "Dreval" - ), - cell_lines_data$grch37$ssm %>% mutate( - Pipeline = "SLMS-3", - Study = NA - ), - all_capture_grch37_ssm_to_bundle %>% mutate( - Pipeline = "SLMS-3", - Study = case_when( - Tumor_Sample_Barcode %in% reddy_data$meta$sample_id ~ "Reddy", - Tumor_Sample_Barcode %in% schmitz_data$meta$sample_id ~ "Schmitz", - Tumor_Sample_Barcode %in% chapuy_data$meta$sample_id ~ "Chapuy", - Tumor_Sample_Barcode %in% golub_data$meta$sample_id ~ "NCI_Golub" - ) - ) -) - - -# Add aSHM mutations for the already released samples -grch37_ashm <- get_ssm_by_regions( - these_samples_metadata = sample_data$meta, - regions_bed = GAMBLR.utils::create_bed_data( - GAMBLR.data::grch37_ashm_regions, - fix_names = "concat", - concat_cols = c("gene","region"),sep="-" - ), - streamlined = FALSE, - basic_columns = FALSE -) %>% - select( - any_of(c(colnames(sample_data$grch37$maf), maf_columns_to_keep)) - ) - -grch37_ashm <- grch37_ashm %>% - filter(Tumor_Sample_Barcode %in% sample_data$meta$Tumor_Sample_Barcode) - -grch37_ashm <- grch37_ashm %>% mutate(Pipeline = "SLMS-3") - -studies <- bind_rows( - sample_data$grch37$maf %>% - distinct(Tumor_Sample_Barcode, Study), - sample_data$hg38$maf %>% - distinct(Tumor_Sample_Barcode, Study) -) %>% -distinct() - -grch37_ashm <- left_join( - grch37_ashm, - studies -) - -hg38_ashm <- get_ssm_by_regions( - these_samples_metadata = sample_data$meta, - regions_bed = GAMBLR.utils::create_bed_data( - GAMBLR.data::hg38_ashm_regions, - fix_names = "concat", - concat_cols = c("gene","region"),sep="-" - ), - projection = "hg38", - streamlined = FALSE, - basic_columns = FALSE -) %>% - select( - any_of(c(colnames(sample_data$hg38$maf), maf_columns_to_keep)) - ) -hg38_ashm <- hg38_ashm %>% - filter(Tumor_Sample_Barcode %in% sample_data$meta$Tumor_Sample_Barcode) - -hg38_ashm <- hg38_ashm %>% mutate(Pipeline = "SLMS-3") - -hg38_ashm <- left_join( - hg38_ashm, - studies -) - -sample_data$grch37$ashm <- grch37_ashm -sample_data$hg38$ashm <- hg38_ashm - - -# Now add the SLMS-3 calls in both projections for those samples that -# are bundled as publication data -publication_samples_grch37 <- sample_data$grch37$maf %>% - filter(Pipeline == "Publication") %>% - pull(Tumor_Sample_Barcode) %>% - unique %>% sort - -publication_samples_hg38 <- sample_data$hg38$maf %>% - filter(Pipeline == "Publication") %>% - pull(Tumor_Sample_Barcode) %>% - unique %>% sort - -publication_samples <- c( - publication_samples_grch37, - publication_samples_hg38 -) - -sample_data$grch37$maf <- get_ssm_by_samples( - these_samples_metadata = get_gambl_metadata() %>% - filter(sample_id %in% publication_samples), - basic_columns = FALSE) %>% - filter( - Hugo_Symbol %in% all_lymphoma_genes - ) %>% - mutate(Pipeline = "SLMS-3") %>% - left_join( - ., - studies - ) %>% - select(colnames(sample_data$grch37$maf)) %>% - bind_rows( - ., - sample_data$grch37$maf - ) - -sample_data$hg38$maf <- get_ssm_by_samples( - these_samples_metadata = get_gambl_metadata() %>% - filter(sample_id %in% publication_samples), - projection = "hg38", - basic_columns = FALSE) %>% - filter( - Hugo_Symbol %in% all_lymphoma_genes - ) %>% - mutate(Pipeline = "SLMS-3") %>% - left_join( - ., - studies - ) %>% - select(colnames(sample_data$hg38$maf)) %>% - bind_rows( - ., - sample_data$hg38$maf - ) - - -setwd("~/my_dir/repos/GAMBLR.data/") - -# Add data from Arthur paper -arthur_maf <- read_tsv( - "inst/extdata/studies/DLBCL_Arthur.maf.gz" -) -sample_data$grch37$maf <- bind_rows( - sample_data$grch37$maf, - arthur_maf %>% - mutate( - Pipeline = "strelka", - Study = "Arthur" - ) -) - -arthur_meta <- read_xlsx( - "inst/extdata/studies/DLBCL_Arthur.xlsx", - sheet = 1 -) %>% filter(`WGS data` == 1) - - -arthur_meta <- gambl_metadata %>% - filter( - patient_id %in% arthur_meta$`Case ID`, - seq_type == "genome", - ! grepl("tumor", sample_id) - ) %>% - mutate( - sample_id = patient_id, - Tumor_Sample_Barcode = patient_id, - cohort = "DLBCL_Arthur", - reference_PMID = pmids$Arthur_DLBCL - ) - -sample_data$meta <- bind_rows( - sample_data$meta, - arthur_meta -) -# Add data from Hilton trios paper -trios_samples <- read_xlsx( - "inst/extdata/studies/DLBCL_Hilton.xlsx" -) %>% -drop_na(DNAseq_sample_id) - -setwd("/projects/rmorin/projects/gambl-repos/gambl-kdreval/") - -trios_meta <- get_gambl_metadata() %>% - filter( - seq_type %in% c("genome", "capture"), - sample_id %in% trios_samples$DNAseq_sample_id - ) %>% - select(any_of(colnames(sample_data$meta))) %>% - mutate( - cohort = "DLBCL_Hilton", - reference_PMID = pmids$Hilton_DLBCL - ) - -sample_data$meta <- bind_rows( - sample_data$meta, - trios_meta -) - -### begin metadata fixing -# This preserves the original cohort column and ensures there are no duplicates -# in the metadata -fix <- sample_data$meta -fix <- fix %>% rename(study = cohort) - -fix <- left_join( - fix, - get_gambl_metadata() %>% - select(sample_id, seq_type, cohort) -) - -fix <- fix %>% filter(!is.na(study)) - -fix <- distinct(fix) - -sample_data$meta <- fix -### end metadata fixing - -# trios grch37 ssm -genome_trios_ssm_grch37 <- get_ssm_by_samples( - these_samples_metadata = trios_meta %>% - filter(seq_type == "genome"), - basic_columns = FALSE, - subset_from_merge = TRUE -) %>% - filter(Hugo_Symbol %in% all_lymphoma_genes) %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" - ) %>% - select(all_of(colnames(sample_data$grch37$maf))) - -capture_trios_ssm_grch37 <- get_ssm_by_samples( - these_samples_metadata = trios_meta %>% - filter(seq_type == "capture"), - basic_columns = FALSE, - subset_from_merge = TRUE -) %>% - filter(Hugo_Symbol %in% all_lymphoma_genes) %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" - ) %>% - select(all_of(colnames(sample_data$grch37$maf))) +diag_summary_maf(fl_data$ssm_to_bundle, "fl_data$ssm_to_bundle (post relabel/enrichment/distinct -- final)", study_lookup = data.frame(sample_id = fl_data$meta_to_bundle$sample_id, study = "FL_Dreval")) -trios_ssm_grch37 <- bind_rows( - genome_trios_ssm_grch37, - capture_trios_ssm_grch37 +# grch37 Publication rows (Dreval FL + Reddy's original variants). +grch37_publication_rows <- bind_rows( + fl_data$ssm_to_bundle %>% mutate(Pipeline = "Publication"), + reddy_original_maf %>% relabel_to_sample_id("Reddy") %>% mutate(Pipeline = "Publication") ) -# trios hg38 ssm -genome_trios_ssm_hg38 <- get_ssm_by_samples( - these_samples_metadata = trios_meta %>% - filter(seq_type == "genome"), - basic_columns = FALSE, - subset_from_merge = TRUE, - projection = "hg38" -) %>% - filter(Hugo_Symbol %in% all_lymphoma_genes) %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" +diag_summary_maf(grch37_publication_rows, "grch37_publication_rows (Dreval FL + Reddy original, final)") + +setwd(PKG_ROOT) + + +# ============================================================================ +# Phase 4: one consolidated SLMS-3 pull for every sample in the bundle, +# instead of the previous per-cohort blocks (all-capture, cell lines, Hilton +# trios, Arthur) each independently deciding which samples to pull for -- +# that pattern silently duplicated mutation rows whenever a sample belonged +# to more than one cohort block. Pulling once per (sample, genome_build), +# keyed off the now-complete deduplicated sample_data$meta, makes +# duplication structurally impossible regardless of how many cohorts a +# sample belongs to. +# ============================================================================ +# Cell lines are excluded here -- they get their own separate, genome-wide +# pull below (not restricted to the lymphoma gene panel), so they must not +# also go through this panel-restricted pull or they'd be pulled twice. +# +# sample_data$meta is used ONLY to decide which sample_ids belong in this +# pull -- never passed directly to a GAMBLR.results call. Cohorts built by +# hand (BL_Thomas, DLBCL_Thomas, FL_Dreval, DLBCL_cell_lines) don't have +# every column a fresh get_gambl_metadata() pull would (e.g. unix_group is +# NA for all of them), which silently caused get_ssm_by_regions() to miss +# real coding-classified calls for those samples (confirmed: a raw pull for +# one such sample had real Missense_Mutation/Nonsense_Mutation/Silent rows +# that never made it into the assembled maf table). Re-fetching complete, +# live metadata for exactly this sample_id set avoids that entirely. +all_slms3_meta <- sample_data$meta %>% + filter(seq_type %in% c("genome", "capture"), + ! sample_id %in% cell_lines_data$meta$sample_id) + +diag_summary(all_slms3_meta, "all_slms3_meta (pre live-metadata refetch)", study_col = "study") +pre_refetch_ids <- all_slms3_meta$sample_id + +all_slms3_meta <- get_gambl_metadata() %>% + filter(sample_id %in% all_slms3_meta$sample_id) + +# The single most important check in this script: if get_gambl_metadata()'s +# live source is missing sample_ids that sample_data$meta has (e.g. an +# externally-integrated cohort not fully registered there), this filter +# silently drops them with no error -- every study's SLMS-3 coverage for +# those samples goes to zero downstream. diag_missing() catches this by name +# instead of leaving it to be inferred later from a low mutation count. +diag_missing(pre_refetch_ids, all_slms3_meta$sample_id, "all_slms3_meta live-metadata refetch (sample_ids dropped because live get_gambl_metadata() doesn't have them)") +# study isn't in all_slms3_meta itself -- join sample_study just for this +# diagnostic so the breakdown is by study, not left un-grouped. Drop the +# existing (unrelated) Tumor_Sample_Barcode column first: get_gambl_metadata() +# already returns one alongside sample_id, so renaming sample_id would +# otherwise collide with it. +diag_summary_maf(all_slms3_meta %>% select(-any_of("Tumor_Sample_Barcode")) %>% rename(Tumor_Sample_Barcode = sample_id), "all_slms3_meta (post live-metadata refetch -- what Phase 4 actually pulls for)") + +slms3_grch37 <- bind_rows( + time_it("SLMS-3 genome grch37", pull_data(all_slms3_meta %>% filter(seq_type == "genome"))), + time_it("SLMS-3 capture grch37", pull_data(all_slms3_meta %>% filter(seq_type == "capture"))) +) %>% mutate(Pipeline = "SLMS-3") +print("Done collecting grch37 SLMS-3") +diag_summary_maf(slms3_grch37, "slms3_grch37 (panel-restricted pull, pre cell-line union)") + +slms3_hg38 <- bind_rows( + time_it("SLMS-3 genome hg38", pull_data(all_slms3_meta %>% filter(seq_type == "genome"), "hg38")), + time_it("SLMS-3 capture hg38", pull_data(all_slms3_meta %>% filter(seq_type == "capture"), "hg38")) +) %>% mutate(Pipeline = "SLMS-3") +print("Done collecting hg38 SLMS-3") +diag_summary_maf(slms3_hg38, "slms3_hg38 (panel-restricted pull, pre cell-line union)") + +# Cell lines get their own, separate, genome-WIDE SNV pull -- not restricted +# to the lymphoma gene panel like every other cohort above. This matches +# their original pre-refactor behaviour, which was lost when they were first +# folded into the panel-restricted consolidated pull above (confirmed via +# compare_bundle_changes.R: cell lines showed ~50-65k "lost" rows per sample +# against the old bundle, correctly diagnosed as a real scope reduction, not +# a bug worth working around). +cell_lines_ssm_grch37 <- time_it("cell lines get_ssm_by_samples grch37", { + get_ssm_by_samples( + these_samples_metadata = cell_lines_data$meta, + basic_columns = FALSE ) %>% - select(all_of(colnames(sample_data$hg38$maf))) - -capture_trios_ssm_hg38 <- get_ssm_by_samples( - these_samples_metadata = trios_meta %>% - filter(seq_type == "capture"), - basic_columns = FALSE, - subset_from_merge = TRUE, - projection = "hg38" -) %>% - filter(Hugo_Symbol %in% all_lymphoma_genes) %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" + # get_ssm_by_samples() doesn't stamp maf_seq_type the way + # get_ssm_by_regions() does (see all_slms3_meta pulls above) -- + # cell_lines_data$meta is unconditionally get_gambl_metadata(seq_type_filter + # = "genome") (see its construction above), so "genome" is correct here, + # not a guess. + mutate(maf_seq_type = "genome") %>% + select(all_of(all_cols)) %>% mutate(Pipeline = "SLMS-3") +}) + +cell_lines_ssm_hg38 <- time_it("cell lines get_ssm_by_samples hg38", { + get_ssm_by_samples( + these_samples_metadata = cell_lines_data$meta, + projection = "hg38", + basic_columns = FALSE ) %>% - select(all_of(colnames(sample_data$hg38$maf))) - -trios_ssm_hg38 <- bind_rows( - genome_trios_ssm_hg38, - capture_trios_ssm_hg38 -) - -sample_data$grch37$maf <- bind_rows( - sample_data$grch37$maf, - trios_ssm_grch37 -) - -sample_data$hg38$maf <- bind_rows( - sample_data$hg38$maf, - trios_ssm_hg38 -) - -regions_bed <- create_bed_data( + mutate(maf_seq_type = "genome") %>% + select(all_of(all_cols)) %>% mutate(Pipeline = "SLMS-3") +}) + +message(sprintf("[DIAG] cell_lines_ssm_grch37: %d rows, %d distinct Tumor_Sample_Barcode (expect 5)", nrow(cell_lines_ssm_grch37), n_distinct(cell_lines_ssm_grch37$Tumor_Sample_Barcode))) +message(sprintf("[DIAG] cell_lines_ssm_hg38: %d rows, %d distinct Tumor_Sample_Barcode (expect 5)", nrow(cell_lines_ssm_hg38), n_distinct(cell_lines_ssm_hg38$Tumor_Sample_Barcode))) + +slms3_grch37 <- bind_rows(slms3_grch37, cell_lines_ssm_grch37) +slms3_hg38 <- bind_rows(slms3_hg38, cell_lines_ssm_hg38) +print("Done collecting cell line SLMS-3 (genome-wide)") +diag_summary_maf(slms3_grch37, "slms3_grch37 (final, incl. cell lines)") +diag_summary_maf(slms3_hg38, "slms3_hg38 (final, incl. cell lines)") + + +# ============================================================================ +# Phase 5: one consolidated aSHM pull, using the now-complete sample_data$meta +# (includes Arthur + Hilton from the start). This replaces both the old main +# aSHM block (which ran before Arthur/Hilton were added to sample_data$meta, +# so neither got any aSHM coverage from it) and Hilton's separate dedicated +# aSHM block (which existed only to compensate for that gap). Arthur gets +# aSHM coverage for the first time as a result. +# ============================================================================ +# sample_data$meta is used only to fix the sample_id set -- see the same +# rationale next to all_slms3_meta above for why a fresh get_gambl_metadata() +# pull is used for the actual GAMBLR.results call instead of sample_data$meta +# directly. +ashm_pull_meta <- get_gambl_metadata() %>% + filter(sample_id %in% sample_data$meta$sample_id) + +diag_missing(sample_data$meta$sample_id, ashm_pull_meta$sample_id, "ashm_pull_meta live-metadata refetch (sample_ids dropped because live get_gambl_metadata() doesn't have them)") +diag_summary_maf(ashm_pull_meta %>% select(-any_of("Tumor_Sample_Barcode")) %>% rename(Tumor_Sample_Barcode = sample_id), "ashm_pull_meta (what Phase 5 actually pulls for)") + +regions_bed_grch37 <- create_bed_data( grch37_ashm_regions, fix_names = "concat", - concat_cols = c("gene", "region"), - sep = "-" + concat_cols = c("gene", "region"), sep = "-" ) -trios_ashm_grch37 <- get_ssm_by_regions( - these_samples_metadata = trios_meta, - regions_bed = regions_bed, - streamlined = FALSE, - basic_columns = FALSE -) - -trios_ashm_grch37 <- trios_ashm_grch37 %>% - filter(Tumor_Sample_Barcode %in% trios_meta$Tumor_Sample_Barcode) - -trios_ashm_grch37 <- trios_ashm_grch37 %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" +grch37_ashm <- time_it("grch37_ashm get_ssm_by_regions", { + get_ssm_by_regions( + these_samples_metadata = ashm_pull_meta, + regions_bed = regions_bed_grch37, + streamlined = FALSE, + basic_columns = FALSE ) %>% - select(all_of(colnames(sample_data$grch37$maf))) - - -sample_data$grch37$ashm <- bind_rows( - sample_data$grch37$ashm, - trios_ashm_grch37 -) %>% distinct - + select(any_of(all_cols)) +}) +grch37_ashm <- grch37_ashm %>% + filter(Tumor_Sample_Barcode %in% sample_data$meta$Tumor_Sample_Barcode) %>% + mutate(Pipeline = "SLMS-3") -regions_bed <- create_bed_data( +regions_bed_hg38 <- create_bed_data( hg38_ashm_regions, fix_names = "concat", - concat_cols = c("gene", "region"), - sep = "-" + concat_cols = c("gene", "region"), sep = "-" ) -trios_ashm_hg38 <- get_ssm_by_regions( - these_samples_metadata = trios_meta, - regions_bed = regions_bed, - projection = "hg38", - streamlined = FALSE, - basic_columns = FALSE -) -trios_ashm_hg38 <- trios_ashm_hg38 %>% - filter(Tumor_Sample_Barcode %in% trios_meta$Tumor_Sample_Barcode) - -trios_ashm_hg38 <- trios_ashm_hg38 %>% - mutate( - Pipeline = "SLMS-3", - Study = "Hilton" +hg38_ashm <- time_it("hg38_ashm get_ssm_by_regions", { + get_ssm_by_regions( + these_samples_metadata = ashm_pull_meta, + regions_bed = regions_bed_hg38, + projection = "hg38", + streamlined = FALSE, + basic_columns = FALSE ) %>% - select(all_of(colnames(sample_data$hg38$maf))) + select(any_of(all_cols)) +}) +hg38_ashm <- hg38_ashm %>% + filter(Tumor_Sample_Barcode %in% sample_data$meta$Tumor_Sample_Barcode) %>% + mutate(Pipeline = "SLMS-3") -sample_data$hg38$ashm <- bind_rows( - sample_data$hg38$ashm, - trios_ashm_hg38 -) %>% distinct +grch37_ashm <- grch37_ashm %>% distinct() +hg38_ashm <- hg38_ashm %>% distinct() -setwd("~/my_dir/repos/GAMBLR.data/") +diag_summary_maf(grch37_ashm, "grch37_ashm (final, pre-merge into maf)") +diag_summary_maf(hg38_ashm, "hg38_ashm (final, pre-merge into maf)") + +print("done extracting aSHM mutations from GAMBLR.results") -# Add data from Reddy paper -reddy_original_maf <- read_tsv( - "inst/extdata/studies/reddy_original_variants_with_VAF.maf.gz" -) %>% -select(any_of(colnames(sample_data$grch37$maf))) +# ============================================================================ +# Phase 6: final maf assembly -- Publication rows + the consolidated SLMS-3 +# pull + the aSHM-region pull, per genome build. No separate ashm table/slot +# -- maf and ashm had identical schemas, independent mutation_id sequences, +# and independent write-time dedup, with no code anywhere querying ashm in +# isolation (confirmed: GAMBLR.open's only reference is get_ssm_from_db()'s +# now-removed include_ashm toggle, which always merged it onto maf anyway). +# Pooled here unconditionally instead -- one table, one dedup pass in +# write_mutations_db.R. No `Study =` in any mutate() anywhere in this +# script; cohort membership lives entirely in sample_data$sample_study. +# ============================================================================ sample_data$grch37$maf <- bind_rows( - sample_data$grch37$maf, - reddy_original_maf %>% - mutate( - Pipeline = "Publication", - Study = "Reddy" - ) + grch37_publication_rows, + slms3_grch37, + grch37_ashm ) -usethis::use_data( +sample_data$hg38$maf <- bind_rows( + hg38_publication_rows, + slms3_hg38, + hg38_ashm +) + +# The other money diagnostic: per-study coding-classified counts, mirroring +# what get_all_coding_ssm() showed missing on the GSC -- printed at build +# time now instead of only being discoverable in a separate later session. +coding_classes <- c("Missense_Mutation", "Nonsense_Mutation", "Frame_Shift_Del", + "Frame_Shift_Ins", "In_Frame_Del", "In_Frame_Ins", + "Splice_Site", "Splice_Region", "Nonstop_Mutation", + "Translation_Start_Site") +diag_summary_maf(sample_data$grch37$maf, "sample_data$grch37$maf (final, all Variant_Classifications)") +diag_summary_maf(sample_data$grch37$maf %>% filter(Variant_Classification %in% coding_classes), "sample_data$grch37$maf (final, coding-classified only)") +diag_summary_maf(sample_data$hg38$maf, "sample_data$hg38$maf (final, all Variant_Classifications)") +diag_summary_maf(sample_data$hg38$maf %>% filter(Variant_Classification %in% coding_classes), "sample_data$hg38$maf (final, coding-classified only)") + +print("done extracting all mutations in lymphoma genes with GAMBLR.results") + + +# --- Persist the assembled data ------------------------------------------- +# Previously this bundled the multi-GB `sample_data.rda`, and later a +# lightweight `sample_metadata` object shipped inside the package itself. +# That bundled-metadata object is retired: GAMBLR.open::get_gambl_metadata() +# (its only real consumer) now queries the sample_meta table in +# gambl_mutations.db directly instead, so metadata and mutation data always +# come from the same build with no separate reinstall-triggered refresh +# cadence to drift out of sync. This is built straight from the in-memory +# `sample_data` above, so there is no sample_data.rda round-trip. + +print("Starting sqlite build") +source("data-raw/write_mutations_db.R") +write_mutations_db( sample_data, - overwrite = TRUE, - compress = "xz" + out_db = "gambl_mutations.db", + source_desc = paste0("assemble_bundled_data.R @ ", format(Sys.Date())) ) +# During the transition you may still want the legacy monolithic object +# (e.g. until other GAMBLR packages are audited off sample_data). Flip to TRUE +# to also write data/sample_data.rda. +WRITE_LEGACY_SAMPLE_DATA <- FALSE +if (WRITE_LEGACY_SAMPLE_DATA) { + usethis::use_data(sample_data, overwrite = TRUE, compress = "xz") +} + library(data.tree) tree <- FromListSimple(sample_data) diff --git a/data-raw/build_mutations_db.R b/data-raw/build_mutations_db.R new file mode 100644 index 0000000..53f7262 --- /dev/null +++ b/data-raw/build_mutations_db.R @@ -0,0 +1,26 @@ +# build_mutations_db.R +# +# ONE-TIME adapter: build gambl_mutations.db from an EXISTING sample_data.rda, +# without re-running the (GSC-dependent) assembly. Use this only to convert a +# pre-existing bundled object. +# +# The ONGOING build writes the DB directly from source at the end of +# data-raw/assemble_bundled_data.R via the same write_mutations_db() helper, so +# there is no sample_data.rda round-trip in the normal pipeline. +# +# Usage: +# Rscript data-raw/build_mutations_db.R [SAMPLE_DATA_RDA] [OUT_DB] + +source("data-raw/write_mutations_db.R") + +args <- commandArgs(trailingOnly = TRUE) +rda <- if (length(args) >= 1) args[[1]] else "data/sample_data.rda" +out <- if (length(args) >= 2) args[[2]] else "gambl_mutations.db" + +message("Loading ", rda, " ...") +e <- new.env() +load(rda, envir = e) +sample_data <- get("sample_data", envir = e) + +write_mutations_db(sample_data, out_db = out, + source_desc = paste("one-time adapter from", rda)) diff --git a/data-raw/build_reference_db.R b/data-raw/build_reference_db.R new file mode 100644 index 0000000..3a3e62c --- /dev/null +++ b/data-raw/build_reference_db.R @@ -0,0 +1,156 @@ +# build_reference_db.R +# +# Proof-of-concept: build a single normalized SQLite reference database from the +# canonical LLMPP curated gene lists (NOT the versioned lymphoma_genes_*.rda +# objects). This one file is intended to become the shared backing store for +# GAMBLR.data accessors, the Lymphopedia wiki, and the Lymphopedia JSON/SQL API. +# +# Usage: +# Rscript data-raw/build_reference_db.R [LLMPP_CURATED_DIR] [OUT_DB] +# +# Design notes: +# * The per-entity curated TSVs have heterogeneous schemas. We normalize the +# common core into one long `gene_entity` table (one row per gene x entity) +# and stash the entity-specific columns as JSON in `extra` so nothing is lost. +# * "Which version" becomes a column (source_version), not a filename +# convention -- this collapses the lymphoma_genes_*_v0.1/v0.2/v_latest soup. + +suppressMessages({ + library(dplyr) + library(readr) + library(tidyr) + library(DBI) + library(jsonlite) +}) + +args <- commandArgs(trailingOnly = TRUE) +curated_dir <- ifelse(length(args) >= 1, args[[1]], + path.expand("~/git/LLMPP/resources/curated")) +out_db <- ifelse(length(args) >= 2, args[[2]], + file.path("inst", "extdata", "gambl_reference.db")) + +stopifnot(dir.exists(curated_dir)) +dir.create(dirname(out_db), recursive = TRUE, showWarnings = FALSE) + +# LLMPP source version -> stamped into every row so the data is self-describing +llmpp_version <- tryCatch( + system(paste("git -C", shQuote(dirname(dirname(curated_dir))), + "describe --always --dirty"), intern = TRUE), + error = function(e) NA_character_ +) +if (length(llmpp_version) == 0 || is.na(llmpp_version[1])) { + llmpp_version <- paste0("LLMPP-", format(Sys.Date())) +} else { + llmpp_version <- paste0("LLMPP@", llmpp_version[1]) +} + +# entity -> curated file +entity_files <- c( + DLBCL = "dlbcl_genes.tsv", + BL = "bl_genes.tsv", + FL = "fl_genes.tsv", + PMBL = "pmbl_genes.tsv", + MCL = "mcl_genes.tsv", + MZL = "mzl_genes.tsv" +) + +# pull a column by any of several candidate names (case-insensitive); NA if absent +pick <- function(df, candidates) { + hit <- names(df)[tolower(names(df)) %in% tolower(candidates)] + if (length(hit) == 0) return(rep(NA_character_, nrow(df))) + as.character(df[[hit[1]]]) +} + +core_cols <- c("gene", "tier", "ashm", "qc", "citekey", "pmid", "mutation_effect") + +harmonize_one <- function(entity, file) { + df <- suppressWarnings(read_tsv(file, col_types = cols(.default = "c"), + name_repair = "unique_quiet")) + core <- tibble( + gene = pick(df, c("Gene", "gene", "Hugo_Symbol")), + entity = entity, + tier = pick(df, c("Tier")), + ashm = pick(df, c("aSHM")), + qc = pick(df, c("QC")), + citekey = pick(df, c("citekey")), + pmid = pick(df, c("PMID")), + mutation_effect = pick(df, c("MutationEffect")) + ) |> filter(!is.na(gene), gene != "") + # keep everything else as JSON so entity-specific fields are not lost + used <- c("Gene", "gene", "Hugo_Symbol", "Tier", "aSHM", "QC", "citekey", + "PMID", "MutationEffect") + extra_df <- df[, !tolower(names(df)) %in% tolower(used), drop = FALSE] + extra_df <- extra_df[, names(extra_df) != "" & !grepl("^\\.\\.\\.", names(extra_df)), + drop = FALSE] + core$extra <- if (ncol(extra_df) == 0) NA_character_ else + vapply(seq_len(nrow(core)), function(i) + toJSON(as.list(extra_df[i, , drop = FALSE]), auto_unbox = TRUE, na = "null"), + character(1)) + core +} + +# symbol -> ensembl map, reused from GAMBLR.data's own gencode_to_symbol object, +# so gene_format = "ensembl" keeps working from the DB. One ensembl id per symbol. +sym2ens <- local({ + e <- new.env(); load(file.path("data", "gencode_to_symbol.rda"), envir = e) + get("gencode_to_symbol", envir = e) |> + dplyr::filter(!is.na(hgnc_symbol), hgnc_symbol != "") |> + dplyr::distinct(hgnc_symbol, .keep_all = TRUE) |> + dplyr::select(gene = hgnc_symbol, ensembl_gene_id) +}) + +gene_entity <- bind_rows( + Map(function(e, f) harmonize_one(e, file.path(curated_dir, f)), + names(entity_files), entity_files) +) |> + left_join(sym2ens, by = "gene") |> + mutate( + tier = suppressWarnings(as.integer(tier)), + source_version = llmpp_version + ) |> + relocate(ensembl_gene_id, .after = gene) |> + relocate(source_version, .after = last_col()) + +# --- second reference table: aSHM regions (both genome builds) ---------------- +read_ashm <- function(file, build) { + df <- suppressWarnings(read_tsv(file, col_types = cols(.default = "c"), + name_repair = "unique_quiet")) + tibble( + gene = pick(df, c("gene")), + chrom = pick(df, c("chr_name", "chrom", "chr")), + start = as.integer(df[[2]]), + end = as.integer(df[[3]]), + region = pick(df, c("region")), + regulatory_comment = pick(df, c("regulatory_comment")), + build = build, + source_version = llmpp_version + ) |> filter(!is.na(gene)) +} +ashm_dir <- file.path(curated_dir, "aSHM") +ashm_regions <- bind_rows( + read_ashm(file.path(ashm_dir, "somatic_hypermutation_locations_GRCh37.txt"), "grch37"), + read_ashm(file.path(ashm_dir, "somatic_hypermutation_locations_hg38.txt"), "hg38") +) + +# --- write the database ------------------------------------------------------- +if (file.exists(out_db)) file.remove(out_db) +con <- dbConnect(RSQLite::SQLite(), out_db) +on.exit(dbDisconnect(con)) + +dbWriteTable(con, "gene_entity", as.data.frame(gene_entity), overwrite = TRUE) +dbWriteTable(con, "ashm_regions", as.data.frame(ashm_regions), overwrite = TRUE) +dbWriteTable(con, "build_info", data.frame( + key = c("source", "source_version", "built_at", "builder"), + value = c("LLMPP resources/curated", llmpp_version, + format(Sys.time(), tz = "UTC", usetz = TRUE), + "data-raw/build_reference_db.R") +), overwrite = TRUE) + +dbExecute(con, "CREATE INDEX idx_ge_gene ON gene_entity(gene)") +dbExecute(con, "CREATE INDEX idx_ge_entity ON gene_entity(entity)") +dbExecute(con, "CREATE INDEX idx_ge_tier ON gene_entity(tier)") +dbExecute(con, "CREATE INDEX idx_ashm_gene ON ashm_regions(gene)") + +cat(sprintf("Wrote %s\n gene_entity : %d rows (%d genes, %d entities)\n ashm_regions: %d rows\n version : %s\n", + out_db, nrow(gene_entity), dplyr::n_distinct(gene_entity$gene), + dplyr::n_distinct(gene_entity$entity), nrow(ashm_regions), llmpp_version)) diff --git a/data-raw/compare_bundle_changes.R b/data-raw/compare_bundle_changes.R new file mode 100644 index 0000000..4ee4150 --- /dev/null +++ b/data-raw/compare_bundle_changes.R @@ -0,0 +1,172 @@ +# compare_bundle_changes.R +# +# Compares the OLD bundled sample_data.rda against a NEW gambl_mutations.db +# build at the level of individual events (not just row totals), for each of +# three data types: SNV (maf), CNV (seg), SV (bedpe). For each type: +# - counts, per sample_id, how many events are NEW (gained: present in the +# new build but not the old) and how many are LOST (missing: present in +# the old build but not the new) +# - writes one example gained-event row per affected sample to +# /{snv,cnv,sv}_gained_examples.log, and one example lost-event +# row per affected sample to /{snv,cnv,sv}_lost_examples.log -- +# for quick manual inspection without dumping every gained/lost row. At +# most one row per sample that has >=1 gained (or lost) event, so at most +# as many rows as there are samples -- in practice fewer, since not every +# sample is affected. +# - also writes the full per-sample gained/lost counts to +# /{snv,cnv,sv}_persample_counts.tsv for further analysis. +# +# "Same event" is determined by an exact match on a type-specific key (sample +# id + genome_build + genomic coordinates), NOT by row position or count -- +# two rows only count as "the same" if every key column matches exactly. +# +# Usage: Rscript data-raw/compare_bundle_changes.R [gambl_mutations.db] [old sample_data.rda] [output dir] + +suppressMessages({library(DBI); library(dplyr)}) + +args <- commandArgs(trailingOnly = TRUE) +db <- if (length(args) >= 1) args[[1]] else "gambl_mutations.db" +old_rda <- if (length(args) >= 2) args[[2]] else "data/sample_data.rda" +outdir <- if (length(args) >= 3) args[[3]] else "." +dir.create(outdir, recursive = TRUE, showWarnings = FALSE) + +con <- DBI::dbConnect(RSQLite::SQLite(), db) +on.exit(DBI::dbDisconnect(con), add = TRUE) + +e <- new.env() +load(old_rda, envir = e) +old_sd <- get("sample_data", envir = e) + +# Generic per-sample gained/lost comparison for one data type. +# old_df, new_df : data frames (already build-stamped with genome_build if +# the type spans both builds) +# sample_col : name of the sample-id column for this type +# key_cols : columns (including sample_col and genome_build, if +# present) that together define a unique "event" -- an +# exact match on all of these means "the same event" +compare_events <- function(old_df, new_df, sample_col, key_cols, label) { + key_of <- function(df) { + # ":::" separator avoids ambiguous concatenation, e.g. chrom="1", + # start=23456789 vs chrom="12",start=3456789 both -> "123456789" if + # pasted with no separator at all + do.call(paste, c(as.list(select(df, all_of(key_cols))), sep = ":::")) + } + old_df$.key <- key_of(old_df) + new_df$.key <- key_of(new_df) + + gained <- new_df %>% filter(!.key %in% old_df$.key) + lost <- old_df %>% filter(!.key %in% new_df$.key) + retained <- old_df %>% filter(.key %in% new_df$.key) + + gained_per_sample <- gained %>% count(.data[[sample_col]], name = "n_gained") + lost_per_sample <- lost %>% count(.data[[sample_col]], name = "n_lost") + retained_per_sample <- retained %>% count(.data[[sample_col]], name = "n_retained") + + # n_retained (present in both old and new) is joined in separately -- it's + # looked up only for the samples already in per_sample (those with >=1 + # gained or lost event), not full_join'd, since retained_per_sample on its + # own covers essentially every sample including untouched ones. + per_sample <- full_join(gained_per_sample, lost_per_sample, by = sample_col) %>% + left_join(retained_per_sample, by = sample_col) %>% + mutate( + n_gained = coalesce(n_gained, 0L), + n_lost = coalesce(n_lost, 0L), + n_retained = coalesce(n_retained, 0L) + ) %>% + arrange(desc(n_gained + n_lost)) + + cat(sprintf( + "\n=== %s: %d gained rows, %d lost rows (%d samples w/ gains, %d samples w/ losses) ===\n", + label, nrow(gained), nrow(lost), + sum(per_sample$n_gained > 0), sum(per_sample$n_lost > 0) + )) + print(head(per_sample, 10)) + + list(per_sample = per_sample, gained = gained, lost = lost) +} + +# One example row per affected sample, from either the gained or the lost +# data frame returned by compare_events() -- whichever is passed in. +write_examples <- function(event_df, sample_col, path) { + if (nrow(event_df) == 0) { + message("no rows -- skipping ", path) + return(invisible(NULL)) + } + examples <- event_df %>% + group_by(.data[[sample_col]]) %>% + slice_head(n = 1) %>% + ungroup() %>% + select(-.key) + write.table(examples, path, sep = "\t", quote = FALSE, row.names = FALSE) + message("wrote ", nrow(examples), " example rows to ", path) +} + +write_counts <- function(per_sample, path) { + write.table(per_sample, path, sep = "\t", quote = FALSE, row.names = FALSE) + message("wrote per-sample counts to ", path) +} + +# --- SNV: maf (which now also absorbs what used to be a separate ashm pull, +# deduplicated against it at write time -- see write_mutations_db.R) -------- +# .source records which table/pull each row came from on the OLD side so +# gained/lost examples can still be traced back to the right code path. +# The new side no longer has that distinction to make (maf and ashm were +# merged into one table with no retained provenance column -- deliberately; +# nothing consumed it), so .source is uniformly "maf" there. Not part of the +# comparison key either way: a row that moved from one source to the other +# between builds still counts as retained if its sample/position match. +new_snv <- dbGetQuery(con, "SELECT * FROM maf") %>% mutate(.source = "maf") +old_snv <- bind_rows( + old_sd$grch37$maf %>% mutate(genome_build = "grch37", .source = "maf"), + old_sd$hg38$maf %>% mutate(genome_build = "hg38", .source = "maf"), + old_sd$grch37$ashm %>% mutate(genome_build = "grch37", .source = "ashm"), + old_sd$hg38$ashm %>% mutate(genome_build = "hg38", .source = "ashm") +) + +# Pipeline == "strelka" (Arthur's old, unfiltered raw-flat-file dump) is +# intentionally and entirely gone from the new bundle -- that source file +# isn't read anymore. Confirmed (twice: 08-15460, DO52686) that its +# Tumor_Sample_Barcode values are the raw file's own non-standard IDs, not +# GAMBL's real sample_id, so every one of these rows is expected to show up +# as "lost" regardless of anything else in the pipeline. Excluded here so +# that expected, already-understood disappearance doesn't drown out +# genuinely new findings on every run. +old_snv <- old_snv %>% filter(is.na(Pipeline) | Pipeline != "strelka") + +snv_key <- c("Tumor_Sample_Barcode", "genome_build", "Chromosome", "Start_Position", "End_Position") +snv_cmp <- compare_events(old_snv, new_snv, "Tumor_Sample_Barcode", snv_key, "SNV (maf, incl. aSHM-region pull)") +write_examples(snv_cmp$gained, "Tumor_Sample_Barcode", file.path(outdir, "snv_gained_examples.log")) +write_examples(snv_cmp$lost, "Tumor_Sample_Barcode", file.path(outdir, "snv_lost_examples.log")) +write_counts(snv_cmp$per_sample, file.path(outdir, "snv_persample_counts.tsv")) + +# --- CNV: seg ----------------------------------------------------------------- +new_cnv <- dbGetQuery(con, "SELECT * FROM seg") +old_cnv <- bind_rows( + old_sd$grch37$seg %>% mutate(genome_build = "grch37"), + old_sd$hg38$seg %>% mutate(genome_build = "hg38") +) +cnv_key <- c("ID", "genome_build", "chrom", "start", "end") +cnv_cmp <- compare_events(old_cnv, new_cnv, "ID", cnv_key, "CNV (seg)") +write_examples(cnv_cmp$gained, "ID", file.path(outdir, "cnv_gained_examples.log")) +write_examples(cnv_cmp$lost, "ID", file.path(outdir, "cnv_lost_examples.log")) +write_counts(cnv_cmp$per_sample, file.path(outdir, "cnv_persample_counts.tsv")) + +# --- SV: bedpe ------------------------------------------------------------------ +new_sv <- dbGetQuery(con, "SELECT * FROM bedpe") +old_sv <- bind_rows( + old_sd$grch37$bedpe %>% mutate(genome_build = "grch37"), + old_sd$hg38$bedpe %>% mutate(genome_build = "hg38") +) +sv_key <- c("tumour_sample_id", "genome_build", "CHROM_A", "START_A", "CHROM_B", "START_B") +sv_cmp <- compare_events(old_sv, new_sv, "tumour_sample_id", sv_key, "SV (bedpe)") +write_examples(sv_cmp$gained, "tumour_sample_id", file.path(outdir, "sv_gained_examples.log")) +write_examples(sv_cmp$lost, "tumour_sample_id", file.path(outdir, "sv_lost_examples.log")) +write_counts(sv_cmp$per_sample, file.path(outdir, "sv_persample_counts.tsv")) + +cat("\n=== summary ===\n") +cat(sprintf("SNV: %d samples with gains, %d samples with losses\n", + sum(snv_cmp$per_sample$n_gained > 0), sum(snv_cmp$per_sample$n_lost > 0))) +cat(sprintf("CNV: %d samples with gains, %d samples with losses\n", + sum(cnv_cmp$per_sample$n_gained > 0), sum(cnv_cmp$per_sample$n_lost > 0))) +cat(sprintf("SV: %d samples with gains, %d samples with losses\n", + sum(sv_cmp$per_sample$n_gained > 0), sum(sv_cmp$per_sample$n_lost > 0))) diff --git a/data-raw/public_sample_meta_columns.txt b/data-raw/public_sample_meta_columns.txt new file mode 100644 index 0000000..3501eb6 --- /dev/null +++ b/data-raw/public_sample_meta_columns.txt @@ -0,0 +1,61 @@ +# Columns from live GAMBLR.results::get_gambl_metadata() to carry into the +# public sample_meta table -- and therefore into the publicly-distributed +# gambl_mutations.db (a GitHub Release asset anyone can download and query +# directly with SQL), and into GAMBLR.open::get_gambl_metadata()'s output. +# +# One column name per line. Lines starting with # are comments; blank +# lines are ignored. Add a line to include a new column in future builds -- +# no code changes needed, just edit this file and re-run +# assemble_bundled_data.R. +# +# sample_id, seq_type, study, and reference_PMID are always included +# separately (required for the pipeline's own internal joins/filters) and +# don't need to be listed here. +# +# IMPORTANT: this is a whitelist for a PUBLIC, downloadable database file, +# not just a filter on one R function -- anyone can query gambl_mutations.db +# directly, bypassing GAMBLR.open entirely. Only add a column here if it's +# genuinely safe to distribute publicly (e.g. NOT biopsy_id, data_path, +# fastq_data_path, fastq_link_name, link_name, library_id, res_id, +# unix_group, or anything else that's an internal tracking/file-path +# identifier rather than sample-level scientific metadata). When in doubt, +# ask before adding. +# +# Confirmed list, 2026-07-27. age_group and hiv_status were deliberately +# left out: neither is actually populated by the current live +# get_gambl_metadata() (age_group is only computed by the old, unused +# og_get_gambl_metadata(), gated behind an optional with_outcomes +# parameter; hiv_status was never implemented anywhere -- its own roxygen +# entry says {TODO}). Add them here once/if that derivation exists in the +# current function. +patient_id +pathology +biopsy_res +genome_build +pairing_status +Tumor_Sample_Barcode +cohort +COO_consensus +DHITsig_consensus +EBV_status_inf +ffpe_or_frozen +fl_grade +bcl2_ba +bcl2_cn +bcl6_ba +bcl6_cn +myc_ba +myc_cn +lymphgen +lymphgen_cnv_noA53 +lymphgen_no_cnv +lymphgen_with_cnv +lymphgen_wright +normal_sample_id +sex +time_point +TotalDuplicatedreads +TotalReads +TotalUniquelyMapped +TotalUnmappedreads +transformation diff --git a/data-raw/release_mutations_db.R b/data-raw/release_mutations_db.R new file mode 100644 index 0000000..0072152 --- /dev/null +++ b/data-raw/release_mutations_db.R @@ -0,0 +1,94 @@ +# release_mutations_db.R +# +# Admin-only: uploads a freshly-built gambl_mutations.db as a GitHub Release +# asset, so GAMBLR.data::gambl_mutations_db() can auto-download it for new +# installs (see R/download_gambl_mutations_db.R). Run this AFTER +# assemble_bundled_data.R / write_mutations_db() has produced a fresh +# gambl_mutations.db, from the GAMBLR.data repo root. This script only +# uploads -- it does not build the database itself. +# +# Requires: +# - the piggyback package: install.packages("piggyback") +# - a GitHub PAT with write access to this repo, either already configured +# (gitcreds::gitcreds_set()) or set via the GITHUB_PAT/GITHUB_TOKEN +# environment variable +# +# The release tag defaults to "data-v", matching what +# download_gambl_mutations_db() looks for by default. Keep these in lockstep: +# bump DESCRIPTION's Version whenever the DB schema changes (as it did +# several times in the sample_study/variant_pipeline refactor), and cut a +# matching release here, so a given package version always downloads data +# built for its own schema rather than whatever the newest release happens +# to be. +# +# Usage: Rscript data-raw/release_mutations_db.R [path to .db] [tag] + +args <- commandArgs(trailingOnly = TRUE) +db_path <- if (length(args) >= 1) args[[1]] else "gambl_mutations.db" +pkg_version <- read.dcf("DESCRIPTION")[, "Version"] +tag <- if (length(args) >= 2) args[[2]] else paste0("data-v", pkg_version) +repo <- "morinlab/GAMBLR.data" + +stopifnot(file.exists(db_path)) +if (!requireNamespace("piggyback", quietly = TRUE)) { + stop("Install piggyback first: install.packages(\"piggyback\")") +} + +# The tag is derived from DESCRIPTION's Version on disk, but +# download_gambl_mutations_db() on every *other* install resolves the same +# tag from utils::packageVersion("GAMBLR.data") -- i.e. whatever's actually +# committed and installed elsewhere. If the version bump used to build this +# tag isn't committed and pushed, every install (including a future +# reinstall of this exact checkout) will keep resolving the OLD version and +# never find this release at all. +git_diff_status <- suppressWarnings(system2("git", c("diff", "--quiet", "--", "DESCRIPTION"), stdout = FALSE, stderr = FALSE)) +git_staged_status <- suppressWarnings(system2("git", c("diff", "--quiet", "--cached", "--", "DESCRIPTION"), stdout = FALSE, stderr = FALSE)) +if (git_diff_status != 0 || git_staged_status != 0) { + stop("DESCRIPTION has uncommitted changes. Commit and push the version bump ", + "to \"", pkg_version, "\" BEFORE running this script, or every install ", + "(including your own) will keep resolving a different version and ", + "never find the \"", tag, "\" release you're about to create.", + call. = FALSE) +} +head_tag_check <- system2("git", c("log", "-1", "--format=%H"), stdout = TRUE, stderr = FALSE) +message("Releasing from commit ", head_tag_check, " -- make sure this has been pushed to origin.") + +message(sprintf("Releasing %s to %s @ %s (%.0f MB)", + db_path, repo, tag, file.info(db_path)$size / 1024^2)) + +# create the release/tag if it doesn't already exist +existing_releases <- piggyback::pb_releases(repo = repo) +if (!tag %in% existing_releases$tag_name) { + message("Creating new release tag: ", tag) + piggyback::pb_new_release(repo = repo, tag = tag) +} + +# GitHub's list-releases endpoint can lag a moment after pb_new_release() +# creates a tag -- pb_upload() re-fetches that list to confirm the tag +# exists, and can spuriously fail with "Release not found" if it queries +# before GitHub's side catches up. piggyback's own pb_upload() has a +# Sys.sleep()-and-retry fallback for this, but only when interactive(), +# which Rscript never is -- so retry it here instead. +upload_with_retry <- function(max_tries = 5) { + for (i in seq_len(max_tries)) { + ok <- tryCatch({ + piggyback::pb_upload(db_path, repo = repo, tag = tag) + TRUE + }, error = function(e) { + if (grepl("not found", conditionMessage(e), ignore.case = TRUE) && i < max_tries) { + message("Release not yet visible via GitHub's API (attempt ", i, "/", max_tries, + "), retrying in ", 2 * i, "s...") + Sys.sleep(2 * i) + FALSE + } else { + stop(e) + } + }) + if (isTRUE(ok)) return(invisible()) + } +} +upload_with_retry() +message(sprintf( + "Done. Users on GAMBLR.data v%s will now auto-download gambl_mutations.db from this release.", + pkg_version +)) diff --git a/data-raw/test_gambl_db.R b/data-raw/test_gambl_db.R new file mode 100644 index 0000000..079c885 --- /dev/null +++ b/data-raw/test_gambl_db.R @@ -0,0 +1,129 @@ +# test_gambl_db.R +# +# Simple, dependency-light sanity checks for a freshly built gambl_mutations.db. +# Run after building the DB (on the GSC side) to confirm the output is usable by +# the GAMBLR.open accessors. +# +# Usage: +# Rscript data-raw/test_gambl_db.R [OUT_DB] [SAMPLE_DATA_RDA] +# +# Exits non-zero if any critical check fails. The optional SAMPLE_DATA_RDA +# enables a regression check (DB row counts vs the in-memory frames). + +suppressMessages({ library(DBI); library(RSQLite) }) + +args <- commandArgs(trailingOnly = TRUE) +db <- if (length(args) >= 1) args[[1]] else "gambl_mutations.db" +rda <- if (length(args) >= 2) args[[2]] else "data/sample_data.rda" + +fails <- 0L +check <- function(cond, msg) { + status <- if (isTRUE(cond)) "PASS" else { fails <<- fails + 1L; "FAIL" } + cat(sprintf(" [%s] %s\n", status, msg)) +} + +stopifnot(file.exists(db)) +con <- dbConnect(SQLite(), db, flags = SQLITE_RO) +on.exit(dbDisconnect(con)) +tbls <- dbListTables(con) +n <- function(t) dbGetQuery(con, sprintf("SELECT COUNT(*) n FROM %s", t))$n +cols <- function(t) dbListFields(con, t) +builds_in <- function(t) sort(dbGetQuery(con, sprintf("SELECT DISTINCT genome_build g FROM %s", t))$g) + +cat("== tables present ==\n") +for (t in c("maf","seg","bedpe","sample_meta","sample_study","variant_pipeline","build_info")) + check(t %in% tbls, sprintf("table %s exists", t)) +check(!("ashm" %in% tbls), "table ashm does NOT exist (merged into maf)") + +cat("\n== row counts (all > 0) ==\n") +for (t in c("maf","seg","bedpe","sample_meta","sample_study","variant_pipeline")) { + cnt <- if (t %in% tbls) n(t) else 0 + check(cnt > 0, sprintf("%-11s has %d rows", t, cnt)) +} + +cat("\n== both genome builds present ==\n") +for (t in c("maf","seg","bedpe")) + if (t %in% tbls) + check(all(c("grch37","hg38") %in% builds_in(t)), + sprintf("%-6s has builds: %s", t, paste(builds_in(t), collapse=", "))) + +cat("\n== required columns present (what the accessors filter on) ==\n") +req <- list( + maf = c("mutation_id","Chromosome","Start_Position","End_Position","Tumor_Sample_Barcode", + "Pipeline","Hugo_Symbol","Variant_Classification","t_alt_count","genome_build"), + seg = c("ID","chrom","start","end","CN","genome_build"), + bedpe = c("tumour_sample_id","CHROM_A","START_A","CHROM_B","START_B","VAF_tumour","SCORE","FILTER","genome_build"), + sample_meta = c("sample_id","Tumor_Sample_Barcode","seq_type","study"), + sample_study = c("sample_id","study","study_id","reference_PMID"), + variant_pipeline = c("mutation_id","Pipeline") +) +for (t in names(req)) if (t %in% tbls) { + missing <- setdiff(req[[t]], cols(t)) + check(length(missing) == 0, + sprintf("%-11s columns present%s", t, + if (length(missing)) paste0(" (MISSING: ", paste(missing, collapse=", "), ")") else "")) +} + +cat("\n== indexes present ==\n") +# NB: idx_maf_gene was asserted here previously but was never actually +# created by write_mutations_db.R (maf has no Hugo_Symbol index -- gene +# queries resolve to a region first, see gambl_mutations_db()'s docs), so +# this check always failed regardless of any other change. Fixed here to +# assert indexes that actually exist. +idx <- dbGetQuery(con, "SELECT name FROM sqlite_master WHERE type='index'")$name +for (i in c("idx_maf_pos","idx_maf_sample","idx_maf_pipe","idx_maf_mutation_id", + "idx_seg_sample", + "idx_bedpe_sample","idx_bedpe_pos_a","idx_bedpe_pos_b", + "idx_study_sample","idx_study_study", + "idx_vp_pipe","idx_vp_mutation_id")) + check(i %in% idx, sprintf("index %s", i)) + +cat("\n== pipelines present ==\n") +# Pipeline is normalized to lowercase at write time (see write_mutations_db.R) +# so it can be matched with a plain, indexed equality instead of every query +# needing to wrap the column in LOWER()/tolower(). +pipes <- dbGetQuery(con, "SELECT DISTINCT Pipeline p FROM maf")$p +for (p in c("slms-3","publication")) + check(p %in% pipes, sprintf("Pipeline '%s' present (all: %s)", p, paste(pipes, collapse=", "))) + +cat("\n== spot query: MYC locus, grch37, slms-3 returns variants ==\n") +myc <- dbGetQuery(con, "SELECT COUNT(*) n FROM maf WHERE genome_build='grch37' + AND Chromosome='8' AND Start_Position>128723128 AND Start_Position<128774067 + AND Pipeline='slms-3'")$n +check(myc > 0, sprintf("MYC slms-3 variants found: %d", myc)) + +cat("\n== build_info counts match actual table counts ==\n") +if ("build_info" %in% tbls) { + bi <- dbGetQuery(con, "SELECT key, value FROM build_info") + getbi <- function(k) as.numeric(bi$value[bi$key == k]) + for (t in c("maf","seg","bedpe")) + check(isTRUE(getbi(paste0("n_", t)) == n(t)), + sprintf("build_info n_%s == COUNT(%s)", t, t)) + check(isTRUE(getbi("n_variant_pipeline") == n("variant_pipeline")), + sprintf("build_info n_variant_pipeline == COUNT(variant_pipeline)")) +} + +cat("\n== regression vs sample_data.rda (optional) ==\n") +if (file.exists(rda)) { + e <- new.env(); load(rda, envir = e); sd <- get("sample_data", envir = e) + rda_rows <- function(t) sum(vapply(c("grch37","hg38"), + function(b) { x <- sd[[b]][[t]]; if (is.null(x)) 0L else nrow(x) }, integer(1))) + for (t in c("seg","bedpe")) { + expect <- rda_rows(t) + check(expect == n(t), sprintf("%-6s DB rows (%d) == sample_data rows (%d)", t, n(t), expect)) + } + # maf no longer has an exact-equality invariant against the legacy .rda: + # the DB's maf now also absorbs what used to be a separate ashm pull, + # deduplicated against it at write time (see write_mutations_db.R), so the + # new count is expected to be somewhere between the old maf alone (nothing + # was ever removed from maf's own content) and old maf + old ashm combined + # (the ceiling, only reached if the two pulls never overlapped at all). + old_maf <- rda_rows("maf") + old_ashm <- rda_rows("ashm") + check(n("maf") >= old_maf && n("maf") <= old_maf + old_ashm, + sprintf("maf DB rows (%d) within [old maf (%d), old maf+ashm (%d)]", + n("maf"), old_maf, old_maf + old_ashm)) +} else cat(" [skip] no", rda, "(fresh GSC build has no legacy .rda — expected)\n") + +cat(sprintf("\n==== %s: %d check(s) failed ====\n", if (fails==0) "ALL PASSED" else "FAILURES", fails)) +if (fails > 0) quit(status = 1) diff --git a/data-raw/write_mutations_db.R b/data-raw/write_mutations_db.R new file mode 100644 index 0000000..3020c9c --- /dev/null +++ b/data-raw/write_mutations_db.R @@ -0,0 +1,201 @@ +# write_mutations_db.R +# +# Reusable build-time helper: given the in-memory `sample_data` list (as +# assembled by data-raw/assemble_bundled_data.R), write the normalized +# `gambl_mutations.db` SQLite database with indexes matching the GAMBLR.open +# accessors. This is the ONGOING path that replaces bundling `sample_data.rda`: +# assemble_bundled_data.R -> sample_data (in memory) -> write_mutations_db() +# +# `source()` this file (it defines the function; it does not run anything). +# Requires the DBI and RSQLite packages (Imports of GAMBLR.data). + +write_mutations_db <- function(sample_data, + out_db = "gambl_mutations.db", + builds = c("grch37", "hg38"), + source_desc = "assemble_bundled_data.R") { + stopifnot(is.list(sample_data), !is.null(sample_data$meta)) + + if (file.exists(out_db)) file.remove(out_db) + con <- DBI::dbConnect(RSQLite::SQLite(), out_db) + on.exit(DBI::dbDisconnect(con)) + # speed up the bulk load + DBI::dbExecute(con, "PRAGMA journal_mode = OFF") + DBI::dbExecute(con, "PRAGMA synchronous = OFF") + + # sample metadata (single table, no genome build). reference_PMID is + # dropped here -- it's single-valued (one row per sample) and can't + # correctly represent a sample belonging to more than one study (each with + # its own PMID); sample_study is the source of truth for that now. + sample_meta_out <- as.data.frame(sample_data$meta) + sample_meta_out$reference_PMID <- NULL + DBI::dbWriteTable(con, "sample_meta", sample_meta_out, overwrite = TRUE) + + # sample_study: many-to-many bridge table (sample_id, study) replacing the + # old inline Study column on maf rows -- guarded so older sample_data + # objects that predate this table (e.g. via the legacy .rda adapter) don't + # hard-crash a build that simply won't have cohort-membership data. + n_sample_study <- 0L + if (!is.null(sample_data$sample_study) && nrow(sample_data$sample_study)) { + DBI::dbWriteTable(con, "sample_study", as.data.frame(sample_data$sample_study), + overwrite = TRUE) + n_sample_study <- nrow(sample_data$sample_study) + } + + # Natural mutation-call key, reused from GAMBLR.open::get_ssm_by_region() + # (which already applies the same distinct() at read time) as a write-time + # safety net: no cohort-specific pull block should ever again be able to + # introduce duplicate mutation rows for a sample that's pulled more than + # once, regardless of the reason. + # maf_seq_type is part of the key so a genome-called and capture-called + # mutation for the same sample that happen to land on the same + # position/allele stay two distinct rows instead of collapsing into one -- + # dplyr::distinct() treats NA (publication-pipeline rows, which don't carry + # per-row seq_type) as matching other NAs, so this doesn't change dedup + # behaviour for rows that predate maf_seq_type, only adds separation where + # it's actually known. + dedup_keys <- list( + maf = c("Tumor_Sample_Barcode", "Chromosome", "Start_Position", "End_Position", "Tumor_Seq_Allele2", "maf_seq_type") + ) + + # variant_pipeline: many-to-many bridge table (mutation_id, Pipeline), same + # pattern as sample_study (sample_id, study). A single Pipeline column on + # maf can only hold one value per variant, but the same real mutation + # can legitimately be produced by more than one pipeline (e.g. a cohort's + # own published/curated maf and our independent SLMS-3 recall both calling + # the same position) -- when that happens the dedup below collapses them + # to one row, silently making the call invisible to a tool_name-filtered + # query for the other pipeline even though the variant itself is still + # fully present in the table. Captured here, before dedup, so no + # pipeline's claim to a variant is ever lost regardless of which row + # "wins" as maf's single representative copy. + # + # mutation_id is a surrogate integer key assigned to each deduped maf + # row (unique across both genome builds), and variant_pipeline references + # it as a plain foreign key instead of repeating the 5-column natural key + # on every row. Matching pre-dedup pipeline claims back to their surviving + # row still requires the natural key -- that cost doesn't disappear, it + # just moves here, to a one-time write-time join, instead of every + # read-time query in get_ssm_from_db() having to join on a 5-column + # composite (two of them text) instead of a single indexed integer column. + # (No elem column here -- that only ever existed to disambiguate maf's and + # ashm's independent mutation_id spaces; maf and ashm were merged into one + # table with one mutation_id sequence, so there's nothing left to + # disambiguate.) + variant_pipeline_rows <- list() + + # genome-build-stamped frames; append per build to keep peak memory low + write_element <- function(elem) { + first <- TRUE + total <- 0L + dk <- dedup_keys[[elem]] + next_id <- 1L + for (b in builds) { + x <- sample_data[[b]][[elem]] + if (is.null(x) || !nrow(x)) next + x <- as.data.frame(x) + x$genome_build <- b + # Normalized to lowercase here, once, at the write boundary -- so + # every reader (get_ssm_from_db(), raw SQL) can match Pipeline with a + # plain equality against a plain index, instead of every query needing + # to wrap the column in LOWER()/tolower() (which a normal b-tree index + # can't be used to satisfy). + if ("Pipeline" %in% names(x)) x$Pipeline <- tolower(x$Pipeline) + + if (!is.null(dk)) { + before <- nrow(x) + deduped <- dplyr::distinct(x, dplyr::across(dplyr::all_of(dk)), .keep_all = TRUE) + if (nrow(deduped) < before) { + message(sprintf("[write_mutations_db] %s/%s: dropped %d duplicate row(s)", + elem, b, before - nrow(deduped))) + } + deduped$mutation_id <- seq.int(next_id, length.out = nrow(deduped)) + next_id <- next_id + nrow(deduped) + + if ("Pipeline" %in% names(x)) { + id_lookup <- deduped %>% dplyr::select(dplyr::all_of(dk), mutation_id) + variant_pipeline_rows[[length(variant_pipeline_rows) + 1]] <<- x %>% + dplyr::select(dplyr::all_of(dk), Pipeline) %>% + dplyr::distinct() %>% + dplyr::inner_join(id_lookup, by = dk) %>% + dplyr::select(mutation_id, Pipeline) + } + x <- deduped + } + DBI::dbWriteTable(con, elem, x, append = !first, overwrite = first) + total <- total + nrow(x) + first <- FALSE + } + total + } + counts <- vapply(c("maf", "seg", "bedpe"), write_element, integer(1)) + + n_variant_pipeline <- 0L + if (length(variant_pipeline_rows)) { + variant_pipeline <- dplyr::bind_rows(variant_pipeline_rows) + DBI::dbWriteTable(con, "variant_pipeline", as.data.frame(variant_pipeline), overwrite = TRUE) + n_variant_pipeline <- nrow(variant_pipeline) + } + + # indexes mirroring how the GAMBLR.open accessors query the data. + # wrapped in try() so a build missing an optional table/column is non-fatal. + # + # idx_maf_pipe/idx_vp_pipe are plain indexes on Pipeline -- safe because + # Pipeline is normalized to lowercase above, at write time, so readers + # (get_ssm_from_db(), raw SQL) can match it with a plain equality instead + # of wrapping the column in LOWER()/tolower(), which a normal b-tree index + # can't be used to satisfy. idx_*_mutation_id give the variant_pipeline + # semi-join (in get_ssm_from_db()'s tool_name resolution) an index to + # actually use on both sides of the join -- without one, that join has to + # scan maf in full for every query with a non-NULL tool_name (the + # default), which is most of them. + idx <- c( + "CREATE INDEX idx_maf_pos ON maf(genome_build, Chromosome, Start_Position)", + "CREATE INDEX idx_maf_sample ON maf(Tumor_Sample_Barcode)", + "CREATE INDEX idx_maf_pipe ON maf(Pipeline)", + "CREATE INDEX idx_maf_seqtype ON maf(maf_seq_type)", + "CREATE INDEX idx_maf_mutation_id ON maf(mutation_id)", + "CREATE INDEX idx_seg_sample ON seg(genome_build, ID)", + "CREATE INDEX idx_seg_pos ON seg(genome_build, chrom, start)", + # bedpe previously only had an index on tumour_sample_id alone -- every + # get_sv_from_db()/get_manta_sv() call also filters genome_build first, + # 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 regardless of how selective the actual filters were. + # tumour_sample_id leads its own index (more selective than genome_build + # alone) with genome_build appended so sample-restricted queries resolve + # in one lookup; separate genome_build-led indexes on each breakpoint end + # let SQLite's OR-optimization use a different index per side of a + # region search. + "CREATE INDEX idx_bedpe_sample ON bedpe(tumour_sample_id, genome_build)", + "CREATE INDEX idx_bedpe_pos_a ON bedpe(genome_build, CHROM_A, START_A)", + "CREATE INDEX idx_bedpe_pos_b ON bedpe(genome_build, CHROM_B, START_B)", + "CREATE INDEX idx_meta_sample ON sample_meta(sample_id)", + "CREATE INDEX idx_meta_barcode ON sample_meta(Tumor_Sample_Barcode)", + "CREATE INDEX idx_study_sample ON sample_study(sample_id)", + "CREATE INDEX idx_study_study ON sample_study(study)", + "CREATE INDEX idx_vp_pipe ON variant_pipeline(Pipeline)", + "CREATE INDEX idx_vp_mutation_id ON variant_pipeline(mutation_id)" + ) + for (stmt in idx) try(DBI::dbExecute(con, stmt), silent = TRUE) + + # self-describing provenance / expected-counts table (used by the tests) + DBI::dbWriteTable(con, "build_info", data.frame( + key = c("source", "built_at", "builder", + "n_maf", "n_seg", "n_bedpe", "n_samples", "n_sample_study", "n_variant_pipeline"), + value = c(source_desc, + format(Sys.time(), tz = "UTC", usetz = TRUE), + "write_mutations_db", + counts[["maf"]], counts[["seg"]], + counts[["bedpe"]], nrow(sample_data$meta), n_sample_study, n_variant_pipeline), + stringsAsFactors = FALSE + ), overwrite = TRUE) + + DBI::dbExecute(con, "VACUUM") + DBI::dbExecute(con, "ANALYZE") + + message(sprintf("Wrote %s (%.0f MB) maf=%d seg=%d bedpe=%d samples=%d", + out_db, file.info(out_db)$size / 1024^2, + counts[["maf"]], counts[["seg"]], + counts[["bedpe"]], nrow(sample_data$meta))) + invisible(out_db) +} diff --git a/data/gambl_metadata.rda b/data/gambl_metadata.rda deleted file mode 100644 index d830dac..0000000 Binary files a/data/gambl_metadata.rda and /dev/null differ diff --git a/envs/r.yaml b/envs/r.yaml index c6af660..d256efa 100644 --- a/envs/r.yaml +++ b/envs/r.yaml @@ -149,6 +149,7 @@ dependencies: - r-rmarkdown=2.25 - r-roxygen2=7.2.3 - r-rprojroot=2.0.3 + - r-rsqlite=2.3.7=r43h0d4f4ea_0 - r-rstudioapi=0.15.0 - r-rversions=2.1.2 - r-sass=0.4.7 diff --git a/inst/extdata/gambl_reference.db b/inst/extdata/gambl_reference.db new file mode 100644 index 0000000..8f44e11 Binary files /dev/null and b/inst/extdata/gambl_reference.db differ diff --git a/man/download_gambl_mutations_db.Rd b/man/download_gambl_mutations_db.Rd new file mode 100644 index 0000000..4eb8796 --- /dev/null +++ b/man/download_gambl_mutations_db.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/download_gambl_mutations_db.R +\name{download_gambl_mutations_db} +\alias{download_gambl_mutations_db} +\title{Download the GAMBL mutations database.} +\usage{ +download_gambl_mutations_db( + dest_path = file.path(tools::R_user_dir("GAMBLR.data", "cache"), "gambl_mutations.db"), + tag = paste0("data-v", utils::packageVersion("GAMBLR.data")), + repo = "morinlab/GAMBLR.data", + overwrite = FALSE +) +} +\arguments{ +\item{dest_path}{Where to save the downloaded file. Defaults to the +standard cache location, \code{tools::R_user_dir("GAMBLR.data", "cache")}.} + +\item{tag}{GitHub release tag to download from. Defaults to +\code{paste0("data-v", utils::packageVersion("GAMBLR.data"))}.} + +\item{repo}{GitHub repo in \code{"owner/repo"} form.} + +\item{overwrite}{Re-download even if \code{dest_path} already exists. Default +\code{FALSE}.} +} +\value{ +\code{dest_path}, invisibly. +} +\description{ +Downloads \code{gambl_mutations.db} from this package's GitHub +Release assets, if not already present. Normally called automatically by +\code{\link[=gambl_mutations_db]{gambl_mutations_db()}} the first time it's needed; call it directly to +pre-fetch the file (e.g. before a batch job on a machine with no other +network access) or to force a re-download. +} +\details{ +Requires the \code{piggyback} package (Suggests, not a hard +dependency of GAMBLR.data -- most calls to \code{\link[=gambl_mutations_db]{gambl_mutations_db()}} never +need to download anything). Downloads are tagged per package version by +default (\code{"data-v"}), so a given install always fetches data +built for its own schema, not whatever happens to be the newest release +-- see \code{data-raw/release_mutations_db.R} for the corresponding upload +step (admin-only, run manually after a fresh build; bump +\code{DESCRIPTION}'s \code{Version} and cut a matching release whenever the schema +changes, so the two stay in lockstep). +} diff --git a/man/gambl_metadata.Rd b/man/gambl_metadata.Rd deleted file mode 100644 index 03d44c1..0000000 --- a/man/gambl_metadata.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{gambl_metadata} -\alias{gambl_metadata} -\title{GAMBL Metadata} -\format{ -\subsection{\code{gambl_metadata}}{ - -A data frame with 4785 rows and 27 columns -\describe{ -\item{sample_id}{Sample identifier.} -\item{patient_id}{Patient identifier.} -\item{pathology}{Pathology.} -\item{seq_type}{Sample sequencing type.} -\item{genome_build}{Genome build the sample coordinates are in reference to.} -\item{pairing_status}{Matched or unmatched.} -\item{Tumor_Sample_Barcode}{Sample ID in another column, needed for certain functions.} -\item{age_group}{Sample age group.} -\item{compression}{The compression available for a particular sample.} -\item{bam_available}{Boolean.} -\item{pathology_rank}{Pathology rank.} -\item{cohort}{Sample cohort} -\item{COO_consensus}{COO consensus.} -\item{DHITsig_consensus}{DHIT signature consensus.} -\item{EBV_status_inf}{EBV status.} -\item{ffpe_or_frozen}{FFPE or frozen.} -\item{fl_grade}{FL grade.} -\item{hiv_status}{Sample HIV status.} -\item{lymphgen}{Lymphgen.} -\item{lymphgen_cnv_noA53}{Lymphgen with CNV no A53.} -\item{lymphgen_no_cnv}{Lymphgen no CNV.} -\item{lymphgen_with_cnv}{Lymphgen with CNV} -\item{lymphgen_wright}{Lymphgen Wright.} -\item{molecular_BL}{Molecualr BL.} -\item{normal_sample_id}{Normal sample ID} -\item{sex}{Female or Male} -\item{time_point}{Smaple timepoint.} -} -} -} -\usage{ -gambl_metadata -} -\description{ -A data frame with metadata for a collection of GAMBL samples. -} -\keyword{internal} diff --git a/man/gambl_mutations_db.Rd b/man/gambl_mutations_db.Rd new file mode 100644 index 0000000..4bf26e0 --- /dev/null +++ b/man/gambl_mutations_db.Rd @@ -0,0 +1,120 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gambl_mutations_db.R +\name{gambl_mutations_db} +\alias{gambl_mutations_db} +\title{Connect to the GAMBL mutations database.} +\usage{ +gambl_mutations_db(db_path = NULL, auto_download = TRUE) +} +\arguments{ +\item{db_path}{Optional explicit path to the .db file.} + +\item{auto_download}{When no existing copy is found, download it +automatically via \code{\link[=download_gambl_mutations_db]{download_gambl_mutations_db()}}. Default \code{TRUE}; +can also be disabled globally with +\code{options(GAMBLR.data.auto_download = FALSE)}.} +} +\value{ +A read-only DBIConnection. +} +\description{ +Returns a (cached) read-only DBI connection to \code{gambl_mutations.db}, +the SQLite database holding the large sample-level frames (MAF, seg, +bedpe) built from the in-memory \code{sample_data} object assembled by +\code{data-raw/assemble_bundled_data.R}, via the shared \code{write_mutations_db()} +helper (\code{data-raw/write_mutations_db.R}). +} +\details{ +Unlike the small bundled \code{\link[=gambl_reference_db]{gambl_reference_db()}}, this file is large +(several hundred MB) and is NOT shipped inside the package. It is +distributed as a GitHub Release asset (see \code{\link[=download_gambl_mutations_db]{download_gambl_mutations_db()}}) +and cached locally. The path is resolved in this order: +\enumerate{ +\item the \code{db_path} argument; +\item \code{getOption("GAMBLR.data.mutations_db")}; +\item the \code{GAMBLR_MUTATIONS_DB} environment variable; +\item the user cache dir, \code{tools::R_user_dir("GAMBLR.data", "cache")}; +\item \code{gambl_mutations.db} in the working directory (development). +} +If none of these exist and \code{auto_download} is \code{TRUE} (the default), the +file is downloaded automatically into the user cache dir via +\code{\link[=download_gambl_mutations_db]{download_gambl_mutations_db()}}. Set \code{auto_download = FALSE}, or +\code{options(GAMBLR.data.auto_download = FALSE)} to disable this globally +(e.g. offline/CI environments) and get an explicit error instead. +} +\section{Table reference}{ +Both genome builds are stacked into the same table (not split across +separate tables); filter on \code{genome_build} (\code{"grch37"} / \code{"hg38"}) to select +one. \code{maf} holds both coding-panel and aSHM-region mutations, pooled +unconditionally at write time (a mutation found by both pulls -- common +for genes that are both aSHM targets and coding-panel genes, e.g. BCL2, +BCL6, MYC, PIM1 -- collapses to one row, not two; see +\code{\link[=get_ssm_from_db]{get_ssm_from_db()}}).\tabular{llll}{ + Table \tab Rows (typical) \tab Grain \tab Key columns \cr + \code{maf} \tab ~3.7M \tab one somatic mutation call \tab \code{mutation_id} (surrogate integer key, unique within \code{maf} across both genome builds -- assigned at write time, after deduplication; \code{variant_pipeline} references it as a foreign key), \verb{Hugo_Symbol, Chromosome, Start_Position, End_Position, Tumor_Sample_Barcode, Variant_Classification, HGVSp_Short, t_alt_count, n_alt_count, Pipeline (lowercase, e.g. }"slms-3"\code{/}"publication"\verb{ -- normalized at write time so it can be matched with a plain, indexed equality), genome_build, maf_seq_type} (\code{"genome"}/\code{"capture"}, so a Tumor_Sample_Barcode with both isn't indistinguishable once merged; \code{NA} for publication-pipeline rows, which don't carry per-row seq_type provenance -- see \code{\link[=get_ssm_from_db]{get_ssm_from_db()}}'s \code{seq_type} parameter) (+ ~35 more MAF-standard columns) \cr + \code{seg} \tab ~126k \tab one copy-number segment \tab \verb{ID, chrom, start, end, LOH_flag, log.ratio, CN, genome_build} \cr + \code{bedpe} \tab ~900 \tab one Manta structural-variant breakpoint pair \tab \verb{CHROM_A, START_A, END_A, CHROM_B, START_B, END_B, manta_name, SCORE, STRAND_A, STRAND_B, tumour_sample_id, normal_sample_id, VAF_tumour, DP, pair_status, FILTER, genome_build} \cr + \code{sample_meta} \tab ~3.3k \tab one sample \tab \verb{patient_id, sample_id, Tumor_Sample_Barcode, seq_type, pathology, cohort, study, ...} (its own \code{genome_build} column records the sample's native alignment build, unrelated to the per-row build stamp used in the other tables) \cr + \code{sample_study} \tab varies \tab one (sample, study) membership fact -- many-to-many, a sample belonging to N studies is N rows \tab \verb{sample_id, study, study_id, reference_PMID}. \code{study_id} is that study's own identifier for the sample where it differs from GAMBL's \code{sample_id} (e.g. a paper's own case/patient ID); NA where no study-specific ID has been sourced. \cr + \code{variant_pipeline} \tab varies \tab one (mutation, Pipeline) fact -- many-to-many, a variant independently called by N pipelines is N rows \tab \verb{mutation_id, Pipeline}. \code{mutation_id} is a foreign key to \code{maf.mutation_id} -- not a natural key, so this table doesn't repeat the variant's own position/sample columns. \cr + \code{build_info} \tab 9 \tab key/value \tab provenance (\code{source}, \code{built_at}, \code{builder}) and expected row counts (\code{n_maf}, \code{n_seg}, \code{n_bedpe}, \code{n_samples}, \code{n_sample_study}, \code{n_variant_pipeline}) used by the \code{data-raw/test_gambl_db.R} regression checks \cr +} + + +\code{maf} does NOT carry a \code{Study} column. Cohort/study membership is +tracked once per sample in \code{sample_study}, not once per mutation row -- +a mutation row's \code{Study} couldn't represent a sample belonging to more +than one study anyway, and a per-row tag caused the same sample's +mutations to be pulled and duplicated once per cohort-specific code path +that claimed it (see \code{GAMBLR.data::get_ssm_from_db()}'s \code{this_study} +parameter for how to filter by study post-refactor). + +Similarly, \code{maf} only ever carries ONE \code{Pipeline} value per variant, +even though the same real mutation can be independently produced by more +than one pipeline (e.g. a cohort's own published/curated maf and our +separate SLMS-3 recall both calling the same position) -- when that +happens, \code{write_mutations_db()}'s write-time dedup keeps one row +arbitrarily (whichever pipeline's data was assembled first), which would +otherwise make the variant invisible to a query for the \emph{other} pipeline +even though it's still fully present in the table. \code{variant_pipeline} +records every pipeline that actually produced each variant, independent +of which one "won" as \code{maf}'s single representative row; see +\code{GAMBLR.data::get_ssm_from_db()}'s \code{tool_name} parameter, which resolves +against this table rather than \code{maf}'s own \code{Pipeline} column. +\subsection{Join keys}{ + +\code{maf} and \code{bedpe} (via \code{tumour_sample_id}) key to +\code{sample_meta$Tumor_Sample_Barcode} / \code{sample_meta$sample_id}; \code{seg$ID} keys +to \code{sample_meta$sample_id}; \code{sample_study$sample_id} keys to +\code{sample_meta$sample_id}; \code{variant_pipeline$mutation_id} keys to +\code{maf.mutation_id}. +} + +\subsection{Indexes}{ + +\code{maf}: \verb{(genome_build, Chromosome, Start_Position)}, +\code{Tumor_Sample_Barcode}, \code{Pipeline}, \code{maf_seq_type}, \code{mutation_id}. No index on +\code{Hugo_Symbol}: gene-restricted queries resolve the gene to a region first +(the same logic used to populate this table) and filter on +\verb{(genome_build, Chromosome, Start_Position)} instead -- see +\code{GAMBLR.data::get_ssm_from_db()}. \code{seg}: \verb{(genome_build, ID)}, +\verb{(genome_build, chrom, start)}. \code{bedpe}: \verb{(tumour_sample_id, genome_build)}, \verb{(genome_build, CHROM_A, START_A)}, \verb{(genome_build, CHROM_B, START_B)} (one per breakpoint end, so a region search that OR's +both ends can use a different index per side). +\code{sample_meta}: \code{sample_id}, \code{Tumor_Sample_Barcode}. \code{sample_study}: +\code{sample_id}, \code{study}. \code{variant_pipeline}: \code{mutation_id}, \code{Pipeline}. +} +} + +\examples{ +\dontrun{ +con <- gambl_mutations_db() +DBI::dbListTables(con) + +# coding mutations in TP53, grch37, one pipeline (Pipeline is stored +# lowercase -- see the Pipeline column note above) +dplyr::tbl(con, "maf") |> + dplyr::filter(genome_build == "grch37", Hugo_Symbol == "TP53", + Pipeline == "slms-3") |> + dplyr::collect() +} +} diff --git a/man/gambl_reference_db.Rd b/man/gambl_reference_db.Rd new file mode 100644 index 0000000..03efd5d --- /dev/null +++ b/man/gambl_reference_db.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_genes_db.R +\name{gambl_reference_db} +\alias{gambl_reference_db} +\title{Connect to the bundled GAMBLR reference database.} +\usage{ +gambl_reference_db(db_path = NULL) +} +\arguments{ +\item{db_path}{Optional explicit path to the .db file. Defaults to the copy +installed with the package; falls back to \code{inst/extdata/gambl_reference.db} +during development.} +} +\value{ +A DBIConnection. +} +\description{ +Returns a (cached) read-only DBI connection to the normalized +SQLite reference database built from the LLMPP curated gene lists. This is the +shared backing store used by \code{\link[=get_genes]{get_genes()}} (and intended for the Lymphopedia +wiki / data API), replacing the versioned \verb{lymphoma_genes_*} .rda objects. +} +\details{ +The database is rebuilt from source by \code{data-raw/build_reference_db.R} +and is small (well under 1 MB), so — unlike \code{\link[=gambl_mutations_db]{gambl_mutations_db()}} — it is +bundled inside the package (\code{inst/extdata/gambl_reference.db}). +} +\section{Table reference}{ +\tabular{lll}{ + Table \tab Grain \tab Key columns \cr + \code{gene_entity} \tab one row per gene x lymphoma entity \tab \verb{gene, ensembl_gene_id, entity, tier, ashm, qc, citekey, pmid, mutation_effect, extra} (JSON of entity-specific fields), \code{source_version} \cr + \code{ashm_regions} \tab one row per gene x aSHM region x genome build \tab \verb{gene, chrom, start, end, region, regulatory_comment, build, source_version} \cr + \code{build_info} \tab key/value \tab provenance: \code{source} (LLMPP resources/curated), \code{source_version}, \code{built_at}, \code{builder} \cr +} + + +\code{gene_entity$entity} takes values like \code{"DLBCL"}, \code{"BL"}, \code{"FL"}, \code{"PMBL"}, +\code{"MCL"}, \code{"MZL"}; \code{tier} is \code{1} (high-confidence), \code{2} (lower-confidence), or +\code{3} (retired). \code{ashm_regions$build} is \code{"grch37"} or \code{"hg38"}. +} + +\examples{ +\dontrun{ +con <- gambl_reference_db() +dplyr::tbl(con, "gene_entity") + +# high-confidence DLBCL genes +dplyr::tbl(con, "gene_entity") |> + dplyr::filter(entity == "DLBCL", tier == 1) |> + dplyr::collect() +} +} diff --git a/man/get_cn_segments_from_db.Rd b/man/get_cn_segments_from_db.Rd new file mode 100644 index 0000000..9118654 --- /dev/null +++ b/man/get_cn_segments_from_db.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_cn_segments_from_db.R +\name{get_cn_segments_from_db} +\alias{get_cn_segments_from_db} +\title{Query copy-number segments from the mutations database.} +\usage{ +get_cn_segments_from_db(projection = "grch37", sample_ids = NULL, con = NULL) +} +\arguments{ +\item{projection}{Genome build / \code{genome_build} column value. Default "grch37".} + +\item{sample_ids}{Optional character vector of segment \code{ID}s (sample ids) to keep.} + +\item{con}{Optional DBI connection (defaults to \code{\link[=gambl_mutations_db]{gambl_mutations_db()}}).} +} +\value{ +A data frame of seg rows (the \code{genome_build} helper column is dropped; +chromosome-prefix handling is left to the caller, matching the bundled seg). +} +\description{ +Data-access helper that retrieves copy-number segments from the +\code{seg} table of \code{gambl_mutations.db}. Backs the GAMBLR.open \code{get_cn_segments} +and \code{get_sample_cn_segments} accessors. +} diff --git a/man/get_genes.Rd b/man/get_genes.Rd index 247de33..7be853a 100644 --- a/man/get_genes.Rd +++ b/man/get_genes.Rd @@ -7,6 +7,7 @@ get_genes( entities = c("DLBCL", "MCL", "BL"), curated_only = TRUE, + tier = NULL, gene_format = "symbol", version = "_latest" ) @@ -15,29 +16,45 @@ get_genes( \item{entities}{Optional vector specifying one or more lymphoma entities e.g. MCL, DLBCL, BL.} -\item{curated_only}{Specify FALSE to retrieve all genes or leave default -for the curated subset.} +\item{curated_only}{When TRUE (default) return the curated set (Tier 1 and +Tier 2) and drop retired Tier 3 genes. Set FALSE to include all tiers. +Ignored when \code{tier} is supplied.} + +\item{tier}{Optional integer vector of confidence tiers to keep, e.g. +\code{tier = 1} for high-confidence genes only, or \code{tier = c(1, 2)}. When +supplied it takes precedence over \code{curated_only}.} \item{gene_format}{Specify what to return as output. Can be one of: * "symbol" (the default): list of gene symbols -* "ensembl": list of ENSEMBLE IDs -* "data.frame": data frame with column Gene and per-entity gene status} +* "ensembl": list of Ensembl IDs +* "data.frame": one row per gene with a \verb{_Tier} column per +requested entity (the gene's tier there, or NA if absent). Filter +these columns to focus on tiers, e.g. +\code{dplyr::filter(df, DLBCL_Tier == 1, FL_Tier == 1) |> dplyr::pull(Gene)}.} -\item{version}{Specify which version to return. Currently supported versions -are 0.0 (legacy version from original GAMBLR), 0.1, and _latest. The -latter will always point to the highest numeric version of the genes.} +\item{version}{Data version to return. Defaults to "_latest", the current +LLMPP-sourced version in the reference database. Legacy numeric versions +(0.0, 0.1) are no longer bundled.} } \value{ -A character vector of gene symbol or Ensembl IDs or a data frame. +A character vector of gene symbols or Ensembl IDs, or a data frame. } \description{ Retrieve gene names from bundled lymphoma gene lists. } \details{ -Complete lists of genes described as significantly mutated in large -lymphoma studies have been curated and provided with this package. +Gene lists are read from the bundled normalized reference database +(\code{inst/extdata/gambl_reference.db}), which is built from the LLMPP curated +lists. Confidence tier and data version are columns in that database, so the +previous per-entity, per-version \code{.rda} objects are no longer required. } \examples{ -all_dlbcl_genes <- get_genes(entities = "DLBCL", curated_only = FALSE) +# high-confidence (Tier 1) DLBCL genes +dlbcl_tier1 <- get_genes(entities = "DLBCL", tier = 1) + +# genes that are Tier 1 in DLBCL, FL and BL simultaneously +gene_df <- get_genes(entities = c("DLBCL", "FL", "BL"), gene_format = "data.frame") +shared_tier1 <- dplyr::filter(gene_df, DLBCL_Tier == 1, FL_Tier == 1, BL_Tier == 1) + all_curated_genes <- get_genes() } diff --git a/man/get_ssm_from_db.Rd b/man/get_ssm_from_db.Rd new file mode 100644 index 0000000..04db7c4 --- /dev/null +++ b/man/get_ssm_from_db.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_ssm_from_db.R +\name{get_ssm_from_db} +\alias{get_ssm_from_db} +\title{Query SSMs from the mutations database.} +\usage{ +get_ssm_from_db( + projection = "grch37", + sample_ids = NULL, + tool_name = "slms-3", + coding_only = FALSE, + include_silent = TRUE, + min_read_support = 0, + this_study = NULL, + regions = NULL, + seq_type = NULL, + con = NULL +) +} +\arguments{ +\item{projection}{Genome build / \code{genome_build} column value. Default "grch37".} + +\item{sample_ids}{Optional character vector of \code{Tumor_Sample_Barcode} to keep.} + +\item{tool_name}{Pipeline to keep, matched case-insensitively (the value +you pass is lowercased before comparing; the stored Pipeline column is +itself normalized to lowercase at write time, so this is a plain, +indexed equality check rather than a function-wrapped one). Default +"slms-3"; set NULL to skip the Pipeline filter. Resolved via the +\code{variant_pipeline} join table (see \code{\link[=gambl_mutations_db]{gambl_mutations_db()}}'s schema +docs) rather than \code{maf}'s own Pipeline column, so a variant +independently called by more than one pipeline is still found by a +query for either one, even though its single maf row can only +carry one Pipeline value. Falls back to filtering the column directly +against DBs built before variant_pipeline existed.} + +\item{coding_only}{When TRUE, keep only coding \code{Variant_Classification}s.} + +\item{include_silent}{When FALSE (and \code{coding_only}), drop Silent mutations.} + +\item{min_read_support}{Keep only variants with \code{t_alt_count} >= this value.} + +\item{this_study}{Optional single study to restrict to, matched against +\code{sample_study$study} (not a \code{maf} column -- see +\code{\link[=gambl_mutations_db]{gambl_mutations_db()}}'s schema docs). Resolves to the set of sample_ids +belonging to that study, intersected with \code{sample_ids} if both are +supplied; composes with \code{tool_name}/Pipeline exactly as before (e.g. +\code{this_study="Reddy"} alone returns Reddy's SLMS-3 recall by default, +while \verb{this_study="Reddy", tool_name="publication"} returns Reddy's +as-published rows only).} + +\item{regions}{Optional data frame with columns \code{chrom}, \code{start}, \code{end}; +rows are OR-ed and applied as a single combined query (not one query per +region), so passing hundreds of regions (e.g. one per gene in a panel) +doesn't cost hundreds of separate round trips.} + +\item{seq_type}{Optional character vector to restrict rows to one or more +\code{maf_seq_type} values (\code{"genome"}/\code{"capture"}), so a sample_id with both +can be queried separately instead of getting both back merged together. +Default \code{NULL} (no filter, matching prior behaviour). Rows predating +\code{maf_seq_type} (e.g. publication-pipeline rows, which don't carry +per-row seq_type provenance) have \code{NA} here and are excluded whenever +this filter is set, same as any other \code{NA} under a \code{\%in\%} filter. Silently +ignored against a \code{gambl_mutations.db} built before this column existed.} + +\item{con}{Optional DBI connection (defaults to \code{\link[=gambl_mutations_db]{gambl_mutations_db()}}).} +} +\value{ +A data frame of MAF rows (the \code{genome_build}/\code{mutation_id} helper +columns are dropped). +} +\description{ +Shared data-access helper that retrieves simple somatic mutations +(SSM, MAF format) from \code{gambl_mutations.db}, pushing the common filters down to +indexed SQL. This is the single query path behind the GAMBLR.open SSM +accessors (get_coding_ssm, get_ssm_by_samples, get_ssm_by_region(s), ...). +} diff --git a/man/get_sv_from_db.Rd b/man/get_sv_from_db.Rd new file mode 100644 index 0000000..bd8f74c --- /dev/null +++ b/man/get_sv_from_db.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_sv_from_db.R +\name{get_sv_from_db} +\alias{get_sv_from_db} +\title{Query structural variants (SV) from the mutations database.} +\usage{ +get_sv_from_db( + projection = "grch37", + sample_ids = NULL, + min_vaf = 0, + min_score = 0, + pass_only = FALSE, + region = NULL, + con = NULL +) +} +\arguments{ +\item{projection}{Genome build / \code{genome_build} column value. Default "grch37".} + +\item{sample_ids}{Optional character vector of \code{tumour_sample_id} to keep.} + +\item{min_vaf}{Minimum \code{VAF_tumour}. Default 0.} + +\item{min_score}{Minimum Manta \code{SCORE}. Default 0.} + +\item{pass_only}{When TRUE, keep only \code{FILTER == "PASS"} breakpoints.} + +\item{region}{Optional list with \code{chrom}, \code{start}, \code{end}; keeps breakpoints +whose A- or B-end falls in the range (positions compared numerically).} + +\item{con}{Optional DBI connection (defaults to \code{\link[=gambl_mutations_db]{gambl_mutations_db()}}).} +} +\value{ +A data frame of BEDPE rows (the \code{genome_build} helper column is dropped). +} +\description{ +Data-access helper that retrieves Manta SV breakpoints (BEDPE +format) from the \code{bedpe} table of \code{gambl_mutations.db}, pushing the common +filters down to SQL. Backs the GAMBLR.open \code{get_manta_sv} accessor. +}