RSS Fine-mapping with GWAS Summary Statistics

RSS Fine-mapping with GWAS Summary Statistics#

End-to-end SuSiE-RSS fine-mapping over GWAS summary statistics, one fit per study per LD block.

Overview#

Fine-mapping asks which variants in a region are causal, not merely correlated with a causal one. Doing that properly needs individual-level genotypes, which for GWAS are usually unavailable. SuSiE-RSS fits the same sum-of-single-effects model to what can be obtained – per-variant z-scores plus an LD matrix from a reference panel – and returns credible sets and posterior inclusion probabilities rather than a single lead SNP per locus (Zou et al., 2022).

Because the LD reference stands in for the study own genotypes, results are only as good as the match between the two. Allele flips and z-score outliers relative to the panel are the usual source of spurious credible sets, which is why the QC options exist: SLALOM or DENTIST screening flags suspicious variants, and RAISS can impute variants that are missing from the study but present in the panel. The LD reference must be genotype-backed (PLINK2/PLINK1/GDS/VCF); a precomputed .cor.xz correlation matrix is not supported here, because the QC steps need the genotypes themselves.

The module takes a list of LD reference blocks and a list of GWAS summary statistics files and runs four SoS steps, which chain in the order below:

Workflow

What it does

generate_manifest

Pair each study with the LD blocks it covers

generate_gwas_sumstats

Harmonise, QC and optionally impute summary statistics per block

gwas_fine_mapping

Fit SuSiE-RSS per block

gwas_rss_plot

Diagnostic plots for the fitted regions

generate_manifest resolves --gwas-meta or --gwas-tsv-list against --region-list or --regions into one manifest TSV with a row per study x region. generate_gwas_sumstats then builds one QC-d GwasSumStats RDS per study x LD block, harmonising variants to the LD-panel ID convention. gwas_fine_mapping dispatches to pecotmr::fineMappingPipeline(), and gwas_rss_plot renders a PIP plot unless --no-plot is given. The results feed TWAS and PRS weight construction downstream.

This chain replaces the legacy two-step get_analysis_regions + univariate_rss workflow, which called the now-removed rss_analysis_pipeline().

When to run it. After GWAS summary statistics and a genotype-backed LD reference panel are in place, and before TWAS or PRS weight construction, which read the GwasFineMappingResult objects written here. This is the summary-statistics-only route: with individual-level genotypes and phenotypes, use mnm_regression instead.

Input#

GWAS inputs – supply either form:

  • --gwas-meta input/rss_analysis/protocol_example.rss_mwe.gwas_meta.tsv TSV with columns study_id, path (relative to the meta file, or absolute) and an optional column_mapping naming a per-study YAML. One row per study.

    study_id	path	column_mapping
AD_Bellenguez_2022	protocol_example.gwas_sumstats.chr22.tsv.gz	protocol_example.gwas_column_mapping.fixed.yml

    In this file path resolves to input/rss_analysis/protocol_example.gwas_sumstats.chr22.tsv.gz and column_mapping to input/rss_analysis/protocol_example.gwas_column_mapping.fixed.yml.

  • --gwas-tsv-list S1=path1 ...: explicit STUDY=PATH pairs; no per-study column mapping is possible this way.

Each summary statistics file is tab-delimited, with chrom, pos, A1, A2 as the first four columns. A column-mapping YAML lets other layouts be read: it must resolve chrom, pos, A1, A2, either z or both beta and se, and n_sample or the pair n_case / n_control (for a case-control study, give one or the other and fill the unused field with 0). var_y is optional.

Region inputs – supply either form:

  • --region-list: BED-like TSV of LD blocks.

  • --regions chr:start-end ...: explicit regions.

LD reference and genome:

  • --ld-meta (required): LD-reference meta TSV with columns #chr, start, end, path, where path is a genotype prefix. One LD block per row, or a single row with start=end=0 for a whole chromosome. The reference must be genotype-backed (PLINK2/PLINK1/GDS/VCF). The protocol ships input/ld_reference/protocol_example.ld_meta_file.tsv, whose prefixes point into input/ld_reference/chr*/:

    #chr	start	end	path
