Skip to content

Expose maf, skip_regions, skip_analysis_pip_cutoff in rss_analysis #1386

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
gaow merged 1 commit into from
Jul 8, 2026
+55 −23
Merged

Expose maf, skip_regions, skip_analysis_pip_cutoff in rss_analysis #1386

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 27 additions & 18 deletions code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,51 +8,54 @@
"source": [
"# RSS Fine-mapping with GWAS Summary Statistics\n",
"\n",
"End-to-end SuSiE-RSS fine-mapping over GWAS summary statistics, per study × per LD block, using\n",
"End-to-end SuSiE-RSS fine-mapping over GWAS summary statistics, per study \u00d7 per LD block, using\n",
"the current pecotmr S4 API (`GwasSumStats` / `summaryStatsQc()` / `fineMappingPipeline()`). Four SoS\n",
"steps:\n",
"\n",
"1. **`[generate_manifest]`** — `gwas_rss_manifest.R` resolves `--gwas-meta` + `--gwas-tsv-list` against\n",
" `--region-list` + `--regions` into one manifest TSV (one row per study × region).\n",
"2. **`[generate_gwas_sumstats]`** — `gwas_sumstats_construct.R` builds one `GwasSumStats` RDS per\n",
" study × LD block, running `summaryStatsQc()` (optional SLALOM/DENTIST QC + RAISS imputation;\n",
"1. **`[generate_manifest]`** \u2014 `gwas_rss_manifest.R` resolves `--gwas-meta` + `--gwas-tsv-list` against\n",
" `--region-list` + `--regions` into one manifest TSV (one row per study \u00d7 region).\n",
"2. **`[generate_gwas_sumstats]`** \u2014 `gwas_sumstats_construct.R` builds one `GwasSumStats` RDS per\n",
" study \u00d7 LD block, running `summaryStatsQc()` (optional SLALOM/DENTIST QC + RAISS imputation;\n",
" harmonization re-keys variants to the LD-panel id convention).\n",
"3. **`[gwas_fine_mapping]`** — `fine_mapping.R` dispatches to\n",
"3. **`[gwas_fine_mapping]`** \u2014 `fine_mapping.R` dispatches to\n",
" `pecotmr::fineMappingPipeline(gwasSumStats, methods = ...)` (SuSiE-RSS).\n",
"4. **`[gwas_rss_plot]`** (opt-out via `--no-plot`) — `gwas_rss_plot.R` renders a PIP plot.\n",
"4. **`[gwas_rss_plot]`** (opt-out via `--no-plot`) \u2014 `gwas_rss_plot.R` renders a PIP plot.\n",
"\n",
"This replaces the legacy two-step workflow (`get_analysis_regions` + `univariate_rss`, which called the\n",
"removed `rss_analysis_pipeline()`). The LD reference must be **genotype-backed** (PLINK2/PLINK1/GDS/VCF);\n",
"a precomputed `.cor.xz` LD matrix is not supported for the `GwasSumStats` path.\n",
"\n",
"## Inputs\n",
"\n",
"- `--gwas-meta` — TSV with columns `study_id`, `path` (relative to the meta file or absolute), and\n",
"- `--gwas-meta` \u2014 TSV with columns `study_id`, `path` (relative to the meta file or absolute), and\n",
" optional `column_mapping` (per-study YAML). One row per study.\n",
"- `--gwas-tsv-list S1=path1 ...` — explicit STUDY=PATH pairs (alternative to `--gwas-meta`; no\n",
"- `--gwas-tsv-list S1=path1 ...` \u2014 explicit STUDY=PATH pairs (alternative to `--gwas-meta`; no\n",
" per-study column mapping).\n",
"- `--region-list` — BED-like TSV of LD blocks; and/or `--regions chr:start-end ...`.\n",
"- `--ld-meta` — LD-reference meta TSV (`#chr`, `start`, `end`, `path` -> genotype prefix). One LD block\n",
"- `--region-list` \u2014 BED-like TSV of LD blocks; and/or `--regions chr:start-end ...`.\n",
"- `--ld-meta` \u2014 LD-reference meta TSV (`#chr`, `start`, `end`, `path` -> genotype prefix). One LD block\n",
" per row (or one row with `start=end=0` for a whole chromosome).\n",
"\n",
"## QC knobs (forwarded to `gwas_sumstats_construct.R` -> `summaryStatsQc()`)\n",
"\n",
"- `--qc-method` — `none` (default), `slalom`, or `dentist`.\n",
"- `--impute` — enable RAISS sumstat imputation.\n",
"- `--qc-args '{\"key\":value,...}'` — JSON passthrough for any other `summaryStatsQc()` kwarg (e.g.\n",
"- `--qc-method` \u2014 `none` (default), `slalom`, or `dentist`.\n",
"- `--impute` \u2014 enable RAISS sumstat imputation.\n",
"- `--qc-args '{\"key\":value,...}'` \u2014 JSON passthrough for any other `summaryStatsQc()` kwarg (e.g.\n",
" `{\"mafCutoff\":0.01,\"nCutoff\":0,\"alleleFlipKriging\":true}`).\n",
"- `--maf` \u2014 MAF cutoff (`mafCutoff`). Default **0.0025**; 0 disables.\n",
"- `--skip-regions chr:start-end \u2026` \u2014 drop variants overlapping these window(s) (`skipRegion`).\n",
"- `--skip-analysis-pip-cutoff` \u2014 skip a region whose single-effect max PIP is below this (`pipCutoffToSkip`). Default **0.025**; 0 disables.\n",
"\n",
"## Fine-mapping knobs (forwarded to `fine_mapping.R`)\n",
"\n",
"- `--methods` (default `susie`), `--coverage` (0.95), `--secondary-coverage` (`0.7,0.5`).\n",
"- `--min-abs-corr` (0.8), `--median-abs-corr` (off; OR-logic purity), `--pip-cutoff` (0.025).\n",
"- `--L` (20) / `--L-greedy` (5) — SuSiE fit defaults; `--method-args '{\"susie\":{...}}'` overrides per key.\n",
"- `--L` (20) / `--L-greedy` (5) \u2014 SuSiE fit defaults; `--method-args '{\"susie\":{...}}'` overrides per key.\n",
"\n",
"## Outputs\n",
"\n",
"- `{cwd}/{manifest_name}.manifest.tsv`\n",
"- `{cwd}/sumstats/{study}.{chrN_S_E}.gwas_sumstats.rds`\n",
"- `{cwd}/fine_mapping/{study}.{chrN_S_E}.gwas_finemap.rds` — one `GwasFineMappingResult` per (study, block).\n",
"- `{cwd}/fine_mapping/{study}.{chrN_S_E}.gwas_finemap.rds` \u2014 one `GwasFineMappingResult` per (study, block).\n",
"- `{cwd}/plots/{study}.{chrN_S_E}.pip_plot.png` (unless `--no-plot`).\n"
]
},
Expand Down Expand Up @@ -90,7 +93,7 @@
"kernel": "SoS"
},
"source": [
"This fine-maps one study × one LD block and writes the manifest, the QC'd `GwasSumStats`, the\n",
"This fine-maps one study \u00d7 one LD block and writes the manifest, the QC'd `GwasSumStats`, the\n",
"`GwasFineMappingResult`, and a PIP plot under `output/rss_analysis/`.\n",
"\n",
"#### With SLALOM QC + RAISS imputation + per-method tuning:"
Expand Down Expand Up @@ -147,6 +150,9 @@
"parameter: qc_method = 'none' # none | slalom | dentist\n",
"parameter: impute = False\n",
"parameter: qc_args = '' # JSON object spliced into summaryStatsQc()\n",
"parameter: maf = 0.0025 # MAF cutoff (summaryStatsQc mafCutoff); 0 disables\n",
"parameter: skip_regions = [] # chr:start-end window(s) whose variants are dropped (skipRegion)\n",
"parameter: skip_analysis_pip_cutoff = 0.025 # skip a region below this max PIP (pipCutoffToSkip); 0 disables\n",
"# --- fine-mapping knobs (forwarded to fine_mapping.R) ---------------\n",
"parameter: methods = 'susie'\n",
"parameter: coverage = 0.95\n",
Expand Down Expand Up @@ -219,6 +225,9 @@
" --genome ${genome} \\\n",
" ${('--column-mapping ' + _jobs['column_mapping']) if _jobs['column_mapping'] else ''} \\\n",
" --qc-method ${qc_method} \\\n",
" --maf ${maf} \\\n",
" --pip-cutoff-to-skip ${skip_analysis_pip_cutoff} \\\n",
" ${('--skip-region ' + ','.join(str(x) for x in skip_regions)) if skip_regions else ''} \\\n",
" ${'--impute' if impute else ''} \\\n",
" ${('--qc-args ' + repr(qc_args)) if qc_args else ''} \\\n",
" --output ${_output}"
Expand Down Expand Up @@ -299,4 +308,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
}
}
33 changes: 28 additions & 5 deletions code/script/pecotmr_integration/gwas_sumstats_construct.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,13 @@
# (default), "slalom", "dentist".
# --impute Flag: run RAISS sumstat imputation in
# summaryStatsQc(impute = TRUE) against the LD sketch.
# --maf MAF cutoff (summaryStatsQc mafCutoff). Default 0.0025.
# --skip-region Comma-separated chr:start-end window(s) whose variants
# are dropped (summaryStatsQc skipRegion).
# --qc-args Optional JSON object of extra named kwargs spliced
# into summaryStatsQc().
# into summaryStatsQc(). May not set a key already
# controlled by a dedicated flag (mafCutoff / skipRegion /
# pipCutoffToSkip / zMismatchQc / impute).
# --output Output RDS path

