RSS LD Sketch Pipeline

RSS LD Sketch Pipeline#

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 \(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.

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

  • \(\hat{R}\) : (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).

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-file input/rss_ld_sketch/protocol_example.ld_blocks.bed (regions to sketch, tab-separated with columns chr, start, end in 0-based half-open coordinates; the toy file holds 3 chr22 blocks)

#chr	start	end
chr22	16000000	20000000
chr22	30000000	34000000

  • --vcf-base input/rss_ld_sketch together with --vcf-prefix protocol_example.genotype. (directory of bgzipped, tabix-indexed VCFs named {vcf_prefix}{chr}.*.bgz; the toy cohort is protocol_example.genotype.chr22.bgz, 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, or process_block fails 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 to output.)

  • --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-matrix output/rss_ld_sketch/W_B50.npy (the shared projection matrix, passed to process_block after generate_W has 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.npy (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)

float32 (60, 50)
[[-0.1402  0.1288  0.0365 -0.1945 -0.0747]
 [-0.1671 -0.1341  0.2251 -0.1030  0.0038]]

  • output/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 before merge_chrom concatenates them.)

  • output/rss_ld_sketch/chr22/protocol_example.chr22.pgen (with .pvar, .psam, .afreq) (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.)

  • 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	T

  • output/rss_ld_sketch/chr22/protocol_example.chr22.psam (the B sketch pseudo-samples, named S1..S{B})

#FID	IID	SEX
S1	S1	NA
S2	S2	NA

  • output/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.069429

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).

The generate_W step declares its output as W_B<B>.rds, but the example on disk is W_B50.npy. The code is authoritative; the example data predates the current naming.

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.npy
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 input/rss_ld_sketch/protocol_example.ld_blocks.bed \
    --chrom 22 \
    --vcf-base input/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 into one per-chromosome PLINK2 pgen.#

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 -h

usage: 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 plink2

Workflow 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}"

    # Remove PLINK merge intermediates immediately after merge
    rm -f "${final_prefix}-merge.pgen" "${final_prefix}-merge.pvar" "${final_prefix}-merge.psam"

    # Step 3: Concatenate .afreq
    first=1
    for f in "${chrom_dir}"/*/*.afreq; do
        if [ "${first}" -eq 1 ]; then
            cat "${f}" > "${final_prefix}.afreq"
            first=0
        else
            tail -n +2 "${f}" >> "${final_prefix}.afreq"
        fi
    done

R: expand = "$[ ]"

    meta_files <- list.files("$[output_dir]/$[_chroms]",
                             pattern = "[.]meta$", recursive = TRUE,
                             full.names = TRUE)
    if (length(meta_files) > 0) {
      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, "=")
        vals <- setNames(as.integer(sapply(kv, `[`, 2)), sapply(kv, `[`, 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 $[_chroms] ===\n")
      print(data.frame(value = unlist(summary), row.names = names(summary)))
    }

bash: expand = "$[ ]"

    # Step 5: Cleanup block intermediates
    chrom_dir="$[output_dir]/$[_chroms]"
    final_prefix="${chrom_dir}/$[cohort_id].$[_chroms]"

    rm -f "${final_prefix}_pmerge_list.txt"
    for block_dir in "${chrom_dir}"/*/; do
        rm -rf "${block_dir}"
    done

Troubleshooting#

Symptom

Cause

Fix

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.

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.