Marker-QTL Unit Congruency

πŸ“ Theory: Theory overview

Objective

We developed a theory called the Marker-QTL Unit Congruence Theory, where:

Genomic prediction accuracy is maximised when the marker unit is aligned with the biological QTL unit.

In other words, if the QTL effects actually live at the haplotype-block level, microhaplotype (MH) markers are expected to outperform single-SNP markers, and vice versa. To test this, we ran a proof-of-concept simulation under both QTL scenarios and compared four genomic prediction models across two marker representations.

Study design

A forward simulation over a fully crossed grid.

Factor Levels
QTL scenario QTL@SNP (true genetic effects at single-SNP level), QTL@MH (effects at haplotype-block level)
Marker representation SNP (z-score encoding via construct_snp_matrix(encoding="zscore")), MH (multi-allelic \(\mathbf{W}_{\alpha h}\) coding, Da 2015, via construct_wah_matrix)
Model BayesR, BayesA, GBLUP, GWABLUP
Trait type Continuous, binary (threshold at the training-set liability median)
Total \(2 \times 2 \times 4 \times 2 = 32\) model fits

Population parameters: \(n_{\text{total}} = 300\) (training 200, test 100), 50 haplotype blocks Γ— 2 SNPs/block = 100 SNPs total, 10 QTL, target heritability \(h^2 = 0.3\). SNP MAF \(\sim\) Unif(0.1, 0.5). In both scenarios the microhaplotype matrix is derived from the same SNP data via phasing and \(\mathbf{W}_{\alpha h}\) encoding, so any difference in performance reflects the marker representation alone, not the underlying genotypes.

We ran MCMC for 20,000 iterations with a 10,000 burn-in, thinning every 10 samples, and seed 123. Prior hyperparameters follow the masbayes defaults, where for BayesA (\(\nu = 4.5\)) and BayesR (\(\pi = (0.90,\, 0.05,\, 0.03,\, 0.02)\), \(\text{variance\_class} = (0,\, 0.01,\, 0.1,\, 1)\)).

Metrics and interpretation

Let \(\hat{u}_i = \widehat{\text{GEBV}}_i\) for individual \(i\), \(y_i\) = observed phenotype, \(g_i\) = true breeding value (TBV), \(\hat{z}_i\) = posterior mean liability (binary trait), \(\Phi(\cdot)\) = standard normal CDF.

Metric Formula Range Interpretation
\(r_{\text{train}}\) \(\displaystyle r_{\text{train}} = \frac{\mathrm{Cov}(\hat{u},\, y)}{\sqrt{\mathrm{Var}(\hat{u})\,\mathrm{Var}(y)}}\) on the training set. For binary traits, replace \(y\) with \(\hat{z}\). \([-1,\, 1]\) Training fit. Closer to 1 is better, but values > 0.95 can signal overfitting.
\(r_{\text{test},y}\) \(\displaystyle r_{\text{test},y} = \frac{\mathrm{Cov}(\hat{u}_{\text{test}},\, y_{\text{test}})}{\sqrt{\mathrm{Var}(\hat{u}_{\text{test}})\,\mathrm{Var}(y_{\text{test}})}}\) \([-1,\, 1]\) Predictive ability on held-out data. 1 = ideal.
\(r_{\text{test},g}\) \(\displaystyle r_{\text{test},g} = \frac{\mathrm{Cov}(\hat{u}_{\text{test}},\, g_{\text{test}})}{\sqrt{\mathrm{Var}(\hat{u}_{\text{test}})\,\mathrm{Var}(g_{\text{test}})}}\) \([-1,\, 1]\) Accuracy against TBV β€” the primary metric, available only in simulation. Closer to 1 is better.
bias Continuous: slope \(\hat{b}\) of \(y_i = a + b\,\hat{u}_i + \varepsilon_i\), \(\displaystyle \hat{b} = \frac{\mathrm{Cov}(\hat{u},\, y)}{\mathrm{Var}(\hat{u})}\). Binary: slope \(\hat{b}\) of \(y_i = a + b\,\Phi(\hat{u}_i) + \varepsilon_i\). around 1 \(\hat{b} = 1\) is ideal; \(\hat{b} < 1\) indicates over-dispersion (predictions too extreme); \(\hat{b} > 1\) indicates under-dispersion (predictions too compressed).
\(\hat{h}^2_{\text{post}}\) \(\displaystyle \hat{h}^2_{\text{post}} = \frac{\hat{\sigma}^2_g}{\hat{\sigma}^2_g + \hat{\sigma}^2_e}\) with variance components taken from the posterior. \([0,\, 1]\) Posterior heritability estimate. Compare with the target \(h^2 = 0.3\).
AUC \(\displaystyle \text{AUC} = \Pr\!\left(\hat{u}_i > \hat{u}_j \,\middle|\, y_i = 1,\; y_j = 0\right)\) \([0.5,\, 1]\) Binary only. 0.5 = random, 1.0 = perfect. > 0.7 is good, > 0.9 is very good. Rank-invariant.

