Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

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, EE, hg2h^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))1/(p(1-p)) where pp is the minor allele frequency (MAF), the expected χ2\chi^2 statistic of variant jj is:

E[χj2∣ℓj]  =  N h2 ℓjM  +  Na  +  1(1)E[\chi^2_j \mid \ell_j] \;=\; \frac{N\,h^2\,\ell_j}{M} \;+\; N a \;+\; 1 \quad (1)
(1)

where NN is the sample size; MM is the number of SNPs, so that h2/Mh^2/M is the average heritability per SNP; aa measures the contribution of confounding biases such as cryptic relatedness and population stratification; and ℓj=∑krjk2\ell_j = \sum_k r^2_{jk} is the LD Score of variant jj, which measures the amount of genetic variation tagged by jj. 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 ajCa_{jC} denote the value of annotation CC at SNP jj:

  • Binary annotation (e.g. an indicator for “in enhancer”, “in exon”, “in cell-type-specific peak”): ajC∈{0,1}a_{jC} \in \{0, 1\}.

  • Continuous annotation (e.g. gene-specificity score, conservation score, continuous epigenomic signal): ajC∈Ra_{jC} \in \mathbb{R}.

Under a polygenic model the per-SNP heritability for SNP jj is

Var(βj)  =  ∑CajC τC\mathrm{Var}(\beta_j) \;=\; \sum_C a_{jC}\, \tau_C
(2)

and the expected χ2\chi^2 statistic of SNP jj is

E[χj2∣aj]  =  N∑CτC ℓ(j,C)  +  Na  +  1(2)E[\chi^2_j \mid \mathbf{a}_j] \;=\; N \sum_C \tau_C\, \ell(j, C) \;+\; N a \;+\; 1 \quad (2)
(3)

where ℓ(j,C)=∑kakC rjk2\ell(j, C) = \sum_k a_{kC}\, r^2_{jk} is the partitioned LD score of SNP jj with respect to annotation CC, and aa measures confounding bias. Equation (2) allows joint estimation of all τC\tau_C via a (computationally simple) multiple regression of χj2\chi^2_j against ℓ(j,C)\ell(j, C).

Interpretation of τC\tau_C:

  • Binary CC: τC\tau_C is the additive increase in per-SNP heritability for SNPs in category CC, on top of the contributions from any other annotations they belong to.

  • Continuous CC: τC\tau_C is the additive change in per-SNP heritability per unit increase in the value of annotation CC.

Tau estimation and enrichment analysis

Notation

For each annotation CC:

  • πCh2\pi^{h^2}_C = proportion of trait heritability hg2h^2_g assigned to annotation CC.

  • πCM\pi^{M}_C = proportion of (effective) SNPs in annotation CC. For binary annotations this is MC/MrefM_C / M_{\mathrm{ref}}; for continuous annotations it is the share of total annotation weight in CC.

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:

  • τC\tau_C and its standard error — (polyfun).

  • πCh2\pi^{h^2}_C and πCM\pi^{M}_C — (polyfun).

  • EC=πCh2/πCME_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 hg2h^2_g — (polyfun).

  • The 200-block jackknife delete-values of τC\tau_C — (polyfun).

Quantities from the post-processing layer (pecotmr)

From the polyfun outputs above plus the reference panel, the post-processing layer computes:

  • sdCsd_C — per-annotation standard deviation over MAF >> cutoff SNPs — (pecotmr: compute_sldsc_annot_sd).

  • MrefM_{\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).

  • τC∗\tau^*_C point estimate and per-block τC∗\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 τC∗\tau^*_C, ECE_C, or EnrichStat across traits — (pecotmr: metaSldscRandom).

The top-level entry point pecotmr::sldscPostprocessingPipeline orchestrates all of the above.

Standardized tau (τ∗\tau^*) — (pecotmr)

τC\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)

τC∗  =  τC⋅sdC⋅Mrefhg2\tau^*_C \;=\; \tau_C \cdot \frac{sd_C \cdot M_{\mathrm{ref}}}{h^2_g}
(4)