suppressPackageStartupMessages({
Expand Down Expand Up @@ -87,6 +92,12 @@ parser <- add_argument(parser, "--n-control",
parser <- add_argument(parser, "--pip-cutoff-to-skip",
help = "Skip a study whose single-trait max PIP is below this cutoff (summaryStatsQc pipCutoffToSkip); 0 disables, <0 uses 3/n_variants.",
type = "numeric", default = 0)
parser <- add_argument(parser, "--maf",
help = "Minor-allele-frequency cutoff (summaryStatsQc mafCutoff); drop variants below it. 0 disables.",
type = "numeric", default = 0.0025)
parser <- add_argument(parser, "--skip-region",
help = "Comma-separated chr:start-end window(s) whose overlapping variants are dropped (summaryStatsQc skipRegion). Empty = none.",
type = "character", default = "")
parser <- add_argument(parser, "--qc-method",
help = "LD-mismatch QC: 'none' (default), 'slalom', or 'dentist'",
type = "character", default = "none")
Expand Down Expand Up @@ -163,12 +174,21 @@ qc_extra <- if (nzchar(argv$qc_args) && argv$qc_args != "." &&
} else {
list()
}
# Reject collisions between explicit flags and --qc-args.
clash <- intersect(names(qc_extra), c("zMismatchQc", "impute"))
# Reject collisions between explicit flags and --qc-args. Dedicated flags win;
# passing the same key via --qc-args is an error so behavior is unambiguous.
clash <- intersect(names(qc_extra),
c("zMismatchQc", "impute", "mafCutoff", "skipRegion",
"pipCutoffToSkip"))
if (length(clash) > 0L)
stop("--qc-args sets ", paste(clash, collapse = ", "),
" which is also controlled by a dedicated flag (--qc-method / ",
"--impute). Pass it via the dedicated flag.")
"--impute / --maf / --skip-region / --pip-cutoff-to-skip). Pass it via ",
"the dedicated flag.")

# --skip-region: comma-separated chr:start-end windows -> character vector
# (NULL when unset so summaryStatsQc's skipRegion default applies).
skip_region_vec <- splitCsv(argv$skip_region)
if (length(skip_region_vec) == 0L) skip_region_vec <- NULL

qc_method <- match.arg(argv$qc_method, c("none", "slalom", "dentist"))

Expand Down Expand Up @@ -334,7 +354,10 @@ gss_out <- if (argv$skip_qc) {
qc_call_args <- c(list(gss,
zMismatchQc = qc_method,
impute = isTRUE(argv$impute),
pipCutoffToSkip = argv$pip_cutoff_to_skip),
pipCutoffToSkip = argv$pip_cutoff_to_skip,
mafCutoff = argv$maf),
if (!is.null(skip_region_vec))
list(skipRegion = skip_region_vec) else list(),
qc_extra)
do.call(summaryStatsQc, qc_call_args)
}
Expand Down
Loading