zoukankan      html  css  js  c++  java
  • 【GS模型】使用R包sommer进行基因组选择的GBLUP和RRBLUP分析?

    简介

    R包sommer内置了C++,运算速度还是比较快的,功能也很丰富,可求解各种复杂模型。语法相比于lme4包也要好懂一些。
    建议查看文档:vignette("v1.sommer.quick.start")
    image.png

    混合线性模型关键在于协方差结构的建立,有以下几类:

    • 复合对称(Compound Symmetry,CS),所有方差相等,所有协方差也相等,对应于单变量方法。但是对于不同尺度的变量是无意义的.
    • 方差组分(Variance Components),每个方差都不相同,并且全部协方差等于0。如果变量完全独立,并且彼此测量尺度不同,才是有意义的模式。
    • 非结构化(Unstructured),没有模式,每个方差和协方差完全不同,彼此间没有关系,对应于多变量方法。协方差结构有很多种,只有在特定的统计条件下才有意义。

    调用模型并不难,难的是在理解的基础上如何随心所欲地应用。

    GS示例代码

    • 预处理
    library(sommer)
    
    data(DT_wheat)
    DT <- DT_wheat
    GT <- GT_wheat
    dim(GT)
    GT[1:10,1:10]
    
    colnames(DT) <- paste0("X",1:ncol(DT))
    DT <- as.data.frame(DT)
    DT$id <- as.factor(rownames(DT))
    rownames(GT) <- rownames(DT)
    
    # check NA
    which(is.na(GT))
    which(is.na(DT))
    
    set.seed(12345)
    y.trn <- DT
    
    #制造1/5的缺失值,作为验证集
    vv <- sample(rownames(DT),round(nrow(DT)/5))
    y.trn[vv,"X1"] <- NA
    y.trn[1:5,1:5]
    
    • GBLUP
    #######-----------GBLUP--------------------------------------
    # GBLUP pedigree-based approach
    K <- A.mat(GT) # additive relationship matrix
    K[1:5,1:5]
    colnames(K) <- rownames(K) <- rownames(DT)
    #test first trait X1
    system.time(
      ans <- mmer(X1~1, 
                random=~vs(id,Gu=K),
                rcov=~units,
                data=y.trn,verbose = FALSE) # kinship based
    )
    ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
    rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
    
    cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
    

    image.png

    • RRBLUP
    #######-----------rrBLUP--------------------------------------
    system.time(
      ans2 <- mmer(X1~1,
                   random=~vs(list(GT)),
                   rcov=~units,
                   data=y.trn,verbose = FALSE) # kinship based
    )
    u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
    rownames(u) <- rownames(GT)
    cor(u[vv,],DT[vv,"X1"]) # same correlation
    

    image.png

    两者结果相差不大(如果去掉随机种子,循环运行的结果相差还是很大的),运算时间相差比较大。

    Ref: 协方差矩阵,协方差结构

  • 相关阅读:
    在Windows10上安装Linux子系统
    《ln--软连接》--linux命令
    【转】最详细的Log4J使用教程
    7 str字符串的内置函数
    6.1 range创建数值序列
    6 int数值类型的内置函数
    linux shell 字符串操作(长度,查找,截取,替换)详解
    awk 字符串处理函数
    Shell字符串截取
    5 python基本数据类型
  • 原文地址:https://www.cnblogs.com/jessepeng/p/14207126.html
Copyright © 2011-2022 走看看