Builds a compact stochastic sketch of a genotype panel for use as an LD reference in summary-statistic fine-mapping.
Overview¶
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.
Rather than storing the full genotype matrix G (n × p), the sketch computes U = WᵀG (B × p) using a random projection matrix . The approximate LD matrix by the Johnson–Lindenstrauss lemma. G is never stored.
Matrix dimensions:
G : (n × p) — n individuals × p variants
W : (n × B) — projection matrix, generated once per cohort
U : (B × p) — stochastic genotype sample = WᵀG, stored in pgen
: (p × p) — approximate LD matrix, computed on-the-fly by SuSiE-RSS from U
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. 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:GAT and chr22:pos:DEL:GAT); 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.
Summary-statistic fine-mapping needs an LD matrix, and storing the full genotype matrix for a whole-genome panel is expensive. The sketch keeps a random projection of the genotypes instead: small enough to distribute, but sufficient to reconstruct the LD structure that SuSiE-RSS actually uses.
When to run it. Once per reference panel, before any rss_analysis run that needs an
LD reference. The output is reference data, not a per-study result.
Input¶
--ld-block-filetests/fixtures/rss_ld_sketch/protocol_example.ld_blocks.bed(regions to sketch, tab-separated with columnschr,start,endin 0-based half-open coordinates; the toy file holds 3 chr22 blocks)
#chr start end
chr22 16000000 20000000
chr22 30000000 34000000--vcf-baseinput/rss_ld_sketchtogether with--vcf-prefix protocol_example.genotype.(directory of bgzipped, tabix-indexed VCFs named{vcf_prefix}{chr}.*.bgzor{vcf_prefix}{chr}.*.vcf.gz; the toy cohort isprotocol_example.genotype.chr22.vcf.gz, 60 individuals)
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE_001 ...--n-samples 60(number of individuals in the VCF. Must match the VCF sample count, orprocess_blockfails with a W shape mismatch.)--B 50(number of sketch (pseudo-)samples to project down to. Sets the second dimension of W.)--cwd output/rss_ld_sketch(working directory for logs and job files. Defaults tooutput.)--output-dir output/rss_ld_sketch(directory the sketch products are written to)--seed 999(random seed for the projection matrix, so a sketch can be reproduced)--W-matrixoutput/rss_ld_sketch/W_B50.rds(the shared projection matrix, passed toprocess_blockaftergenerate_Whas written it)--chrom chr22(chromosome to process)--cohort-id protocol_example(tag used in the per-block and merged output filenames)
Output¶
output/rss_ld_sketch/W_B50.rds(shared projection matrix written once bygenerate_W, shape n_samples x B;process_blockreads it back for every LD block so all blocks share one projection)
R matrix, double (60, 50) # saveRDS/readRDS; one row per sample, one column per sketch dimoutput/rss_ld_sketch/chr22/chr22_<start>_<end>/<cohort_id>.chr22_<start>_<end>.dosage.gz(process_block) (the per-LD-block sketch: B pseudo-samples x the variants in that block, gzipped dosage text. One directory per block, written beforemerge_chromconcatenates them.)output/rss_ld_sketch/chr22/protocol_example.chr22.pgen(with.pvar,.psam,.afreq) (per-chromosome PLINK2 genotype sketch assembled bymerge_chrom: B pseudo-samples x p variants. The.pgenis binary; the three companion files describe it.)output/rss_ld_sketch/chr22/protocol_example.chr22.pvar(variant information)
#CHROM POS ID REF ALT
22 16073625 chr22:16073625:G:T G T
22 16102937 chr22:16102937:C:T C Toutput/rss_ld_sketch/chr22/protocol_example.chr22.psam(the B sketch pseudo-samples, namedS1..S{B})
#FID IID SEX
S1 S1 NA
S2 S2 NAoutput/rss_ld_sketch/chr22/protocol_example.chr22.afreq(allele frequencies and the sketch value range per variant)
#CHROM ID REF ALT ALT_FREQS OBS_CT U_MIN U_MAX
chr22 chr22:16073625:G:T G T 0.108333 120 -2.588204 2.248544
chr22 chr22:16102937:C:T C T 0.308333 120 -1.849294 2.069429output/rss_ld_sketch/chr22/protocol_example.chr22.event_id.tsv(mirror-pair indel mapping, written bymerge_chromonly when such pairs exist; records the original ID alongside the canonical event ID)
ID CHROM POS REF ALT event_id event_type event_pos event_ref event_alt
chr22:17250000:C:CGAT 22 17250000 C CGAT chr22:17250000:INS:GAT INS 17250000 GAT
chr22:17250000:CGAT:C 22 17250000 CGAT C chr22:17250001:DEL:GAT DEL 17250001 GATOnly 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 three-base event at chr22:17250000 represented as both an insertion (C→CGAT) and a deletion (CGAT→C), each assigned a canonical event ID. Read the two rows as one event seen from opposite sides: the first is an insertion of GAT after 17250000, so the event keeps the anchor position; the second is a deletion of GAT beginning at 17250001, so the event position is the anchor plus one. The same rule covers a single-base pair at a different site: at chr22:16600000, A→AT and AT→A become chr22:16600000:INS:T and chr22:16600001:DEL:T — an insertion of T after 16600000, and a deletion of T beginning at 16600001. SNPs, equal-length substitutions, symbolic alleles, and non-mirror indels are omitted and retain their standard IDs.
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).
Minimal Working Example¶
Run the three workflows in order. generate_W builds the shared projection matrix once; process_block sketches each LD block; merge_chrom assembles the per-chromosome pgen.
Timing: ~10-20 min (on toy dataset)
Step 1. Generate the projection matrix W (run once per cohort; --n-samples must equal the VCF sample count).¶
Timing: TBD (on toy dataset)
sos run pipeline/rss_ld_sketch.ipynb generate_W \
--n-samples 60 \
--output-dir output/rss_ld_sketch \
--B 50 \
--seed 123 \
--cwd output/rss_ld_sketch<path>/.pixi/envs/python/lib/python3.12/site-packages/sos/targets.py:22: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
INFO: Running generate_W:
INFO: generate_W is completed.
INFO: generate_W output: output/rss_ld_sketch/W_B50.rds
INFO: Workflow generate_W (ID=w68f63c60d8da4b5e) is executed successfully with 1 completed step.
Step 2. Process all LD blocks for the chromosome — read the VCF, filter variants, and write per-block dosage sketches U = WᵀG.¶
Timing: TBD (on toy dataset)
sos run pipeline/rss_ld_sketch.ipynb process_block \
--ld-block-file tests/fixtures/rss_ld_sketch/protocol_example.ld_blocks.bed \
--chrom 22 \
--vcf-base output/rss_ld_sketch/rss_ld_sketch \
--vcf-prefix protocol_example.genotype. \
--output-dir output/rss_ld_sketch \
--W-matrix output/rss_ld_sketch/W_B50.rds \
--B 50 \
--cohort-id protocol_example \
--cwd output/rss_ld_sketch<path>/.pixi/envs/python/lib/python3.12/site-packages/sos/targets.py:22: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
INFO: Running process_block:
3 LD blocks queued
INFO: process_block (index=0) is completed.
INFO: process_block (index=1) is completed.
INFO: process_block (index=2) is completed.
INFO: process_block output: output/rss_ld_sketch/chr22/chr22_16000000_20000000/protocol_example..chr22_16000000_20000000.dosage.gz output/rss_ld_sketch/chr22/chr22_30000000_34000000/protocol_example..chr22_30000000_34000000.dosage.gz... (3 items in 3 groups)
INFO: Workflow process_block (ID=w8c1e759d62203ef6) is executed successfully with 1 completed step and 3 completed substeps.
Step 3. Merge the per-block dosage sketches, reconcile allele frequencies, and relabel mirror-pair indels with canonical event IDs.¶
Timing: TBD (on toy dataset)
sos run pipeline/rss_ld_sketch.ipynb merge_chrom \
--output-dir output/rss_ld_sketch \
--cohort-id protocol_example \
--chrom 22 \
--cwd output/rss_ld_sketch<path>/.pixi/envs/python/lib/python3.12/site-packages/sos/targets.py:22: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
INFO: Running merge_chrom:
PLINK v2.0.0-a.6.9LM 64-bit Intel (29 Jan 2025) cog-genomics.org/plink/2.0/
(C) 2005-2025 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to output/rss_ld_sketch/chr22/protocol_example..chr22.log.
Options in effect:
--make-pgen
--out output/rss_ld_sketch/chr22/protocol_example..chr22
--pmerge-list output/rss_ld_sketch/chr22/protocol_example..chr22_pmerge_list.txt pfile
--sort-vars
Start time: Tue Jun 23 09:52:54 2026
191527 MiB RAM detected, ~187643 available; reserving 95763 MiB for main
workspace.
Using up to 32 threads (change this with --threads).
--pmerge-list: 3 filesets specified.
--pmerge-list: 50 samples present.
--pmerge-list: Merged .psam written to
output/rss_ld_sketch/chr22/protocol_example..chr22-merge.psam .
--pmerge-list: 3 .pvar files scanned.
Concatenation job detected.
Concatenating... 673/673 variants complete.
Results written to
output/rss_ld_sketch/chr22/protocol_example..chr22-merge.pgen +
output/rss_ld_sketch/chr22/protocol_example..chr22-merge.pvar .
50 samples (0 females, 0 males, 50 ambiguous; 50 founders) loaded from
output/rss_ld_sketch/chr22/protocol_example..chr22-merge.psam.
673 variants loaded from
output/rss_ld_sketch/chr22/protocol_example..chr22-merge.pvar.
Note: No phenotype data present.
Writing output/rss_ld_sketch/chr22/protocol_example..chr22.pvar ... done.
Writing output/rss_ld_sketch/chr22/protocol_example..chr22.psam ... done.
Writing output/rss_ld_sketch/chr22/protocol_example..chr22.pgen ... done.
End time: Tue Jun 23 09:52:54 2026
=== Filter Summary for chr22 ===
value
n_total 6157.0
n_passed 673.0
n_multiallelic 0.0
n_monomorphic 4701.0
n_all_na 0.0
n_high_msng 101.0
n_low_maf 0.0
n_low_mac 682.0
pct_dropped 89.1
INFO: merge_chrom is completed.
INFO: merge_chrom output: output/rss_ld_sketch/chr22/protocol_example..chr22.pgen
INFO: Workflow merge_chrom (ID=w8e5a670551e06660) is executed successfully with 1 completed step.
Command Interface¶
sos run pipeline/rss_ld_sketch.ipynb -husage: sos run pipeline/rss_ld_sketch.ipynb
[workflow_name | -t targets] [options] [workflow_options]
workflow_name: Single or combined workflows defined in this script
targets: One or more targets to generate
options: Single-hyphen sos parameters (see "sos run -h" for details)
workflow_options: Double-hyphen workflow-specific parameters
Workflows:
generate_W
process_block
merge_chrom
Global Workflow Options:
--cwd output (as path)
--modular-script-dir code/script (as path)
Directory holding the modular analysis scripts
(code/script)
--job-size 1 (as int)
--walltime '24:00:00'
--mem 32G
--numThreads 8 (as int)
Sections
generate_W:
Workflow Options:
--n-samples VAL (as int, required)
Generate projection matrix $W \sim N(0, 1/\sqrt{n})$,
shape (n x B). Run ONCE before processing any
chromosome; W depends only on n and B, and all
chromosomes reuse the same W so per-chromosome sketches
are mergeable.
--output-dir VAL (as str, required)
--B 10000 (as int)
--seed 123 (as int)
process_block:
Workflow Options:
--ld-block-file VAL (as str, required)
--chrom 0 (as int)
--vcf-base VAL (as str, required)
--vcf-prefix VAL (as str, required)
--cohort-id 'ADSP.R5.EUR'
--output-dir VAL (as str, required)
--W-matrix VAL (as str, required)
--B 10000 (as int)
--maf-min 0.0005 (as float)
--mac-min 5 (as int)
--msng-min 0.05 (as float)
--sample-list ''
merge_chrom:
Workflow Options:
--chrom 0 (as int)
--output-dir VAL (as str, required)
--cohort-id VAL (as str, required)
--plink2-bin plink2Workflow implementation¶
[global]
parameter: cwd = path("output")
# Directory holding the modular analysis scripts (code/script)
parameter: modular_script_dir = path('code/script')
parameter: job_size = 1
parameter: walltime = "24:00:00"
parameter: mem = "32G"
parameter: numThreads = 8
cwd = path(f'{cwd:a}')
[generate_W]
# Generate projection matrix $W \sim N(0, 1/\sqrt{n})$, shape (n x B).
# Run ONCE before processing any chromosome; W depends only on n and B, and all
# chromosomes reuse the same W so per-chromosome sketches are mergeable.
parameter: n_samples = int
parameter: output_dir = str
parameter: B = 10000
parameter: seed = 123
input: []
output: f'{output_dir}/W_B{B}.rds'
task: trunk_workers = 1, trunk_size = 1, walltime = '00:05:00', mem = '4G', cores = 1
bash: expand = "${ }", stdout = f'{_output:n}.stdout', stderr = f'{_output:n}.stderr'
Rscript ${modular_script_dir}/reference_data/rss_ld_sketch.R \
--step generate_w \
--n-samples ${n_samples} \
--B ${B} \
--seed ${seed} \
--output ${_output}
[process_block]
parameter: ld_block_file = str
parameter: chrom = 0
parameter: vcf_base = str
parameter: vcf_prefix = str
parameter: cohort_id = "ADSP.R5.EUR"
parameter: output_dir = str
parameter: W_matrix = str
parameter: B = 10000
parameter: maf_min = 0.0005
parameter: mac_min = 5
parameter: msng_min = 0.05
parameter: sample_list = ""
# Build the LD-block list from the BED (chr1..chr22; optional single-chrom filter).
blocks = []
with open(ld_block_file) as _fh:
for _line in _fh:
if _line.startswith("#") or not _line.strip():
continue
_p = _line.split()
_c = _p[0]
if not (_c.startswith("chr") and _c[3:].isdigit()):
continue
_cnum = int(_c[3:])
if not (1 <= _cnum <= 22):
continue
if chrom != 0 and _cnum != chrom:
continue
blocks.append({"chr": _c, "start": int(_p[1]), "end": int(_p[2])})
del _fh
stop_if(len(blocks) == 0, msg = f"No blocks found for chrom={chrom} in {ld_block_file}")
input: for_each = "blocks"
output: f'{output_dir}/{_blocks["chr"]}/{_blocks["chr"]}_{_blocks["start"]}_{_blocks["end"]}/{cohort_id}.{_blocks["chr"]}_{_blocks["start"]}_{_blocks["end"]}.dosage.gz'
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads
bash: expand = "${ }"
Rscript ${modular_script_dir}/reference_data/rss_ld_sketch.R \
--step process_block \
--vcf-base ${vcf_base} \
--vcf-prefix ${vcf_prefix} \
--chrom ${_blocks["chr"]} \
--block-start ${_blocks["start"]} \
--block-end ${_blocks["end"]} \
--w-matrix ${W_matrix} \
--cohort-id ${cohort_id} \
--output-dir ${output_dir} \
--B ${B} \
--maf-min ${maf_min} \
--mac-min ${mac_min} \
--msng-min ${msng_min} \
--sample-list "${sample_list}"
[merge_chrom]
parameter: chrom = 0
parameter: output_dir = str
parameter: cohort_id = str
parameter: plink2_bin = "plink2"
import os, glob
if chrom != 0:
chroms = [f"chr{chrom}"]
else:
chroms = sorted(set(
os.path.basename(_d)
for _d in glob.glob(os.path.join(output_dir, "chr*"))
if os.path.isdir(_d)
))
input: for_each = "chroms"
output: f"{output_dir}/{_chroms}/{cohort_id}.{_chroms}.pgen"
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads
bash: expand = "$[ ]"
set -euo pipefail
shopt -s nullglob
chrom_dir="$[output_dir]/$[_chroms]"
final_prefix="${chrom_dir}/$[cohort_id].$[_chroms]"
merge_list="${chrom_dir}/$[cohort_id].$[_chroms]_pmerge_list.txt"
# Step 1: Convert each block dosage.gz -> sorted per-block pgen
> "${merge_list}"
files=("${chrom_dir}"/*/*.dosage.gz)
if [ ${#files[@]} -eq 0 ]; then
echo "No dosage files found in ${chrom_dir}" >&2
exit 1
fi
for dosage_gz in "${files[@]}"; do
block_dir=$(dirname "${dosage_gz}")
block_tag=$(basename "${block_dir}")
prefix="${block_dir}/$[cohort_id].${block_tag}_tmp"
map_file="${block_dir}/$[cohort_id].${block_tag}.map"
psam_file="${block_dir}/$[cohort_id].${block_tag}.psam"
meta_file="${block_dir}/$[cohort_id].${block_tag}.meta"
B=$(grep "^B=" "${meta_file}" | cut -d= -f2)
printf '#FID\tIID\n' > "${psam_file}"
for i in $(seq 1 ${B}); do
printf 'S%d\tS%d\n' ${i} ${i} >> "${psam_file}"
done
$[plink2_bin] \
--import-dosage "${dosage_gz}" format=1 noheader \
--psam "${psam_file}" \
--map "${map_file}" \
--make-pgen \
--out "${prefix}_unsorted" \
--silent
$[plink2_bin] \
--pfile "${prefix}_unsorted" \
--make-pgen \
--sort-vars \
--out "${prefix}" \
--silent
rm -f "${prefix}_unsorted.pgen" "${prefix}_unsorted.pvar" "${prefix}_unsorted.psam"
echo "${prefix}" >> "${merge_list}"
done
# Step 2: Merge all per-block pgens -> one per-chrom pgen
$[plink2_bin] \
--pmerge-list "${merge_list}" pfile \
--make-pgen \
--sort-vars \
--out "${final_prefix}"
# Step 3: reconcile .afreq, relabel mirror-pair indels with event IDs, summarize filters, and cleanup.
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]"
Troubleshooting¶
| Symptom | Cause | Fix |
|---|---|---|
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. |
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. |
No passing variants in chrXX | Filters removed everything (small toy cohort) | Widen --maf-min / --mac-min / --msng-min, or choose blocks with more variants. |
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. |
| Region query returns nothing | Missing tabix index | Run tabix -p vcf file.bgz to create the .tbi. |