SNP heritability and functional enrichment

Overview

estimateH2() estimates SNP heritability from GWAS summary statistics and an LD reference. Four estimators are available, differing in what they model and what they need:

`method`	Approach	LD input
`"sldsc"`	Stratified LD score regression (Finucane 2015)	`LdScore`
`"gldsc"`	Generalized LDSC — GLS with full block LD (Xiong 2024)	`LdScore` + matrices
`"lder"`	LD eigenvalue regression (Song 2022)	`LdEigen`
`"hdl"`	High-definition likelihood (Ning 2020)	`LdEigen`

Supplying annotations makes any of them stratified, adding a per-annotation enrichment table.

There is a second, separate track in the package: when stratified LDSC has already been run externally by polyfun, sldscPostprocessingPipeline() postprocesses those outputs. Both are covered below.

library(pecotmr)
data(ldScoreExample, ldEigenExample)
ldScoreExample

## LdScore: 1300 SNPs, 1 LD score columns
##   Block LD: stored matrices
##   Reference N: 20000, In-sample: FALSE

Both bundled references describe the same 1300 variants across 14 chr22 LD blocks, so the four methods can be compared on one dataset. Their LD is that of a reference panel of 20000 people simulated from the chr22 LD of the de-identified protocol_example panel. That panel itself is a sketch of 60 people: far too few for these estimators, which correct for the noise in a reference’s LD, so the examples carry the LD of a large panel drawn from it instead (inst/scripts/build_h2_examples.R).

Variants must line up

Every estimator reads the z-scores positionally, indexed by the LD reference’s own variant order — nothing is matched on variant ID along the way. estimateH2() therefore checks up front that the summary statistics, the annotations and the reference describe the same variants, in the same order, on the same genome build, and refuses to run otherwise.

The practical consequence: build the summary statistics from the reference rather than alongside it.

# A GWAS simulated on this reference so the true h2 is known: z is drawn
# from the standard model z ~ N(0, N R diag(v) R + R), block by block.
M <- length(ldScoreExample)
nGwas <- 10000
h2True <- 0.4
perSnpVar <- rep(h2True / M, M)

set.seed(1)
z <- numeric(M)
for (block in ldMatrixList(ldScoreExample)) {
    idx <- block$snpIdx
    p <- length(idx)
    sigma <- nGwas * (block$R %*% diag(perSnpVar[idx], p) %*% block$R) +
        block$R
    e <- eigen((sigma + t(sigma)) / 2, symmetric = TRUE)
    z[idx] <- as.vector(e$vectors %*% (sqrt(pmax(e$values, 0)) * rnorm(p)))
}

gr <- as(ldScoreExample, "GRanges")
S4Vectors::mcols(gr) <- S4Vectors::DataFrame(
    SNP = names(ldScoreExample),
    A1 = S4Vectors::mcols(ldScoreExample)$A1,
    A2 = S4Vectors::mcols(ldScoreExample)$A2,
    Z = z,
    N = as.integer(nGwas)
)
gwas <- GwasSumStats(
    studyName = "trait1",
    entry = list(gr),
    genome = unique(unname(GenomeInfoDb::genome(ldScoreExample))),
    ldReference = readGenotypes(
        system.file("extdata", "toy_ref.bed", package = "pecotmr")
    )
)

Estimating heritability

methods <- c("sldsc", "gldsc", "lder", "hdl")
estimates <- lapply(methods, function(m) {
    ref <- if (m %in% c("sldsc", "gldsc")) ldScoreExample else ldEigenExample
    suppressWarnings(estimateH2(gwas, ref, method = m))
})
names(estimates) <- methods

data.frame(
    method = methods,
    h2 = round(vapply(estimates, heritabilityEstimate, numeric(1)), 3),
    se = round(vapply(estimates, heritabilityEstimateSe, numeric(1)), 3),
    row.names = NULL
)

##   method    h2    se
## 1  sldsc 0.235 0.051
## 2  gldsc 0.364 0.084
## 3   lder 0.397 0.041
## 4    hdl 0.416 0.046

