ASReml-R 4.2 & ASRgenomics 速查表

library(asreml)
asreml.options(
  threads    = 8,       # 并行线程 (v4.2)
  workspace  = "8gb",   # REML 求解内存
  pworkspace = "4gb",   # predict() 内存
  maxit      = 30,      # 默认 13
  ai.sing    = TRUE     # AI 矩阵奇异时继续迭代
)
asreml.options()        # 查看当前设置

# 1. 不拟合，仅导出参数表
sv <- asreml(yield ~ Variety, random = ~ Rep,
             data = nin89, start.values = TRUE)
iv <- sv$vparameters.table
print(iv)   # Component, Value, Constraint

# 2. 按名称修改
iv$Value[iv$Component == "Rep"] <- 0.5

# 3. 传回重新拟合
mod <- asreml(yield ~ Variety, random = ~ Rep,
              data = nin89, G.param = iv, R.param = iv)

# 场景 A: 已知参数，固定方差直接解 BLUP —— 按名称赋值
sv <- asreml(DBH_Adj ~ 1, random = ~ vm(Genotype, ginv),
             residual = ~ idv(units), data = dat1,
             start.values = TRUE)
iv <- sv$vparameters.table
print(iv$Component)   # 先确认组分名称
iv$Value[iv$Component == "vm(Genotype, ginv)"] <- 1.1427
iv$Value[iv$Component == "units!units"]        <- 1.3724
iv$Constraint <- "F"
mod_blup <- asreml(DBH_Adj ~ 1, random = ~ vm(Genotype, ginv),
                   residual = ~ idv(units), data = dat1,
                   G.param = iv, R.param = iv, maxit = 1)

# 场景 B: 两阶段 MET，固定残差尺度 = 1 (wt = 1/SE^2)
sv_wt <- asreml(blue ~ Site, random = ~ fa(Site, 2):Variety,
                weights = wt, data = stage1_blues,
                start.values = TRUE)
iv_wt <- sv_wt$vparameters.table
r_row <- iv_wt$Component == "units!R"
iv_wt$Value[r_row] <- 1; iv_wt$Constraint[r_row] <- "F"
mod_stage2 <- asreml(blue ~ Site, random = ~ fa(Site, 2):Variety,
                     weights = wt, data = stage1_blues,
                     G.param = iv_wt, R.param = iv_wt)

# 场景 C: 约束父本方差 = 母本方差 (NCII)
mod_unc <- asreml(length ~ rep,
                  random = ~ female + male + female:male,
                  data = pine)                       # 无约束模型
sv_nc <- asreml(length ~ rep,
                random = ~ female + male + female:male,
                data = pine, start.values = TRUE)
gam <- sv_nc$vparameters.table
# 第 1 列: 约束组号(相同 = 相等); 第 2 列: 比例系数
M <- cbind(c(1, 1, 2, 3), c(1, 1, 1, 1))
dimnames(M) <- list(gam$Component, NULL)  # 行数须等于组分数
mod_eq <- asreml(length ~ rep,
                 random = ~ female + male + female:male,
                 vcc = M, data = pine)
lrt.asreml(mod_unc, mod_eq, boundary = FALSE)   # H0: σ²m = σ²f

ainv <- ainverse(ped)   # Animal, Sire, Dam

# 模型 1: 不用系谱，个体作普通随机因子 -> 总重复力
m_rep <- asreml(date ~ 1 + Age, random = ~ Animal,
                residual = ~ idv(units), data = dat)
vpredict(m_rep, rep ~ V1 / (V1 + V2))

# 模型 2: 基础动物模型 (仅适用于单次测量)
m_her <- asreml(date ~ 1 + Age, random = ~ vm(Animal, ainv),
                residual = ~ idv(units), data = dat)

# 模型 3: 重复力动物模型 (加性 + 永久环境)
m_full <- asreml(date ~ 1 + Age,
                 random = ~ vm(Animal, ainv) + ide(Animal),
                 residual = ~ idv(units), data = dat)
summary(m_full)$varcomp  # V1 vm, V2 ide, V3 units!units, V4 units!R

vpredict(m_full, h2  ~ V1 / (V1 + V2 + V3))
vpredict(m_full, pe2 ~ V2 / (V1 + V2 + V3))
vpredict(m_full, rep ~ (V1 + V2) / (V1 + V2 + V3))

lrt.asreml(m_full, m_her, boundary = TRUE)  # PE 是否显著

# M1: Dam 作普通随机因子 (母体总效应, 不分遗传/环境)
m1 <- asreml(bwt ~ 1 + sex, random = ~ vm(Animal, ainv) + Dam,
             residual = ~ idv(units), data = gryphon)

# M2: 直接加性 + 母体遗传, Cov(a,m) = 0
m2 <- asreml(bwt ~ 1 + sex,
             random = ~ vm(Animal, ainv) + vm(Dam, ainv),
             residual = ~ idv(units), data = gryphon)

