diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index 2646a3baa..8315e4130 100644 --- a/code/SoS/reference_data/rss_ld_sketch.ipynb +++ b/code/SoS/reference_data/rss_ld_sketch.ipynb @@ -19,27 +19,7 @@ "kernel": "SoS" }, "source": [ - "## Overview\n", - "\n", - "This pipeline generates a stochastic genotype sample **U = WᵀG** from whole-genome sequencing VCF files and stores it as a PLINK2 pgen file for use as an LD reference panel with SuSiE-RSS fine-mapping.\n", - "\n", - "Rather than storing the full genotype matrix G (n × p), the sketch computes U = WᵀG (B × p) using a random projection matrix $W \\sim N(0, 1/\\sqrt{n})$. The approximate LD matrix $R = U^T U / B \\approx G^T G / n$ by the Johnson–Lindenstrauss lemma. G is never stored.\n", - "\n", - "**Matrix dimensions:**\n", - "- G : (n × p) — n individuals × p variants\n", - "- W : (n × B) — projection matrix, generated once per cohort\n", - "- U : (B × p) — stochastic genotype sample = WᵀG, stored in pgen\n", - "- $\\hat{R}$ : (p × p) — approximate LD matrix, computed on-the-fly by SuSiE-RSS from U\n", - "\n", - "The workflow has three steps run in order: `generate_W` (build the projection matrix), `process_block` (read VCF per LD block and write per-block dosage sketches), and `merge_chrom` (merge per-block dosages into one per-chromosome pgen).\n", - "\n", - "Summary-statistic fine-mapping needs an LD matrix, and storing the full genotype matrix for\n", - "a whole-genome panel is expensive. The sketch keeps a random projection of the genotypes\n", - "instead: small enough to distribute, but sufficient to reconstruct the LD structure that\n", - "SuSiE-RSS actually uses.\n", - "\n", - "**When to run it.** Once per reference panel, before any `rss_analysis` run that needs an\n", - "LD reference. The output is reference data, not a per-study result." + "## Overview\n\nThis pipeline generates a stochastic genotype sample **U = WᵀG** from whole-genome sequencing VCF files and stores it as a PLINK2 pgen file for use as an LD reference panel with SuSiE-RSS fine-mapping.\n\nRather than storing the full genotype matrix G (n × p), the sketch computes U = WᵀG (B × p) using a random projection matrix $W \\sim N(0, 1/\\sqrt{n})$. The approximate LD matrix $R = U^T U / B \\approx G^T G / n$ by the Johnson–Lindenstrauss lemma. G is never stored.\n\n**Matrix dimensions:**\n- G : (n × p) — n individuals × p variants\n- W : (n × B) — projection matrix, generated once per cohort\n- U : (B × p) — stochastic genotype sample = WᵀG, stored in pgen\n- $\\hat{R}$ : (p × p) — approximate LD matrix, computed on-the-fly by SuSiE-RSS from U\n\nThe workflow has three steps run in order: `generate_W` (build the projection matrix), `process_block` (read VCF per LD block and write per-block dosage sketches), and `merge_chrom` (merge per-block dosages into one per-chromosome pgen). During `merge_chrom`, allele-frequency rows are reconciled to the final sorted `.pvar`; missing or duplicate IDs stop the workflow, while extra block rows are reported and excluded. During `merge_chrom`, the step also detects indels that form an exact REF/ALT-swap (“mirror”) pair at the same position — the ambiguous cases where insertion vs. deletion anchoring can flip effect-allele orientation during association harmonization. Only these mirror-pair indels are assigned a canonical directional event ID (e.g. `chr22:pos:INS:T` and `chr22:pos:DEL:T`); the event ID replaces the variant ID in `.pvar` and `.afreq`, and the original↔event mapping is recorded in a `.event_id.tsv`. All other variants — SNPs, equal-length substitutions, and non-mirror indels — keep their standard IDs. If a panel contains no mirror pairs, no `.event_id.tsv` is written and `.pvar`/`.afreq` are left fully standard. The `.pgen` genotype data is never modified.\n\nSummary-statistic fine-mapping needs an LD matrix, and storing the full genotype matrix for\na whole-genome panel is expensive. The sketch keeps a random projection of the genotypes\ninstead: small enough to distribute, but sufficient to reconstruct the LD structure that\nSuSiE-RSS actually uses.\n\n**When to run it.** Once per reference panel, before any `rss_analysis` run that needs an\nLD reference. The output is reference data, not a per-study result." ] }, { @@ -61,7 +41,7 @@ "```\n", "\n", "- `--vcf-base` **`input/rss_ld_sketch`** together with `--vcf-prefix protocol_example.genotype.`\n", - "(directory of bgzipped, tabix-indexed VCFs named `{vcf_prefix}{chr}.*.bgz`; the toy cohort is `protocol_example.genotype.chr22.bgz`, 60 individuals)\n", + "(directory of bgzipped, tabix-indexed VCFs named `{vcf_prefix}{chr}.*.bgz` or `{vcf_prefix}{chr}.*.vcf.gz`; the toy cohort is `protocol_example.genotype.chr22.vcf.gz`, 60 individuals)\n", "\n", "```\n", "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMPLE_001\t...\n", @@ -73,7 +53,7 @@ "- `--cwd output/rss_ld_sketch` (working directory for logs and job files. Defaults to `output`.)\n", "- `--output-dir output/rss_ld_sketch` (directory the sketch products are written to)\n", "- `--seed 999` (random seed for the projection matrix, so a sketch can be reproduced)\n", - "- `--W-matrix` **`output/rss_ld_sketch/W_B50.npy`**\n", + "- `--W-matrix` **`output/rss_ld_sketch/W_B50.rds`**\n", "(the shared projection matrix, passed to `process_block` after `generate_W` has written it)\n", "- `--chrom chr22` (chromosome to process)\n", "- `--cohort-id protocol_example` (tag used in the per-block and merged output filenames)" @@ -88,13 +68,11 @@ "source": [ "## Output\n", "\n", - "- **`output/rss_ld_sketch/W_B50.npy`**\n", + "- **`output/rss_ld_sketch/W_B50.rds`**\n", "(shared projection matrix written once by `generate_W`, shape n_samples x B; `process_block` reads it back for every LD block so all blocks share one projection)\n", "\n", "```\n", - "float32 (60, 50)\n", - "[[-0.1402 0.1288 0.0365 -0.1945 -0.0747]\n", - " [-0.1671 -0.1341 0.2251 -0.1030 0.0038]]\n", + "R matrix, double (60, 50) # saveRDS/readRDS; one row per sample, one column per sketch dim\n", "```\n", "\n", "- **`output/rss_ld_sketch/chr22/chr22__/.chr22__.dosage.gz`** (`process_block`)\n", @@ -127,9 +105,15 @@ "chr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n", "```\n", "\n", - "These feed SuSiE-RSS fine-mapping: load with a metadata TSV (one row per chromosome, columns `#chrom start end path`, `path` = pgen prefix). Use the X (genotype) interface for `susie_rss(z, X=X)` or the R (correlation) interface for `susie_rss(z, R=R)`.\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv`** (mirror-pair indel mapping, written by `merge_chrom` only when such pairs exist; records the original ID alongside the canonical event ID)\n", + "```\n", + "ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\n", + "chr22:16500000:A:AT 22 16500000 A AT chr22:16500000:INS:T INS 16500000 T\n", + "chr22:16500000:AT:A 22 16500000 AT A chr22:16500001:DEL:T DEL 16500001 T\n", + "```\n", + "Only indels that form an exact REF/ALT-swap pair at the same position are included — the cases where insertion vs. deletion anchoring is ambiguous. The example above shows one such pair: the same single-base event at chr22:16500000 represented as both an insertion (`A`→`AT`) and a deletion (`AT`→`A`), each assigned a canonical event ID. SNPs, equal-length substitutions, symbolic alleles, and non-mirror indels are omitted and retain their standard IDs.\n", "\n", - "The `generate_W` step declares its output as `W_B.rds`, but the example on disk is `W_B50.npy`. The code is authoritative; the example data predates the current naming." + "These feed SuSiE-RSS fine-mapping: load with a metadata TSV (one row per chromosome, columns `#chrom start end path`, `path` = pgen prefix). Use the X (genotype) interface for `susie_rss(z, X=X)` or the R (correlation) interface for `susie_rss(z, R=R)`." ] }, { @@ -189,7 +173,7 @@ " import pkg_resources\n", "INFO: Running \u001b[32mgenerate_W\u001b[0m: \n", "INFO: \u001b[32mgenerate_W\u001b[0m is \u001b[32mcompleted\u001b[0m.\n", - "INFO: \u001b[32mgenerate_W\u001b[0m output: \u001b[32moutput/rss_ld_sketch/W_B50.npy\u001b[0m\n", + "INFO: \u001b[32mgenerate_W\u001b[0m output: \u001b[32moutput/rss_ld_sketch/W_B50.rds\u001b[0m\n", "INFO: Workflow generate_W (ID=w68f63c60d8da4b5e) is executed successfully with 1 completed step.\n" ] } @@ -267,7 +251,7 @@ "kernel": "SoS" }, "source": [ - "### Step 3. Merge the per-block dosage sketches into one per-chromosome PLINK2 pgen." + "### Step 3. Merge the per-block dosage sketches, reconcile allele frequencies, and relabel mirror-pair indels with canonical event IDs.\n" ] }, { @@ -572,7 +556,7 @@ }, "outputs": [], "source": [ - "[merge_chrom]\nparameter: chrom = 0\nparameter: output_dir = str\nparameter: cohort_id = str\nparameter: plink2_bin = \"plink2\"\n\nimport os, glob\n\nif chrom != 0:\n chroms = [f\"chr{chrom}\"]\nelse:\n chroms = sorted(set(\n os.path.basename(_d)\n for _d in glob.glob(os.path.join(output_dir, \"chr*\"))\n if os.path.isdir(_d)\n ))\n\ninput: for_each = \"chroms\"\noutput: f\"{output_dir}/{_chroms}/{cohort_id}.{_chroms}.pgen\"\ntask: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads\nbash: expand = \"$[ ]\"\n\n set -euo pipefail\n shopt -s nullglob\n\n chrom_dir=\"$[output_dir]/$[_chroms]\"\n final_prefix=\"${chrom_dir}/$[cohort_id].$[_chroms]\"\n merge_list=\"${chrom_dir}/$[cohort_id].$[_chroms]_pmerge_list.txt\"\n\n # Step 1: Convert each block dosage.gz -> sorted per-block pgen\n > \"${merge_list}\"\n files=(\"${chrom_dir}\"/*/*.dosage.gz)\n if [ ${#files[@]} -eq 0 ]; then\n echo \"No dosage files found in ${chrom_dir}\" >&2\n exit 1\n fi\n for dosage_gz in \"${files[@]}\"; do\n block_dir=$(dirname \"${dosage_gz}\")\n block_tag=$(basename \"${block_dir}\")\n prefix=\"${block_dir}/$[cohort_id].${block_tag}_tmp\"\n map_file=\"${block_dir}/$[cohort_id].${block_tag}.map\"\n psam_file=\"${block_dir}/$[cohort_id].${block_tag}.psam\"\n meta_file=\"${block_dir}/$[cohort_id].${block_tag}.meta\"\n\n B=$(grep \"^B=\" \"${meta_file}\" | cut -d= -f2)\n printf '#FID\\tIID\\n' > \"${psam_file}\"\n for i in $(seq 1 ${B}); do\n printf 'S%d\\tS%d\\n' ${i} ${i} >> \"${psam_file}\"\n done\n\n $[plink2_bin] \\\n --import-dosage \"${dosage_gz}\" format=1 noheader \\\n --psam \"${psam_file}\" \\\n --map \"${map_file}\" \\\n --make-pgen \\\n --out \"${prefix}_unsorted\" \\\n --silent\n\n $[plink2_bin] \\\n --pfile \"${prefix}_unsorted\" \\\n --make-pgen \\\n --sort-vars \\\n --out \"${prefix}\" \\\n --silent\n\n rm -f \"${prefix}_unsorted.pgen\" \"${prefix}_unsorted.pvar\" \"${prefix}_unsorted.psam\"\n echo \"${prefix}\" >> \"${merge_list}\"\n done\n\n # Step 2: Merge all per-block pgens -> one per-chrom pgen\n $[plink2_bin] \\\n --pmerge-list \"${merge_list}\" pfile \\\n --make-pgen \\\n --sort-vars \\\n --out \"${final_prefix}\"\n\n # Remove PLINK merge intermediates immediately after merge\n rm -f \"${final_prefix}-merge.pgen\" \"${final_prefix}-merge.pvar\" \"${final_prefix}-merge.psam\"\n\n # Step 3: Concatenate .afreq\n first=1\n for f in \"${chrom_dir}\"/*/*.afreq; do\n if [ \"${first}\" -eq 1 ]; then\n cat \"${f}\" > \"${final_prefix}.afreq\"\n first=0\n else\n tail -n +2 \"${f}\" >> \"${final_prefix}.afreq\"\n fi\n done\n\nR: expand = \"$[ ]\"\n\n meta_files <- list.files(\"$[output_dir]/$[_chroms]\",\n pattern = \"[.]meta$\", recursive = TRUE,\n full.names = TRUE)\n if (length(meta_files) > 0) {\n fields <- c(\"n_total\", \"n_passed\", \"n_multiallelic\", \"n_monomorphic\",\n \"n_all_na\", \"n_high_msng\", \"n_low_maf\", \"n_low_mac\")\n stats <- do.call(rbind, lapply(meta_files, function(f) {\n lines <- grep(\"^n_\", readLines(f), value = TRUE)\n kv <- strsplit(lines, \"=\")\n vals <- setNames(as.integer(sapply(kv, `[`, 2)), sapply(kv, `[`, 1))\n as.data.frame(as.list(vals[fields]))\n }))\n totals <- colSums(stats, na.rm = TRUE)\n summary <- data.frame(t(totals))\n summary$pct_dropped <- round(100 * (1 - summary$n_passed / summary$n_total), 1)\n cat(\"\\n=== Filter Summary for $[_chroms] ===\\n\")\n print(data.frame(value = unlist(summary), row.names = names(summary)))\n }\n\nbash: expand = \"$[ ]\"\n\n # Step 5: Cleanup block intermediates\n chrom_dir=\"$[output_dir]/$[_chroms]\"\n final_prefix=\"${chrom_dir}/$[cohort_id].$[_chroms]\"\n\n rm -f \"${final_prefix}_pmerge_list.txt\"\n for block_dir in \"${chrom_dir}\"/*/; do\n rm -rf \"${block_dir}\"\n done" + "[merge_chrom]\nparameter: chrom = 0\nparameter: output_dir = str\nparameter: cohort_id = str\nparameter: plink2_bin = \"plink2\"\n\nimport os, glob\n\nif chrom != 0:\n chroms = [f\"chr{chrom}\"]\nelse:\n chroms = sorted(set(\n os.path.basename(_d)\n for _d in glob.glob(os.path.join(output_dir, \"chr*\"))\n if os.path.isdir(_d)\n ))\n\ninput: for_each = \"chroms\"\noutput: f\"{output_dir}/{_chroms}/{cohort_id}.{_chroms}.pgen\"\ntask: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads\nbash: expand = \"$[ ]\"\n\n set -euo pipefail\n shopt -s nullglob\n\n chrom_dir=\"$[output_dir]/$[_chroms]\"\n final_prefix=\"${chrom_dir}/$[cohort_id].$[_chroms]\"\n merge_list=\"${chrom_dir}/$[cohort_id].$[_chroms]_pmerge_list.txt\"\n\n # Step 1: Convert each block dosage.gz -> sorted per-block pgen\n > \"${merge_list}\"\n files=(\"${chrom_dir}\"/*/*.dosage.gz)\n if [ ${#files[@]} -eq 0 ]; then\n echo \"No dosage files found in ${chrom_dir}\" >&2\n exit 1\n fi\n for dosage_gz in \"${files[@]}\"; do\n block_dir=$(dirname \"${dosage_gz}\")\n block_tag=$(basename \"${block_dir}\")\n prefix=\"${block_dir}/$[cohort_id].${block_tag}_tmp\"\n map_file=\"${block_dir}/$[cohort_id].${block_tag}.map\"\n psam_file=\"${block_dir}/$[cohort_id].${block_tag}.psam\"\n meta_file=\"${block_dir}/$[cohort_id].${block_tag}.meta\"\n\n B=$(grep \"^B=\" \"${meta_file}\" | cut -d= -f2)\n printf '#FID\\tIID\\n' > \"${psam_file}\"\n for i in $(seq 1 ${B}); do\n printf 'S%d\\tS%d\\n' ${i} ${i} >> \"${psam_file}\"\n done\n\n $[plink2_bin] \\\n --import-dosage \"${dosage_gz}\" format=1 noheader \\\n --psam \"${psam_file}\" \\\n --map \"${map_file}\" \\\n --make-pgen \\\n --out \"${prefix}_unsorted\" \\\n --silent\n\n $[plink2_bin] \\\n --pfile \"${prefix}_unsorted\" \\\n --make-pgen \\\n --sort-vars \\\n --out \"${prefix}\" \\\n --silent\n\n rm -f \"${prefix}_unsorted.pgen\" \"${prefix}_unsorted.pvar\" \"${prefix}_unsorted.psam\"\n echo \"${prefix}\" >> \"${merge_list}\"\n done\n\n # Step 2: Merge all per-block pgens -> one per-chrom pgen\n $[plink2_bin] \\\n --pmerge-list \"${merge_list}\" pfile \\\n --make-pgen \\\n --sort-vars \\\n --out \"${final_prefix}\"\n\n # Step 3: reconcile .afreq, relabel mirror-pair indels with event IDs, summarize filters, and cleanup.\n Rscript $[modular_script_dir]/reference_data/rss_ld_sketch.R --step merge_chrom --chrom-dir \"$[output_dir]/$[_chroms]\" --final-prefix \"$[output_dir]/$[_chroms]/$[cohort_id].$[_chroms]\"\n" ] }, { @@ -586,7 +570,7 @@ "\n", "| Symptom | Cause | Fix |\n", "|---|---|---|\n", - "| `No VCF files for chrXX in {vcf_base}` | VCF naming or extension mismatch | Files must end in `.bgz` and be named `{vcf_prefix}{chr}.*.bgz`; check `--vcf-base` and `--vcf-prefix`. |\n", + "| `No VCF files for chrXX in {vcf_base}` | VCF naming or extension mismatch | Files must end in `.bgz` or `.vcf.gz` and be named `{vcf_prefix}{chr}.*`; check `--vcf-base` and `--vcf-prefix`. |\n", "| `W shape mismatch` | `--n-samples` or `--B` differs from the W used | Re-run `generate_W` with the same `--n-samples` and `--B`, and pass that `W_B{B}.rds` to `process_block`. |\n", "| `No passing variants in chrXX` | Filters removed everything (small toy cohort) | Widen `--maf-min` / `--mac-min` / `--msng-min`, or choose blocks with more variants. |\n", "| `No blocks found for chrom=XX` | `--chrom` does not match any BED rows | Ensure the BED `chr` column matches (e.g. `chr22`) and `--chrom` is the matching number. |\n", @@ -630,4 +614,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index c96f06fea..117b15253 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -61,7 +61,7 @@ do_process_block <- function(argv) { if (is.na(argv$vcf_base)) stop("process_block needs --vcf or --vcf-base") all_files <- list.files(argv$vcf_base, full.names = TRUE) pat <- paste0("^", argv$vcf_prefix, chrom, "[:.]") - vcf_files <- sort(all_files[grepl("\\.bgz$", all_files) & + vcf_files <- sort(all_files[grepl("\\.(bgz|vcf\\.gz)$", all_files) & grepl(pat, basename(all_files))]) if (!length(vcf_files)) stop(sprintf("No VCF files for %s in %s", chrom, argv$vcf_base)) } @@ -164,9 +164,206 @@ do_process_block <- function(argv) { cat(sprintf("Wrote %d variants -> %s\n", n_passed, basename(paste0(pfx, ".dosage.gz")))) } + +# ---- pure: mirror-pair canonical event IDs --------------------------------- +# Given variant records (parallel CHROM/POS/ID/REF/ALT vectors), compute the +# canonical directional event-ID mapping for same-position REF/ALT-swap +# ("mirror") indel pairs -- the ambiguous cases where INS vs DEL anchoring can +# flip effect-allele orientation. Returns a data.frame (ID/CHROM/POS/REF/ALT/ +# event_id/event_type/event_pos/event_ref/event_alt) with ONLY mirror-pair rows +# (0 rows if none). No file I/O; do_merge_chrom owns reading/writing/substitution. +canonicalize_mirror_event_ids <- function(chrom, pos, id, ref, alt) { + empty_map <- data.frame( + ID = character(0), CHROM = character(0), POS = integer(0), + REF = character(0), ALT = character(0), event_id = character(0), + event_type = character(0), event_pos = integer(0), + event_ref = character(0), event_alt = character(0), + check.names = FALSE, stringsAsFactors = FALSE) + if (!length(id)) return(empty_map) + canonical_event <- function(chrom, pos, id, ref, alt) { + ref <- toupper(ref); alt <- toupper(alt); pos <- as.integer(pos) + chrom <- if (grepl("^chr", id)) sub(":.*$", "", id) else paste0("chr", sub("^chr", "", chrom)) + # Only length-changing alleles receive an internal event ID. + if (nchar(ref) == nchar(alt)) + return(c(event_id = id, event_type = "UNCHANGED", event_pos = pos, + event_ref = ref, event_alt = alt)) + if (grepl("[<>*]|\\[|\\]", ref) || grepl("[<>*]|\\[|\\]", alt) || + grepl(",", alt, fixed = TRUE)) + return(c(event_id = paste(chrom, pos, ref, alt, sep = ":"), + event_type = "SYMBOLIC", event_pos = pos, + event_ref = ref, event_alt = alt)) + r <- strsplit(ref, "", fixed = TRUE)[[1]] + a <- strsplit(alt, "", fixed = TRUE)[[1]] + prefix <- 0L + while (length(r) && length(a) && r[1] == a[1]) { + r <- r[-1]; a <- a[-1]; prefix <- prefix + 1L + } + while (length(r) && length(a) && tail(r, 1) == tail(a, 1)) { + r <- head(r, -1); a <- head(a, -1) + } + rr <- paste(r, collapse = ""); aa <- paste(a, collapse = "") + event_pos <- pos + prefix + if (!nzchar(rr) && nzchar(aa)) { + event_pos <- event_pos - 1L; type <- "INS" + id <- paste(chrom, event_pos, type, aa, sep = ":") + } else if (nzchar(rr) && !nzchar(aa)) { + type <- "DEL"; id <- paste(chrom, event_pos, type, rr, sep = ":") + } else { + type <- "SUB"; id <- paste(chrom, event_pos, type, rr, aa, sep = ":") + } + c(event_id = id, event_type = type, event_pos = event_pos, + event_ref = rr, event_alt = aa) + } + event <- t(mapply(canonical_event, chrom, pos, id, ref, alt, SIMPLIFY = TRUE)) + event_map <- data.frame(ID = id, CHROM = chrom, POS = pos, REF = ref, ALT = alt, + event, check.names = FALSE, stringsAsFactors = FALSE) + event_map <- event_map[event_map$event_type %in% c("INS", "DEL"), , drop = FALSE] + # MIRROR-ONLY: keep only indels whose exact ref/alt swap exists at the same + # position (the sign-flip cases). Non-mirror indels stay standard. + nk <- paste(event_map$CHROM, event_map$POS, event_map$REF, event_map$ALT, sep = ":") + sk <- paste(event_map$CHROM, event_map$POS, event_map$ALT, event_map$REF, sep = ":") + event_map[sk %in% nk, , drop = FALSE] +} + +# ---- merge_chrom helpers --------------------------------------------------- + +# Shared reader for the tab tables plink2 writes (.pvar/.afreq). +read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", + stringsAsFactors = FALSE) + +# Gather per-block .afreq rows, validate against the merged .pvar, reorder to +# match, and atomically install the reconciled .afreq. Returns the pvar table. +reconcile_afreq <- function(final_prefix, afreq_files) { + pvar <- read_tab(paste0(final_prefix, ".pvar")) + afreq <- do.call(rbind, lapply(afreq_files, read_tab)) + if (!"ID" %in% names(pvar) || !"ID" %in% names(afreq)) + stop("Both .pvar and .afreq must contain an ID column") + if (anyNA(pvar$ID) || anyNA(afreq$ID)) stop("Missing variant ID detected") + if (anyDuplicated(pvar$ID)) stop("Duplicate IDs in final .pvar") + if (anyDuplicated(afreq$ID)) stop("Duplicate IDs across block .afreq files") + + missing_ids <- setdiff(pvar$ID, afreq$ID) + extra_ids <- setdiff(afreq$ID, pvar$ID) + message("AFREQ_RECONCILE pvar=", nrow(pvar), " block_afreq=", nrow(afreq), + " missing=", length(missing_ids), " extra=", length(extra_ids)) + if (length(missing_ids)) { + writeLines(missing_ids, paste0(final_prefix, ".afreq.missing_ids.txt")) + stop("Final .afreq is incomplete; see .afreq.missing_ids.txt") + } + if (length(extra_ids)) + writeLines(extra_ids, paste0(final_prefix, ".afreq.extra_ids.txt")) + + afreq <- afreq[match(pvar$ID, afreq$ID), , drop = FALSE] + stopifnot(identical(as.character(afreq$ID), as.character(pvar$ID))) + afreq_tmp <- paste0(final_prefix, ".afreq.tmp") + write.table(afreq, afreq_tmp, sep = "\t", quote = FALSE, + row.names = FALSE, col.names = TRUE) + if (!file.rename(afreq_tmp, paste0(final_prefix, ".afreq"))) + stop("Could not atomically install reconciled .afreq") + pvar +} + +# Given a mirror-pair event_map (nrow > 0), write the .event_id.tsv sidecar and +# substitute the event IDs into the final .pvar and .afreq (atomic, guarded). +apply_event_id_substitution <- function(final_prefix, event_map) { + duplicated_event <- duplicated(event_map$event_id) | + duplicated(event_map$event_id, fromLast = TRUE) + if (any(duplicated_event)) { + write.table(event_map[duplicated_event, ], + paste0(final_prefix, ".event_id.collisions.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE) + stop("Canonical event-ID collision detected; see .event_id.collisions.tsv") + } + event_tmp <- paste0(final_prefix, ".event_id.tsv.tmp") + write.table(event_map, event_tmp, sep = "\t", quote = FALSE, + row.names = FALSE, col.names = TRUE) + if (!file.rename(event_tmp, paste0(final_prefix, ".event_id.tsv"))) + stop("Could not atomically install event-ID mapping") + + id2event <- setNames(event_map$event_id, event_map$ID) + pv <- read_tab(paste0(final_prefix, ".pvar")) + hit <- pv$ID %in% names(id2event) + pv$ID[hit] <- id2event[pv$ID[hit]] + if (anyDuplicated(pv$ID)) stop("Event-ID substitution produced duplicate IDs in .pvar") + pv_tmp <- paste0(final_prefix, ".pvar.tmp") + write.table(pv, pv_tmp, sep = "\t", quote = FALSE, row.names = FALSE, col.names = TRUE) + if (!file.rename(pv_tmp, paste0(final_prefix, ".pvar"))) + stop("Could not atomically install event-ID .pvar") + + af <- read_tab(paste0(final_prefix, ".afreq")) + hitf <- af$ID %in% names(id2event) + af$ID[hitf] <- id2event[af$ID[hitf]] + if (anyDuplicated(af$ID)) stop("Event-ID substitution produced duplicate IDs in .afreq") + af_tmp <- paste0(final_prefix, ".afreq.tmp") + write.table(af, af_tmp, sep = "\t", quote = FALSE, row.names = FALSE, col.names = TRUE) + if (!file.rename(af_tmp, paste0(final_prefix, ".afreq"))) + stop("Could not atomically install event-ID .afreq") +} + +# Tally the per-block .meta filter counts and print the per-chromosome summary. +summarize_block_filters <- function(chrom_dir) { + meta_files <- list.files(chrom_dir, pattern = "[.]meta$", recursive = TRUE, + full.names = TRUE) + if (!length(meta_files)) return(invisible(NULL)) + fields <- c("n_total", "n_passed", "n_multiallelic", "n_monomorphic", + "n_all_na", "n_high_msng", "n_low_maf", "n_low_mac") + stats <- do.call(rbind, lapply(meta_files, function(f) { + lines <- grep("^n_", readLines(f), value = TRUE) + kv <- strsplit(lines, "=", fixed = TRUE) + vals <- setNames(as.integer(vapply(kv, function(x) x[[2]], character(1))), + vapply(kv, function(x) x[[1]], character(1))) + as.data.frame(as.list(vals[fields])) + })) + totals <- colSums(stats, na.rm = TRUE) + summary <- data.frame(t(totals)) + summary$pct_dropped <- round(100 * (1 - summary$n_passed / summary$n_total), 1) + cat("\n=== Filter Summary for ", basename(chrom_dir), " ===\n", sep = "") + print(data.frame(value = unlist(summary), row.names = names(summary))) +} + +# Remove the plink2 merge intermediates and per-block working directories. +cleanup_merge_intermediates <- function(final_prefix, block_dirs) { + unlink(c(paste0(final_prefix, "_pmerge_list.txt"), + paste0(final_prefix, "-merge.pgen"), + paste0(final_prefix, "-merge.pvar"), + paste0(final_prefix, "-merge.psam"))) + unlink(block_dirs, recursive = TRUE) +} + +# ---- merge_chrom ----------------------------------------------------------- +do_merge_chrom <- function(argv) { + chrom_dir <- argv$chrom_dir + final_prefix <- argv$final_prefix + block_dirs <- list.dirs(chrom_dir, recursive = FALSE, full.names = TRUE) + afreq_files <- unlist(lapply(block_dirs, function(d) { + list.files(d, pattern = "[.]afreq$", full.names = TRUE) + }), use.names = FALSE) + if (!length(afreq_files)) stop("No block .afreq files found under ", chrom_dir) + + # Reconcile per-block .afreq against the merged .pvar (installs the .afreq). + pvar <- reconcile_afreq(final_prefix, afreq_files) + + # Canonical event IDs for mirror-pair indels (pure transform; defined above). + chrom_col <- if ("#CHROM" %in% names(pvar)) "#CHROM" else "CHROM" + event_map <- canonicalize_mirror_event_ids(pvar[[chrom_col]], pvar$POS, + pvar$ID, pvar$REF, pvar$ALT) + + # Only when there ARE mirror pairs do we emit a mapping and relabel IDs. + # Panels with no ambiguous indels (e.g. R4) stay fully standard. + if (nrow(event_map) > 0) { + apply_event_id_substitution(final_prefix, event_map) + message("EVENT_ID_SUBSTITUTION mirror_variants=", nrow(event_map)) + } else { + message("EVENT_ID_SUBSTITUTION mirror_variants=0 (no ambiguous indels; panel left standard)") + } + + summarize_block_filters(chrom_dir) + cleanup_merge_intermediates(final_prefix, block_dirs) +} + # ---- CLI ------------------------------------------------------------------- -p <- arg_parser("RSS LD random-projection sketch (generate_w / process_block)") -p <- add_argument(p, "--step", help = "generate_w | process_block") +p <- arg_parser("RSS LD random-projection sketch (generate_w / process_block / merge_chrom)") +p <- add_argument(p, "--step", help = "generate_w | process_block | merge_chrom") p <- add_argument(p, "--n-samples", help = "generate_w: total sample size n") p <- add_argument(p, "--B", help = "sketch dimension B", default = "10000") p <- add_argument(p, "--seed", help = "generate_w: RNG seed", default = "123") @@ -184,12 +381,20 @@ p <- add_argument(p, "--mac-min", help = "process_block: min MAC", default = "5" p <- add_argument(p, "--msng-min", help = "process_block: max missingness", default = "0.05") p <- add_argument(p, "--sample-list", help = "process_block: optional sample subset file", default = "") p <- add_argument(p, "--output", help = "generate_w: output W .rds") -argv <- parse_args(p) - -if (identical(argv$step, "generate_w")) { - do_generate_w(argv) -} else if (identical(argv$step, "process_block")) { - do_process_block(argv) -} else { - stop("--step must be 'generate_w' or 'process_block'") +p <- add_argument(p, "--chrom-dir", help = "merge_chrom: chromosome output directory", default = "") +p <- add_argument(p, "--final-prefix", help = "merge_chrom: final pgen/pvar/afreq prefix", default = "") +# Run the CLI only when executed directly (Rscript), not when sourced for +# unit testing of the pure helpers above. +if (sys.nframe() == 0L) { + argv <- parse_args(p) + + if (identical(argv$step, "generate_w")) { + do_generate_w(argv) + } else if (identical(argv$step, "process_block")) { + do_process_block(argv) + } else if (identical(argv$step, "merge_chrom")) { + do_merge_chrom(argv) + } else { + stop("--step must be 'generate_w', 'process_block', or 'merge_chrom'") + } } diff --git a/tests/fixtures/rss_ld_sketch/expected/afreq_deterministic.tsv b/tests/fixtures/rss_ld_sketch/expected/afreq_deterministic.tsv index e3c0476e6..2225c8fd7 100644 --- a/tests/fixtures/rss_ld_sketch/expected/afreq_deterministic.tsv +++ b/tests/fixtures/rss_ld_sketch/expected/afreq_deterministic.tsv @@ -27,6 +27,8 @@ chr22 chr22:16400343:T:C T C 0.125000 120 chr22 chr22:16427771:T:C T C 0.475000 120 chr22 chr22:16439458:A:G A G 0.133333 120 chr22 chr22:16442170:G:A G A 0.050000 120 +chr22 chr22:16500000:A:AT A AT 0.450000 120 +chr22 chr22:16500000:AT:A AT A 0.383333 120 chr22 chr22:16549730:A:G A G 0.058333 120 chr22 chr22:16560460:ACAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAGGAAAGCCAGCCTGCCAAGCCAGCCAAGCCAGCCAGCCACTCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGCCAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCCAGCCAGCCAGCCCACACAGCCAAGC:* ACAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAGGAAAGCCAGCCTGCCAAGCCAGCCAAGCCAGCCAGCCACTCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGCCAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCCAGCCAGCCAGCCCACACAGCCAAGC * 0.449153 118 chr22 chr22:16560533:TCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGC:CCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGC TCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGC CCAAGCCAGCCAAGTCAGCGAGCCAACCAAGCCAGCCAACTCAGCCAGCCACCTAAGCCAGCCAAGCCAGCCAGCCAGCCAAGCAAGCCAGCCAGCCAAGCCAGCCAAGCCAGCCAAGCCAGGCAGCCAGC 0.050000 120 diff --git a/tests/fixtures/rss_ld_sketch/expected/event_id.tsv b/tests/fixtures/rss_ld_sketch/expected/event_id.tsv new file mode 100644 index 000000000..c205ce295 --- /dev/null +++ b/tests/fixtures/rss_ld_sketch/expected/event_id.tsv @@ -0,0 +1,3 @@ +ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt +chr22:16500000:A:AT 22 16500000 A AT chr22:16500000:INS:T INS 16500000 T +chr22:16500000:AT:A 22 16500000 AT A chr22:16500001:DEL:T DEL 16500001 T diff --git a/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz.tbi b/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz.tbi deleted file mode 100644 index bb13bb2ae..000000000 Binary files a/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz.tbi and /dev/null differ diff --git a/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz b/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz similarity index 94% rename from tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz rename to tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz index d80543e90..8522a7f72 100644 Binary files a/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz and b/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz differ diff --git a/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz.tbi b/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz.tbi new file mode 100644 index 000000000..bd3f8d1b9 Binary files /dev/null and b/tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz.tbi differ diff --git a/tests/helpers/canonicalize_mirror_check.R b/tests/helpers/canonicalize_mirror_check.R new file mode 100644 index 000000000..f88ff0746 --- /dev/null +++ b/tests/helpers/canonicalize_mirror_check.R @@ -0,0 +1,37 @@ +#!/usr/bin/env Rscript +# Unit check for canonicalize_mirror_event_ids (pure function). +# Sources the worker (CLI is guarded) and exercises the mirror-pair rule on a +# tiny in-memory fixture. Prints "UNIT_OK" on success; stop()s on any mismatch. +args <- commandArgs(trailingOnly = TRUE) +worker <- args[1] +source(worker) + +# Fixture: one mirror pair (INS A>AT + its swap AT>A) at 100; one non-mirror +# indel (C>CA, no swap) at 200; one SNP (G>T) at 300. +chrom <- c("22","22","22","22") +pos <- c(100L,100L,200L,300L) +id <- c("chr22:100:A:AT","chr22:100:AT:A","chr22:200:C:CA","chr22:300:G:T") +ref <- c("A","AT","C","G") +alt <- c("AT","A","CA","T") + +em <- canonicalize_mirror_event_ids(chrom, pos, id, ref, alt) + +# Only the mirror pair (2 rows) should survive. +stopifnot(nrow(em) == 2) +stopifnot(all(em$ID %in% c("chr22:100:A:AT","chr22:100:AT:A"))) +# INS keeps pos, DEL is pos+1; single-base event "T". +ins <- em[em$event_type == "INS", ] +del <- em[em$event_type == "DEL", ] +stopifnot(nrow(ins) == 1, nrow(del) == 1) +stopifnot(ins$event_id == "chr22:100:INS:T") +stopifnot(del$event_id == "chr22:101:DEL:T") +# Non-mirror indel and SNP must be absent. +stopifnot(!any(grepl("200", em$event_id))) +stopifnot(!any(em$ID == "chr22:300:G:T")) + +# Empty input -> 0-row result (no error). +em0 <- canonicalize_mirror_event_ids(character(0), integer(0), + character(0), character(0), character(0)) +stopifnot(nrow(em0) == 0) + +cat("UNIT_OK\n") diff --git a/tests/notebooks/reference_data/test_rss_ld_sketch.py b/tests/notebooks/reference_data/test_rss_ld_sketch.py index 853d5ed8e..a28192520 100644 --- a/tests/notebooks/reference_data/test_rss_ld_sketch.py +++ b/tests/notebooks/reference_data/test_rss_ld_sketch.py @@ -3,7 +3,7 @@ Worker: code/script/reference_data/rss_ld_sketch.R (generate_w / process_block). Fixtures (tests/fixtures/rss_ld_sketch), from the MWE: - protocol_example.genotype.chr22.bgz(.tbi) 60-sample chr22 VCF + protocol_example.genotype.chr22.vcf.gz(.tbi) 60-sample chr22 VCF protocol_example.ld_blocks.bed 3 blocks (16-20M / 30-34M / 44-48M) expected/afreq_deterministic.tsv MWE reference .afreq, columns 1-6 @@ -26,7 +26,7 @@ def _run_blocks(run_r, repo_root, tmp_path): """generate_w then process_block over the 3 fixture blocks; return tmp_path.""" worker = repo_root / WORKER - vcf = repo_root / FIX / "protocol_example.genotype.chr22.bgz" + vcf = repo_root / FIX / "protocol_example.genotype.chr22.vcf.gz" W = tmp_path / "W_B50.rds" p = run_r(worker, ["--step", "generate_w", "--n-samples", "60", "--B", B, "--seed", "123", "--output", W]) @@ -34,14 +34,14 @@ def _run_blocks(run_r, repo_root, tmp_path): for bs, be in BLOCKS: p = run_r(worker, ["--step", "process_block", "--vcf", vcf, "--chrom", "chr22", "--block-start", bs, "--block-end", be, "--w-matrix", W, - "--cohort-id", "protocol_example.", "--output-dir", tmp_path, "--B", B]) + "--cohort-id", "protocol_example", "--output-dir", tmp_path, "--B", B]) assert p.returncode == 0, p.stdout + p.stderr return tmp_path def _block_file(tmp_path, bs, be, ext): tag = f"chr22_{bs}_{be}" - return tmp_path / "chr22" / tag / f"protocol_example..{tag}.{ext}" + return tmp_path / "chr22" / tag / f"protocol_example.{tag}.{ext}" def test_rss_process_block_afreq_deterministic(run_r, repo_root, tmp_path): @@ -83,7 +83,7 @@ def test_rss_process_block_via_sos(run_r, run_sos, repo_root, tmp_path): p = run_sos(repo_root / NB, "process_block", { "ld_block_file": repo_root / FIX / "protocol_example.ld_blocks.bed", "chrom": 22, "vcf_base": repo_root / FIX, "vcf_prefix": "protocol_example.genotype.", - "cohort_id": "protocol_example.", "output_dir": tmp_path, "W_matrix": W, "B": B, + "cohort_id": "protocol_example", "output_dir": tmp_path, "W_matrix": W, "B": B, "modular_script_dir": repo_root / "code/script"}) assert p.returncode == 0, p.stdout + p.stderr rows = [] @@ -100,10 +100,48 @@ def test_rss_merge_chrom(run_r, run_sos, repo_root, tmp_path): # + concatenated afreq + the base-R filter summary (ported off data.table). _run_blocks(run_r, repo_root, tmp_path) p = run_sos(repo_root / NB, "merge_chrom", { - "chrom": 22, "output_dir": tmp_path, "cohort_id": "protocol_example.", + "chrom": 22, "output_dir": tmp_path, "cohort_id": "protocol_example", "modular_script_dir": repo_root / "code/script"}) assert p.returncode == 0, p.stdout + p.stderr chrom_dir = tmp_path / "chr22" - assert (chrom_dir / "protocol_example..chr22.pgen").exists() - afreq = (chrom_dir / "protocol_example..chr22.afreq").read_text().splitlines() - assert afreq[0].startswith("#CHROM") and len(afreq) - 1 == 673 # all variants merged + assert (chrom_dir / "protocol_example.chr22.pgen").exists() + afreq = (chrom_dir / "protocol_example.chr22.afreq").read_text().splitlines() + assert afreq[0].startswith("#CHROM") and len(afreq) - 1 == 675 # all variants merged (incl. injected mirror pair) + + +def test_rss_merge_chrom_mirror_event_ids(run_r, run_sos, repo_root, tmp_path): + # The fixture cohort carries one same-position ref/alt-swap indel pair + # (A/AT and AT/A at 16500000). merge_chrom must relabel that pair with + # canonical event IDs and emit the .event_id.tsv sidecar; every other + # variant keeps its standard chr:pos:ref:alt ID. + _run_blocks(run_r, repo_root, tmp_path) + p = run_sos(repo_root / NB, "merge_chrom", { + "chrom": 22, "output_dir": tmp_path, "cohort_id": "protocol_example", + "modular_script_dir": repo_root / "code/script"}) + assert p.returncode == 0, p.stdout + p.stderr + chrom_dir = tmp_path / "chr22" + + # the mirror pair is relabeled in the .pvar: INS keeps pos, DEL is pos+1 + pvar = (chrom_dir / "protocol_example.chr22.pvar").read_text().splitlines() + ids_at_site = [ln.split("\t")[2] for ln in pvar + if not ln.startswith("#") and ln.split("\t")[1] in ("16500000", "16500001")] + assert "chr22:16500000:INS:T" in ids_at_site + assert "chr22:16500001:DEL:T" in ids_at_site + + # the .event_id.tsv sidecar matches the expected mapping byte-for-byte + got = (chrom_dir / "protocol_example.chr22.event_id.tsv").read_text() + exp = (repo_root / FIX / "expected" / "event_id.tsv").read_text() + assert got == exp + + + +def test_canonicalize_mirror_event_ids_unit(run_r, repo_root): + # Direct unit test of the pure canonicalize_mirror_event_ids() function: + # sources the worker and exercises the mirror-pair rule on a tiny in-memory + # fixture -- no pipeline run. A mirror pair is relabeled (INS keeps pos, + # DEL is pos+1); a non-mirror indel and a SNP are excluded; empty input + # yields an empty (0-row) mapping. + checker = repo_root / "tests/helpers/canonicalize_mirror_check.R" + p = run_r(checker, [repo_root / WORKER]) + assert p.returncode == 0, p.stdout + p.stderr + assert "UNIT_OK" in p.stdout