From f000a7d760b4941af94691842f8da85530feaeaa Mon Sep 17 00:00:00 2001 From: Yining97 Date: Sun, 19 Jul 2026 23:55:37 -0400 Subject: [PATCH] rss_analysis: effective sample size for case/control GWAS construct.R attaches per-variant n_case/n_control (N_CASE/N_CONTROL) and no longer requires a precomputed N when they are present; --effective-n (default on) forwards to summaryStatsQc(effectiveN=). Study-level fallback: gwas_rss_manifest.R carries n_case/n_control from --gwas-meta into the manifest and the construct step forwards them as --n-case/--n-control. effective_n exposed in rss_analysis.ipynb. --- .../mnm_methods/rss_analysis.ipynb | 6 +++++- .../pecotmr_integration/gwas_rss_manifest.R | 19 ++++++++++++++++--- .../gwas_sumstats_construct.R | 14 +++++++++++--- 3 files changed, 32 insertions(+), 7 deletions(-) diff --git a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb index 549f65ab5..43515e68a 100644 --- a/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/rss_analysis.ipynb @@ -154,6 +154,7 @@ "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", @@ -221,7 +222,7 @@ "[generate_gwas_sumstats]\n", "# Fan out over the manifest's rows: one (study, region) GwasSumStats per\n", "# row. Manifest columns: study_id, gwas_tsv, column_mapping, chr, start,\n", - "# end, region_id, gwas_tsv_basename.\n", + "# n_case, n_control (optional), end, region_id, gwas_tsv_basename.\n", "import csv\n", "jobs = list(csv.DictReader(open(f\"{cwd}/{manifest_name}.manifest.tsv\"), delimiter='\\t'))\n", "input: for_each = 'jobs'\n", @@ -242,6 +243,9 @@ " ${'--impute' if impute else ''} \\\n", " ${('--qc-args ' + repr(qc_args)) if qc_args else ''} \\\n", " ${'--allele-flip-kriging' if allele_flip_kriging else ''} \\\n", + " --effective-n ${'TRUE' if effective_n else 'FALSE'} \\\n", + " ${('--n-case ' + _jobs['n_case']) if _jobs.get('n_case','') not in ('', 'NA') else ''} \\\n", + " ${('--n-control ' + _jobs['n_control']) if _jobs.get('n_control','') not in ('', 'NA') else ''} \\\n", " --output ${_output}" ] }, diff --git a/code/script/pecotmr_integration/gwas_rss_manifest.R b/code/script/pecotmr_integration/gwas_rss_manifest.R index 5ec7621e3..39ee0bdb3 100644 --- a/code/script/pecotmr_integration/gwas_rss_manifest.R +++ b/code/script/pecotmr_integration/gwas_rss_manifest.R @@ -10,6 +10,7 @@ # # Inputs (zero or more of each input pair; at least one of each pair): # --gwas-meta Optional. Columns: study_id, path, +# n_case, n_control (optional; study-level fallback), # column_mapping (optional). Relative # paths resolve against the meta file's # own directory. @@ -24,7 +25,7 @@ # --output Output manifest path. # # Output TSV columns: -# study_id, gwas_tsv, column_mapping, chr, start, end, region_id, +# study_id, gwas_tsv, column_mapping, n_case, n_control, chr, start, end, region_id, # gwas_tsv_basename # # `region_id` is the SoS-safe sanitised region label (`:` and `-` replaced @@ -61,6 +62,8 @@ source(file.path(.d, "manifest_common.R")) studies <- data.frame(study_id = character(0), gwas_tsv = character(0), column_mapping = character(0), + n_case = numeric(0), + n_control = numeric(0), stringsAsFactors = FALSE) seenStudies <- character(0) @@ -75,7 +78,11 @@ if (nzchar(argv$gwas_meta) && argv$gwas_meta != ".") { paste(missing, collapse = ", "), " (got: ", paste(names(meta), collapse = ", "), ").") metaDir <- dirname(normalizePath(argv$gwas_meta)) - hasCm <- "column_mapping" %in% names(meta) + hasCm <- "column_mapping" %in% names(meta) + # Optional study-level case/control counts: the fallback source of effective + # sample size when the sumstats has no per-variant n_case/n_control columns. + hasNca <- "n_case" %in% names(meta) + hasNco <- "n_control" %in% names(meta) for (i in seq_len(nrow(meta))) { sid <- as.character(meta$study_id[[i]]) tsv <- as.character(meta$path[[i]]) @@ -83,12 +90,15 @@ if (nzchar(argv$gwas_meta) && argv$gwas_meta != ".") { cm <- if (hasCm) as.character(meta$column_mapping[[i]]) else "" if (!is.na(cm) && nzchar(cm) && !startsWith(cm, "/")) cm <- file.path(metaDir, cm) + nca <- if (hasNca) suppressWarnings(as.numeric(meta$n_case[[i]])) else NA_real_ + nco <- if (hasNco) suppressWarnings(as.numeric(meta$n_control[[i]])) else NA_real_ if (sid %in% seenStudies) stop("Duplicate study_id in --gwas-meta: ", sid) seenStudies <- c(seenStudies, sid) studies <- rbind(studies, data.frame(study_id = sid, gwas_tsv = tsv, column_mapping = if (is.na(cm)) "" else cm, + n_case = nca, n_control = nco, stringsAsFactors = FALSE)) } } @@ -108,7 +118,8 @@ for (item in tsvItems) { seenStudies <- c(seenStudies, sid) studies <- rbind(studies, data.frame(study_id = sid, gwas_tsv = tsv, - column_mapping = "", stringsAsFactors = FALSE)) + column_mapping = "", n_case = NA_real_, n_control = NA_real_, + stringsAsFactors = FALSE)) } if (nrow(studies) == 0L) stop("No GWAS inputs supplied (give --gwas-meta and/or --gwas-tsv-list).") @@ -163,6 +174,8 @@ for (i in seq_len(nrow(studies))) { study_id = studies$study_id[[i]], gwas_tsv = studies$gwas_tsv[[i]], column_mapping = studies$column_mapping[[i]], + n_case = studies$n_case[[i]], + n_control = studies$n_control[[i]], chr = chr, start = s, end = e, region_id = region_id, gwas_tsv_basename = tools::file_path_sans_ext( diff --git a/code/script/pecotmr_integration/gwas_sumstats_construct.R b/code/script/pecotmr_integration/gwas_sumstats_construct.R index 1304c46db..25562f84c 100644 --- a/code/script/pecotmr_integration/gwas_sumstats_construct.R +++ b/code/script/pecotmr_integration/gwas_sumstats_construct.R @@ -111,6 +111,9 @@ parser <- add_argument(parser, "--impute", parser <- add_argument(parser, "--allele-flip-kriging", help = "Enable kriging allele-flip QC (summaryStatsQc alleleFlipKriging): sign-flip allele-switched z (logLR>2 & |z|>2) in place, retaining the variants. Off by default.", flag = TRUE) +parser <- add_argument(parser, "--effective-n", + help = "TRUE/FALSE. For case/control GWAS (n_case + n_control columns), use the effective sample size 4/(1/n_case + 1/n_control) as N (summaryStatsQc effectiveN). Default TRUE; FALSE keeps the raw N column (or total when absent).", + type = "character", default = "TRUE") parser <- add_argument(parser, "--qc-args", help = "JSON object of extra named kwargs for summaryStatsQc()", type = "character", default = "") @@ -185,12 +188,12 @@ qc_extra <- if (nzchar(argv$qc_args) && argv$qc_args != "." && # 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", "alleleFlipKriging")) + "pipCutoffToSkip", "alleleFlipKriging", "effectiveN")) if (length(clash) > 0L) stop("--qc-args sets ", paste(clash, collapse = ", "), " which is also controlled by a dedicated flag (--qc-method / ", "--impute / --maf / --skip-region / --pip-cutoff-to-skip / ", - "--allele-flip-kriging). Pass it via the dedicated flag.") + "--allele-flip-kriging / --effective-n). 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). @@ -286,11 +289,15 @@ gss_out <- if (argv$skip_qc) { message("--skip-qc set; serialising raw GwasSumStats without summaryStatsQc().") gss } else { + effN <- as.logical(argv$effective_n) + if (is.na(effN)) + stop("--effective-n must be TRUE or FALSE (got: ", argv$effective_n, ")") qc_call_args <- c(list(gss, zMismatchQc = qc_method, impute = isTRUE(argv$impute), pipCutoffToSkip = argv$pip_cutoff_to_skip, - mafCutoff = argv$maf), + mafCutoff = argv$maf, + effectiveN = effN), if (isTRUE(argv$allele_flip_kriging)) list(alleleFlipKriging = TRUE) else list(), if (!is.null(skip_region_vec)) @@ -318,6 +325,7 @@ if (!argv$skip_qc) { seg <- c( if (!is.null(a$variantsIn) && !is.null(a$variantsOut)) sprintf("variants %s->%s", a$variantsIn, a$variantsOut), + if (!is.null(a$nSource) && !is.na(a$nSource)) paste0("N=", a$nSource), if (!is.null(cf$mafDropped) && cf$mafDropped > 0L) paste0("maf-drop ", cf$mafDropped), if (!is.null(cf$nDropped) && cf$nDropped > 0L) paste0("N-drop ", cf$nDropped), if (!is.null(a$krigingFlipped)) paste0("kriging-flip ", a$krigingFlipped),