From 912fd8bb40866cc3aedfe68b62c36bd02ea7d510 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Fri, 14 Aug 2026 14:39:43 -0400 Subject: [PATCH 1/6] Harden RSS LD sketch chromosome merging --- code/SoS/reference_data/rss_ld_sketch.ipynb | 143 +++++++++++++- code/script/reference_data/rss_ld_sketch.R | 201 +++++++++++++++++++- 2 files changed, 331 insertions(+), 13 deletions(-) diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index 2646a3baa..a387a998a 100644 --- a/code/SoS/reference_data/rss_ld_sketch.ipynb +++ b/code/SoS/reference_data/rss_ld_sketch.ipynb @@ -21,17 +21,17 @@ "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", + "This pipeline generates a stochastic genotype sample **U = W\u1d40G** 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", + "Rather than storing the full genotype matrix G (n \u00d7 p), the sketch computes U = W\u1d40G (B \u00d7 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\u2013Lindenstrauss 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", + "- G : (n \u00d7 p) \u2014 n individuals \u00d7 p variants\n", + "- W : (n \u00d7 B) \u2014 projection matrix, generated once per cohort\n", + "- U : (B \u00d7 p) \u2014 stochastic genotype sample = W\u1d40G, stored in pgen\n", + "- $\\hat{R}$ : (p \u00d7 p) \u2014 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", + "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). 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. The step also writes an additional directional indel event map without changing the original `.pgen`, `.pvar`, or `.afreq` allele representation.\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", @@ -127,6 +127,14 @@ "chr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n", "```\n", "\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv`** (additional one-to-one event mapping; original variant fields are preserved)\n", + "```\n", + "ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\n", + "chr22:16111042:C:CA 22 16111042 C CA chr22:16111042:INS:A INS 16111042 A\n", + "chr22:16262334:GATTCCATTCGATGATTCCATTCA:G 22 16262334 GATTCCATTCGATGATTCCATTCA G chr22:16262335:DEL:ATTCCATTCGATGATTCCATTCA DEL 16262335 ATTCCATTCGATGATTCCATTCA\n", + "```\n", + "Only length-changing alleles receive canonical `INS` or `DEL` event IDs. SNPs and equal-length substitutions retain their original IDs with `event_type = UNCHANGED`.\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", "\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." @@ -210,7 +218,7 @@ "kernel": "SoS" }, "source": [ - "### Step 2. Process all LD blocks for the chromosome — read the VCF, filter variants, and write per-block dosage sketches U = WᵀG.\n" + "### Step 2. Process all LD blocks for the chromosome \u2014 read the VCF, filter variants, and write per-block dosage sketches U = W\u1d40G.\n" ] }, { @@ -267,7 +275,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 write directional indel event IDs.\n" ] }, { @@ -353,6 +361,27 @@ " --cwd output/rss_ld_sketch" ] }, + { + "cell_type": "markdown", + "id": "event-id-mode", + "metadata": {}, + "source": [ + "### Optional: Regenerate indel event IDs from an existing `.pvar`\n", + "This additive mode reads `.pvar` and writes `.event_id.tsv`. It does not modify `.pgen`, `.pvar`, or `.afreq`.\n", + "\n", + "**Timing**: TBD (on toy dataset)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "event-id-command", + "metadata": {}, + "outputs": [], + "source": [ + "sos run pipeline/rss_ld_sketch.ipynb event_id --final-prefix output/rss_ld_sketch/chr22/protocol_example.chr22\n" + ] + }, { "cell_type": "markdown", "id": "32c022be", @@ -572,7 +601,84 @@ }, "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]\n", + "parameter: chrom = 0\n", + "parameter: output_dir = str\n", + "parameter: cohort_id = str\n", + "parameter: plink2_bin = \"plink2\"\n", + "\n", + "import os, glob\n", + "\n", + "if chrom != 0:\n", + " chroms = [f\"chr{chrom}\"]\n", + "else:\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", + "\n", + "input: for_each = \"chroms\"\n", + "output: f\"{output_dir}/{_chroms}/{cohort_id}.{_chroms}.pgen\"\n", + "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads\n", + "bash: 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, write indel 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" ] }, { @@ -592,6 +698,23 @@ "| `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", "| Region query returns nothing | Missing tabix index | Run `tabix -p vcf file.bgz` to create the `.tbi`. |" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9b262eb", + "metadata": {}, + "outputs": [], + "source": [ + "[event_id]\n", + "# Generate directional minimal-event IDs for an existing chromosome sketch.\n", + "parameter: final_prefix = str\n", + "input: f\"{final_prefix}.pvar\"\n", + "output: f\"{final_prefix}.event_id.tsv\"\n", + "task: trunk_workers = 1, trunk_size = 1, walltime = \"00:30:00\", mem = \"4G\", cores = 1\n", + "bash: expand = \"${ }\"\n", + " Rscript ${modular_script_dir}/reference_data/rss_ld_sketch.R --step event_id --final-prefix \"${final_prefix}\"\n" + ] } ], "metadata": { diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index c96f06fea..b93858140 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -164,9 +164,198 @@ do_process_block <- function(argv) { cat(sprintf("Wrote %d variants -> %s\n", n_passed, basename(paste0(pfx, ".dosage.gz")))) } + +# ---- 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) + + read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", + stringsAsFactors = FALSE) + 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") + + # Create a directional minimal-event identifier without changing VCF alleles. + 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) + } + chrom_col <- if ("#CHROM" %in% names(pvar)) "#CHROM" else "CHROM" + event <- t(mapply(canonical_event, pvar[[chrom_col]], pvar$POS, pvar$ID, + pvar$REF, pvar$ALT, SIMPLIFY = TRUE)) + event_map <- data.frame(ID = pvar$ID, CHROM = pvar[[chrom_col]], + POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, + event, check.names = FALSE, stringsAsFactors = FALSE) + 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") + + # Summarize block-level filtering before removing intermediates. + meta_files <- list.files(chrom_dir, pattern = "[.]meta$", recursive = TRUE, + full.names = TRUE) + if (length(meta_files)) { + 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(" +=== Filter Summary for ", basename(chrom_dir), " === +", sep = "") + print(data.frame(value = unlist(summary), row.names = names(summary))) + } + + # Cleanup happens only after all final outputs and summaries are complete. + 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) +} + +# ---- event_id -------------------------------------------------------------- +do_event_id <- function(argv) { + if (!nzchar(argv$final_prefix)) stop("--final-prefix is required for event_id") + final_prefix <- argv$final_prefix + read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", + stringsAsFactors = FALSE) + pvar <- read_tab(paste0(final_prefix, ".pvar")) + # Create a directional minimal-event identifier without changing VCF alleles. + 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) + } + chrom_col <- if ("#CHROM" %in% names(pvar)) "#CHROM" else "CHROM" + event <- t(mapply(canonical_event, pvar[[chrom_col]], pvar$POS, pvar$ID, + pvar$REF, pvar$ALT, SIMPLIFY = TRUE)) + event_map <- data.frame(ID = pvar$ID, CHROM = pvar[[chrom_col]], + POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, + event, check.names = FALSE, stringsAsFactors = FALSE) + 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") +} + + # ---- 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 / event_id)") +p <- add_argument(p, "--step", help = "generate_w | process_block | merge_chrom | event_id") 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 +373,18 @@ 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") +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 = "") 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 if (identical(argv$step, "event_id")) { + do_event_id(argv) } else { - stop("--step must be 'generate_w' or 'process_block'") + stop("--step must be 'generate_w', 'process_block', 'merge_chrom', or 'event_id'") } From 5aad2564181100a6e663926756a76ffed192fb93 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Mon, 17 Aug 2026 21:14:38 -0400 Subject: [PATCH 2/6] Restrict LD-sketch event IDs to indels --- code/SoS/reference_data/rss_ld_sketch.ipynb | 42 +----------- code/script/reference_data/rss_ld_sketch.R | 73 ++------------------- 2 files changed, 6 insertions(+), 109 deletions(-) diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index a387a998a..4047d34b1 100644 --- a/code/SoS/reference_data/rss_ld_sketch.ipynb +++ b/code/SoS/reference_data/rss_ld_sketch.ipynb @@ -127,13 +127,13 @@ "chr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n", "```\n", "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv`** (additional one-to-one event mapping; original variant fields are preserved)\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv`** (indel-only mapping written automatically by `merge_chrom`; original variant fields are preserved)\n", "```\n", "ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\n", "chr22:16111042:C:CA 22 16111042 C CA chr22:16111042:INS:A INS 16111042 A\n", "chr22:16262334:GATTCCATTCGATGATTCCATTCA:G 22 16262334 GATTCCATTCGATGATTCCATTCA G chr22:16262335:DEL:ATTCCATTCGATGATTCCATTCA DEL 16262335 ATTCCATTCGATGATTCCATTCA\n", "```\n", - "Only length-changing alleles receive canonical `INS` or `DEL` event IDs. SNPs and equal-length substitutions retain their original IDs with `event_type = UNCHANGED`.\n", + "Only length-changing alleles that canonicalize to `INS` or `DEL` are included. SNPs, equal-length substitutions, and symbolic alleles are omitted from this mapping.\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", "\n", @@ -361,27 +361,6 @@ " --cwd output/rss_ld_sketch" ] }, - { - "cell_type": "markdown", - "id": "event-id-mode", - "metadata": {}, - "source": [ - "### Optional: Regenerate indel event IDs from an existing `.pvar`\n", - "This additive mode reads `.pvar` and writes `.event_id.tsv`. It does not modify `.pgen`, `.pvar`, or `.afreq`.\n", - "\n", - "**Timing**: TBD (on toy dataset)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "event-id-command", - "metadata": {}, - "outputs": [], - "source": [ - "sos run pipeline/rss_ld_sketch.ipynb event_id --final-prefix output/rss_ld_sketch/chr22/protocol_example.chr22\n" - ] - }, { "cell_type": "markdown", "id": "32c022be", @@ -698,23 +677,6 @@ "| `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", "| Region query returns nothing | Missing tabix index | Run `tabix -p vcf file.bgz` to create the `.tbi`. |" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b9b262eb", - "metadata": {}, - "outputs": [], - "source": [ - "[event_id]\n", - "# Generate directional minimal-event IDs for an existing chromosome sketch.\n", - "parameter: final_prefix = str\n", - "input: f\"{final_prefix}.pvar\"\n", - "output: f\"{final_prefix}.event_id.tsv\"\n", - "task: trunk_workers = 1, trunk_size = 1, walltime = \"00:30:00\", mem = \"4G\", cores = 1\n", - "bash: expand = \"${ }\"\n", - " Rscript ${modular_script_dir}/reference_data/rss_ld_sketch.R --step event_id --final-prefix \"${final_prefix}\"\n" - ] } ], "metadata": { diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index b93858140..802978f16 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -245,6 +245,7 @@ do_merge_chrom <- function(argv) { event_map <- data.frame(ID = pvar$ID, CHROM = pvar[[chrom_col]], POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, event, check.names = FALSE, stringsAsFactors = FALSE) + event_map <- event_map[event_map$event_type %in% c("INS", "DEL"), , drop = FALSE] duplicated_event <- duplicated(event_map$event_id) | duplicated(event_map$event_id, fromLast = TRUE) if (any(duplicated_event)) { @@ -289,73 +290,9 @@ do_merge_chrom <- function(argv) { unlink(block_dirs, recursive = TRUE) } -# ---- event_id -------------------------------------------------------------- -do_event_id <- function(argv) { - if (!nzchar(argv$final_prefix)) stop("--final-prefix is required for event_id") - final_prefix <- argv$final_prefix - read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", - stringsAsFactors = FALSE) - pvar <- read_tab(paste0(final_prefix, ".pvar")) - # Create a directional minimal-event identifier without changing VCF alleles. - 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) - } - chrom_col <- if ("#CHROM" %in% names(pvar)) "#CHROM" else "CHROM" - event <- t(mapply(canonical_event, pvar[[chrom_col]], pvar$POS, pvar$ID, - pvar$REF, pvar$ALT, SIMPLIFY = TRUE)) - event_map <- data.frame(ID = pvar$ID, CHROM = pvar[[chrom_col]], - POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, - event, check.names = FALSE, stringsAsFactors = FALSE) - 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") -} - - # ---- CLI ------------------------------------------------------------------- -p <- arg_parser("RSS LD random-projection sketch (generate_w / process_block / merge_chrom / event_id)") -p <- add_argument(p, "--step", help = "generate_w | process_block | merge_chrom | event_id") +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") @@ -383,8 +320,6 @@ if (identical(argv$step, "generate_w")) { do_process_block(argv) } else if (identical(argv$step, "merge_chrom")) { do_merge_chrom(argv) -} else if (identical(argv$step, "event_id")) { - do_event_id(argv) } else { - stop("--step must be 'generate_w', 'process_block', 'merge_chrom', or 'event_id'") + stop("--step must be 'generate_w', 'process_block', or 'merge_chrom'") } From 2738bca10736448191e9b16b475e150a3f2a0f23 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Tue, 18 Aug 2026 15:37:21 -0400 Subject: [PATCH 3/6] Restrict LD-sketch event IDs to mirror-pair indels Narrow event-ID assignment from all indels to only indels that form an exact REF/ALT-swap pair at the same position (the ambiguous cases where INS vs DEL anchoring can flip effect-allele orientation). Non-mirror indels, SNPs, and equal-length substitutions retain standard IDs. merge_chrom now substitutes the canonical event ID into .pvar and .afreq for the mirror set, and writes no .event_id.tsv when a panel has no mirror pairs (leaving .pvar/.afreq fully standard). Notebook narrative updated to match. Validated against deployed R5 EUR chr1 (reproduces 48014 mirror event IDs exactly). --- code/SoS/reference_data/rss_ld_sketch.ipynb | 158 +------------------- code/script/reference_data/rss_ld_sketch.R | 58 +++++-- 2 files changed, 51 insertions(+), 165 deletions(-) diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index 4047d34b1..5d99c6b9b 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\u1d40G** 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 \u00d7 p), the sketch computes U = W\u1d40G (B \u00d7 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\u2013Lindenstrauss lemma. G is never stored.\n", - "\n", - "**Matrix dimensions:**\n", - "- G : (n \u00d7 p) \u2014 n individuals \u00d7 p variants\n", - "- W : (n \u00d7 B) \u2014 projection matrix, generated once per cohort\n", - "- U : (B \u00d7 p) \u2014 stochastic genotype sample = W\u1d40G, stored in pgen\n", - "- $\\hat{R}$ : (p \u00d7 p) \u2014 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). 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. The step also writes an additional directional indel event map without changing the original `.pgen`, `.pvar`, or `.afreq` allele representation.\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\u1d40G** 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 \u00d7 p), the sketch computes U = W\u1d40G (B \u00d7 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\u2013Lindenstrauss lemma. G is never stored.\n\n**Matrix dimensions:**\n- G : (n \u00d7 p) \u2014 n individuals \u00d7 p variants\n- W : (n \u00d7 B) \u2014 projection matrix, generated once per cohort\n- U : (B \u00d7 p) \u2014 stochastic genotype sample = W\u1d40G, stored in pgen\n- $\\hat{R}$ : (p \u00d7 p) \u2014 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 (\u201cmirror\u201d) pair at the same position \u2014 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\u2194event mapping is recorded in a `.event_id.tsv`. All other variants \u2014 SNPs, equal-length substitutions, and non-mirror indels \u2014 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." ] }, { @@ -86,58 +66,7 @@ "kernel": "SoS" }, "source": [ - "## Output\n", - "\n", - "- **`output/rss_ld_sketch/W_B50.npy`**\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", - "```\n", - "\n", - "- **`output/rss_ld_sketch/chr22/chr22__/.chr22__.dosage.gz`** (`process_block`)\n", - "(the per-LD-block sketch: B pseudo-samples x the variants in that block, gzipped dosage text. One directory per block, written before `merge_chrom` concatenates them.)\n", - "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pgen`** (with `.pvar`, `.psam`, `.afreq`)\n", - "(per-chromosome PLINK2 genotype sketch assembled by `merge_chrom`: B pseudo-samples x p variants. The `.pgen` is binary; the three companion files describe it.)\n", - "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pvar`** (variant information)\n", - "\n", - "```\n", - "#CHROM\tPOS\tID\tREF\tALT\n", - "22\t16073625\tchr22:16073625:G:T\tG\tT\n", - "22\t16102937\tchr22:16102937:C:T\tC\tT\n", - "```\n", - "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.psam`** (the B sketch pseudo-samples, named `S1`..`S{B}`)\n", - "\n", - "```\n", - "#FID\tIID\tSEX\n", - "S1\tS1\tNA\n", - "S2\tS2\tNA\n", - "```\n", - "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.afreq`** (allele frequencies and the sketch value range per variant)\n", - "\n", - "```\n", - "#CHROM\tID\tREF\tALT\tALT_FREQS\tOBS_CT\tU_MIN\tU_MAX\n", - "chr22\tchr22:16073625:G:T\tG\tT\t0.108333\t120\t-2.588204\t2.248544\n", - "chr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n", - "```\n", - "\n", - "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv`** (indel-only mapping written automatically by `merge_chrom`; original variant fields are preserved)\n", - "```\n", - "ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\n", - "chr22:16111042:C:CA 22 16111042 C CA chr22:16111042:INS:A INS 16111042 A\n", - "chr22:16262334:GATTCCATTCGATGATTCCATTCA:G 22 16262334 GATTCCATTCGATGATTCCATTCA G chr22:16262335:DEL:ATTCCATTCGATGATTCCATTCA DEL 16262335 ATTCCATTCGATGATTCCATTCA\n", - "```\n", - "Only length-changing alleles that canonicalize to `INS` or `DEL` are included. SNPs, equal-length substitutions, and symbolic alleles are omitted from this mapping.\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", - "\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." + "## Output\n\n- **`output/rss_ld_sketch/W_B50.npy`**\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```\nfloat32 (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```\n\n- **`output/rss_ld_sketch/chr22/chr22__/.chr22__.dosage.gz`** (`process_block`)\n(the per-LD-block sketch: B pseudo-samples x the variants in that block, gzipped dosage text. One directory per block, written before `merge_chrom` concatenates them.)\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pgen`** (with `.pvar`, `.psam`, `.afreq`)\n(per-chromosome PLINK2 genotype sketch assembled by `merge_chrom`: B pseudo-samples x p variants. The `.pgen` is binary; the three companion files describe it.)\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pvar`** (variant information)\n\n```\n#CHROM\tPOS\tID\tREF\tALT\n22\t16073625\tchr22:16073625:G:T\tG\tT\n22\t16102937\tchr22:16102937:C:T\tC\tT\n```\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.psam`** (the B sketch pseudo-samples, named `S1`..`S{B}`)\n\n```\n#FID\tIID\tSEX\nS1\tS1\tNA\nS2\tS2\tNA\n```\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.afreq`** (allele frequencies and the sketch value range per variant)\n\n```\n#CHROM\tID\tREF\tALT\tALT_FREQS\tOBS_CT\tU_MIN\tU_MAX\nchr22\tchr22:16073625:G:T\tG\tT\t0.108333\t120\t-2.588204\t2.248544\nchr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n```\n\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```\nID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\nchr22:36543767:A:AT 22 36543767 A AT chr22:36543767:INS:T INS 36543767 T\nchr22:36543767:AT:A 22 36543767 AT A chr22:36543768:DEL:T DEL 36543768 T\n```\nOnly indels that form an exact REF/ALT-swap pair at the same position are included \u2014 the cases where insertion vs. deletion anchoring is ambiguous. The example above shows one such pair: the same single-base event at chr22:36543767 represented as both an insertion (`A`\u2192`AT`) and a deletion (`AT`\u2192`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\nThese 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\nThe `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." ] }, { @@ -275,7 +204,7 @@ "kernel": "SoS" }, "source": [ - "### Step 3. Merge the per-block dosage sketches, reconcile allele frequencies, and write directional indel event IDs.\n" + "### Step 3. Merge the per-block dosage sketches, reconcile allele frequencies, and relabel mirror-pair indels with canonical event IDs.\n" ] }, { @@ -580,84 +509,7 @@ }, "outputs": [], "source": [ - "[merge_chrom]\n", - "parameter: chrom = 0\n", - "parameter: output_dir = str\n", - "parameter: cohort_id = str\n", - "parameter: plink2_bin = \"plink2\"\n", - "\n", - "import os, glob\n", - "\n", - "if chrom != 0:\n", - " chroms = [f\"chr{chrom}\"]\n", - "else:\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", - "\n", - "input: for_each = \"chroms\"\n", - "output: f\"{output_dir}/{_chroms}/{cohort_id}.{_chroms}.pgen\"\n", - "task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads\n", - "bash: 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, write indel 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" + "[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" ] }, { @@ -715,4 +567,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 802978f16..687079e20 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -246,19 +246,53 @@ do_merge_chrom <- function(argv) { POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, event, check.names = FALSE, stringsAsFactors = FALSE) event_map <- event_map[event_map$event_type %in% c("INS", "DEL"), , drop = FALSE] - 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") + + # 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 <- event_map[sk %in% nk, , drop = FALSE] + + # 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) { + 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") + + message("EVENT_ID_SUBSTITUTION mirror_variants=", nrow(event_map)) + } else { + message("EVENT_ID_SUBSTITUTION mirror_variants=0 (no ambiguous indels; panel left standard)") } - 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") # Summarize block-level filtering before removing intermediates. meta_files <- list.files(chrom_dir, pattern = "[.]meta$", recursive = TRUE, From 6184b60047b041951dcb994bac1fd3e43b503616 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 19 Aug 2026 09:54:36 -0400 Subject: [PATCH 4/6] Add mirror-pair fixture + test; accept .vcf.gz in VCF discovery Test coverage for the mirror-pair event-ID feature: - process_block VCF discovery now accepts .vcf.gz in addition to .bgz (rss_ld_sketch.R); both extensions validated. - Augment the rss_ld_sketch test fixture with one same-position REF/ALT-swap indel pair (A/AT + AT/A at chr22:16500000, ~40% AF, well-called) injected into the existing 60-sample cohort; convert fixture to .vcf.gz. - Regenerate expected/afreq_deterministic.tsv (673->675 variants) and add expected/event_id.tsv (the mirror-pair mapping). - Update test_rss_ld_sketch.py: .bgz->.vcf.gz, single-dot cohort-id, 673->675, and a new test_rss_merge_chrom_mirror_event_ids asserting the INS:T/DEL:T relabeling and the .event_id.tsv sidecar. - Notebook narrative: name the .vcf.gz fixture, note both extensions accepted, show the mirror pair in the Output example, correct stale W_B50.npy -> .rds. All 5 tests in test_rss_ld_sketch.py pass. --- code/SoS/reference_data/rss_ld_sketch.ipynb | 61 ++++++++++++++++-- code/script/reference_data/rss_ld_sketch.R | 2 +- .../expected/afreq_deterministic.tsv | 2 + .../rss_ld_sketch/expected/event_id.tsv | 3 + .../protocol_example.genotype.chr22.bgz.tbi | Bin 20170 -> 0 bytes ...=> protocol_example.genotype.chr22.vcf.gz} | Bin 304997 -> 305100 bytes ...protocol_example.genotype.chr22.vcf.gz.tbi | Bin 0 -> 20172 bytes .../reference_data/test_rss_ld_sketch.py | 43 +++++++++--- 8 files changed, 94 insertions(+), 17 deletions(-) create mode 100644 tests/fixtures/rss_ld_sketch/expected/event_id.tsv delete mode 100644 tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.bgz.tbi rename tests/fixtures/rss_ld_sketch/{protocol_example.genotype.chr22.bgz => protocol_example.genotype.chr22.vcf.gz} (94%) create mode 100644 tests/fixtures/rss_ld_sketch/protocol_example.genotype.chr22.vcf.gz.tbi diff --git a/code/SoS/reference_data/rss_ld_sketch.ipynb b/code/SoS/reference_data/rss_ld_sketch.ipynb index 5d99c6b9b..8315e4130 100644 --- a/code/SoS/reference_data/rss_ld_sketch.ipynb +++ b/code/SoS/reference_data/rss_ld_sketch.ipynb @@ -19,7 +19,7 @@ "kernel": "SoS" }, "source": [ - "## Overview\n\nThis pipeline generates a stochastic genotype sample **U = W\u1d40G** 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 \u00d7 p), the sketch computes U = W\u1d40G (B \u00d7 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\u2013Lindenstrauss lemma. G is never stored.\n\n**Matrix dimensions:**\n- G : (n \u00d7 p) \u2014 n individuals \u00d7 p variants\n- W : (n \u00d7 B) \u2014 projection matrix, generated once per cohort\n- U : (B \u00d7 p) \u2014 stochastic genotype sample = W\u1d40G, stored in pgen\n- $\\hat{R}$ : (p \u00d7 p) \u2014 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 (\u201cmirror\u201d) pair at the same position \u2014 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\u2194event mapping is recorded in a `.event_id.tsv`. All other variants \u2014 SNPs, equal-length substitutions, and non-mirror indels \u2014 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." + "## 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." ] }, { @@ -41,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", @@ -53,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)" @@ -66,7 +66,54 @@ "kernel": "SoS" }, "source": [ - "## Output\n\n- **`output/rss_ld_sketch/W_B50.npy`**\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```\nfloat32 (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```\n\n- **`output/rss_ld_sketch/chr22/chr22__/.chr22__.dosage.gz`** (`process_block`)\n(the per-LD-block sketch: B pseudo-samples x the variants in that block, gzipped dosage text. One directory per block, written before `merge_chrom` concatenates them.)\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pgen`** (with `.pvar`, `.psam`, `.afreq`)\n(per-chromosome PLINK2 genotype sketch assembled by `merge_chrom`: B pseudo-samples x p variants. The `.pgen` is binary; the three companion files describe it.)\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pvar`** (variant information)\n\n```\n#CHROM\tPOS\tID\tREF\tALT\n22\t16073625\tchr22:16073625:G:T\tG\tT\n22\t16102937\tchr22:16102937:C:T\tC\tT\n```\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.psam`** (the B sketch pseudo-samples, named `S1`..`S{B}`)\n\n```\n#FID\tIID\tSEX\nS1\tS1\tNA\nS2\tS2\tNA\n```\n\n- **`output/rss_ld_sketch/chr22/protocol_example.chr22.afreq`** (allele frequencies and the sketch value range per variant)\n\n```\n#CHROM\tID\tREF\tALT\tALT_FREQS\tOBS_CT\tU_MIN\tU_MAX\nchr22\tchr22:16073625:G:T\tG\tT\t0.108333\t120\t-2.588204\t2.248544\nchr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n```\n\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```\nID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt\nchr22:36543767:A:AT 22 36543767 A AT chr22:36543767:INS:T INS 36543767 T\nchr22:36543767:AT:A 22 36543767 AT A chr22:36543768:DEL:T DEL 36543768 T\n```\nOnly indels that form an exact REF/ALT-swap pair at the same position are included \u2014 the cases where insertion vs. deletion anchoring is ambiguous. The example above shows one such pair: the same single-base event at chr22:36543767 represented as both an insertion (`A`\u2192`AT`) and a deletion (`AT`\u2192`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\nThese 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\nThe `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." + "## Output\n", + "\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", + "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", + "(the per-LD-block sketch: B pseudo-samples x the variants in that block, gzipped dosage text. One directory per block, written before `merge_chrom` concatenates them.)\n", + "\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pgen`** (with `.pvar`, `.psam`, `.afreq`)\n", + "(per-chromosome PLINK2 genotype sketch assembled by `merge_chrom`: B pseudo-samples x p variants. The `.pgen` is binary; the three companion files describe it.)\n", + "\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.pvar`** (variant information)\n", + "\n", + "```\n", + "#CHROM\tPOS\tID\tREF\tALT\n", + "22\t16073625\tchr22:16073625:G:T\tG\tT\n", + "22\t16102937\tchr22:16102937:C:T\tC\tT\n", + "```\n", + "\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.psam`** (the B sketch pseudo-samples, named `S1`..`S{B}`)\n", + "\n", + "```\n", + "#FID\tIID\tSEX\n", + "S1\tS1\tNA\n", + "S2\tS2\tNA\n", + "```\n", + "\n", + "- **`output/rss_ld_sketch/chr22/protocol_example.chr22.afreq`** (allele frequencies and the sketch value range per variant)\n", + "\n", + "```\n", + "#CHROM\tID\tREF\tALT\tALT_FREQS\tOBS_CT\tU_MIN\tU_MAX\n", + "chr22\tchr22:16073625:G:T\tG\tT\t0.108333\t120\t-2.588204\t2.248544\n", + "chr22\tchr22:16102937:C:T\tC\tT\t0.308333\t120\t-1.849294\t2.069429\n", + "```\n", + "\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", + "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)`." ] }, { @@ -126,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" ] } @@ -147,7 +194,7 @@ "kernel": "SoS" }, "source": [ - "### Step 2. Process all LD blocks for the chromosome \u2014 read the VCF, filter variants, and write per-block dosage sketches U = W\u1d40G.\n" + "### Step 2. Process all LD blocks for the chromosome — read the VCF, filter variants, and write per-block dosage sketches U = WᵀG.\n" ] }, { @@ -523,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", diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index 687079e20..35a775b47 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)) } 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 bb13bb2ae9ccdee482a8b371f4408fe11e4a4f9a..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 20170 zcmZ^qdpMK-AOAb6RL-WHLUOirKsl@&E2m0D5jCIECz+PX7+WRT0qck&3!#J3Qr#LO z+4vMShl;w*VJch9EH%?MesAC3f4}Scb#-0->2kGu_x*mp9*^hiHl;0HG;i+D?2>t2 z-V5f*Jjd8c?FVB8n4CKBXOtFjK`=fZWkl73ylf8l@%icws~q-#pKR6N_7_in#BH4Z zzHX$baaGXeEphA#cfx=o&*PuiPl+>DC3(zBe04vyTpqI&M9?XjIoR z{OMWL-{ci5%NhAcjN~rUZHp8=C@XJa_mjsV-7bR&6B;UOziIdJprO|cx55WWujY|T z@;$m9pLlK9e!9p0wFUE|DqZYNu@cE1lyJLogj-7mt68cex>9$3@wVIe3m6H)yM8x> z%(3ae2Zyw3Cv=LNhRKjf*0+n}hEsh0CP-(g&MlYD=QnJsli$W5_wGpxfRdM^-G#Dz zU(OUxaq*f!jpec5S8BnJdVWDp!Jt0vRR%&DBHNi=(8y2A8DB_ZvhFhdF^04`JbOEo zY~tv557x-lwdfRosi0gH$r=i|);Pkb1$qU`qaW+Eo?ra4%7BJJCGE`)fYj>RM>@r2 zPbgB6?8p>n5smytF%V5Yy_jN3zjoV?PE!U_9r7j< zIdNImx7WU|R`p4pv^c+T1Z6MUlBG&oz}(0$_AC$2#6;intJ(p@^45O72Y3I_jq4P* zO_GyDvYHfbDvt1sf6j- znggBUn=NFCNS2MFnx@^6|J)_;XQ^a?QXPInY=it3ruyk_<#i0>=`JPs%2)v_^s~|p zs+BZ@B8)`w$lxL-I`NXU3{EuG@$MAwh~;$P6yuX<1D2{|i<3RS;a9FfE`~R=eMJZq zO)~ep3nw?A8@k16p46-0`!Rb{zHGxre!pwIK_+Iv*gN~D;!WiuN~chE$(Jj{5pp)p z-7^GVr#uOwaqofbAc%C@zNG}_z0@7<6#s##%Y?F(;#W~D)gc!rExyU-@W8Q=+wt39mqJC_!mm zQ4W{+N4jdul zq1|pZaAz?Fwu3vX=aj8}IUMmqiv}KqRti@nyQss}#wmEAoNG8j~&+|2J! zY>?~mXRuepKaUz$zua%{4>69I5yk2Ac&x+rA&~!o&B42Hk66#HORTkud{Zck_vgIB zDN1$FdMq|fVA}GDzb&Rl5?59mmmc~)@6l{cLYrY0B2Ec@f=4VamCreVqp zk!(4IyRF}JBI#ca%fwX4t4DjV(og)$Y(SYdcI-Yqn7$&RFZi@j@h18l!mv}ki|}j2 zy?R0tOAZpst|f9SaD-?L71%y2T{6 zF5Dpw@}$L!WNlh%E=!*p^P@Ju>@CL3Zo1Xm#8anV-8CkhTuQojWXuOjZgV&%OKS1fgetIL?z?9Poa}p>@&tpHU7-A2ff{sYO{TQW4;i5Y^|&U z9+qhh`~+CZ0BGD1lg$%Fx1Vn1mDdc>K*5uOGK<}obfR~`$tvs9`bQa2KkIGv2W zGP0~+@Y1yaBl$&^l)$0s+F4yE1zhe=Xw!?D%@ivQW-T~i8T72 zz7UzW*GMekRvy=H?GW#}J}Sv-oNBnhd%rff9F~k|wswlcXf#xZ_>LWUSR|{v!!;Gj z+*F)pIKsnY0w4OOeoU$BTZtk54TJL*oCOC2x-32^+Ouc)?Lc8xitUFo*w7P^-z6@K zBDV=;9}>9VP1Mp_+_0AXs4eYMt@%&NVLN|amkx1&Ijv44vsuPnh*MNkI2ky?n|lH? zx?%s8SNH-<$kn~>yP>f5bu)fY{s-htFjQ%6ZZbJCNXpZt6G$y_M9bp4uz#5L>|NL; zOh2SctaFN5DU^*l%mvQODiApANEfL+pIBds&&8PbhPi`V(bzF+Vw#+PXKmCyxGqE& zNoeFwD%VAUFD+Q7`fvT~n0Ec@>H21dy=JwdzB%zRBB4WU`I71)lHEhmpnBfmBuOg) z9L%}2WM$Na1f_R<0OV3<-(Lo|9?_C?i^Uj>Bi)W7rVptxt z5iPvQMG4a2w4alUO=qLh{Ggsyou|;}c_Nva85iLxnOL@$n_tY62?cfxg*^9!AdZ!O z*if3oQ?DmkNgeo})pR*Hwzln03xVj7yT?Na(cAyF2@3OLzAZ}uJ> z-Jm!7%XRS9MaFGo`-zaqB*m30=rbU{PC6~(7uyvXT*C~k*}K9QYNb0|u7F3&^_n}x zZqF#kL^5^?_Y6)^wHJ+NsU}W3E#McUs|>DUs$GM#eW2>wwtf|`Rl9Z>xGp2Ky&~DK zHqJ>LA#V#B$6|jw;dB|L#$X=E-UXTVt#2uXt%?vaifvL&7p^N#ag%_KWU(<{?|z@( zQU=%k&>rp<*Cf)wY58t3=Mqlw5`i{nd9=hiE#pUh%PF{lVaS35{h(PDL(_Q?F(?>n zrL3>G3y&Pv4($>j8GQ*G94{?Ykk3qW{w64HF%q+%X*C0e3Ezzx#n$ zHu)awlrK@2;P@?G_fUj5R(d4fN7uHaGJ!a$_FSHX&m1e_q<$JPO4n-SF z82HX^aA>kstCon1CFl@VddCJP zf=^5+>Mg_!pmvOTL&=Ab!~lr*b*&71q|*pR%8|ijqQ0VAJhY0uP$+B3;#>gv=^{wR z5t=pnT=)&^KgtR7SIN?vn=|#Y5TuE%G!n)!w-Q!+p-=i<|EBr* zh^J7N$LCUTigpQi7fzA;Lhuhu_2nRzI%VfB+kt)9W86RVRDPXzTq)ASuHnyoN@}5F zqEkYPhNg^^f4#?>74u9~Y*!j)R*|qgN{x6&#UqgG19h>vQxVGa6<9}pKPuOijv?Lh zDUu`_Np5;YZih(i){gh!>v?*JY4+mT@V${^19ZkCN5>L4;e}Rys__tg8x5Sfhxpv@ z{UH+vyt#Tf0)LkP3>Llr%oOpk^#XV77JgJJyY@0hQsW;9rcIQ^lS%pTA4ke4_}UPr zWJ96MFNBkWQuR$w$i5IImqCZb$eEwKRz{lDcsbmf zjhN~b>r_$|BH6x8b0OeT!^tplQT~aPUKxm((4ldgiKr$NQj3mn(vPf^Zsa$VHQ>1; zh$!Yw{3nIROqiQ=4c{{+8-E9!gi-ZFyrzsG{e-`1O*98BKreBhsSx35+z= zm2`<$B56}XnGv0%z!ADN1Pdir@{TJK?M3n{7@l!(nh)e}W8ql=tAEyHbcoq&$RX*g z1g(Y?o=~>oF?WQ>Xi7-ZjFv|>yG?0FcZpjYb)u&yrFzZW9-M+K;kw`mKi&!=Sn1*r zfg6i0vX`=3SJR&bDhIL6h$>xFmw4bY0@W#Qg=mcQXn7x7z`!Y9Hw&V1ivM)bPw`_h z^7q$4B6tk@f<%aV{#Uel{^fgc+h%<*Fi#z(l!#<9F=qo#(Z32E&PtCAaWdi;r(ZEB z#tgXb2RwyRMrH4WOpDg~mBRE89gA*pQY>Y!P-aKt-p3KHzCdqgvH9Daw(*P4<{M;S z7zS>~Wb7&-`yoHnL2BVBjnr7P0!zsfa#eiq7+=kBq+VweIST{JRwtZdM?}U zwSH}nBaELxA7rJ+`8XN!8@3b|@Qnwq>}%fv(L<0w?!l<GF_(m~$t#Bqar$k5%l=ZLyTzlI$=^k?icOs5I7N*w+E=m=uW|rPm)=91W4f{0 zaIQkHc!Uyu6N8*NT=iolvwGQ!hkoZjd3=yhxT06G*oVxS>-nCRH^XTdQsJR0v1+yH z?T{*O2zA>JS&+om4L9s(l`$GP+L8pTt$yT$MB`~5`4+c+B zNRbZQDc<5p5wnp)Bh^%6p{xz!PhC;Ck)@Z=<=`8(PWk?>poh8mF~ZpX~i#I|EQOSp{u;cf*=IJMpgp5)oVAF z!-XVGMVFW~yYd`7Pthqp-b_)6WDoap=i>;oC&7&E5pY@xo+*qoR+CROYXlQEZ~rbM z0k);*pkCHrkNFAv0*w41jfAD_=%JKCiiJ>C_n31Pr|8ZVShGBi#bP%D1MtZ8GNyX- z9%V7c)cLSc0F+R;Hg{67#PfKDJY(vByF?@uZY*(G)%(7-d@xJtjWQ@kvT(hcT`v8#d=EOtvc zwqnZC{auK_o~3GCoo2!(VlTL6VNfrFiUOc8*b(#?F6)O+}r$wiD!d)-GTS81RE-2I`(NQdG|7?q%n><~`{ zPl7=yJeC^(hO|k+aTdG%fOI*(e;u}6G(YN(Vqt@G7XH1&Q2tpdNUHWvi_wxU;d?S~ z;ICs);XCSmA!OD1{tCD?NmtS#cB&&Mi)7|UxP)Q`{#QRRVII zaxzF{_vX@hrrqK=FG`3|cIFtj5~nCwD%iyGuzmU!jm z+BI^uNG7tG%Z{$+AUixH;ZpFuHERe6EVg`|bfL<0c9&8(sb1|C4R}H$VlChaiM=-9 zi5PpW*-{v;Mr3efjfeOzDNQ2TyE?8Wj?mkT-pNXzl%T<_|8tqsW}q*g8-QE?Ofg{6 zu>PaOPVqaNt$RIL&7jro{1tXKQ4o9dkNAF{$?{son%_ly%R8WMej@%JAl#*Nh$fn{R zZ1GVmv{P(mNO=I*C!Vtzr?@o^Odh*wUJE!ao%aWZKv{)$CxF3ttiS2+U+Ev%0L5bW z?*$ZLBE1AC!Zdj|pokup?TIoN*{@mCEzTJyZxzbkPjQOC7(wSqczznt}CTYi)5}$c_L0xC5CrLOust|g*VyQn zz%{|zv@UT&1l(&9Q=jiy|0*2rq4=M6_$>82kC;7J z1!hVE`_E-?cebWshgi`{TPc#Q)8YEz6sDn^798QR2(2SQ;d{5DLjg}D=271ZW&7>8 zKX8f{FWg8E{R6?0sC>L0mA(3F- z0HPngZ1|^YzU1j)<%dyYQbVwj50uZhUR41PBSZ(zL%++%>$G6Vp*Kkfy@ z;rZ|CnK&`NZKW5DXfu*Q!qOboJaE^Py0X#$pe9VCVB^ z1SPc{ig{7vlv;4e{MsEE2tmt2h1l}}LfNgooN{9I>+7wW`JH0D7;2kH#uIV$aEhPp zARRn{%|JTzzpx|aVo0!eq&Gx@=LjQq&DtkG#IDro=@xgk&hf?DVNOnQ0zQA2U_DFK zzXDq`b)f#wJ}hlSuRlJg_9li=bFlA6qUmh>o!UDXUZZ#4(DY@Aa&>CM(8v@G1cb3a zmjPyhRxTe$__|2|MzC*|QUpu&BV8KCQhncqMS`P+UCSF;Iy3d|ujU|;l^$d*oyRwE ztioTxL|>0`|29B(d2eG<0bA(m&^yGw@2LJFSyc{qKqwPBa0fgFCU(VhKsbE|p8GLQWV?@uvwGCy&_;(j10GL*hO8XN5z5dKFMk3Dei*O_M7n8~v5S(N-9# z2ZoTx&J^%1tZt-{y2Pf5>jI4i0a~<~qn|DsYnI>NJtv#|^ybjwZZ42b$liTGHc7tP z0ohb%stsgQqp?;qu$^@3kWe<$%ALR|Sc&LGfMkw4ne&Mww+$$mf#T;Wf9CZEl&KN0X@`WeJ9ju$#k2CsVze_)N10A?s_Q`#iR6L|3~x(h`!3K`&7H}e zn>c;NDE!B*fGgOa+yPgZBr^+m;L1aO?UIz-*nVrE?6czafwHfWQFh=I@WMH)_-;0b z6+$Y&3VKqA{cgxV%^Xlgp@$BwTdas8TMA|43EXS_kkppJ?ZPP@b)fXsC+t{kFIVOY zev~5xBt*z{D!46@HqGHg*fZ*P4~F#ijx-d%;o(z|54^pRARne4n}K}DaM60yAs&jL z9T3UZ$8%5N2!#oppFj#=&?_aA_~Ds8Q$BIj4Wxs~W}hmc&P**om%ti1$LmG1ejN@d zki-nmU!0=jU-Uv2`^`pX$;iTDVoom^6aAuCA)9=ba4gbe6`pL24h7?`T@Vq3;<;yBz3{cPa|uIrxp(TL(#C~_+;8^XvS>35{8TPG&;m# z5{jCKhP16DD{0xp6`odol&c@;byZr^vA)S=>ra zA~9D^58_Vt=TEpgaMf{y*3SZ8Fp9WIwfRxZO8l46HM5fuQfkkZ;xRSpLtb3HZ2veo z?A#;jgQ4gn=KT|C&B+-31VBwPDg=(uEDj7F>|xN%rhi!sn%Uy1mw29Jt2~+q1eD@8 z5Ktf+t)ymr&&}T;)k5H(=s9aWE8rMkE#RMcXljT^wh=4e%Aaw(Ais#=-Po753u+BO zHk87f{?P%4!mD4@1fk4=%YlHELy*3Y?t@~Lhd~vg@nkLWE(kI$Ae}mI(->(pg7U94 z5Dt=5DhP*TmuaL<@!l3{llDWnE)bE?EzXLinh9l?QO*@G@|Mk!!_>ck9QwyhdIF$< z6o<>@FoWB_&Gz*fusVPTMFjioK2StNA#*?x(Re#5+XotXzYY`;lH)p{h!}oI21Ud_ zf4C<>sptoai1Zy+PD}U=rCA2oG0_{Nvwfk2ujVZkaM@|ZaH{8^Lo!tgj9)BQ4k8jp zFJP%A7dn}NR`>zP1k`L`U?6B_uAVu2e%JrV>K>TIPNO%nRJ&80wt_kFMZrZ(^}vBZ zP#FHStSE(JR%(WJi+Ne(e}%HM&D`feLo?CaSnNR`Cl7Fd-v+Ayq-WdPcS4O3jwecC zUbWsPf>K%+!2zzRkB>e`LN-X)1vx-g3f>mYtC^iN`kwDi0{{VBYY`~X$!JGV=eavU zt*cSRK7bx70~`T8B;3;n^zdjS8PLOqTigXW!l%ADdIMFFbmWfpfijvyLF`-1fOm#NJ!{&Ow-7ABf>pvnsxVZTk2PG1~g z+Z~)_HQQzNH=4ix>De=JREx6)2E@t@q9F`H>4q|6P&N|GKO@{d&~>6YOn;*z;M^9n=X12#*KVs`FgxRK}ur+av&&u z=DLa#`v5ugyKMnWEcf(4g)7=b9stwduel00h07&^TcMge<2fKz+rpWTQ)q~}psL7S zHCI)Pn+VpjRF8;)qoC8B5@@p619+@8KPr~$nvIE8TnC}h`1A(06i%$t^am5-6t%3r zaY*X}6{NylSYW{Ou5Y!F5*U^8%VC0W9?88c0t5q5a@`derz_^7hIyplGnXNxE$~V` zX_f#Q7zeejQU;pCHfvdRiwzIcLeo_EFHZ#f!95TywPN((?;R6t9zZrE-f#tyd&{mO z;5={He7FZk@DQ<`VuC)kMkI^1<81ABnRxby`>FxyE;@)c;YY>OT^SgSSJ9xMm`&Op zDNNUEOga(?%ECGeM_@HPwevf~$;IS=^kPA16tyjFl>kOb`-HN+$540Vy!yMa8&@l= zOB~cqegt4k#QoG@qRg)ntYNX^9I>GGY&c(=i)l1n_oNbzLFj|#xSK{X1j{JLIYBtZ zo1cOcERR!hSdil<3Ti)1>NIZg>iaQ5L$b{oh?c#D`MUKw~CaAuuC(7TXnqPTKUk zw>!{)5hFVERMXc4JuOlhjM77pfO^Q&rvuUPj@&Ae#p`i(a0F?~oJwkU6zrPQm#Jh> zO^a0AZV+`pwdT6%pZ|bvy5Y3#)D8#{ZCgB*BAv+~4<$5BySh{S!rQ10F_^?kKe-c_q++)se7D2_e|W=~0dHJO(!er{JUuVgc*fNJCicsjX5Z;8p2u1N*CgQ6^O*kLH0$ z03HiO0`OQM5+3gak)RrK2a!NLN&zy7{_5b000?oYb*fu@9Ga6!uV5~aNn;*BCVBkv z1TrZqz7fbI(=*$FOltMC^aBHZkycr!n1ZBsiexvZxCS@^%_>`mA2oKitVlfTF|%|8t4PL{O7Q~0-i+Z zUI$gVI&GIo#yQUU$0Sxc`V4&*Ajxs3Oom(Wi;9A4nCh^AKwrou&bp!kCdFG;l)#D? z2+%`pB~Sn~XlHPK0|lprUczFJSUXwpi_aAmfDN{-q3wYXFLG{Qj8Gj@B$V+Cxu-$3 zDMp`Td9*t2gitp98`A!Jt?Cy0zazgA%1lSeuR!p<=epw*mCooWR=TghlPN#y?3Drr zhO{LRln|z-7NCSs%RyuaWm{0x%=%zyfDY#kkcAJ>S|Iq$nfF^wQH@^qeh{+XQNI#) zVQ2?*h;@(~fF8+v`aVVt1Ve)1E)Wb;WE&6+B(env1`C8%8`!v@(Y~j7$X|U0y|m(m zA3BBQ;p5L-z$daF0cJ?(46XX6T8-Mcy=rWl61K+vUt5c4 z30rg{z_!zE^5gp1L8}PrW46Zd`#MT|TGa4VDm5MmKo#d4_$XRv0yt3zPBOtXdg@tk z=j3w9mY}h3iZ^N>^xV3|BR^;Zcz5Y(8pjcA+pR>0NeFYe;!+)2KkmqM;}<7a_hw*7 zNqfdbiZ>e3!!zzEsc$l;A8Zk3p-q8m*(3$_HFTo?*21v>2q+!dq)L#}OanDUof8NO z!qB?88seBOs3AO3^yg{_^DIz9nBE985=^RVgzhog4Gmb^t}25y^7OO1#c$tHEx>2+ zen_=oGX~v$Qp<$0ea4(l(5V`80YF^+APASh_;QTERRYVenE<~;5GQb@0J@2Dda1NY zu!v{pK70mSI^~ikw#I7nO=6$9g8j1eSFT_k@a4gx-B8wmIgm>Ue6!#(*i=Tpf_g2& z5n6@>*`z1Y{rDo+?#b1Sf682Gn1PvnMWAs>xBCD#CWx9NsRPYA>G~j);j$anCWC1t zPvWC1wXb_J!j!uscR~YpR*t_svqK)Q18KoC*2>4DI2jmffJinN!Z{8&aYS$yOo)cq z#lZQ>n7`&f`Q4BOyR;cRDgL@DV$ZIk5nX1{8!Gg!x4&p2@de= z>JKM@O1q)t0;sgxBY}$FINS5!!aRpnU zK481agRtDmgzvfWrU3)fc-1d(JA~MYTz?NX{YPJ?Q(S$PLJ`SYGP&za)Z7!dI@$0W zUS2H7#w6@J5Eu+4Jaz=zO^i>7&@S=!!{i4-*+Vh+BTn%w7rh2_k_aaUKCwKjAQwYN zZ3m^q$V?PaP9#(?AdA-GwSX*YLSiC8bB5hi3sz6Z?VUhVcDwT-!rWk^wWg~GPzH|-c17GBS= zH%oPJ9aD?nZ+o@307K7@JPj(TotDn!u+DcxA=t%QPJIt}mP&IF$`0yr^T2>AHUc6e z#Ab6+*!jn7l0feo14bM*1j9or(}Cabr~!5tzxSZ` z%ZR<_i^uYne9sX>AfSpL)&T((-4y@?l$zq;*!Xlt7HgesB9&vo*CnJ5ymq;iino6a zEL6rxbONYNuxN5RSVnBVP98`-GMo}h1-zm7caArvBIh(zdKA!5QPpJEj*-R0x@&UK z^4(>+{x_8w5@!ZHRD``Y@KD2uIS+MqlzWUA{YK6B7#KJ%TXD2-ggd;d@?PWh9NUnX9#b{hqT^pmqj!>Yy92QzkEKYO^u)uh%8glpL#ZJc*l&lw;q;gn8raRRoW-F))F#j#&s2YfA-~R#a zMU{ znwmdq(YwWbHL6uw?r^L>mDC|lu_9ZkiqD7boCJ%-AzRlbaepexu*h=4ImvF%=Nt;kV5xl1CwK8?^g*LR?o zJSGnUz0@B|2YSgAw`Wc-xm$sE0bGtCu6K!(qsecDGU@~Gq~#mHw#@nrpqJ={VL&gD zj`#qIn6kG86p`_h3iQ$$J_qO}N(0bKD%EzNmx@!b0KKH)whM&7z^(Nl1Y+CuLF=Hd zPXVpNFpXPr`qdg46l9ej-pMlN3M^eBU{W7%l-6 z;W4@2i812PP}#+5nh%n|?xVgsa2tly#Uk1J$6TN-PFgWF`3?8#!V567u3`3q^kpG# z(V%x|jJDCAR5l@>AVBf3H+Qah@TlZmGf`7^_%e0*4H6cpr$bbGlt0(u_dqXbM zYxBxsUgIH@^Xs!X#qEhB>t|!xz@eaNKr{8do!e;2wFjChJsQKjCQ(Q3-&qw1h1i+{ zn261rgNZg4`4hN)D$Y-^P>O2QnM;Dplr7+R;xjmh`>mwrn>b&=&_hQ9HA?DxFu4G9 z)cFC{pra03r$5>)hWAlc3S}}8_cD&~WGzSm;-mW@1(35h0ZgExf=~ZUU(tBL_w)~5 z)9@#~Lhzc1Ar<7nf%jan7dx;i1$%y%`QjPXo^qt}&Hw)|>9W{cCBdE&l4=4MNsru^ z{qV7qLv)(ysonfzw!NFl8N_nH+fR?*&blR3Y9^nUZEV(FJo`R)hoH?|aPb=YJ{8@$ zMDQS2aODa5Nt9rG$87BM$)~+3wfHml3gU0~CfD@-{cpA%-^$$|`)YjF+DfwO_w7;{lPtMHVi7=;v7na^?tE!#$Ib(yFq{)-g*$>H<9LwXJ!0%k}a#F+# z4&@K`L(r^S^b1YG3oY<2lA%n~ZKc-rT;0(SJyL8h#nE zV%uvuTktnJ5oNYuwz{S1Vh&|(3tcrZz3oDQXL;?$_p|L|t2xiNa=)j{xbc4IdH5C?kAc%u!Aimj{Z-1_4z=l z#rXJ!o7q0e@awcE`@5I=xlhc~aNp5fzQVfoy!p=UD;>30XQ4)#aw31Ygef4dg)m0@;O^ zhKF46hZdKg!l$lyGV6BY=b&QRD^2gT3zN0Qt=-f&fD~K{?sS_U2$L( z>uor*-n%fW>frR(-Q)XKrd8+sOzqvKKE>mHCjQMUJ-7Ab!R(vm24@YAYn(Zf zshO}Bx6|-&qcYKiV5wvTXN&gXCa*jzBYpF`IWW~nZv3osnDilTcyGX^aUJZnC#CCC zYfoL+(70%zdV9?$gZlZ2&ge_4KYYzmk3939Xg}%0sNSx$NFz#ACGbZ_9NH!c+n=^!~Ed{P&f_8ODq! zuc`J`E!OM@+w0FN-~9R}q2pNHv#1*hTl|k*G+`A6P`{l!JoJ+<9KO8Cqq9pLY0SCW z-}{Y5Shzf}oR($6 zjmvyPJ7d@HM#OM;uZYT!8FqH8IQ|3wdX0zCnagfpjlD+7Ts;?A1z@*t4&K>(qs?H? z`DWtll%742_fEchKMTO0!vBd6*Ez1! z5vSeo$@eOAgL9LA1oMMd;Ss0L3%iw_O2ZXopX-YvZ)U1hFaKkD*H_PTMlbpK^69!W zff_$;RQK17v=eurPzO7H+*fyjei8Ki>*I( z+_s$hv&g6-1Bt`FGg6*08{gt{Rcp=dqa!DV@QJu*KJkV{Kh6__U2z&O>xLgMbagi7 zeH*%JcR%tHW(yuNd2|>%Vs`@<({i+#TczH9K+7ltfohounc$ypkI44V!OBx3Pf=sI za*D<6ZOc*?Y}iK~F5Mv3bzV34YZKmfHLYgD3(fAH6NYWZjd}_3t*7J0_&Eok6bN!| z5v_BZWHl#KR;%ycQ@3K|h2Yn!K+~cVwk~^g;#-F)Z+o{It5D7VSkwzH#MD3i5^bPX zdd_-%=J4m8O7|s@m4KXhF$G^YPl71>8ee}nH`W9_ZL#_B8>&_Kku1}O?TRIvPjj0# zJ|BBMH)>eW~^$QW}xmH<&TMsQ_`o%5ixL(JLwD6Y~3-s-|-~^=!TB z@z)_98(P)lGzMQ(@3osZe6eG*u75Fy)H*o4uKs#f#PW&oWu8c`ZwT|Le$1M^b&GEj zo>mVxoE<)Z?LO(tixVIYTNhta{AS`!WCo>w^wThx8xIee82&!{I$2ve&L|-083#Q5 zzNEi4;~o8}cx$T5vHTU5rmHW>ra~@dRUCb`3e(MP{cG@i9>UTo3bBGRQScX35{ zS%U#LSe|OSKK)^WHXOhr)Z<1-^k*^r*zfwN)uP#-amLM zbAG3J_eZCgHBnn$KRpzFGIvB@U>0UHczMaV>l;Ho<w>#Hd*7-pd_nETP+9E;g&w@Aa(uk4HxYZU+C=-{C$ivO)exQ1*Km(HKgvihtZ z)4RVd+Eb`Mju-q|($+5+O*Z(uY_zcBySj6F>vb?{Q|qmAbJ(5UYp&01cUoM(bwRA{ zG1|uBJLKVy*vciIjIUqt#>y}aD}$rw&GHP(Zf>nWHbu2Pg>BoD-o!?n$uwOvzJt>3 zHDg(yrV}rTk@(~9FUmhS*yd9Ea7lH~(BCEYC!5o+!Owhl|KW=E;UA&it>jc^Fwg4D zFc-q-N{(*`Uhe-`Cp=xhS%*2aKtHOv`IYNaO~cuimF2Dq)1oli$WNye_ft#$yV~q@ zOt+XpB_wq=E%nan>{xUp&?&`}iq!r!#)K@^#QXF3uIumDOqbj@xJ;)Hb)6hrO=D7A zw=K5oemb8Ex*thROr5zJ@GL*9O{jmZk?U3x|l!vaE<>MbYBZuWN>%WHIoiR^? zmY0n&Q||2HJ*XLWE~!O*7ISa!EQ#8@ZQb2A%hS94yC3!5zOeh2!7ZA%I`f<75_jhi1K`A*Leek)qKBPMPmgI8*Q_I7gA)+lOy zSpQ*LKi8%FzG6Y;0kkXv_qh4T)Jd6^)yun~8tlut(|K>SHlB2v2e%Gc9=Or%ot&)v zL4l#6MW~NUIWuY-K((70!(WCA^r7$14r*@f#vax`F{{@0%W2ee4Pv-Aw)e5oz~|u+ zKit3Efn-_MfBUgFgf>QY!?$gE+cZ}OVy_v?qO*pA-%JPeDs24@l*Q|udS3dXx7-9> zIzuRXWxFJ$i6qqV!jbzngq(#OUy5;AhQ^Yr8D+96Igu#aqu3jazjsGP!JNc;pN$QMWzb0Jp4=RmqSh*N|ZUbgnQ&Owjs6Ooj zrFLdbG5$R2Qy9`-#dJd|>)dF-a! z-_K2H7HKD`>hJd+w|NqoxiZ&a_+lPL-lkJ1E9DW3QWr;`7@l#nW(*t}UeLFn%%F9Q znibo$FfDtB<``>$|KM$;^&$XJIh?t+tqDiaI95t5BkmC%Us#im4&C}vL@ggE;&%IZR zZrcb7I!|tzA070@l7-)E`YGp3D|(Ez#?l!3X9@0W_DYkJBeVl`dtBzJNt)}-gz37E z>BGg=B?IgB2kZTPW&e?V0e#3oZd!kXzGmvoNvGnJi)X(Ze9C4-Y+&kT+GEw7+%(UHiGkH#Yj z%{FdK&+Quwax@0758uzb)P}r!V9of>)@ycNjkgM<26uL6rhghj@fY&9jaTbFxycPo zaj!U^=vR7kKao386N!f7qL=)ze|92`tTS^KcKEDnobCLw ztnAE{jSugaT-MiLEtuZwN^~4|ZcS=?%JBGbw&h4vZwzvG=G}7>`{$Gzu4VoJKxbvE$B$8}v8y zSgB&MeL(Y9L$AL0_Bh`-CY@s@$R`d|t%&CYTakFJ%tlo!3}d2_E1+*3k^d|4#t72uSzg>|z}C z6epGAUkSbQ7T*(VA3D-mszanNgXg?>v>_Z2q2L4_*(_z^D-NVy^M z=w-P=Yi^?h(q>XV#xNKWE3p=n+AZG8wrQvuW z1|82}&N*yM7V2EWVHvnE2fO9tuM(_!jp;QQ@&Ri;;<3+|*@pXiwmlNLNe!PeaMjK+O37%Vo(vEl_5UYYP@eLNV;F%woDt@h8q@5&g=-1ghGXLAz za(HqOe*Htp8!fD>poSVAP)A2?+-iVrM!0h}J~zRe7HDLLF>bit6Qg}mDF_!tU{a)z zmxxgpaKu&Yn2N$b z8RNB4MHl_`h3hQD3|s8KQmAi-4QuhP6OM33Z&!Tli8|hx6NK_7Fy<8A_z#BVVO@z( z=Q&=j#NS`OEzAc>^u($@_&yE`5^z8o8lMCG=|z;jj;C(ow=4`Ozlf5g)0e zzZUk;!8bFoa5fG%#efBvWrMR^aI-5G{n=Sd4@~ex>+P8CgYmm@{9a7khjJ14-~uM4 zq2vwB$inx&y{3_mc89R)2xfl98Q(DX2M%o))}@S%I^y-vVh+}sVWcHKFGAh7*!3gs zX+q7f*dqIIo=7J(%+MB&^>DBe7R^RWGaPT#+44mgXN`ZZ!9sgHw*fD1#`cTYH3!r4 zguD+zdB;V7_tLh{IT~qskI&Z z4@1vre0l~|<8f*Vs^7u4`Ka(3|E?C|KjEGx^pGFvD$++0_b6fCStvIT7cRmpD{$a? z%-x8M&e&@^uJJ+r-KcZ~B_eTk6jsMz^I0s+M4LjP&OCl(KZ}=-$ln~Xi|>fE77+K`+UN~&8XOpl0Cca{r~otLW3dr zX&BBLhhHY3#Uz|Jt+UDom~DYB3vugxw0n$BW%ximHdN$i8C*ORCy&J7(fDQzPFKcv zs@S4|cKWDlhH2)wdI4rE!p7xj=70;FG07Dty5TEN9PEdh2l3z`lszi+8HI;VVC_jf zb_PG4L-$0~zlIyG`;6orh4PvIJ}vN^RHmn z3{1<$j9lFP45eS-;97jzj1{f;d(tk8Eboq0eQ^Cy95ez?jK-nU&_f4z842+g_$?NX z=|$cWX+Gk4p2$baxN$0)>SF(ScxnkM+2HCm7`XvsH{<9nn6MSQY)A2O?DH0{h`)c> zG3WhIP6qcXqSh#^AB$TwuvH6f4Y1l$sAGpw4(Je%q36-;GWr+c%15~9m5^VDPxPW+ zc05;Oyk&+_R`|dgeOKV+U>tZ9BO}l&4)t^tYeY5}VAfomVTp|kaoc)~cSIFu9JE!4 z+k+4Grp#$S%ObR+is6Q#XS!58!QV%$NLIEJnT z_~aP|RATuXA+8Fm>(Q+NeOj>Vcbw2wp}ix&8+!J~BpGy-$96URW`gcZany2rYKMUi zSm=e`JA`~cOo>2iF$)Rtj_+l>tBCDds5}iPTH^wlSu!G@$zqRXxMe4n2BX$dG#|5B zLF8gpj5fvCWw>J%c3F#QftYa+kB8&Wlei)QKNX<)7qk|07$@>}Z@edi`Xf zc)$&Fw&5>-RFA>O33%-yI^4q{BOJ{-?)wA`3&5Gdm>Y{8=kR9+PP&7~zGIPu$C8fc z(HA%C;8%StnTbcNFu)EE2BK*Z-YLcI(i_)`oU4fswQ;^VmRVzpJ+6qs#0#jJf!lJ> zzYv!^$7`>#wizvy#5Z-!qYBDv;amgsnt{4DDB*#|JFs;pD*0p2aIA^K6Q}TOJQknF zfiJN66&fl=dW%#Tjeo1+er>!x0~_|BY9!`Zrtj{UXKPemiPs(R{w7S=f{U3P!JaeKf001A02m}BC000301^_}s0stET0{{R300000006RGhAjX9 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 d80543e90868f6f2ed74103ef57e5e834b7fac7a..8522a7f72b3e0e6bc7bd2be78785617b4948b92e 100644 GIT binary patch delta 16668 zcmX_HRa9JE&&Db4P~3HJDQ=}uq&N(NI}~@9DQ?9I0~9UA3Jfy1Q{3HM2Pj(H;qUuh z{1>^&Icx2`vy(hI$#dZ2grgIL)l33lH>{H!)f6ZjY2?DFdTE~%?gjnK6jmr|YfKEr z4d-d`!t*N`QOI|=EZxYcAlI0f@NgE5d|y<0PszjM893bRM??}K%H#8`^!ukrdzSrB zkEsB%;?B1gSA3r+Ex*#Ek-c|a7u`4hx{kVkkt^CX(PjDdRKOSaq##lB%JYyHjGD-& zW_w#H65zRhf7OLc%6QWi3^P6Me0`j`@4U-9l?I1wa(?K4-ITVwUpUkLGe`@NiqgKG zy3?RB+~s3d$R0?2tNkfgK2jgUuq(9R#f6er=BQ$4YsMtkp_w!BtfD9WGSDmdpo@Ri z3>29~{LMBbxjGG+iWCE4N9DfMVO^+eB5!4(%TP3gd|TksCvMa|5FHQ; zh-|i>t7dKop&j});6IfU5M(t&ky1Z zEh^-DllHBjDG%vV$Kh9aQIsYl(En}eyJbsjg+2m+0EZf?c@~61K|`S2ktLnMqqoe0 zS?H8y+y=8tw5{o?s|@TE&+RigQT=p1P_?G-=sMP_-x6yz)_^Ne+hAG~$9?+rcFbqTAsprna zw%givE04y(Qr=&B#8g?TK`sMZWmC|e4u`Q5A7a$smWG=d$3m6DA$gTh!Q z1l@z|yI9%GP`H&h#WOeZHD7i8vS~@vod*Kno*N7S0{O9eow$S#pPg!I3c8`kzl!}$ zD)m3jTv=d3mlZexX|WzMoFw~677u30!=Vr;lEa@vu?ovq5%9-EgR$Ya3D4j80t+>x zaLcFc6F~7m@MV9tWTr@#xjt^ww2Np4v0U!y$wlDF+o0TMW2h=Pz>}4uzRYdNmfd>FyYBV!U{Df|8|y>`jYL_ zFwFed&D}32DU4FzsB|jqtrB3L;`xMN#d!0pCdz>l*}VZPCY6c4j27Z6 zUgcx3Nu;K||IhldVb=)S-65Ms)$1TRdxtbKoK|`uA@=ahQXzdk4In1KsrLS7{;~^| z^IG^E9ZI@K)rpn2*Bp4!k}HbnB2%B8erl#O<$FnwT2+1~jdl}Mpe|?YWJOx$a{ZH> zEI)o{Y~CiSbeFaw?brU|g_^@z;a?kU*Jm*%iR$vUpd|_1#Ge6ombK@j7OQ`;$EU|V-)!}Sec!Hzhm$6 zaU|v|C+1JDv&Xg+-nLlX-(9@NIQr8W{O-2WczB}R@#?)q@W|JvKnR)qcQW~9w)0^hWeE?Ki~M*`fg(Jiq<6y( z+Hy{D0`=4Jx6>7Ur#De>(=2Dx68O~gm{Y2~VbKflT>rsQYP~9dVSm4(>RFQF0(~Zl z6(>hc&b)ROFV#7qcm-DTqA*2|s*B0t@T5@aq~$0!6suqmOmJVR=Q#TTOgq(eZvWG# z|Hw@2Z4y#DO^rbBw#x9B{7_Acp0_1nGuk-}7E#Dh7b zqSqH6jqiCtJj$*D+t8lNf^TNuKMq{X`4^~48Eu%(`G1uB`jKkJ|6B{dwvEo_Ua%uD zhjr4TztOBg_zalP#E4s-0fuw>Z@sIin#jLh?GdlRB6AMOm>c~t^r_H6Gf8UE>-tKx z-Y17F*+b7c4?(&y+c0FsL(qt#cgodHOH-PigOi?^gr2r)YY%-t8&9g~XzM<@#%Asp zC7i6{EO0E#?{LaiPDq%!=uzVTJ^r_ZZ;u*)^sBG_J4a^Yq}`zVe# zVjjX>SLsY1w0aYu+hF#rEA0Z~0xuKfs(>@Jferehzt9t5t%E zK_fMhEQG!iepAU61wT;u_9@#ZJn6nmQBUx@_=kT&2=fpMguzB(eWh`9(Z}+kPker) zaT0IJLcNXxh({JKFVNSI+85i+X-DUp0ouZN~Ivhp<2&0wD z8YND#O!ln*HaBfr3f%F#=gNZnI4*A#*ezMh-8{5YuKulD7&(K>?6>{3My7~$c$bi8 zj`MrUunNqYLJ{s{j-zG%ugVw>rv~D*D0G&qwP=V?%*m%@$q_kce!@9gG9krh8dX(*B^p&{mDm4=! z+fCz;F=Fk15Z@mQAr1Qgt$WK9zb35DJ}#Q~SjQsJyvyJ*YHpxR+H^=@CZpZBC>f-epk^A^AJG`4Ns2Zzv7YVqd6eZ}Xg?w$QA+#=Su-!Jrz7Vzi6xlwz6s<>4ARLmY zt4+cURMowbQlW<$gp>GOY*W7G;XI8kv}Z4T)is`%7q*wRJ_dq!=G5rn->_Kp?i+vi&zJ{@U+E$xiylf5Fs!A_pkN1Tj$ zlkk)BgW%^oOhj186a*^9*JE0mE?H=6Spapw@w5O?&yF=r^iO~e-2LY;u z9s;TT&YoJ41O@bL86tDG4}6P$F_nXWdnj9VlN$C0gm3(lKZFe@$?>CmSqq|u68BdW|5VPREq#MQ4SmEM46za$720z4C6F zwC}1mW7PG_^yj+Cu~UO5I+d?dy|jVv0g@gb8CfbKKTJGMg!8CWbGc6USCwi4;Ug56 zy=(P6uEZt=;Dqt5jmO*V^;c_9S`H5j;_beoAo_;_j0kfI?)=};Smi`6a33)$hj_o) zjju`a#=4bf&~i*9Ufxv@LeY78>t%z4Bfu`MnVvH!bjHy^(oHriaFU7_Ed307UO*^LV{)oC0go8v+;yLg#` z!;5}204O>#jz;iNEaiTO2<|u(ckR3uhYqexu`9&qguB5zuZi%c=40BWe+@T0%$6m+ zr(yqKYjoe1BR8MM+IjHU2W046RDP38L^)Gh#Z6jYch0Yj4*JEbM8l+2C=L~0hnxli z)*+}M@Z{)pE;<$zCA~6{&`7zOXz2Kz+@l+OXGEji^|_{!aVu3Ab|$k*)%BR9i3|Ur ze10p|*`CtS?A5cc)QjdO^- zam=`EedPyqOjWiy!yMixzr}lg@f?j2h&=&gHu(1Dl^m5{^n8Be%w05rKNfvKO47_~yj9&3^umUzUBRhY)jqwH@WPrw^{Z zjK4oF31R}z2u-XX!P|yDp1Pqp;?|F@6YZH>9}+#hRQ62-E#ne4yxM0zu6RGu)4TBN zfWc+%Mk!WbaK(PePYQz;%w^~JhA+EhZBIwZtr9K zo!AS}dqC|b1hr)Bjdw+58?E-FwpjLEZq_-u1=)R;8*6v{`%JqCy-djCW=q{j10kdL z8El2gus+o~t9s7l{G;q+hV8TF*xFPhR5U7yQzRY@nH@L&c-t}w;}@@avM8GsczPM- zcnC%@59w|Bh96z=KjXLx~R`9x;^XUK?x;ChVpCoop_ z{P-iAf|l76Q>iR<_{W}_GwhWnvYSi&cNa&L+}tx@TLRqWJ6d4IWrBwog>c7c@g+X?O5lK6PK3mx&Tn`NK zkTdT8dL-g!oHL{z2}QNMnCbTND)#PLKw9pXA9l}8dO{3$sp zd&{gqfx|PaF^3$R_HVilWj7DIwVSuM0vlH*Y}&L|(@?z{DK_1Jh}({Pds?0PL3CX% zCi6e$NdA~4`Rf5{&^aX?Mv_0w7C0BQOLusecFsj1l#e3Lg-WI=NQsg1i0r&0adQrQ zq|TgKEtOZce&k5Ei3ye9TFp)za z|GzRpKzKz%crS;9P>FGpBbjo^MmS)3S|H*}ieJ@Ok2r_8G&e^g+Lhwiwjtfb@pdI1 zrG_&=;^LFE;j$a0l7jTh8lB9+?Z&I(Q212;Tg%0>&e!CN4O@YMH~0)PY_G^yTi*y7 z!#&l0wX%Sl-ojLy1GAQE`x$2COke_Ppk+cj-AM)cQXkLRl7cNVta zyxfokJ2m(Jlktz7O=`-AR)rkZc?&4zdYasxZD17-_trtpAeG?zm|*gl z6&h#aF>m*{UMI3fKWe(wT(Qjs+Rg<0kvuK*UwgSDWk-?=<(^@utESf}vlJr5gl0v9 zCY8zcqe=dxhshQGR#F<*DeA{E&G`qH-VaIxLnHaAflS=LHdLWva|R!cREgGtdUJAo zZ;ie0V#e+b@8kPRuVzIcxfhXI!>kSFJeU^hJJjzc*R>c7$K|XW*aX(nzK?sBQ=!|& zBi#2!Phir<#_d|-JuV(_eqcq?c+xMDTh=cSpvc}Joaq=t2qwx}`I9}LVxh=Q=c1M4 zufY|b4%oV04tN&&HotGs*f7;ne;*T(tbQD0n^LXFjhK_Np-2y|H*)@I`gT9Hj_qp$ zwna(Cuq*-H7V%S0X)c>gXvXjhU!6SoIob7|fm>dvGuKm9>HM4Uu0H%+d0WOOU{_!z z2rc?ON-Sjrk{x~Tz+T4YNyUCv9c-o^HUp=vOoqoOzn9`enn9N1BD^SjHX(VVGEt() zWz1Qa#!VXAsrm%29~Qf+)8q3P{!LggRY4ALX_oPs=o{rjCnf)5KR%2Bzz48G>43(- zNXnD@2wfqM(HQx(=(%i&>?KK=(z*Tb;maXLV6J*DC?3R~&j(-@gdNrrcDwbKWstjD z$ei1EDRaaHR%i6u|F+-$d|qt7mx>ggp@qrE{N)X-5)X_|&-^8k5&_{2;(JAQL?%V{ zUVOP|b_GSg$T#fsh8a40T~Tkfza&|x9eaSkDEny+59`M$>x%+B29u|#N-Cg)ScT=N z02d12I_Nv178;sN{atAK?D!`=lYfUmj_b7_x6e!!(gXGtW(b+Vo93e-J1q;}(p*NA zC0e9y;9$rAwC!S6W!knV+(S{XuQ2WE3&eex9yrYYQBoNhXdtp1);CHYjv`upksj#Q zHL6Y}a7lElLz@5+~5RKwp*Z%iT(mDKGAm{7B56?nAT z=spi<4C535ll6SLs~2>`bTf9<<5fs&|TA-6NO;Jm5Pfe`d<4L zbpr^hn0aO{`0tMsio_~umF#YlvGm+S2FC{V$K5FV-K8{Xr@ugmXlmWN)kO<4#tq4T ztDIk`Dk5q<8xCCBJoYELa6clhpVdfnPPE@V^HcA$RqhLS@s^|qS%`rZ*x~n{7Mi#A zV1~LozK+^yJ=KodL-7FiCy)iP11XcDxM56C)Z)`$k7SW|DyOcFX{&G-_C)uR7o4gi znUi>)Ad}mV_74hboip{p+Go#7%)Y52*F}7lG1~otZo?Ty!&URTb%zJWqtoaDnI$z= zHhw*$7uff`il&D8Q${cvKouv=&iSxC_)z+^mVZ_47alZ2(ANL=u)joAWrpBc>BO*F z@cYZ$SS2YDGwWueXb4V~ZblM(;?b3e+fy=;3e$wh1HX06ZbABDuBxnwJSIUBUmkj-JR$&Rj& zlHR{|v2Bz;l*NwILP$xjeCjKRd`F{=gXMbYN?*xbu(O7P1x8WLPgpNsu5@xERX!wI zw`?yIqniiOim!}Sa*C7!oIvGa?HQu*--IxOnu}!%po!u1+!x3#(M6ZV^yQBrpqgLZ zse(|95ChKp#@5{;7c;}Px~*T6xs0sjLh*n*dUCZso|Nq zQTC5%ih3SEt*z^KWlAUIh=BbfS!i~8(=qJ)`*gX-#Iflm ziRWRE54j7wsoOwDK%@F|yf}GJD90MX@pe<$>^sPv)Zi`JrjrN1$=07Ul{*&8>GRN< zdEA9p1MwW9RFdo;?H)3Rnhb+O8vAq`rX^gqOSA4~mWin$w$n^@s-{+9`6>={k z0*;zKK%qN=*!*VJ@$U|Y%Q6J9r)}X@0*;2CC>p3?HJ=}}gxF1RKz50i0etKx*ZO|t zBAK`vpE5gUt>Mj{pT%T^Jci*#;e2U*0+o4WBk*AI)_d@Db3ppEByCh;ryF;#v-`pM zmi_dXL1XCWShYlm`|`g;VCN5?j)v$(G6GG=bBV5A0V+uG)6q_e6?RLRs|EfRlDT!M zlP0;PjN`dx>81qdsBQ~{K3Ln;GAjgE?yS?dwVzA4>e$72>I^!V=C3QpwVDwi;oY7e3Z?y1v+Ncler?PHzuN(oCc8{)(4l1Ui5Gx zLzQBR8s|j(NUv^UX zlE|el&P8K4I)j`IcWD}`mupT-e(&7~uH6bV`kRO*yVQSfsIyi4BYsdhU!OQzS9p#0 zRKtbyeSzWO#8*lZa-=uz9HaryYVIBjPV!?a5qUZ0SM-6%FuXM0V|877!c9MA8V8U8 z{59;96CoG>NP$Nnc-X1#mkzZJ+kZl=1C|{~~v2KODg!5JZwUbwVALQ zI!SYDuIT`?m1b_q8lLRhnP}VJ1<9LJg%d>MqB*{SmQ4k*!Lg2dJ_W_A$ zP3+M!8~AUI?UlI2NkFGrM10`6jInD;8*5ejf>0GJR0HtuQHC~+yal-43ZC;5+NOP6 zB?cv?JJ*wUtj1IYO6CmU>>sz|`WUVFKZAb-5dvvK-HFN>+(#$PBi0%IlYtW_psYR9 zOIft&9@9_7QTu*n2PX_bMZMD?vJ<^HFG6L9?#CpAC(7$?_pOg>)$7ze|^Pii38fMQsAWX~vX znuPP1$zX8#2QnGI_ht8z*GzzjI9@Aq=I!AE$Igim(#Tr?2B#&Iv-&VW!<*!&7P%|{ zW0E44nfWcDAaQb7Brc#zh_9RB0tJ`_3FY|BVIllkCIACCA)R~1k=G7gmhJJeC6Qt8 z=J=u|!9HVY{pC=X7heSY+rD;sGvq**K>weSu)lwhij!DcFoz>b{*xdY55j*pRm*>X z1M6Z)1J&K=O*@+m>Eph4y`=jk{9}nS8?%?wi$dB`A7~}Oug4Y_R z26!7I`!DQj2)2S73EvMx%SuZ#C14}!7KDp6x6>-@TEXn{^POub7@e4KCf=F$PgyeG zZBmNmtT9xhz{3}8m1Dyps6YMWN8Ls>dz1IkN$cu#Gbg8kcJokNK*7OJHv&JY+MeO> zn-~Dw(Ea;S-0~mDbLBte2HA6i-3hCrNJLk?T*l;z#&#KG2Xeuzt(Y`3D=KFsMqa;( z@iZJ=p^D1>S>^#SaQa^K~NLC{J(-CG>pIusn zUiD^#VK2j-IN9YlX!?C(&xlOeT0m z`Bc^K1|OSWajk=?JHzDzQ4>}_jyld~4coMLJ=xVYqL9^6U8XcypZjdjRY6+X2P7R7 zSJyLzb&II41Hmb?LhrbYJuGz!paJvEo~E&2&%SP#hoSwBJHR4pu(v+Mcu`|10_PrS zO!Kts?#hUgvV*qS<)oOhZr?X#zAZe`=w~j}Lm&rhvveHwuV2kWL(ck3|JTG~&M_e?0k#Jq(mRtfB@3@MZ76myc=Zk)kk2JF+(R< z4^d5 z4bx~i6i5&=G{GyByP_^g*@V8tO1n$BwZAl*NO4PD4ec7XuE8werl!RjfG5WJ%HY6 z9DU0}63U^`s19SE_qOyQCqAz1#Ahj)Kt@Ellnx*T_$JptX&cZ_7Ay+U9O14ze&3Pu zNPfE7;bpWZ-aHm|Kshlc*E1&9zs%C&9kC=4oeLW+PFqSnGhh*$pNJcIOV>786R$h^ zXMQ|`wj$BbH}D-auY%KwOhv?R&z9=T*|CVE*RyjP2A^Y#;g`*GKxUmW>GB_BRy^+) zOsn*-;1=|bI0dh{116#8*gM~kY`ZVwC<`Zo9_>nWG{iI}`YkRa5{l@q2CZLE%GIje zQsdveBM;>6kwk6_!ImKgi?B-pb~_lba~^ZTI^B=|yj01aBDRN6tEry6Ea^F%<9aL| z&4@{#ZU+M#9G;?Exa3nm2x9|;<8fi`Gi(ma;A*G@caMDPcRXysS-S$X?KTs%shYU9 zZqSDePB zk?)&7+_e1`Q+}lP^cb{27(DPgnTF&Nn zgI&HkRhFmAL@n~bY*5v-v!%MP$z#Zw3sjEKBho<1bS$`8#iR-j794Wwd=*$jy)N52 zd+%+_auSU-N;4V097*xTYe%^vdK+K`iccxjXGCmG8Pul-3kQc&W1IfoJAR??X0ucw z@sPb0AU4^y<21yBiEN}+0N=9a))gdo?+u70kz78nw40NS!&>)gJpagY_wi)3+>%5z z!O^&c2l*CQsCDU~d_0|0-1jGKr@i6&Jhe!U>yL2_O$jWLzi_o1*D5n=e40yh{c$LmBKSJDu9q*2v? z_Cj({0>_~x`4V38_oi(7n(CD-wKZ-f$AF>R4ev*n($#abg%tWpnQH5?xXX_SQiC^s zLw2A@PbrF34*{yCTCIoCo*qhb0%i{#=l3?y=rDEfBi7Frq>88)Ax{j%1~}Edht^GTlw_w_(@blagfAv6 z&Jt#Ud+qt4qtS-A@~KG_k&xS_TI~;0fJ}%cBX^j5Dk=NG@FuAp!AP{zC%$>R9K4WU zv|EId&Ptsdu95BMRAdKh>uayfy){82$MuG-k0ZW4l;h{@>k>P4SpkzHB#QyBYLZe< zK@ZK-pOkX@>3)34e(OoWk=;hPAxQV5+%i-RoZKf0xGT)4_AAQau2zCJQ(G^_E7r-# zVF~J-Kx;)!K7~EeQZ(3SEEO-m^XAA*w|F5xr2LNJ5n;fU=i@1?RQA_7xD++R8|8JN z9~ZxS>T^7Pf`Wsci@@PK#MnFN&#RrcKpUN=&zYL4gr`ciHOk1-Fq5;j+;5C_st{>S_p*q|mAhSz!( zmR5T86Cd>!Q2W_8;yGsfRJSIcc&ZI|y&3gD3>9M6c~h?Fi5vgS0QwtO{Ubu?=S=`f zUyDG|q*w_7;3%7An=YKBN$DN|5q5zXdO$LUIZB-5#q+b;9`uYwI0eE%jQ{WBg;gHs zB&XH=NQ6+v^0-8#DXgDTTnI51v8JMeYxfPyci$K3MUKe^ij)>(N^qA!NnalLxnj0+ zqqk+2fmE{>?_BQb7Tj3^$GZ=Bn4tYpI4D23u`6x&xR7@S&DO(&xto+`IuPnJ+0>hp zw5LStUa9f_^Rek^hIV@ku2ND>Lc^_uG%PfHE;8&{Ei%pCP*7=87;4IyOJ2w#YmDV? z*Zu5u!l^*Zz!YCA_>mE0LSYrD3T-A$)`#6d8O-`F2!5p16czHQ>gmsnSI^<_14yPn&I&a@4>&Ju0Pj~?Z)PJTd z3~&;$FlXFl=X7cnvpdi9;I4@_QVY1y+umUZ-eaW#efd&Odj}!dUQ6_7 zZXGq??C+aM7Q%JnDab?JGXZ>b~J zut0gX@=SvQK%y|IkNvO`%V#P2!{8N1Ql{HQH*UN$VaY_mwMW-%W6K)0+nB9WbB%Poy$CB86>Tcyle=_jWW zXO728cB}M58sC~&TGahPFiQ67*Ie;1H1x+?lMtuOFH$dPT#7pyV@)wa@XQa}-Y7~c zw{@Olr5(}PB0Ske^>r*@$y>n_-;nf0(k_+0c2;DdvX~4s)lv&+rZz- z`c(t~>-Ig9FNo!4lhD<2#lF|>O7mzQUlYCz%1PI`Bgq+tQ%Z2OXOD;H6I$A|R_4D; zn=j(l=d@VLq7Of-H1HXn5RwI?_XEgS8WM$wfsS z@b?$rWN)T@9Qg4QyabyZ@q9>-rnoPJ8wF4x=&oMHAe_SuTW;~Ep%lC46n`WRkhJbQ z=OP9~MTJ%|sJ~aJbyxn0zohG^S~PE%Md-6QJ^E>4bk@l9!`|t#5-{6E6Z*9t7sfv? zRs$Ta`)z2;cnCF?_db|`N+=$=Kr@%<@ zjBvma2b$Q$fc%P>7byO1BtDPoL1Sy&EGys=lD97Q0UFxf{y}w@?w@TB;^m<~g=Y)J z7D^y(`}g3;k6Tx{!a!N{c>4lXtM-2_cLlrN;kQM2#YSlIQx9cF;$wy8+Zh-JPP|CV z6z$-(ai^f`T04{4orwZbrk4C+F@MM*y?*f{N=r%j*)aL$G|!W*MlQpcnCeV}q?4hx zvq1Yw6KIyF?gk$b02b`> zyoTp=l$DBR$<7~Qju_LMihP`N2AI@v6c23a^AqzMmJl&`oAF8zaT9YCZPC?y`v25H z&xYsODUiic=Ig!|Ab7V0@-FfR@wO%X`J$UDMaLscB$^id8&l%u`du=sW=C&qKm$7% z$gz3wjb&cT6lI~sM*qACITnIkhcL`B@0H|@<>P!d>X=Es*MD4DyK-R}dCONQ!mko) z9X-yNMlWU`fuAsler7fMIW|YvDv^8)(`xqSA40C$P*m4U)Bn5G+MoLDG3n06eed2m zJ+cCCw{2Kwj_n@a3jX|KIVqIBXFEz5A-Sv<4TroS-2b4

B-dtN%`JRyz&AUGv&uOBc&t^e$l^U()md=LGEEF+K z-#;^nXQM;gzWpoul+ki)yn8uh1~_DXsTjwC{K*ndgTh=v{bvIHi%p&h3IkW*@7(Rb zFEs0e8tn1H{(gj_ABIg<`L()T%gS)goMU;nJ7*SfrBm!q|L?PwpVWn zN2XwIj;*g-KVq66ci^T9urZ;qba|+cH4?!EHNQo+fLj!1Pi}Do1rk^@L`-sz(nz z2|+4PGWqwTgI6p<*-?wvCD?}RiRe^Pwmh{z&i;}Yh6oe38a+m;g@FDnP-(e^;oxTT z?8p7?4k3!VXxkWX7X7k|j7ufry7_xD#N+FU@sgMq7TN5HLCG@2UKG5V@66kLrs>tQ zxv0^&n|hOxJ>8)SF-TGL69DN{{hrNAygpz7%^lDxKI+-NKQt@zimmcVEc4=6L-vs5 z)$ff@wq1pa(Nhzr{FPh_`KAXciW%F8$l z3!Hc(uw%d6aPsTiaDQ4t(UXX216`{JTba z)j>yo5)f1hw^?D*4dkiWrcAbOU0o!&hr*W6eBovM@2~JfL=a zEU4}u9HTojGEa9H;p8zgYWRD!4b>hNXO*OMwV}q~)5r8?OrTx6Ouk2U!GJ}tmxs1$E#nmwXkdmJ^S~QY6k)K}I-S^jV4tX=2jer9jBkaV zl;UpGkYpqulhgh=&G4kt0tZCDx)(&NLG(63yWjJAxGWF(@GN+EugS(%Ou+%AM?1$(G z88E}x@ZG@u*t4|)k9RkkLdG-#J2B*YkdY%tJ|Zzk5EB{}Sew%7c!AElT;i7+HqYo! ze+-V8lbaK`&+>|xV+K{C8|Qv!aqFD997_X0*m;u0dq$PD6(CIejhPro=aWGMl>7ko z{|n8omT1a=QM9F4a8m~2UNCu&(9#F z5##|Z(2@H+9CcjI3^XWBf&e*SiUUxIY|RGxa?OU_Ne0oaNru}<0gD|GwDx>FjCi+x zin@7GOZ+J@+c)wTA`6QcGl&!wRU_Z`ldwg*Ml3vl#tVaAbQgraVe6YWd1cj7%)YIm z0M8D#?-7>6LXzed5k&dB&p$g3m!0VFGKxh>bc4+1q7uu->Gk=$YpRO>5qMA-7d!@^ zH<8LU(P{@g6#GR#D!koqn&Y6~>qv^=KI&#YxkP^*>)yDELodBN+E<=do5;J><`FJz zT=K}0_}p0+nQCEv~vHhL0R;HVrLaxCG09vwWGhi_WEuzFQ&N!5qZ|$ z_vrGCDIWCE7l3{-h*#jiBGkxMlmg;22jd zP$|9d3mEs_&ma?Gzbplkzc&Agd$V@zF!NotF8>1G$=cQP26G8?f}8yJ)52uD@@JN6t$Qm}uO{-7Na@(~U>6E5(lr!D@?nS0poij02*_RxYZoJz#HN^SuC-hSG80E=8hwBk$6YetQW?H%Gho--6zp zlA5Mwi-sGaJmT#k9Zdx-0b;fjkGWFBe(%(!th!!h)HBJh3;m*5n{y<(!6~M#-ge$% zfX9fyfGb-|7wL$`%aLG77uuVpVjdjF+jmwskB`2Om3$v}Pby>D)`WB-pKEL&IU%<2 zf4f@vPNwWxAUN+*Y)W>tpo!)z7V}ld$h7RdpL4wpGi!5gpetN5OH;8` zR|XFh2u0_w8ld?F56Sdut9xmXUzX)^ZZ)VBgp*0M8F1}{_6`>_Lp(+Og_h=&~^k;MOxqT zSltZFJ9C(#XyfBkZLX@Q-CLxR@vnS^$Y4sKswF47y9(IOPl7(c6H;C zhx~Nu!f`6ErG5EPp=T0@hjJk8&2C)4MY~vN&x+Oe6ENFfp;g63_=((LM~Ui7)Z(iQ zV3_m8zEFBZV5s$e8*@$Q!r>m2@+inhDa5Cr=^0fgSKkl|K1up2lw4ZDbEE_viT#DU+7*KAwYe$bu4dPr)$T++hP8uXlH2ld_{~;*f9W=dC-KO` zWI5oBzWv8HaD(1)%;vpUME+7oQhiu3W4dV0y+G2v_V+kGbe=8KoL8>j9?N7DGC*|= zl5L&eGIveGtHh!K!FYR-{AcW!>-YY&1Gf^&+jq+uR43iFm+Sosfg1f;U%L`XQ-9?h z$n1OLR>b*I^Nr)|<`!u(WznPoG-=gWaVtih=+jpD`e^Q1 zj5^_;^BQ*#Pam)qQ5UM9mCe#9=bHiKoUGlJ*Gj#82AXDCM`(8rW-4*Io^WMf-bIUl zET|^b9at5M9opCB5>n^5T9kZnc70{t@1Q@!dnoqmn%K{$4N^Zm+O_#wXm815t(cj`SO+Sm@dUnmYGo?^_Xj1wKh({LnL{O>KvMa-KPP)QR+b-hZ9I^VL2bOc}0GU^!8Kvqd}r_dYsN z1G8kiA!$U&=}ZC`52CHRX`4yoG@xbL4CqQ?Lw$CWV^*Cr2{AgvH} z-hy=)^U&O%lOpdKl_L+SE7l$3mc!BVj%|1XgMZvVT?FSry6 zsYvZ^jinfl=Uh=Y;oac*+ALOG{n}f4<@6WtpDT(F-+R4D@u8#fp>rI4r>O*ihik7& zUIEKfz^VcZ1=^k{(Dtr^ShhTP{h_PedTUfTu_vLb>Q>r-DuZ12Q<~6d%sE^8kG4N; m5fUxpA6Iez2y#}29u&6Whm9;7vl!cEoZ4lOwpgyxET??LwEgbt5|6^ zk)}2N@pq((9y^^A=RGK9sy+NOR{laR2>$TLII?G= z4G>M4;U=OtKPKPHei-Oe`SNlVasfUUU4D6f*QFI2`LZc$aW{Sd?xTD|=p zqvWj^M){D8>l6N`X7Cfh(5y219Bf7QN$jBRX=5sr(yCeU#;&XX#Hkp6sp z7^sM`{z3FkfdvsYQ5CtE#u{~aAG4;qIWW93R3mr_{qR>kQT-$tm;sm_^e(SCR)bGX zHWq1eCK6KX=L{J=3@j*bZujHGnEwkpJ!!3?2*<%(etKGcshdzq5mh`7`E!=j&Jj?M z&xx8C*H54I_F~#9`MrW$n>5UMPR*{iZ|E|d4p^v}PTDO9n7PW9{_l^Zc0ymjKYh#P zz^8^fw$WFBo(8so>N=u4p$*6Uph~L&(Ts;^HIw`iSn^1Xt$T{&^=`7^w*ys*E(^GU ze(wNg6{25+53qGdtk9DirsBrT4qg7}7kQm^PaHKd5jPj=58wJ-u>cLx?F3)kC~h3; zH9x!Ixn&vIGScyOhOTEx`lR?NSHH(d z-{V=@ln;i}K+@2@J>0BjoxRuni6I}|%ft;W z{-h%3RvY7cXUg1fW-D1@%S$RF-V7P!2v8U0*pUy6XV?AQgHO#lXgC-e(TuNkz;r%O zs-+bN*{;H>%1z1Q1^#7;0*M##DCJvLeay5Q2AiV&sS=Bs9Ra;xW$BeaRn`Da{5JJN z!bMl`%eMk68M;tl9_JqXOvQXxE}iK5}IWHoPW+(1@>bHrg)?Jx^ra4tcIMO zG=6z;xLezc8JkmE!|ngn>dc{!h3j)Uh9%%rSV2yWg>j=*-Hp;RNTGucxrvJ%OK) z>j0;9aDBnzZ70#u$A$vkpP7IZX?VRUPp5Y?16vYy>Qx;yJVU?_Sg@~DAi{! zR5T6`pe!j%hs@9R;88-3)OH6N0(;2nlfcUK?nlxScenDKJw7Zu4< zp+r0j3m3L9RRvm#GNxGeLvty1@S) zy;>mWLQ6*#l?tUIpm{`K@Z@{746xnqQY?&R(89`hc zN&L}9FYjO6B~U9dhC?NxGX43%%*oDE?(PAyUO(Hl)4lZgVF|}hMaAH$GIUp^!bwH& z#WmxmVCrvNkUfuSsWoT}-Mlvrho>4->e)B`~&ususb0<)o)vH}%_9C$Gt3 z!D68pL!hz#QD33nc=RL*OwxCk!Vx}6B&A()KIA6f>PR1+1M%cG`M~=dcwT3iUEd#T zdDX9T1vGYh={Y|sYp6JIw7c2?me$o-g3HbojrmTf=KelBc^tKN$Vm8_k)G%0^jI*% zfDiurHp8AvNe30`iwos_ydO? zE1Q?z9W6kDLc2!+W0KOTzx7-0jFU zWgc_FoBmv^Mu~dh%OM-kkK>$}^d6FhAV7oI$j*r5AH{qfwcqTA^`7o3Ht-8Z$K(DZ zy!G8gV}QL+bg)7oKPkJ?$9pEJ5Qbm5blY2OE4SSIx$#Rv>lD(v^{`By&^g~`R1mm)Ebs}jR8e~Jc$ee~Rhh5`uxHGA=eX#{G_DTu+E z^}jm);g z&6n`K1CKkrlPmqlqZ%P}VTa`fKrjyLUy5L)=Dc-%)=TcrRAvEJS z75^Dq4s0kX;JgbB@MbjzW5WO!N(jIvta)pwCYy-_-j5suxzx)X8Snuzkd2VEzdCf8 ze)s0CkwH8g94L~#V%{cBwDBKKyDziXfR(g3*!Lo zo+S;`qE)r7>*7*}Z|Y=zxB=FI1?L(%Sizq|&0 zEJiQ}>OuE@es|ZGwC!IQz&_T+uK$=vc$d>n8O&hoDHR6$oeUv}bDl9Ud3PLX8V}+H zU6KDd9Aj|>GAqUR|AtN%hEr4_x`Lo4rFI|JRQ5|=&9eU^W7EZ?L!&Ya6o zu(EV!12A4WKvO5VlOdnhu}}fhS^~!lI1ZTHnDeZvb@KBa0FE%8F9qX42?SYS27e{g zmj;Fd#Pn{gfnN!$NY~eo&H3asF@SM8VleuKIysCrswKPO9Zt$BE<(@VAyBa^nQM>b zg`hYLhQnNM79T9&PNHSO_7_n*P^)Q&-L=nBO;i#Vi4+AQK~G)m5(LoIO8bC=tqcF9 zE5u6e{$uZ>>(xbeeF!m{k4&#Z5-{FtH;m9Qv(EygO^e!8#%UV?RTWaLri}{tkDW}o z`|w`@-l1PULYWq*?HqT!qN&AV>-D9Ybi}k(`$W7M@k=j%^;^uLXjrveRoj?6YS>~| z7d`M5hCT9mf!YeuAPY8uD@qbl*aw-5Rg`$|%CXWd>4R4y_6G$W`11B9QLa!X9)}wb z?2US)*HR5*#*a!-*;^<7F6m{v;zvuJmYPbt)}>_Y{UVzDk&2Cs%0|@M&8}0GwlaYM`~Mxgy_j7~Se@YRhNO%|BCsDejw}Sh zyvkj!VrOT5uB5mGNzM}B&Q7?%xMOt1k{T*dtqAwa%qf zwn~W3=&v5!_0hHX|LgqF(@CAnlk^U;r?nMNx@=>2O z1EE#hG5YL{jDd_=1*6T?QJL( zMBpagGsw+4gLVZM0x#i8?qJm0cI$XR{T~g#VDzcpwG5C}VG0Ax8&(cFm;*X{>L=VnH~+51@-PY6E94{}nv1mW#1e<-K|WF+|Z7B>&Ej@8N4lum!ezR=q_$S;kL ziNvcRyQ5!Y(qK53(|B)Y-$*v){aimSA2f{WgiNJOv%$jav&)mF*<&%~{OJT~r_GOe z$+5ALIHGMmtS8Iw-hBQnkyT219~GUVC^T2{J|jN;2gX3HWD0nHwY0?zwc*l1(65oE zh|+KO{V&NH%wMcxdqRoCJ0dP9t}D`uxyiPK@q1X`NTUdw7nCASeuho7+y<(c7ckjg zJ2|~anYJo_+KJEmQ`2{pE+m)8;jLwy+Sqc56y2>bJR$Jt^`O7}Ac}QSeX&5leoP)_ z8e-0q^H-&Z0iPJ`VCb+}=AGYDHJvQ;G>WPnC@udYNJYrx2bM;p67Um=}PgJ!Sp4; znN$YYg0qcu`{0d~ON96Lm9K9YPbl>=_0cR_ayKFbIzEVlKk>aIm*s9ZB^7?Q`xx3v zL18@oSB@C?tRw0`Ws@<#6j*8vRVhPb+0Ojx+z%OF68LH3aURkB!q5S5Tn|~&NUZbL-sG%4T<8}!I@qi{ZCH`46M>f_3VLF(h zQ=EfGGwld@-F#|_1>W`g+tS*7D#)q-BKZ8zeh3#oUZUINrRZ1Bt@7MN4VBgR&mc97 z!%BSS`mNF+68D@fA6YT?hWU^>a@$as7BEJ~pBK^44A;u3-WCB@hW_$j^Fs5_TREsu zv5(_{C1iaOdB|q~HebANX~AED@atO26~mjsZ(0nQyaoccOy8veHCJOe|8Y_%dk>G^ zeyFN>9IeFLP^B4d>;{Yu`Pl0%TIht~v^31y#Fc3eRV&R~R@~-g`;62)7w;VZ)0s_| zkR>FS3G~Mgk#;@f9DJvqTscUN&t+lr_IewHaUgG`Cx_Y#T)}Fr{2eZlDIkKp=ZA&| zCQJ$Z8te=|Ks$1!wwma)pd&r$6)tQ+Qno7c!KY*SE%;Bk)Z5al-RxPJ93-Pjg@(JVDXr_?#f!h)ydjkAykEfo+Ix|49toyCQ(m{ z!(nhzQMDM>qY2H%FA?z<<431B}`76ez*>APQ zmKtTXl=0$34PxNB4+ak9{pc5U16CxMC}X-0y_V8X88DBsDVS!prP7RkmN(%e_}OQKTt9CN`4llz_Q|W^>fxyU@75Uc4-hj>xzn)N)GL`*LYZ+V3=h-ZJQ)a{F#pI z#ZY9ljK|)XY9NYd&q{TDqLP^!$<(9Dlcm2q$V6=}=$Vy6V*Pn_Vx}gqvf_l}Gn)f= zy?v;lkTomSv`619kYhoh=4&)M80+qoICY#5^IS*xCfW<7OY^J$`PS;&uQ%>PMWEK% zwl(4q3bOW_y}eO!T>=3oDq|N1~aS947T<^~1xBC#*d%aiB0ZuTVC9$V>c zT6n2`edpRrG1Bj6NR*_Vc1c`J`_Wq0C=5+&(&bF=@k%r+Qy1m!)0;O`5lO*Mh?wB| z_eU9D;txLJ1{+@&@$(a5XL%hUtE2qX(b3s{$JBK8^wI_FnzyC@ZL+WxYj&Tm8l@x4 z<>}rc?st09@e&*NJ}-(>`Tun~c*nRYGAhm6cW***vc02mz6vt+p3`*oh1v|6*fk%& zF~{O~Z|EWa5|ERH#+fxQWRjFI&jG%Uc)9_Dk#8ejNS-f}ey8{K-cmg>*}U98`QRt7 zGS{5IS?#2*_vdoYtQ_8fAJ97Lb&X$bmyS@EiNS?Sm0*6{DqyCC`Lm!%ep3b9rWgAn z45vwNySSa9<599BveyRdbGWHclhkH?#8N9Jwyg@q=|Rj7nNYY&38MN1c*JDpD&;_< zSzP6~;$GeogIw)>rRqC=J>_q0!@tRk;1W%r zFak3y%jn7r!;U#O$CO*WgdTf$9?QICw$wc>)8$n>Bv{~hOZNCkRYhl*Waf4(v21fW zTsjFUw!?~Mp6bzoP-#DJ7NcIXAgn&A^QsUZYiEL7JOt)v*d)ynY=Q)1%r4Ag&Gb!2 z*l0gb)k6H5EIYNEYKG@GwqHqrjjKhkr!R4{yV?!)6U}vZVL`DfhhduWY;x=f@$tHH zG+=@yRT9E}0R}d-Jhs4jmO@9>ce4CQV24s7wBOW)j@}<^CVq6#>$>(96z|-`m$R_w zp9c#=TlMA0|5?O%o3@O<`5r3fkoS$6?bq!*IbLhqm3KU!Yx8S+?BT6${U$faByHen zP%QOTzitis)^P6*zoZ}{XnL~m>sdos2Nf|RG9>1#M8`c3}++yx7U;`u+vFP zr#Gr4V2ES{|3%PSIy_{b_obt<)HGV4Y7n<7wT8>Q{VOkvJ_g7`#G`lU;tw7a2f)B% zTJVbz=mY%99}OblmQen{tj`9LGl|6_j0a(Wl*rju)ko7+{jk;KX$s$A)68BWfNAsJ zou0jU^M*zYm4(6Y9ZL|g5h`aNDoZ~$P2W9g{E+Z1wlNyE>PVQ+G0`^`4tXsBz;_ma znF*!V@iY82pe90z;yU6`?{?`GzH4-SBdabzeg40NSsVdy~Co=oCsupOuXl;+Qe z8tq`K^4q~zz4@0MVUA|N0j`^-PD0u$`;*hVcEbc&|Gqh$)ZfUcOfuEWV(iG$ zp{+Sln&-x;Ks92XITV)y$KW33-o~Tzu`!Wy>hpxPR=Y3qk4j6`fy$~r~FuZD= zXFpyBQv(A4Uf?Lr8G5~yzVKCw{i@Xbs?ZsFj529%Ba=vpv)5=~nlfVuns|H-pI*z( zA=zodGRC^42N>h-UrXo={cEJl^EQ8EtRuOH)z%R9T^sErAS(1){u#20(DJa6(?kTW zstlq;$ejdUYYD(7n#o?ry!vuA7E0xgpdtWK67V#Kw$lE=LA9gqmIgia)r$b2D-eT4 zCh$aa8!$x508})&i^#QTUz1TTXl|*Pv91CcmY#%M<)}|uit;*@u7bAekBnt~^JWE& z_@Nq{+&@u25I+S!n`#pfooSy|;|;VP_6n83ZZa|=8}#z5LLtA+U4QVFn3?ttWygZ= zElp8JJIEKo)`?>~dz$^t!|n=q;X=_XEk65h%YTRo)%U^-hz2wY>PvcW#M?i=u}h0T zLsJ3&=>AL+A5D}{=i2x@x44&{4O@SA6hC~F@>NC^u69u!MDz|7qIKsC=d1gQ7Z|zO zam2Y7M*Wet&Xwk)TEZPM7qdE@5g-m!_^rQO3_M%Vm`%$is-9(g2W^rJ8dixWs* zXRgUVonVm|Jt9?+s=A-0GzZtHh9D-zvkQlA<3{Lq2n&7+RpiIYU>|Us;$*NYjc7^V z`XZ%;rdR%4W{|BYO3DU+Dl6~yX=%4!>evt+%68z3W^Jo7u<#8!3#%r25*$p(=m=g38(aUuF zaP7r6_>~-eczJx=rA12o($%4Mdtmajptc{DFxy)?Iz((SGgulb0Pn;FqxV{N4{0%B z1I#i8C8DigL23rZ`gt-67C!$M@YxtC%Gppx%OJ0np{#6|cU%T0CO2JA#91U>)FOHeJG(El?7TmLY3m;OS%!v-)~cMl*&52^XOd)v! ziVh9dOoy-6cW)iC&x)xnQT1z`Q~*PYz9^?>h$SepGL*-E#+1PWZa`3Vz~k411l^?? zFOoL{KU?Yk^LZStiyGVgrO= zxDwJZS&DHqP353bwzO5Du9AM36jv_^pbr2HC(|C?t8!e;8LJH4=C?MI!n6p4aO!8v z!Cn`+u!Dv5*IeLZuWRO?%pJsUZq9AX=D)!nVr{>CAk__AZl(>sXq)=HPL0UL+KB1} z2zX9(WQ)r7rEpkD$o{P7t-j!H5&zrtzurQOrQD^>g{2pR?4niQsfSnPISTyl3h-rL z`0ykjKPKdfDK>D~yk$sLN+rw}=V6nT$ZO;S3}eWoYW9)`lY@xOeGQ})NJ_BTO!C3B zqx+SUxE!QJu{+#CFxWL_RrR7!>iA*h`H*qj`OrgpH=K#-+7lETsB|+ctPWCQr*@t` zxBIEtMdzur*A@W94Bf;G$A~77^S9%+=HR;+%ROP}>@Ej^+EBmf-IkmDV;8)c)V)7E z7g+DxgD~r`%CRw7?ls1Tn&3Fe3VH=nQ4!@`Vp$$hVa34 zB|VvqLm=B*OP1wEPG0GHPnd}X8{-oZ1eG~34DGian%f`wGJ$WQpf5kgV87}xt8ixp z7rUnoO&YI4AsdPH{VpK6oQaiUqN)kJQVrO3yV{OE!gT^lP4zz#6e7-n1iALA*Ur$Y6T4CCfpsMPQ7;(EXny|ccL4YE@Swn^;)C)7R7MHSk7Yuy>j=tA^D@Mk)Z;g^5j#Isk+{= zWLVaZUvK;+dxlDaA!>)XqE>UNnsJDm?1vj>K9yQBKY8&TFAY{P3H20!=!w+KRTbsyx=8vh`nN zU85|!*q@8ihn-+x-kr1OOl`SMKtGYM0Ob2Hd+p&@-`_15Shj0|>=zl!JK~p7yE{1> zJ2b{NczP|})|8UtAXDRE>k%wAJc@(eK9djF`sNG$A9tMX%tpp%2ZV)1$!5oaiS-c- z>B1KP%uW^8bp2+q@8{uP#C$7kVyCwRN{j61N0g&tIsPH(cWEC)k~T6CR>hF1!$s5a zvj`l2-~^97jdlGHxORTIp}z53?Zhl9bbUF5*4_SUI=r#~XHdC0bU^AdTBTPKY2|at zf+1ve60qfB6D6`l$bd;>DA$F*A-wH#DCA}p8@-ieJpd$D3?$oS}LbqcF_yBf$Zg;x4{f z>cPYHt- z;d~>8s^f#jA4ht8wj#zD(=!6j@nQ1e%Pp;)qN6bXR90GK%j!upLy$$_t}TA7y{rp; z!pR;nV5c!}?#Xwc-M7y<^uYBHnMvpXHIbxDx^e_;UMJZ*5f0C&8ytLvW8z_Z)dXcN z(F%}%Uk1n1%3Y4xN`~eDmFM`Gsv9%{koKP9!2+y6NW158Ngkl!AmXy;+8NO}FV@Qn z9qXrco4{e6_tS9QCEPbW8qkx|UkZMB;0$6VE0yiu`2I>opbzV>flL zr%!XRax$J;n+~BUUQ(L|Tm=3_@!?13&fznuE3>fzfs?o!85+kKPhwdJDE(vz2#CZ} zS)7nHyfCI148MEcTCGIionLXnYQCgT+9Xn3c?kcy>=>9xww-JhA6<$)o1F9mGwIP) zAIqp4mZXKick(2qqbn0nsUA|YWUHyTP_*o7wi|+}>t>K^kb%YrwhMrcG=`l$U5saC zG#{$bI!NrNho7HdG|6-ZH(yhPWS05Eml(^i?r?r)2J@~A-kqztdug6|M!%pB zt@Isb^H@%ispbhclu+dQlif>apCO(dAhS}Z>Bjg7#{ z(6|x^Su6=4x6$ql91KD*k9jo?7(Y-@F&ToXrmw0W%)mQ#9GD)s&oV_&NPAO6%j?xW z?Jtbgj_h^_8LIY(acGh=*0r7F;u3?-{xkgAR~d>sJ64J2-d6lkzWH4$*VC z$^Az#YN<+^)I}>FkzmZbL1=+99BQM;Ict%j;$KX)f&0Z)u4BzU*yog-=-VfP-Cu`{(+kA92KH`z&m%A4m5Pw(4@OQ(XAIN_rA2uRz#*;C&(wjsBL z2nc!>Nw8jL|bw((RYOl^Q)Gd^6WS z2XZ}Fc9y=a4bmqArY-Nhvs)BYgWN6@*g9+JZbkWA{jm4{5dZq7(l6_k_Pqh;Tdx^G z0o6cdYHhca)1+mBy-F%X1NpZ;!I{;ZbpC(YP7Ma;ImuRpGgBI*LdXW*-(`d|UxynL z==OZ%A5<9O#8e=she~r)n1J%(e^sERbF}w$bbNb$`@OW`8D=0YXr-$bRIJ~~97q{Z z4S>?30fHeI{Qb?LTw6hgaOqA@xb*L3$Oou_bb?HMxOF0n$vDU`YosV0+P2%|*_P)8 zmf|vkV+_&iNL6XNdsUM{>O(aM!noYs9={kksEO#vr0#2*DOM>9FUA(RK6LhbptL;9 zEAQr$M|KJmQG@cn233Yj?%G%Xc!%zbLs(nMW7(nk9VVoPR_26URXs*Ri}@AgPX@ul z#s*MPhw4wj3`h~HkEQT|zu5F82Jj`XC)arinGjB4; zo=KD$82y^*PG<-zSbpBw~g69f@hjemR# zWTisTRuX{2c}=Q3WgahiO%M@Oz-sjkRx9m+LJ$mP$WfwgHiRwU2jV4L15s=rI$#MO z`jLxlClP`S|He4mm#9-t8xd*YPm-ELJ~-|6$>VVWIezY!S-$bxx$cM=tJFEbC&qfZ z!h)XOTIaHQ+EV5b%bt(21lUk&4QejWCmx`^!^`1~XrP8Of{u~1c>QuXco3$k1C%Bs zu!A+v%PZuD{z?-sf~^A{%l{myYI8!(-9YdWYOk zLSWB(e8d4go&*8tQ+&6O!{eUG0Faagh7mylZ`lq|3eGqlod}`+E;ldL{l{bv84KIj zKh|k~HWOu2nEDYx?3M#;;If(`aU`sC;PovFp)0`rKl6uTAog)#$C9{K&CF7aLbz|) z*bZnverOpCq8S~MEBS6XFJ;PWWK!D$i*(JGiLh`HW@lC$Mw(Ft&^kP~n-SRb!g8Y| z`~0>1=P=U;Jdwr#V^zp7*cJak^=wO`4sMFp5vZwH*a)5rJ};|{lf((Kurc>2d1C=2 zT0ncPt;7;mRZDB<(qATJ_-c~B`P6i>dK6o02%Nb014$i709YWZIyl29<)YZrZ(+qq z<^|J8+9hPB_fyGTDk64n&*U|#YsC;x^PJQC%hBg!GdC2aAK=6#vbu&j0#O3sn-T*e zkIhR#*U%#{qRTphyE^u$Vk%3YAK*;`>Kyh)+?8!xzS@YDE< zX$AZ{b_zmDg~YMpK*2y#&wog7R6;Lr$Aa>LA7W}(gEq91lpKxTwYPHm*OJ-gutiC|{8 z#%XflHz3H6@$e>M!BXPiOa7(T#k+XCuMAFpR3yqC*v(jdK^_i3^a`BmOLS5$Fqv4> z>Z%_}$G~9$5qT7bApmxTCo9--Ve5@Q9Das=a!*1TP#wt7^J1j)L-zv)M;S;lAqEg+ zd?;h53vvJVcd34dV6N2&Do~`V_K}5+%BsAn4b@R}eD>_klsrNKXM6}{uvup|YXTSr zzBwYDB0^4%j+=`t(^+i z5$gq;tMLc`1`2ARydZTe)s8I?%q<#{H{I3!DnFQ0E>e=SSqkN!3&`7<>J)BaYT~8_Cz(<5Cy!z@si+t^84| zyor&J(V$LEXjy)Os^wFi+kIg|FZ5UI%0ADKOzu-n|253Ge!%t%=k8KhmD8JHG9h-3 zif>+;O*8Xe1-?z^Z!jm3Zu}&+8y6j-mcuEgmL{jrge^ZfR~5`|iKsoXNJxu$b47(% zIAo)D4KYV8PelwW6)C{L_2roZ!O9}TV2b`Gt~g$5 zB{#E$y%hC4+%=}A+7=wS_HIY?fUyQvdyTl4nqrf61!#iM9Xp1tQd88kakx^9BK z`c-uz$>Vvg(NChNu?Oaw`^;Qrm`;9j7y*#i6>7C^H~0tpJ)&oBW|dgej6#bgUYjM=J?^%X`U$XnL0i0`2_e|y^bX%w zI?rr-c6_j^V;>jVh1P*J@Xt5s{9*w%gsY`pY|hKbM;lzs*)+V(cbUAd_YQ?V<)2@% zhCFR9kmb+QEFI9cz-N|aPPL48W>LTQfT=u2cQbLdsZ5yB!BVW-lMwfw*fNdc_nY_! zd(!pO`7r15&OL~OAlp()br~2x%Aexxfp-7*H$wU%#MUVy#QKBIogd|?ekWx|kSmB9 zAW=SewhYEI|1j$NvV;#Huy>Xt(y`$Pu+po6xi2I14Iya1F~M2ll5GhQx^Xk*8~xy9 z>^iq#bH~~9K&(oh!qCMfxr;N5O;Pt1Rc1TTjmwg!Z;ZBjAkq|x%hC&MJR@bxjRDfy z`~A9&S`)B@UWEpBBdWCeC%gzRhnjeca*KOnp4>g-zh%?MW~`7ei<+U3AmXHV7ro7L zi495je+NS*kg1#xL-f#NBffh{vP&O%yrUkKR}h}-rmfH1p+fjk%xRkVbR2n_l@ikUaTKb|iJ<;}H=^RsTuk)$5#0tD{-Bn!QM;DxV(uLhbLUfC6$nG`G+}qwlg`1s~Sohx)`yCo7l!f=E z%tb@6%oDb0>leK$;z5EHP4>Fn7n8(ZQx#I%-42~YeZ<%RB0BglOt-332jQ;Xv-6)F z-qHu20gxE;;|HmqOqYHd2{<+3T;TH&*z57=F0d!{^HtZ97Ki|mz*%ExC@LKh$_8pf zbulMl_@JQ6(g-#-wRhk2=GxL;xUDv_>`&fn*SYGer-g)JgtpNbn-V3n_3nF)(hz)} zU+=tgY)54!=XYQpe>h4{LzER+@AOAR!TdnMKQO-vu_>_K2x!(jP;XIX1BVDrV=&%) zDYYdgC)-aKD#>pUI&q{ z@NeS%G$+qeQ2w$ML493bK#JqiEobt1!;ISM~sA%6`Rj&r6cfvv;;fplDkCp@R z(jP!16v9R(&Zp@OQ!rNd<%xNzvH#@7Xq+lX!G6jy2+^A3F;%EdW!iFRhTTiuj<$Y7 zP=ysKea6-Vy$ zZBcw}J#L>IUJoO!*u)nB?OZbzEW z#~BuTOdP#KO+=!=WO45Az^-vDsHC2_hy19 zL?;S`H$Prb2DL92VE-4065bzXw8W)-!Vk8Bf7i9<*L_DHzZr?)#KXPvdLN&#v#|zK4<$kQPL$7Jti| zEuP*C3ovS7+VX0C2# z!{pg?;L%j691SilS`ag~L;r{EJZ`Ut`$S$na&D6_mK}oO18YvWpf!jd8T};~GxRni zaK7xET8Ojt)5Si+Aa(sxgWH$V(Au2VgSh~?))-@f7hVdO?o_{&QH$divW-zgyJhpv z&?vGT9va(kp^^7cFrhGz4Fr@nr_>z*VMGbx*NB@EnDW*V)M^6xWW&EPm%l!IH$J)_ zLT+l)HexV2sgA2=qG$MuI74K{>*{C^TA+Tg4AbB!7%okU`&}!RlQNSqGe5%)t%Cl^ zmgeftmOjQX>+N@N_LLx_S07?7`(|aF@Kw?#k`}+qy}Sw!T!^Ud1J)Z8Mv;rfX=UXz zo8i&et+Y%nrIXTTwq4nQAqzCWnqx&%uO+P0gcm3lobXAQG_AUnY;PoVOq|k?wbF8> ztkmf@4@=I1rW!B`mL7W`F!;3o7Fy+G8na8yV*G(6CvpP_2ui*a+i8X&ppd$IktsY- zJk#luo#jCDM&;t|no3N<2{sitsqbzk24%ZNB%=0Gua9C8}5>oJT z-8Qd`jH~x~aMo}JqZ+R9{>C2kr)%Lw6vac7ZA2dEasbO+G|0xRAKmZl$}$$q4n(&5 zw%$LA)hO^ymg?QFU9lW4s6{%W)gkf#h;MKh*Uzwb&U{36a(>}*2QJD-{1AARnm1KaSRt_+u1mf6qtCcV^> zI1UYUxkLXj6whKktjlj|+A_v6j4ULLSAd*==&xDD>bg0w9mUaA%^D^O0qymwF4s)K z!T5(KfpJpdOx|y^0zu$#!Y3<#%))ZHPhqXU&^$9dawR{@8$X%e^S1k9NUayfhsyR{F_76sZjOSX zl?bb4<1Gcnv%i#SKT=c(oBWLa-n-=-28hr!ceE&?bJ@_+_Vd`~K4G@(lwiVvuun_> z(KHZg3Sk{-7~d}}E%q+Sl{dmV{D%Hx09;=8-|jo^~>Owu1XY;jV@$R7w)^+EWN9W-of>Q%Zb}|Yz)VJWvx-a7(G7H+?KscD8LPeTV zEBS!;ZeOL z*Tad$ zTTfyHscOO1EU6U4YfzV?BK1em`J8(OLP=gQ*JTf!u{Ikmsbzm<2i(CKuBE&+~jOtwR5>#rHP z;ns{DA>bbG*F|~XvH!3};;z5>T<69ohZ!J-uRpF)jcbgnt6H=Ew(UP9s+2P&9|;mH z-L)5_!Hx62{bqs|9EL#+B??8iU>D23tt@RSb-a3{Oa;K2`ctilEPcq=&Wq&+cEzVB z+(*WA$jHvbc~+=aE~jT$Q0heD1tsK(vha7Gt+K7S7%xWi3Kl9XXTv2Hg3rZ4ue7{E-Z>sUuygszW| zQB^$M?r8Hxj#jaIPYkBlB=?lhnAv5{s%m4buaWeIGnfd{z^*aZhKW{mt~Snw1uJlX z_uD?6sD1$gIhM>?mN}NM_kBDTrM}9)QHRO$lbU114q@$$2b(c(FTwjTp=oe-%c@jI z)m@0@TtB#HXhzj{O-G~WU`;nnV1JRn=NB&CZ2LxbIx>jdB%_Zn z9svzQh@NU^3i;>9?^x0$bS<>)-15(lj*UDS&F#VuwvfY8q(^DX?E9r2Y)O4N(AU3b zCoDaY^Ip`e*GF<*__)_7oti3~7p870)eVzPN0Q)YE95TZIL8JDSkhN+W1tpXscF-> zcdLF-0LLIyT@ve-D}al3ypp!JgI2E`ca2&v7md>$1m}H8d_JFiNOSW3a0>Yof%qk; zK!x$zEbM4LiouIt0LX9QhsXTS%~ILh#It7FuDZ6}m3Bn+et9|%)D2j@lFjg!%8X)v zxQEzpr0pO~a2XwMaa=N{8Fro- z+T(+ND7>eQSAd4EAxJoNh9YJnwn<}NF7jdzKlUDQ=v^i1=Y5{i(V!&dVt6-~=z8deEAh~6DgWJys`lIo>7 z6&oGY918Dd#=J9PX30!5`#$}C|9!9P+vTe3>blCjUa#ll{~mkXWX5Hp;Y53g@B!|cuHE5&9&d(yPcS)&?k}m;p#@gf>Lei!8yp3H&-+ouVYYrx z@eR4)cmC^K;;EN=dh`(D>Got!#EAL7;lcD<>y4blUK}%7SjgQVv)C~IORU^+`d_Qq zo@0j1E7d|y)+<@p5C5He>Ty|lqbleHS&^{qs*C8y597bf=bNJ?z1M%0IcJfnb6(mWCXtR)Hwvuxg@gwz2THgO7MkL25paDf-Bvb&g2`;p{L1>BZ+1 zBc99eyLwqTKHfju2TGbVYc1s>GQJw_ez%-JlD@!VnW=%xk*j$NmUj!j4d+_qAj%$a zO;73uc8hA%=%>e|rXH-d1bN*9tS2)SZ-s>MiW@5n?%=Az{o40I{HG=-A8{>ubZR<9 zhWhl~DE?IF7gi*JlpBJ@F%@C@NOc~%sjwg)C*B&??gd$FG+EgxN-U?%jY*AT=qpC0 zC_4KJ__5bv12_@Rij&m^Z5xd1gKhmF%n$p#DQ8xVSlj2Kk-$yPK0=U}gTXB})8FogI@NcN3^EQTt3Jn!MtQ7qSAJaYRtc zJ}BPK+~pBhEl-QlDUy7mbH=1lFbhr~F_{7%CMwuMa=UGVEBBgkLc(&FYVe89l5&-~ zG^S$eoWziaKK@EpjPrLu&kv5XLjAN{yF^e7jWH^Xnr3kbq_;W(20?yxod83S|6_=K zQ83X$Wb?4+m<^nG>?&sJ_v5X(I8Iv7S`SF|o~C=Z$caP8bc$4~sK-a8Ka*H52&6n! ztRqwL7lw4^H7azb@8S5$fb0WM?^?}WCoW;SMeACq|BOoS1+YIfB2v7s>t*0_Ll4Yn z;;KTtv%R4hD}?8)Lf~LcM~+H!*Rcx-@_$FLam)r>EYhCW_}s472XY+aXUVx8B;koSVec3G5L4L{qHB zq%+3sAOh*(ht8T61$S|XjXsCmA&iY(*&}XY9&CGu$TbhXy;GF>loDQ_Ci$w$t|Z7G z3}bgO6x&{}s!AqXs^vz>9o zkN;xvOob?=dO5gTnoV91BhosrjN5Go|J)(638nHc@`jH0&`yj=4A{qHnNQ}Uoh-GP8$Td8%uIFLf!E`?SiW)WpEfYQeZasC$~|b3 z^C#sUCg#7{`3B;gt>6?B6`-Y1OsNy(MXv;UOhsAFEA>a8ttk47?z|$TFv# zEB>QnL3A1Zm_>t)Nq5262MMHl2J0h%^g&0U&E%0WE?JXJ4$5Co9P@Ei%E(ATs=o2u zzQ|t*FtIFTo;`etL$t6j=h~2O$L_yWF(##HvJwe$K^6NRL4J!UNNZvYHnHLA9U_Sq zRTYIJ3e5y2Ow{Q;l2tskOM|RVsE^*9D^ukazq(KS#zwBGqUJ_LN?pVl)(fmtH<5br!L<6G(UL1>s4bg14Ps&X|M0H-igSKKLzuA+b`xC>&b4meB^XW3cBN}@1JV9Fuae@&1e zdBsAzVPWzW2a*MthiwgTFb{XS&m)oVgvxs12h}WXjLd*#v0oQ+(S_Tj3C0xnE~RQr zdioj5fk0|n540~8hX6CN_%sE~L}Nv|>@rTg!)<l5flI;P{fD`I|WM zT*&RnZaCJA+0f^V+{Qy+ zcmk#(r`8)xMZ}wc>|dPsF>$&;N>$%bfRwHlvh^DMoLRm|pz=S9_2{^&duIj^OIyph zPBHLdFujBt!(AesX_{nIYKLIS2&8Yl*nLcJPm!8D7ty_f0-U;T+!>xJ+{;> z<*mwFP+NL_FHqro z?0QPL++xk~i~UniV%Qf6^5O)n3gfe^W*WQr$bQ>Wy_)ADU&8^G!zT2zz4G32lZDaOb##O-rCq`$D z=Pz!SI5AP4XR2?`B2&&;RoApFPC5wg{ zm6D>_U^F(a2BWc|`zRQVF0uE(Xym*No)}NO%%2OM@Pe{<2x&Q2-$+}vLj)h7@ViAb zuIKm!`P0MV$rX+trq`+Kdm2oL+3_CQSskK)%~W@=dLmgv1o@bwU=_2$<2?T1Y)s0& zHF#~_z~(H+A7rE4!qEB4I0pKNfjdNw2r=-1#LsLNMiX?E{Xt+Js>fO7%+@zIo%x*F zA^K9cI1lf*Qf`h)&sMNM667^s1lrYynXLX#KF#c_G>zVrYyGu^D`#s9yF>^WZ8{}q z7N9&`!zH9>>MpXu~@p(G)g%_hj7g$T5OtOj`Dt#~dK6*BrPEQPIs z-@rqYZqydxjQ8B(Pw&pjXV^S{$n~kz3GEa$T%dwcoEF2{X#^(J7W^}g+PpvoZ^#=M zq}CSW$bbC$2PN;L$r=t_KwV!y9lctNkM0F3al^ zTwpfr^~9_42AEf5HoS#|b$AP&QGG=%9f#-$?H6&{B4%OaY2mPcBYd_)RN6)*jY+A2 z>{J3taz_9*zPfP1(p>$8i#fy&v7JofAus{Dz<|@U8u)^eP^#CcR5O;%ACp!&(8%E8 z@1lE+Nv}AuqQN9cjYRO$hHep?g@|JLGvW=L@hsZeV{wZ!J#q<1!HeUp?aHFhzsovV zHED=uuDLF!Q)CfH2Pm>N9H5B2%s}A5Y&a4i(c(3h-6euC8WLJBm{vpZ_tpPWs46S2 z0Zm{aj{;4|H?|m<5)QsMr#qxJ4rajV(5N)Hg3Sg0Z`*{nxI>h$PDvV*{>^8T2&5h_ zK`Jo&MM=o~xF`uT{sJfP4QwS1Z5wi51T3nPbxi2otlws_PzKLkrmz!WwM)5q5U~zkV=J3-c z?p>ls8)+rr=-IJ42qf+&-EqgCyT(HR# ztYW^frZR3UN6WNR6sk-82i%nkwk3g7xKjW==1W6t5Hl5C zXNoJYn7PHZoU01y0!cW|r2|Qbxd$Xc{!cOXJQMZRy7e{=^TOlML5SQGc*q0F8fXE> z_1<&J=~3y;B=#$?U0PuGFdK%{kSN~3NJGK}!KW z@X*g*gO%FYA!zOI)S;Uei_BIV$X*5%bPTN3ih)Ql0=A<6cNy%}j7AFB_ei}ok1^?q zCK?+&Q2bT|>=>SjN(%;yb=t$Su&(h{7!1btr@}zlA8G3g^^RNSm2mT=I!%wby=!$p zcZsgWP=1a|W7OEc2UMq+g{(CMc}^i4 zERq~Yvd<7mH)YsFW~%r3>Xp2OUFl)n^6f5WP#7JD3G%M=g~SuqKrNzYHO9L{xtWxi zQE6Wb>+ixeQGmp%*&xf%#%+&Yc)vr15dD-sjbv(4_MqiIyD{bLvUkr99?%03QfsiQW%57&1k71Y%gs^I39(wo7J_b$<_b}KC z#bleL`9r;5ZP!2KBL34n6ky|AclZcybeHv7=IT5B6&UNntJ6!&G=C z;v3qG$%p;>_d-y-8SWvM^h(>IQ^fc|DIAmDRASlPEk^mj=DVZomw6|w?o9Sp*MlYVw^&j^|$5s zz;nAq8}#XIDaC?o5%jhO`>9Yq)nZgS$79_jkiHKHteC0qG$cECMi;LVTc?*ZPyx>4 z)9ck8e4N2(inJ;&=eovgAUZ^|e$&ol($N^U4;Y)%f)mV!Z^4o^ya7wMsjA<3SJ*4^*!n;E9FdTmgyM@WD)?qA*q- zoFAQ5t3vtD?}yayA=1jYVq0zHLvE?PuJQqw`&=iTecW&;>J8;Rpg1SC27y%Ijy=Fk z9X4J>iQjmD5{=`p07{@Io%!w%qUrto11^WJ-O?%QHKq)WN%8*dRRp<0h}{oPW2_{W zshD#@YVllb{>S^SeCy9t9JWNR;Vsx+k*(tmY^8>T}Aw;`{`F*ak9bx{E+YQ%-b&Ar}=+rT( zj|ba~Kw_WA9^oXU_$DGZ@QfPs3W{**|L#5H28mB^DSOD}5Oe}NMKgc5KL>Ne1TX|@ z!T>Nt&kAAJpSC06w@* zYWjdjNS(3;NMaNHQN6L?8He_0OsbX7-bNt(TPJ`B@+$@;tf<%!hn4d8dNDrSc8c4uKbm{BtGU^Yq^-b`HJF~yi6P1pihA*|%n&p&n$xQ2o zy%6WEX^v!1Gsi`%?V{n37n%YVxZVN!aJ~BATUYAvnDo7gPo&zpD`I$)_R*pPz)?r| z{dOF_Ib85dx;u((1-wajR~O)Q?v9+N;uy|;MI)Rxp`T89x2P0O-99RfAhE*!KYs3z+QGAiF^L{7`VQW+d z+-Y5)4tF-vfI9GZz?D-eEf%k}fI8d>p|#ELCSIKp;F+mE&PuNGZKn2c*gzemvjVFs zwfd8kg5&$h|FV{t@LW#a1OgHNB3*b8vN(9%=FBDk zl*U0dExP}h6!{ha2I=x9$!^}jlLA?@d^`HZ6*-c`@RLr&jv`{EC0a^@*3l`BL z!4IH?L!+(deXckXX4)l+fEoJ11_AO&Q-I%;iqQmqQ`)OndKtK?BYwbd@{e%O_J(@n zH+hzFjUnyQZqarw6aZ zFEla=*yC_l*D#ix<$UlRF6?< z1A_GcxJeKeDdrL%aj>hIir~@N3>DjhP`;W?-UIGuGi@KR7u}prU^ZO1D2ZSyyv?gC zh1=06inyE`s;vP8BAiYE0^#zVMI(@=Jh9**j`RUMG4S|x7~C*;3?zmjtEvrojkcM6#W;?E*Th~ZJ{xi2 zA=lJGt7$roF7&az@PKP!rsLEN_6lm$sMP&58wfw+5Z0l}%Q*#I?F7!s9nDglkZZewTnwbRTCA5lHD_KqpWawgH_ec6{bY z!-@Yv8~8%TpDlq*2;1#)O1Nrp?Y1rviB56-o5HfW4Qrb-V7UkBwvm+YjlXGFfI+5D zfi+A{RTWZij!IA3EgJQ8QUDQ(b6x-$nO*jF84mr}Aq;LL3JeFX2qOFQ|t_7YCL`bgxcFoO|h1K6Q0 zUJKqW@_^_*qtf0w)+2N;@m5t5v}y_t@>d7JHm2gtHvDd0;XS;U&9=di-#XM z(E&hsH?Zsg0K8cW2tu)N3=l*b`D!iC5V7a3WICbEf5B|wmqM3Qe@i$qlhdp4==Xp^ zBkd;3xbjL(^7L*B<7gY$;dZ&dQiVml1oC9wB3J({gxX2*?py# zYG(A?<&p!{1jq!`ss&^sP>TX&V$01%6H2F>;sx@G7ufE2! z@^E&8^Zw(KTfq4_g#;AgQd>{V#2H8V0)4<#SQJkuu-Y7SIV4PT^e#A}LnMlz!N;WV zvn)UnbCZG-Rk8D&0fA0ctg|Fn5X|T$etFUl+_k9~$sS%>Vl}a4dL6$pp#C6a+;2MY zfGaoDj_(##NIBh$%J32O91Q5`qifwESNcE1pi9Tcq(`dR=YaO@ zUrg}*ycZMvo!1r<{M%@d;5X%mg9JbKnVBbeR(o|yJ49;w6!0mhohSjoOyRLkfuJcF zO8}m|IdUb>=yZkN=jo;KjZxXYQ2ZX7R^VUv>w^4&VN3z}!%u&gf6UG6#9vF5NwKe?1vy-_=?@eL}AsDyLm21bMbqCuIKy5#)PX3-1kofdF zW74}OPrt0*`q`Bekhc3(!EKy6+_T*iiodoA5Co%5H-T6riCzt22+}Vzzz+(~HH-YP z`B?$5Gpy4KK=@YV-VPGjk9MFYz@sf%PPZPhfaTN`4{8FWe|}(1n2PA_AbuEV(gC

v8d19{9G@Wb|5Eqb>oz>@~z2d`)lKgcg11o1=a0-^fVth@7Pt?ISBh3Snz98gNP zY9QHvuZQK}5cos3{*Z-}UC4dz_9!izE>Ub1?e?fN(%}xm1v6Eb#I{B!5!oFZ@>IjV z3pX0s)`wlgRT%~Ve!x&|d#2JTe5zH?pE@CW)}FcnuL zb7rP~i>qG2OG|j(N5i?k4U6-Il&|;4xkJ6-=5J=?YmG>yz-7Wf7~nF^I+F=DtEtCZoo)`y~AAxiw0vp3re1uigIaVniuSN@1Ccs>2(uEj6?hwarGoXiYvW;|lF30t;Imjx=aW?DkbD`I;rgClq z0jAm|+9Ia`*qe2wq>o7(hgl{B(le{Y{1zCyTbYVxtyU_I!#rdFJf42joIeYhS*|v~ zUrNW(OHt|)TrSHPxEJGSYP7RbLWyTrgYem3;K*zsca@M-48)3c@ev1#_Ci(GR(=n; zj1UcrWPVHjPA!XW(fCnXV2Xk`_EiwXROGt|{9URf6nFONMofzQ0Dcb-9eu+w9cPRR zEIJ55MHXKlaFgsb$YorH9V~{#f9F%BEtEA*?Js1%1qM?T&<7MSm|S=Zy;k_kv%&K{ ziFgxU<48^|7&Frc06y@OH`_~+Y5vJp_A|RF@d{Q#Dc39W^m&2n^c9>~BQVkra`3Zy z4pzU1nxTI$PVw%rf`PR^AvnWKjYQ+u@zPum0ETF@-1HoP8(g=yTNL9%)gP6v2xea+ zkp9aPm@`xV9LCoOEz)9ZYuOw)=Dd5qHzX!+!Ig4_sT$CtTe79rU~qj^QvGp9Nhq@+ zQ5z5EH99@5ZR4y!Xa@F!*pp^AS<2^y;Q zjTq6*B^tjYX={^#2c)B^yw6p=r2&9J=uK?{Uf)mlFgR|Rg4N8_lNOSPKy7t-tU!?T zrwL92v1Ol=tmUDTUlQ|hm_4`UAYr!E5Hp4<;22de3^Fy-hyYa2>0N(JAf_f5N=+5qEk@QOuqJg}yMa+g1$uzuP|0(RQ!~A+W0LPm$+1xM z3cVs6M#pq%NeDj7b%AcGal#R0Vr0f%JYSF#pjPYXB@TYJveQ zDC>Pd4pH7`19FIpiNzd3E{1~~0_nGyLtKbqfgIur9pn)5-*S*cH2ez#=Ls?Io$UwZ z%2UTXkq zNafO=j7sBBED%k1US717I?n-HX(07sK@QG;XUHLch*4w(qUjiW*miK0dbj!! z5-w()4xoiLPdcE5%fqZqpb7p>;I$1y%vI>AOD8{c~PlI4bvpTy<)ipE_1k+7D zdu@H7tWyZ!Dpe&dx=LS|be($F;Vd~_2jwbRlfGChyC1NaOmlOwI>2YB{D=)^D!vn{ z?+Qy&FNXmSxM3gY_;?Lif{qV)U~|v!WLgu?yr+z-Du$1P<-?DXI3`8^e}VObLToC> z-T!)Flf4l=9~IeFm6NnvKY0oK@d^y ziB%K3%07-@eZ(!qZhNc^qgU22iH&8Z`eQ&8k#@GWuV>nrgYu38uH0SQHz4#S9fCQ5 ziMKR?l|zufO%@~>`8!YVtA_E?PT%YMS8im4E$m~DYgaY61Aj^Q!xZ>S7(gvCe8^&n=8(E&+>N&+=KHEjs<69y7(hc%0dw#pXPagd933C@9^Vepe+mDob| z=7-<_vw`RYg6)MU5`G6St+v#$RlW;dLU-(u?`o_ql$r9-HIO5a!0u23;FlVh`AyX~ zcdvt4TqL)I=4xAjLu5g}IVN>Y(<}$wf>3%O*li10pnUg*%Q{Vvx1{lZu0994fw%(4 zVa9mAgMrQ0KPj#8*WfXoq6iEP9QTX?8pt3}#;n}~vr{q#Pzrg^0U!EuM-C&G{|WJR zc4qLW4^0gujwe_d16h(MIxG)xX9G}T3udsYRJO3UL ztm8mru!KF@AkPM6%feBqitk%j@DN<-5$G~g*-48mF}((`gi%U>- z2F^Snwoe;~cSSk-L5yOHVz8UoYWss|GM`>PDy^)Wb7m?=T!4lY4>Ey3vrIXKQ# zkh6IFeU>vEF<-;@P{-QX>GF_u?9-g z`yE!*E$WP*IDk6oSyrJ7d@ANFJEGW?C>s!*Wj5e3lJz|FvFpTqoIes>4@TiVYak=~ zzAzvoppB^xm}MHMi2U+1L9{rYC{hu?Kq4I=c_LP8(*7G>rOduL8<7%D1f~<(;iU|S zy~?Y4(gwZHaLhTP(v|D zZJv?+q4DEO(>x$D>>ywZc(z66BQC>S3*-{)K>E;_w08{T5~R69i@8Lc5lAG8%kJxe z8UPxdy$^yaniZWQRL&w&WLf~EAXQ8Qq)_A@KHjRTmyRQE2+lqTh0+nN54bgEI?>%C z(;IaBs8rvIT}vPh?7)Ki<+c(DmLu!Bf=jsgu)T*y@$a_Z(gd+Y#b+vrB|L9~SVBJX z6ng}ufKecoK$qPCsYF$O2uLLmXAvNk*nU9^q!OV_>Jwn7DA+pyG8r#|g?cobKpBnB)}8ty3Xh`~l`iFE-@;5En(!`Y#679&6i4{3v5 zQ)oG51rX@%?0p0|l)!2xkS2d(HNnl_Q*Fsx5L^n&#wmx@t7~VCoqysNb8N|-up7Ae zDQ{a3P@B48zy1Zt){xi2T7{ZM+iPS%3lX_tJ;IH}HML=&t#!m5umwEcGGuD*O4b?E zkcZsR8EqR7ypPask4f!(Sk?s6qeLJeDA~!yhSSBmFrM6`F*!G^MQ-FWnjLl%mz&~N z=>fI9wds*IUS*`2fEc@~5k7u#KtL8h9$!{D%bC87S!34&yh;yk;6K?3=$eN)Ee{s*cGW1Fh&cniMazAqX*Pa-oEqKBHTAabd8;p318?l29Djlb}zSSYmCg_Cl4nu!+{!Ew&@JKgJ?J zMEv3^HXir`q1&p}!FqhV4;01kMoj?4;%ou1F|MtxDyx4NAF2Up!+3|b>fgqx(_dI_ z1X4|g;21Nte={BwEKfdj{3iC1y=%#a@zBazNAT2L^Z>l!T5FN>XJ*)9I;&QD(2|)N z9Dsi{dpyN!gG7~?dej4--L_O&TISdzG#scZmX&hgDWX+)1D;Xp6~|&6V^4U|;LMdI zh7G`uX`%jki+2YJ^5H9j|1nd!5b&ZJ4zB=S)Iuc<1lj`?8IE7TM`+Fy5SXa55%^4g z=~N$zb!UJfiE(FzFQDgh0|Xo+@_f4^eiyItc(r2|j{L&AXnZbRc_FaK3*t0^DKv^P zBxw#VbM=}(wGFVxg3m^fUMz(A;E_CZ>9yLcIQ8y;NN>n%^6RXuz7_TdmftCoXHq|d zKA|fGD9MY!#b8^m0MUiwzYX};y{0a*yCA!0BDhD6Cm5@fd?JOQOU8%ylE^JiC|!_V zq++Q0W70t}%Qr2zxm(qEe)l_&UkqLs=we(4MEal~)wttnEjVq)KK+yU7S|Zvqz+NM z54Gv9VyQUNy&pLA6*eH>_IFs!x7`b=O`tZ(qyq*nPGvQK?QcQgf^lt3c?rVpstwrG zzQdO}kGKg}b<}|;cNLUm6$$o|7fo>gn_kZQpqNUl1wSY>hiQA}61hVZgQa*ia0W5o zwQ@T}QL2=-)STfMFNy`Aa~SJpn*Y1ckDT*2y$d9eUau8^RqooJ)E|O+BW%&YH=2Z3b%CNfz1P`=_(j9( zh$>ORNZ!KV;8yPj8Plu>9&**pbmBWj0pm1iOv?X3PXeZ`;vzLkQG6+Bqy4~iT84VK zEl?Cw1k9MoXsa^NFzym1<c=odF7S3ODC{nrvu*qjrDyLyU zqS?uY6G*c=1U|s&`uvgtP~jGe-3#6?;0S;R+(xKgvr?9aQ-5mJ2yin=E8XIM_sU$4 znh(NM7mKi1Qx+(~7XJ&<>xR$EKzhA^uL0?G@o7je6UX@N2a2%9>q*rYFAV|j8?z`B z;@?!wQK`cU7U;r0!7g@TWwga?;+#9kCKg%_gY+8x_k~{8|Nj;Y<~My2Azp+7gqZ!Y z1t5f~(IP^8=PV+`x?c8vf}95j?*TNV&~_OvJG2iX-wrwC4YhsRVqVH^vag9{4Y3d!yT>|pD7Mo3_&Ap zC+5;Ikv6&|TvJ0`dY7oqkIHslFMCf1?+K)e6M+mA3rfPergf`kjsP8S_`4~ujLT`F zI_-M54Fq-Pgciz}hAvS{H08vY)PTD9&cK^&Y%()7etq>7u^NYgsRRM2wPgrs+2m>1 zfLe+ZH5jy$Q`t`k+NNBa;1xhrGil#Z>t#F!YvX{0_~+7i$M1q{^m|mEbKF z1kAq@^b>qG^*jQl0{D2R=xjG-W=uNahE-uILQDWnpx->~qu^XaJmV%(*ER(oj1zF) zNB8O$c8OGf(E&;5_%4zJ-ya)rX|?If1-Zr}#C-Gr`-_};)?FdSpBW|B1Q8`s50ziR zTTqPw%HCSDAIkRbPWB*@ZLpgC<=}!-yoC-FGuPDeU~b{ddN;x6D}s3}_VI6Q$wdKp z!)fpF>e>>%6sqa>84dwr%7s* z3pIbAFlSyFBYmN~SGbvy`+!OvR#qXgC-E1C$kPY@QGQMGV6z{w{+wXn{mF{h$)YB* zbFL|!FgszSS6aERHOG@+PHS~x*EGicbjU?^li{|-Kc-x_D^B94LR#@o!*2rD{@}G%}scg4B%?|pPg*m^FRQzR2X8UvoMh3pf zd3S!(5`QZSTduHHWh{J&%wSi0u%EG6S01s5KiS}IAwnWM=pL(vCtENQbYB(>4Pi4n z=Ld-RpJo#2A$+fMb^jWC^GeBw_0-JMnXd$j>U0zc9Sq z2Cww0K1q3Dr4py-*S{gSNXx9%#w^5U%SkI`{J>+|v^dM+vkUP(XBF|!V{KfD`bC`)L&0;JdQv=%npqn$5bw+Q*jo zrTLL7Rkob!*IH#?T1)p{dw=EXW!vIDY3<#mxzcU*+D&xI;B`7>@AaLI**EB?G@37| zZT)a+In`79@MSlR@a5qa?b-Qc*Z)Eu{l5OR{!i7{IX^On(Ov$xq@qIjOWD7~DRS~E z(`a;4>su;u*M^zybCzne{~e2ZUU@fvN`JMV-MgG|;>$IK6_U>+y+2Ema!d01eTOA; zpQFdu_CLI8y*cuJ=If1*r%f_(+5)G}9^Vhc)HE;TwBA?0mZ!GuH||TDrshz#aUK$% z%JG(A_fyE~e?6bx)5@zjHoEyxSk6G~Iilh7zD$u*N5_(A-|8<0@kdJ;uuS4$drN`A zCDRRhvYo|mDNFN_m8*0XOtR}tBWF`5PPM%^`gQ}nt(dr1vqmVu>tu(t^|6~B&+a%% z{iTL&ync7VFl8v??B!a^WjsN^jDg>$mY5YU_vk)IUg7e*xuC93sgcyH zCb17M|8rk7&L8>ca@@teV4H2l$6K!hn#v@9Fw*9Bes9gwrY{j^9zV@-Q9D;sR_S?E z_VHo~2CkGN?~X_RF|ys^S-tx?p?SurcA675ob?`a|A^T}TBid#3nXw&*!r(SgY9X7*9I($B(ef6O}Llu}a^Jj$Z+#dE^qF|NWzD?P+bR4?4Q4CsZBZZ{=4< z-2-(rM)vgo!d&>OUORKSa#9^B@vC(P|BLEK zf2}A$eyI%j^)RW^djq4Q^wfsT4TBJ=v(xL8>iuNLOwVKT)0Nt#NlVMnZui75o~_kc zpsM1AhuvNodVKtENGibxRwu6fdW6N^9B$n1F7&WO+rH`@B$}R z33B`zWx;Fn{)>*+q&@8@-AKB7Sh*X0=DqNqBWXE*vi-e8UP&=c*m_`tTv)2yH=k?% zRDaD`)0F?wo~cT1!^{3+(T+!p-&2INeyS6ED5!C*rvqKdi}%bv+HxB z_9n0P)ER8wH($Rh(WkRq%UF47Jaa27K71N$V;jG6O>euXoDs7@>Ey zXvOY0->Q6j{F7CQ>v{{b!8dA_AKUuI`@%BbT&w-1^CK54d#HFi%pbJ5|6nO@hdWun z^55LR$QR#L8s;;ItI*jj7*ZCnY=q%14=VAqKN;10eVGX-sQ{hx{<~f<-7fL*+;ZqH z?x^+m!6k2;O6&CVNQ0+cBR9}OTdKFOtYC%Ryqfo^##qv<<#C!DEI0#LRmqFMRRx8and7J)iPR zlVCj!?=o@D9a%Wj0@vM1+vUsXWol0cZ_gqpRpw9Y=*G9gG98Bl`_E)fPe{_!PmNgr zM}2pdn1axpzEO=*fnuW1uMg24WN*7QafE9X{6$+FHW(aXkg!i5(Y}uL`S>Zul?F`d z%ud^9`>$>MHr6UqvAkouEq~*a?Dc1e7c|$WUL6|W+P0)h{($D%xCZ{?pU0Z?aQvye zmkN(`Y=*_4HEdh!9?9-6i{Z2NUl%A$dtCO_tgIx{u050a>F zZ6_!X&hK2K^W>IQz(MN%C)&4>x2yWi&wsqWUwfnG!OV`Cgo|SF(Waii#bk7>SmDe` zb0N9LGE#fSx=swexvCOXM;Lo)szJJXEYi{2sLfwkt=P8fR!Q9MJ;C?!I?s>vZ#3N& zz3rqg+vl&IiyHmQ*5gYvnS)D0za_bEi4XeIVjtWmTWUdA=c9dRz;Fu-k3oF9?di1J zV5^?P`u7h?^tNkdU2}V;LY&&P(Q|ON^b~1Lx?gnE=1*oN-75%9nRyO}=Xass(z|An zRcpUkZx!i;)VJoxZ=F!tm9}5VkxBH8#BXff{}WC2_3QH_66U$!P7r7z#(Y><3SW76Px*Bz_&Q*Z8DNUlG( zE`Q+tz%25Qru~!Vn5uIp4AH0G&G?4>fz2LhY}d5<4lh4AG_zodpkr1ZpD{P6HzZH) zr)v~xUq0IQLb}7#EV?_~wlQ0EI^d07q4NE@pN{g~Un|7Uz@z z7xl`4wbgXTR?isNp?kSJ_?qO{YDi@x&YctGxTm%z~8tBX#C|L-LP2HcKB-d68;c9O^MQHae1tLE*U} zMT1$t-1^AIo|z9~m+vxg&-x$1zpZxdyalX$9TImy{PYe}MM+8Iyb zu6e~Zx%>YvCf(A2zX?7!qdjrr(qD7OkL%@CEvvkaE0ZTuSC_`($Gn2#-D8o-cw~UG zq4SWA`;K#ecf@tI+%$LO8<1@GT17-cYwZlG`&TE)27|~!t6E{lPi>$|ET(Vz}RSn zs@gh9@oM|6EqI@EZ11@9@4{kUonvE)Wq60b6CPjaGF$lU!2XQpuN!rS%CPf6SJ1Ol zLEvg!1%z1YIWf0k<_9vU*E{fGp~+Fb;i_aAcEv~~yWdyVsm4vI-n`^{((tvd{iL;z zhR>xsTNk-82tf^3-9o2`OVjL{|ADqjR(#rZOgvj#wn6P^$FR-UcyU~8?k<&oPHLUn zqwSTukGooi)w1tlCtP&=R|7&cKd(64wKVkX-nK7&|K(h046A`9s+{XEs!v?fpk6u? zur^!wudmx}LJ31{AFskPFWpV~emdiEnts+T!ILaY47Jlv_C)6x<98=pK8!E)AVNW7Jn~ z$Wa|`*_g7RJMALKsCxW;xZ}t_C$wshELF{64Go|ByxjI4QLOZ3%{7;gRapXJq z&`wLZtMGVJqeCyfzSFqEQ!rN@_U!VHAH*nC{^)}em%E%+vah6Q#R$SPCT0h|WiA5$ zX?@&UjkRv2sj?}{r;i3-)||R?nY&%m=~2RuYrl9Y=Vvo><5Y64d2i0RO-Ag0TO*d{ zS_+CZPd&L3WsYtA->q+ZZKAZdzq@}&?^upp_xQ|t=T{ql_Lk69-o))sK+GrUuXm;N z1e{7*R&itN?3#g3Ui^r@qcki3fC^pm&gk}}`9Cm~x`y96$@k%(ypK;W*^$>4pYG&z z+-!Z&d;FQl&oj2HB%WNK!=Tr-pS!nC8bPlqI#&*t_Fp@q+OczGeA@$(dD`3bjVEf$ zmX5oqjbGPqgrW^qI?wCb8(s^xe!O$3{fud^e!a1=AZyy-+#P-CzgJ$S$&gEYlZI(3k6}4$I%Nwx z(wyiz8^XxJl)B3+&YrS+Gy3p@5Ow&7K@DlR?Q+;2?KbNnn-k6s>q$=^++OL@{wpZTJL+|>- z8Aq$C8;$HpkISeNU213LgKDQMZ?=7KJtrXkDa`a*`+owC2yyp1`@smMmSMOfzT1c) z{xHk@8$9MK&68=iYrKvdZ3SKNi-v>g!<*54tn_r<~4ZeSm=LX1T ziF_uDgN9)2F!WZ&`(sdL90qG+ydkcbi*FX;K_|>{Mh#D#=t_D>nS>tlx;B0+CKLFntbwu)x$+c+mk(-0=T5!dvHkv z?u$i-Q<#^4&8b42^EfOE7vy2LLj3gzYhGez9Ugdxb?@=WC(Qndd&Ji~5V=VnSC7Uz z6}&kK-PAE`8aC>poF2v)V7W0SEJSGs+!us3!I=CCi(2v2cTDfD_*|r&ByQ-}*_$%| z+u8DXd?0=qEaZ(6j*mliHQYBD9kg+aKK?Ys9kcO?F&Z--55@RlPgIia@;zV<{NZ_Eoug<}|d0y7!{vk zS|i%OM}t-j`T5`eN`Hx5D~XML@V*T0R6zA1*nJ8v*TAVdn5&11#yHmm7g*w?#kh0{ zu33S9uSXvzoU#Rt{Bc+?dhNk6q4+QsV^h&P6WcH0(`=l49p4)kNOlvw=NaQ68yvnC ze>!5S3(oY%)*vk1hyG!heh4cMR{B~i>;ADooL)~ zQaIl!R7=Flizt?Zw%1WBU#L@zx2rL+4&C3POEXUBTCLhquLs&mV!aIN4aGl3q0&UW zKNjAN-^w+pN{5pu)!3gEbvJQ>b}OV z?=kQrYJ9;~IZZQ>j%t{tEo@K6L55f|8!b$6yk%!A7h=3M{@FTh_gt&dM*t2)~0p{glnnBgib%vy-e%h1#wt(-8$1t+@V3r`&6hZ_5Fe<;cw7W$0F&|_GC z9FLsB4+-d&jCz-G!xg-lg|64}OCD}5z;U;6d#~lOB5nHNpUPIfFa;@p#{ zorbF~VoNdJy^EvoxbY#B|J6?hfGCx9o%Us#GB*S zI6N}_*L9I*O8SK&@2lWOO*GNP{`2s}VpO)pRjV;-1IBH}Q68AI1-op+?v>c*HD2s) zRoXG<{ZL*8cMnCak=Qr}J=C#H3zzC+t%Xp>4x{bSJ`uytpy>tlFToY}ap4OgzX2ak z7pv^JuSR&?6r(Nit~L5D#|t4Sdl;i4(JLOO>jc(`Y|_V^xj4fDn-}2L^_b{@jqL4!C3^ivNYuUa07cg}X6+KS~@y zmm+-l7z3)Y@|6%*gSCz5+Jru>*!3Gu=o--8k>3qH`(ug>x+q|~8h$lKw6o$y6mkOr-iKD%7d^;W-ub?kd z*BbY^V%}E#<&TqN@j((^K8NZ8{T)U`zkcQo3LZ97oeAH^fEE*g)Wz|)CX zdIn{mV#^CO82Zawq~a+2M-}&K Date: Wed, 19 Aug 2026 11:53:51 -0400 Subject: [PATCH 5/6] Extract mirror-pair canonicalization into a pure, testable function merge_chrom conflated two responsibilities: assembling the per-chromosome pgen, and canonicalizing mirror-pair indel IDs. Separate them: - Extract canonicalize_mirror_event_ids(chrom, pos, id, ref, alt) as a pure function that computes the mirror-pair event-ID mapping from variant records and returns it (0-row data.frame if no pairs). No file I/O. - do_merge_chrom now calls it, then owns reading/writing/ID-substitution as before. Output is byte-identical (5 pipeline tests unchanged; full SoS pipeline reproduces the same 675-variant / INS:T+DEL:T result). - Guard the CLI with sys.nframe()==0L so the worker is sourceable for unit tests. - Add a direct unit test (test_canonicalize_mirror_event_ids_unit + helper R script) exercising the rule on an in-memory fixture in ~1s: mirror pair relabeled, non-mirror indel and SNP excluded, empty input -> empty mapping. Surfaced and fixed a latent empty-input edge case. --- code/script/reference_data/rss_ld_sketch.R | 132 ++++++++++-------- tests/helpers/canonicalize_mirror_check.R | 37 +++++ .../reference_data/test_rss_ld_sketch.py | 13 ++ 3 files changed, 126 insertions(+), 56 deletions(-) create mode 100644 tests/helpers/canonicalize_mirror_check.R diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index 35a775b47..54c3bae71 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -165,46 +165,21 @@ do_process_block <- function(argv) { } -# ---- 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) - - read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", - stringsAsFactors = FALSE) - 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") - - # Create a directional minimal-event identifier without changing VCF alleles. +# ---- 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)) @@ -239,19 +214,60 @@ do_merge_chrom <- function(argv) { c(event_id = id, event_type = type, event_pos = event_pos, event_ref = rr, event_alt = aa) } - chrom_col <- if ("#CHROM" %in% names(pvar)) "#CHROM" else "CHROM" - event <- t(mapply(canonical_event, pvar[[chrom_col]], pvar$POS, pvar$ID, - pvar$REF, pvar$ALT, SIMPLIFY = TRUE)) - event_map <- data.frame(ID = pvar$ID, CHROM = pvar[[chrom_col]], - POS = pvar$POS, REF = pvar$REF, ALT = pvar$ALT, + 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 <- event_map[sk %in% nk, , drop = FALSE] + event_map[sk %in% nk, , drop = FALSE] +} + +# ---- 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) + + read_tab <- function(f) read.delim(f, check.names = FALSE, comment.char = "", + stringsAsFactors = FALSE) + 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") + + # 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. @@ -346,14 +362,18 @@ p <- add_argument(p, "--sample-list", help = "process_block: optional sample sub p <- add_argument(p, "--output", help = "generate_w: output W .rds") 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 = "") -argv <- parse_args(p) +# 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'") + 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/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 dd4cbc723..a28192520 100644 --- a/tests/notebooks/reference_data/test_rss_ld_sketch.py +++ b/tests/notebooks/reference_data/test_rss_ld_sketch.py @@ -132,3 +132,16 @@ def test_rss_merge_chrom_mirror_event_ids(run_r, run_sos, repo_root, tmp_path): 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 From 488934d2a871d84418acf69891332e30c56ca566 Mon Sep 17 00:00:00 2001 From: Anak Empawi Date: Wed, 19 Aug 2026 12:05:13 -0400 Subject: [PATCH 6/6] Modularize do_merge_chrom into a short orchestrator Follow-through on the canonicalization extraction: pull the remaining inline concerns out of do_merge_chrom so each step is a named, focused function and the merge step reads as a sequence of calls. - reconcile_afreq(): gather per-block .afreq, validate against .pvar, reorder, atomically install; returns pvar. - apply_event_id_substitution(): write the .event_id.tsv sidecar and substitute event IDs into .pvar/.afreq (collision- and duplicate-guarded, atomic). - summarize_block_filters(): tally per-block .meta counts and print the summary. - cleanup_merge_intermediates(): remove plink2 merge temps and block dirs. - Hoist read_tab() to a file-level helper shared by the above. do_merge_chrom drops from ~113 lines to a ~25-line orchestrator: reconcile -> canonicalize -> (substitute | note-standard) -> summarize -> cleanup. Behavior byte-identical: 6 tests pass and the full SoS pipeline reproduces the same 675-variant / INS:T+DEL:T result with identical reconcile/substitution/ filter-summary messages. --- code/script/reference_data/rss_ld_sketch.R | 165 ++++++++++++--------- 1 file changed, 93 insertions(+), 72 deletions(-) diff --git a/code/script/reference_data/rss_ld_sketch.R b/code/script/reference_data/rss_ld_sketch.R index 54c3bae71..117b15253 100644 --- a/code/script/reference_data/rss_ld_sketch.R +++ b/code/script/reference_data/rss_ld_sketch.R @@ -225,18 +225,15 @@ canonicalize_mirror_event_ids <- function(chrom, pos, id, ref, alt) { event_map[sk %in% nk, , drop = FALSE] } -# ---- 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) +# ---- 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) - 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)) @@ -263,6 +260,88 @@ do_merge_chrom <- function(argv) { 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" @@ -272,72 +351,14 @@ do_merge_chrom <- function(argv) { # 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) { - 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") - + 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-level filtering before removing intermediates. - meta_files <- list.files(chrom_dir, pattern = "[.]meta$", recursive = TRUE, - full.names = TRUE) - if (length(meta_files)) { - 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(" -=== Filter Summary for ", basename(chrom_dir), " === -", sep = "") - print(data.frame(value = unlist(summary), row.names = names(summary))) - } - - # Cleanup happens only after all final outputs and summaries are complete. - 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) + summarize_block_filters(chrom_dir) + cleanup_merge_intermediates(final_prefix, block_dirs) } # ---- CLI -------------------------------------------------------------------