【原】群體遺傳三劍客第二篇：Admixture群體結(jié)構(gòu)分析

育種數(shù)據(jù)分析 2025-03-17 發(fā)布于河南

展開全文

大家好，我是鄧飛。

之前計(jì)劃寫群體結(jié)構(gòu)三劍客的博文，寫了兩篇了：

搞起來(lái)！群體遺傳三劍客：PCA、Admixture、進(jìn)化樹

群體遺傳三劍客第一篇：分組和不分組的PCA分析，添加解釋百分比

今天介紹一下群體結(jié)構(gòu)分析最常用的軟件：admixture。注意，admixture需要在Linux中使用，不能在windows系統(tǒng)，windows的子系統(tǒng)Linux也是不能使用的。

軟件介紹

群體結(jié)構(gòu)分析，一般使用的軟件就是STRUCTURE，但是STREUTURE運(yùn)行速度極慢，admixture憑借其運(yùn)算速度，成為了主流的分析軟件。下面介紹一下admixture的使用方法。

官方網(wǎng)址：

http://dalexander./admixture/download.html

軟件不太好下載，官網(wǎng)不太好打開，官網(wǎng)上面一共四個(gè)文件：

- admixture的Linux版本

- admixture的Mac版本

- admixture的pdf文檔

- admixture的示例數(shù)據(jù)

公眾號(hào)回復(fù)：admixture，獲得上面四個(gè)文件的下載鏈接。

另外，還可以使用conda進(jìn)行軟件安裝.

conda install admixture

安裝完成之后, 鍵入admixture, 顯示如下信息, 說(shuō)明安裝成功

(base) [dengfei@localhost test]$ admixture****                   ADMIXTURE Version 1.3.0                  ********                    Copyright 2008-2015                     ********           David Alexander, Suyash Shringarpure,            ********                John  Novembre, Ken Lange                   ********                                                            ********                 Please cite our paper!                     ********   Information at www.genetics.ucla.edu/software/admixture  ****
Usage: admixture <input file> <K>See --help or manual for more advanced usage.

1，示例數(shù)據(jù)解壓

下載之后, 上傳到Linux系統(tǒng)，然后解壓:

tar zxvf hapmap3-files.tar.gz

查看解壓后的文件:

(base) [dengfei@localhost admixture]$ lshapmap3.bed  hapmap3.bim  hapmap3.fam  hapmap3-files.tar.gz  hapmap3.map

2 admixture支持的格式

plink的bed文件或者ped文件
EIGENSTRAT軟件的.geno格式
注意:
如果你的數(shù)據(jù)格式是plink的bed文件, 比如a.bed, 那么你應(yīng)該包含a.bim, a.fam
如果你的數(shù)據(jù)格式是plink的ped文件, 比如b.ped, 那么你應(yīng)該包括b.map

3 選擇合適的分群數(shù)目k值

這里你要有一個(gè)k值, 如果你不知道你的群體能分為幾個(gè)類群, 可以做一個(gè)測(cè)試, 比如從1~7分別分群, 然后看他們的cv值哪個(gè)小，用哪個(gè)k值。

4 運(yùn)行k=3的admixture

注意, 這里的名稱為hapmap3.bed, 而不是hapmap3(不像plink那樣不加后綴), 而且沒(méi)有--file 參數(shù), 直接加plink的bed文件

admixture hapmap3.bed 3

運(yùn)算結(jié)果:

(base) [dengfei@localhost admixture]$ admixture hapmap3.bed 3****                   ADMIXTURE Version 1.3.0                  ********                    Copyright 2008-2015                     ********           David Alexander, Suyash Shringarpure,            ********                John  Novembre, Ken Lange                   ********                                                            ********                 Please cite our paper!                     ********   Information at www.genetics.ucla.edu/software/admixture  ****
Random seed: 43Point estimation method: Block relaxation algorithmConvergence acceleration algorithm: QuasiNewton, 3 secant conditionsPoint estimation will terminate when objective function delta < 0.0001Estimation of standard errors disabled; will compute point estimates only.Size of G: 324x13928Performing five EM steps to prime main algorithm1 (EM)     Elapsed: 0.318    Loglikelihood: -4.38757e+06    (delta): 2.87325e+062 (EM)     Elapsed: 0.292    Loglikelihood: -4.25681e+06    (delta): 1307623 (EM)     Elapsed: 0.29    Loglikelihood: -4.21622e+06    (delta): 40582.94 (EM)     Elapsed: 0.29    Loglikelihood: -4.19347e+06    (delta): 22748.25 (EM)     Elapsed: 0.29    Loglikelihood: -4.17881e+06    (delta): 14663.1Initial loglikelihood: -4.17881e+06Starting main algorithm1 (QN/Block)     Elapsed: 0.741    Loglikelihood: -3.94775e+06    (delta): 2310582 (QN/Block)     Elapsed: 0.74    Loglikelihood: -3.8802e+06    (delta): 67554.63 (QN/Block)     Elapsed: 0.852    Loglikelihood: -3.83232e+06    (delta): 47883.84 (QN/Block)     Elapsed: 1.01    Loglikelihood: -3.81118e+06    (delta): 21138.25 (QN/Block)     Elapsed: 0.903    Loglikelihood: -3.80682e+06    (delta): 4354.366 (QN/Block)     Elapsed: 0.85    Loglikelihood: -3.80474e+06    (delta): 2085.657 (QN/Block)     Elapsed: 0.856    Loglikelihood: -3.80362e+06    (delta): 1112.588 (QN/Block)     Elapsed: 0.908    Loglikelihood: -3.80276e+06    (delta): 865.019 (QN/Block)     Elapsed: 0.852    Loglikelihood: -3.80209e+06    (delta): 666.66210 (QN/Block)     Elapsed: 1.015    Loglikelihood: -3.80151e+06    (delta): 579.4911 (QN/Block)     Elapsed: 0.908    Loglikelihood: -3.80097e+06    (delta): 548.15612 (QN/Block)     Elapsed: 0.961    Loglikelihood: -3.80049e+06    (delta): 473.56513 (QN/Block)     Elapsed: 0.855    Loglikelihood: -3.80023e+06    (delta): 258.6114 (QN/Block)     Elapsed: 0.959    Loglikelihood: -3.80005e+06    (delta): 179.94915 (QN/Block)     Elapsed: 1.011    Loglikelihood: -3.79991e+06    (delta): 146.70716 (QN/Block)     Elapsed: 0.903    Loglikelihood: -3.79989e+06    (delta): 13.194217 (QN/Block)     Elapsed: 1.01    Loglikelihood: -3.79989e+06    (delta): 4.6074718 (QN/Block)     Elapsed: 0.85    Loglikelihood: -3.79989e+06    (delta): 1.5001219 (QN/Block)     Elapsed: 0.851    Loglikelihood: -3.79989e+06    (delta): 0.12891620 (QN/Block)     Elapsed: 0.851    Loglikelihood: -3.79989e+06    (delta): 0.0018298321 (QN/Block)     Elapsed: 0.851    Loglikelihood: -3.79989e+06    (delta): 4.33805e-05Summary:Converged in 21 iterations (21.788 sec)Loglikelihood: -3799887.171935Fst divergences between estimated populations:    Pop0    Pop1    Pop0    Pop1    0.163    Pop2    0.073    0.156    Writing output files.

