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つのファイルを合わせて、サンプルの血統情報と遺伝子型、遺伝子座の情報を提供します。
このファイルの中身をテキストエディタで見ると、以下のようになっています。
# 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より引用、改変
# 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でプロット
これが最も一般的な方法だと思われます。
いくつか方法がありますが、例えば以下のサイトが非常に参考になります。
引用文献・サイト
- Alexander, D.H., Lange, K. 2011. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12: 246 https://doi.org/10.1186/1471-2105-12-246
- ADMIXTUREのマニュアル：http://dalexander.github.io/admixture/admixture-manual.pdf
- PLINKのマニュアル：https://www.cog-genomics.org/plink/1.9/
- CLUMPAK：http://clumpak.tau.ac.il/
- StructureSelector：https://lmme.qdio.ac.cn/StructureSelector/