HIBLUP 速查表 v2 · 详细版

• 从 https://www.hiblup.com/download 下载对应平台的可执行文件（学术免费，商用需联系开发团队），Linux 下先 chmod +x hiblup。
• 所有分析的命令结构都相同：./hiblup <功能开关> <输入文件> <模型参数> [过滤参数] --threads n --out 前缀。
• 每个功能由一个开关决定（如 --single-trait），HIBLUP 再根据你给的输入文件自动决定用哪个模型（见第 3 节）。
• 各文件中个体、SNP 的顺序不需要对齐，软件自动按 ID 匹配；表型文件中没有、但系谱/基因型中有的个体也会得到育种值。

1. 外部准备：用 PLINK/Beagle 等做基因型质控和填充，合并成一套 PLINK 二进制文件（HIBLUP 不做填充，只读一套 bfile）。
2. 群体结构检查：--pca、--cluster、--allele-freq --pop-class，决定是否需要混合 GRM 或把群体作为固定效应。
3. 构建关系矩阵一次：--make-xrm 生成 GA/GD 或 HA/HD，后续都用 --xrm 复用。
4. 单性状模型：逐性状 --single-trait，看 .anova 精简固定效应，看 .vars 得遗传力。
5. 多性状模型：--multi-trait 得遗传相关，必要时把单性状结果作为 --vc-priors 初值。
6. 附加输出：需要可靠性加 --r2，需要标记效应加 --snp-effect。
7. 应用：用 .rand 中的 GEBV 排名选种；用 .snpeff 对新个体 --pred 打分；用 --mating 做选配。
8. 验证：用交叉验证评估预测准确性（见第 9 节）。

过滤参数对所有涉及个体或 SNP 的功能都有效。

与其他软件不同，HIBLUP 的表型文件同时包含表型和所有环境因子（协变量、固定效应、环境随机效应），模型中的效应都通过“列号”引用。

id    sex  season  day  weight  loc  dam     tr1      tr2      tr3
Ind1  F    Winter  100  2.9     l3   Ind902  73.9776  21.3252  35.7737
Ind2  M    Summer  100  1.5     l2   Ind902  94.1994  19.1162  54.2453

• 必须有表头。
• 缺失值符号：NA Na . - NaN NAN nan na N/A n/a <NA>。
• 默认分析第 2 列，用 --pheno-pos 指定。
• 遗传随机效应不要通过 --rand 1 添加：提供 --pedigree、--bfile 或 --xrm 后自动拟合。

规则：性状内用逗号，性状间用空格，某性状没有该效应写 0，顺序与 --pheno-pos 一一对应。

--pheno-pos 8 9 10
--dcovar    2,3  0  2,3,6     # tr1: sex+season; tr2: 无; tr3: sex+season+loc
--qcovar    4,5  5  4
--rand      6,7  7  0

互作项用冒号：2:3 = sex×season（详见 4.6）。

id     father  mother
Ind2   NA      NA
Ind5   Ind1    NA
Ind17  Ind11   Ind5

• 三列依次为个体、父亲、母亲；未知亲本用缺失符号。
• 不要求按出生顺序排列，软件自动排序整理。
• 系谱中出现但未单独成行的亲本（如上例 Ind1）也会被处理。

• 只接受 PLINK 二进制 .bed/.bim/.fam，前缀传给 --bfile。
• 只能读一套文件：按染色体分开的需先 plink --merge-list 合并。
• 编码：A1A1=2，A1A2=1，A2A2=0，A1 为 .bim 第 5 列；所有 SNP 效应都以 A1 为效应等位基因。
• 不做填充：缺失基因型被当作杂合（编码 1），务必提前填充。
• 选配时 .fam 第 5 列须标性别（1=公，2=母）。

其他格式转换：

# VCF -> 二进制
plink2 --vcf demo.vcf --make-bed --out demo
# PED/MAP -> 二进制
plink2 --ped demo.ped --map demo.map --make-bed --out demo
# HapMap -> VCF（TASSEL，先排序再导出）
run_pipeline.pl -SortGenotypeFilePlugin -inputFile demo.hmp.txt -outputFile demo.sort.hmp.txt -fileType Hapmap
run_pipeline.pl -fork1 -h demo.sort.hmp.txt -export -exportType VCF -runfork1

HIBLUP 拟合的一般模型可写成：