chr1	16103	2888443	chr1/protocol_example.LD.chr1
chr1	2888443	4320284	chr1/protocol_example.LD.chr1

  • --genome (default GRCh38): genome build label carried into the output objects.

QC knobs, forwarded to gwas_sumstats_construct.R and on to summaryStatsQc():

  • --qc-method (default none): none, slalom or dentist.

  • --impute / --no-impute (default off): enable RAISS imputation of variants present in the panel but missing from the study.

  • --qc-args: JSON passthrough for any other summaryStatsQc() argument, e.g. {"mafCutoff":0.01,"nCutoff":0,"alleleFlipKriging":true}.

  • --maf (default 0.0025; 0 disables): MAF cutoff, i.e. mafCutoff.

  • --allele-flip-kriging / --no-allele-flip-kriging (default off) and --effective-n / --no-effective-n (default on): individual summaryStatsQc() switches.

  • --skip-regions chr:start-end ...: drop variants overlapping these windows, i.e. skipRegion.

  • --skip-analysis-pip-cutoff (default 0.025; 0 disables): skip a region whose single-effect max PIP falls below this, i.e. pipCutoffToSkip.

Fine-mapping knobs, forwarded to fine_mapping.R:

  • --methods (default susie), --coverage (0.95), --secondary-coverage (0.7,0.5).

  • --min-abs-corr (0.8), --median-abs-corr (off; sets OR-logic purity), --pip-cutoff (0.025).

  • --L (20) and --L-greedy (5): SuSiE fit defaults.

  • --method-args: JSON overriding any of the above per method, e.g. {"susie":{"L":10}}.

  • --ser-fallback / --no-ser-fallback (default on), --r-finite, --r-mismatch (eb_mix), --rss-control: numerical safeguards for the RSS likelihood.

Other:

  • --no-plot: skip gwas_rss_plot.

  • --manifest-name (default gwas_rss): output prefix for the manifest file.

  • --cwd (default output), --modular-script-dir (default code/script), --container, --job-size, --walltime, --mem, --numThreads.

Output#

  • gwas_rss.manifest.tsv – one row per study x region, recording which summary statistics file and column mapping were paired with which LD block. Every later step reads it; the gwas_rss prefix comes from --manifest-name. The example command produces a single row, since it runs one study over one LD block, and gwas_tsv / column_mapping hold absolute paths, abbreviated here:

    study_id	gwas_tsv	column_mapping	chr	start	end	region_id	gwas_tsv_basename
AD_Bellenguez_2022	.../input/rss_analysis/protocol_example.gwas_sumstats.chr22.tsv.gz	.../input/rss_analysis/protocol_example.gwas_column_mapping.fixed.yml	chr22	49355984	50799822	chr22_49355984_50799822	protocol_example.gwas_sumstats.chr22

  • sumstats/<study>.<region>.gwas_sumstats.rds – a GwasSumStats object per study x block: the QC-d, optionally imputed summary statistics, with variant IDs harmonised to the LD panel. The step reports how many variants survived into the object:

    Wrote QC'd GwasSumStats for study 'AD_Bellenguez_2022' over chr22:49355984-50799822 (7609 variants in)

  • fine_mapping/<study>.<region>.gwas_finemap.rds – a GwasFineMappingResult object per study x block, holding the SuSiE-RSS fit under MethodNames (susie by default): credible sets, PIPs, single-effect posterior summaries, and the QC provenance of the variants that went in:

    Wrote fineMapping result for GwasSumStats (1 row(s))

  • plots/<study>.<region>.pip_plot.png – diagnostic PIP plot, unless --no-plot was given.

  • a .stdout and a .stderr log beside each of the targets above.

Everything is written under --cwd. <study> is the study_id from the GWAS meta file and <region> is the LD block written as chrN_start_end, which for the example command above are AD_Bellenguez_2022 and chr22_49355984_50799822 – one file of each kind.

Minimal Working Example#

Default fine-mapping, no summary-statistics QC#

This fine-maps one study x one LD block and writes the manifest, the QC-d GwasSumStats, the GwasFineMappingResult, and a PIP plot under output/rss_analysis/.

Timing: ~2 min (on toy dataset)