Results

Show full R script
library(masbayes)
library(masreml)

# ── CONFIG ───────────────────────────────────────────────────────────────────
config <- list(
  seed            = 42,
  n_total         = 300,
  n_train         = 200,
  n_test          = 100,
  n_blocks        = 50,
  n_snp_per_block = 2,
  h2_target       = 0.3,
  n_qtl           = 10,
  mcmc = list(
    n_iter = 20000L,
    n_burn = 10000L,
    n_thin = 10L,
    seed   = 123L
  ),
  bayesr = list(
    pi_vec         = c(0.90, 0.05, 0.03, 0.02),
    variance_class = c(0, 0.01, 0.1, 1),
    a0_e           = 5,
    a0_g           = 5
  ),
  bayesa = list(
    nu   = 4.5,
    a0_e = 10
  )
)

set.seed(config$seed)
n_total         <- config$n_total
n_train         <- config$n_train
n_test          <- config$n_test
n_blocks        <- config$n_blocks
n_snp_per_block <- config$n_snp_per_block
n_snp_total     <- n_blocks * n_snp_per_block

# ── 1. Generate SNP genotype ─────────────────────────────────────────────────
maf <- runif(n_snp_total, 0.1, 0.5)
geno_snp_all <- matrix(0L, nrow = n_total, ncol = n_snp_total)
for (j in 1:n_snp_total)
  geno_snp_all[, j] <- as.integer(rbinom(n_total, 2, maf[j]))

# ── 2. Phase SNP -> haplotype ────────────────────────────────────────────────
hap_all <- matrix(0L, nrow = n_total, ncol = n_snp_total * 2)
for (j in 1:n_snp_total) {
  for (i in 1:n_total) {
    g  <- geno_snp_all[i, j]
    h1 <- if (g == 2) 1L else if (g == 0) 0L else as.integer(rbinom(1, 1, 0.5))
    h2 <- as.integer(g - h1)
    hap_all[i, 2*j-1] <- h1 + 1L
    hap_all[i, 2*j  ] <- h2 + 1L
  }
}
storage.mode(hap_all) <- "integer"

# ── 3. Reorder hap columns by block ─────────────────────────────────────────
hap_cols_per_block <- n_snp_per_block * 2
hap_reordered <- matrix(0L, nrow = n_total, ncol = n_snp_total * 2)
col_out <- 1
for (b in 1:n_blocks) {
  for (j in ((b-1)*n_snp_per_block + 1):(b*n_snp_per_block)) {
    hap_reordered[, col_out]   <- hap_all[, 2*j-1]
    hap_reordered[, col_out+1] <- hap_all[, 2*j  ]
    col_out <- col_out + 2
  }
}
storage.mode(hap_reordered) <- "integer"

# ── 4. Encode MH per block ───────────────────────────────────────────────────
idx_train <- 1:n_train
idx_test  <- (n_train+1):n_total

encode_hap <- function(mat)
  apply(mat, 1, function(x) sum(x * 3^(seq_along(x)-1)))

