Batched cross-trait LD-score-regression genetic correlation for the UKB-PPP database.
One external trait is correlated against every protein. Because every protein's z-score is defined on the SAME shared SNP set (see ppp_ldsc_context), we process K proteins at once and run the three regressions each genetic correlation needs:
rg = gcov(trait, protein) / sqrt(h2_trait * h2_protein)
Following the GenomicSEM convention (see genomic_sem_ldsc.run_ldsc), each DIAGONAL heritability is estimated on that trait's own kept SNP set, and only the OFF-DIAGONAL covariance on the pair's intersection. So h2_trait -- its point estimate, its per-block leave-one-out delete values, and n_bar -- depends only on the trait and is computed exactly ONCE (estimate_trait_context) and reused for every protein. Only h2_protein and gcov are per-protein.
Per-protein exclusions (missing variants, the chi-square filter, and cis variants) are encoded as ZERO regression weights over the shared row set rather than by dropping rows, so all proteins keep the same contiguous jackknife blocks -- exactly as in batched_ldsc_h2.
Weighting (the easy bug): the univariate h2 weight is omega = het * oc (the PRODUCT), because estimate_h2 multiplies BOTH design and response by sqrt(het * oc). The covariance weight is omega_cov = oc * (sqrt(het_trait) + sqrt(het_protein))^2, from estimate_cov's weights_cov = (initial_w_trait + initial_w_protein) with initial_w = sqrt(het * oc); the sum-to-1 normalization is a common scalar that cancels in the 2x2 slope solve. The heteroskedasticity aggregate-h2 that shapes het is computed over the SAME set the regression runs on: own kept set for a diagonal h2, the pair intersection for the covariance.
The rg standard error is a delete-block jackknife on the ratio: leave-one-block-out slopes for gcov and h2_protein are combined with the once-computed h2_trait deletes at the same block index to form per-block rg values, and rg_se is their jackknife spread. This is the classic ldsc rg SE.
Classes:
-
BatchedRgResult–Per-protein cross-trait rg outputs; all arrays are length K (protein order preserved).
-
ExactRgResult– -
TraitLdscContext–The trait-only LDSC quantities, computed once and reused for every protein.
Functions:
-
batched_rg–Estimate the trait-vs-protein genetic correlation for K proteins sharing the SNP set.
-
estimate_trait_context–Estimate the trait heritability once over the shared SNP set and package the trait-only
-
exact_rg_single–Reference single-protein genetic correlation via the repo's exact GenomicSEM port: h2 on
Per-protein cross-trait rg outputs; all arrays are length K (protein order preserved). h2_trait and n_bar_trait are trait-level scalars shared by every protein.
Attributes:
-
gcov(ndarray) – -
gcov_intercept(ndarray) – -
h2_protein(ndarray) – -
h2_trait(float) – -
n_bar_protein(ndarray) – -
n_bar_trait(float) – -
n_snps(ndarray) – -
rg(ndarray) – -
rg_p(ndarray) – -
rg_se(ndarray) – -
rg_z(ndarray) –
Attributes:
-
gcov(float) – -
gcov_intercept(float) – -
h2_protein(float) – -
h2_trait(float) – -
rg(float) –
The trait-only LDSC quantities, computed once and reused for every protein.
z / keep / n are parallel (S,) arrays over the shared context SNP set: z is the signed trait z-score (index-effect-allele oriented, zeroed off the kept set), keep marks the trait's kept SNPs (present and passing the chi-square filter), and n is the per-SNP trait sample size. h2 / n_bar are the trait heritability and its N_bar; h2_delete holds the (n_blocks,) per-block leave-one-out h2_trait values used by the rg jackknife.
Attributes:
batched_rg(
trait_ctx: TraitLdscContext,
z_protein: ndarray,
ld: ndarray,
n_protein: ndarray,
m: float,
*,
n_blocks: int = DEFAULT_N_BLOCKS,
exclude: ndarray | None = None,
) -> BatchedRgResult
Estimate the trait-vs-protein genetic correlation for K proteins sharing the SNP set.
trait_ctx: the once-computed trait LDSC context (estimate_trait_context). z_protein: (S, K) signed protein z-score (= BETA/SE, index-effect-allele oriented), NaN where a protein lacks the variant. ld: (S,) LD score, genome-sorted. n_protein: (K,) per-protein sample size. m: total reference-SNP count. exclude: optional (S, K) boolean; True drops that variant for that protein (e.g. cis).
Source code in mecfs_bio/build_system/task/ppp_ldsc/batched_ldsc_rg.py
288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 | |
estimate_trait_context(
z_trait: ndarray,
n_trait: ndarray,
ld: ndarray,
m: float,
*,
n_blocks: int = DEFAULT_N_BLOCKS,
) -> TraitLdscContext
Estimate the trait heritability once over the shared SNP set and package the trait-only quantities the batched rg kernel reuses for every protein.
z_trait: (S,) signed trait z-score aligned to the context (index) effect allele, NaN where the trait lacks the variant. n_trait: (S,) trait sample size (NaN where absent). ld: (S,) LD score, genome-sorted for contiguous blocks. m: total reference-SNP count.
Source code in mecfs_bio/build_system/task/ppp_ldsc/batched_ldsc_rg.py
exact_rg_single(
z_trait: ndarray,
n_trait: float,
z_protein: ndarray,
ld: ndarray,
n_protein: float,
m: float,
*,
n_blocks: int = DEFAULT_N_BLOCKS,
exclude: ndarray | None = None,
) -> ExactRgResult
Reference single-protein genetic correlation via the repo's exact GenomicSEM port: h2 on each trait's own kept set and gcov on the intersection (dropping, not zero-weighting, the filtered/cis SNPs). This is what the batched kernel's point estimates are validated against. z_trait/z_protein/ld are genome-sorted (S,) arrays; n_trait/n_protein are scalar.