Independent list of variants using LD clumping

Prunes a reference genotype panel down to a set of approximately independent variants, one LD block at a time.

Overview¶

Many downstream methods assume the variants they are handed are approximately independent. A reference genotype panel is not: neighbouring variants are correlated through linkage disequilibrium, so a raw variant list counts the same signal several times over. This step runs PLINK LD clumping (--indep-pairwise) within each LD block and then merges the survivors into a single list, which is the form mashr and similar analyses expect.

Working one block at a time keeps each clumping job small and lets the blocks run in parallel. The merge step afterwards also rewrites the variant IDs into one consistent format.

Synthetic-data note. The minimal working example runs on a small synthetic chr22 PLINK panel, protocol_example.ld_genotype.chr22, of 60 samples and roughly 18k variants, built from the toy genotype VCF. It is for demonstration only and is not real individual-level data.

When to run it. Once against the reference panel, before any analysis that needs a set of independent variants.

Input¶

--genotype-list: a tab-separated file with no header. Each row is an LD block label in the first column followed by one or more PLINK bfile prefixes. Every prefix needs its .bed, .bim and .fam present together, otherwise the run stops with Error: Failed to open ... .bed. Example tests/fixtures/ld_prune_reference/protocol_example.ld_genotype.list:
```
chr22_10516173_17414263  tests/fixtures/ld_prune_reference/genotype/protocol_example.ld_genotype.chr22.bed
chr22_17414263_19608021  tests/fixtures/ld_prune_reference/genotype/protocol_example.ld_genotype.chr22.bed
chr22_19608021_22004675  tests/fixtures/ld_prune_reference/genotype/protocol_example.ld_genotype.chr22.bed
```
--window, --shift and --r2: the PLINK --indep-pairwise parameters, 50, 10 and 0.001 by default. An --r2 this strict is aggressive on a small panel; if a block comes back with an empty .prune.in, loosen --r2 to something like 0.2 or widen --window.
--cwd: the directory outputs are written to.
--modular-script-dir: the code/script directory holding the modular analysis scripts.

Output¶

{cwd}/LD_pruning/{block}.prune.in - one file per LD block, listing the variants clumping kept, one per line. Example tests/fixtures/ld_prune_reference/prune_in/chr22_19608021_22004675.prune.in, 273 lines:
```
chr22:10761227_G_A
chr22:11312810_G_A
chr22:11328856_G_T
chr22:11334723_T_C
```
PLINK writes the discarded variants to {block}.prune.out and a run log to {block}.log beside it.
{cwd}/LD_pruning/LD_pruned_variants.txt - the merged list across every block, with the variant IDs rewritten from PLINK’s underscore separator to colons. Example tests/fixtures/ld_prune_reference/expected/LD_pruned_variants.txt, 820 lines and 5 columns:
```
chrom  pos       alt  ref  variant_id
chr22  10761227  G    A    chr22:10761227:G:A
chr22  11312810  G    A    chr22:11312810:G:A
chr22  11328856  G    T    chr22:11328856:G:T
```

Minimal Working Example¶

LD_pruning runs both numbered steps in order: LD_pruning_1 clumps each block, then LD_pruning_2 merges the survivors. LD_pruning_2 runs in R and needs data.table, dplyr and stringr in the active environment.

Timing: ~5-10 min (on toy dataset)

sos run pipeline/ld_prune_reference.ipynb LD_pruning \
    --genotype-list tests/fixtures/ld_prune_reference/protocol_example.ld_genotype.list \
    --cwd output/ld_pruned

Command Interface¶

sos run pipeline/ld_prune_reference.ipynb -h

usage: sos run code/SoS/reference_data/ld_prune_reference.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:
  LD_pruning

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 1h
  --mem 80G
  --numThreads 10 (as int)

Sections
  LD_pruning_1:
    Workflow Options:
      --genotype-list VAL (as path, required)
      --window 50 (as int)
      --shift 10 (as int)
      --r2 0.001 (as float)
  LD_pruning_2:

Workflow implementation¶

The [global] section declares parameters (output directory, modular_script_dir pointing at code/script, and resource options). [LD_pruning_1] reads the genotype list and runs PLINK --indep-pairwise per block; [LD_pruning_2] merges all *.prune.in files and reformats the variant IDs into a single table via the reference_data/ld_prune_reference.R worker.

[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 = "1h"
parameter: mem = "80G"
parameter: numThreads = 10

[LD_pruning_1]
parameter: genotype_list = path
parameter: window = 50
parameter: shift = 10
parameter: r2 = 1e-3

import csv
# genotype_list: header row + one row per block (block_id, bfile...); pandas read
# used header=0, so the first row is treated as a header here too.
with open(genotype_list) as _fh:
    _rows = [r for r in csv.reader(_fh, delimiter="\t") if r]
_rows = _rows[1:]
input_blocks = [r[0] for r in _rows]
input_files = [r[1:] for r in _rows]
del _rows, _fh

input: input_files, group_by = 1, group_with = "input_blocks"
output: prune = f'{cwd}/{step_name[:-2]}/{_input_blocks}.prune.in'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    plink2 \
        --bfile  ${_input:anr} \
        --indep-pairwise  ${window} ${shift} ${r2} \
        --rm-dup force-first \
        --out ${_output["prune"]:annr} \
        --threads ${numThreads}

[LD_pruning_2]
input: group_by = "all"
output: f'{cwd}/{step_name[:-2]}/LD_pruned_variants.txt'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/reference_data/ld_prune_reference.R \
        --input ${_input:r,} \
        --output ${_output:r}