hap_block_all    <- matrix(0L, nrow = n_total, ncol = n_blocks * 2)
allele_freq_list <- list(haplotype=c(), allele=c(), freq=c())

for (b in 1:n_blocks) {
  cols    <- ((b-1)*hap_cols_per_block + 1):(b*hap_cols_per_block)
  hap_sub <- hap_reordered[, cols]
  h1_id   <- encode_hap(hap_sub[, seq(1, hap_cols_per_block, 2), drop=FALSE])
  h2_id   <- encode_hap(hap_sub[, seq(2, hap_cols_per_block, 2), drop=FALSE])
  hap_block_all[, 2*b-1] <- h1_id
  hap_block_all[, 2*b  ] <- h2_id
  tbl     <- table(c(h1_id[idx_train], h2_id[idx_train]))
  freqs   <- as.numeric(tbl) / sum(tbl)
  alleles <- as.integer(names(tbl))
  allele_freq_list$haplotype <- c(allele_freq_list$haplotype,
                                   rep(paste0("block_", b), length(alleles)))
  allele_freq_list$allele    <- c(allele_freq_list$allele, alleles)
  allele_freq_list$freq      <- c(allele_freq_list$freq, freqs)
}
storage.mode(hap_block_all) <- "integer"
colnames_block <- paste0("block_", rep(1:n_blocks, each = 2))

# ── 5. Construct W matrices ──────────────────────────────────────────────────
wah_train  <- construct_wah_matrix(
  hap_block_all[idx_train,], colnames_block, allele_freq_list, NULL, TRUE)
W_mh_train <- wah_train$W_ah
ref_struct  <- list(allele_info=wah_train$allele_info,
                    dropped_alleles=wah_train$dropped_alleles)
W_mh_test  <- construct_wah_matrix(
  hap_block_all[idx_test,], colnames_block, NULL, ref_struct, TRUE)$W_ah

snp_train_obj <- construct_snp_matrix(geno_snp_all[idx_train, ],
                                      encoding = "zscore")
W_snp_train   <- snp_train_obj$W
p_snp_tr      <- snp_train_obj$freq
sd_snp_tr     <- snp_train_obj$sd

W_snp_test    <- construct_snp_matrix(geno_snp_all[idx_test, ],
                                      encoding = "zscore",
                                      ref_freq = p_snp_tr,
                                      ref_sd   = sd_snp_tr)$W
W_snp_all     <- construct_snp_matrix(geno_snp_all,
                                      encoding = "zscore",
                                      ref_freq = p_snp_tr,
                                      ref_sd   = sd_snp_tr)$W
storage.mode(W_snp_train) <- "double"
storage.mode(W_snp_test)  <- "double"
storage.mode(W_snp_all)   <- "double"

wah_all    <- construct_wah_matrix(
  hap_block_all, colnames_block, allele_freq_list, NULL, TRUE)
W_mh_all   <- wah_all$W_ah

# ── 6. Simulate y (continuous and binary) ───────────────────────────────────
simulate_y <- function(W_source, idx_tr, label, h2_target = 0.3,
                       n_qtl = 10, type = "snp") {
  n_col     <- ncol(W_source)
  beta_true <- rep(0, n_col)
  qtl_idx   <- sample(n_col, n_qtl)
  if (type == "snp") {
    raw <- rgamma(n_qtl, shape=0.4, scale=1) * sample(c(-1,1), n_qtl, replace=TRUE)
    beta_true[qtl_idx] <- raw / sqrt(sum(raw^2))
  } else {
    raw <- rnorm(n_qtl)
    beta_true[qtl_idx] <- raw / sqrt(sum(raw^2))
  }
  tbv_all  <- as.vector(W_source %*% beta_true)
  tbv_mean <- mean(tbv_all[idx_tr])
  tbv_sd   <- sd(tbv_all[idx_tr])
  tbv_std  <- (tbv_all - tbv_mean) / tbv_sd
  sg <- var(tbv_std[idx_tr])
  se <- sg * (1 - h2_target) / h2_target
  y_cont <- tbv_std + rnorm(length(tbv_std), 0, sqrt(se))
  threshold <- median(y_cont[idx_tr])
  y_bin  <- as.numeric(y_cont > threshold)
  h2_obs <- sg / (sg + se)
  list(y_cont=y_cont, y_bin=y_bin, g=tbv_std, sigma2_g=sg, sigma2_e=se, h2=h2_obs)
}

