BLUE和BLUP的比较及使用R语言计算它们

听我说一句：这是我的学习笔记，仅供参考

现在的分析喜欢计算表型值

目前就只有三个答案计算表型值：

1 平均值

2 BLUE值 (最佳线性无偏估计, 固定因子)

3 BLUP值 (最佳线性无偏预测, 随机因子)

1 各自的原理

1.1 BLUE

植物育种，感觉多年多点的数据，其实应该用BLUE，但是目前大多数都用BLUP了。在论坛的讨论如下，

论坛来源：https://www.researchgate.net/post/Does-any-one-have-an-idea-of-which-one-BLUE-or-BLUP-to-use-for-a-GWAS-analysis-of-a-trait-in-wheat-eg-resistance-to-rust

问题：

回答1：

回答2：

论坛说用BLUE好

1.2 BLUP

动物育种用这个多

参考文章1：https://www.omicsclass.com/article/1380

参考文章2：https://cloud.tencent.com/developer/article/1771941

BULP的优点：

• 能有效地充分利用所有亲属的信息
• 能更有效地校正环境效应
• 能校正由于非随机交配（选择交配）造成的偏差
• 能考虑不同群体及不同世代的遗传差异（注意：不同群体联合评定时，群体间有一定的遗传联系）
• 能提供个体育种值的最精确的无偏估计值

1.3 BLUE值和BLUP值

BLUP —— Best Liner Unbiased Prediction 最佳线性无偏预测

BLUE —— 最佳线性无偏估计

• **「BLUE」**值，相当于是对混合线性模型中固定因子的估算
• **「BLUP」**值，相当于是对混合线性模型中随机因子的预测

BLUE值一般是矫正的表型值，尺度和表型值一致，如果是多个重复或者多年多点的数据，可以将其代替平均值进行相关GS和GWAS的分析。

BLUP值一般用于品种排名，品种选种时的依据。

2 计算BLUP育种值（多环境无重复）

使用R语言计算BLUP。

多环境：意味着实验在多个不同的环境中进行，比如不同的地点、年份、栽培条件等。

无重复：在每个环境中，基因型样本没有重复。这意味着同一个基因型在某个特定环境下只出现一次，没有进行重复实验。

2.1 读入数据

数据格式

 lines env   y                 
  <chr> <chr> <chr>             
1 sanple1 env1  29.2              
2 sanple2 env1  25.399999999999995
3 sanple3 env1  NA                
4 sanple4 env1  27.1              
5 sanple5 env1  28.633333333333336
6 sanple6 env1  29.933333333333334
7 sanple1 env2  29              
8 sanple2 env3  25
9 sanple3 env4  27                
10 sanple4 env5  NA              
11 sanple5 env6  28.6
12 sanple6 env7  29.93

getwd()
library(lme4)
library(readxl)
data <- read_excel("blup_DATA.xlsx",sheet = 1,col_names = TRUE)
head(data)
# 使用 sapply 查看数据框每列的类
sapply(data, class)

2.2 转换因子—数据清洗

#因子是一种用于分类数据的特殊类型，通常用于表示分类变量。
# 将 data$lines 列转换为因子类型
data$lines=factor(data$lines)
# 将 data$env 列转换为因子类型
data$env=factor(data$env)

##报错发现没有将y列转换成数值
#将 y 列转换为数值类型：
data$y = as.numeric(data$y)
# 删除含有NA值的行
data <- na.omit(data)

#再次查看一下数据类型
sapply(data, class)

2.3 BLUP的计算

用lmer进行BLUP分析,在lmer中 1|env 表示把env当作随机效应，我们把env和lines当作随机效应。

blp=lmer(y~(1|env)+(1|lines),data=data)
blp

计算结果

> summary(blp)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ (1 | env) + (1 | lines)
   Data: data

REML criterion at convergence: 17441.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0656 -0.5439 -0.0272  0.5378  5.7567 

Random effects:
 Groups   Name        Variance Std.Dev.
 lines    (Intercept) 17.3720  4.1680  
 env      (Intercept)  0.2514  0.5014  
 Residual             11.8308  3.4396  
Number of obs: 3039, groups:  lines, 615; env, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  28.2340     0.2871   98.33

2.4 计算遗传力

遗传力的计算公式

#遗传力的计算
H = 17.3720/(17.3720 + 11.8308/5)
H

数据来源

2.5 获取BLUP育种值

ranef函数获取随机效应值。

blups返回一个list，包含env和lines的随机效应值即BLUP。

blp@beta为整体均值。

#获取BLUP值
blups= ranef(blp)
names(blups)
lines=blups$lines+blp@beta
res=data.frame(id=rownames(lines),blup=lines)