y = Xb + Σ Wk·ck + Σ Zj·gj + e
    │      │          │       └ 残差 e
    │      │          └ 遗传随机效应（PA/PD、GA/GD、HA/HD 及其互作），来自 --pedigree / --bfile / --xrm
    │      └ 环境随机效应（--rand，如场地、母体、永久环境）
    └ 固定效应（--dcovar）与协变量（--qcovar），截距 mu 自动加入

HIBLUP 基于 V 矩阵（表型方差协方差矩阵）而不是 MME 稀疏左手项求解，因此：基因组数据（LHS 稠密）很高效；纯系谱或高度重复记录时反而不如 DMU、BLUPF90、ASReml。

• --xrm demo.GA 等价于 --bfile demo --add；--xrm demo.HA 即一步法。
• rrBLUP 与 GBLUP 理论等价：GBLUP 命令加 --snp-effect 即得 rrBLUP 的 SNP 效应。
• --pheno-pos 给 1 个位置为单性状，多个位置为多性状（--mme 据此自动切换）。

建议先建好、再用 --xrm 复用：逐性状分析时每次都重建矩阵很慢，复用结果完全相同。

(1) 系谱矩阵 A

# --add 加性 A（--dom 得显性）；--add-inv 同时输出 A 逆（给其他软件用）
./hiblup --make-xrm --pedigree demo.ped --add --add-inv --threads 32 --out demo
# -> demo.PA.bin  demo.PA.id  demo.PAinv.sparse.bin  demo.PAinv.sparse.id

(2) 基因组矩阵 G

# --add --dom 同时得 GA、GD；--code-method 编码方法（见下表）；--step 每块 SNP 数，控制内存
./hiblup --make-xrm --bfile demo --add --dom --code-method 1 --step 10000 --threads 32 --out demo
# -> demo.GA.bin  demo.GA.id  demo.GD.bin  demo.GD.id

G 的计算公式：G = ZZ' / [tr(ZZ')/n]，Z 为 n×m 的编码后基因型矩阵，不同方法只是 Z 的编码方式不同：

同一参数也影响基因组近交系数和亲缘系数（--ibc --rc --code-method）。

(3) G 的变体

(4) 一步法 H 矩阵：同时给系谱和基因型即自动构建。

./hiblup --make-xrm --pedigree demo.ped --bfile demo \
         --alpha 0.05 --add --dom --out demo
# -> demo.HA.*  demo.HD.*；用 --add-inv / --dom-inv 得 H 逆（*.sparse.*）

# 用自己算好的 GRM（须为二进制）构建 H
./hiblup --make-xrm --pedigree demo.ped --xrm yourGRM --alpha 0.05 --add --out demo

• 构建过程：先把 G 线性校正到与 A22 尺度一致（使 G 对角线均值、非对角线均值分别等于 A22 的），得到 Gadj；再 Gw = (1−α)·Gadj + α·A22，然后按 Aguilar 2010 / Christensen 2010 的公式组成 H。
• --alpha 越大，Gw 中系谱成分越多，可改善 G 的奇异性；示例使用 0.05。
• H 的构建同样支持 --snp-weight、--pop-class。

(5) 环境随机效应的矩阵

./hiblup --make-xrm --pheno demo.phe --rand 6,7 --out demo
# -> demo.loc.sparse.*  demo.dam.sparse.*

用途：把环境效应作为 --xrm 输入多性状模型，从而估计环境效应在性状间的相关（默认只估遗传和残差协方差）。

(6) 格式转换与跨软件

# 二进制 -> 文本方阵（demo.txt + demo.id.txt）；加 --triangle 输出下三角 i j k
./hiblup --trans-xrm --xrm demo.GA --out demo
# 文本 -> 二进制（ijk 格式需加 --triangle）
./hiblup --trans-xrm --xrm-txt demo.txt --xrm-id demo.id.txt --out demo
# GCTA/LDAK 的 GRM 转文本
./hiblup --trans-xrm --xrm gcta.grm --out hiblup
# HIBLUP 的稠密 GRM 给 GCTA 用
awk '{print $1,$1}' demo.GA.id > gcta.grm.id && cp demo.GA.bin gcta.grm.bin

• 文件名带 .sparse. 的矩阵不能直接给 GCTA/LDAK。
• 求逆报 “matrix is not invertible”：矩阵非正定，加 --ridge-value 小值。

(7) 基因型转码