set.seed(config$seed)
sc_snp <- simulate_y(W_snp_all, idx_train, "QTL@SNP",
                     h2_target = config$h2_target,
                     n_qtl     = config$n_qtl,
                     type      = "snp")
set.seed(config$seed)
sc_mh  <- simulate_y(W_mh_all,  idx_train, "QTL@MH",
                     h2_target = config$h2_target,
                     n_qtl     = config$n_qtl,
                     type      = "mh")

# ── 7. Bayesian model fitting ────────────────────────────────────────────────
mcmc_p <- config$mcmc

run_scenario <- function(sc, W_tr, W_te, y_tr, y_te, g_te,
                         marker_label, trait_type = "continuous") {
  y_train <- y_tr
  y_test  <- y_te
  resp    <- if (trait_type == "binary") "binary" else "gaussian"
  mtype   <- if (identical(marker_label, "SNP")) "snp" else "multiallelic"
  rows    <- list()
  se_init <- if (trait_type == "binary") 1.0 else sc$sigma2_e
  sg_init <- if (trait_type == "binary") 1.0 else sc$sigma2_g

  for (model in c("BayesR", "BayesA")) {
    result <- tryCatch({
      if (model == "BayesR") {
        res <- run_bayesr(
          w=W_tr, y=y_train, 
          marker_type   = mtype,
          pi_vec        = config$bayesr$pi_vec,
          sigma2_e_init = se_init,
          sigma2_ah     = sg_init,
          prior_params  = list(
            a0_e           = config$bayesr$a0_e,
            a0_g           = config$bayesr$a0_g,
            variance_class = config$bayesr$variance_class
          ),
          mcmc_params   = mcmc_p,
          method        = "mcmc",
          response_type = resp,
          fold_id       = 0L)
      } else {
        res <- run_bayesa(
          w=W_tr, y=y_train, 
          marker_type   = mtype,
          nu            = config$bayesa$nu,
          sigma2_g      = sg_init,
          sigma2_e_init = se_init,
          prior_params  = list(a0_e = config$bayesa$a0_e),
          mcmc_params   = mcmc_p,
          method        = "mcmc",
          response_type = resp,
          fold_id       = 0L)
      }

      gebv_tr  <- res$GEBV
      pred_te  <- predict(res, W_te, y_test)
      gebv_te  <- pred_te$GEBV
      r_test_y <- cor(gebv_te, y_test)
      r_test_g <- cor(gebv_te, g_te)

      if (trait_type == "binary" && !is.null(res$z_hat) && is.numeric(res$z_hat)) {
        r_train       <- cor(gebv_tr, res$z_hat)
        bias_test     <- pred_te$metrics$bias
        sigma2_g_post <- var(gebv_tr)
        h2_post       <- sigma2_g_post / (sigma2_g_post + 1.0)
      } else {
        r_train       <- cor(gebv_tr, y_train)
        bias_test     <- pred_te$metrics$bias
        sigma2_g_post <- var(gebv_tr)
        sigma2_e_post <- mean(res$sigma2_e_samples)
        h2_post       <- sigma2_g_post / (sigma2_g_post + sigma2_e_post)
      }

      auc <- if (trait_type == "binary") pred_te$metrics$AUC else NA

      list(status="OK", r_train=round(r_train,3),
           r_test_y=round(r_test_y,3), r_test_g=round(r_test_g,3),
           bias=round(bias_test,3), h2=round(h2_post,3),
           auc=round(auc,3), p=ncol(W_tr))
    }, error = function(e)
      list(status=paste("ERROR:", conditionMessage(e)),
           r_train=NA, r_test_y=NA, r_test_g=NA,
           bias=NA, h2=NA, auc=NA, p=ncol(W_tr)))

    rows[[model]] <- data.frame(
      Trait=trait_type, Marker=marker_label, Model=model,
      Status=result$status, r_train=result$r_train,
      r_test_y=result$r_test_y, r_test_g=result$r_test_g,
      bias=result$bias, h2_post=result$h2, AUC=result$auc,
      p=result$p, stringsAsFactors=FALSE)
  }
  do.call(rbind, rows)
}

