【原】GCTA PCA分析cookbook

育种数据分析 2021-11-18

展开全文

GCTA软件介绍系列

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")))

结果:

后记1, 使用示例数据`b.ped`和`b.map`使用gcta64做PCA分析

看完gcta, 发现plink也可以构建G矩阵, 也可以进行PCA分析, 本数据使用plink的解决方案:

将ped文件, 转化为bed文件

plink --file b --make-bed --out test

生成test.bed, test.bim,test.fam三个文件

构建G矩阵grm

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

生成三个文件:

grm.grm.bin  grm.grm.id  grm.grm.N.bin

生成PCA, 数目为3

gcta64 --grm grm --pca 3

生成两个文件:

gcta.eigenval  gcta.eigenvec

作图

dd=read.table("gcta.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")))

结果:

后记2, 使用示例数据`b.ped`和`b.map`使用`plink`做PCA分析

看完gcta, 发现plink也可以构建G矩阵, 也可以进行PCA分析, 本数据使用plink的解决方案:

只用一行代码, 就可以生成PCA的数据, 比gcta64简单太多了.

plink --file b --pca 3

比较两个数据的结果, 可以看出, plink和gcta64结果一致.

对PCA作图:

结果一致, 因为plink调用的是gcta64的算法, 构建G矩阵, 构建PCA.

福利1
计算gcta64或者plink可以构建矩阵, asreml也支持下三角的G矩阵或者G逆矩阵, 问题来了, 两者怎么联系到一起呢? 这样asreml就可以愉快的进行GBLUP的分析了.

福利2
之前的博客中有提到利用H矩阵构建PCA分析, 那么如何操作呢?

欲听后事如何, 请听下回分解.

公众号后台回复:plink, 获得测试数据:b.ped和b.map, 用于本次分析.

如果您对于数据分析，对于软件操作，对于数据整理，对于结果理解，有任何问题，欢迎联系我。

作者其它博文:

生物统计：

主要包括试验设计，生物统计中的数据分析，育种中的数据分析，相关的文献解读。

1，用R语言生成增广试验设计

4，如何对增广试验数据进行分析

5，如何对数据进行汇总统计(R语言)

6，关于联合方差分析的讨论-1

7，农业统计分析系列1-软件包介绍

8，农业统计分析系列2-试验设计

9，Excel中的数据透视功能处理农业数据

10，进军机器学习--序言

11, 植物育种中全基因组选择是成熟的方法么?

12, 不同试验设计遗传力的计算方法

13, 农业大数据时代的几个案例

14, 农业试验中如何分析单因素方差分析

15, P-rep designs 文献解析及实现方法

16, 文献阅读：林木中遗传参数评估

17, 育种4.0世代的到来个人应该准备什么

18, 农业试验设计中田间种植图的绘制方法

数量遗传：

主要是动物数量遗传，动物育种中应用比较广泛，无论是基于系谱的动物模型，近交系数，亲缘关系系数，配合力，育种值，还是单性状模型，重复力模型，多性状模型等相关知识。

10, DMU遗传参数评估cookbook pdf

17，单性状动物模型矩阵形式计算BLUP值

18, 文献阅读: ABLUP-GBLUP-SSGBLUP模拟数据比较

19, 文献阅读：林木中遗传参数评估

20, 遗传变异系数怎么计算

21, 测定日模型及随机模型介绍

22, Admixture使用说明文档cookbook

编程语言：

包括Python，R语言，Julia，Perl语言，Linux的Shell语言，主要是我平时学习时的一些笔记和总结。

1，R，Julia以及Python共享数据

2，Python生物统计---笔记1

3，Python学生物统计---笔记2

4，Python学生物统计---笔记3

5，Python学生物统计---数据导入笔记4

6，Python学生物统计---可视化---笔记5

7，Python学生物统计---T检验笔记6

8，Python学生物统计---方差分析笔记7

9，shiny学习笔记1---上传数据

10，shiny学习笔记2-下载数据

11，shiny学习笔记3--生成html报告

12，data.table学习笔记1

13，data.table学习笔记2

14, R语言与独孤九剑以及Python与降龙十八掌

15, snakemake 学习笔记1

16, snakemake 学习笔记2

17, 远程访问服务器 jupyter的设置方法

18, WOX 糙快猛的实现方法

19, R语言中如何写入xlsx的不同sheet表格

20, 几种加快R语言运算的方法

21, 如何批量安装R语言包

22, 如何高效的在服务器和本地进行上传和下载文件

23, 如何优雅的使用markdown写博客

基因组选择：

育种数据分析中，表型选择，方差分析，混合线性模型的BLUP育种值是学科的枝干，MAS，GWAS是花苞, GS则是盛开的花朵，其依赖于常规的数量遗传理论，但青出于蓝而胜于蓝，具有光明的前景，由于GS的应用，分子育种的落地又大大提前了一步。现在GS在动物育种中，特别是牛，猪，鸡，羊中正在大规模落地，以后再玉米，水稻，小麦，大豆的应用也将落地。冬天来了，春天还会远么？这个章节有文献解析，SNP数据清洗，G矩阵及H矩阵构建，模拟数据，软件使用，理论介绍等等。