集団構造の解析を行うにあたってよく利用されるツールとして、STRUCTUREやADMIXTUREがあります。STRUCTUREでなく、ADMIXTUREが使用されることが多くなったが、各々の計算手法は異なることから両方で計算できるに越したことはないであろう、と思いとりあえずやってみたのでSTRUCTURE v2.3.4の解析手法について備忘録的にまとめておきます。
とはいっても、基本的にはipyrad analysis toolsのサイト上にSTRUCTUREに関して書かれているページに書かれた手法をそのままさらっと書くといった内容。
STRUCTUREは、K(集団数)の値を与え、各集団がHardy–Weinberg平衡・連鎖平衡の状態であると仮定して、各集団における対立遺伝子頻度と各個体のクラスター所属(集団の割合)を推定します(Hubisz et al.，2009)。また、クラスターへの割り当ては、マルコフ連鎖モンテカルロ法（MCMC）を用いてパラメータ空間全体を積分しているため、本解析の前にパラメータ空間を絞り込むための計算がある程度の回数必要となります。
この記事ではipyradの出力ファイルである.snp.hdf5をinputファイルとしています。
どうもこのファイルは.vcfからvcf to hdf5 toolなるものを使えば変換できるらしいですが、やってないのでここには書かない。あと、どうも全ゲノムデータの場合はld_block_size引数の指定が必須らしいですが、他のデータセットの場合どうなるかなどはいまいちわかりません(やってないから)。STRUCTUREは長い間使われているのでipyrad analysis tools以外の手法が色々充実しているであろうと思います。ほかの手法をお探し下さい。
環境
macOS 14.2.1
ipyrad 0.9.95
STRUCTURE 2.3.4
解析
ipyrad環境を立ち上げてpythonに移動する。
#python環境に移動
conda activate ipyrad-env
python3
移動したのち、ipyrad.analysisからstructureを呼び出して集団、個体の指定やフィルタリングの設定を行う。
これらを指定した後、探索空間を絞り込むためのMCMCの回数と本解析の回数を設定し、struct run()でstructureを動かす。
from ipyrad.analysis import structure
#集団ラベルと個体ラベルの指定
imap = {
"pop1": ["ind1", "ind2", "ind3", "ind4"],
"pop2": ["ind5", "ind6", "ind7",],
"pop3": ["ind8", "ind9"],
"pop4": ["ind10", "ind11", "ind12"],
"pop5": ["ind13", "ind14", "ind15", "ind16"],
}
#各集団でとある割合以上の共有がされていないSNPをフィルタリング(0.5で50%)
minmap = {i: 0.5 for i in imap}
#mincovは各個体でとある割合以上の共有がされていないSNPをフィルタリング
#.snps.hdf5のパス
struct = ipa.structure(
workdir="保存/先の/パス"
data="/home/deren/Downloads/ref_pop2.snps.hdf5",
imap=imap,
minmap=minmap,
mincov=0.9,
)
#burninがMCMCの反復回数。足りない場合、結果が収束しない
#numrepが本解析の計算回数。
struct.mainparams.burnin = 2000000
struct.mainparams.numreps = 5000000
#nrepsは出力される結果の数。kpopは集団数
struct.run(nreps=10, kpop=[2, 3, 4, 5, 6, 7, 8], auto=True)
#evanno法で計算されたΔKの値を出力。
etable = struct.get_evanno_table([2, 3, 4, 5])
etable
#出力されたtableを保存
etable.to_csv("保存/先の/パス/choosingK_table.txt", sep="/t", index=FALSE)
というような感じ。
結果の出力
STRUCTURE側でも図の出力はできるようですが、ここではADMIXTUREと同様にRパッケージのpophelperを使用する方法でいきます。
pophelperを使用する方法については、少し調べてみるとHatena blogでSTRUCTURE→pophelperで図の出力までの工程を書いてくださっている記事があるのでそちらも参考にするとよいと思います。
ここではfor構文を使用して、出力されたSTRUCTURE 10回分の出力結果をすべて図にしています。また、 K=2~8 までの図をすべてつなげた図を出したかったので、そのようなスクリプトにしています。
library(pophelper)
# ベースディレクトリの指定
basedir <- "your/data/pass/STRUCTUREの出力ファイルがあるところ"
# rep = 0 〜 9（10個分の出力データ） に対応する図を出力
for (i in 0:9) {
pattern <- paste0("-", i, "_f$")
# 該当する rep のファイルを取得
files_rep <- list.files(
path = basedir,
pattern = pattern,
full.names = TRUE
)
# ファイルが存在する場合のみ処理
if (length(files_rep) > 0) {
# データ読み込み
slist_rep <- readQ(files_rep)
# 図の出力。ここをいじると出力する図について色々いじれる
#sorting = "label"でインプットした順番通りに個体が並ぶ
p <- plotQ(
slist_rep,
imgoutput = "join",
imgtype = "pdf",
exportpath = basedir,
outputfilename = paste0("STRUCTURE_plot_rep", i),
sortind = "label",
sharedindlab = FALSE,
showindlab = FALSE
)
cat("Finished rep", i, "\n")
} else {
cat("No files found for rep", i, "\n")
}
}
以上です。
間違いなどあるような気がしているので、もしよければお教えいただければ嬉しく思います。
引用文献・サイト
Hubisz et al.，2009. Inferring weak population structure with the assistance of
sample group information.Molecular Ecology Resources 9, 1322–1332 https://doi.org/10.1111/j.1755-0998.2009.02591.x
STRUCTUREのサイト：https://web.stanford.edu/group/pritchardlab/structure.html
ipyrad analysis toolsのサイト：https://ipyrad.readthedocs.io/en/master/API-analysis/index.html
Gai memo【SNP 解析】STRUCTUREのきれいな図作成への道 by R：https://nemunemu-nyanko.hatenablog.com/entry/2021/01/27/135253