# ── 7b. GBLUP via masreml ────────────────────────────────────────────────────
rownames(geno_snp_all) <- as.character(1:n_total)
G_snp_full <- build_G_snp(geno_snp_all, ref_W = geno_snp_all[idx_train, ])
rownames(hap_block_all) <- as.character(1:n_total)
G_mh_full <- build_G_mh(
  mh_list = hap_block_all,
  ref_mh  = hap_block_all[idx_train, ],
  ids     = as.character(1:n_total)
)
train_ids_ch <- as.character(idx_train)
test_ids_ch  <- as.character(idx_test)

run_gblup_scenario <- function(sc, marker_label, G_full, y_tr, y_te, g_te,
                                trait_type = "continuous") {
  G_tr <- G_full[train_ids_ch, train_ids_ch]
  y_named <- setNames(y_tr, train_ids_ch)
  fit <- tryCatch(
    masreml(y = y_named, G = list(g = G_tr),
            trait = trait_type, method = "auto"),
    error = function(e) NULL
  )
  if (is.null(fit)) {
    return(data.frame(Trait=trait_type, Marker=marker_label, Model="GBLUP",
                      Status="ERROR", r_train=NA, r_test_y=NA, r_test_g=NA,
                      bias=NA, h2_post=NA, AUC=NA, p=NA, stringsAsFactors=FALSE))
  }
  pred <- predict(fit, G_full = list(g = G_full),
                  train_ids = train_ids_ch, test_ids = test_ids_ch)
  gebv_tr  <- fit$total_gebv + fit$fixed_effects[1]
  gebv_te  <- pred$GEBV      + fit$fixed_effects[1]
  h2_post  <- as.numeric(fit$varcomp$h2["g"])
  ev <- evaluate_prediction(
          gebv        = gebv_te,
          y           = y_te,
          h2          = h2_post,
          tbv         = g_te,
          fitted_prob = if (trait_type == "binary") pred$prob else NULL
        )
  r_train  <- cor(gebv_tr, y_tr)
  data.frame(Trait=trait_type, Marker=marker_label, Model="GBLUP",
             Status="OK",
             r_train=round(r_train,3), r_test_y=ev$r_test_y,
             r_test_g=ev$r_test_g, bias=round(ev$bias,3),
             h2_post=round(h2_post,3), AUC=ev$AUC,
             p=nrow(G_tr), stringsAsFactors=FALSE)
}

gblup_results <- list()
for (sc_name in c("QTL@SNP", "QTL@MH")) {
  sc   <- if (sc_name == "QTL@SNP") sc_snp else sc_mh
  g_te <- sc$g[idx_test]
  for (trait in c("continuous", "binary")) {
    y_tr <- if (trait == "binary") sc$y_bin[idx_train] else sc$y_cont[idx_train]
    y_te <- if (trait == "binary") sc$y_bin[idx_test]  else sc$y_cont[idx_test]
    r_snp <- run_gblup_scenario(sc, "SNP", G_snp_full, y_tr, y_te, g_te, trait)
    r_mh  <- run_gblup_scenario(sc, "MH",  G_mh_full,  y_tr, y_te, g_te, trait)
    r_snp$Scenario <- r_mh$Scenario <- sc_name
    gblup_results[[paste(sc_name, trait)]] <- rbind(r_snp, r_mh)
  }
}
gblup_final <- do.call(rbind, gblup_results)
gblup_final <- gblup_final[, c("Scenario","Trait","Marker","Model",
                                "r_train","r_test_y","r_test_g",
                                "bias","h2_post","AUC","p","Status")]

