HEART tests the association between each variant and a phenotype in a target population. It borrows information from a larger source population whose genetic effects may differ. This vignette walks through the two steps of the method on simulated data.
library(HEART)
sim <- heart_simulate(n_tar = 3000, n_src = 20000, J = 3000, seed = 2026)
table(sim$truth$config)
#>
#> 00 01 10 11
#> 2732 110 30 128The configurations are 00 (null in both populations),
01 (source-specific signal, so the target null holds),
10 (target-specific signal) and 11 (shared
signal).
Step 1: target-only, auxiliary and decorrelated estimators
fit <- heart_assoc(sim$Y, sim$G, sim$D, sim$X, aux = c("pooled", "source"))
head(fit[, c("snp", "beta_tar", "p_tar", "beta_pooled", "rho_pooled",
"beta_de_pooled", "p_de_pooled")])
#> snp beta_tar p_tar beta_pooled rho_pooled beta_de_pooled
#> 1 snp1 -0.02314106 0.60962163 -0.012256548 0.09230360 -0.010120545
#> 2 snp2 0.05525188 0.18973649 0.002104585 0.12017348 -0.004535225
#> 3 snp3 -0.05259263 0.54332479 -0.008172641 0.06696513 -0.004650768
#> 4 snp4 0.04210665 0.37582215 -0.014303461 0.13383141 -0.019938654
#> 5 snp5 -0.08442441 0.02783308 0.001396242 0.14270131 0.013443716
#> 6 snp6 -0.08371033 0.26529879 0.031850898 0.13492930 0.043145875
#> p_de_pooled
#> 1 0.4785980
#> 2 0.7565946
#> 3 0.8371745
#> 4 0.2565701
#> 5 0.3557711
#> 6 0.1170738The pooled estimator is not valid for the target null. Its p-values
are small for the source-specific variants (01), whose
target effect is zero:
Step 2: HEART test
res <- heart_test(fit$p_tar, fit$p_de_pooled, snp = fit$snp, alpha = 0.05)
res
#> HEART: Heterogeneity-aware Empowerment via Adaptive Robust Transfer
#> tests: 3000
#> four-group proportions: pi00 = 0.8815, pi01 = 0.0807, pi10 = 0.0263, pi11 = 0.0115
#> tail parameters: alpha1 = 0.129, C1 = 0.838, alpha2 = 0.056, C2 = 0.964
#> discoveries at level 0.05: FDR = 84, FWER = 16 (target-only BH = 24, Bonferroni = 4)Compare the discoveries with the truth:
signal <- sim$truth$beta_tar != 0
bh <- p.adjust(fit$p_tar, "BH") <= 0.05
rbind(
HEART = c(discoveries = sum(res$results$reject_fdr),
FDP = mean(!signal[res$results$reject_fdr]),
power = mean(res$results$reject_fdr[signal])),
`Target-only BH` = c(sum(bh), mean(!signal[bh]), mean(bh[signal]))
)
#> discoveries FDP power
#> HEART 84 0.03571429 0.5126582
#> Target-only BH 24 0.04166667 0.1455696Inspecting the null distribution
heart_null_cdf() returns the components of the estimated
composite null distribution: F00, Q̂, F̂01, N̂ and
F̂0.
nc <- res$null_cdf
head(nc)
#> log_t F00 Q F01 N F0
#> 1 -14.234983 1.768597e-09 1.020412e-08 9.374244e-08 9.127820e-09 9.486187e-09
#> 2 -11.313273 5.056040e-08 2.917139e-07 2.679895e-06 2.609447e-07 2.711897e-07
#> 3 -9.978783 2.338409e-07 1.349171e-06 1.239446e-05 1.206864e-06 1.254247e-06
#> 4 -9.743951 3.061689e-07 1.766476e-06 1.622813e-05 1.580153e-06 1.642191e-06
#> 5 -9.709737 3.184294e-07 1.837214e-06 1.687798e-05 1.643430e-06 1.707953e-06
#> 6 -9.572134 3.729018e-07 2.151499e-06 1.976523e-05 1.924565e-06 2.000125e-06