add RSS TWAS/finemapping notebook

code/SoS/mnm_analysis/mnm_methods/qtl_rss_analysis.ipynb

@@ -0,0 +1,110 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "# RSS Fine-mapping and TWAS Weights with QTL Summary Statistics\n\nFine-maps a cis window and learns TWAS weights from cis-QTL **summary statistics** plus an\nLD reference panel, without individual-level genotypes."
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Overview\n\n`mnm_regression.ipynb` fits a cis window from individual-level data: it needs the genotype\nmatrix and the phenotype for every sample. That is not always possible -- summary statistics\nare shareable where genotypes are not, and a scan that has already run need not be repeated.\n\nThis module takes the other route. It reads a cis-QTL nominal association table (effect,\nstandard error and the variant id per SNP), pairs it with an LD reference panel, and hands\nthe resulting `QtlSumStats` to the same two pipelines the individual-level route uses:\n\n- **`qtl_rss_fine_mapping`** runs SuSiE-RSS (`susieR::susie_rss`) over the window. The\n  regression-with-summary-statistics likelihood replaces the individual-level one; credible\n  sets and PIPs mean what they always did.\n- **`qtl_rss_twas_weights`** learns predictive weights from the same object. Not every\n  method has a summary-statistics implementation: `mrash`, `lasso`, `scad`, `mcp`,\n  `l0learn`, `mrmash` and `dpr_gibbs` do; `enet` and the `bayes_*` family are\n  individual-level only and will be rejected here.\n\n**The LD panel is the load-bearing input.** RSS reconstructs the joint fit from marginal\nstatistics and LD, so the panel must cover the window's variants and must be on the same\nallele orientation. `summaryStatsQc()` harmonizes the two and reports what it corrected --\nread that line. If it says it sign- or strand-flipped everything, the alleles were declared\nwrong, not fixed (see `--variant-id-alleles` below).\n\n**When to run it.** After a cis scan (`TensorQTL.ipynb`), in place of\n`mnm_regression.ipynb`'s `susie_twas` when only summary statistics are available."
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Input\n\n- `--sumstats` **`output/cis/example.cis_qtl.pairs.tsv.gz`**\n(the cis-QTL nominal table from `TensorQTL.ipynb`: one row per (trait, variant) with\n`bhat`/`sebhat`, `pvalue`, `af` and `n`. A `z` column is used if present, otherwise the Wald\nz `bhat/sebhat` is derived.)\n\n- `--ld-sketch` **`tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed`**\n(the LD reference panel: a genotype path/prefix, or a per-chromosome LD meta file. Its\nvariants must cover the window.)\n\n- `--study` / `--context` / `--trait` -- the tuple this collection describes. `--trait` is\nthe gene id, and is also what the trait filter matches.\n\n- `--trait-column` (default `molecular_trait_id`) -- a cis scan writes *every* gene into one\ntable, so the rows for `--trait` are selected by this column. Set it empty for a file that\nalready holds one gene.\n\n- `--variant-id-alleles` (`none` | `A2A1` | `A1A2`) -- where the alleles come from when the\ntable has no `A1`/`A2` columns, as a cis scan typically does not: they sit inside the variant\nid. **The order cannot be inferred from the string** -- a `.pvar` writes `REF:ALT` (pecotmr's\ncanonical `A2:A1`) while a PLINK `.bim` writes `A1:A2` -- so declare which one your ids use.\nA real `A1`/`A2` column always wins. Declaring it wrong is silent: QC will \"correct\" the\napparent mismatch by flipping every variant against the panel.\n\n- `--genome` (default `GRCh38`), `--region`, `--n-sample`, `--column-mapping` -- optional.\n\nQC knobs are forwarded to `summaryStatsQc()`: `--maf`, `--mac`, `--imiss`,\n`--z-mismatch-qc`, `--pip-cutoff-to-skip`, and `--skip-qc` for diagnostics."
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Output\n\n- `{cwd}/sumstats/{study}.{context}.{trait}.qtl_sumstats.rds` -- the `QtlSumStats`: the\nwindow's variants with `SNP`/`A1`/`A2`/`Z`/`N` (plus `BETA`/`SE`/`P`/`AF` when supplied), the\nLD panel attached as the `ldSketch`, and a `qcInfo` audit of what QC did.\n- `{cwd}/fine_mapping/{...}.qtl_rss_finemap.rds` -- a `QtlFineMappingResult`: credible sets\nand PIPs from SuSiE-RSS, the same class the individual-level route produces.\n- `{cwd}/twas_weights/{...}.qtl_rss_twas_weights.rds` -- a `TwasWeights` collection.\n\nEach step also writes `.stdout` / `.stderr` beside its output. The QC line in the sumstats\nlog is worth reading every time:\n\n```\n[study/context/gene] QC summary: 200 in -> 200 out | corrected: sign-flip 0, strand-flip 0\n```"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Minimal Working Example\n\nRuns on the committed chr22 toy data: the TensorQTL nominal table for 16 genes, with the\nsame 49-sample genotypes used as the LD reference (in-sample LD -- fine for a smoke test,\noptimistic for real inference, where a separate reference panel belongs)."
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "sos run pipeline/qtl_rss_analysis.ipynb qtl_rss \\\n    --cwd output/qtl_rss \\\n    --sumstats tests/fixtures/tensorqtl/expected/cis_qtl.pairs.tsv.gz \\\n    --ld-sketch tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \\\n    --study test_study --context context1 --trait ENSG00000283047 \\\n    --variant-id-alleles A1A2 \\\n    --methods susie --twas-methods lasso -j1"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Command Interface"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "sos run pipeline/qtl_rss_analysis.ipynb -h"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "## Workflow implementation"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "[global]\nparameter: cwd = path('output')\nparameter: modular_script_dir = path('code/script')\n# --- the (study, context, trait) this run describes ------------------\nparameter: study = str\nparameter: context = str\nparameter: trait = str\n# --- inputs ----------------------------------------------------------\nparameter: sumstats = path\nparameter: ld_sketch = path\nparameter: genome = 'GRCh38'\nparameter: region = ''\nparameter: n_sample = -1.0            # study-level total N; <0 = take it from the file\nparameter: column_mapping = ''\n# A cis scan writes every gene into one table; empty = the file holds one trait.\nparameter: trait_column = 'molecular_trait_id'\n# Where the alleles live when there is no A1/A2 column: none | A2A1 | A1A2.\n# Declaring this wrong is silent -- see the Input section.\nparameter: variant_id_alleles = 'none'\n# --- QC knobs (forwarded to summaryStatsQc) --------------------------\nparameter: maf = 0.0\nparameter: mac = 0.0\nparameter: imiss = 1.0\nparameter: z_mismatch_qc = 'none'     # none | slalom | dentist\nparameter: pip_cutoff_to_skip = 0.0\nparameter: skip_qc = False\n# --- fine-mapping knobs (forwarded to fine_mapping.R) ----------------\nparameter: methods = 'susie'\nparameter: coverage = 0.95\nparameter: secondary_coverage = '0.7,0.5'\nparameter: min_abs_corr = 0.5\nparameter: pip_cutoff = 0.025\nparameter: L = 10\nparameter: L_greedy = 'none'   # 'none'/'off' = greedy off (default); a positive int enables greedy-L\nparameter: ser_fallback = True\nparameter: r_mismatch = 'none'        # none | eb | eb_mix\nparameter: method_args = ''           # JSON {token: {kwarg: value}}\n# --- TWAS-weight knobs (forwarded to twas_weights.R) -----------------\n# Summary-statistics implementations only: mrash, lasso, scad, mcp, l0learn,\n# mrmash, dpr_gibbs. enet and the bayes_* family are individual-level only.\nparameter: twas_methods = 'lasso'\nparameter: seed = 999\n# --- cluster resources -----------------------------------------------\nparameter: job_size = 1\nparameter: walltime = '5h'\nparameter: mem = '16G'\nparameter: numThreads = 1\nparameter: container = ''\nparameter: entrypoint = ''\n\nprefix = f'{study}.{context}.{trait}'"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "[qtl_rss_1, generate_qtl_sumstats]\n# Read the cis-QTL nominal table, restrict it to this trait, attach the LD panel\n# and run summaryStatsQc -> one QtlSumStats for the window.\noutput: f'{cwd:a}/sumstats/{prefix}.qtl_sumstats.rds'\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\nbash: expand = '${ }', stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint\n    Rscript ${modular_script_dir}/pecotmr_integration/qtl_sumstats_construct.R \\\n        --sumstats ${sumstats:a} \\\n        --study ${study} \\\n        --context ${context} \\\n        --trait ${trait} \\\n        --trait-column '${trait_column}' \\\n        --variant-id-alleles ${variant_id_alleles} \\\n        --ld-sketch ${ld_sketch:a} \\\n        --genome ${genome} \\\n        ${('--region ' + region) if region else ''} \\\n        ${('--n-sample ' + str(n_sample)) if n_sample >= 0 else ''} \\\n        ${('--column-mapping ' + column_mapping) if column_mapping else ''} \\\n        --maf ${maf} \\\n        --mac ${mac} \\\n        --imiss ${imiss} \\\n        --z-mismatch-qc ${z_mismatch_qc} \\\n        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \\\n        ${'--skip-qc' if skip_qc else ''} \\\n        --output ${_output}"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "[qtl_rss_2, qtl_rss_fine_mapping]\n# SuSiE-RSS over the window: the same fine_mapping.R the individual-level and\n# GWAS routes use, dispatching on the QtlSumStats class.\ninput: f'{cwd:a}/sumstats/{prefix}.qtl_sumstats.rds'\noutput: f'{cwd:a}/fine_mapping/{prefix}.qtl_rss_finemap.rds'\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\nbash: expand = '${ }', stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint\n    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n        --qtl-sumstats ${_input} \\\n        --methods ${methods} \\\n        --coverage ${coverage} \\\n        --secondary-coverage ${secondary_coverage} \\\n        --min-abs-corr ${min_abs_corr} \\\n        --pip-cutoff ${pip_cutoff} \\\n        --L ${L} \\\n        --L-greedy ${L_greedy} \\\n        --ser-fallback ${'TRUE' if ser_fallback else 'FALSE'} \\\n        --r-mismatch ${r_mismatch} \\\n        ${('--method-args ' + repr(method_args)) if method_args else ''} \\\n        --seed ${seed} \\\n        --output ${_output}"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "[qtl_rss_3, qtl_rss_twas_weights]\n# RSS TWAS weights from the same collection. A QtlSumStats already spans one\n# window, so no --gene-id / --region selection applies.\ninput: f'{cwd:a}/sumstats/{prefix}.qtl_sumstats.rds'\noutput: f'{cwd:a}/twas_weights/{prefix}.qtl_rss_twas_weights.rds'\ntask: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\nbash: expand = '${ }', stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint\n    Rscript ${modular_script_dir}/pecotmr_integration/twas_weights.R \\\n        --qtl-sumstats ${_input} \\\n        --methods ${twas_methods} \\\n        --seed ${seed} \\\n        --output ${_output}"
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "SoS",
+   "language": "sos",
+   "name": "sos"
+  },
+  "language_info": {
+   "codemirror_mode": "sos",
+   "file_extension": ".sos",
+   "mimetype": "text/x-sos",
+   "name": "sos",
+   "nbconvert_exporter": "sos_notebook.converter.SoS_Exporter",
+   "pygments_lexer": "sos"
+  },
+  "sos": {
+   "kernels": [
+    [
+     "SoS",
+     "sos",
+     "sos",
+     "",
+     ""
+    ]
+   ],
+   "version": ""
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

code/script/pecotmr_integration/fine_mapping.R

@@ -16,6 +16,9 @@
 #
 # GWAS — one call per per-block GwasSumStats RDS (each carrying its own
 # z-scores + LD sketch; no gene/region concept):
+#   --qtl-sumstats <RDS>           pecotmr::QtlSumStats (QTL RSS mode, from
+#                                   qtl_sumstats_construct.R): SuSiE-RSS over a
+#                                   cis window using the collection's ldSketch
 #   --gwas-sumstats <RDS>           pecotmr::GwasSumStats (per LD block,
 #                                   typically from gwas_sumstats_construct.R)
 #
@@ -64,6 +67,9 @@ parser <- add_argument(parser, "--qtl-dataset",
 parser <- add_argument(parser, "--gwas-sumstats",
                        help = "Path to a GwasSumStats RDS (GWAS mode)",
                        type = "character", default = "")
+parser <- add_argument(parser, "--qtl-sumstats",
+                       help = "Path to a QtlSumStats RDS (QTL RSS mode; from qtl_sumstats_construct.R)",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--gene-id",
                        help = "Trait identifier (QTL gene mode); mutually exclusive with --region",
                        type = "character", default = "")
@@ -275,10 +281,11 @@ if (!is.null(seed_val)) set.seed(seed_val)
 
 has_qtl  <- nzchar(argv$qtl_dataset)
 has_gwas <- nzchar(argv$gwas_sumstats)
-if (has_qtl && has_gwas)
-  stop("--qtl-dataset and --gwas-sumstats are mutually exclusive; pass exactly one.")
-if (!has_qtl && !has_gwas)
-  stop("Specify either --qtl-dataset (QTL mode) or --gwas-sumstats (GWAS mode).")
+has_qss  <- nzchar(argv$qtl_sumstats)
+if (sum(has_qtl, has_gwas, has_qss) > 1L)
+  stop("--qtl-dataset, --gwas-sumstats and --qtl-sumstats are mutually exclusive; pass exactly one.")
+if (!has_qtl && !has_gwas && !has_qss)
+  stop("Specify --qtl-dataset (QTL individual-level), --qtl-sumstats (QTL RSS) or --gwas-sumstats (GWAS RSS).")
 
 # Optional context restriction (QTL mode): NULL = all contexts in the dataset.
 contexts_arg <- if (nzchar(argv$contexts) && argv$contexts != ".")
@@ -336,9 +343,12 @@ if (!is.null(median_abs_corr)) cs_args$medianAbsCorr <- median_abs_corr
 # fineMappingResult). Added only when supplied, for pecotmr-version tolerance.
 if (!is.null(fmr_obj)) cs_args$fineMappingResult <- fmr_obj
 
-if (has_gwas) {
-  # ----- GWAS mode -------------------------------------------------------
-  gss <- readRDS(argv$gwas_sumstats)
+if (has_gwas || has_qss) {
+  # ----- RSS mode (GwasSumStats or QtlSumStats) --------------------------
+  # The SuSiE-RSS knobs below are pipeline arguments of the summary-statistics
+  # path, which pecotmr documents as "QtlSumStats / GwasSumStats only", so the
+  # two RSS inputs share this branch; fineMappingPipeline dispatches on class.
+  gss <- readRDS(if (has_gwas) argv$gwas_sumstats else argv$qtl_sumstats)
   ser_fallback <- as.logical(argv$ser_fallback)
   if (is.na(ser_fallback))
     stop("--ser-fallback must be TRUE or FALSE (got: ", argv$ser_fallback, ")")
@@ -368,8 +378,10 @@ if (has_gwas) {
            '\'{"check_prior":true,"mismatch_estimator":"map"}\'.')
     gwas_args$rssControl <- rc
   }
+  if (!is.null(seed_val)) gwas_args$seed <- seed_val
   res <- do.call(fineMappingPipeline, c(list(gss), gwas_args))
-  label <- paste0("GwasSumStats '", basename(argv$gwas_sumstats), "'")
+  label <- if (has_gwas) paste0("GwasSumStats '", basename(argv$gwas_sumstats), "'")
+           else paste0("QtlSumStats '", basename(argv$qtl_sumstats), "'")
 } else {
   # ----- QTL mode --------------------------------------------------------
   has_gene   <- nzchar(argv$gene_id)