# ── 7c. GWABLUP via masreml ──────────────────────────────────────────────────
run_gwablup_scenario <- function(sc, marker_type, y_tr, y_te, g_te,
                                  G_full, geno_train, geno_all,
                                  trait_type = "continuous") {
  train_ids <- train_ids_ch
  test_ids  <- test_ids_ch
  y_named   <- setNames(y_tr, train_ids)
  G_tr      <- G_full[train_ids, train_ids]
  comp_name <- if (marker_type == "SNP") "snp_add" else "mh_add"

  result <- tryCatch({
    fit_tr <- masreml(y = y_named, G = list(g = G_tr),
                      trait = trait_type, method = "auto")
    if (marker_type == "SNP") {
      markers_tr  <- list(snp_add = geno_train)
      ref_markers <- list(snp_add = geno_train)
    } else {
      hap_tr      <- hap_block_all[idx_train, ]
      rownames(hap_tr) <- train_ids
      markers_tr  <- list(mh_add = hap_tr)
      ref_markers <- list(mh_add = hap_tr)
    }
    gwas_tr <- run_gwas(
      markers     = markers_tr,
      y           = y_named,
      masreml_fit = fit_tr,
      ref_markers = ref_markers
    )
    fit_wa <- gwablup(
      y           = y_named,
      markers     = markers_tr,
      gwas_result = gwas_tr,
      trait       = trait_type,
      ref_markers = ref_markers
    )
    g_full_named        <- list(G_full)
    names(g_full_named) <- comp_name
    pred <- predict(fit_wa,
                    G_full    = g_full_named,
                    train_ids = train_ids,
                    test_ids  = test_ids)
    gebv_tr <- fit_wa$total_gebv + fit_wa$fixed_effects[1]
    gebv_te <- pred$GEBV         + fit_wa$fixed_effects[1]
    h2_post <- as.numeric(fit_wa$varcomp$h2[comp_name])
    ev <- evaluate_prediction(
            gebv        = gebv_te,
            y           = y_te,
            h2          = h2_post,
            tbv         = g_te,
            fitted_prob = if (trait_type == "binary") pred$prob else NULL
          )
    list(status="OK",
         r_train  = round(cor(gebv_tr, y_tr), 3),
         r_test_y = ev$r_test_y,
         r_test_g = ev$r_test_g,
         bias     = ev$bias,
         h2_post  = round(h2_post, 3),
         auc      = ev$AUC,
         p        = length(fit_wa$total_gebv))
  }, error = function(e)
    list(status=paste("ERROR:", conditionMessage(e)),
         r_train=NA, r_test_y=NA, r_test_g=NA,
         bias=NA, h2_post=NA, auc=NA, p=NA))

  data.frame(Trait=trait_type, Marker=marker_type, Model="GWABLUP",
             Status=result$status,
             r_train=result$r_train, r_test_y=result$r_test_y,
             r_test_g=result$r_test_g, bias=result$bias,
             h2_post=result$h2_post, AUC=result$auc,
             p=result$p, stringsAsFactors=FALSE)
}

geno_train <- geno_snp_all[idx_train, ]
rownames(geno_train) <- train_ids_ch
storage.mode(geno_train) <- "double"

gwablup_results <- list()
for (sc_name in c("QTL@SNP", "QTL@MH")) {
  sc   <- if (sc_name == "QTL@SNP") sc_snp else sc_mh
  g_te <- sc$g[idx_test]
  for (trait in c("continuous", "binary")) {
    y_tr <- if (trait == "binary") sc$y_bin[idx_train] else sc$y_cont[idx_train]
    y_te <- if (trait == "binary") sc$y_bin[idx_test]  else sc$y_cont[idx_test]
    r_snp <- run_gwablup_scenario(sc, "SNP", y_tr, y_te, g_te,
                                   G_snp_full, geno_train, geno_snp_all, trait)
    r_mh  <- run_gwablup_scenario(sc, "MH",  y_tr, y_te, g_te,
                                   G_mh_full,  geno_train, geno_snp_all, trait)
    r_snp$Scenario <- r_mh$Scenario <- sc_name
    gwablup_results[[paste(sc_name, trait)]] <- rbind(r_snp, r_mh)
  }
}
gwablup_final <- do.call(rbind, gwablup_results)
gwablup_final <- gwablup_final[, c("Scenario","Trait","Marker","Model",
                                    "r_train","r_test_y","r_test_g",
                                    "bias","h2_post","AUC","p","Status")]

