用于gwas的blup与遗传力计算

1，数据准备

（使用文件夹下以测量性状命名的tsv文件）

2.代码

library(lme4)
library(data.table)

# 设置工作目录（根据需要调整）
setwd("E:/pyProjects/zc/zc_gwas2")

# 创建result文件夹（如果不存在）
if (!dir.exists("blup_result")) {
  dir.create("blup_result")
}

# 定义性状列表
traits <- c("TSN", "FTN", "TKW", "NGS", "FSN", "AL")
# traits <- "FSN"

# 遍历每个性状
for (trait in traits) {
  # 导入表型数据
  phe <- fread(paste0('traitplus/', trait, ".csv"), data.table = FALSE)

  # 准备数据
  TRAIT <- as.numeric(phe[[trait]])
  LINE <- as.factor(phe$Cul)
  LOC <- as.factor(phe$loc)
  REP <- as.factor(phe$Rep)

  # 建立模型
  blup_model <- lmer(TRAIT ~ (1 | LINE) +
    (1 | LOC) +
    (1 | REP / LOC) +
    (1 | LINE:LOC))


  # 获取结果摘要
  result <- summary(blup_model)

  # 准备计算数据
  result_g <- as.list(result$ngrps)
  print(result_g)
  vg <- as.numeric(VarCorr(blup_model)$LINE[1])
  vge <- if ("LINE:LOC" %in% names(VarCorr(blup_model))) {
    as.numeric(VarCorr(blup_model)$`LINE:LOC`[1])
  } else {
    0
  }
  ve <- as.numeric(VarCorr(blup_model)$LOC[1])
  L <- as.numeric(result_g$LOC)
  R <- as.numeric(max(phe$Rep))

  # 计算广义遗传力
  h <- vg / (vg + vge / L + ve / (L * R))

  # 提取BLUP值
  blup_values <- ranef(blup_model)$LINE
  blup_df <- data.frame(LINE = rownames(blup_values), BLUP = blup_values[[1]])

  # 将结果保存到文件
  write.csv(blup_df, file = paste0("blup_result/", trait, "_BLUP.csv"), row.names = FALSE)

  # 将遗传力保存到文件
  write.table(data.frame(Trait = trait, Heritability = h),
              file = "blup_result/Heritability.csv",
              append = TRUE,
              sep = ",",
              row.names = FALSE,
              col.names = !file.exists("blup_result/Heritability.csv"))

  # 打印进度
  cat(paste("Finished processing", trait, "\n"))
}

cat("All traits have been processed. Results are saved in the 'result' folder.\n")

 BLUP的模型解释：

模型解释：
 

• TRAIT：这是因变量，即我们要研究的性状。
• (1 | LINE)：这表示品系（基因型）的随机效应。每个品系都有一个与之相关的随机效应。
• (1 | LOC)：这表示位置（环境）的随机效应。每个位置都有一个与之相关的随机效应。
• (1 | REP / LOC)：这表示重复嵌套在位置中的随机效应。每个位置内的每个重复都有一个随机效应。
• (1 | LINE:LOC)：这表示品系和位置之间的交互作用，即基因型与环境的交互作用（G×E）。

从数学角度来看，这个模型可以表示为：

Y_ijkl = μ + L_i + E_j + R_k(j) + (LE)_ij + ε_ijkl

其中：

• Y_ijkl 是观察到的表型值
• μ 是总体平均值
• L_i 是第i个品系的随机效应，假设 L_i ~ N(0, σ²_L)
• E_j 是第j个环境的随机效应，假设 E_j ~ N(0, σ²_E)
• R_k(j) 是第j个环境中第k个重复的随机效应，假设 R_k(j) ~ N(0, σ²_R)
• (LE)_ij 是品系和环境的交互作用，假设 (LE)_ij ~ N(0, σ²_LE)
• ε_ijkl 是随机误差，假设 ε_ijkl ~ N(0, σ²_ε)

计算方法：

1. 模型使用最大似然估计（ML）或限制最大似然估计（REML）来估计方差组分。
2. BLUP（最佳线性无偏预测）值通过解混合模型方程来获得： [X'R^(-1)X X'R^(-1)Z] [b] = [X'R^(-1)y] [Z'R^(-1)X Z'R^(-1)Z+G^(-1)] [u] [Z'R^(-1)y] 其中，X和Z是设计矩阵，R是残差协方差矩阵，G是随机效应的协方差矩阵，b是固定效应，u是随机效应。
3. 广义遗传力计算如下： h² = σ²_L / (σ²_L + σ²_LE/n_loc + σ²_ε/(n_loc * n_rep)) 其中，n_loc是位置数，n_rep是重复数。

 遗传力的计算方法

遗传力的生物学含义：

遗传力是一个重要的遗传学参数，用于描述一个性状的表型变异中由遗传因素引起的比例。它反映了一个群体中，某个性状的遗传变异在总变异中所占的比例。遗传力的值介于0到1之间：
- 接近1表示该性状主要受遗传因素控制
- 接近0表示环境因素对该性状的影响更大

在育种中，遗传力高的性状更容易通过选择来改良。

广义遗传力的计算方法：

在这个脚本中，我们计算的是广义遗传力（broad-sense heritability）。计算公式如下：

h² = Vg / (Vg + Vge/L + Ve/(L*R))

其中：
- h² 是广义遗传力
- Vg 是遗传变异（基因型变异）
- Vge 是基因型与环境的交互作用变异
- Ve 是环境变异（包括误差）
- L 是位置（环境）数量
- R 是重复数量

让我们逐步解析这个公式：

1. Vg（遗传变异）：
   在模型中，这对应于 LINE 随机效应的方差组分。
   vg <- as.numeric(VarCorr(blup_model)$LINE[1])
 

2. Vge（基因型与环境交互作用变异）：
   这对应于 LINE:LOC 交互作用的方差组分。
   vge <- if ("LINE:LOC" %in% names(VarCorr(blup_model))) {
     as.numeric(VarCorr(blup_model)$`LINE:LOC`[1])
   } else {
     0
   }
 

3. Ve（环境变异）：
   在这个模型中，我们使用 LOC 随机效应的方差组分来代表环境变异。
   ve <- as.numeric(VarCorr(blup_model)$LOC[1])
 

4. L 和 R：
   L <- as.numeric(result_g$LOC)
   R <- as.numeric(max(phe$Rep))
 

最后，将这些组分代入公式计算遗传力：
h <- vg / (vg + vge / L + ve / (L * R))

这个公式的数学含义：

- 分子 Vg 代表了纯粹的遗传变异。
- 分母是总变异，包括：
  - 遗传变异 (Vg)
  - 平均每个环境的基因型与环境交互作用变异 (Vge/L)
  - 平均每个环境和重复的环境变异 (Ve/(L*R))

通过这种方式，我们可以估计在多环境试验中，遗传因素对表型变异的贡献程度。这个估计考虑了不同环境和重复的影响，使得遗传力的估计更加准确和有代表性。

高遗传力意味着表型更多地受基因控制，选择效率会更高；低遗传力则表明环境因素的影响更大，可能需要更精确的表型评估或更大的群体规模来进行有效选择。