./hiblup --trans-geno --bfile demo --add --out x             # x.geno.A.txt（0/1/2）
./hiblup --trans-geno --bfile demo --dom --out x             # x.geno.D.txt（杂合=1，纯合=0）
./hiblup --trans-geno --bfile demo --blupf90 --add --out x   # BLUPF90 格式（ID + 连续基因型字符串）

# 方式一：一步到位（每次都会重建矩阵）
./hiblup --single-trait \
         --pheno demo.phe --pheno-pos 8 \
         --dcovar 2,3 --qcovar 4,5 --rand 6,7 \
         --pedigree demo.ped --bfile demo \
         --add --dom --threads 32 --out demo

# 方式二（推荐）：先建矩阵，再复用
./hiblup --make-xrm --pedigree demo.ped --bfile demo --add --dom --out demo
./hiblup --single-trait \
         --pheno demo.phe --pheno-pos 8 \
         --dcovar 2,3 --qcovar 4,5 --rand 6,7 \
         --xrm demo.HA,demo.HD --add --dom \
         --vc-method AI --ai-maxit 30 \
         --threads 32 --out demo

• 构建 XRM 支持的选项（加权 G、混合 G、H 矩阵等）在这条命令里都能直接用。
• 输出：demo.vars、demo.beta、demo.anova、demo.rand（解读见第 6 节）。
• 多个性状逐个分析时，只改 --pheno-pos 和 --out，矩阵复用即可。

# 表型文件中同一 ID 可出现多行（多次记录）
./hiblup --single-trait \
         --pheno demo.repeat.phe --pheno-pos 2 \
         --rand 1 \
         --bfile demo --threads 32 --out demo

• --rand 1：把第 1 列个体 ID 作为环境随机效应，即永久环境效应；遗传效应仍由 --bfile 拟合。
• 结果中 ID 行（永久环境）方差 + 遗传方差，占表型方差的比例即重复力 r = (Va + Vpe) / Vp。
• 性能提示：重复次数少（如每个体 2 次）HIBLUP 更快；重复次数多（如每个体 10 次）建议 BLUPF90、ASReml 等 MME 类软件。

./hiblup --multi-trait \
         --pheno demo.phe --pheno-pos 8 9 10 \
         --qcovar 4,5 5 4 \
         --dcovar 2,3 0 2 \
         --rand 6,7 7 0 \
         --xrm demo.GA,demo.GD \
         --vc-method AI --ai-maxit 30 \
         --threads 32 --out demo

本例每个性状的模型：

• 可用于一个群体内多个性状，或把不同群体的同一性状当作不同性状估跨群体遗传相关。
• 支持不平衡数据（v1.5.0 及以上），性状数量不限。
• 默认估计残差协方差；如需不估残差协方差，加 --ignore-cove。
• 只估遗传效应和残差的协方差，至少要有一个遗传随机效应；环境效应相关见 4.1(5)。
• 输出：demo.vars、demo.covars，以及每个性状的 demo.<trait>.beta/.anova/.rand。
• 先用单性状结果作为 --vc-priors、--covc-priors 初值，可明显加快多性状收敛。

适用于方差组分已知（文献值或已估计）、只需快速计算育种值的场景，如例行遗传评估。

# 单性状：loc 方差、GA 方差、残差方差
./hiblup --mme --pheno demo.phe --pheno-pos 8 --rand 6 \
         --xrm demo.GA --vc-priors v1,v2,v3 --pcg --threads 32 --out demo

# 多性状：加性 + 显性
./hiblup --mme --pheno demo.phe --pheno-pos 8 9 10 --rand 6,7 7 0 \
         --bfile demo --add --dom \
         --vc-priors t1_v1,t1_v2,t1_v3,t1_v4,t1_v5 t2_v1,t2_v2,t2_v3,t2_v4 t3_v1,t3_v2,t3_v3 \
         --covc-priors cov1,cov2,cov3 cov4,cov5,cov6 cov7,cov8,cov9 \
         --pcg --threads 32 --out demo

--vc-priors 的顺序：每个性状内依次为 环境随机效应（按 --rand 顺序）→ 遗传效应（按 --xrm 或 add、dom 顺序）→ 残差；性状内逗号、性状间空格。

--covc-priors 的顺序：按“遗传效应 1、遗传效应 2 … 残差”分块（块数 = 遗传矩阵数 + 1），块内为性状间协方差矩阵的下三角，按列依次填写，块内逗号、块间空格。3 个性状时每块为 cov(1,2), cov(1,3), cov(2,3)：