sos run pipeline/rss_analysis.ipynb \
    generate_manifest+generate_gwas_sumstats+gwas_fine_mapping+gwas_rss_plot \
    --cwd output/rss_analysis --modular-script-dir code/script \
    --gwas-meta input/rss_analysis/protocol_example.rss_mwe.gwas_meta.tsv \
    --regions chr22:49355984-50799822 \
    --ld-meta input/ld_reference/protocol_example.ld_meta_file.tsv

SLALOM QC, RAISS imputation and per-method tuning#

Timing: ~2 min (on toy dataset)

sos run pipeline/rss_analysis.ipynb \
    generate_manifest+generate_gwas_sumstats+gwas_fine_mapping+gwas_rss_plot \
    --cwd output/rss_analysis --modular-script-dir code/script \
    --gwas-meta input/rss_analysis/protocol_example.rss_mwe.gwas_meta.tsv \
    --regions chr22:49355984-50799822 \
    --ld-meta input/ld_reference/protocol_example.ld_meta_file.tsv \
    --qc-method slalom --impute --qc-args '{"mafCutoff":0.01}' \
    --min-abs-corr 0.5 --method-args '{"susie":{"L":10}}'

Command Interface#

sos run pipeline/rss_analysis.ipynb -h

usage: sos run pipeline/rss_analysis.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_manifest
  generate_gwas_sumstats
  gwas_fine_mapping
  gwas_rss_plot

Global Workflow Options:
  --cwd output (as path)
  --modular-script-dir code/script (as path)
  --gwas-meta . (as path)
                        --- per-study GWAS inputs (use either or both)
                        ---------------------
  --gwas-tsv-list  (as list)
  --region-list . (as path)
                        --- region inputs (use either or both)
                        -----------------------------
  --regions  (as list)
  --ld-meta VAL (as path, required)
                        --- LD reference + genome
                        ------------------------------------------
  --genome GRCh38
  --qc-method none
                        --- QC knobs (forwarded to gwas_sumstats_construct.R)
                        --------------
  --[no-]impute (default to False)
  --qc-args ''
  --maf 0.0025 (as float)
  --skip-regions  (as list)
  --skip-analysis-pip-cutoff 0.025 (as float)
  --[no-]allele-flip-kriging (default to False)
  --[no-]effective-n (default to True)
  --methods susie
                        --- fine-mapping knobs (forwarded to fine_mapping.R)
                        ---------------
  --coverage 0.95 (as float)
  --secondary-coverage '0.7,0.5'
  --min-abs-corr 0.8 (as float)
  --median-abs-corr ''
  --pip-cutoff 0.025 (as float)
  --L 20 (as int)
  --L-greedy 5 (as int)
  --method-args ''
  --[no-]ser-fallback (default to True)
                        GWAS SuSiE-RSS fine-mapping: SER fallback is ON by
                        default (reproduces the single-panel notebook): fit
                        finite-sample R + EB LD-mismatch SuSiE-RSS and fall back
                        to the single-effect (SER) result for regions susieR
                        flags as unreliable. ser_fallback and r_mismatch are
                        orthogonal: ser_fallback = False only disables the SER
                        fallback (EB LD-mismatch stays on); pass r_mismatch =
                        'none' too for the pre-change plain multi-effect fit.
                        Requires susieR >= 0.16.6 (which moved R_mismatch_method
                        / check_prior into susie_rss_control(); pass those and
                        other control settings via rss_control).
  --r-finite ''
  --r-mismatch 'eb_mix'
  --rss-control ''
  --[no-]no-plot (default to False)
                        --- plot opt-out
                        ---------------------------------------------------
  --manifest-name 'gwas_rss'
                        --- output prefix for the manifest file
                        ----------------------------
  --container ''
                        --- infrastructure
                        -------------------------------------------------
  --job-size 1 (as int)
  --walltime 2h
  --mem 16G
  --numThreads 1 (as int)

Sections
  generate_manifest:
  generate_gwas_sumstats:
  gwas_fine_mapping:
  gwas_rss_plot:

Workflow implementation#

