Genetic correlation from summary statistics

Genetic-correlation analysis uses the gcorr R package.

Input boundary

Genome-wide score workflows also accept gact’s list of named z and n matrices, or direct ldsc(z=..., n=..., ldscores=...) input. Matrix summaries must already share allele orientation. Marker names are matched once to their common LD-score reference in reference order; sample sizes are unchanged. No allele metadata or BED loading is needed for this route.

Checked and regional inputs use an ordinary data frame with the qgg column names rsids, chr, pos, ea, nea, stat, and n. Here stat is the signed Z statistic. Use gbase::checkStat() when a Glist reference is available and alignment or QC is required. Marker order, coordinates and allele orientation must already agree across traits. The genome-wide fitters never intersect, flip or rescale traits. The BED HESS preparation described below explicitly matches regional markers and swaps alleles to the reference; it does not infer strand complements. The high-level ldsc() workflow obtains LD scores, annotation scores and reference metadata from a Glist; the focused interfaces also accept those quantities explicitly. A single data frame uses the label "trait"; pass a uniquely named list for multiple traits.

library(gbase)
library(gcorr)

markers <- paste0("rs", 1:20)
one <- data.frame(rsids = markers, chr = "1", pos = 1:20,
  ea = "A", nea = "G", stat = rep(2, 20), n = 10000)

Genome-wide methods

With a Glist containing stored LD scores and the output of gact::getMarkerStat(), the ordinary workflows are:

fit <- gcorr(Glist, stat, method = "ldsc")
fit <- gcorr(Glist, stat, method = "gnova")
fit <- gcorr(Glist, stat, method = "sumher")
fit$estimates

# Stored scores from gact can also be supplied directly.
fit <- ldsc(z = stat$z, n = stat$n, ldscores = ldscores, what = "rg")

The Glist workflow defaults to free LDSC intercepts and an inclusive squared-Z cutoff of 80. GNOVA estimates its overlap intercept using LDSC unless supplied through overlap_intercepts; that preliminary fit adds to workflow time. SumHer constructs one equal-marker baseline with multiplicative inflation and Glist’s reference sample size. This is an illustrative architecture; use gcorr_sumher() for other components. All three create delete blocks when omitted. Matrix inputs use 200 balanced contiguous blocks in LD-score order; checked frames retain chromosome-aware blocks. GNOVA requires constant N per trait. Direct ldsc() retains its previous filtering and intercept defaults.

The recommended LDSC interface uses the same stat contract as glma::vegas():

ordinary <- ldsc(stat = one, Glist = Glist, what = "h2")

multiple <- ldsc(
  stat = list(height = height_stat, bmi = bmi_stat),
  Glist = Glist,
  what = "rg"
)

partitioned <- ldsc(stat = one, Glist = Glist,
  sets = list(coding = coding_markers, regulatory = regulatory_markers),
  residual = TRUE, what = "h2")

grouped <- ldsc(stat = one, Glist = Glist,
  sets = list(region1 = region1_markers, region2 = region2_markers),
  set_model = "grouped", residual = TRUE, what = "h2")

Marker sets are named marker-ID vectors and may overlap. The default set_model = "stratified" uses standard annotation-specific LD scores and accounts for LD crossing set boundaries. set_model = "grouped" instead fits the qgg-style grouped construction by masking the stored ordinary LD-score vector with focal-marker membership; it avoids a sparse-LD scan and is most appropriate when cross-boundary LD is limited. Glist supplies the ordinary and annotation-specific LD scores and reference sample size. Deterministic contiguous delete blocks are generated when blocks is omitted. The result records the selected set model, retains the full fit, and adds convenient h2 and rg components.

Set names are matched once to integer marker positions and reused. Disjoint grouped sets keep one integer membership and one LD score per marker throughout fitting, rather than a dense matrix with one column per set. The explicit partitioned interface also accepts gcorr_set_scores(marker_ids, ldscores, sets) for disjoint sets covering all markers. These indices follow the supplied statistic/reference order; physical BED indices remain a separate Glist mapping. Overlapping and stratified annotation scores retain the existing matrix input.

