Skip to content

update for sldsc processing #1388

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
danielnachun merged 1 commit into from
Jul 8, 2026
Merged

update for sldsc processing #1388

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
122 changes: 39 additions & 83 deletions code/SoS/enrichment/sldsc_enrichment.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,13 +50,14 @@
"\n",
"**Stage 2 - pecotmr post-processing.** A single `pecotmr::sldscPostprocessingPipeline` call consumes all polyfun outputs: it extracts $\\tau$, $E$, $h^2_g$, EnrichStat p-value and per-block jackknife $\\tau$ values; computes $sd_C$ and $M_{\\mathrm{ref}}$ over the regression's MAF-cutoff SNP set; standardizes $\\tau \\to \\tau^*$ for single and joint modes; auto-detects binary vs continuous annotations; and runs a DerSimonian-Laird random-effects meta-analysis across traits, producing three meta tables ($\\tau^*$ cross-type comparable, $E$ within-binary, EnrichStat within-binary). Output is an R list with `per_trait` and `meta` entries.\n",
"\n",
"**Stage 3 - subset meta-analysis.** `pecotmr::meta_sldsc_random` re-runs the meta on a trait subset without re-running the regression (lightweight, interactive):\n",
"**Stage 3 - subset meta-analysis.** `pecotmr::metaSldscRandom` re-runs the meta on a trait subset without re-running the regression (lightweight, interactive):\n",
"\n",
"```r\n",
"res <- readRDS(\"sldsc_results.rds\")\n",
"neuro <- c(\"AD_GWAX\", \"PD_meta\", \"ALS_meta\")\n",
"meta_neuro_taustar <- pecotmr::meta_sldsc_random(\n",
" res$per_trait[neuro], category = \"my_target_anno\", quantity = \"tau_star\"\n",
"view_neuro <- pecotmr:::.sldscViewForMeta(res$per_trait[neuro], \"single\")\n",
"meta_neuro_taustar <- pecotmr::metaSldscRandom(\n",
" view_neuro, category = \"my_target_anno\", quantity = \"tauStar\"\n",
")\n",
"```\n"
]
Expand Down Expand Up @@ -178,7 +179,7 @@
"- Whether each annotation is binary or continuous — **(pecotmr: `is_binary_sldsc_annot`)**.\n",
"- $\\tau^*_C$ point estimate and per-block $\\tau^*_C$ — **(pecotmr: `standardize_sldsc_trait`)**.\n",
"- EnrichStat point estimate and its standard error (formula below) — **(pecotmr: `standardize_sldsc_trait`)**.\n",
"- DerSimonian-Laird random-effects meta-analysis of $\\tau^*_C$, $E_C$, or EnrichStat across traits — **(pecotmr: `meta_sldsc_random`)**.\n",
"- DerSimonian-Laird random-effects meta-analysis of $\\tau^*_C$, $E_C$, or EnrichStat across traits — **(pecotmr: `metaSldscRandom`)**.\n",
"\n",
"The top-level entry point `pecotmr::sldscPostprocessingPipeline` orchestrates all of the above.\n",
"\n",
Expand Down Expand Up @@ -234,7 +235,7 @@
"\n",
"### Meta-Analysis across Traits (Random Effects) — (pecotmr)\n",
"\n",
"DerSimonian-Laird random-effects meta-analysis of per-annotation estimates across traits, implemented in `pecotmr::meta_sldsc_random` (which delegates the numerics to `rmeta::meta.summaries(..., method = \"random\")`):\n",
"DerSimonian-Laird random-effects meta-analysis of per-annotation estimates across traits, implemented in `pecotmr::metaSldscRandom` (which delegates the numerics to `rmeta::meta.summaries(..., method = \"random\")`):\n",
"\n",
"$$\\hat\\theta_{\\mathrm{meta}} \\;=\\; \\frac{\\sum_i w_i\\, \\hat\\theta_i}{\\sum_i w_i}, \\qquad SE_{\\mathrm{meta}} \\;=\\; \\sqrt{\\frac{1}{\\sum_i w_i}}, \\qquad w_i \\;=\\; \\frac{1}{SE_i^2 + \\hat\\sigma^2}$$\n",
"\n",
Expand Down Expand Up @@ -544,16 +545,17 @@
"kernel": "SoS"
},
"source": [
"## Step 4. `Subset Meta-Analysis (`pecotmr::meta_sldsc_random`)` (optional)\n",
"## Step 4. `Subset Meta-Analysis (`pecotmr::metaSldscRandom`)` (optional)\n",
"\n",
"The default meta in Step 2 pools all traits the user supplied. To re-run the meta on a subset (e.g., neurodegenerative traits only, or autoimmune traits only) without re-running the regression layer:\n",
"\n",
"\n",
"```r\n",
"res <- readRDS(\"sldsc_results.rds\")\n",
"neuro <- c(\"AD_GWAX\", \"PD_meta\", \"ALS_meta\")\n",
"meta_neuro_taustar <- pecotmr::meta_sldsc_random(\n",
" res$per_trait[neuro], category = \"my_target_anno\", quantity = \"tau_star\"\n",
"view_neuro <- pecotmr:::.sldscViewForMeta(res$per_trait[neuro], \"single\")\n",
"meta_neuro_taustar <- pecotmr::metaSldscRandom(\n",
" view_neuro, category = \"my_target_anno\", quantity = \"tauStar\"\n",
")\n",
"```\n",
"\n",
Expand Down Expand Up @@ -688,6 +690,8 @@
"parameter: annotation_name = str\n",
"parameter: python_exec = \"python\" # e.g. \"/home/you/.conda/envs/polyfun/bin/python\"\n",
"parameter: polyfun_path = path # e.g. \"/home/you/tools/polyfun\"\n",
"# Directory holding the modular R wrapper scripts (code/script).\n",
"parameter: modular_script_dir = path('code/script')\n",
"\n",
"# MAF cutoff for sLDSC. Default 0.05 per sLDSC recommendation (rare-variant LD is unstable\n",
"# and HapMap3-style regression weights are common-variant by construction).\n",
Expand Down Expand Up @@ -1304,12 +1308,14 @@
"outputs": [],
"source": [
"[postprocess]\n",
"# Post-processing of polyfun outputs via pecotmr::sldscPostprocessingPipeline.\n",
"# Reads .results / .log / .part_delete for all traits in `traits_file`, both\n",
"# single-target and (when present) joint-target runs, computes Gazal-style\n",
"# Post-processing of polyfun outputs via pecotmr::sldscPostprocessingPipeline,\n",
"# driven by the sldsc_postprocess.R wrapper. The wrapper reads .results / .log /\n",
"# .part_delete for all traits in `traits_file` (both single-target and, when\n",
"# present, joint-target runs), the target .annot.gz files, and the reference\n",
"# .frq files, bundles them into a pecotmr SldscData object, computes Gazal-style\n",
"# tau*, EnrichStat with back-solved jackknife SE, and runs the default\n",
"# DerSimonian-Laird random-effects meta across all supplied traits. Writes\n",
"# one RDS containing per-trait tables and three meta tables (tau*, E, EnrichStat).\n",
"# DerSimonian-Laird random-effects meta across all supplied traits. Writes one\n",
"# RDS with per-trait tables and the meta tables (tauStar, enrichment, enrichstat).\n",
"\n",
"parameter: traits_file = path() # text file: one trait sumstats filename per line\n",
"parameter: heritability_cwd = path() # parent directory of [get_heritability] outputs (contains <annotation_name>_single_<i>/ subdirs and optionally <annotation_name>_joint/)\n",
Expand All @@ -1324,52 +1330,18 @@
"output: f\"{cwd:a}/{annotation_name}.sldsc_postprocess.rds\"\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads\n",
"\n",
"R: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
" library(pecotmr)\n",
"\n",
" traits <- readLines(\"${traits_file}\")\n",
" target_cats <- c(${\",\".join('\"%s\"' % c for c in target_categories)})\n",
" target_lab <- c(${\",\".join('\"%s\"' % c for c in target_categories_label)})\n",
"\n",
" # Auto-detect single-target and joint-target output directories.\n",
" her_root <- \"${heritability_cwd}\"\n",
" all_subdirs <- list.dirs(her_root, recursive = FALSE)\n",
" single_pattern <- paste0(\"^\", \"${annotation_name}\", \"_single_([0-9]+)$\")\n",
" joint_name <- paste0(\"${annotation_name}\", \"_joint\")\n",
" single_dirs <- all_subdirs[grepl(single_pattern, basename(all_subdirs))]\n",
" single_indices <- as.integer(sub(single_pattern, \"\\\\1\", basename(single_dirs)))\n",
" single_dirs <- single_dirs[order(single_indices)]\n",
" joint_dir <- file.path(her_root, joint_name)\n",
" has_joint <- dir.exists(joint_dir)\n",
"\n",
" message(sprintf(\"Detected %d single-target dirs%s\",\n",
" length(single_dirs),\n",
" if (has_joint) \"; joint-target dir present\" else \"; no joint-target dir\"))\n",
"\n",
" # Build per-trait prefix maps. Each trait's polyfun output is at <dir>/<trait>\n",
" # (polyfun appends .results / .log / .part_delete).\n",
" trait_single_prefixes <- lapply(traits, function(t) file.path(single_dirs, t))\n",
" names(trait_single_prefixes) <- traits\n",
"\n",
" if (has_joint) {\n",
" trait_joint_prefix <- setNames(file.path(joint_dir, traits), traits)\n",
" } else {\n",
" trait_joint_prefix <- setNames(rep(NA_character_, length(traits)), traits)\n",
" }\n",
"\n",
" res <- sldscPostprocessingPipeline(\n",
" trait_single_prefixes = trait_single_prefixes,\n",
" trait_joint_prefix = trait_joint_prefix,\n",
" target_anno_dir = \"${target_anno_dir}\",\n",
" frqfile_dir = \"${frqfile_dir}\",\n",
" plink_name = \"${plink_name}\",\n",
" maf_cutoff = ${maf_cutoff},\n",
" target_categories = if (length(target_cats) > 0) target_cats else NULL,\n",
" target_labels = if (length(target_lab) > 0) target_lab else NULL\n",
" )\n",
"\n",
" saveRDS(res, \"${_output[0]}\")\n",
" message(\"S-LDSC post-processing complete; results written to ${_output[0]}\")\n"
"bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
" Rscript ${modular_script_dir}/pecotmr_integration/sldsc_postprocess.R \\\n",
" --traits-file ${traits_file} \\\n",
" --heritability-cwd ${heritability_cwd} \\\n",
" --annotation-name ${annotation_name} \\\n",
" --target-anno-dir ${target_anno_dir} \\\n",
" --frqfile-dir ${frqfile_dir} \\\n",
" --plink-name ${plink_name} \\\n",
" --maf-cutoff ${maf_cutoff} \\\n",
" ${('--target-categories ' + ','.join(target_categories)) if target_categories else ''} \\\n",
" ${('--target-categories-label ' + ','.join(target_categories_label)) if target_categories_label else ''} \\\n",
" --output ${_output}\n"
]
},
{
Expand All @@ -1382,7 +1354,8 @@
"source": [
"[meta_subset]\n",
"# Optional: re-run random-effects meta on a user-defined subset of traits, using\n",
"# the cached per-trait standardized results from [postprocess]. No regression rerun.\n",
"# the cached per-trait standardized results from [postprocess] via the\n",
"# sldsc_meta_subset.R wrapper. No regression rerun.\n",
"\n",
"parameter: postprocess_rds = path() # output of [postprocess]\n",
"parameter: subset_traits_file = path() # text file: one trait id per line, subset of those passed to [postprocess]\n",
Expand All @@ -1396,29 +1369,12 @@
"output: f\"{cwd:a}/{annotation_name}.{subset_name}.meta.rds\"\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads\n",
"\n",
"R: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
" library(pecotmr)\n",
"\n",
" res <- readRDS(\"${postprocess_rds}\")\n",
" subset_traits <- readLines(\"${subset_traits_file}\")\n",
" target_cats <- c(${\",\".join([f'\"{c}\"' for c in target_categories])})\n",
" if (length(target_cats) == 0) target_cats <- res$params$target_categories\n",
"\n",
" subset_per_trait <- res$per_trait[subset_traits]\n",
"\n",
" # Map wide names (tau_star_single/joint) to bare names meta_sldsc_random expects.\n",
" view_single <- pecotmr:::.sldsc_view_for_meta(subset_per_trait, \"single\")\n",
" view_joint <- pecotmr:::.sldsc_view_for_meta(subset_per_trait, \"joint\")\n",
"\n",
" out <- list(\n",
" tau_star_single = setNames(lapply(target_cats, function(c) meta_sldsc_random(view_single, c, \"tau_star\")), target_cats),\n",
" tau_star_joint = setNames(lapply(target_cats, function(c) meta_sldsc_random(view_joint, c, \"tau_star\")), target_cats),\n",
" enrichment = setNames(lapply(target_cats, function(c) meta_sldsc_random(view_single, c, \"enrichment\")), target_cats),\n",
" enrichstat = setNames(lapply(target_cats, function(c) meta_sldsc_random(view_single, c, \"enrichstat\")), target_cats)\n",
" )\n",
"\n",
" saveRDS(out, \"${_output[0]}\")\n",
" message(\"Subset meta complete; results written to ${_output[0]}\")"
"bash: expand = \"${ }\", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'\n",
" Rscript ${modular_script_dir}/pecotmr_integration/sldsc_meta_subset.R \\\n",
" --postprocess-rds ${postprocess_rds} \\\n",
" --subset-traits-file ${subset_traits_file} \\\n",
" ${('--target-categories ' + ','.join(target_categories)) if target_categories else ''} \\\n",
" --output ${_output}\n"
]
}
],
Expand Down
72 changes: 72 additions & 0 deletions code/script/pecotmr_integration/sldsc_meta_subset.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
#!/usr/bin/env Rscript
# sldsc_meta_subset.R
#
# Re-run the DerSimonian-Laird random-effects meta on a user-defined subset
# of traits, reusing the cached per-trait standardized results from
# sldsc_postprocess.R. No regression rerun -- this only re-meta's the
# already-standardized per-trait tables via pecotmr::metaSldscRandom().
#
# Inputs:
# --postprocess-rds <RDS> Output of sldsc_postprocess.R.
# --subset-traits-file <txt> One trait id per line; a subset of the traits
# passed to sldsc_postprocess.R.
# --target-categories c1,c2 Target annotation names to meta on. When empty,
# uses params$target_categories from the RDS.
# (If sldsc_postprocess.R ran with display
# labels, that field already holds the labels.)
# --output <RDS> Subset meta result RDS.

