遗传力又叫遗传率，分为广义遗传力（board heritability）和狭义遗传力（narrow sense heritability）。广义遗传力一般用H2表示，分子是所有的遗传方差，分母是总方差（包括遗传方差、环境方差和误差方差）。相应的，狭义遗传力用h2表示，分子是加性效应的遗传方差，分母是总方差。
（1）完全双列杂交设计
模型为：
y = Xβ + Zu1 + Zu2 + ZuS + ε这里，X是固定效应的发生率矩阵，β是固定效应向量，Z是随机效应的发生率矩阵，u1和u2分别是杂种优势群GCA1和GCA2的随机效应向量，uS是SCA的随机效应向量，ε是残差向量。GCA和SCA的BLUP可以用于预测杂交。
计算代码为：
library(sommer) # 加载sommer包，第一次安装请先运行install.packages("sommer")
data(DT_cornhybrids) # 加载数据集
DT <- DT_cornhybrids # 杂交数据矩阵赋值给DT变量
DTi <- DTi_cornhybrids # 亲本信息矩阵赋值给DTi变量，计算遗传力的时候用不到
GT <- GT_cornhybrids # 数值化的基因型矩阵赋值给GT变量，计算遗传力时用不到
modFD <- mmer(Yield~Location, # 用地点做固定效应
random=~GCA1+GCA2+SCA, # GCA1、GCA2和SCA做随机效应
rcov=~units, #
data=DT, verbose = FALSE) # 数据集使用DT，不显示迭代过程
(suma <- summary(modFD)$varcomp) # 获得方差组分并显示
Vgca <- sum(suma[1:2,1]) # 两个GCA方差的和
Vsca <- suma[3,1] # SCA方差的
Ve <- suma[4,1] # 环境方差的
Va = 4*Vgca # 4个环境×GCA方差 ？？？为什么不直接用GCA而要×4？
Vd = 4*Vsca # 4个环境×SCA方差 ？？？为什么不直接用SCA？
Vg <- Va + Vd # 遗传方差=加性方差+显性方差
(H2 <- Vg / (Vg + (Ve)) ) # 广义遗传力=遗传方差/总方差，Ve外面的括号是多余的
(h2 <- Va / (Vg + (Ve)) ) # 狭义遗传力=加性方差/总方差
本例中，使用来自2个杂种优势群的40个自交系，每个杂种优势群有20个自交系。因此，共有20×20=400个可能的组合。该模拟数据在4个地点鉴定了100个系的表型数据，每个地点只有1次重复，共计4×100=400条数据。
下面是杂交数据DT的前100行（1600×6）：
head(GT,100)
Location GCA1 GCA2 SCA Yield PlantHeight
1 1 A258 AS5707 A258:AS5707 NA NA
2 1 A258 B2 A258:B2 NA NA
3 1 A258 B99 A258:B99 NA NA
4 1 A258 LH51 A258:LH51 NA NA
5 1 A258 Mo44 A258:Mo44 NA NA
6 1 A258 NC320 A258:NC320 NA NA
7 1 A258 Oh40B A258:Oh40B 129.0082 1.632108
8 1 A258 W117Ht A258:W117Ht 139.5254 1.613334
9 1 A258 W182BN A258:W182BN 122.8658 1.688060
10 1 A258 W611S A258:W611S NA NA
11 1 A258 A641 A258:A641 166.6094 2.010321
12 1 A258 B10 A258:B10 159.3310 1.947324
13 1 A258 B119 A258:B119 NA NA
14 1 A258 B14A A258:B14A NA NA
15 1 A258 CM174 A258:CM174 NA NA
16 1 A258 H91 A258:H91 NA NA
17 1 A258 N7A A258:N7A NA NA
18 1 A258 PHG86 A258:PHG86 NA NA
19 1 A258 R226 A258:R226 NA NA
20 1 A258 W610S A258:W610S 147.6366 1.564184
21 1 B118 AS5707 B118:AS5707 NA NA
22 1 B118 B2 B118:B2 126.9928 2.027177
23 1 B118 B99 B118:B99 NA NA
24 1 B118 LH51 B118:LH51 NA NA
25 1 B118 Mo44 B118:Mo44 135.6515 1.862859
26 1 B118 NC320 B118:NC320 NA NA
27 1 B118 Oh40B B118:Oh40B NA NA
28 1 B118 W117Ht B118:W117Ht NA NA
29 1 B118 W182BN B118:W182BN NA NA
30 1 B118 W611S B118:W611S 130.1698 1.629009
31 1 B118 A641 B118:A641 156.6514 2.130744
32 1 B118 B10 B118:B10 NA NA
33 1 B118 B119 B118:B119 NA NA
34 1 B118 B14A B118:B14A NA NA
35 1 B118 CM174 B118:CM174 NA NA
36 1 B118 H91 B118:H91 NA NA
37 1 B118 N7A B118:N7A 145.3416 1.525327
38 1 B118 PHG86 B118:PHG86 NA NA
39 1 B118 R226 B118:R226 NA NA
40 1 B118 W610S B118:W610S 165.0516 1.603079
41 1 B97 AS5707 B97:AS5707 119.5537 2.031154
42 1 B97 B2 B97:B2 NA NA
43 1 B97 B99 B97:B99 NA NA
44 1 B97 LH51 B97:LH51 121.0262 1.728354
45 1 B97 Mo44 B97:Mo44 141.8828 2.007985
46 1 B97 NC320 B97:NC320 NA NA
47 1 B97 Oh40B B97:Oh40B 121.4198 1.555613
48 1 B97 W117Ht B97:W117Ht 141.3195 1.525166
49 1 B97 W182BN B97:W182BN 129.2736 1.558707
50 1 B97 W611S B97:W611S NA NA
51 1 B97 A641 B97:A641 NA NA
52 1 B97 B10 B97:B10 NA NA
53 1 B97 B119 B97:B119 NA NA
54 1 B97 B14A B97:B14A 146.8906 1.884599
55 1 B97 CM174 B97:CM174 NA NA
56 1 B97 H91 B97:H91 NA NA
57 1 B97 N7A B97:N7A NA NA
58 1 B97 PHG86 B97:PHG86 NA NA
59 1 B97 R226 B97:R226 NA NA
60 1 B97 W610S B97:W610S NA NA
61 1 C102 AS5707 C102:AS5707 NA NA
62 1 C102 B2 C102:B2 134.6043 1.918804
63 1 C102 B99 C102:B99 NA NA
64 1 C102 LH51 C102:LH51 148.6455 2.105199
65 1 C102 Mo44 C102:Mo44 NA NA
66 1 C102 NC320 C102:NC320 NA NA
67 1 C102 Oh40B C102:Oh40B NA NA
68 1 C102 W117Ht C102:W117Ht NA NA
69 1 C102 W182BN C102:W182BN NA NA
70 1 C102 W611S C102:W611S NA NA
71 1 C102 A641 C102:A641 NA NA
72 1 C102 B10 C102:B10 NA NA
73 1 C102 B119 C102:B119 NA NA
74 1 C102 B14A C102:B14A NA NA
75 1 C102 CM174 C102:CM174 NA NA
76 1 C102 H91 C102:H91 NA NA
77 1 C102 N7A C102:N7A NA NA
78 1 C102 PHG86 C102:PHG86 NA NA
79 1 C102 R226 C102:R226 136.8265 1.684973
80 1 C102 W610S C102:W610S NA NA
81 1 LH61 AS5707 LH61:AS5707 NA NA
82 1 LH61 B2 LH61:B2 NA NA
83 1 LH61 B99 LH61:B99 NA NA
84 1 LH61 LH51 LH61:LH51 NA NA
85 1 LH61 Mo44 LH61:Mo44 NA NA
86 1 LH61 NC320 LH61:NC320 NA NA
87 1 LH61 Oh40B LH61:Oh40B NA NA
88 1 LH61 W117Ht LH61:W117Ht 126.6400 1.607918
89 1 LH61 W182BN LH61:W182BN NA NA
90 1 LH61 W611S LH61:W611S NA NA
91 1 LH61 A641 LH61:A641 NA NA
92 1 LH61 B10 LH61:B10 NA NA
93 1 LH61 B119 LH61:B119 NA NA
94 1 LH61 B14A LH61:B14A NA NA
95 1 LH61 CM174 LH61:CM174 NA NA
96 1 LH61 H91 LH61:H91 NA NA
97 1 LH61 N7A LH61:N7A NA NA
98 1 LH61 PHG86 LH61:PHG86 138.8772 1.652411
99 1 LH61 R226 LH61:R226 NA NA
100 1 LH61 W610S LH61:W610S NA NA下面是亲本信息数据DTi（40×4）：
DTi
Genotype Group Yield PlantHeight
1 A258 NSS 88.20679 0.9272210
2 A634 SSS 93.31265 1.0226329
3 A641 SSS 109.00473 1.0095320
4 A680 SSS 98.81964 1.1096407
5 AS5707 NSS 80.65517 1.3279382
6 B10 SSS 95.84411 1.3262020
7 B105 SSS 100.65840 1.1560672
8 B118 NSS 111.36915 1.4949915
9 B119 SSS 86.68451 0.9925757
10 B121 SSS 104.95121 0.6704880
11 B14A SSS 114.69768 1.0709695
12 B2 NSS 108.69982 1.1506838
13 B84 SSS 92.18648 1.2417499
14 B97 NSS 84.54403 1.2251734
15 B99 NSS 88.14189 1.3141365
16 C102 NSS 99.86213 1.3150211
17 CM174 SSS 107.16554 1.0458131
18 H105W SSS 94.94959 0.9217209
19 H91 SSS 105.75978 0.8894477
20 LH51 NSS 81.85276 1.2395905
21 LH61 NSS 109.39862 0.9673939
22 LP5 SSS 103.10739 0.9972820
23 Mo44 NSS 108.44876 0.8871609
24 N7A SSS 103.19176 1.4401722
25 NC262 NSS 99.52945 0.7187014
26 NC320 NSS 104.42666 1.2892791
27 NC344 NSS 96.68368 0.8644697
28 Oh40B NSS 102.63953 0.9376076
29 Pa91 NSS 124.03104 1.1982550
30 PHG71 SSS 83.70858 1.1174897
31 PHG86 SSS 90.85945 1.2136497
32 PHW17 SSS 104.01616 1.1334628
33 R226 SSS 95.31983 1.3587364
34 Tzi25 SSS 91.65843 1.0268473
35 W117Ht NSS 116.60487 0.9949743
36 W153R NSS 113.74265 1.5680403
37 W182BN NSS 89.71262 1.4102404
38 W23 NSS 93.20116 1.2686134
39 W610S SSS 106.17348 0.9484871
40 W611S NSS 92.02427 1.1937574基因型矩阵前6行和前6列（40×40）：
GT[1:6,1:6]
A258 A634 A641 A680 AS5707 B10
A258 2.23285528 -0.3504778 -0.04756856 -0.32239362 -0.07776163 -0.01257374
A634 -0.35047780 1.4529169 0.45203869 -0.02293680 -0.43538636 0.19984929
A641 -0.04756856 0.4520387 1.96940221 -0.09896791 -0.28059417 0.02019641
A680 -0.32239362 -0.0229368 -0.09896791 1.65221984 -0.33095920 0.12259252
AS5707 -0.07776163 -0.4353864 -0.28059417 -0.33095920 2.36536453 -0.18705982
B10 -0.01257374 0.1998493 0.02019641 0.12259252 -0.18705982 1.78689265
modFD <- mmer(Yield~Location,
random=~GCA1+GCA2+SCA,
rcov=~units,
data=DT, verbose = FALSE)
(suma <- summary(modFD)$varcomp)
VarComp VarCompSE Zratio Constraint
GCA1.Yield-Yield 0.000000 16.50337 0.0000000 Positive
GCA2.Yield-Yield 7.412226 18.94200 0.3913116 Positive
SCA.Yield-Yield 187.560303 41.59428 4.5092817 Positive
units.Yield-Yield 221.142463 18.14716 12.1860656 Positive
Vgca <- sum(suma[1,1])
Vsca <- suma[2,1]
Ve <- suma[3,1]
Va = 2*Vgca
Vd = 2*Vsca
Vg <- Va + Vd
(H2 <- Vg / (Vg + Ve) )
(h2 <- Va / (Vg + (Ve)) )
> (H2 <- Vg / (Vg + Ve) )
[1] 0.7790856
> (h2 <- Va / (Vg + (Ve)) )
[1] 0.02961832由于该数据主要模拟显性效应，因此，加性方差非常小，导致了h2很小。
（2）不考虑自交的半双列杂交
本例中有7个亲本，共有7×(7-1)/2=21个杂交组合。实验设计是2个重复的完全随机设计（completely randomized disign，CRD）。双亲的GCA和SCA的方差组分可以被估计，
模型为：
y = Xβ + Zug + Zus + ε代码为：
data("DT_halfdiallel") # 加载数据集
DT <- DT_halfdiallel # 将数据集赋值给DT变量
DT$femalef <- as.factor(DT$female) # 增加femalef列，作为因子
DT$malef <- as.factor(DT$male) # 增加malef列，作为因子
DT$genof <- as.factor(DT$geno) # 增加genof列，作为因子
#### model using overlay
modh <- mmer(sugar~1, # 斜率为固定效应，预测糖
random=~vs(overlay(femalef,malef)) # 亲本的设计矩阵+基因型向量为随机效应
+ genof,
data=DT, verbose = FALSE) # 数据集是DT，不显示迭代
(suma <- summary(modh)$varcomp)
summary(modh)$varcomp
Vgca <- sum(suma[1,1])
Vsca <- suma[2,1]
Ve <- suma[3,1]
Va = 2*Vgca
Vd = 2*Vsca
Vg <- Va + Vd
(H2 <- Vg / (Vg + Ve) )
(h2 <- Va / (Vg + (Ve)) )
VarComp VarCompSE Zratio Constraint
u:femalef.sugar-sugar 5.507899 3.5741151 1.541052 Positive
genof.sugar-sugar 1.815784 1.3629575 1.332238 Positive
units.sugar-sugar 3.117538 0.9626094 3.238632 Positive这里，vs()是构造方差协方差矩阵的主函数，里面包含了overlay()辅助函数。overlay()函数的作用是构造综合考虑两个变量的发生率矩阵，对于每个因子（本例中共7个），无论是父本还是母本，只要出现即为1，否则为0。
本例中发生率矩阵如下：
overlay(DT$femalef,DT$malef)
1 2 3 4 5 6 7
1 1 1 0 0 0 0 0
2 1 1 0 0 0 0 0
3 1 0 1 0 0 0 0
4 1 0 1 0 0 0 0
5 1 0 0 1 0 0 0
6 1 0 0 1 0 0 0
7 1 0 0 0 1 0 0
8 1 0 0 0 1 0 0
9 1 0 0 0 0 1 0
10 1 0 0 0 0 1 0
11 1 0 0 0 0 0 1
12 1 0 0 0 0 0 1
13 0 1 1 0 0 0 0
14 0 1 1 0 0 0 0
15 0 1 0 1 0 0 0
16 0 1 0 1 0 0 0
17 0 1 0 0 1 0 0
18 0 1 0 0 1 0 0
19 0 1 0 0 0 1 0
20 0 1 0 0 0 1 0
21 0 1 0 0 0 0 1
22 0 1 0 0 0 0 1
23 0 0 1 1 0 0 0
24 0 0 1 1 0 0 0
25 0 0 1 0 1 0 0
26 0 0 1 0 1 0 0
27 0 0 1 0 0 1 0
28 0 0 1 0 0 1 0
29 0 0 1 0 0 0 1
30 0 0 1 0 0 0 1
31 0 0 0 1 1 0 0
32 0 0 0 1 1 0 0
33 0 0 0 1 0 1 0
34 0 0 0 1 0 1 0
35 0 0 0 1 0 0 1
36 0 0 0 1 0 0 1
37 0 0 0 0 1 1 0
38 0 0 0 0 1 1 0
39 0 0 0 0 1 0 1
40 0 0 0 0 1 0 1
41 0 0 0 0 0 1 1
42 0 0 0 0 0 1 1
attr(,"variables")
[1] "DT$femalef" "DT$malef"
发表评论