#保存BLUP育种值
write.table(res,file="data_blup_result.txt",row.names = F,quote = F,sep="\t")

2.6 描述性统计

#画一个柱状图看看
hist(lines[,1],col="black",border="white",xlab="BLUP of lines",main="")
#描述统计读入数据
PL_BLUP <- read.table("data_blup_result.txt",header = T,sep = "\t")
#install.packages("fBasics")
library(fBasics)
#描述统计
summary(data_BLUP$X.Intercept.)
var(data_BLUP$X.Intercept.) #方差
sd(data_BLUP$X.Intercept.) #标准差
skewness(data_BLUP$X.Intercept.) #偏度
kurtosis(data_BLUP$X.Intercept.) #计算峰度

在r中没有直接函数可以调用，但是有两个包可以使用：moments、fBasics

这两个包区别是：峰度moments没有减3，fBasics减3

3 计算BLUE育种值

该计算方式是昨年写的，又复习了一遍

3.1 读入数据

Line:品种名 Rep：重复（基因型重复） Year：种植年份 Loc：地理位置 Trait：表型数据

从上往下垒：

Line	 Rep	Year	Loc	Trait
C005	3	2018	SZ	110.8
C007	3	2018	SZ	101.6
C009	3	2018	SZ	98.2
C015	3	2018	SZ	94
C018	3	2018	SZ	101
C019	3	2018	SZ	100.2
C020	3	2018	SZ	NA
C021	3	2018	SZ	124.2
C023	3	2018	SZ	94.4
C005	3	2019	BJ	110
C007	3	2019	BJ	101.63
C009	3	2019	BJ	98.35
C015	3	2019	BJ	94.89
C018	3	2019	BJ	101
C019	3	2019	BJ	100
C020	3	2019	BJ	102
C021	3	2019	BJ	NA
C023	3	2019	BJ	94.8

3.2 计算BLUE和遗传力

getwd()
library(lme4)
options(lmerControl=list(check.nobs.vs.rankZ ="warning",
                         check.nobs.vs.nlev = "warning",
                         check.nobs.vs.nRE= "warning",
                         check.nlev.gtreq.5 = "warning",
                         check.nlev.gtr.1= "warning"))
phe <- read.table("你的数据",header = T,
                  sep = "\t")
trait <- phe$Trait
line <-  phe$Line
loc <-  phe$Loc
Rep <- phe$Rep
year <- phe$Year
TRAIT <- as.numeric(trait)
LINE <-  as.factor(line)
LOC <- as.factor(loc)
REP <-  as.factor(Rep)
YEAR <-  as.factor(year)
blue <-  lmer(TRAIT~ line + (1|loc) + (1|loc:Rep) + (1|year),data = phe)
as.data.frame(fixef(blue))
summary(blue)
#install.packages("lsmeans")
library(lsmeans)
re = lsmeans(blue,"line") 
re
write.table(re,"re_BLUE.txt",sep = "\t")

#这个方法也能直接计算遗传力：
blue_h <- lmer(TRAIT~(1|LINE)+(1|LOC)+
               (1|YEAR)+(1|REP%in%LOC:YEAR)+(1|LINE:LOC)+(1|LINE:YEAR))
summary(blue_h)
#计算公式：H2=Vg/{Vg+(VGL/L)+(VGY/Y)+(VE/LYR)}，
#其中H2为广义遗传力，
#Vg、VGL、VGY、VE分别为遗传方差、基因型与地点互作方差、基因型与年份互作方差、残差方差，
#L、Y、R为地点数、年份数、重复数。

##line是遗传方差，LINE:LOC 是基因型与地点互作方差Vgl，
##LINE:YEAR是品种与年份互作方差Vgy，Residual是残差方差Ve，
##以上都是取第一列的数值，重复数R，年份数Y，环境个数L

  (1|YEAR)+(1|REP%in%LOC:YEAR)+(1|LINE:LOC)+(1|LINE:YEAR))

summary(blue_h)
#计算公式：H2=Vg/{Vg+(VGL/L)+(VGY/Y)+(VE/LYR)}，
#其中H2为广义遗传力，
#Vg、VGL、VGY、VE分别为遗传方差、基因型与地点互作方差、基因型与年份互作方差、残差方差，
#L、Y、R为地点数、年份数、重复数。

##line是遗传方差，LINE:LOC 是基因型与地点互作方差Vgl，
##LINE:YEAR是品种与年份互作方差Vgy，Residual是残差方差Ve，
##以上都是取第一列的数值，重复数R，年份数Y，环境个数L