# M3: + 母体永久环境 ide(Dam)
m3 <- asreml(bwt ~ 1 + sex,
             random = ~ vm(Animal, ainv) + vm(Dam, ainv) + ide(Dam),
             residual = ~ idv(units), data = gryphon)

# M4: 允许 Cov(a,m) != 0
m4 <- asreml(bwt ~ 1 + sex,
  random = ~ str(~ vm(Animal, ainv) + vm(Dam, ainv),
                 ~ us(2, init = c(1, -0.1, 0.5)):vm(Animal, ainv)) +
             ide(Dam),
  residual = ~ idv(units), data = gryphon)

# V1 = Va, V2 = Cov(a,m), V3 = Vm, V4 = Vc, V5 = Ve
summary(m4)$varcomp
vpredict(m4, h2_d ~ V1 / (V1+V2+V3+V4+V5))
vpredict(m4, h2_m ~ V3 / (V1+V2+V3+V4+V5))
vpredict(m4, c2   ~ V4 / (V1+V2+V3+V4+V5))
vpredict(m4, r_am ~ V2 / sqrt(V1 * V3))
vpredict(m4, h2_T ~ (V1 + 0.5*V3 + 1.5*V2) / (V1+V2+V3+V4+V5))

init_G <- c(3.0, 0.5, 2.0)   # V1, Cov12, V2
init_R <- c(4.0, 0.2, 2.5)
biv <- asreml(
  fixed    = cbind(bwt, tarsus) ~ trait + trait:sex,
  random   = ~ us(trait, init = init_G):vm(Animal, ainv) +
               us(trait):Dam,
  residual = ~ id(units):us(trait, init = init_R),
  na.action = na.method(y = "include", x = "include"),
  data = gryphon, maxit = 30)

init_G3 <- c(10, 2, 15, 1.5, 3, 20)  # V1,C12,V2,C13,C23,V3
triv <- asreml(
  fixed    = cbind(y1, y2, y3) ~ trait + trait:Line,
  random   = ~ us(trait, init = init_G3):vm(Calf, ainv),
  residual = ~ id(units):us(trait),
  na.action = na.method(y = "include"),
  data = harvey, maxit = 50)

# V1 Va1, V2 CovA, V3 Va2, V4 units!R, V5 Ve1, V6 CovE, V7 Ve2
summary(biv)$varcomp
vpredict(biv, h2_1 ~ V1 / (V1 + V5))
vpredict(biv, h2_2 ~ V3 / (V3 + V7))
vpredict(biv, r_g  ~ V2 / sqrt(V1 * V3))
vpredict(biv, r_e  ~ V6 / sqrt(V5 * V7))
vpredict(biv, r_p  ~ (V2 + V6) / sqrt((V1 + V5) * (V3 + V7)))

# corgh: 直接估计 r_g 及其 SE
biv_cor <- asreml(cbind(t1, t2) ~ trait + trait:Block/Rep,
                  random = ~ corgh(trait):vm(ID, ainv),
                  residual = ~ id(units):us(trait), data = nested)

# 约束模型 Cov_A = 0
biv0 <- asreml(cbind(t1, t2) ~ trait + trait:Block/Rep,
               random = ~ idh(trait):vm(ID, ainv),
               residual = ~ id(units):us(trait), data = nested)
lrt.asreml(biv, biv0, boundary = FALSE)

library(ASRgenomics); library(asreml); library(dplyr)

data(pheno.pine); dat <- pheno.pine
dat[, 1:3] <- lapply(dat[, 1:3], as.factor)

data(geno.pine655)
M_clean <- qc.filtering(M = geno.pine655, maf = 0.05,
             marker.callrate = 0.2, ind.callrate = 0.2,
             na.string = "-9", plots = FALSE)$M.clean
Gmat <- G.matrix(M_clean, method = "VanRaden")$G

# 表型与 G 对齐
chk  <- match.kinship2pheno(K = Gmat, pheno.data = dat,
          indiv = "Genotype", clean = FALSE, mism = TRUE)
dat1 <- droplevels(dat[chk$matchesP, ])

# 稀疏逆矩阵
G_b  <- G.tuneup(G = Gmat, blend = TRUE, pblend = 0.02)$Gb
ginv <- G.inverse(G = G_b, sparseform = TRUE)$Ginv.sparse

modg <- asreml(DBH_Adj ~ 1, random = ~ vm(Genotype, ginv),
               residual = ~ idv(units),
               na.action = na.method(y = "include"), data = dat1)
summary(modg)$varcomp
vpredict(modg, h2 ~ V1 / (V1 + V2))

# 导出 GEBV
blup <- as.data.frame(summary(modg, coef = TRUE)$coef.random)
gebv <- data.frame(
  ID   = sub("^vm\\(Genotype, ginv\\)_", "", rownames(blup)),
  gebv = blup$solution, se = blup$std.error)
