RSS Fine-mapping and TWAS Weights with QTL Summary Statistics

Fine-maps a cis window and learns TWAS weights from cis-QTL summary statistics plus an LD reference panel, without individual-level genotypes.

Overview¶

mnm_regression.ipynb fits a cis window from individual-level data: it needs the genotype matrix and the phenotype for every sample. That is not always possible -- summary statistics are shareable where genotypes are not, and a scan that has already run need not be repeated.

This module takes the other route. It reads a cis-QTL nominal association table (effect, standard error and the variant id per SNP), pairs it with an LD reference panel, and hands the resulting QtlSumStats to the same two pipelines the individual-level route uses:

qtl_rss_fine_mapping runs SuSiE-RSS (susieR::susie_rss) over the window. The regression-with-summary-statistics likelihood replaces the individual-level one; credible sets and PIPs mean what they always did.
qtl_rss_twas_weights learns predictive weights from the same object. Not every method has a summary-statistics implementation: mrash, lasso, scad, mcp, l0learn, mrmash and dpr_gibbs do; enet and the bayes_* family are individual-level only and will be rejected here.

The LD panel is the load-bearing input. RSS reconstructs the joint fit from marginal statistics and LD, so the panel must cover the window’s variants and must be on the same allele orientation. summaryStatsQc() harmonizes the two and reports what it corrected -- read that line. If it says it sign- or strand-flipped everything, the alleles were declared wrong, not fixed (see --variant-id-alleles below).

When to run it. After a cis scan (TensorQTL.ipynb), in place of mnm_regression.ipynb’s susie_twas when only summary statistics are available.

Input¶

--sumstats output/cis/example.cis_qtl.pairs.tsv.gz (the cis-QTL nominal table from TensorQTL.ipynb: one row per (trait, variant) with bhat/sebhat, pvalue, af and n. A z column is used if present, otherwise the Wald z bhat/sebhat is derived.)
--ld-sketch tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed (the LD reference panel: a genotype path/prefix, or a per-chromosome LD meta file. Its variants must cover the window.)
--study / --context / --trait -- the tuple this collection describes. --trait is the gene id, and is also what the trait filter matches.
--trait-column (default molecular_trait_id) -- a cis scan writes every gene into one table, so the rows for --trait are selected by this column. Set it empty for a file that already holds one gene.
--variant-id-alleles (none | A2A1 | A1A2) -- where the alleles come from when the table has no A1/A2 columns, as a cis scan typically does not: they sit inside the variant id. The order cannot be inferred from the string -- a .pvar writes REF:ALT (pecotmr’s canonical A2:A1) while a PLINK .bim writes A1:A2 -- so declare which one your ids use. A real A1/A2 column always wins. Declaring it wrong is silent: QC will “correct” the apparent mismatch by flipping every variant against the panel.
--genome (default GRCh38), --region, --n-sample, --column-mapping -- optional.

QC knobs are forwarded to summaryStatsQc(): --maf, --mac, --imiss, --z-mismatch-qc, --pip-cutoff-to-skip, and --skip-qc for diagnostics.

Output¶

{cwd}/sumstats/{study}.{context}.{trait}.qtl_sumstats.rds -- the QtlSumStats: the window’s variants with SNP/A1/A2/Z/N (plus BETA/SE/P/AF when supplied), the LD panel attached as the ldSketch, and a qcInfo audit of what QC did.
{cwd}/fine_mapping/{...}.qtl_rss_finemap.rds -- a QtlFineMappingResult: credible sets and PIPs from SuSiE-RSS, the same class the individual-level route produces.
{cwd}/twas_weights/{...}.qtl_rss_twas_weights.rds -- a TwasWeights collection.

Each step also writes .stdout / .stderr beside its output. The QC line in the sumstats log is worth reading every time:

[study/context/gene] QC summary: 200 in -> 200 out | corrected: sign-flip 0, strand-flip 0

Minimal Working Example¶

Runs on the committed chr22 toy data: the TensorQTL nominal table for 16 genes, with the same 49-sample genotypes used as the LD reference (in-sample LD -- fine for a smoke test, optimistic for real inference, where a separate reference panel belongs).

sos run pipeline/qtl_rss_analysis.ipynb qtl_rss \
    --cwd output/qtl_rss \
    --sumstats tests/fixtures/tensorqtl/expected/cis_qtl.pairs.tsv.gz \
    --ld-sketch tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --study test_study --context context1 --trait ENSG00000283047 \
    --variant-id-alleles A1A2 \
    --methods susie --twas-methods lasso -j1