• 默认直接求 V 矩阵的逆；--pcg 用预条件共轭梯度迭代求解，很多情况下更快。
• --vc-priors、--covc-priors 在 --single-trait/--multi-trait 中同样可用，作为 REML 起始值。

• G×G 只能通过 --xrm 指定，矩阵为各矩阵的 Hadamard 积再按迹标准化：G = (G1∘G2∘…)/[tr(G1∘G2∘…)/N]。
• 多性状时同样“性状间空格、无则写 0”：

--dcovar 2:3 0 2,3,2:3 --rand 6 7,6:7 0
--rand-gxe 6:demo.GA 6:demo.GA,7:demo.GA 0

(1) 计算 SNP 效应

# 方式一：拟合模型时直接输出（此例为一步法单性状，加性）
./hiblup --single-trait --pheno demo.phe --pheno-pos 8 \
         --bfile demo --pedigree demo.ped --add --snp-effect --threads 32 --out demo

# 方式二：由已有 GEBV 反推（会内部构建 G）
./hiblup --snp-effect --gebv demo.gebv.d.txt --bfile demo --dom --threads 32 --out demo

# 方式二加速：提供已建好的 G
./hiblup --snp-effect --gebv demo.gebv.ad.txt --bfile demo \
         --xrm demo.GA,demo.GD --add --dom --threads 32 --out demo

• 结果写入 demo.snpeff。
• GEBV 若来自加权 G，反推时必须给同样的 --snp-weight。
• 应用：GWAS 式查看大效应位点；作为权重构建加权 G；给新个体打分。

(2) 预测 GEBV / PRS

./hiblup --pred --bfile new_animals --score demo.snpeff --threads 10 --out new
# 同时含加性、显性效应时：加 --add --dom（score 文件须为 6 列）
# 只有显性效应：加 --dom（按 0/1/0 编码）

• 默认按加性编码（0/1/2）打分；功能同 PLINK --score，但速度快数倍。
• 输出 new.bv，每个效应列对应一个预测值列。
• 新个体的 .bim 与 .snpeff 中同名 SNP 的 A1/A2 必须一致。

(3) 基因组选配

./hiblup --mating --bfile demo --score demo.snpeff --threads 32 --out demo

• 对所有公母组合预测后代期望遗传表现，同时利用加性和显性效应（适合杂交生产优化）。
• .fam 第 5 列必须标性别：1=公、2=母，其他视为未知并被剔除。
• .bim 与 .snpeff 的等位基因方向必须一致。

# 等位基因 / 基因型频率（可按群体分列）
./hiblup --bfile demo --allele-freq --pop-class demo.popclass.txt --out demo   # demo.afreq
./hiblup --bfile demo --geno-freq --out demo                                   # demo.gfreq
# 按频率筛选 SNP，例如 0.05 < freq <= 0.1
awk '$4 > 0.05 && $4 <= 0.1 {print $1}' demo.afreq > freq_filter_snp.txt

# 纯合度 / 杂合度
./hiblup --bfile demo --homo --hete --out demo        # demo.homo demo.hete

# 近交系数与亲缘系数（系谱、基因型或已有矩阵均可）
./hiblup --ibc --rc --pedigree demo.ped --out demo
./hiblup --ibc --rc --bfile demo --code-method 1 --out demo
./hiblup --ibc --rc --xrm demo.GA --out demo

# PCA（--npc 输出前 n 个 PC，--step 控制内存）
./hiblup --pca --bfile demo --npc 5 --out demo
./hiblup --pca --xrm demo.GA --npc 5 --out demo

# 无监督聚类（三种输入任选其一）
./hiblup --cluster --bfile demo --npc 10 --threads 16 --out demo
./hiblup --cluster --xrm demo.GA --npc 10 --threads 16 --out demo
./hiblup --cluster --pc demo.pc.txt --npc 5 --threads 16 --out demo
# 关闭 UMAP，直接在 PC 空间做 HDBSCAN（自动选 DBCV 最高的参数）
./hiblup --cluster --bfile demo --umap-off --npc 10 --out demo

• 聚类流程：PCA 提取主要遗传变异 → UMAP 降到二维揭示非线性结构 → 改进 HDBSCAN 按局部密度识别稳定亚群，不需要预设群体数。
• --npc 必须小于分析个体数。
• 用 --pc 直接输入 PC 时没有特征值加权，结果可能与用基因型/XRM 输入不同。

