基因型数据：使用rrblup的A.mat()函数计算得到的G矩阵，也可以使用其他软件得到KINSHIP或G矩阵，二者等价。G矩阵的规模应为n*n（n为数据集的个体数）
表型数据：命名好的unlist表型数据向量或数据框，对应n个个体的数据
变量结构实例：
> str(G)
num [1:260, 1:260] 0.69585 -0.0362 0.01101 -0.01341 0.00942 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:260] "H-R100" "H-R101" "H-R102" "H-R103" ...
..$ : chr [1:260] "H-R100" "H-R101" "H-R102" "H-R103" ...
> str(pheno1)
Named int [1:260] 13 22 10 9 10 16 8 4 15 12 ...
- attr(*, "names")= chr [1:260] "H-R100" "H-R101" "H-R102" "H-R103" ...
library(sommer) #加载sommer包
df <- data.frame(ID = names(pheno1), Y = as.numeric(pheno1)) #将表型转为数据框结构
df$ID <- as.factor(df$ID) #将ID作为因子
fit <- mmer(fixed = Y ~ 1, random = ~ vs(ID, Gu = G), data = df) #使用mmer()函数，其中G为随机变量，Y为固定效应
var.u <- fit$sigma$`u:ID` #提取varU
var.e <- fit$sigma$`units` #提取varE
h2 <- var.u / (var.u + var.e) #计算狭义遗传力
h2
library(rrBLUP) #加载rrBLUP包
fit <- mixed.solve(y = pheno1, K = G) #使用mixed.solve()函数进行计算，y为表型数据，G为G矩阵
sigma_a2 <- fit$Vu #提取varU
sigma_e2 <- fit$Ve #提取varE
h2 <- sigma_a2 / (sigma_a2 + sigma_e2) #计算狭义遗传力
h2
额外补充：表型可以不使用均值替换
在遗传学数据处理中，不建议直接使用均值替换离群值来计算遗传力。线性混合模型（LMM）是通过个体间的“差异”与亲缘关系的“相关性”来估算方差的。均值是分布的中心。如果你把离群值改为均值，相当于在数据中插入了几个“完全没有变异”的点。软件会认为这些原本有差异的个体突然变得一模一样，但这并不是由于遗传引起的，这会扭曲 Vg 和 Ve 的比例，导致估计出的遗传力失去生物学意义。更好的替代方案应为设置为 NA。
因此将 pheno8（leafW）和 pheno10（leafLW）原本的三个填补个体 H-R223、H-R230、H-R233 的表型数据替换为NA后重新进行遗传力的计算：
| sommer | rrBLUP | NA | |
| pheno8 | 0.254436 | 0.254478 | 0.25583 |
| pheno10 | 0.209271 | 0.210463 | 0.215488 |
替换后的结果有少量提升，但提升有限。这说明表型潜在的真实偏离值并不多，因此不论使用哪种方式，对GWAS等的结果影响有限，可以忽略不计。
pheno6(corollaW)的0遗传力问题
# --- 1. 计算变异系数 (CV) ---
# CV = (标准差 / 平均值) * 100%
x <- as.numeric(pheno6)
cv_val <- (sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)) * 100
cv_val
[1] 14.42179
# --- 2. 检查极端离群值 (Outliers) ---
# 使用 IQR 方法：超过 [Q1 - 1.5*IQR, Q3 + 1.5*IQR] 的值为离群值
quartiles <- quantile(x, probs=c(.25, .75), na.rm = TRUE)
IQR_val <- IQR(x, na.rm = TRUE)
lower_bound <- quartiles[1] - 1.5 * IQR_val
upper_bound <- quartiles[2] + 1.5 * IQR_val
outliers <- x[x < lower_bound | x > upper_bound]
outliers
[1] 31.8 32.0 32.0 37.0 12.0
# --- 3. 确定是否为正态分布 (Normality Test) ---
# Shapiro-Wilk 检验：p > 0.05 说明符合正态分布
shapiro_test <- shapiro.test(x)
shapiro_test
Shapiro-Wilk normality test
data: x
W = 0.97832, p-value = 0.0005341
核心问题：显著的非正态性与离群值
• 离群值的干扰：CW(pheno6)识别出了 5 个离群值（31.8, 32.0, 32.0, 37.0, 12.0）。在混合模型中，方差计算是对差值的平方进行求和。这些极端值（尤其是 37.0 和 12.0）会产生巨大的残差变异。计算中，软件会将这些无法由 G 矩阵解释的剧烈波动全部归类到残差方差 σe2 中。由于分母变大，遗传力 h2自然趋于 0。
• 不符合正态分布 (p < 0.001)：Shapiro-Wilk 检验的 p 值远小于 0.05，说明CW(pheno6)的分布具有显著的偏态或峰度异常。线性混合模型（LMM）的前提假设是残差服从正态分布，当数据严重偏离正态时，模型对 σa2 的估计往往会失效。
• 正常的变异系数：14.42% 的变异系数说明该性状本身是有变异的（通常 CV > 10% 被认为具有一定的变异程度）。
pheno6_fixed <- pheno6
outliers_values <- c(31.8, 32.0, 32.0, 37.0, 12.0)
pheno6_fixed[pheno6_fixed %in% outliers_values] <- NA
fit <- mixed.solve(y = pheno6_fixed, K = G)
sigma_a2 <- fit$Vu
sigma_e2 <- fit$Ve
h2 <- sigma_a2 / (sigma_a2 + sigma_e2)
h2
[1] 1.000033e-09