LDSC, partitioned LDSC, GNOVA and SumHer check marker alignment in R, then pass numeric statistics and LD scores with integer jackknife block membership to native code. Constant sample sizes use one value per trait. Marker names stay in R; filtering returns retained row indices. Nested results retain estimates, uncertainty and fit summaries without per-marker IDs, weights or fitted means.

The explicit focused functions are:

• gcorr_ldsc() and gcorr_ldsc_partitioned();
• gcorr_gnova() for ordinary or binary-annotation GNOVA;
• gcorr_sumher() for architecture components and optional categories.

The package-owned gcorr(..., method=) dispatcher retains the focused interfaces alongside the Glist workflow.

scores <- setNames(seq(1, 2, length.out = 20), markers)
blocks <- setNames(rep(paste0("b", 1:5), each = 4), markers)

explicit <- gcorr_ldsc(one, scores, blocks, marker_count = 50000)
gnova <- gcorr(one, blocks = blocks, ld_scores = scores, method = "gnova")

Partitioned LDSC requires annotation LD scores, weighting LD scores and annotation sums. Annotated GNOVA additionally requires the matching binary annotation matrix and reference-panel sample size. SumHer requires uncorrected architecture LD scores, component weight sums and separate tagging values. These quantities are not inferred from one another.

Regional methods

With independent LD-region sets already saved by gbase::gprep(task="eigen", sets=sets), HESS uses the same Glist/stat interface for one or two traits:

fit <- gcorr(Glist, stat, method="hess", nthreads=4)
fit$estimates

The stored region definitions and reference samples are reused. Named Z/N matrices must already be oriented to Glist’s a1 alleles; checked frames retain coordinate and allele validation. Missing statistic markers are omitted from stored regions and empty regions are dropped, with counts and names recorded in fit$projection$reference_selection. Spectra alone cannot project Z, so HESS recomputes the required vectors per region, projects both traits together, then discards the vectors. The result keeps compact moments in fit$projection. The overlap inputs shared_sample_size and phenotypic_correlation both default to zero; override them for overlapping samples. Fits retain the existing joint correction, spectral and uncertainty defaults.

Alternatively, prepare allele-labelled frames and explicit independent LD-region sets. For example, the EUR boundaries can define the sets once:

ids <- intersect(rownames(stat$z), rownames(stat$n))
ids <- intersect(ids, unlist(Glist$rsids, use.names = FALSE))
stat_hess <- setNames(lapply(c("GWAS1", "GWAS2"), function(id)
  gact::getMarkerStat(GAlist=GAlist, studyID=id, rsids=ids)), c("GWAS1", "GWAS2"))
boundaries <- read.table(
  "https://psoerensen.github.io/gact/Document/data/EUR-fourier-ac125e47.bed",
  header = TRUE)
boundaries$chr <- sub("chr", "", boundaries$chr, fixed = TRUE)
markers <- stat_hess[[1]]
membership <- integer(nrow(markers))
for (chr in unique(boundaries$chr)) {
  rows <- which(as.character(markers$chr) == chr)
  regions_in_chr <- which(boundaries$chr == chr)
  intervals <- boundaries[regions_in_chr, ]
  at <- findInterval(markers$pos[rows], intervals$start)
  covered <- at > 0L
  covered[covered] <- markers$pos[rows[covered]] < intervals$stop[at[covered]]
  membership[rows[covered]] <- regions_in_chr[at[covered]]
}
regions <- split(markers$rsids[membership > 0L], membership[membership > 0L])
saveRDS(regions, "EUR_LD_region_marker_sets.rds")

The boundaries use start-inclusive, stop-exclusive coordinates, matching the CAD/T2D example. Omit regions containing no retained reference markers. These sets jointly define the modeled universe; do not use whole chromosomes as individual dense HESS regions.