# ── 8. Run all Bayesian combinations ─────────────────────────────────────────
all_results <- list()
for (sc_name in c("QTL@SNP", "QTL@MH")) {
  sc   <- if (sc_name == "QTL@SNP") sc_snp else sc_mh
  g_te <- sc$g[idx_test]
  for (trait in c("continuous", "binary")) {
    y_tr_use <- if (trait == "binary") sc$y_bin[idx_train] else sc$y_cont[idx_train]
    y_te_use <- if (trait == "binary") sc$y_bin[idx_test]  else sc$y_cont[idx_test]
    r_mh  <- run_scenario(sc, W_mh_train, W_mh_test,
                          y_tr_use, y_te_use,
                          g_te, "MH", trait)
    r_snp <- run_scenario(sc, W_snp_train, W_snp_test,
                          y_tr_use, y_te_use,
                          g_te, "SNP", trait)
    res        <- rbind(r_mh, r_snp)
    res$Scenario <- sc_name
    all_results[[paste(sc_name, trait)]] <- res
  }
}

# ── 9. Combined results ──────────────────────────────────────────────────────
final <- rbind(do.call(rbind, all_results), gblup_final, gwablup_final)
final <- final[, c("Scenario","Trait","Marker","Model",
                   "r_train","r_test_y","r_test_g",
                   "bias","h2_post","AUC","p","Status")]