interpreted as the additive change in per-SNP heritability associated with a 1 standard deviation increase in annotation CC, divided by the average per-SNP heritability across all SNPs. τC∗\tau^*_C is dimensionless and comparable across annotations and across traits. In a joint multi-annotation regression it is the independent contribution of annotation CC after controlling for overlapping effects of the others.

Here sdCsd_C is the standard deviation of annotation CC across reference SNPs (MAF >> cutoff), MrefM_{\mathrm{ref}} is the count of those SNPs, and hg2h^2_g is the trait heritability. Applying the same scaling to each of the 200 jackknife blocks yields per-block τC∗\tau^*_C values; their sample variance gives the jackknife standard error

SEjackknife(τC∗)  =   (B−1)2B Varb(τC,(b)∗) SE^{\text{jackknife}}(\tau^*_C) \;=\; \sqrt{\,\tfrac{(B-1)^2}{B}\, \mathrm{Var}_b(\tau^*_{C,(b)})\,}
(5)

with B=200B = 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 CC differs from outside it (Finucane et al. 2015):

EnrichStatC  =  hg2Mref ⁣[ πCh2πCM  −  1−πCh21−πCM ]\text{EnrichStat}_C \;=\; \frac{h^2_g}{M_{\mathrm{ref}}} \!\left[\, \frac{\pi^{h^2}_C}{\pi^{M}_C} \;-\; \frac{1 - \pi^{h^2}_C}{1 - \pi^{M}_C} \,\right]
(6)

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:

∣ZC∣  =  Φ−1 ⁣(1−pC2),SEEnrichStatC  =  ∣EnrichStatC∣∣ZC∣.|Z_C| \;=\; \Phi^{-1}\!\left(1 - \tfrac{p_C}{2}\right), \qquad SE_{\text{EnrichStat}_C} \;=\; \frac{|\text{EnrichStat}_C|}{|Z_C|}.
(7)

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), πCM=MC/Mref\pi^{M}_C = M_C / M_{\mathrm{ref}} has a direct interpretation and ECE_C reads as “the category explains ECE_C-fold more heritability than its share of SNPs.” The within-type headline quantities are therefore ECE_C and the EnrichStat p-value; τC∗\tau^*_C is reported alongside.

For a continuous annotation (e.g. gene-specificity score, conservation score, continuous epigenomic signal), ECE_C depends on the scale of the annotation: rescaling the annotation by a constant changes ECE_C even though the underlying biology is unchanged. The within-type headline quantities are therefore τC∗\tau^*_C and its p-value; ECE_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}\{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-coefficients outputs 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-coefficients to polyfun for this reason.

Cross-type comparison: always use τC∗\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 τC∗\tau^*_C. The standardization in Gazal et al. (2017) was designed for exactly this purpose: sdC=p(1−p)sd_C = \sqrt{p(1-p)} for a binary annotation (where pp is the proportion in the category) and sdC=sd_C = empirical standard deviation for a continuous annotation, so the resulting τC∗\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. ECE_C does not have this property and must not be compared across types.

The pipeline emits both ECE_C and τC∗\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 EE are conditional on the other annotations in the model; joint τC∗\tau^*_C is reported as the independent contribution of annotation CC 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 NN single-target outputs, the joint output, or both in one invocation. With both defaults the post-processing layer reads all N+1N+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")):

θ^meta  =  ∑iwi θ^i∑iwi,SEmeta  =  1∑iwi,wi  =  1SEi2+σ^2\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}
(8)

where θ^i\hat\theta_i is the per-trait estimate and SEiSE_i its standard error:

  • For τC∗\tau^*_C meta: SEiSE_i is the jackknife SE from the per-block τC∗\tau^*_C values.

  • For ECE_C meta: SEiSE_i is the polyfun-reported Enrichment_std_error.

  • For EnrichStat meta: SEiSE_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 ECE_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.