For joint HESS with limited disk space, prepare compact moments directly from the existing BED reference:

projection <- gcorr_prepare_hess(
  list(CAD=cad_stat, T2D=t2d_stat), Glist, regions=ld_regions, nthreads=4)
cad <- gcorr_fit_hess(projection, "CAD")
pair <- gcorr_fit_hess(projection, c("CAD","T2D"),
  shared_sample_size=n_shared, phenotypic_correlation=r_phenotype)
saveRDS(projection, "regional-hess-projections.rds")
pair$diagnostics$conditioning
pair$diagnostics$inference$traits
pair$diagnostics$inference$negative_regions

gcorr_prepare_hess() matches region names to physical BED indices once. Descriptor validation retains one component’s full BIM metadata at a time; the numerical reader uses the validated integer sample and marker rows without rebuilding chromosome-wide metadata or marker-name maps. Preparation-only buffers are released before numerical decoding. It computes centered, unit-norm genotypes, decomposes the smaller marker- or sample-space Gram and projects all supplied traits through the same retained components. Large genotype, Gram and eigenvector buffers are discarded after each region. The serializable object retains integer physical BED indices, marker counts, ranks, spectra, trait sample sizes and quadratic/cross moments. Its spectral controls are fixed; prepare again to change truncation or summary statistics. The subsequent native fit receives only regional counts, ranks and moments and reuses the existing joint correction and conditional uncertainty formulas. Schema-1 RDS projections with marker names remain readable; new schema-2 projections store integer membership.

The compact preparation defaults to maximum_components="auto": a common per-region cap of min(50, floor(min(N)/(2*number_of_regions))). Its conservative total-rank bound is at most half the smallest supplied sample size, so the heritability correction denominator is at least 0.5. The policy uses sample sizes and region counts before examining trait projections. If even one component per region exceeds this budget, automatic preparation stops before reading the reference. Numeric caps, including zero for all eligible components, remain explicit overrides. Their conservative bound warns before preparation if it exceeds the half-sample budget; fits additionally warn on actual rank fractions above 0.5 and retain the native singularity guard. Older schema-1 projections remain readable and keep their original ranks. The policy and actual conditioning are recorded separately. A smaller rank can lose signal: this safeguard does not establish unbiased estimates or statistical calibration. The independent prepared-LD fitters retain their existing numeric defaults.

Genomic-control adjustments are explicit fit-time inputs. For example, after establishing the factor applied to a study:

cad <- gcorr_fit_hess(projection, "CAD",
  genomic_control=c(CAD=applied_gc_factor),
  genomic_control_provenance=c(CAD="Study methods: source of applied factor"))
cad$input_scaling

The factor multiplies Z-squared; a two-trait fit requires both named factors and scales its cross moment by their geometric mean. Every nonunit adjustment requires named provenance, either a documented source or an explicitly labelled sensitivity assumption. Each fit starts from unchanged projection moments, so repeated fitting cannot accumulate adjustments and needs no new eigensystem. The default factor is one. No correction is inferred from an LDSC intercept, and supplied sample sizes are unchanged. This follows the explicit-factor approach in the HESS guidance; an assumed factor is not a validated correction.

Regions must be disjoint and collectively define the modeled genetic universe. Each region must contain markers from one chromosome. One constant positive sample size per trait is required. Preparation matches the supplied regional markers and explicitly swaps alleles to BIM allele-1 dosage when necessary; strand complements and sample-size transformations are not inferred. Two-trait fits require a shared-sample count and phenotypic correlation, which are different from a cross-LDSC intercept. Binary-trait liability conversion is not performed. Negative estimates remain visible, with reasons for missing uncertainty.

diagnostics$conditioning reports total rank, rank/N, correction denominators and amplification. diagnostics$inference reports negative regional h2 by trait, total h2 outside [0,1], indefinite regional genetic matrices, invalid total residual matrices and missing-SE counts.