# 1) 参考群 LD 矩阵（输出 .bin + .info，--write-txt 可读）
./hiblup --ld --bfile ref --window-bp 1e6 --threads 16 --out ref_ldm

# 2) LD score
./hiblup --ldscore --bfile ref --window-bp 1000000 --threads 10 --out ref     # ref.ldsc

# 3) LD score 回归：一个 sumstat 估 h2；多个 sumstat 同时估 h2 与两两 rG
./hiblup --ldreg --sumstat t1.ma --lds ref.ldsc --out x
./hiblup --ldreg --sumstat t1.ma t2.ma t3.ma --lds ref.ldsc --out x

# 4) SBLUP：用基因型或预计算 LD 矩阵
./hiblup --sblup --sumstat t1.ma --bfile ref --window-bp 1e6 --h2 0.3234 --threads 16 --out x
./hiblup --sblup --sumstat t1.ma --ldm ref_ldm --h2 0.3234 --threads 10 --out x

# 5) 多性状 SBLUP：h2 每性状一个，rg 按下三角 1-2 1-3 1-4 2-3 2-4 3-4
./hiblup --sblup --sumstat t1.ma t2.ma t3.ma t4.ma --ldm ref_ldm \
         --h2 0.3234 0.1256 0.6345 0.3536 \
         --rg 0.1336 0.5567 0.2345 0.8454 0.3446 0.4633 --threads 10 --out x

# 6) 用估得的 SNP 效应给个体打分
./hiblup --pred --bfile target --score x.snpeff --out target

• --h2、--rg 可来自 REML（有个体数据时）或 LD score 回归。
• 用 --ldm 做 SBLUP 要求 SNP 满足哈迪-温伯格平衡，否则效应有偏、预测变差。
• 窗口内 SNP 超过约 1 万时加 --pcg。
• 不要删减 LD score 文件中的 SNP 去对齐 sumstat：总 SNP 数对回归很关键。
• 使用公开 LD score 且无 maf 列：若指定了 --M，maf 列可填任意值。

HE 细节：

• 有固定效应时先把表型对其回归，再对残差做 HE 回归。
• 默认不输出随机效应解，需要加 --he-pred；此时可再加 --pcg，即 “HE+PCG” 策略，用于超大数据。
• 样本量不足时 HE 估计可能有偏、不稳定，中小数据坚持用 REML。

以下为通用经验，非官方规定：

1. 看 LOG 中每次迭代的 loglike 是否单调上升并趋于稳定；震荡或达到上限说明未收敛。
2. 检查表型：极端值、单位、缺失编码是否在 HIBLUP 认可的列表内。
3. 检查固定效应是否有水平记录过少或与随机效应混杂（.anova 中 Df 异常）。
4. 换 --vc-method EMAI 或 EM，或用单性状结果作为 --vc-priors。
5. 多性状先两两拟合，再合并；去掉估计接近 0 的随机效应。
6. 某方差组分趋于 0 或 SE 远大于估计值，说明数据不支持该效应，考虑从模型中删除并用 LRT 验证（第 7 节）。

单性状示例（一步法加性 + 显性，loc、dam 为环境随机效应）：

      Var      Var_SE   h2      h2_SE
loc   12.2127  8.9109   0.1035  0.0682
dam   10.6268  4.8605   0.0901  0.041
HA    59.2007  11.7641  0.5019  0.0807
HD    28.9946  8.9294   0.2458  0.083
e     6.9232   1.6957   0.0587  0.0161

怎么算、怎么读：

• 表型方差 Vp = 各行 Var 之和 = 12.21 + 10.63 + 59.20 + 28.99 + 6.92 ≈ 117.96。
• 狭义遗传力 h² = V(HA) / Vp = 59.20 / 117.96 ≈ 0.502，即 HA 行的 h2。
• 显性方差占比 d² = V(HD) / Vp ≈ 0.246；广义遗传力 H² ≈ (59.20 + 28.99) / 117.96 ≈ 0.75。
• loc 行 h2 = 0.10 表示场地效应解释 10% 的表型方差；dam 行即母体（共同环境）效应比例 c²。
• 所有行的 h2 之和约为 1。

经验判读：

多性状时行名加性状前缀，如 tr1_demo.GA、tr2_e，每个性状单独计算 h2。