Zmeta  =  θ^metaSEmeta,p  =  2 Φ(−∣Zmeta∣)Z_{\mathrm{meta}} \;=\; \frac{\hat\theta_{\mathrm{meta}}}{SE_{\mathrm{meta}}}, \qquad p \;=\; 2\,\Phi(-|Z_{\mathrm{meta}}|)
(9)

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 (default output): work directory for the analysis.

  • --modular-script-dir (default code/script): directory holding the modular R wrapper scripts.

  • --maf-cutoff (default 0.05): MAF cutoff, per the S-LDSC recommendation, since rare-variant LD is unstable and HapMap3-style regression weights are common-variant by construction. 0 opts 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 (columns path, 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.gz form: 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-snps writes .l2.ldscore.gz; without one, compute_ldscores writes .l2.ldscore.parquet. It also changes the LD-score column name that postprocess matches on, so --target-categories may need to be passed explicitly.

  • --ld-wind-kb (default 0) / --ld-wind-cm (default 1.0): LD window. Using the cM window requires genetic-map positions in the .bim fourth column; switch to --ld-wind-kb if it holds 0.

  • --compute-single / --no-compute-single and --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 (default 3) 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 by make_annotation_files_ldscore. Each (trait, target directory) pair is one polyfun call.

  • --baseline-ld-dir and --baseline-name (default baseline_chr): baseline LD scores and their per-chromosome file prefix.

  • --frqfile-dir and --plink-name (default ADSP_chr): PLINK .frq allele frequencies, used to enforce the MAF cutoff.

  • --weights-dir and --weight-name (default weights_chr): LD-score regression weights.

  • --n-blocks (default 200): 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 (default 0.6), --min-maf (default 0.001), --chi2-cut (default 30): 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: the get_heritability output directory to read.

  • --target-categories and --target-categories-label: which .results Category rows are the target annotation, and how to label them. Auto-detected unless --snp-list changed 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-file and --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 -- from make_annotation_files_ldscore: polyfun .annot.gz files plus per-annotation LD scores (.l2.ldscore.parquet or .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> -- from get_heritability: the S-LDSC regression output, one set per trait per target directory. The .results Category column carries the annotation name with a _<ref-ld-index> suffix.

  • <annotation_name>.sldsc_postprocess.rds -- from postprocess: 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 -- from meta_subset: a re-meta of the cached .sldsc_postprocess.rds over 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, tests/fixtures/sldsc_enrichment/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 output/sldsc_ldscore/colocboost_test_annotation_path.txt \
  --reference_anno_file output/sldsc_ldscore/reference_annotation0.txt \
  --genome_ref_file output/sldsc_ldscore/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 tests/fixtures/sldsc_enrichment/sumstats_test_all.txt \
  --sumstat_dir output/sldsc_ldscore/sldsc \
  --baseline_ld_dir output/sldsc_ldscore/sldsc \
  --weights_dir output/sldsc_ldscore/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 tests/fixtures/sldsc_enrichment/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 tests/fixtures/sldsc_enrichment/expected/sldsc_postprocess.rds \
  --subset_traits_file tests/fixtures/sldsc_enrichment/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}
References
  1. Finucane, H. K., Bulik-Sullivan, B., Gusev, A., Trynka, G., Reshef, Y., Loh, P.-R., Anttila, V., Xu, H., Zang, C., Farh, K., Ripke, S., Day, F. R., Purcell, S., Stahl, E., Lindstrom, S., Perry, J. R. B., Okada, Y., Raychaudhuri, S., Daly, M. J., … Price, A. L. (2015). Partitioning heritability by functional annotation using genome-wide association summary statistics. Nature Genetics, 47(11), 1228–1235. 10.1038/ng.3404
  2. Gazal, S., Finucane, H. K., Furlotte, N. A., Loh, P.-R., Palamara, P. F., Liu, X., Schoech, A., Bulik-Sullivan, B., Neale, B. M., Gusev, A., & Price, A. L. (2017). Linkage disequilibrium–dependent architecture of human complex traits shows action of negative selection. Nature Genetics, 49(10), 1421–1427. 10.1038/ng.3954
Enrichment and Validation
GREGOR enrichment analysis
xQTL Modifier Score
scEEMS Model Training