The true value is 0.4. Three of the four land close to it; sldsc comes in low, and that is worth understanding rather than tuning away.

sldsc is attenuated because the reference is small. LD score regression is a genome-scale method — real analyses regress over roughly a million variants. On 1300 the intercept is not identifiable: an unweighted OLS on this same data returns an intercept several times its true value of 1, and the heritability absorbs the difference. lder is the least sensitive to reference size, which is why it is the default.

All four are checked against their reference implementations — ldsc, HDL, LDER and gldsc — on identical inputs: HDL agrees exactly, LDER to machine precision, gLDSC to ~3e-8 and S-LDSC to ~1e-6. The small residuals are a ridge on the weighted solve, which keeps an exactly singular design from returning NA. Agreement holds even where the reference implementations misbehave, such as an S-LDSC intercept below 1 or a negative partitioned heritability.

One deliberate difference: LDER here clamps negative eigenvalues before forming its regression weights. A real, rank-deficient LD block decomposes with a tail of tiny negative values that are floating-point noise, and the original LDER implementation takes sqrt() of the resulting negative weight and aborts. The clamp changes nothing wherever the original runs.

hdl has a reference-size limit, and it tells you when it hits it. It is the only one of the four that models reference-panel noise explicitly: its likelihood subtracts a term scaling as 1/nRef from the modelled variance. Once the reference is small relative to the GWAS, that term outgrows the signal, the modelled variance is driven to its floor across most eigen-directions, and the optimiser walks the intercept out to the top of its range — returning a number that reflects that bound rather than the data. This reference’s 20000 people are plenty: even at N = 1e5 HDL recovers about 0.40 for a true 0.4. Against a 1,000-person panel at that N it drifts to about 0.57 and warns that the fit sits on a bound. If HDL disagrees sharply with the others, check for that warning and compare nSamples() against ldPanelNSamples().

A standard error of zero, or an estimate pinned at a boundary, almost always means the reference is too small rather than that the trait has no heritability. The block jackknife needs at least two LD blocks and wants many more; estimateH2() refuses a single-block reference outright.

Building an LD reference

The bundled references were built with buildLdEigen() / buildLdScore(), which turn an LD reference into the statistic each estimator needs. The reference comes per LD block: a genotype panel together with the blocks to cut it into, or one LdMatrix per block. One block has no jackknife.

panel <- readGenotypes(plink2Prefix = "/path/to/reference.chr22")
blocks <- read.delim("/path/to/ld_blocks.tsv") # chrom, start, end

ldEigen <- buildLdEigen(panel, ldBlocks = blocks, genome = "hg38")
ldScore <- buildLdScore(panel, ldBlocks = blocks, genome = "hg38")

# Or from LD matrices you built, one LdMatrix per block:
ldScore <- buildLdScore(list(ldBlock1, ldBlock2, ldBlock3), genome = "hg38")

nSamples (the reference panel size) and the genome build default to what the reference records: a panel’s sample count (for a stochastic sketch, the effective size of its LD), or the nSamples of each LdMatrix. The panel size matters beyond HDL: LD scores sum ldsc’s unbiased r², r^2 - (1 - r^2) / (n - 2), so n should be the number of samples the LD was estimated from. buildLdScore(keepLdMatrices = FALSE) drops the per-block matrices, which shrinks the object by roughly an order of magnitude — but only unstratified "sldsc" can use the result, since "gldsc" and stratified "sldsc" both form their per-annotation scores from each block’s LD. An LdScore or LdEigen built from a genotype panel stores none of it: it keeps the panel and reads each block’s LD when an estimator reaches that block, so every method runs on it while only one block’s LD is in memory at a time.

Stratified estimation and enrichment

Passing an AnnotationMatrix partitions heritability across annotations. The annotation rows must line up with the reference, exactly as the summary statistics do.

