diff --git a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb index 88c18435..6a3b5856 100644 --- a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb @@ -8,18 +8,18 @@ "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", @@ -27,32 +27,35 @@ "\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" ] }, @@ -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:" @@ -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", @@ -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}" @@ -299,4 +308,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file diff --git a/code/script/pecotmr_integration/gwas_sumstats_construct.R b/code/script/pecotmr_integration/gwas_sumstats_construct.R index 077057d1..4da48f83 100644 --- a/code/script/pecotmr_integration/gwas_sumstats_construct.R +++ b/code/script/pecotmr_integration/gwas_sumstats_construct.R @@ -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({ @@ -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") @@ -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")) @@ -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) }