admixture只有Linux版本，若是想在windows下使用，建議使用R包LEA，功能相似。詳見以前的一篇博文：能夠作structure的R語言包：LEA，本次推送的次條也介紹了這個軟件包。
軟件介紹
基因組選擇中， 有時候測量了不少家系，若是想看一下這些家系的分類狀況，可使用軟件對其進行分羣。通常使用的軟件就是STRUCTURE，可是STREUTURE運行速度極慢，admixture憑藉其運算速度，成爲了主流的分析軟件。下面介紹一下admixture的使用方法。javascript
官方網址
Admixturecss
http://software.genetics.ucla.edu/admixture/download.html
html
軟件安裝
使用conda進行軟件安裝.java
conda install admixture
安裝完成以後, 鍵入admixture, 顯示以下信息, 說明安裝成功nginx
(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. 快速起步
1.1 下載示例數據
wget http://software.genetics.ucla.edu/admixture/hapmap3-files.tar.gz
下載完成以後, 解壓:windows
tar zxvf hapmap3-files.tar.gz
查看解壓後的文件:ruby
(base) [dengfei@localhost admixture]$ lshapmap3.bed hapmap3.bim hapmap3.fam hapmap3-files.tar.gz hapmap3.map
或者在官網上, 下載示例數據: hapmap3-files.tar.gzbash
1.2 admixture支持的格式
plink的bed文件或者ped文件微信
EIGENSTRAT軟件的
.geno格式
注意:markdown若是你的數據格式是plink的bed文件, 好比
a.bed, 那麼你應該包含a.bim,a.fam若是你的數據格式是plink的ped文件, 好比
b.ped, 那麼你應該包括b.map
1.3 選擇合適的分羣數目k值
這裏你要有一個k值, 若是你不知道你的羣體能分爲幾個類羣, 能夠作一個測試, 好比從1~7分別分羣, 而後看他們的cv值哪一個小, 用那個k值.
1.4 運行k=3的admixture
注意, 這裏的名稱爲hapmap3.bed, 而不是hapmap3(不像plink那樣不加後綴), 並且沒有--file 參數, 直接加plink的bed文件
admixture hapmap3.bed 3
運算結果:
[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.
會生成兩個文件:P,Q
hapmap3.3.P hapmap3.3.Q
1.5 運算admixture時, 添加偏差信息
在命令彙總增長一個參數:-B, 速度會變慢喔.
admixture -B hapmap3.bed 3
會生成三個文件:P,Q,Se
1.6 若是你的SNP數據量很大, 跑的很慢
在選擇最佳k值時, 能夠將SNP分爲子集, 好比50k snp分爲50個子集, 每一個子集1k SNP, 那麼根據子集選擇最佳K值, 而後根據最佳的K值去跑全部的SNP
1.7 多線程
若是你有多個線程(processors), 能夠添加參數-jn, n爲線程的個數, 好比你想用4個線程跑:
admixture hapmap3.bed 3 -j 4
2. 參考信息
2.1 如何選擇合適的K值
能夠同時運行多個程序, 每一個程序不一樣的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
這樣跑完以後, 會生成幾個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(交叉驗證的偏差)值:
grep -h CV *.out
[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時, CV error最小
2.2 如何繪製Q的圖表
使用R語言
ta1 = read.table("hapmap3.3.Q")head(ta1)barplot(t(as.matrix(ta1)),col = rainbow(3), xlab = "Individual", ylab = "Ancestry", border = NA)
2.3 我須要根據LD去掉一些SNP麼?
admixture不考慮LD的信息, 若是你想這麼作, 可使用plink
好比, 這裏根據plink 的bed文件進行LD的篩選
plink --bfile hapmap3 --indep-pairwise 50 10 0.1
這裏的過濾參數的意思是:
50, 滑動窗口是50
10, 每次滑動的大小是10
0.1 表示R方小於0.1
而後將其轉化爲bed文件:
plink --bfile hapmap3 --extract plink.prune.in --make-bed --out prunedData
結果輸出過濾後的文件爲:
prunedData.bed prunedData.bim prunedData.fam
使用過濾後的文件, 重新運行admixture:
for K in 1 2 3 4 5 ; do admixture --cv prunedData.bed $K | tee log${K}.out;done
[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)
3. 其它
其它見官方的pdf文檔
若是您對於數據分析，對於軟件操做，對於數據整理，對於結果理解，有任何問題，歡迎聯繫我。
本文分享自微信公衆號 - 育種數據分析之放飛自我（R-breeding）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。