【原】GCTA PCA分析cookbook

育種數(shù)據(jù)分析 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. 參數(shù)說明

5.1 輸入輸出文件

輸入文件:

--bfile test: 類似plink的參數(shù)格式. 支持binary文件(test.fam,test.bim,test.bed)
--dosage-mach test.mldose test.mlinfo支持其它數(shù)據(jù)格式

輸出文件:

--out result: 類似plink的--out參數(shù), 定義輸出文件名

5.2 數(shù)據(jù)清洗

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

計算近交系數(shù)

--ibc: 會用三種方法計算近交系數(shù).

示例:

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的數(shù)目為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, 使用示例數(shù)據(jù)`b.ped`和`b.map`使用gcta64做PCA分析

看完gcta, 發(fā)現(xiàn)plink也可以構建G矩陣, 也可以進行PCA分析, 本數(shù)據(jù)使用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, 數(shù)目為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, 使用示例數(shù)據(jù)`b.ped`和`b.map`使用`plink`做PCA分析

看完gcta, 發(fā)現(xiàn)plink也可以構建G矩陣, 也可以進行PCA分析, 本數(shù)據(jù)使用plink的解決方案:

只用一行代碼, 就可以生成PCA的數(shù)據(jù), 比gcta64簡單太多了.

plink --file b --pca 3

比較兩個數(shù)據(jù)的結果, 可以看出, plink和gcta64結果一致.

對PCA作圖:

結果一致, 因為plink調用的是gcta64的算法, 構建G矩陣, 構建PCA.

福利1
計算gcta64或者plink可以構建矩陣, asreml也支持下三角的G矩陣或者G逆矩陣, 問題來了, 兩者怎么聯(lián)系到一起呢? 這樣asreml就可以愉快的進行GBLUP的分析了.

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

欲聽后事如何, 請聽下回分解.

公眾號后臺回復:plink, 獲得測試數(shù)據(jù):b.ped和b.map, 用于本次分析.

如果您對于數(shù)據(jù)分析，對于軟件操作，對于數(shù)據(jù)整理，對于結果理解，有任何問題，歡迎聯(lián)系我。

作者其它博文:

生物統(tǒng)計：

主要包括試驗設計，生物統(tǒng)計中的數(shù)據(jù)分析，育種中的數(shù)據(jù)分析，相關的文獻解讀。

1，用R語言生成增廣試驗設計

2，P-rep designs 文獻解析及實現(xiàn)方法

4，如何對增廣試驗數(shù)據(jù)進行分析

5，如何對數(shù)據(jù)進行匯總統(tǒng)計(R語言)

6，關于聯(lián)合方差分析的討論-1

7，農業(yè)統(tǒng)計分析系列1-軟件包介紹

8，農業(yè)統(tǒng)計分析系列2-試驗設計

9，Excel中的數(shù)據(jù)透視功能處理農業(yè)數(shù)據(jù)

10，進軍機器學習--序言

11, 植物育種中全基因組選擇是成熟的方法么?

12, 不同試驗設計遺傳力的計算方法

13, 農業(yè)大數(shù)據(jù)時代的幾個案例

14, 農業(yè)試驗中如何分析單因素方差分析

15, P-rep designs 文獻解析及實現(xiàn)方法

16, 文獻閱讀：林木中遺傳參數(shù)評估

17, 育種4.0世代的到來個人應該準備什么

18, 農業(yè)試驗設計中田間種植圖的繪制方法

數(shù)量遺傳：

主要是動物數(shù)量遺傳，動物育種中應用比較廣泛，無論是基于系譜的動物模型，近交系數(shù)，親緣關系系數(shù)，配合力，育種值，還是單性狀模型，重復力模型，多性狀模型等相關知識。

1，R語言求解混合線性方程組（有系譜）

2，R語言混合線性模型包代碼演示

3，DMU-遺傳參數(shù)評估-學習筆記1

10, DMU遺傳參數(shù)評估cookbook pdf

11，育種中一般和特殊配合力的計算方法

12，為什么要學習數(shù)量遺傳學1--序言

13，2-數(shù)量遺傳學課程介紹

14，3-數(shù)量遺傳學課程介紹-R語言基礎

15，4-數(shù)量遺傳學課程介紹-R挖掘數(shù)據(jù)

16，利用系譜計算近交親緣關系系數(shù)

17，單性狀動物模型矩陣形式計算BLUP值

18, 文獻閱讀: ABLUP-GBLUP-SSGBLUP模擬數(shù)據(jù)比較

19, 文獻閱讀：林木中遺傳參數(shù)評估

20, 遺傳變異系數(shù)怎么計算

21, 測定日模型及隨機模型介紹

22, Admixture使用說明文檔cookbook

編程語言：

包括Python，R語言，Julia，Perl語言，Linux的Shell語言，主要是我平時學習時的一些筆記和總結。

1，R，Julia以及Python共享數(shù)據(jù)

2，Python生物統(tǒng)計---筆記1

3，Python學生物統(tǒng)計---筆記2

4，Python學生物統(tǒng)計---筆記3

5，Python學生物統(tǒng)計---數(shù)據(jù)導入筆記4

6，Python學生物統(tǒng)計---可視化---筆記5

7，Python學生物統(tǒng)計---T檢驗筆記6

8，Python學生物統(tǒng)計---方差分析筆記7

9，shiny學習筆記1---上傳數(shù)據(jù)

10，shiny學習筆記2-下載數(shù)據(jù)

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 糙快猛的實現(xiàn)方法

19, R語言中如何寫入xlsx的不同sheet表格

20, 幾種加快R語言運算的方法

21, 如何批量安裝R語言包

22, 如何高效的在服務器和本地進行上傳和下載文件

23, 如何優(yōu)雅的使用markdown寫博客

基因組選擇：

育種數(shù)據(jù)分析中，表型選擇，方差分析，混合線性模型的BLUP育種值是學科的枝干，MAS，GWAS是花苞, GS則是盛開的花朵，其依賴于常規(guī)的數(shù)量遺傳理論，但青出于藍而勝于藍，具有光明的前景，由于GS的應用，分子育種的落地又大大提前了一步?，F(xiàn)在GS在動物育種中，特別是牛，豬，雞，羊中正在大規(guī)模落地，以后再玉米，水稻，小麥，大豆的應用也將落地。冬天來了，春天還會遠么？這個章節(jié)有文獻解析，SNP數(shù)據(jù)清洗，G矩陣及H矩陣構建，模擬數(shù)據(jù)，軟件使用，理論介紹等等。