狭义遗传力（）：仅考虑加性遗传效应（SNP的线性效应），反映可遗传的变异占表型总变异的比例。
GBLUP模型-rrBLUP包
library(rrBLUP)
#读取样本名
sample.names <- read.table("/root/home/liuqirui/SNPdata/sample_names.txt",header=F)
#读取012格式的SNP数据
geno <- fread("/root/home/liuqirui/SNPdata/260SNP.imputed.txt", header=T)
geno <- geno[,-1]
geno <- sapply(geno, as.numeric)
#矩阵-1转为-1,0,1格式
geno <- geno-1
rownames(geno) <- sample.names[,1]
#使用A.mat()函数计算K矩阵
K <- A.mat(geno)
#读取表型数据
pheno1 <- as.vector(unlist(read.table("/root/home/liuqirui/SNPdata/260sample_flowerNUM.txt",header=F,row.names=1)))
names(pheno1) <- sample.names[,1]
#转为数据框
pheno1_df <- data.frame(ID = rownames(geno), Phenotype = pheno1)
#使用kin.blup()进行计算，这一步为核心步骤
gblup1_fit <- kin.blup(data = pheno1_df, geno = "ID", pheno = "Phenotype", K = K)
#从结果获取遗传方差和环境方差
var_g_1 <- gblup1_fit$Vg
var_e_1 <- gblup1_fit$Ve
#利用狭义遗传力公式计算得到狭义遗传力
h2_1 <- var_g_1 / (var_g_1 + var_e_1)
Bayes模型-BGLR包
贝叶斯模型全部使用BGLR包，差异仅限于参数上分别为“BayesA” “BayesB” “BayesC” “BRR”
读取基因型和表型与上文一致（无需计算K矩阵，A.mat()是rrBLUP的函数）
ETA <- list(list(X=geno, model="BayesA"))
fit <- BGLR(y=pheno1, ETA=ETA, nIter=12000, burnIn=2000, saveAt="BA_")
#提取方差
var_g <- fit$ETA[[1]]$varB
var_e <- fit$varE
h2 <- var_g / (var_g + var_e)
RKHS模型-BGLR包
RKHS模型也使用BGLR包，但需要提前计算核函数，参数设置为“RKHS”
K <- kernelMatrix(rbfdot(sigma = 0.5), as.matrix(geno))
K <- K / mean(diag(K))
ETA <- list(list(K = K, model = "RKHS"))
RKHS_10 <- BGLR(y = pheno10, ETA = ETA, nIter = 12000, burnIn = 2000)
var_g_rkhs1 <- RKHS_1$ETA[[1]]$varU
var_e_rkhs1 <- RKHS_1$varE
h2_rkhs1 <- var_g_rkhs1 / (var_g_rkhs1 + var_e_rkhs1)
BGLR包：nIter和burnIn参数的设置
nIter（总迭代次数）和 burnIn（预烧期迭代次数）是 马尔可夫链蒙特卡洛 (MCMC) 算法的关键参数，直接影响结果的收敛性和可靠性。
一般经验值
小数据集（n < 500, p < 10K）：
nIter = 10,000–20,000
burnIn = 2,000–5,000
大数据集（n > 1K 或 p > 100K）：
nIter = 50,000–100,000
burnIn = 10,000–20,000
LASSO/EN模型-glmnet包
LASSO和EN都使用glmnet包，当alpha = 1为LASSO模型，alpha = 0.5为EN模型
X <- scale(geno)
fit_lasso_1 <- cv.glmnet(X, pheno1, alpha = 1)
#EN fit_lasso_1 <- cv.glmnet(X, pheno1, alpha = 0.5)
y_pred_1 <- predict(fit_lasso_1, X, s = "lambda.min")
h2_1 <- var(y_pred_1) / var(pheno1)