Stratified LD Score Regression (S-LDSC) Enrichment#
Per-SNP heritability enrichment of genome annotations from GWAS summary statistics alone.
Overview#
Heritability is not spread evenly over the genome. This module asks whether a functional category – a chromatin state, a set of regulatory elements, an xQTL annotation – carries more heritability than its share of SNPs would predict, using GWAS summary statistics alone.
The reasoning is short. A variant’s association statistic reflects every causal variant it tags, so under real polygenic signal a SNP that tags more variation has a higher expected chi-square, whereas confounding such as population stratification inflates chi-square independently of LD; regressing chi-square on LD score separates the two (Bulik-Sullivan et al. 2015). Computing a separate LD score per annotation category turns that one slope into one slope per category, each estimating the per-SNP heritability that category contributes (Finucane et al. 2015). Categories may overlap, so every regression conditions on a baseline-LD model – without it an estimate reflects LD-related architecture the annotation merely correlates with.
Two numbers come out, for different reasons. Enrichment, the proportion of heritability over the proportion of SNPs, reads easily but is defined only for binary annotations and is not comparable across annotations of different size. Standardized tau (\(\tau^*\)) is the per-SNP effect of a one-standard-deviation increase in the annotation, conditional on every other annotation in the model; it is comparable across annotations and across traits and extends to continuous annotations (Gazal et al. 2017). Both are reported.
The regression itself comes from polyfun, not the
original LDSC in bulik/ldsc; standardization and the across-trait meta-analysis are done in
pecotmr.
When to run it. Run this module after GWAS summary statistics and the required LD-score annotations are available to estimate annotation-specific heritability enrichment and its uncertainty.
Pipeline stages#
The pipeline runs in three stages.
Stage 1 – polyfun. Three SoS workflows wrap polyfun: make_annotation_files_ldscore converts
target annotations into polyfun .annot.gz and computes LD scores, munge_sumstats_polyfun
preprocesses each GWAS into LDSC format, and get_heritability runs polyfun’s ldsc once per
(trait, target directory) pair.
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; standardizes \(\tau \to \tau^*\) for both single- and joint-tau modes; auto-detects
binary versus continuous annotations; and runs a DerSimonian-Laird random-effects meta-analysis across
traits. The result is an R list with per_trait and meta entries.
Stage 3 – subset meta-analysis. pecotmr::metaSldscRandom re-runs the meta over a trait subset
from the cached .sldsc_postprocess.rds, without re-running the regression.
Model#
Under a polygenic assumption, in which effect sizes for variants are drawn independently from distributions with variance proportional to \(1/(p(1-p))\) where \(p\) is the minor allele frequency (MAF), the expected \(\chi^2\) statistic of variant \(j\) is:
where \(N\) is the sample size; \(M\) is the number of SNPs, so that \(h^2/M\) is the average heritability per SNP; \(a\) measures the contribution of confounding biases such as cryptic relatedness and population stratification; and \(\ell_j = \sum_k r^2_{jk}\) is the LD Score of variant \(j\), which measures the amount of genetic variation tagged by \(j\). A full derivation is given in the Supplementary Note of Bulik-Sullivan et al. (2015); an alternative derivation appears in the Supplementary Note of Zhu and Stephens (2017) AoAS.
Equation (1) computes SNP heritability from summary statistics alone, with no individual-level genotypes.
Let \(a_{jC}\) denote the value of annotation \(C\) at SNP \(j\):
Binary annotation (e.g. an indicator for “in enhancer”, “in exon”, “in cell-type-specific peak”): \(a_{jC} \in \{0, 1\}\).
Continuous annotation (e.g. gene-specificity score, conservation score, continuous epigenomic signal): \(a_{jC} \in \mathbb{R}\).
Under a polygenic model the per-SNP heritability for SNP \(j\) is
and the expected \(\chi^2\) statistic of SNP \(j\) is
where \(\ell(j, C) = \sum_k a_{kC}\, r^2_{jk}\) is the partitioned LD score of SNP \(j\) with respect to annotation \(C\), and \(a\) measures confounding bias. Equation (2) allows joint estimation of all \(\tau_C\) via a (computationally simple) multiple regression of \(\chi^2_j\) against \(\ell(j, C)\).
Interpretation of \(\tau_C\):
Binary \(C\): \(\tau_C\) is the additive increase in per-SNP heritability for SNPs in category \(C\), on top of the contributions from any other annotations they belong to.
Continuous \(C\): \(\tau_C\) is the additive change in per-SNP heritability per unit increase in the value of annotation \(C\).
Tau estimation and enrichment analysis#
Notation#
For each annotation \(C\):
\(\pi^{h^2}_C\) = proportion of trait heritability \(h^2_g\) assigned to annotation \(C\).
\(\pi^{M}_C\) = proportion of (effective) SNPs in annotation \(C\). For binary annotations this is \(M_C / M_{\mathrm{ref}}\); for continuous annotations it is the share of total annotation weight in \(C\).
Quantities from the regression layer (polyfun)#
Solving Equation (2) jointly across annotations, with 200-block genomic jackknife for inference, is performed by polyfun’s ldsc. Each polyfun run yields, per annotation:
\(\tau_C\) and its standard error — (polyfun).
\(\pi^{h^2}_C\) and \(\pi^{M}_C\) — (polyfun).
\(E_C = \pi^{h^2}_C / \pi^{M}_C\) and its standard error — (polyfun).
The p-value of the differential per-SNP heritability test (defined below) — (polyfun), computed internally with the full coefficient covariance matrix.
Also obtained, per run:
The total trait heritability \(h^2_g\) — (polyfun).
The 200-block jackknife delete-values of \(\tau_C\) — (polyfun).
Quantities from the post-processing layer (pecotmr)#
From the polyfun outputs above plus the reference panel, the post-processing layer computes:
\(sd_C\) — per-annotation standard deviation over MAF \(>\) cutoff SNPs — (pecotmr:
compute_sldsc_annot_sd).\(M_{\mathrm{ref}}\) — reference SNP count at the MAF cutoff — (pecotmr:
compute_sldsc_M_ref).Whether each annotation is binary or continuous — (pecotmr:
is_binary_sldsc_annot).\(\tau^*_C\) point estimate and per-block \(\tau^*_C\) — (pecotmr:
standardize_sldsc_trait).EnrichStat point estimate and its standard error (formula below) — (pecotmr:
standardize_sldsc_trait).DerSimonian-Laird random-effects meta-analysis of \(\tau^*_C\), \(E_C\), or EnrichStat across traits — (pecotmr:
metaSldscRandom).
The top-level entry point pecotmr::sldscPostprocessingPipeline orchestrates all of the above.
Standardized tau (\(\tau^*\)) — (pecotmr)#
\(\tau_C\) has units that depend on the scale of the annotation and on the total heritability of the trait, so raw \(\tau\) is not directly comparable across annotations or across traits; the standardized version is computed instead (Gazal et al. 2017)
interpreted as the additive change in per-SNP heritability associated with a 1 standard deviation increase in annotation \(C\), divided by the average per-SNP heritability across all SNPs. \(\tau^*_C\) is dimensionless and comparable across annotations and across traits. In a joint multi-annotation regression it is the independent contribution of annotation \(C\) after controlling for overlapping effects of the others.
Here \(sd_C\) is the standard deviation of annotation \(C\) across reference SNPs (MAF \(>\) cutoff), \(M_{\mathrm{ref}}\) is the count of those SNPs, and \(h^2_g\) is the trait heritability. Applying the same scaling to each of the 200 jackknife blocks yields per-block \(\tau^*_C\) values; their sample variance gives the jackknife standard error $\(SE^{\text{jackknife}}(\tau^*_C) \;=\; \sqrt{\,\tfrac{(B-1)^2}{B}\, \mathrm{Var}_b(\tau^*_{C,(b)})\,}\)\( with \)B = 200$, used as the per-trait input to cross-trait meta-analysis.
Differential per-SNP heritability (“EnrichStat”) — (polyfun + pecotmr)#
To test whether the per-SNP heritability inside annotation \(C\) differs from outside it (Finucane et al. 2015):
The point-estimate p-value of this test is computed by polyfun internally using the full coefficient covariance and reported as Enrichment_p. Its standard error is recovered from the reported p-value:
This per-trait point + SE is the input to cross-trait meta-analysis.
Reporting: binary vs. continuous annotations — (pecotmr)#
The estimation machinery applies to both annotation types, but the headline quantity to report within each type differs.
For a binary annotation (e.g. enhancer indicator, exon, in/out of a cell-type peak), \(\pi^{M}_C = M_C / M_{\mathrm{ref}}\) has a direct interpretation and \(E_C\) reads as “the category explains \(E_C\)-fold more heritability than its share of SNPs.” The within-type headline quantities are therefore \(E_C\) and the EnrichStat p-value; \(\tau^*_C\) is reported alongside.
For a continuous annotation (e.g. gene-specificity score, conservation score, continuous epigenomic signal), \(E_C\) depends on the scale of the annotation: rescaling the annotation by a constant changes \(E_C\) even though the underlying biology is unchanged. The within-type headline quantities are therefore \(\tau^*_C\) and its p-value; \(E_C\) is reported alongside but should not be interpreted for continuous annotations.
The pipeline determines whether an annotation is binary by inspecting whether its values lie in \(\{0, 1\}\) and selects the appropriate within-type headline statistic automatically (pecotmr).
From the official LDSC tutorial (Partitioned Heritability from Continuous Annotations):
“Enrichment is (Prop. heritability) / (Prop. SNPs). These outputs make sense only for binary annotations. Do not try to interpret them for continuous annotations. Using
--print-coefficientsoutputs the regression coefficients and corresponding standard errors and Z score for each annotation. These coefficients measure the additional contribution of one annotation to the model and are interpretable for both binary and continuous annotations.”The pipeline always passes
--print-coefficientsto polyfun for this reason.
Cross-type comparison: always use \(\tau^*_C\) — (pecotmr)#
For an apple-to-apple comparison across binary and continuous annotations — ranking annotations on a single axis, meta-analyzing a mixed set, or reporting a leaderboard that pools both types — use \(\tau^*_C\). The standardization in Gazal et al. (2017) was designed for exactly this purpose: \(sd_C = \sqrt{p(1-p)}\) for a binary annotation (where \(p\) is the proportion in the category) and \(sd_C = \) empirical standard deviation for a continuous annotation, so the resulting \(\tau^*_C\) is dimensionless and has the same interpretation in both cases — additive change in per-SNP heritability per 1 SD increase in the annotation, normalized by the average per-SNP heritability. \(E_C\) does not have this property and must not be compared across types.
The pipeline emits both \(E_C\) and \(\tau^*_C\) for every annotation, with the binary/continuous flag, so callers can pick the right column for the comparison they are making.
Joint analysis — (polyfun runs the regression; pecotmr standardizes both modes)#
For joint analysis (multiple annotations fit together), both \(\tau\) and \(E\) are conditional on the other annotations in the model; joint \(\tau^*_C\) is reported as the independent contribution of annotation \(C\) after controlling for the others. The annotation-prep step exposes two independent toggles, compute_single and compute_joint (both default True), so the user can produce the \(N\) single-target outputs, the joint output, or both in one invocation. With both defaults the post-processing layer reads all \(N+1\) regression outputs per trait and presents single + joint side-by-side. When the joint subset is decided after looking at single-target results (exploratory \(\rightarrow\) conditional workflow), the user runs the annotation-prep step a second time with compute_single=False on the curated subset.
Meta-Analysis across Traits (Random Effects) — (pecotmr)#
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")):
where \(\hat\theta_i\) is the per-trait estimate and \(SE_i\) its standard error:
For \(\tau^*_C\) meta: \(SE_i\) is the jackknife SE from the per-block \(\tau^*_C\) values.
For \(E_C\) meta: \(SE_i\) is the polyfun-reported
Enrichment_std_error.For EnrichStat meta: \(SE_i\) is the back-solved SE from polyfun’s
Enrichment_p.
Binary-annotation enrichment reporting uses a two-channel meta: the effect size and SE come from the meta on \(E_C\) (interpretable on the original enrichment-fold scale), while the p-value comes from the meta on EnrichStat (the appropriate hypothesis test). The pipeline produces a default meta over all supplied traits; users can re-run meta on any subset of traits without re-running the regression layer.
When to run it. Once GWAS summary statistics and an annotation worth testing are in hand. S-LDSC characterises a trait against a whole category; when the question is about a specific variant set, use gregor instead.
Input#
Every path-list input is a two-column index file (#id, path), one row per chromosome, so the
pipeline can fan out per chromosome. Example index files live in input/enrichment/sldsc/.
Required throughout:
--annotation-name(required): prefix for every output directory and file of this run.--cwd(defaultoutput): work directory for the analysis.--modular-script-dir(defaultcode/script): directory holding the modular R wrapper scripts.--maf-cutoff(default0.05): MAF cutoff, per the S-LDSC recommendation, since rare-variant LD is unstable and HapMap3-style regression weights are common-variant by construction.0opts out of MAF filtering (not recommended); any other value requires recomputing LD scores at that cutoff.
make_annotation_files_ldscore – build annotation and LD-score files:
--annotation-file: index of the target annotation files to test. Binary and continuous annotations are both supported and the type is auto-detected per column. Multiple target annotations may be given in one file (columnspath,path1,path2, …); single-target and joint-target analyses are then produced in the same pass.--reference-anno-file: index of the baseline annotations, typically a baseline-LD baseline-LD model in.annot.gzform: v1.1 of Gazal et al. 2017 carries 75 annotations, v2.2 carries 97. The baseline conditions every regression.--genome-ref-file: index of the PLINK-format genome reference.--chromosome: restrict to a subset of chromosomes.--snp-list: HapMap3-style SNP list. Its presence switches the polyfun call: given a list,ldsc --l2 --print-snpswrites.l2.ldscore.gz; without one,compute_ldscoreswrites.l2.ldscore.parquet. It also changes the LD-score column name thatpostprocessmatches on, so--target-categoriesmay need to be passed explicitly.--ld-wind-kb(default0) /--ld-wind-cm(default1.0): LD window. Using the cM window requires genetic-map positions in the.bimfourth column; switch to--ld-wind-kbif it holds 0.--compute-single/--no-compute-singleand--compute-joint/--no-compute-joint(both default on): whether to emit the per-annotation single-target and the joint directories. The joint directory is only emitted for two or more annotations.--score-column(default3) and--is-range/--no-is-range(default off): which column carries the annotation score, and whether the annotation file is interval- or position-based.
With --no-is-range, one row per position and an optional score column, defaulting to 1:
chr pos score
1 10001 1
1 10002 1
With --is-range, one row per interval:
chr start end score
1 10001 20001 1
1 30001 40001 1
get_heritability – run the S-LDSC regression. All LD-derived inputs (baseline LD scores, target LD
scores, regression weights, allele frequencies) must be computed against the same reference panel as
--genome-ref-file; files from different panels (1000G versus ADSP, say) must not be mixed:
--all-traits-file: list of munged summary statistics to run, one per line.--sumstat-dir: directory holding those summary statistics.--target-anno-dir/--target-anno-dirs: the LD-score directory or directories written bymake_annotation_files_ldscore. Each (trait, target directory) pair is one polyfun call.--baseline-ld-dirand--baseline-name(defaultbaseline_chr): baseline LD scores and their per-chromosome file prefix.--frqfile-dirand--plink-name(defaultADSP_chr): PLINK.frqallele frequencies, used to enforce the MAF cutoff.--weights-dirand--weight-name(defaultweights_chr): LD-score regression weights.--n-blocks(default200): jackknife blocks for the standard errors.--all-traits: run every trait found rather than the list.
munge_sumstats_polyfun – optional preprocessing, run before get_heritability:
--sumstats(required): raw GWAS summary statistics to convert to LDSC format.--n: sample size, when the file does not carry one.--min-info(default0.6),--min-maf(default0.001),--chi2-cut(default30): quality and outlier filters.--keep-hla/--no-keep-hla(default off): retain the HLA region.
postprocess – standardize tau and meta-analyse across traits:
--traits-file: traits to include.--heritability-cwd: theget_heritabilityoutput directory to read.--target-categoriesand--target-categories-label: which.resultsCategoryrows are the target annotation, and how to label them. Auto-detected unless--snp-listchanged the column name.
meta_subset – optional re-meta over a trait subset, no regression re-run:
--postprocess-rds: the cached.sldsc_postprocess.rds.--subset-traits-fileand--subset-name(required): the trait subset and its output label.--target-categories: as above.
Cluster resources: --numThreads (default 16), --job-size, --walltime, --mem.
Output#
<annotation_name>/under the LD-score directory – frommake_annotation_files_ldscore: polyfun.annot.gzfiles plus per-annotation LD scores (.l2.ldscore.parquetor.l2.ldscore.gz,.l2.M,.l2.M_5_50). One single-target directory per annotation, plus a joint directory when there is more than one. Recomputable from the annotation and reference panel, and read by every regression.<trait>.<results,log,part_delete,delete>– fromget_heritability: the S-LDSC regression output, one set per trait per target directory. The.resultsCategorycolumn carries the annotation name with a_<ref-ld-index>suffix.<annotation_name>.sldsc_postprocess.rds– frompostprocess: a single object holding per-trait tables (Gazal-style tau*, EnrichStat with back-solved jackknife SE) and three DerSimonian-Laird random-effects meta tables (tau*, E, EnrichStat).<subset_name>tables – frommeta_subset: a re-meta of the cached.sldsc_postprocess.rdsover a user-defined trait subset. Lightweight: no regression is re-run.
Each stage writes into its own --cwd, and leaves a .stdout / .stderr log and the generated shell
or R script beside its targets.
LD scores for the example annotation, under output/sldsc_ldscore/protocol_example_single_1/:
protocol_example_single_1.19.annot.gz
protocol_example_single_1.19.l2.ldscore.parquet
protocol_example_single_1.19.l2.M
protocol_example_single_1.19.l2.ldscore.parquet.log
The regression result, output/sldsc_heritability/protocol_example_single_1/sumstats.parquet.results
(first four columns of each row shown; the file also carries Enrichment, Enrichment_std_error,
Enrichment_p, Coefficient, Coefficient_std_error and Coefficient_z-score):
Category Prop._SNPs Prop._h2 Prop._h2_std_error
ANNOT_0 1.6443e-02 4.1714e-01 2.9297e-01
Coding_UCSC_common_1 5.2671e-03 1.5724e-01 1.0492e-01
postprocess and meta_subset have no example output committed: output/sldsc_postprocess/ holds only
the generated R script and its log, so no .sldsc_postprocess.rds was produced for the toy data.
Minimal Working Example#
Step 1. make_annotation_files_ldscore#
Annotation preparation and S-LDSC regression (polyfun). This step accepts a single annotation file for a single-tau analysis (one annotation as input) or several annotation files for a joint-tau analysis (multiple annotations as input).
Timing: TBD (on toy dataset)
sos run pipeline/sldsc_enrichment.ipynb make_annotation_files_ldscore \
--annotation_file input/enrichment/sldsc/colocboost_test_annotation_path.txt \
--reference_anno_file input/enrichment/sldsc/reference_annotation0.txt \
--genome_ref_file input/enrichment/sldsc/genome_reference_bfile.txt \
--annotation_name protocol_example \
--plink_name reference. --baseline_name annotations. --weight_name weights. \
--cwd output/sldsc_ldscore -j 4
Munge summary statistics (preprocessing, run before Step 2)#
Before estimating heritability, each raw GWAS summary-statistics file must be converted into the LDSC-compatible format consumed by get_heritability. Run munge_sumstats_polyfun once per trait; the munged files are then collected in the directory passed to get_heritability via --sumstat_dir.
# sos run pipeline/sldsc_enrichment.ipynb munge_sumstats_polyfun \
# --sumstats data/polyfun_new/example_data/trait_raw_sumstats.tsv \
# --n 0 \
# --min-info 0.6 \
# --min-maf 0.001 \
# --chi2-cutoff 30 \
# --cwd data/polyfun_new/example_data
Step 2. get_heritability#
Timing: TBD (on toy dataset)
sos run pipeline/sldsc_enrichment.ipynb get_heritability \
--target_anno_dirs output/sldsc_ldscore/protocol_example_single_1 \
--all_traits_file input/enrichment/sldsc/sumstats_test_all.txt \
--sumstat_dir input/enrichment/sldsc \
--baseline_ld_dir input/enrichment/sldsc \
--weights_dir input/enrichment/sldsc \
--plink_name reference. --baseline_name annotations. --weight_name weights. \
--annotation_name protocol_example \
--maf_cutoff 0 --cwd output/sldsc_heritability -j 4
Step 3. Post-processing (pecotmr) and meta-analysis#
Reads every polyfun output under --heritability-cwd and calls
pecotmr::sldscPostprocessingPipeline to produce the per-trait standardized tables and the default
random-effects meta across all traits.
--target-categories-label (same order as --target-categories) gives the target annotations
friendly names in the output – --target-categories ANNOT_1_0 ANNOT_2_0 --target-categories-label quantile_eQTL eQTL makes the target column read quantile_eQTL / eQTL rather than ANNOT_1_0 /
ANNOT_2_0. The original names stay in params$target_categories_orig; omit the flag to keep the
polyfun .results names.
Timing: TBD (on toy dataset)
sos run pipeline/sldsc_enrichment.ipynb postprocess \
--traits_file input/enrichment/sldsc/sumstats_test_all.txt \
--heritability_cwd output/sldsc_heritability \
--target_categories ANNOT_0 --target_categories_label protocol_example_annotation \
--target_anno_dir output/sldsc_ldscore/protocol_example_single_1 \
--annotation_name protocol_example \
--maf_cutoff 0 --cwd output/sldsc_postprocess -j 4
Step 4. Subset Meta-Analysis (pecotmr::metaSldscRandom) (optional)#
Step 3 pools every trait passed to postprocess. To re-meta a subset – neurodegenerative traits
only, say – without re-running the regression or the per-trait standardization, meta_subset works
from the cached .sldsc_postprocess.rds. It is light enough to run interactively:
res <- readRDS("sldsc_results.rds")
neuro <- c("AD_GWAX", "PD_meta", "ALS_meta")
view_neuro <- pecotmr:::.sldscViewForMeta(res$per_trait[neuro], "single")
meta_neuro_taustar <- pecotmr::metaSldscRandom(
view_neuro, category = "my_target_anno", quantity = "tauStar"
)
Timing: TBD (on toy dataset)
sos run pipeline/sldsc_enrichment.ipynb meta_subset \
--postprocess_rds output/sldsc_postprocess/protocol_example.sldsc_postprocess.rds \
--subset_traits_file input/enrichment/sldsc/sumstats_test_category1.txt \
--subset_name category1 --target_categories ANNOT_0 \
--annotation_name protocol_example \
--maf_cutoff 0 --cwd output/sldsc_postprocess -j 4
Command Interface#
sos run pipeline/sldsc_enrichment.ipynb -h
usage: sos run pipeline/sldsc_enrichment.ipynb
[workflow_name | -t targets] [options] [workflow_options]
workflow_name: Single or combined workflows defined in this script
targets: One or more targets to generate
options: Single-hyphen sos parameters (see "sos run -h" for details)
workflow_options: Double-hyphen workflow-specific parameters
Workflows:
make_annotation_files_ldscore
get_heritability
munge_sumstats_polyfun
postprocess
meta_subset
Global Workflow Options:
--cwd output (as path)
Path to the work directory of the analysis.
--annotation-name VAL (as str, required)
Prefix for the analysis output
--modular-script-dir code/script (as path)
Directory holding the modular R wrapper scripts
(code/script).
--maf-cutoff 0.05 (as float)
MAF cutoff for sLDSC. Default 0.05 per sLDSC
recommendation (rare-variant LD is unstable and
HapMap3-style regression weights are common-variant by
construction). Set to 0 to opt out of MAF filtering (NOT
recommended; only use if you understand the
implications). Other values would require recomputing LD
scores at that cutoff.
--annotation-file . (as path)
for make_annotation_files_ldscore workflow:
--reference-anno-file . (as path)
--genome-ref-file . (as path)
--chromosome (as list)
--snp-list . (as path)
--ld-wind-kb 0 (as int)
--ld-wind-cm 1.0 (as float)
--all-traits-file . (as path)
for get_heritability workflow. Note: all LD-derived
inputs (baseline LD scores, target LD scores, regression
weights, allele frequencies) must be computed against
the same reference panel as `genome_ref_file`. Do not
mix files derived from different reference panels (e.g.,
1000G vs ADSP).
--sumstat-dir . (as path)
--target-anno-dir . (as path)
--baseline-ld-dir . (as path)
--frqfile-dir . (as path)
--plink-name 'ADSP_chr'
--weights-dir . (as path)
--baseline-name 'baseline_chr'
--weight-name 'weights_chr'
--n-blocks 200 (as int)
--numThreads 16 (as int)
Number of threads
--job-size 1 (as int)
For cluster jobs, number commands to run per job
--walltime 12h
--mem 16G
Sections
make_annotation_files_ldscore:
Workflow Options:
--[no-]compute-single (default to True)
--- snplist (--snp_list) vs no-snplist: which polyfun
script, output format, column name, and the CM
requirement --- --snp_list given -> ldsc --l2 --print-
snps -> output .l2.ldscore.gz --snp_list absent ->
compute_ldscores -> output .l2.ldscore.parquet
LD-score column name (this is what becomes the .results
"Category" in [get_heritability], with a "_<ref-ld-
index>" suffix appended there): * compute_ldscores
ALWAYS keeps the annot column name(s): single annot
column "ANNOT" -> ldscore column "ANNOT" joint
annot columns "ANNOT_1","ANNOT_2",... ->
"ANNOT_1","ANNOT_2",... * ldsc --l2 has a quirk: with
EXACTLY ONE annotation (n_annot == 1) it HARD-CODES the
ldscore column name to "L2" and DROPS the annot's
original column name. With >=2 annotations it uses
"<annot_name>L2" ("ANNOT_1L2","ANNOT_2L2",...). => a
single-target snplist run reports "L2_0" in .results,
while a single-target no-snplist run reports "ANNOT_0".
[postprocess] auto- detects either; only matters if you
pass --target-categories explicitly. CM column
requirement for snplist: ldsc --l2 --print-snps
requires the target annot to (a) carry a "CM"
(centimorgan) column and (b) line up with the plink .bim
(same SNP set, same row order). This step handles both
internally (normalize_for_ldsc: takes CM from the .bim
4th column, re-expands the annot onto the .bim rows,
filling 0). Therefore the plink .bim files MUST carry
genetic-map (cM) positions when using --ld-wind-cm (the
default); if your .bim has 0 in the cM column, switch to
--ld-wind-kb instead.
--[no-]compute-joint (default to True)
--score-column 3 (as int)
--[no-]is-range (default to False)
get_heritability:
Workflow Options:
--target-anno-dirs paths()
--- about the ".results" Category column and the "_0 /
_1" suffix --- Each (trait, target_dir) pair is ONE
polyfun call; its `ldsc --ref-ld-chr` always gets
exactly two LD-score sources, in this order:
"<target_dir>/<target>." (index 0) ,
"<baseline_dir>/<baseline>" (index 1) With --overlap-
annot, every annotation column in the .results
"Category" is named <ldscore_column_name>_<ref-ld-
index>: index 0 = the target file -> "ANNOT_0" (no-
snplist; compute_ldscores keeps the annot col name) ->
"L2_0" (snplist + single annot; ldsc hard-codes "L2",
see below) -> "ANNOT_1_0","ANNOT_2_0" (no-snplist
joint dir, N>=2 annot cols) ->
"ANNOT_1L2_0","ANNOT_2L2_0" (snplist joint dir, N>=2 ->
"<name>L2") index 1 = the baseline file ->
"base_1","Coding_UCSC_1", ... (the 97 baseline annots)
So in this pipeline the suffix is only ever 0 (target)
or 1 (baseline); it would continue 0,1,2,... only if you
handed `ldsc --ref-ld-chr` more than two sources. (Why
ANNOT_0 vs L2_0: see the [make_annotation_files_ldscore]
header — ldsc's "n_annot == 1 -> column name 'L2'" quirk
vs compute_ldscores keeping the annot column name.)
[postprocess] auto-detects the target Category; if you
instead pass --target-categories, the names must match
this column exactly.
--all-traits (as list)
munge_sumstats_polyfun:
Workflow Options:
--sumstats VAL (as path, required)
--n 0 (as int)
--min-info 0.6 (as float)
--min-maf 0.001 (as float)
--[no-]keep-hla (default to False)
--chi2-cut 30 (as int)
postprocess:
Workflow Options:
--traits-file . (as path)
--heritability-cwd . (as path)
--target-categories (as list)
--target-categories-label (as list)
meta_subset:
Workflow Options:
--postprocess-rds . (as path)
--subset-traits-file . (as path)
--subset-name VAL (as str, required)
--target-categories (as list)
Workflow implementation#
The cells below are the pipeline definition (preserved from the original notebook): the [global] parameter block and the workflow step bodies.
[global]
# Path to the work directory of the analysis.
parameter: cwd = path('output')
# Prefix for the analysis output
parameter: annotation_name = str
# Directory holding the modular R wrapper scripts (code/script).
parameter: modular_script_dir = path('code/script')
# MAF cutoff for sLDSC. Default 0.05 per sLDSC recommendation (rare-variant LD is unstable
# and HapMap3-style regression weights are common-variant by construction).
# Set to 0 to opt out of MAF filtering (NOT recommended; only use if you understand the implications).
# Other values would require recomputing LD scores at that cutoff.
parameter: maf_cutoff = 0.05
# for make_annotation_files_ldscore workflow:
parameter: annotation_file = path()
parameter: reference_anno_file = path()
parameter: genome_ref_file = path() # with .bed
parameter: chromosome = []
parameter: snp_list = path()
parameter: ld_wind_kb = 0 # use kb if the value is provided
parameter: ld_wind_cm = 1.0 # default using ld_wind_cm
# for get_heritability workflow.
# Note: all LD-derived inputs (baseline LD scores, target LD scores, regression weights,
# allele frequencies) must be computed against the same reference panel as `genome_ref_file`.
# Do not mix files derived from different reference panels (e.g., 1000G vs ADSP).
parameter: all_traits_file = path() # txt file, each row contains all GWAS summary statistics name: e.g. CAD_META.filtered.sumstats.gz
parameter: sumstat_dir = path() # Directory containing GWAS summary statistics
parameter: target_anno_dir = path() # Directory containing target annotation files: output of ldscore
parameter: baseline_ld_dir = path() # Directory containing baseline LD score files (computed against our panel)
parameter: frqfile_dir = path() # Directory containing allele frequency files (.frq, our panel)
parameter: plink_name = "ADSP_chr"
parameter: weights_dir = path() # Directory containing LD weights (computed against our panel)
parameter: baseline_name = "baseline_chr" # Prefix of baseline annotation files
parameter: weight_name = "weights_chr" # Prefix of LD weights files
parameter: n_blocks = 200
# Number of threads
parameter: numThreads = 16
# For cluster jobs, number commands to run per job
parameter: job_size = 1
parameter: walltime = '12h'
parameter: mem = '16G'
Make Annotation File#
[make_annotation_files_ldscore]
# Annotation preparation. Takes one annotation_file with N target annotations
# and produces, in one invocation, any combination of:
# - N single-target LD-score directories (when compute_single = TRUE, default)
# - 1 joint LD-score directory containing all N (when compute_joint = TRUE
# and N >= 2, default)
#
# Outputs per chromosome <chr>:
# <cwd>/<annotation_name>_single_<i>/<annotation_name>_single_<i>.<chr>.annot.gz (i in 1..N, when compute_single)
# <cwd>/<annotation_name>_single_<i>/<annotation_name>_single_<i>.<chr>.l2.ldscore.{parquet|gz}
# <cwd>/<annotation_name>_single_<i>/<annotation_name>_single_<i>.<chr>.l2.M
# <cwd>/<annotation_name>_single_<i>/<annotation_name>_single_<i>.<chr>.l2.M_5_50 (when .frq present)
#
# <cwd>/<annotation_name>_joint/<annotation_name>_joint.<chr>.{...} (when compute_joint and N>=2)
#
# Workflows:
# - Workflow A ("all at once"): compute_single=TRUE, compute_joint=TRUE (defaults).
# Produces both, fits the case where you have already chosen the joint set.
# - Workflow B ("exploratory then conditional"):
# Step 1: compute_single=TRUE, compute_joint=FALSE.
# Run on N candidate annotations -> N single-target dirs.
# Inspect single-target results, identify K significant ones.
# Step 2: compute_single=FALSE, compute_joint=TRUE.
# Run on a NEW annotation_file with the K selected annotations
# -> 1 joint dir with the conditional model.
#
# --- snplist (--snp_list) vs no-snplist: which polyfun script, output format,
# column name, and the CM requirement ---
# --snp_list given -> ldsc --l2 --print-snps -> output .l2.ldscore.gz
# --snp_list absent -> compute_ldscores -> output .l2.ldscore.parquet
#
# LD-score column name (this is what becomes the .results "Category" in
# [get_heritability], with a "_<ref-ld-index>" suffix appended there):
# * compute_ldscores ALWAYS keeps the annot column name(s):
# single annot column "ANNOT" -> ldscore column "ANNOT"
# joint annot columns "ANNOT_1","ANNOT_2",... -> "ANNOT_1","ANNOT_2",...
# * ldsc --l2 has a quirk: with EXACTLY ONE annotation (n_annot == 1) it
# HARD-CODES the ldscore column name to "L2" and DROPS the annot's original
# column name. With >=2 annotations it uses "<annot_name>L2"
# ("ANNOT_1L2","ANNOT_2L2",...).
# => a single-target snplist run reports "L2_0" in .results, while a
# single-target no-snplist run reports "ANNOT_0". [postprocess] auto-
# detects either; only matters if you pass --target-categories explicitly.
#
# CM column requirement for snplist: ldsc --l2 --print-snps requires the
# target annot to (a) carry a "CM" (centimorgan) column and (b) line up with
# the plink .bim (same SNP set, same row order). This step handles both
# internally (normalize_for_ldsc: takes CM from the .bim 4th column, re-expands
# the annot onto the .bim rows, filling 0). Therefore the plink .bim files MUST
# carry genetic-map (cM) positions when using --ld-wind-cm (the default);
# if your .bim has 0 in the cM column, switch to --ld-wind-kb instead.
#
parameter: compute_single = True
parameter: compute_joint = True
parameter: score_column = 3
parameter: is_range = False
import csv
import os
if not (compute_single or compute_joint):
raise ValueError("[make_annotation_files_ldscore] at least one of compute_single or compute_joint must be TRUE")
def adapt_file_path(file_path, reference_file):
reference_path = os.path.dirname(reference_file)
if os.path.isfile(file_path):
return file_path
file_name = os.path.basename(file_path)
if os.path.isfile(file_name):
return file_name
file_in_ref_dir = os.path.join(reference_path, file_name)
if os.path.isfile(file_in_ref_dir):
return file_in_ref_dir
file_prefixed = os.path.join(reference_path, file_path)
if os.path.isfile(file_prefixed):
return file_prefixed
raise FileNotFoundError(f"No valid path found for file: {file_path}")
def _read_tsv(p):
with open(p) as _f:
return list(csv.DictReader(_f, delimiter="\t"))
# ---- Parse inputs and determine N (stdlib; no pandas so it runs in the base env) ----
if (str(annotation_file).endswith(('rds', 'tsv', 'txt', 'tsv.gz', 'txt.gz')) and
str(reference_anno_file).endswith('annot.gz')):
# Case 1: direct file paths (single-chromosome run). Multiple target files separated by ','.
target_files_direct = str(annotation_file).split(',')
N_targets = len(target_files_direct)
input_files = [[*target_files_direct, str(reference_anno_file), str(genome_ref_file)]]
input_chroms = [int(x) for x in chromosome] if len(chromosome) > 0 else [0]
else:
# Case 2: txt list with #id and one or more 'path' columns
_tgt = _read_tsv(annotation_file); _ref = _read_tsv(reference_anno_file); _gen = _read_tsv(genome_ref_file)
for _rows in (_tgt, _ref, _gen):
for _r in _rows:
_r["#id"] = str(_r["#id"]).replace("chr", "")
path_columns = [k for k in _tgt[0].keys() if k.startswith("path")]
N_targets = len(path_columns)
for _r in _tgt:
for _c in path_columns:
_r[_c] = adapt_file_path(_r[_c], str(annotation_file))
_ref_by = {_r["#id"]: adapt_file_path(_r["path"], str(reference_anno_file)) for _r in _ref}
_gen_by = {_r["#id"]: adapt_file_path(_r["path"], str(genome_ref_file)) for _r in _gen}
_keep = set(str(c) for c in chromosome) if len(chromosome) > 0 else None
input_chroms = []; input_files = []
for _r in _tgt: # inner-join on #id, preserving target order
_id = _r["#id"]
if _id not in _ref_by or _id not in _gen_by:
continue
if _keep is not None and _id not in _keep:
continue
input_chroms.append(_id)
input_files.append([_r[_c] for _c in path_columns] + [_ref_by[_id], _gen_by[_id]])
# ---- Determine output format ----
use_print_snps = snp_list.is_file()
ldscore_ext = "l2.ldscore.gz" if use_print_snps else "l2.ldscore.parquet"
if ld_wind_kb > 0:
ld_window_param = ld_wind_kb
ld_window_flag = "--ld-wind-kb"
else:
ld_window_param = ld_wind_cm
ld_window_flag = "--ld-wind-cm"
emit_single = compute_single
emit_joint = compute_joint and N_targets >= 2
# ---- Build per-chromosome output list ----
def chrom_outputs(chrom):
outs = []
if emit_single:
for i in range(N_targets):
name = f"{annotation_name}_single_{i+1}"
prefix = f"{cwd:a}/{name}/{name}.{chrom}"
outs += [f"{prefix}.annot.gz", f"{prefix}.{ldscore_ext}", f"{prefix}.l2.M"]
if emit_joint:
name = f"{annotation_name}_joint"
prefix = f"{cwd:a}/{name}/{name}.{chrom}"
outs += [f"{prefix}.annot.gz", f"{prefix}.{ldscore_ext}", f"{prefix}.l2.M"]
return outs
input: input_files, group_by = N_targets + 2, group_with = "input_chroms"
output: chrom_outputs(input_chroms[_index])
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bnn}'
# ---- Step A: write the .annot.gz files for this chromosome (R worker) ----
bash: expand = "${ }", stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'
Rscript ${modular_script_dir}/enrichment/make_annotation.R --step annot \
--targets ${",".join(str(p.absolute()) for p in _input[:-2])} \
--reference-anno ${_input[-2]:a} \
--bfile-prefix ${_input[-1]:na} \
--score-column ${score_column} \
${('--is-range' if is_range else '')} \
${('--use-print-snps' if use_print_snps else '')} \
${('--emit-single' if emit_single else '')} \
${('--emit-joint' if emit_joint else '')} \
--annotation-name ${annotation_name} --cwd ${cwd:a} --chrom ${input_chroms[_index]}
# ---- Step C: run polyfun's LD-score computation for each emitted annotation file (external). ----
bash: expand = "$[ ]", stderr = f'{_output[1]}.stderr', stdout = f'{_output[1]}.stdout'
set -e
chrom="$[input_chroms[_index]]"
run_polyfun() {
local annot="$1"
local out_prefix="$2"
if [ "$[str(use_print_snps)]" = "True" ]; then
ldsc \
--print-snps $[snp_list] \
$[ld_window_flag] $[ld_window_param] \
--out "$out_prefix" \
--bfile $[_input[-1]:nar] \
--yes-really \
--annot "$annot" \
--l2
else
compute_ldscores \
--annot "$annot" \
--bfile $[_input[-1]:nar] \
$[ld_window_flag] $[ld_window_param] \
--out "${out_prefix}.$[ldscore_ext]" \
--allow-missing
fi
}
if [ "$[str(emit_single)]" = "True" ]; then
for i in $(seq 1 $[N_targets]); do
name="$[annotation_name]_single_$i"
annot="$[cwd:a]/$name/$name.$chrom.annot.gz"
prefix="$[cwd:a]/$name/$name.$chrom"
run_polyfun "$annot" "$prefix"
done
fi
if [ "$[str(emit_joint)]" = "True" ]; then
name="$[annotation_name]_joint"
annot="$[cwd:a]/$name/$name.$chrom.annot.gz"
prefix="$[cwd:a]/$name/$name.$chrom"
run_polyfun "$annot" "$prefix"
fi
# ---- Step D: write .l2.M (+ .l2.M_5_50 when .frq present) files (R worker) ----
bash: expand = "${ }", stderr = f'{_output[2]}.stderr', stdout = f'{_output[2]}.stdout'
Rscript ${modular_script_dir}/enrichment/make_annotation.R --step mfiles \
--annotation-name ${annotation_name} --cwd ${cwd:a} --chrom ${input_chroms[_index]} \
${('--emit-single' if emit_single else '')} ${('--emit-joint' if emit_joint else '')} \
--n-targets ${N_targets} --ldscore-ext ${ldscore_ext} \
${('--use-print-snps' if use_print_snps else '')} \
--frq-file ${frqfile_dir}/${plink_name}${input_chroms[_index]}.frq
Calculate Functional Enrichment using Annotations#
[get_heritability]
# Per-trait sLDSC regression via polyfun. Fans out across target_anno_dirs:
# each (trait, target_dir) pair becomes one polyfun invocation. Outputs go to
# <cwd>/<basename(target_dir)>/<trait>.{results,log,part_delete}.
#
# `target_anno_dirs` is the list produced by [make_annotation_files_ldscore]:
# typically N _single_<i> directories plus optionally one _joint directory.
#
# --- about the ".results" Category column and the "_0 / _1" suffix ---
# Each (trait, target_dir) pair is ONE polyfun call; its `ldsc --ref-ld-chr`
# always gets exactly two LD-score sources, in this order:
# "<target_dir>/<target>." (index 0) , "<baseline_dir>/<baseline>" (index 1)
# With --overlap-annot, every annotation column in the .results "Category" is
# named <ldscore_column_name>_<ref-ld-index>:
# index 0 = the target file -> "ANNOT_0" (no-snplist; compute_ldscores keeps the annot col name)
# -> "L2_0" (snplist + single annot; ldsc hard-codes "L2", see below)
# -> "ANNOT_1_0","ANNOT_2_0" (no-snplist joint dir, N>=2 annot cols)
# -> "ANNOT_1L2_0","ANNOT_2L2_0" (snplist joint dir, N>=2 -> "<name>L2")
# index 1 = the baseline file -> "base_1","Coding_UCSC_1", ... (the 97 baseline annots)
# So in this pipeline the suffix is only ever 0 (target) or 1 (baseline); it would
# continue 0,1,2,... only if you handed `ldsc --ref-ld-chr` more than two sources.
# (Why ANNOT_0 vs L2_0: see the [make_annotation_files_ldscore] header — ldsc's
# "n_annot == 1 -> column name 'L2'" quirk vs compute_ldscores keeping the annot
# column name.) [postprocess] auto-detects the target Category; if you instead pass
# --target-categories, the names must match this column exactly.
#
parameter: target_anno_dirs = paths()
parameter: all_traits = []
import os
with open(all_traits_file, 'r') as f:
trait_paths = [os.path.join(sumstat_dir, line.strip()) for line in f if line.strip()]
# Build (trait, target_dir) Cartesian product as parallel flat lists.
input_list = []
target_meta = []
for td in target_anno_dirs:
for t in trait_paths:
input_list.append(t)
target_meta.append(str(td))
input: input_list, group_by = 1, group_with = "target_meta"
output: f"{cwd:a}/{os.path.basename(target_meta[_index])}/{os.path.basename(_input[0])}.log", \
f"{cwd:a}/{os.path.basename(target_meta[_index])}/{os.path.basename(_input[0])}.results"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand = "${ }"
target_dir="${target_meta[_index]}"
target_name="$(basename ${target_meta[_index]})"
trait="$(basename ${_input[0]})"
output_dir="${cwd:a}/$target_name"
mkdir -p "$output_dir"
# MAF cutoff handling. Only 0 (disabled) or 0.05 (sLDSC default) are supported;
# other values would require recomputing LD scores at that cutoff.
frq_file_check="${frqfile_dir}/${plink_name}22.frq"
if [ "${maf_cutoff}" = "0" ] || [ "${maf_cutoff}" = "0.0" ]; then
echo "maf_cutoff = 0: skipping MAF filter (--not-M-5-50)"
frq_option="--not-M-5-50"
elif [ "${maf_cutoff}" = "0.05" ]; then
if [ -f "$frq_file_check" ]; then
echo "maf_cutoff = 0.05: using --frqfile-chr (MAF > 5%)"
frq_option="--frqfile-chr ${frqfile_dir}/${plink_name}"
else
echo "ERROR: maf_cutoff=0.05 requires .frq files for the reference panel,"
echo " but none found at ${frqfile_dir}/${plink_name}*.frq."
echo " Provide .frq files in frqfile_dir, or set maf_cutoff=0 (NOT recommended)."
exit 1
fi
else
echo "ERROR: maf_cutoff=${maf_cutoff} is not supported. Only 0 (no filter) or"
echo " 0.05 (sLDSC default) are accepted. Other values would require"
echo " recomputing LD scores at that cutoff."
exit 1
fi
run_ldsc() {
local extra_args="$1"
ldsc \
--h2 ${sumstat_dir}/$trait \
--ref-ld-chr "$target_dir/$target_name.","${baseline_ld_dir}/${baseline_name}" \
--out "$output_dir/$trait" \
--overlap-annot \
--w-ld-chr ${weights_dir}/${weight_name} \
$frq_option \
--print-coefficients \
--print-delete-vals \
--n-blocks ${n_blocks} \
$extra_args
}
run_ldsc ""
log_file="$output_dir/$trait.log"
# FloatingPointError retry ladder (preserved from original): 30 -> 20 -> 10
for max in 30 20 10; do
if [ -f "$log_file" ] && grep -q "FloatingPointError\|invalid value encountered in sqrt" "$log_file"; then
echo "FloatingPointError detected, retrying with --chisq-max $max..."
run_ldsc "--chisq-max $max"
else
break
fi
done
if [ -f "$log_file" ] && grep -q "FloatingPointError\|invalid value encountered in sqrt" "$log_file"; then
echo "ERROR: FloatingPointError persists for trait $trait at target $target_name even with --chisq-max 10"
echo "This trait may have severe numerical instability issues in the summary statistics."
fi
[munge_sumstats_polyfun]
parameter: sumstats = path
parameter: n = 0
parameter: min_info = 0.6
parameter: min_maf = 0.001
parameter: keep_hla = False
parameter: chi2_cut = 30
input: sumstats
output: f"{_input:n}.munged.parquet"
bash: expand=True, stderr=f'{_output:nn}.stderr', stdout=f'{_output:nn}.stdout'
munge_polyfun_sumstats \
--sumstats {_input} \
--out {_output} \
{'--n {}'.format(n) if n>0 else ''} \
{'--min-info {}'.format(min_info)} \
{'--min-maf {}'.format(min_maf)} \
{'--chi2-cutoff {}'.format(chi2_cut)} \
{'--keep-hla' if keep_hla else ''} \
--remove-strand-ambig
[postprocess]
# Post-processing of polyfun outputs via pecotmr::sldscPostprocessingPipeline,
# driven by the sldsc_postprocess.R wrapper. The wrapper reads .results / .log /
# .part_delete for all traits in `traits_file` (both single-target and, when
# present, joint-target runs), the target .annot.gz files, and the reference
# .frq files, bundles them into a pecotmr SldscData object, computes Gazal-style
# tau*, EnrichStat with back-solved jackknife SE, and runs the default
# DerSimonian-Laird random-effects meta across all supplied traits. Writes one
# RDS with per-trait tables and the meta tables (tauStar, enrichment, enrichstat).
parameter: traits_file = path() # text file: one trait sumstats filename per line
parameter: heritability_cwd = path() # parent directory of [get_heritability] outputs (contains <annotation_name>_single_<i>/ subdirs and optionally <annotation_name>_joint/)
parameter: target_categories = [] # target annotation names. Auto-detected from the joint-run results if empty.
parameter: target_categories_label = [] # optional display names, same order as target_categories;
# when given, every "target" column / tau*-block colname in
# the output RDS is renamed to these (params$target_categories
# holds the labels, params$target_categories_orig the originals).
parameter: target_anno_dir = path() # directory of target .annot.gz files used for sd_C and binary detection (typically the joint dir, since it carries all target columns)
input: traits_file
output: f"{cwd:a}/{annotation_name}.sldsc_postprocess.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
Rscript ${modular_script_dir}/pecotmr_integration/sldsc_postprocess.R \
--traits-file ${traits_file} \
--heritability-cwd ${heritability_cwd} \
--annotation-name ${annotation_name} \
--target-anno-dir ${target_anno_dir} \
--frqfile-dir ${frqfile_dir} \
--plink-name ${plink_name} \
--maf-cutoff ${maf_cutoff} \
${('--target-categories ' + ','.join(target_categories)) if target_categories else ''} \
${('--target-categories-label ' + ','.join(target_categories_label)) if target_categories_label else ''} \
--output ${_output}
[meta_subset]
# Optional: re-run random-effects meta on a user-defined subset of traits, using
# the cached per-trait standardized results from [postprocess] via the
# sldsc_meta_subset.R wrapper. No regression rerun.
parameter: postprocess_rds = path() # output of [postprocess]
parameter: subset_traits_file = path() # text file: one trait id per line, subset of those passed to [postprocess]
parameter: subset_name = str # label used in the output filename
parameter: target_categories = [] # target annotation names to meta on; if empty, uses all from postprocess output
# If [postprocess] was run with --target-categories-label, the cached RDS already
# carries the display names (params$target_categories = the labels), so leave
# --target-categories empty here (or pass the labels, not the original ANNOT_* names).
input: postprocess_rds, subset_traits_file
output: f"{cwd:a}/{annotation_name}.{subset_name}.meta.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
Rscript ${modular_script_dir}/pecotmr_integration/sldsc_meta_subset.R \
--postprocess-rds ${postprocess_rds} \
--subset-traits-file ${subset_traits_file} \
${('--target-categories ' + ','.join(target_categories)) if target_categories else ''} \
--output ${_output}