diff --git a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb index 6fc6db24..b1f52356 100644 --- a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb @@ -133,59 +133,7 @@ "kernel": "SoS" }, "outputs": [], - "source": [ - "[global]\n", - "parameter: cwd = path('output')\n", - "parameter: modular_script_dir = path('code/script')\n", - "# --- per-study GWAS inputs (use either or both) ---------------------\n", - "parameter: gwas_meta = path('.')\n", - "parameter: gwas_tsv_list = [] # list of STUDY=PATH items\n", - "# --- region inputs (use either or both) -----------------------------\n", - "parameter: region_list = path('.')\n", - "parameter: regions = [] # list of chr:start-end strings\n", - "# --- LD reference + genome ------------------------------------------\n", - "parameter: ld_meta = path\n", - "parameter: genome = 'GRCh38'\n", - "# --- QC knobs (forwarded to gwas_sumstats_construct.R) --------------\n", - "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", - "parameter: allele_flip_kriging = False # kriging allele-flip QC: sign-flip switched z (logLR>2 & |z|>2); off by default\n", - "parameter: effective_n = True # case/control: use effective sample size 4/(1/n_case+1/n_control) as N (summaryStatsQc effectiveN); False -> raw N\n", - "# --- fine-mapping knobs (forwarded to fine_mapping.R) ---------------\n", - "parameter: methods = 'susie'\n", - "parameter: coverage = 0.95\n", - "parameter: secondary_coverage = '0.7,0.5'\n", - "parameter: min_abs_corr = 0.8\n", - "parameter: median_abs_corr = '' # empty -> off (OR-logic purity; needs pecotmr step-1)\n", - "parameter: pip_cutoff = 0.025\n", - "parameter: L = 20\n", - "parameter: L_greedy = 5\n", - "parameter: method_args = '' # nested per-method JSON object\n", - "# GWAS SuSiE-RSS fine-mapping: SER fallback is ON by default (reproduces the single-panel notebook):\n", - "# fit finite-sample R + EB LD-mismatch SuSiE-RSS and fall back to the single-effect (SER) result for\n", - "# regions susieR flags as unreliable. ser_fallback and r_mismatch are orthogonal: ser_fallback = False\n", - "# only disables the SER fallback (EB LD-mismatch stays on); pass r_mismatch = 'none' too for the\n", - "# pre-change plain multi-effect fit.\n", - "parameter: ser_fallback = True # fall back to single-effect (SER) when susieR flags R unreliable\n", - "parameter: r_finite = '' # finite-sample R size; empty -> auto (LD-panel sample size)\n", - "parameter: r_mismatch = 'eb' # LD-mismatch correction: 'eb' (default) or 'none'\n", - "parameter: r_mismatch_method = '' # rMismatch estimator 'mle'|'map'; empty -> pecotmr default\n", - "parameter: check_prior = '' # susie_rss check_prior 'TRUE'|'FALSE'; empty -> pecotmr default\n", - "# --- plot opt-out ---------------------------------------------------\n", - "parameter: no_plot = False\n", - "# --- output prefix for the manifest file ----------------------------\n", - "parameter: manifest_name = 'gwas_rss'\n", - "# --- infrastructure -------------------------------------------------\n", - "parameter: container = ''\n", - "parameter: job_size = 1\n", - "parameter: walltime = '2h'\n", - "parameter: mem = '16G'\n", - "parameter: numThreads = 1" - ] + "source": "[global]\nparameter: cwd = path('output')\nparameter: modular_script_dir = path('code/script')\n# --- per-study GWAS inputs (use either or both) ---------------------\nparameter: gwas_meta = path('.')\nparameter: gwas_tsv_list = [] # list of STUDY=PATH items\n# --- region inputs (use either or both) -----------------------------\nparameter: region_list = path('.')\nparameter: regions = [] # list of chr:start-end strings\n# --- LD reference + genome ------------------------------------------\nparameter: ld_meta = path\nparameter: genome = 'GRCh38'\n# --- QC knobs (forwarded to gwas_sumstats_construct.R) --------------\nparameter: qc_method = 'none' # none | slalom | dentist\nparameter: impute = False\nparameter: qc_args = '' # JSON object spliced into summaryStatsQc()\nparameter: maf = 0.0025 # MAF cutoff (summaryStatsQc mafCutoff); 0 disables\nparameter: skip_regions = [] # chr:start-end window(s) whose variants are dropped (skipRegion)\nparameter: skip_analysis_pip_cutoff = 0.025 # skip a region below this max PIP (pipCutoffToSkip); 0 disables\nparameter: allele_flip_kriging = False # kriging allele-flip QC: sign-flip switched z (logLR>2 & |z|>2); off by default\nparameter: effective_n = True # case/control: use effective sample size 4/(1/n_case+1/n_control) as N (summaryStatsQc effectiveN); False -> raw N\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.8\nparameter: median_abs_corr = '' # empty -> off (OR-logic purity; needs pecotmr step-1)\nparameter: pip_cutoff = 0.025\nparameter: L = 20\nparameter: L_greedy = 5\nparameter: method_args = '' # nested per-method JSON object\n# GWAS SuSiE-RSS fine-mapping: SER fallback is ON by default (reproduces the single-panel notebook):\n# fit finite-sample R + EB LD-mismatch SuSiE-RSS and fall back to the single-effect (SER) result for\n# regions susieR flags as unreliable. ser_fallback and r_mismatch are orthogonal: ser_fallback = False\n# only disables the SER fallback (EB LD-mismatch stays on); pass r_mismatch = 'none' too for the\n# pre-change plain multi-effect fit. Requires susieR >= 0.16.6 (which moved R_mismatch_method /\n# check_prior into susie_rss_control(); pass those and other control settings via rss_control).\nparameter: ser_fallback = True # fall back to single-effect (SER) when susieR flags R unreliable\nparameter: r_finite = '' # finite-sample R size; empty -> auto (LD-panel sample size)\nparameter: r_mismatch = 'eb_mix' # LD-mismatch correction: 'none' | 'eb' | 'eb_mix' (residual-mixture EB, default)\nparameter: rss_control = '' # JSON object of susie_rss_control() settings, e.g. '{\"check_prior\":true,\"mismatch_estimator\":\"map\"}'; empty -> defaults\n# --- plot opt-out ---------------------------------------------------\nparameter: no_plot = False\n# --- output prefix for the manifest file ----------------------------\nparameter: manifest_name = 'gwas_rss'\n# --- infrastructure -------------------------------------------------\nparameter: container = ''\nparameter: job_size = 1\nparameter: walltime = '2h'\nparameter: mem = '16G'\nparameter: numThreads = 1" }, { "cell_type": "code", @@ -227,30 +175,7 @@ "kernel": "SoS" }, "outputs": [], - "source": [ - "[gwas_fine_mapping]\n", - "# Per-RDS fan-out: one fine-mapping task per (study, region) GwasSumStats.\n", - "output: f\"{cwd}/fine_mapping/{_input:bnn}.gwas_finemap.rds\", group_by = 1\n", - "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n", - "bash: expand = '${ }', stderr = f\"{_output}.stderr\", stdout = f\"{_output}.stdout\", container = container\n", - " Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n", - " --gwas-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", - " ${('--median-abs-corr ' + str(median_abs_corr)) if str(median_abs_corr) != '' else ''} \\\n", - " ${('--method-args ' + repr(method_args)) if method_args else ''} \\\n", - " --ser-fallback ${'TRUE' if ser_fallback else 'FALSE'} \\\n", - " --r-mismatch ${r_mismatch} \\\n", - " ${('--r-finite ' + str(r_finite)) if str(r_finite) != '' else ''} \\\n", - " ${('--r-mismatch-method ' + r_mismatch_method) if r_mismatch_method else ''} \\\n", - " ${('--check-prior ' + check_prior) if check_prior else ''} \\\n", - " --output ${_output}" - ] + "source": "[gwas_fine_mapping]\n# Per-RDS fan-out: one fine-mapping task per (study, region) GwasSumStats.\noutput: f\"{cwd}/fine_mapping/{_input:bnn}.gwas_finemap.rds\", group_by = 1\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}.stderr\", stdout = f\"{_output}.stdout\", container = container\n Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n --gwas-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 ${('--median-abs-corr ' + str(median_abs_corr)) if str(median_abs_corr) != '' else ''} \\\n ${('--method-args ' + repr(method_args)) if method_args else ''} \\\n --ser-fallback ${'TRUE' if ser_fallback else 'FALSE'} \\\n --r-mismatch ${r_mismatch} \\\n ${('--r-finite ' + str(r_finite)) if str(r_finite) != '' else ''} \\\n ${('--rss-control ' + repr(rss_control)) if rss_control else ''} \\\n --output ${_output}" }, { "cell_type": "code", diff --git a/code/script/pecotmr_integration/fine_mapping.R b/code/script/pecotmr_integration/fine_mapping.R index 3dbfd7c5..3299c466 100644 --- a/code/script/pecotmr_integration/fine_mapping.R +++ b/code/script/pecotmr_integration/fine_mapping.R @@ -119,8 +119,8 @@ parser <- add_argument(parser, "--pip-cutoff-to-skip", # fallback is ON by default (reproduces the single-panel notebook): fit a # finite-sample R + EB LD-mismatch SuSiE-RSS and drop to the single-effect (SER) # result for regions susieR flags as unreliable. These ride on gwas_args (GWAS -# mode only); rFinite/rMismatchMethod/checkPrior forward only when set, so the -# call also runs against a pecotmr that predates them. +# mode only); rFinite forwards only when set. susie_rss_control() settings +# (susieR >= 0.16.6) are passed as a JSON object via --rss-control. parser <- add_argument(parser, "--ser-fallback", help = "GWAS mode: TRUE/FALSE. Fall back to the single-effect (SER) result when susieR flags R unreliable (fineMappingPipeline serFallback). Default TRUE.", type = "character", default = "TRUE") @@ -128,13 +128,10 @@ parser <- add_argument(parser, "--r-finite", help = "GWAS mode: finite-sample R correction size (fineMappingPipeline rFinite). Empty = auto (LD-panel sample size).", type = "character", default = "") parser <- add_argument(parser, "--r-mismatch", - help = "GWAS mode: LD-mismatch correction (fineMappingPipeline rMismatch): 'eb' (default) or 'none'.", - type = "character", default = "eb") -parser <- add_argument(parser, "--r-mismatch-method", - help = "GWAS mode: rMismatch estimator (fineMappingPipeline rMismatchMethod): 'mle' or 'map'. Empty = pecotmr default.", - type = "character", default = "") -parser <- add_argument(parser, "--check-prior", - help = "GWAS mode: TRUE/FALSE for susie_rss check_prior (fineMappingPipeline checkPrior). Empty = pecotmr default.", + help = "GWAS mode: LD-mismatch correction (fineMappingPipeline rMismatch): 'none', 'eb', or 'eb_mix' (residual-mixture EB, susieR >= 0.16.6). Default 'eb_mix'.", + type = "character", default = "eb_mix") +parser <- add_argument(parser, "--rss-control", + help = "GWAS mode: JSON object of susieR::susie_rss_control() settings (e.g. '{\"check_prior\":true,\"mismatch_estimator\":\"map\"}'), forwarded as fineMappingPipeline rssControl. Empty = susie_rss_control() defaults.", type = "character", default = "") # --- Multivariate / joint-fit knobs (QTL mode; mvsusie / fsusie). Each is # opt-in and omitted from the pipeline call when left at its default, so this @@ -317,18 +314,23 @@ if (has_gwas) { ser_fallback <- as.logical(argv$ser_fallback) if (is.na(ser_fallback)) stop("--ser-fallback must be TRUE or FALSE (got: ", argv$ser_fallback, ")") - # GWAS-only SuSiE-RSS knobs on top of the shared cs_args. rFinite / - # rMismatchMethod / checkPrior forward only when set (empty -> pecotmr default). + # GWAS-only SuSiE-RSS knobs on top of the shared cs_args. rFinite forwards + # only when set; --rss-control (JSON) becomes the rssControl named list of + # susie_rss_control() settings. gwas_args <- c(cs_args, list(serFallback = ser_fallback, rMismatch = argv$r_mismatch)) if (nzchar(argv$r_finite)) gwas_args$rFinite <- as.numeric(argv$r_finite) - if (nzchar(argv[["r_mismatch_method"]])) - gwas_args$rMismatchMethod <- argv[["r_mismatch_method"]] - if (nzchar(argv$check_prior)) { - cp <- as.logical(argv$check_prior) - if (is.na(cp)) stop("--check-prior must be TRUE or FALSE (got: ", argv$check_prior, ")") - gwas_args$checkPrior <- cp + if (nzchar(argv$rss_control) && argv$rss_control != "." && + argv$rss_control != "{}") { + rc <- tryCatch(jsonlite::fromJSON(argv$rss_control, simplifyVector = TRUE), + error = function(e) stop( + "--rss-control must be a JSON object (got: ", argv$rss_control, + "). Error: ", conditionMessage(e))) + if (!is.list(rc) || is.null(names(rc)) || any(!nzchar(names(rc)))) + stop("--rss-control must be a JSON object with named fields, e.g. ", + '\'{"check_prior":true,"mismatch_estimator":"map"}\'.') + gwas_args$rssControl <- rc } res <- do.call(fineMappingPipeline, c(list(gss), gwas_args)) label <- paste0("GwasSumStats '", basename(argv$gwas_sumstats), "'")