[global]
parameter: cwd = path('output')
parameter: modular_script_dir = path('code/script')
# --- per-study GWAS inputs (use either or both) ---------------------
parameter: gwas_meta = path('.')
parameter: gwas_tsv_list = []   # list of STUDY=PATH items
# --- region inputs (use either or both) -----------------------------
parameter: region_list = path('.')
parameter: regions = []         # list of chr:start-end strings
# --- LD reference + genome ------------------------------------------
parameter: ld_meta = path
parameter: genome = 'GRCh38'
# --- QC knobs (forwarded to gwas_sumstats_construct.R) --------------
parameter: qc_method = 'none'   # none | slalom | dentist
parameter: impute = False
parameter: qc_args = ''         # JSON object spliced into summaryStatsQc()
parameter: maf = 0.0025         # MAF cutoff (summaryStatsQc mafCutoff); 0 disables
parameter: skip_regions = []    # chr:start-end window(s) whose variants are dropped (skipRegion)
parameter: skip_analysis_pip_cutoff = 0.025   # skip a region below this max PIP (pipCutoffToSkip); 0 disables
parameter: allele_flip_kriging = False   # kriging allele-flip QC: sign-flip switched z (logLR>2 & |z|>2); off by default
parameter: effective_n = True   # case/control: use effective sample size 4/(1/n_case+1/n_control) as N (summaryStatsQc effectiveN); False -> raw N
# --- fine-mapping knobs (forwarded to fine_mapping.R) ---------------
parameter: methods = 'susie'
parameter: coverage = 0.95
parameter: secondary_coverage = '0.7,0.5'
parameter: min_abs_corr = 0.8
parameter: median_abs_corr = ''   # empty -> off (OR-logic purity; needs pecotmr step-1)
parameter: pip_cutoff = 0.025
parameter: L = 20
parameter: L_greedy = 5
parameter: method_args = ''     # nested per-method JSON object
# GWAS SuSiE-RSS fine-mapping: SER fallback is ON by default (reproduces the single-panel notebook):
# fit finite-sample R + EB LD-mismatch SuSiE-RSS and fall back to the single-effect (SER) result for
# regions susieR flags as unreliable. ser_fallback and r_mismatch are orthogonal: ser_fallback = False
# only disables the SER fallback (EB LD-mismatch stays on); pass r_mismatch = 'none' too for the
# pre-change plain multi-effect fit. Requires susieR >= 0.16.6 (which moved R_mismatch_method /
# check_prior into susie_rss_control(); pass those and other control settings via rss_control).
parameter: ser_fallback = True    # fall back to single-effect (SER) when susieR flags R unreliable
parameter: r_finite = ''          # finite-sample R size; empty -> auto (LD-panel sample size)
parameter: r_mismatch = 'eb_mix'  # LD-mismatch correction: 'none' | 'eb' | 'eb_mix' (residual-mixture EB, default)
parameter: rss_control = ''       # JSON object of susie_rss_control() settings, e.g. '{"check_prior":true,"mismatch_estimator":"map"}'; empty -> defaults
# --- plot opt-out ---------------------------------------------------
parameter: no_plot = False
# --- output prefix for the manifest file ----------------------------
parameter: manifest_name = 'gwas_rss'
# --- infrastructure -------------------------------------------------
parameter: container = ''
parameter: job_size = 1
parameter: walltime = '2h'
parameter: mem = '16G'
parameter: numThreads = 1

[generate_manifest]
# Resolve --gwas-meta / --gwas-tsv-list x --region-list / --regions into a
# single manifest TSV via gwas_rss_manifest.R. Downstream steps fan out
# over its rows; no Python parsing in this notebook.
input: None
output: f"{cwd}/{manifest_name}.manifest.tsv"
task: trunk_workers = 1, trunk_size = 1, walltime = '15m', mem = '2G', cores = 1, tags = f"{step_name}_{_output:bn}"
bash: expand = '${ }', stderr = f"{_output}.stderr", stdout = f"{_output}.stdout", container = container
    Rscript ${modular_script_dir}/pecotmr_integration/gwas_rss_manifest.R \
        ${('--gwas-meta ' + str(gwas_meta)) if gwas_meta.is_file() else ''} \
        ${('--gwas-tsv-list ' + ' '.join(str(x) for x in gwas_tsv_list)) if gwas_tsv_list else ''} \
        ${('--region-list ' + str(region_list)) if region_list.is_file() else ''} \
        ${('--regions ' + ' '.join(str(x) for x in regions)) if regions else ''} \
        --output ${_output}

