rss_analysis: effective sample size for case/control GWAS

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}"
    ]
   },

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 <TSV>           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 <TSV>              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,20 +78,27 @@ 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]])
     if (!startsWith(tsv, "/")) tsv <- file.path(metaDir, tsv)
     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(

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),