集団構造解析ソフトAdmixtureの使い方

形態などのデータからそれなりに多様性があることがわかっている個体群について、各個体の全ゲノム的な多型データが得られているとする。

このとき、 1. その集団を遺伝的に同じ起源を持つグループ（＝祖先集団）に分けた場合、どのように分けられるか？ 2. 祖先集団の交雑により成立した個体はあるか？

を推定する解析を 集団構造解析 と呼ぶ（あくまで自分の分野の場合）。

集団構造解析に用いられるツールとしては、ひと昔前にはSTRUCTUREがもっぱら使われていたが、近年はより計算が早いとされるAdmixtureによる解析が主流になっている。

自分の研究でも活用させてもらったが、データを読み込む際につまずいたので備忘録的に使い方（の注意点）を書いておく。

インストール

公式サイトからDLして展開する。残念ながらWindows版はない模様。admixtureのアイコンをクリックするとターミナルが開くので、あとはマニュアル通りで基本的には問題ない。

ファイルの準備

ただし、読み込みファイルの形式には注意が必要。マニュアルには、.bed, .ped, .genoの3種類の形式が読めると書いているが、.pedもしくは .genoファイルを使いたい場合は、同じディレクトリに対応する.mapファイルを置いておく必要がある。また、.pedファイルは見本通りにエクセルとテキストエディタで作成しても認識されず、必ずPLINKというソフトフェアで保存したものでなければならない。

解析ファイルの作り方

解析ファイルの作成にはいったんTASSELを挟むとうまくいった。

1.PLINKをインストールして任意のフォルダに展開。

2.手持ちのジェノタイプデータをエクセル等で作成してhapmap形式（.txt）で保存。形式はこんな感じ

AlleleID	alleles	chrom	pos	strand	assembly#	center	protLSID	assayLSID	panelLSID	QCcode	ind.1	ind.2	ind.3
marker1	G/A	1	1	NA	NA	NA	NA	NA	NA	NA	GG	GG	GG
marker2	C/T	1	2	NA	NA	NA	NA	NA	NA	NA	CC	CC	CC
marker3	A/G	1	3	NA	NA	NA	NA	NA	NA	NA	AA	AA	AA
marker4	G/C	1	4	NA	NA	NA	NA	NA	NA	NA	GG	GG	GG

3.TASSELで保存したファイルを開き、「Save as...」からPLINK形式を選んで保存。.plkファイルができていることを確認。

4.PLINKを起動し、以下のコマンドを打つ

plink --file 入力ファイル名.plk --recode12 --out 出力ファイル名

5.出力したmapファイルと.pedファイルを、Admixtureを展開したのと同じフォルダに入れれば準備完了

解析

解析の流れとしては、

適切な祖先集団数（K）の推定：マニュアル2.1章の通り
指定したKの値でAdmixture解析、Q値出力
各個体のQ値をプロット：Rのbarplotが一般的

という感じ。具体的なコマンドは以下のブログに詳しく解説されている。 https://indo-european.eu/human-ancestry/admixture-ancestry-components-r-statistics-plink-convertf-bed-and-ped-files/

ついでながら以下にggplot2でbarplotを書く場合のスクリプトを置いておく。ご参考まで。

#K=2
q.k2 <- read_csv(file="acc_name.csv") %>>% # read accessions name
  cbind(read_table(file="acc.plk_fmt.2.Q", col_names = c("q1", "q2"))) %>>% # add Q value for each accession
  dplyr::arrange(q1, q2, species)

#K=3
q.k3 <- read_csv(file="acc_name.csv") %>>%
    cbind(read_table(file="acc.plk_fmt.3.Q", col_names = c("q1", "q2", "q3"))) %>>%
    #note: round the values to properly order the accessions
    dplyr::mutate_if(is.numeric, funs(round(., digits=5))) %>>%
    dplyr::arrange(species, q1, q2, q3)  

#barplot by ggplot2
#K=2
    qg.k2 <- q.k2 %>>%
      tidyr::gather(key = "Q_number", value = "Q_portion", -c("species", "accession")) %>>%
      dplyr::mutate(K_value = "K = 2")

    #K=3
    qg.k3 <-  q.k3 %>>%
      tidyr::gather(key = "Q_number", value = "Q_portion", -c("species", "accession")) %>>%
      dplyr::mutate(K_value = "K = 3")

    qg <- rbind(qg.k2, qg.k3)

    qg %>>%
      ggplot(aes(y = Q_portion, x = accession, fill = Q_number)) +
      geom_bar(stat="identity") +
      labs(x="", y="Q value")+
      theme(
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1),
        legend.position = "left") +
      scale_x_discrete(limits = lab.order) +
      scale_fill_manual(values = c("blue3", "forestgreen", "darkorange3"),
                        name = "")+
      facet_wrap(~K_value, nrow = 2, strip.position = "left")