會(huì)生成兩個(gè)文件:P,Q

hapmap3.3.P  hapmap3.3.Q

5 運(yùn)算admixture時(shí), 添加誤差信息

在命令匯總增加一個(gè)參數(shù):-B, 速度會(huì)變慢.

admixture -B hapmap3.bed 3

會(huì)生成三個(gè)文件:P,Q,Se

6 如果你的SNP數(shù)據(jù)量很大, 跑的很慢

在選擇最佳k值時(shí), 可以將SNP分為子集, 比如50k snp分為50個(gè)子集, 每個(gè)子集1k SNP, 那么根據(jù)子集選擇最佳K值, 然后根據(jù)最佳的K值去跑所有的SNP

7 多線程

如果你有多個(gè)線程(processors), 可以添加參數(shù)-jn, n為線程的個(gè)數(shù), 比如你想用4個(gè)線程跑:

admixture  hapmap3.bed 3 -j 4

8 如何選擇合適的K值

可以同時(shí)運(yùn)行多個(gè)程序, 每個(gè)程序不同的k值, 比如, 想要k值選擇1,2,3,4,5, 可以寫為:

 for K in 1 2 3 4 5; do admixture --cv hapmap3.bed $K | tee log${K}.out; done

這樣跑完之后, 會(huì)生成幾個(gè)out文件,

hapmap3.1.P  hapmap3.1.Q  hapmap3.2.P  hapmap3.2.Q  hapmap3.3.P  hapmap3.3.Q  hapmap3.4.P  hapmap3.4.Q  hapmap3.5.P  hapmap3.5.Q log1.out  log2.out  log3.out  log4.out  log5.out

使用grep查看*out文件的cv error(交叉驗(yàn)證的誤差)值:

grep -h CV  *.out

(base) [dengfei@localhost admixture]$ grep -h CV *outCV error (K=1): 0.55248CV error (K=2): 0.48190CV error (K=3): 0.47835CV error (K=4): 0.48236CV error (K=5): 0.49001

可以看出, K=3時(shí), CV error最小

9 如何繪制Q的圖表

使用R語(yǔ)言

ta1 = read.table("hapmap3.3.Q")head(ta1)barplot(t(as.matrix(ta1)),col = rainbow(3),        xlab = "Individual",        ylab = "Ancestry",        border = NA)

10 我需要根據(jù)LD去掉一些SNP么?

admixture不考慮LD的信息, 如果你想這么做, 可以使用plink

比如, 這里根據(jù)plink 的bed文件進(jìn)行LD的篩選

plink  --bfile hapmap3 --indep-pairwise 50 10 0.1

這里的過(guò)濾參數(shù)的意思是:

50, 滑動(dòng)窗口是50
10, 每次滑動(dòng)的大小是10
0.1 表示R方小于0.1

然后將其轉(zhuǎn)化為bed文件:

plink  --bfile hapmap3 --extract plink.prune.in --make-bed --out prunedData

結(jié)果輸出過(guò)濾后的文件為:

prunedData.bed  prunedData.bim  prunedData.fam

使用過(guò)濾后的文件, 從新運(yùn)行admixture:

for K in 1 2 3 4 5 ; do admixture --cv prunedData.bed $K | tee log${K}.out;done

(base) [dengfei@localhost ld-test]$ grep -h CV *outCV error (K=1): 0.52305CV error (K=2): 0.48847CV error (K=3): 0.48509CV error (K=4): 0.49404CV error (K=5): 0.49828

可以看出K=3, cv error最小, 因此選擇k=3

作圖:

ta1 = read.table("prunedData.3.Q")head(ta1)barplot(t(as.matrix(ta1)),col = rainbow(3),        xlab = "Individual",        ylab = "Ancestry",        border = NA)