育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

GCTA学习1 | 抛砖引玉--初步介绍

已有 1113 次阅读 2022-1-7 21:13 |个人分类:GWAS|系统分类:科研笔记

康德说:世界上只有两样东西是值得我们深深景仰的,一个是我们头上的灿烂星空,另一个是我们内心的崇高道德法则。

论起软件,也有两款软件,让我学习的过程中,不断的有所收获,不断的感慨功能的强大,一个是plink,另外一个是GCTA。

本篇介绍,是2020年的博客,里面的内容有些老了,最近,GCTA又增加了GWAS的模块fast-GWA,速度远超同类软件,更重要的话,windows电脑也可以用,这样混线性模型GEMMA只能在Linux下使用的弊端,GCTA解决了。不用Linux,就不用Linux吧,现在GCTA也可以分析MLM的GWAS模块了!

1. GCTA介绍

在群体遗传中,GCTA中做PCA非常方便, 下面介绍一下GCTA的安装方法.

2. 安装命令

「使用conda自动安装」

conda install -c biobuilds gcta 

「手动安装」

官方地址

说明文档

3. 安装成功测试

这里, 应该键入gcta64, 而不是gcta

(base) [dengfei@localhost bin]$ gcta64
*******************************************************************
* Genome-wide Complex Trait Analysis (GCTA)
* version 1.26.0
* (C) 2010-2016, The University of Queensland
* MIT License
* Please report bugs to: Jian Yang <jian.yang@uq.edu.au>
*******************************************************************
Analysis started: Wed Apr 24 14:07:43 2019

Options:



Error: no analysis has been launched by the option(s).


Analysis finished: Wed Apr 24 14:07:43 2019
Computational time: 0:0:0

显示上面信息, 表明软件安装成功.

4. 功能介绍

\"在这里插入图片描述\"在这里插入图片描述

5. 参数说明

5.1 输入输出文件

「输入文件:」

  • --bfile test: 类似plink的参数格式. 支持binary文件(test.fam,test.bim,test.bed)
  • --dosage-mach test.mldose test.mlinfo支持其它数据格式

「输出文件:」

  • --out result: 类似plink的--out参数, 定义输出文件名

5.2 数据清洗

「ID保留和删除」如果不写, 默认全部使用

  • --keep test.indi.list定义分析个体列表
  • --remove test.indi.list 删除个体列表

「选择SNP」

--chr 1:选择染色体--autosome 选择所有SNP

6. 构建G矩阵

--make-grm 会生成三个文件:\"在这里插入图片描述\"「如何你想在R中读取二进制文件, 可以使用如下代码:」

ReadGRMBin=function(prefix,AllN=F,size=4){
  
  sum_i=function(i){return(sum(1:i))}
  
  BinFileName=paste(prefix,".grm.bin",sep="")
  NFileName=paste(prefix,".grm.N.bin",sep="")
  IDFileName=paste(prefix,".grm.id",sep="")
  id = read.table(IDFileName) # read the ID of the gmatrix
  n=dim(id)[1]
  BinFile=file(BinFileName,"rb")
  grm=readBin(BinFile,n=n*(n+1)/2,what=numeric(0),size=size)  # generate the fack gmatrix
  NFile=file(NFileName,"rb");
  if(AllN==T){
    N=readBin(NFile,n=n*(n+1)/2,what=numeric(0),size=size)
  }else{
    N=readBin(NFile,n=1,what=numeric(0),size=size)
  }
  i=sapply(1:n,sum_i)
  return(list(diag=grm[i],off=grm[i],id=id,N=N))
}

「计算近交系数」

--ibc: 会用三种方法计算近交系数.

示例:

gcta64 --bfile test --autosome --make-grm --out grm

这里:

  • --bfile 读取的是plink的二进制文件
  • --autosome 是利用常染色体上的所有SNP信息, 这个不能省略
  • --make-grm生成grm矩阵
  • --out 生成前缀名

会生成如下三个文件夹:

(base) [dengfei@localhost plink_file]$ ls grm*
grm.grm.bin  grm.grm.id  grm.grm.N.bin

7. 利用构建好的G矩阵, 计算PCA分析

--grm test: 这里的xx是前缀, 它其实包括三个文件:

test.grm.bin,
test.grm.N.bin
test.grm.id

命令:

gcta64 --grm grm --pca 3 --out out_pca

  • --grmgrm文件
  • --pca PCA的数目为3
  • --out 结果输出文件

结果生成两个文件:

(base) [dengfei@localhost plink_file]$ ls out_pca.eigenv*
out_pca.eigenval  out_pca.eigenvec

8. 利用PCA结果画图

在R语言中, 设置好工作路径, 键入如下命令:

dd=read.table("out_pca.eigenvec",header=F)
head(dd)
names(dd) = c("Fid","ID","PC1","PC2","PC3")
plot(dd$PC1,dd$PC2,pch=c(rep(1,112),rep(2,103)),col=c(rep("blue",112),rep("red",103)),main="PCA",xlab="pc1",ylab="pc2")
legend("bottomright",c("TEXT1","TEXT2"),pch=c(rep(1),rep(2)),col=c(rep("blue"),rep("red")))

结果:在这里插入图片描述

9. GCTA最新版介绍

GCTA最新版是 v1.93.3beta2,它的功能也分为三部分:

9.1 第一部分:遗传力、遗传相关、表型预测

9.2 第二部分:GWAS

在这里插入图片描述在这里插入图片描述

9.3 第三部分:GWAS模拟数据、群体遗传、孟德尔随机化

10. 学习一件事物最好的方法,就是对此写一本书

那就开始吧。

比如,我不太理解GWAS中显著性SNP的解释百分比(PVE),写了几篇博客后,就搞懂了:

GWAS分析中SNP解释百分比PVE | 第一篇,SNP解释百分比之和大于1?

GWAS分析中SNP解释百分比PVE | 第二篇,GLM模型中如何计算PVE?

GWAS分析中SNP解释百分比PVE | 第三篇,MLM模型中如何计算PVE?

GWAS分析中SNP解释百分比PVE | 第四篇,MLM模型中手动计算PVE?

下一章,介绍一下GCTA在Windows和Linux的安装。欢迎关注我,一起学习GCTA软件。

欢迎关注我的公众号:育种数据分析之放飞自我。主要分享R语言,Python,育种数据分析,生物统计,数量遗传学,混合线性模型,GWAS和GS相关的知识。




https://wap.sciencenet.cn/blog-2577109-1320007.html

上一篇:R语言箱线图添加显著性--不同水平实现方法
下一篇:GCTA学习2 | 软件下载安装--windows和Linux

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-6-29 09:29

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部