final <- final[order(final$Scenario, final$Trait, final$Marker, final$Model), ]
Combined results: 2 QTL scenarios Γ— 2 traits Γ— 2 markers Γ— 4 models = 32 model fits. Train n = 200, test n = 100, hΒ² target = 0.3, n_QTL = 10.
Scenario Trait Marker Model r_train r_test_y r_test_g bias h2_post AUC p Status
QTL@MH binary MH BayesA 0.899 0.295 0.726 0.595 0.366 0.671 150 OK
QTL@MH binary MH BayesR 0.923 0.310 0.754 0.501 0.521 0.675 150 OK
QTL@MH binary MH GBLUP 0.744 0.128 0.573 1.466 0.182 0.568 200 OK
QTL@MH binary MH GWABLUP 0.744 0.128 0.572 1.460 0.182 0.571 200 OK
QTL@MH binary SNP BayesA 0.841 0.178 0.521 0.435 0.277 0.592 100 OK
QTL@MH binary SNP BayesR 0.786 0.204 0.550 0.645 0.177 0.611 100 OK
QTL@MH binary SNP GBLUP 0.686 0.164 0.497 1.943 0.182 0.600 200 OK
QTL@MH binary SNP GWABLUP 0.687 0.165 0.498 1.946 0.182 0.599 200 OK
QTL@MH continuous MH BayesA 0.784 0.444 0.801 1.071 0.202 NA 150 OK
QTL@MH continuous MH BayesR 0.786 0.461 0.833 0.865 0.289 NA 150 OK
QTL@MH continuous MH GBLUP 0.817 0.290 0.696 1.337 0.365 NA 200 OK
QTL@MH continuous MH GWABLUP 0.754 0.270 0.682 1.524 0.290 NA 200 OK
QTL@MH continuous SNP BayesA 0.751 0.261 0.530 0.737 0.186 NA 100 OK
QTL@MH continuous SNP BayesR 0.732 0.267 0.557 0.651 0.223 NA 100 OK
QTL@MH continuous SNP GBLUP 0.769 0.233 0.525 1.186 0.326 NA 200 OK
QTL@MH continuous SNP GWABLUP 0.776 0.232 0.524 1.084 0.345 NA 200 OK
QTL@SNP binary MH BayesA 0.875 0.338 0.561 0.688 0.347 0.678 150 OK
QTL@SNP binary MH BayesR 0.844 0.339 0.627 0.724 0.308 0.680 150 OK
QTL@SNP binary MH GBLUP 0.725 0.290 0.444 3.244 0.182 0.657 200 OK
QTL@SNP binary MH GWABLUP 0.724 0.289 0.447 3.204 0.182 0.657 200 OK
QTL@SNP binary SNP BayesA 0.866 0.439 0.797 0.892 0.383 0.752 100 OK
QTL@SNP binary SNP BayesR 0.818 0.454 0.894 0.921 0.332 0.763 100 OK
QTL@SNP binary SNP GBLUP 0.708 0.360 0.622 4.038 0.182 0.726 200 OK
QTL@SNP binary SNP GWABLUP 0.709 0.362 0.626 4.027 0.182 0.726 200 OK
QTL@SNP continuous MH BayesA 0.774 0.480 0.831 1.146 0.232 NA 150 OK
QTL@SNP continuous MH BayesR 0.713 0.531 0.939 0.937 0.303 NA 150 OK
QTL@SNP continuous MH GBLUP 0.762 0.392 0.601 2.696 0.265 NA 200 OK
QTL@SNP continuous MH GWABLUP 0.708 0.314 0.518 2.882 0.214 NA 200 OK
QTL@SNP continuous SNP BayesA 0.756 0.562 0.931 1.008 0.336 NA 100 OK
QTL@SNP continuous SNP BayesR 0.692 0.562 0.967 0.994 0.319 NA 100 OK
QTL@SNP continuous SNP GBLUP 0.770 0.499 0.766 2.303 0.336 NA 200 OK
QTL@SNP continuous SNP GWABLUP 0.769 0.531 0.817 2.182 0.337 NA 200 OK

Analysis

Simulation results support the theory of marker-QTL unit congruence. When the QTL is a combined effect of microhaplotype blocks, the accuracy of microhaplotype-based BayesR prediction reaches 0.832 for continuous traits and 0.758 for binary traits, which is 20-30% higher than that of SNPs. This phenomenon shows that when the QTL is a combined effect of multiple microhaplotype blocks, SNPs cannot capture the QTL signal, as it disperses across all blocks.

When the QTL is only a single independent SNP variant, the prediction accuracy of SNP-based BayesR reaches 0.973 for continuous traits and 0.897 for binary traits. Interestingly, although SNPs appear superior in this scenario, microhaplotypes remain quite competitive, with a smaller accuracy gap for continuous traits. This phenomenon shows that microhaplotypes can still capture the QTL signal, even though the QTL is composed of independent SNPs.

In both scenarios, the Bayesian models (BayesR and BayesA) outperform GBLUP and GWABLUP in bias calibration. GBLUP exhibits severe bias inflation, with values ​​exceeding 2.0 to 4.0 under QTL@SNP, as variance is spread uniformly across all markers rather than concentrated at the true QTL. BayesR has more calibrated predictions (bias approaching 1.0) than BayesA, due to its mixture priors that more aggressively downscale non-QTL markers toward zero. GWABLUP offers only marginal improvement over GBLUP, with the GWAS weighting step insufficient to recover the rare QTL signal captured by the full mixture prior.

Note

It should be noted, however, that the genotype data used in this analysis were generated through simulation and may not fully reflect the properties of empirical genotype data. In practice, the number of microhaplotype blocks that can be constructed depends on the density of SNPs within short physical windows (i.e., 150 bp), and regions containing only a single SNP within the defined window length cannot form a valid microhaplotype block. The effective marker density and genomic coverage of a microhaplotype panel may therefore differ considerably across species, sequencing platforms, and target enrichment designs, which could influence the magnitude of the advantages reported here.