Command Interface¶

sos run pipeline/qtl_rss_analysis.ipynb -h

Workflow implementation¶

[global]
parameter: cwd = path('output')
parameter: modular_script_dir = path('code/script')
# --- the (study, context, trait) this run describes ------------------
parameter: study = str
parameter: context = str
parameter: trait = str
# --- inputs ----------------------------------------------------------
parameter: sumstats = path
parameter: ld_sketch = path
parameter: genome = 'GRCh38'
parameter: region = ''
parameter: n_sample = -1.0            # study-level total N; <0 = take it from the file
parameter: column_mapping = ''
# A cis scan writes every gene into one table; empty = the file holds one trait.
parameter: trait_column = 'molecular_trait_id'
# Where the alleles live when there is no A1/A2 column: none | A2A1 | A1A2.
# Declaring this wrong is silent -- see the Input section.
parameter: variant_id_alleles = 'none'
# --- QC knobs (forwarded to summaryStatsQc) --------------------------
parameter: maf = 0.0
parameter: mac = 0.0
parameter: imiss = 1.0
parameter: z_mismatch_qc = 'none'     # none | slalom | dentist
parameter: pip_cutoff_to_skip = 0.0
parameter: skip_qc = False
# --- 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.5
parameter: pip_cutoff = 0.025
parameter: L = 10
parameter: L_greedy = 'none'   # 'none'/'off' = greedy off (default); a positive int enables greedy-L
parameter: ser_fallback = True
parameter: r_mismatch = 'none'        # none | eb | eb_mix
parameter: method_args = ''           # JSON {token: {kwarg: value}}
# --- TWAS-weight knobs (forwarded to twas_weights.R) -----------------
# Summary-statistics implementations only: mrash, lasso, scad, mcp, l0learn,
# mrmash, dpr_gibbs. enet and the bayes_* family are individual-level only.
parameter: twas_methods = 'lasso'
parameter: seed = 999
# --- cluster resources -----------------------------------------------
parameter: job_size = 1
parameter: walltime = '5h'
parameter: mem = '16G'
parameter: numThreads = 1
parameter: container = ''
parameter: entrypoint = ''

prefix = f'{study}.{context}.{trait}'

[qtl_rss_1, generate_qtl_sumstats]
# Read the cis-QTL nominal table, restrict it to this trait, attach the LD panel
# and run summaryStatsQc -> one QtlSumStats for the window.
output: f'{cwd:a}/sumstats/{prefix}.qtl_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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/qtl_sumstats_construct.R \
        --sumstats ${sumstats:a} \
        --study ${study} \
        --context ${context} \
        --trait ${trait} \
        --trait-column '${trait_column}' \
        --variant-id-alleles ${variant_id_alleles} \
        --ld-sketch ${ld_sketch:a} \
        --genome ${genome} \
        ${('--region ' + region) if region else ''} \
        ${('--n-sample ' + str(n_sample)) if n_sample >= 0 else ''} \
        ${('--column-mapping ' + column_mapping) if column_mapping else ''} \
        --maf ${maf} \
        --mac ${mac} \
        --imiss ${imiss} \
        --z-mismatch-qc ${z_mismatch_qc} \
        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \
        ${'--skip-qc' if skip_qc else ''} \
        --output ${_output}

[qtl_rss_2, qtl_rss_fine_mapping]
# SuSiE-RSS over the window: the same fine_mapping.R the individual-level and
# GWAS routes use, dispatching on the QtlSumStats class.
input: f'{cwd:a}/sumstats/{prefix}.qtl_sumstats.rds'
output: f'{cwd:a}/fine_mapping/{prefix}.qtl_rss_finemap.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --qtl-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} \
        --ser-fallback ${'TRUE' if ser_fallback else 'FALSE'} \
        --r-mismatch ${r_mismatch} \
        ${('--method-args ' + repr(method_args)) if method_args else ''} \
        --seed ${seed} \
        --output ${_output}

[qtl_rss_3, qtl_rss_twas_weights]
# RSS TWAS weights from the same collection. A QtlSumStats already spans one
# window, so no --gene-id / --region selection applies.
input: f'{cwd:a}/sumstats/{prefix}.qtl_sumstats.rds'
output: f'{cwd:a}/twas_weights/{prefix}.qtl_rss_twas_weights.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:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/twas_weights.R \
        --qtl-sumstats ${_input} \
        --methods ${twas_methods} \
        --seed ${seed} \
        --output ${_output}