标签: R

  • Visual Studio Code+R

    Visual Studio Code是微软开发的开源、跨平台代码编辑器,号称宇宙最强!使用该编辑器直接编辑和运行R程序比使用RStudio占用资源更低,稳定性也好一些。

    R语言的下载

    进入R语言的官方网站,根据自己的所在地选择最适合的镜像站下载R(https://www.r-project.org/)。

    我选择的是兰州大学开源学会镜像(https://mirror.lzu.edu.cn/CRAN/),根据动画操作下载,版本号可能会有不同。

    如果电脑中已经安装了其他版本的R,建议先卸载,然后再安装。

    R语言的安装

    安装时,安装目录建议去掉横线,有些程序不识别目录中的短横线,容易造成错误。

    其他保持默认即可。

    下载VSCode和安装

    进入Visual Studio Code官网(https://code.visualstudio.com/),下载适合的版本,我这里以Windows 11 64 bit为例。安装保持默认即可。

    VSCode插件安装

    运行VSCode,按快捷键Ctrl+Shift+X调出扩展商店,在搜索栏输入“R”,安装R拓展。

    R语言安装languageserver包

    打开R,输入下列代码安装languageserver包,回车运行。

    install.packages("languageserver")
    

    在VSCode里运行R代码

    在VSCode中,选择菜单栏的文件→打开文件夹,选择一个经常编辑代码的文件夹。并信任该文件夹。

    按Ctrl+N新建一个文件,语言选择R。

    在第一行输入任意R代码,按Ctrl+Enter运行,查看结果。若能看到运行结果,则环境搭建完毕。

    getwd()
    
  • sommer包的GbyE分析

    sommer包可以做单环境模型、主效应模型、对角线模型(DG)、复合对称模型(CS)、复合对称+对角线模型(CS+DG)、非结构化模型(US)、随机回归模型(RR)和协方差结构G×E。当然,还可以有其他自己的组合。

    library(sommer)
    data(DT_example)
    DT <- DT_example  # 41个土豆材料,1个地点,3个年份=3个环境
    A <- A_example   # 系谱矩阵
    # 1. 单地点
    # BLUP
    ansSingle <- mmer(Yield~1,
                      random= ~ vs(Name, Gu=A),
                      rcov= ~ units,
                      data=DT, verbose = FALSE)
    summary(ansSingle)
    
    # 2. 主效应,假设没有基因型与环境互作
    # 环境作为固定效应+BLUP
    ansMain <- mmer(Yield~Env,
                    random= ~ vs(Name, Gu=A),
                    rcov= ~ units,
                    data=DT, verbose = FALSE)
    summary(ansMain)
    
    # 3. 对角线模型
    # 假设存在基因型与环境互作
    # 假设每个环境基因型的方差不同
    # 每个环境单独拟合基因型效应
    # 缺点:没有环境间的协方差(每个环境为不同个体)
    # 环境作为固定效应+单环境的BLUP(环境的对角线矩阵)
    ansDG <- mmer(Yield~Env,
                  random= ~ vs(ds(Env),Name, Gu=A),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansDG)
    
    # 4. 复合对称模型
    # 假设存在基因型与环境互作
    # 主要的基因型方差协方差分量跨所有环境
    # 假设主要的基因型效应被跨环境效应影响
    # 缺点:各环境之间有相同的方差-协方差
    # 环境作为固定效应+BLUP+GEI
    E <- diag(length(unique(DT$Env)))   # 环境的对角线矩阵
    rownames(E) <- colnames(E) <- unique(DT$Env)   #  添加环境名称
    # 克罗内克乘积,矩阵1的每一项×矩阵2的每一项
    # 若矩阵1是m×n,矩阵2是p×q,则克罗内克乘积是mp×nq的方形矩阵
    # make.dimnames表示添加行名列名
    EA <- kronecker(E,A, make.dimnames = TRUE)
    ansCS <- mmer(Yield~Env,
                  random= ~ vs(Name, Gu=A) + vs(Env:Name, Gu=EA),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansCS)
    
    # 5. 复合对称+对角线模型
    # 假设存在基因型与环境互作
    # 假设每个环境都有不同的基因型与环境互作方差
    # 缺点:各环境有相同的协方差
    # 环境作为固定效应+BLUP+每个地点的基因型与环境互作效应
    ansMain <- mmer(Yield~Env,
                    random= ~ vs(Name, Gu=A) + vs(ds(Env),Name, Gu=A),
                    rcov= ~ units,
                    data=DT, verbose = FALSE)
    summary(ansMain)
    
    # 6. 非结构化模型
    # 假设存在基因型与环境互作
    # 非结构化模型可以得到尽可能多的环境到环境组合的协方差
    # 缺点:方差分量很大,方差组分很难收敛,一些方差协方差组分是零,且很难有良好的初始值
    # 环境作为固定效应+单环境BLUP(协方差调整)
    ansUS <- mmer(Yield~Env,
                  random= ~ vs(us(Env),Name, Gu=A),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansUS)
    
    # 7. 随机回归模型
    # 假设环境可以视为连续变量,因此,截距和斜率的方差组分可以被拟合
    # 方差分量的数量将依赖于Legendre多项式拟合的阶数(leg函数)
    GBIT.library("orthopolynom")
    DT$EnvN <- as.numeric(as.factor(DT$Env))
    ansRR <- mmer(Yield~Env,
                  random= ~ vs(leg(EnvN,1),Name),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansRR)
    # 此外,其他结构的方差协方差矩阵可以应用于这个多项式
    GBIT.library("orthopolynom")
    DT$EnvN <- as.numeric(as.factor(DT$Env))
    ansRR <- mmer(Yield~Env,
                  random= ~ vs(us(leg(EnvN,1)),Name),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansRR)
    
    # 8. 基因型环境互作协方差结构
    # 不常用,但1阶自回归和其他协方差结构可以应用于GEI模型
    # 
    E <- AR1(DT$Env) # can be AR1() or CS(), etc.
    rownames(E) <- colnames(E) <- unique(DT$Env)
    EA <- kronecker(E,A, make.dimnames = TRUE)
    ansCS <- mmer(Yield~Env,
                  random= ~ vs(Name, Gu=A) + vs(Env:Name, Gu=EA),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
    summary(ansCS)
    

    sommer使用直接反演(DI)算法,大数据集会非常缓慢。该包比较适合p>n,即随机效应水平大于材料数。1000个材料,2个重复,2000个记录,或者300个材料,100000个标记,sommer包能快速完成。

  • R语言计算方差协方差矩阵

    方差协方差矩阵,多组变量组成的以方差和协方差数值得到的矩阵。要获得方差-协方差矩阵,需要先知道什么是方差。

    方差:是衡量随机变量或一组数据时离散程度的度量,即(样本中的每个值减去均值)平方和/自由度。

    例如,有一组数据,用R手动计算其方差。

    a <- c(2, 2, 1, 5, 4, 1, 2, 5, 2, 2)

    sum((a-mean(a))^2)/(length(a)-1) #手动计算方差

    var(a) #函数计算方差

    可以对比一下,函数计算的方差与手动计算的方差是否一致。

    协方差:是用于衡量两个随机变量的联合变化程度。而方差是协方差的一种特殊情况,即变量与自身的协方差。

    如果我们又有一组数据,用R计算一下协方差。

    b <- c(2, 2, 1, 4, 4, 3, 5, 1, 3, 5)

    sum((a-mean(a))*(b-mean(b)))/(length(a)-1)

    var(a,b)

    可自行查看结果是否一致。

    方差-协方差矩阵,是用不同数据得到的协方差矩阵。

    如果要得到a,b的方差-协方差矩阵,直接使用var()就可以。

    matrixNoHeader <- cbind(a,b)

    var(matrixNoHeader)

    协方差矩阵行列数量相等,等于变量的个数。

    如果是多个数据

    a <- c(2, 2, 1, 5, 4, 1, 2, 5, 2, 2)
    b <- c(2, 2, 1, 4, 4, 3, 5, 1, 3, 5)
    c <- c(1, 5, 5, 5, 3, 3, 2, 1, 4, 2)
    d <- c(4, 2, 5, 1, 5, 5, 4, 3, 5, 1)

    matrixNoHeader <- cbind(a,b,c,d)
    var(matrixNoHeader)

    对角线是每个变量自身的方差,其余部分是协方差。方差协方差矩阵就是这么简单。

  • Win10子系统安装最新版本R语言【失败】

    直接安装R语言一直是老大难的问题,Ubuntu等系统一般会内置R,但是版本较低,不能满足使用需要,升级通常也是失败的。因此,尝试升级或重新安装最新版本的R语言。经过一番折腾,发现此方法行不通。

    替代方案:

    1. 通过anaconda(https://www.anaconda.com/)安装,很简单,百度一下都有。
    2. 安装Microsoft R Open(https://mran.microsoft.com/open),虽然看起来不是R语言的最新版本,但是稳定性非原版可比,根据官方安装说明操作即可。

    以下是安装失败过程:


    首先,我使用的是Ubuntu 20.04 TLS版本。

    1. 输入下面命令:

    cd /etc/apt/
    sudo vi sources.list

    2. 进入vim编辑模式。接下来,按键盘上的[i]键,进入插入编辑模式。将光标移动到最后一行末尾,回车,插入一行,复制下面的代码,然后按鼠标右键。


    deb https://cloud.r-project.org/bin/linux/ubuntu focal-cran40/
    # deb-src https://cloud.r-project.org/bin/linux/ubuntu focal-cran40/

    3. 然后按[ESC]键退出编辑模式,依次输入[:wq!],回车,即可保存退出。

    4. 接下来,输入下面代码升级源信息。


    sudo apt-get update

    5. 若出现如下错误,则需要获得key。

    Err:2 https://cran.itam.mx/bin/linux/ubuntu groovy-cran40/ InRelease

    The following signatures couldn’t be verified because the public key is not available: NO_PUBKEY 51716619E084DAB9

    6. 使用下面代码获得key,注意替换后面的一串数字。


    sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 51716619E084DAB9

    完成后,重复步骤4。

    7. 接下来安装最新版本的R。


    sudo apt-get install r-base

    8. 若出现下列错误,则按照提示操作。

    E: Unmet dependencies. Try ‘apt –fix-broken install’ with no packages (or specify a solution).

    9. 根据提示,输入下列代码。


    sudo apt --fix-broken install

    10. 重复第7步,若仍然安不上,则先卸载当前R。


    sudo apt-get autoremove r-base-core

  • HMP转Flapjack基因型和图谱文件

    【功能】

    将HMP文件转换为Flapjack的genotype和map文件。

    【测试】

    Null

    【更新历史】

    2020年2月23日

    • 自动完成转换,并在输入文件相同目录下生成【flapjackGeno.txt】文件和【flapjackmap.txt】。

    【使用方法】

    在Rstudio中运行下列代码,选择要转换的hmp文件。


    # Hmp converted to flapjack genotype file for import
    rm(list=ls())
    
    ### function section ###
    # the function of Read files
    readFiles <- function(header=TRUE,choose=FALSE,fname){
      if(!"readr" %in% installed.packages()) {
        install.packages("readr")
        library(readr)
      }else{
        library(readr)
      }
      if (choose==TRUE){
        myFileName <- file.choose()
      }else{
        myFileName <- fname
      }
      if (grepl("\\.csv$",myFileName)){
        if (header==TRUE){
          myFile <- read_csv(myFileName, col_names = TRUE)
        }else{
          myFile <- read_csv(myFileName, col_names = FALSE)
        }
      }else{
        if (header==TRUE){
          myFile <- read_tsv(myFileName, col_names = TRUE)
        }else{
          myFile <- read_tsv(myFileName, col_names = FALSE)
        }
      }
      return(myFile)
    }
    
    # the function of replace the letters of hybrids
    repHybrids <- function(myVector){
      myVector <- sub("R","A/G",myVector)
      myVector <- sub("Y","C/T",myVector)
      myVector <- sub("S","C/G",myVector)
      myVector <- sub("W","A/T",myVector)
      myVector <- sub("K","G/T",myVector)
      myVector <- sub("M","A/C",myVector)
      myVector <- sub("N","-",myVector)
      return(myVector)
    }
    ### function section end ###
    
    ### read file ###  
    # choose a file of HMP
    fn <- choose.files()
    
    # set workspace
    setwd(dirname(fn))
    
    # read a file of HMP
    hmpFile <- readFiles(header=T,fname(basename(fn)))
    
    ### read file end ###
    
    ### genotype ###
    # remove useless columns
    myHmpFile <- as.matrix(hmpFile[,c(1,12:ncol(hmpFile))])
    
    # replace the hybrids
    system.time(newHmpFile <- apply(myHmpFile,2,repHybrids))
    
    # check consistancy
    all(myHmpFile[,1] == newHmpFile[,1])
    all(colnames(myHmpFile) == colnames(newHmpFile))
    
    # Transposing
    resultImprotFile <- t(newHmpFile)
    resultImprotFile <- cbind(rownames(resultImprotFile),resultImprotFile)
    resultImprotFile[,1] <- sub("rs#","",rownames(resultImprotFile))
    resultImprotFile <- rbind("",resultImprotFile)
    resultImprotFile[1,1] <- "# fjFile = GENOTYPE"
    
    # output file
    write.table(resultImprotFile,"flapjackGeno.txt",quote = FALSE,sep="\t",row.names=FALSE,col.names = FALSE)
    ### genotype end ###
    
    ### map ###
    # map file
    myMapFile <- as.matrix(hmpFile[,c(1,3,4)])
    myMapFile[,3] <- as.numeric(myMapFile[,3])
    myMapFile <- rbind("",myMapFile)
    myMapFile[1,1] <- "# fjFile = MAP"
    # output file
    write.table(myMapFile,"flapjackMap.txt",quote = FALSE,sep="\t",row.names=FALSE,col.names = FALSE)
    ### map end ###