标签: 遗传力计算

  • sommer包计算双列杂交的遗传力

    遗传力又叫遗传率,分为广义遗传力(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"