# An illustrative annotation over the first half of the reference, simulated
# to carry three times the per-variant heritability of the rest.
inA <- seq_len(M) <= M / 2
perSnpVarStrat <- ifelse(inA, 0.225 / sum(inA), 0.075 / sum(!inA))

set.seed(1)
zStrat <- numeric(M)
for (block in ldMatrixList(ldScoreExample)) {
    idx <- block$snpIdx
    p <- length(idx)
    sigma <- nGwas *
        (block$R %*% diag(perSnpVarStrat[idx], p) %*% block$R) +
        block$R
    e <- eigen((sigma + t(sigma)) / 2, symmetric = TRUE)
    zStrat[idx] <- as.vector(
        e$vectors %*% (sqrt(pmax(e$values, 0)) * rnorm(p))
    )
}
S4Vectors::mcols(gr)$Z <- zStrat
gwasStrat <- GwasSumStats(
    studyName = "trait1",
    entry = list(gr),
    genome = unique(unname(GenomeInfoDb::genome(ldScoreExample))),
    ldReference = readGenotypes(
        system.file("extdata", "toy_ref.bed", package = "pecotmr")
    )
)

annot <- AnnotationMatrix(
    matrix(as.integer(inA), ncol = 1, dimnames = list(NULL, "annotA")),
    as(ldScoreExample, "GRanges"),
    data.frame(name = "annotA", tier = "baseline", type = "binary"),
    genome = unique(unname(GenomeInfoDb::genome(ldScoreExample)))
)

strat <- estimateH2(
    gwasStrat,
    ldScoreExample,
    method = "gldsc",
    annotations = annot
)
as.data.frame(stratifiedHeritabilityEnrichment(strat))[, c(
    "annotation",
    "tau",
    "enrichment",
    "enrichmentSe",
    "propH2",
    "propSnps"
)]

##   annotation          tau enrichment enrichmentSe   propH2 propSnps
## 1     annotA 0.0001337725    1.30079    0.1187998 0.650395      0.5

enrichment is the proportion of heritability divided by the proportion of variants: above 1 means the annotation carries more heritability than its share of the genome. The simulated truth here is 1.5 — three times the per-variant heritability on half the variants puts 75% of it in annotA — and the estimate is directionally right but attenuated, again a consequence of the reference’s size rather than of the annotation.

Annotations come from files via readAnnotations(), which reads BED, BigWig and LDSC .annot formats at a given set of variant positions:

annot <- readAnnotations(
    c(enhancer = "enhancers.bed", conservation = "phastCons.bw"),
    snpRanges = as(ldScoreExample, "GRanges"),
    genome = "hg38"
)

The tier column decides how each annotation is used: "baseline" annotations enter the model and get an enrichment estimate, while "candidate" annotations are scored against the fitted baseline (the sHDL score test, via scoreStats()).

One structural caveat: a set of baseline annotations that exactly partitions the variants is rank-deficient by construction, because the model’s implicit all-ones base annotation is then their sum. Real baseline models — the LDSC baselineLD set, for instance — use overlapping annotations and do not have this problem.

Postprocessing external polyfun runs

The second track handles stratified LDSC that was run outside R by polyfun. The readers load its outputs, SldscData bundles them, and the pipeline does the computation — no file I/O of its own.

sldscDir <- system.file("extdata", "sldsc", package = "pecotmr")
traitDir <- system.file("extdata", "sldsc_trait", package = "pecotmr")

annotTable <- readSldscAnnot(sldscDir)
frq <- readSldscFrq(sldscDir, plinkName = "reference.")

traitNames <- c(traitA = "sumstats", traitB = "sumstats2", traitC = "sumstats3")
runs <- lapply(traitNames, function(t) {
    readSldscTrait(file.path(traitDir, paste0(t, ".parquet")))
})
runs$traitA$categories

## [1] "ANNOT_0"                         "Coding_UCSC_common_1"           
## [3] "Coding_UCSC_lowfreq_1"           "Conserved_LindbladToh_common_1" 
## [5] "Conserved_LindbladToh_lowfreq_1" "Repressed_Hoffman_common_1"     
## [7] "Repressed_Hoffman_lowfreq_1"     "base_1"