suppressPackageStartupMessages({
library(argparser)
library(pecotmr)
})

parser <- arg_parser("Re-run sLDSC random-effects meta on a subset of traits")
parser <- add_argument(parser, "--postprocess-rds",
help = "Output RDS of sldsc_postprocess.R",
type = "character")
parser <- add_argument(parser, "--subset-traits-file",
help = "Text file: one trait id per line",
type = "character")
parser <- add_argument(parser, "--target-categories",
help = "Comma-separated target annotation names",
type = "character", default = "")
parser <- add_argument(parser, "--output",
help = "Output subset-meta RDS path",
type = "character")
argv <- parse_args(parser)

res <- readRDS(argv$postprocess_rds)

subsetTraits <- readLines(argv$subset_traits_file)
subsetTraits <- subsetTraits[nzchar(trimws(subsetTraits))]

targetCats <- if (nzchar(argv$target_categories)) {
trimws(strsplit(argv$target_categories, ",", fixed = TRUE)[[1L]])
} else res$params$target_categories

missingTraits <- setdiff(subsetTraits, names(res$per_trait))
if (length(missingTraits) > 0L)
stop("--subset-traits-file names traits absent from --postprocess-rds: ",
paste(missingTraits, collapse = ", "))

subsetPerTrait <- res$per_trait[subsetTraits]

# Map the wide per-trait columns (tauStarSingle/tauStarJoint, ...) to the bare
# names metaSldscRandom() expects.
viewSingle <- pecotmr:::.sldscViewForMeta(subsetPerTrait, "single")
viewJoint <- pecotmr:::.sldscViewForMeta(subsetPerTrait, "joint")

out <- list(
tau_star_single = setNames(lapply(targetCats, function(cat)
metaSldscRandom(viewSingle, cat, "tauStar")), targetCats),
tau_star_joint = setNames(lapply(targetCats, function(cat)
metaSldscRandom(viewJoint, cat, "tauStar")), targetCats),
enrichment = setNames(lapply(targetCats, function(cat)
metaSldscRandom(viewSingle, cat, "enrichment")), targetCats),
enrichstat = setNames(lapply(targetCats, function(cat)
metaSldscRandom(viewSingle, cat, "enrichstat")), targetCats))

dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE)
saveRDS(out, argv$output)
message("Subset meta complete; results written to ", argv$output)
Loading