Item             COVar    COVar_SE  r          r_SE
tr1:tr2_demo.GA  46.6899  8.56852   0.747524   0.171364
tr1:tr3_demo.GA  16.3376  5.96333   0.454749   0.185599
tr1:tr2_demo.GD  26.9273  4.65452   0.796922   0
tr1:tr2_e        0.24039  1.17428   0.024379   0.119333

• 行名 tr1:tr2_demo.GA：性状 1 与 2 在加性效应上的协方差和相关；_e 行为残差协方差/相关。
• 相关 r = COVar / √(Var₁·Var₂)。例：46.69 / √(61.35 × 63.59) ≈ 0.748（Var 取自 .vars 中 tr1_demo.GA、tr2_demo.GA）。
• 加性遗传相关 rA（GA/HA/PA 行）用于判断间接选择和制定选择指数：rA 高（如 0.75）说明两性状受相似的加性基因控制，选一个会带动另一个。
• 显著性快速判断：r ± 2·r_SE 不跨 0 即可认为显著；tr1:tr2_e 的 0.024 ± 0.24 跨 0，残差相关不显著。
• r_SE = 0（如 tr1:tr2_demo.GD）通常说明该相关估计触到参数边界或 SE 未能计算，应谨慎对待该值。
• 表型相关可由 rP ≈ (ΣCOVar) / √(Vp₁·Vp₂) 计算。

Levels         Estimation  SE
mu             32.5285     3.4758
sex_F          6.2933      1.7593
sex_M          26.2353     1.7634
season_Autumn  18.0905     0.9753
day            0.1547      0.0574
bornweight     1.5703      0.4614

• mu：总体截距。
• 因子_水平（如 sex_F、season_Autumn）：离散固定效应各水平的估计；比较同一因子不同水平之间的差值最有意义，例：公母差约 26.24 − 6.29 ≈ 19.9。
• 数量协变量（day、bornweight）：回归系数，表示协变量增加 1 个单位时表型的变化量，例：出生重每增加 1 单位，性状增加约 1.57。
• 多性状时为 demo.<trait>.beta。

Factors     Df   SumSq       MeanSq      F         Pr(>F)
sex         1    49710.3845  49710.3845  7180.311  <2e-16
season      3    24515.0941  8171.698    1180.344  <2e-16
day         1    925.3352    925.3352    133.658   <2e-16
bornweight  1    437.3432    437.3432    63.171    1.30E-14
e           493  3413.1141   6.9232

• Df：离散因子 = 水平数 − 1（season 4 个季节 → 3），协变量 = 1。
• Pr(>F)：F 检验 p 值，p < 0.05 说明该效应显著，应保留；不显著的可考虑从模型中删除后重跑。
• e 行 MeanSq 等于残差方差（6.9232，与 .vars 中 e 一致）。

• 每行一个个体，列为每个随机效应的解：
• PA / GA / HA 列：加性育种值（EBV/GEBV），包括没有表型、只有系谱或基因型的个体——这正是候选群体选种用的值。
• PD / GD / HD 列：显性离差，不遗传给后代，选种通常只看加性列；杂交组合优化时才用。
• 环境效应列（loc、dam 等）：该个体所在水平的效应估计。
• residuals：残差。
• 多性状为 demo.<trait>.rand；多性状育种值可按经济权重组合成选择指数。
• 育种值是相对值（以群体为基准的离差），只在同一次分析内比较排名有意义。

加 --r2 后每个效应多出 3 列：

ID       GA        GA_SE    GA_PEV   GA_R2
IND1001  6.66127   4.27597  18.2839  0.774767
IND1003  2.08201   4.62091  21.3528  0.736963

• 准确性 = √R²，例 √0.775 ≈ 0.88。
• R² 越接近 1 越可靠：本身有表型、或有很多有表型近亲的个体高；无表型、亲缘远的候选个体低。
• 实用：两个个体 GEBV 之差小于约 2 倍 SE 时，排名差异不确定。
• --r2 计算往往比 REML 本身还慢，只在需要时开启。

• 先核对 LOG 中打印的模型方程，确认每个效应被识别为预期的类型（例如字符型协变量误放进 --qcovar 会出错）。
• 看迭代过程中 loglike 是否收敛；最后一次迭代的 loglike 用于似然比检验。