[generate_gwas_sumstats]
# Fan out over the manifest's rows: one (study, region) GwasSumStats per
# row. Manifest columns: study_id, gwas_tsv, column_mapping, chr, start,
# n_case, n_control, n_sample (optional), end, region_id, gwas_tsv_basename.
# n_sample is the study-level total N: the level-4 fallback used when a study
# has no per-variant n_case/n_control and no per-variant N column.
import csv
jobs = list(csv.DictReader(open(f"{cwd}/{manifest_name}.manifest.tsv"), delimiter='\t'))
input: for_each = 'jobs'
output: f"{cwd}/sumstats/{_jobs['study_id']}.{_jobs['region_id']}.gwas_sumstats.rds"
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}.stderr", stdout = f"{_output}.stdout", container = container
    Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_construct.R \
        --study ${_jobs['study_id']} \
        --gwas-tsv ${_jobs['gwas_tsv']} \
        --ld-block ${_jobs['chr']}:${_jobs['start']}-${_jobs['end']} \
        --ld-meta ${ld_meta} \
        --genome ${genome} \
        ${('--column-mapping ' + _jobs['column_mapping']) if _jobs['column_mapping'] else ''} \
        --qc-method ${qc_method} \
        --maf ${maf} \
        --pip-cutoff-to-skip ${skip_analysis_pip_cutoff} \
        ${('--skip-region ' + ','.join(str(x) for x in skip_regions)) if skip_regions else ''} \
        ${'--impute' if impute else ''} \
        ${('--qc-args ' + repr(qc_args)) if qc_args else ''} \
        ${'--allele-flip-kriging' if allele_flip_kriging else ''} \
        --effective-n ${'TRUE' if effective_n else 'FALSE'} \
        ${('--n-case ' + _jobs['n_case']) if _jobs.get('n_case','') not in ('', 'NA') else ''} \
        ${('--n-control ' + _jobs['n_control']) if _jobs.get('n_control','') not in ('', 'NA') else ''} \
        ${('--n-sample ' + _jobs['n_sample']) if _jobs.get('n_sample','') not in ('', 'NA') else ''} \
        --output ${_output}

[gwas_fine_mapping]
# Per-RDS fan-out: one fine-mapping task per (study, region) GwasSumStats.
output: f"{cwd}/fine_mapping/{_input:bnn}.gwas_finemap.rds", group_by = 1
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}.stderr", stdout = f"{_output}.stdout", container = container
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --gwas-sumstats ${_input} \
        --methods ${methods} \
        --coverage ${coverage} \
        --secondary-coverage ${secondary_coverage} \
        --min-abs-corr ${min_abs_corr} \
        --pip-cutoff ${pip_cutoff} \
        --L ${L} \
        --L-greedy ${L_greedy} \
        ${('--median-abs-corr ' + str(median_abs_corr)) if str(median_abs_corr) != '' else ''} \
        ${('--method-args ' + repr(method_args)) if method_args else ''} \
        --ser-fallback ${'TRUE' if ser_fallback else 'FALSE'} \
        --r-mismatch ${r_mismatch} \
        ${('--r-finite ' + str(r_finite)) if str(r_finite) != '' else ''} \
        ${('--rss-control ' + repr(rss_control)) if rss_control else ''} \
        --output ${_output}

[gwas_rss_plot]
stop_if(no_plot, '--no-plot set; skipping PIP plot step.')
output: f"{cwd}/plots/{_input:bnn}.pip_plot.png", group_by = 1
task: trunk_workers = 1, trunk_size = job_size, walltime = '30m', mem = '4G', cores = 1, tags = f"{step_name}_{_output:bn}"
bash: expand = '${ }', stderr = f"{_output}.stderr", stdout = f"{_output}.stdout", container = container
    Rscript ${modular_script_dir}/pecotmr_integration/gwas_rss_plot.R \
        --input ${_input} \
        --output ${_output}