Compact fits default to uncertainty="hess". This uses the pinned original HESS software variance systems for heritability and covariance, including their full regional LD rank rather than the retained projection rank. The diagonal-plus-rank-one systems are solved in linear storage and arithmetic. Genome-wide correlation uses the original delete-one-region jackknife, centered on the full-sample ratio of summed regional contributions. It does not refit the coupled equations after deletion. Invalid delete-region ratios remain missing and invalidate the SE; none is replaced with zero. An rg outside [-1,1] retains a descriptive jackknife SE but has no Wald Z/P. diagnostics$genome_rg_jackknife records the replicates and their validity.

Point estimates retain the documented row-rank moment equations. The original software instead uses a column-rank convention, so unequal ranks can cause small point differences. Original variance parity is assessed at identical point inputs. This convention does not supply a joint sampling covariance or regional rg SEs. Select uncertainty="moment" for the previous signed conditional fourth-moment covariance and regional delta SEs, or "strict" for its admissibility checks. Prepared-LD fitters retain their strict default. The h2 test uses the directional upper tail, including for negative estimates. These SEs remain approximate: changing their calculation does not restore discarded signal or establish the scale of genomic-control-adjusted statistics.

memory_limit_bytes bounds a conservative estimate of numerical workspace for one region; reference metadata, summary inputs, R and allocator overhead are additional. The default smaller-Gram dimension limit is 2048. This route needs neither stored sparse LD nor a disk cache of eigenvectors. It retains separate reader-setup and per-region decode/product/eigen/projection clocks.

gcorr_prepare_regions() reads named marker regions directly from Glist sparse LD through one native reader. It validates each complete eigensystem once; the resulting session-bound preparation is accepted by:

• gcorr_hess() and gcorr_rho_hess();
• gcorr_supergnova();
• gcorr_lava(), gcorr_lava_partial() and gcorr_lava_regression();
• gcorr_hdll().

regions <- gcorr_prepare_regions(Glist,
  regions = list(locus = markers[1:4]))
regional <- gcorr_hess(one, regions, maximum_components = 0,
  minimum_eigenvalue = 0)

Named signed-LD matrices plus their reference sample sizes remain supported when a Glist is not available.

The summary data frame must contain every regional marker. LAVA sampling correlations, binary-trait metadata and simulation seeds are explicit. SUPERGNOVA exposes covariance rather than correlation. HDL-L accepts exactly two continuous traits per fit. Joint HESS/rho-HESS requires disjoint prepared regions and reports its correction explicitly.

Capability table

Family	Focused entry point	Essential explicit inputs
Convenient LDSC workflow	ldsc()	checked stat; Glist or named ordinary LD scores; optional named marker sets using stratified or grouped scores
Ordinary LDSC	gcorr_ldsc()	named LD scores, delete blocks, marker universe
Partitioned LDSC	gcorr_ldsc_partitioned()	annotation/weight scores, sums, blocks, reference metadata
GNOVA	gcorr_gnova()	ordinary scores or binary annotation inputs; pair overlap intercepts
SumHer	gcorr_sumher()	architecture scores and sums, tagging, blocks
HESS/rho-HESS	gcorr(Glist, stat, method="hess"); gcorr_hess(), gcorr_rho_hess(); gcorr_prepare_hess() plus gcorr_fit_hess()	prepared regional LD or compact BED projections; compact fits default to no sample overlap
SUPERGNOVA	gcorr_supergnova()	two traits, regional LD, overlap intercept and variance
LAVA	gcorr_lava*()	regional LD, sampling correlation, optional seeded inference
HDL-L	gcorr_hdll()	two traits, regional LD and shared sample count

Result and uncertainty conventions

Focused fits return labelled estimate tables, method-specific matrices and diagnostics, explicit assumptions and native software versions. Missing uncertainty remains missing; it is never replaced by zero. Joint annotation or architecture modelling uses the retained full block-jackknife covariance. Grouped SEM/GGM and contrasts are described in covariance models.