# x.ldsr.h2
Item  Intercept  Intercept_SE  h2        h2_SE       h2_Pval
demo  1.08285    0.011433      0.122826  0.00393122  2.71554e-214

# x.ldsr.rg
Item         CovG       CovG_SE     Intercept  Intercept_SE  rG        rG_SE      rG_Pval
demo1:demo2  0.0252414  0.00716841  0.0166108  0.00819384    0.141956  0.0403148  0.000429607

自动完成：看 .anova 最后一列 Pr(>F)（见 6.4）。

HIBLUP 不自动检验随机效应，需手动做：

1. 跑完整模型（含待检验效应，如 HD），在 LOG 最后一次 REML 迭代处记下 loglike = L1。
2. 跑简化模型（去掉该效应，其他完全相同），记下 L0。
3. 在 R 中计算：

stat   <- 2 * (L1 - L0)
pvalue <- pchisq(stat, df = 1, lower.tail = FALSE)
pvalue

• 一次去掉一个方差组分时 df = 1。
• 补充（统计学通用做法，非教程内容）：方差组分在参数边界 0 处检验时，上述 p 值偏保守，常用 p/2（0.5χ²₀ + 0.5χ²₁ 混合分布）。
• 示例判读：stat > 3.84 即 p < 0.05，该效应显著，应保留。

# 1) 显著 SNP 单独建 GRM（snp.txt：一列 SNP ID，无表头）
./hiblup --make-xrm --bfile demo --extract snp.txt --threads 32 --out grm1   # grm1.GA.*
# 2) 其余 SNP 建 GRM
./hiblup --make-xrm --bfile demo --exclude snp.txt --threads 32 --out grm2   # grm2.GA.*
# 3) 两个 GRM 同时作为随机效应
./hiblup --single-trait --pheno demo.phe --pheno-pos X --xrm grm1.GA,grm2.GA --threads 32 --out vc

• 计算：PVE = V(grm1) / [V(grm1) + V(grm2) + V(e)]，三个方差都在 vc.vars 中。
• 也就是 vc.vars 中 grm1.GA 行的 h2；grm2.GA 行的 h2 则为其余多基因背景的贡献。

• 基因型缺失会被当作杂合：HIBLUP 不做填充，先用 Beagle 等填充。
• 只读一套 bfile：分染色体的文件先合并。
• 表型第 1 列必须是 ID，且必须有表头；字符型变量不能放进 --qcovar。
• 多性状参数个数必须与性状数一致：没有的效应写 0，否则会错位。
• 等位基因方向：所有 SNP 效应针对 .bim 中的 A1；跨数据集打分、选配前核对 A1/A2。
• 选配：.fam 第 5 列性别未标注的个体会被剔除。

• --xrm 只接受二进制 .bin + .id；文本矩阵用 --trans-xrm 转换（Q5）。
• 不要把逆矩阵传给 --xrm：HIBLUP 直接用关系矩阵本身。
• 矩阵无法求逆：加 --ridge-value；一步法可适当调整 --alpha。
• 带 .sparse. 的矩阵不能直接给 GCTA/LDAK。

• 内存不够（Q2）：调小 --step；或加 --float-prec（内存约减半，个别情况损失精度）。
• 逐性状重复建矩阵很慢：先 --make-xrm，再 --xrm 复用。
• 超大数据（>10 万个体）：--vc-method HE --he-pred --pcg；或先在子集上用 REML 估参数，再 --mme --pcg 求解全数据。
• 纯系谱、高度重复记录：改用 DMU / BLUPF90 / ASReml。

HIBLUP 教程没有内置交叉验证命令，下面是基于 HIBLUP 的常用做法（通用方法，供参考）：

1. 把有表型个体随机分为 k 份（常用 5 折，可重复多次）。
2. 每次将一份个体的表型设为 NA（作为验证集），其余为参考群；基因型/系谱保持完整。
3. 用相同模型跑 --single-trait（或参数已知时 --mme），从 .rand 取验证集个体的 GEBV。
4. 计算指标：

• “校正表型”指表型减去固定效应（可用 .beta 计算），或直接用原始表型的近似。
• 比较不同模型（PBLUP vs GBLUP vs SSGBLUP、有无显性、加权与否）时，使用同一组折分。
• 你的环境中有 mintyr-pipelines 技能，可按性状 × 组别做 k 折划分并直接导出 HIBLUP 输入文件，需要时可以让我用它批量准备。