write.csv(gebv, "gebv_gblup.csv", row.names = FALSE)

library(ASRgenomics); library(AGHmatrix); library(asreml)

data(ped.pine)
A <- AGHmatrix::Amatrix(data = ped.pine)   # n_a x n_a

G <- G.matrix(M_clean, method = "VanRaden")$G

# 1. 检查并匹配 G 与 A
chk <- match.G2A(A = A, G = G, clean = TRUE, ord = TRUE,
                 mism = TRUE, RMdiff = TRUE)

# 2. 将 G 对齐到 A22，再求完整逆
G_al <- G.tuneup(G = chk$Gclean, A = chk$Aclean,
                 align = TRUE, sparseform = FALSE)$Gb
Ginv <- G.inverse(G = G_al, sparseform = FALSE)$Ginv

# 3. 构建 H^-1 (λ 同时设定 τ = ω，范围 0–1)
hinv <- H.inverse(A = A, Ginv = Ginv, lambda = 0.90,
                  sparseform = TRUE)
attr(hinv, "INVERSE")   # 应为 TRUE

# 4. 拟合：已测序与未测序个体同时得到 GEBV
ss <- asreml(DBH_Adj ~ 1, random = ~ vm(Genotype, hinv),
             residual = ~ idv(units), workspace = "2gb",
             na.action = na.method(y = "include"), data = dat)
summary(ss)$varcomp
vpredict(ss, h2 ~ V1 / (V1 + V2))
summary(ss)$aic
blup_ss <- summary(ss, coef = TRUE)$coef.random

set.seed(2026); k <- 5
dat1$fold <- sample(rep(1:k, length.out = nrow(dat1)))
pa <- b1 <- numeric(k)
term <- "vm(Genotype, ginv)_"

for (i in 1:k) {
  dcv <- dat1
  val <- which(dcv$fold == i)
  dcv$y_train <- dcv$DBH_Adj
  dcv$y_train[val] <- NA                 # 掩盖验证集

  m <- asreml(y_train ~ 1, random = ~ vm(Genotype, ginv),
              residual = ~ idv(units),
              na.action = na.method(y = "include"),
              data = dcv, trace = FALSE)

  cr   <- summary(m, coef = TRUE)$coef.random
  gebv <- cr[paste0(term, dcv$Genotype[val]), "solution"]
  yv   <- dat1$DBH_Adj[val]
  pa[i] <- cor(yv, gebv, use = "complete.obs")
  b1[i] <- coef(lm(yv ~ gebv))[2]
}
cat("PA:", mean(pa), "| Acc:", mean(pa) / sqrt(h2_full),
    "| b1:", mean(b1), "\n")

# G.predict: G 为参考群 + 候选群的完整方阵
gy <- as.matrix(gebv_ref)       # 一列 BLUP
rownames(gy) <- ref_ids
g_cond <- G.predict(G = G_full, gy = gy)
head(g_cond)                     # predicted.value

library(asreml); library(ASRgwas)

pre_c <- pre.gwas(
  pheno.data = dat, indiv = "ID", resp = "amoebic_load",
  geno.data = geno.salmon, Q.method = "K", method = "VanRaden",
  maf = 0.05, marker.callrate = 0.2, ind.callrate = 0.2,
  heterozygosity = 0.8, Fis = 0.98, impute = TRUE)
pre_c$plot.scree               # 确定 npc

gw <- gwas.asreml(
  pheno.data = pre_c$pheno.data, resp = "amoebic_load", gen = "ID",
  Kinv = pre_c$Kinv, Q = pre_c$Q, npc = 7,
  geno.data = pre_c$geno.data, map.data = pre_c$map,
  pvalue.thr = 5e-04, workspace = "2gb", P3D = TRUE)

gw$heritability                # 零模型 h2
gw$gwas.sel                    # 超过阈值的 SNP
qq.plot(gwas.table = gw$gwas.all)
manhattan.plot(gwas.table = gw$gwas.all, padding = 2e6)

pre_b <- pre.gwas(
  pheno.data = dat, indiv = "ID", resp = "ipn_mortality",
  geno.data = geno.salmon, Q.method = "K", method = "VanRaden",
  maf = 0.05, marker.callrate = 0.2, ind.callrate = 0.2,
  heterozygosity = 0.8, Fis = 0.98, impute = TRUE)

gwB <- gwas.asreml(
  pheno.data = pre_b$pheno.data, resp = "ipn_mortality", gen = "ID",
  Kinv = pre_b$Kinv, Q = pre_b$Q, npc = 3,
  cov = "tag_wgt", fixedf = "sex", randomf = "family",
  family = "binomial",
  geno.data = pre_b$geno.data, map.data = pre_b$map,
  pvalue.thr = 5e-04, threads = 1, P3D = FALSE)

gwB$call                       # 内部 asreml 公式
summary(gwB$mod)$varcomp