Those categories are a real baseline model — Coding_UCSC, Conserved_LindbladToh and Repressed_Hoffman, each split into common and low-frequency, plus base — with ANNOT_0 as the target annotation under test.

sldscData <- SldscData(
    annotData = annotTable,
    frqData = frq,
    traits = lapply(runs, function(r) list(single = list(r)))
)
result <- sldscPostprocessingPipeline(
    sldscData,
    targetCategories = "ANNOT_0"
)
as.data.frame(result$meta$enrichment)

##    target isBinary singleMean singleSe   singleP nTraits
## 1 ANNOT_0     TRUE   20.46742  9.91854 0.3909615       3

The pipeline standardizes each trait’s tau to Gazal tau*, back-solves a jackknife standard error for the enrichment statistic, and pools traits by DerSimonian-Laird random-effects meta-analysis. sldscSubsetMeta() re-runs the meta over a subset of traits without refitting anything:

subset <- sldscSubsetMeta(result, subsetTraits = c("traitB", "traitC"))
unlist(subset$enrichment$ANNOT_0)

##                mean                  se                   p             nTraits 
##  "18.2661713350886"  "11.9393819547916" "0.126039184244747"                 "2" 
##         traitsUsed1         traitsUsed2                tau2 
##            "traitB"            "traitC"                 "0"

Each polyfun run may be a single target run or a joint one; supplying both lets the pipeline report the target’s enrichment from the single fit alongside its coefficient from the joint fit.

The postprocessing steps

sldscPostprocessingPipeline() is a wrapper over steps that are exported individually, which is what you want when auditing a surprising number:

computeSldscMRef(sldscData, mafCutoff = 0.05) # the M_ref that scales tau*

## [1] 319

isBinarySldscAnnot(sldscData) # binary vs continuous

## ANNOT 
##  TRUE

computeSldscAnnotSd() gives the per-annotation standard deviation that turns tau into Gazal tau*, standardizeSldscTrait() does that standardisation for one trait’s run, and metaSldscRandom() is the DerSimonian-Laird pooling across traits. Running them yourself is the way to see which of the three a strange tau* came from.

annotData(), frqData(), annotCols() and sldscResults() read the corresponding parts of an SldscData back out.

LD scores from an existing reference

computeLdScores() derives per-variant LD scores from an LdEigen or LdScore, optionally stratified by annotation:

scores <- computeLdScores(ldEigenExample)
head(round(as.vector(scores[, 1]), 3))

## [1] 3.146 3.436 2.655 3.313 2.464 3.370

This is the same quantity buildLdScore() stores: each variant’s sum, over its block, of ldsc’s unbiased r² rather than the raw sample r², which overstates LD on a finite panel. From an LdEigen the raw sum comes out of the eigendecomposition as sum_i V[j,i]^2 d[i]^2, and the correction applies to that sum exactly, so the two agree — worth knowing when you have one object and need the other’s input.

Bridging the two tracks

h2EstimateToSldscTrait() converts an H2Estimate into the same list shape readSldscTrait() produces, so an in-package estimate can be meta-analysed alongside polyfun runs:

names(h2EstimateToSldscTrait(strat))

##  [1] "categories"               "tau"                     
##  [3] "tauSe"                    "enrichment"              
##  [5] "enrichmentSe"             "enrichmentP"             
##  [7] "propH2"                   "propSnps"                
##  [9] "h2g"                      "annotationJackknifeCoefs"
## [11] "nBlocks"

Next steps

Molecular QTL enrichment — the complementary per-variant QTL/GWAS enrichment analysis.
QC for summary statistics — preparing the GWAS side.

Session info

sessionInfo()

## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.34.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] pecotmr_0.9.3
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.3.0                   bitops_1.1-0               
##   [3] httr2_1.3.0                 rlang_1.3.0                
##   [5] magrittr_2.0.5              otel_0.2.0                 
##   [7] matrixStats_1.5.0           susieR_0.16.6              
##   [9] compiler_4.5.3              RSQLite_3.53.3             
##  [11] systemfonts_1.3.2           vctrs_0.7.3                
##  [13] quadprog_1.5-8              stringr_1.6.0              
##  [15] pkgconfig_2.0.3             crayon_1.5.3               
##  [17] fastmap_1.2.0               dbplyr_2.6.0               
##  [19] backports_1.5.1             XVector_0.50.0             
##  [21] Rsamtools_2.26.0            rmarkdown_2.32             
##  [23] tzdb_0.5.0                  UCSC.utils_1.6.1           
##  [25] ragg_1.5.2                  purrr_1.2.2                
##  [27] bit_4.6.0                   xfun_0.61                  
##  [29] Rfast_2.1.5.2               MultiAssayExperiment_1.36.2
##  [31] cachem_1.1.0                GenomeInfoDb_1.46.2        
##  [33] jsonlite_2.0.0              blob_1.3.0                 
##  [35] DelayedArray_0.36.1         reshape_0.8.10             
##  [37] tictoc_1.2.1                BiocParallel_1.44.0        
##  [39] irlba_2.4.1                 parallel_4.5.3             
##  [41] R6_2.6.1                    bslib_0.12.0               
##  [43] stringi_1.8.9               RColorBrewer_1.1-3         
##  [45] GenomicRanges_1.62.1        jquerylib_0.1.4            
##  [47] numDeriv_2016.8-1.1         Rcpp_1.1.2                 
##  [49] Seqinfo_1.0.0               SummarizedExperiment_1.40.0
##  [51] knitr_1.52                  readr_2.2.0                
##  [53] IRanges_2.44.0              BiocBaseUtils_1.12.0       
##  [55] Matrix_1.7-6                tidyselect_1.2.1           
##  [57] abind_1.4-8                 yaml_2.3.12                
##  [59] codetools_0.2-20            curl_8.0.0                 
##  [61] metafor_5.0-1               lattice_0.23-1             
##  [63] tibble_3.3.1                plyr_1.8.9                 
##  [65] withr_3.0.3                 Biobase_2.70.0             
##  [67] S7_0.2.2                    evaluate_1.0.5             
##  [69] archive_1.1.13              desc_1.4.3                 
##  [71] RcppParallel_5.1.11-2       BiocFileCache_3.0.0        
##  [73] Biostrings_2.78.0           filelock_1.0.3             
##  [75] pillar_1.11.1               MatrixGenerics_1.22.0      
##  [77] metadat_1.6-0               checkmate_2.3.3            
##  [79] stats4_4.5.3                generics_0.1.4             
##  [81] vroom_1.7.1                 mathjaxr_2.0-0             
##  [83] S4Vectors_0.48.1            hms_1.1.4                  
##  [85] ggplot2_4.0.3               scales_1.4.0               
##  [87] glue_1.8.1                  tools_4.5.3                
##  [89] fs_2.1.0                    grid_4.5.3                 
##  [91] PlinkMatrix_1.0.0           tidyr_1.3.2                
##  [93] nlme_3.1-170                cli_3.6.6                  
##  [95] zigg_0.0.2                  textshaping_1.0.5          
##  [97] mixsqp_0.3-54               S4Arrays_1.10.1            
##  [99] dplyr_1.2.1                 gtable_0.3.6               
## [101] sass_0.4.10                 digest_0.6.39              
## [103] BiocGenerics_0.56.0         SparseArray_1.10.10        
## [105] htmlwidgets_1.6.4           farver_2.1.2               
## [107] memoise_2.0.1               htmltools_0.5.9            
## [109] pkgdown_2.2.1               lifecycle_1.0.5            
## [111] httr_1.4.9                  bit64_4.8.6                
## [113] MASS_7.3-66

pecotmr authors

2026-10-22