LoginSignup
3
2

More than 1 year has passed since last update.

@YF_bio

populationsから出力したファイルでADMIXTURE解析をする

Last updated at Posted at 2020-12-24

RAD-seqやMIG-seqでシーケンスした多型データから集団構造を調べたい際に、大抵の人はクラスタリング解析をすると思います。

ADMIXTUREはSTRUCTUREと比較して計算が非常に高速なのが特徴で、多くの論文で使われている手法です。SNP検出ソフトであるStacksから出力したファイルは、PLINKで編集しないとADMIXTUREで解析できないため、その方法を備忘録的に残しておきます。

環境

macOS 11.1
Stacks 2.54
PLINK 1.90b6.21
ADMIXTURE 1.3.0

ADMIXTUREの導入

以下のサイトからダウンロードします。
Mac OSとLinuxに対応しています。コンパイル等は必要ありません。
M1チップ搭載のMacの場合、ターミナルをRosettaを使用して開くように設定しておきます。

ディレクトリを指定せずにADMIXTUREを起動するために、実行ファイルをcp /Downloads/admixture_macosx-1.3.0/admixture /usr/local/binでコピーしてパスを通しておきます。

populationsからの出力

Stacksのプログラムの一つであるpopulationsはカタログ化された遺伝子座から、Popmapファイルで指定された個体ごとに一定の基準で遺伝子座およびSNPを出力してくれます。
Stacksの使い方については、以下をご参考ください。

populationsでのジェノタイプ出力の際に、PLINK形式のファイルを出力するオプション**--plink**を付けます。

ターミナル
populations -P ./stacks -M ./popmap.txt -O ./populations --write-single-snp --plink -R 0.6 --min-mac 2 --max-obs-het 0.75 -t 16

これで以下の2つのファイルが出力されます。

populations.plink.ped
populations.plink.map

2つのファイルを合わせて、サンプルの血統情報と遺伝子型、遺伝子座の情報を提供します。
このファイルの中身をテキストエディタで見ると、以下のようになっています。

populations.plink.ped
# Stacks v2.54;  PLINK v1.07; December 18, 2020
Lake	Channa_sp_1	0	0	0	0	C	C	T	T	A...
Lake	Channa_sp_2	0	0	0	0	C	C	T	T	A...
River	Channa_sp_3	0	0	0	0	C	C	C	T	A...
River	Channa_sp_4	0	0	0	0	C	C	C	T	A...
...

1列目　pop ID
2列目　individual ID
3列目　父親のID（不明の場合は0）
4列目　母親のID（不明の場合は0）
5列目　性別（1 = male, 2 = female, 0 = unknown）
6列目　形質値（1 = control, 2 = case, -9/0 = missing）

File format referenceより引用、改変

populations.plink.map
# Stacks v2.54;  PLINK v1.07; December 18, 2020
un	1_6	0	7
un	2_51	0	150
un	20_27	0	533
un	28_17	0	905
un	30_1	0	1063
un	36_25	0	1262
un	38_14	0	1433
...

1列目　染色体番号
2列目　変異サイトの識別子
3列目　遺伝子間距離
4列目　変異サイトの座標
File format referenceより引用、改変

このうち、MAPファイルの染色体番号がunだとADMIXTUREで解析ができないので、PLINKでいじります。

PLINKでの変換

以下のコマンドでファイルを変換します。

ターミナル
plink --file populations.plink --allow-extra-chr 0 --recode 12 --missing -out Admixture-input

#--fileでpopulationsから出力したファイルを指定
#--allow-extra-chr 0 で染色体番号を0として扱う
#PED形式で出力する場合：--recode 12
#BED形式で出力する場合：--make-bed
#--missingで個体ごと、遺伝子座ごとにミッシングの数を計算
#-outで出力ファイルの名前を指定

処理は秒で終わります。

ADMIXTUREでの解析

PLINKで変換したファイルをインプットとしてADMIXTURE解析行います。
基本のコマンドは以下の通りです。

ターミナル
admixture -C 15 -s time --cv Admixture-input.ped 3 log3.out

#インプットファイル、解析するK値、ログファイル名を指定します
#-Cは解析の終了基準を指定します
#-Cが正の整数の場合は繰り返し回数、Cが1未満の場合は対数尤度の増加率が指定の数値を下回るまで繰り返すように設定できます
#-sはシード値です。timeにしておけば現在時刻をシード値として解析を始めます。
#--cvは各KにおいてCross-validation error値を算出してくれます

これだとKの値ごとに回さないといけないため、大変面倒くさいです。
そこでfor文で予め指定したKについて、自動で解析してもらいます。

ターミナル
for K in {1..12}; do admixture -C 0.0001 -s time --cv Admixture-input.ped $K | tee log${K}.out; done

解析が終わるとK値ごとに.Qファイル、.Pファイル、log.outファイルが生成されています。
図示の際にはQファイルを使用します。

最適クラスター数Kを選ぶ基準となるCross-validationは各解析のログファイルに残っています。
いちいち確認するのは面倒なので、以下のコマンドでまとめて表示しテキストファイルに保存します。

ターミナル
grep -h CV log*.out > CV-error.txt
cat CV-error.txt

CV error (K=2): 0.27856
CV error (K=3): 0.19952
CV error (K=4): 0.20640
CV error (K=5): 0.18858
CV error (K=6): 0.20056
CV error (K=7): 0.19954
CV error (K=8): 0.19560
CV error (K=9): 0.20710
CV error (K=10): 0.21673
CV error (K=11): 0.20096
CV error (K=12): 0.22981

結果の可視化

便利なウェブツールがいくつもあるので、簡単に結果を可視化できます。

CLUMPAK (URL)

CLUMPAKは解析したすべてのK値について、まとめて図示してくれる超便利ウェブツールです。
ただし、論文で使えるような綺麗な図は出してくれないので、あくまで予備的なものです。

まず、ADMIXTUREから出力された.Qファイルをフォルダにまとめて圧縮します。

ターミナル
mkdir Qmatrix ; cp *.Q Qmatrix/ ; zip -r Qmatrix.zip Qmatrix/

#Qmatrizフォルダを作成し、拡張子が.Qのファイルをそこに移し、それを圧縮しています。

圧縮したQファイルと1列だけのPopmapをウェブサイトでポン入れしてSubmitすると、図が出力されます。
ちょっと時間がかかります。

StructureSelector (URL)

StructureSelectorはcross-validation error以外の指標によって最適なK値を推定してくれるウェブツールです。
推定方法については詳しくないのでウェブサイトか論文を参照してください。

こちらはADMIXTUREから出力されたすべてのファイルをまとめて圧縮したものを使います。
他はCLUMPAKと同様です。

Rでプロット

これが最も一般的な方法だと思われます。
いくつか方法がありますが、例えば以下のサイトが非常に参考になります。

引用文献・サイト

3
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
2