使用的包：
data.table：用于读取大文件
sommer：用于进行最核心的线性求解计算
lme4：用于部分计算的支持
实例：使用杜鹃花花序花朵数量性状计算基因型遗传力
library(data.table)
library(lme4)
library(sommer)
flowerNUM <- read.table("260sample_flowerNUM.txt",header=F) #读取表型数据
leafW <- read.table("260sample_leafW_new.txt",header=F)
colnames(flowerNUM)[1:2]<-c("ID","phenotype")
colnames(leafW)[1:2]<-c("ID","phenotype")#处理后的flowerNUM和leafW应该是含有ID和phenotype两个数据列的数据框（data.frame）flowerNUM$ID <- leafW$ID
#flowerNUM的样本号似乎有缺失，为了保证效率使用leafW的样本ID直接替换
G <- fread("G.012.imputed.GAPIT.3.txt",header=T) #读取基因型数据
G <- G[,-1] #读取的基因型数据第一列为样本名，需要删除
rownames(G) <- leafW$ID #使用leafW的样本名进行替换
G <- apply(G,2,as.numeric) #因为fread()函数读取后的数据结构并不与read.table()一致，需要手动处理
A <- A.mat(G) #使用sommer包的A.mat()函数进行加性亲缘矩阵的计算
rownames(A)<- colnames(A)<- leafW$ID #对A矩阵的样本进行命名，此时A矩阵应该是一个260*260的可逆矩阵
phenotype <- as.numeric(flowerNUM$phenotype)#从flowerNUM数据框中提取一个表型数据的向量作为固定效应
model_gblup <- mmer( phenotype ~1, ~ vsr(sample.ID, A), ~ units, flowerNUM, tolParInv = 1e-6)#mmer()函数的各部分分别为：固定效应、随机效应、表型数据框、防止奇异参数
#计算可能因矩阵奇异出现以下结果，此时增加tolParInv参数的值直到完成计算为止：
model_gblup.flowerNUM <- mmer( phenotype ~1, ~ vsr(sample.ID, A), ~ units, flowerNUM, tolParInv = 1e-3)#model_gblup.flowerNUM中的结果可以通过两种办法获得：手动计算和自动计算
#手动计算（提取数据后计算得到）
var_components <- summary(model_gblup.flowerNUM)$varcomp
sigma_g2 <- var_components$VarComp[1]
sigma_e2 <- var_components$VarComp[2]
heritability <- sigma_g2 /(sigma_g2 + sigma_e2)#此外还需要计算三个值的标准误差 SD
sigma_g2_SD <- var_components$VarCompSE[1]
sigma_e2_SD <- var_components$VarCompSE[2]
heritability_SE <-sqrt((1/(sigma_g2 + sigma_e2) - sigma_g2 / (sigma_g2 + sigma_e2)^2)^2* sigma_g2_SD^2+(-sigma_g2 /(sigma_g2 + sigma_e2)^2)^2* sigma_e2_SD^2)#使用vpredict()函数计算（同时得到计算结果和标准误差）
h2_vpredict <- vpredict(model_gblup, h2 ~ V1 / (V1 + V2))
Vg_vpredict <- vpredict(model_gblup, Vg ~ V1)
Ve_vpredict <- vpredict(model_gblup, Ve ~ V2)#推荐直接使用vpredict()函数进行计算，使用模型估计的方差成分，引入了更加复杂的模型
#如果想要输出数字，可以用该指令（根据上面的结果替换变量名即可）
print(paste0("Vg(SD):", sigma_g2,"(", sigma_g2_SD,")"))
print(paste0("Ve(SD):", sigma_e2,"(", sigma_e2_SD,")"))
print(paste0("Heritability(SD):", heritability,"(",heritability_SE,")"))