IMa2 (Isolation with Migration analysis 2)
Hey J. 2010. Isolation with migration models for more than two populations. Mol Biol Evol 27:905–920.
1. ソフトウエアの入手先
https://bio.cst.temple.edu/~hey/software/software.htm
2. Infileの準備
example fileに従って作成する（末尾に別記）。
3. プログラムの起動
プログラムを開始するにはコマンドプロンプトを立ち上げ、IMa2のフォルダを指定する必要がある。
たとえば、C:Program files\IMa2にIMa2.exeといった具合。
4. コマンド
IMa2では、プログラムの実行に際して、一行の命令文を入力する。
命令文には、集団サイズ、移入率、分岐年代の探索範囲、Heatingの具合、どのパラメータについて結果を表示するか等を記述する。
Ex) IMa2.exe -i infile.txt -o outfile.txt -q1000 -m0.1 -t30 -b100000 -l50000 -d100 -s123 -p2 -p7 -hfg -hn40 -ha0.975 -hb0.75
この場合、
IMa2.exe （IMa2.exeを以下の条件で実行する）
-i infile.txt （インファイルは同じフォルダ内のinfile.txt）
-o outfile.txt （アウトファイルはoutfile.txtという名前でフォルダ内に作成）
-q 1000 （集団サイズqの事前分布は0<q<1000、一様分布）
-m 0.1 (移入率の事前分布は0<m<0.1、一様分布)
-t 30 (分岐年代の事前分布は0<t<30、一様分布)
-b 100000 （burninの試行数：統計に加えないデータの数：解析を始めてすぐの値は初期値の影響を強く受けるため）
-l 50000 （試行の長さ）
-d 100 （サンプリングの頻度） *-l 50000, -d 100なら、解析の長さは50,000*100=5,000,000サイクル
- s 123 （解析スタート時に「タネ」として使用する乱数）
-p2 （表示オプション: TMRCAのプロットを表示する）
-p7 （表示オプション: migration eventのヒストグラムを表示する）
-hfg （Heatingを使用する）
-hn40 （Heated chainを40本使用）
-ha, -hb (Heatingの度合いを決定するパラメータ）
というコマンド内容になっている。
5. 結果の解釈
outfileのうち.tiファイルにはパラメータの変遷が格納されており、結果が記述されているのはoutfile。
Outfileの中盤にあるHISTOGRAMSの項目で、出力された結果が確認できる。
HISTOGRAM GROUP 1: MARGINAL DISTRIBUTION VALUES AND HISTOGRAMS OF PARAMETERS IN MCMC
の項目において、t (mutation scaled population splitting time)、u (relative mutation rate)およびK (transition/transversion ratio)の分布が得られる。
Tを実際の年代、U1をローカス1の1塩基あたりの塩基置換速度、Lを塩基配列の長さとすると、以下の式でtを実際の歴年代に変換できる。
T = t / (U1 * L)
また、複数のローカスを使っている場合、各ローカスの進化速度の調和平均を算出する必要がある。
IMa2では、各ローカスの進化速度の比を算出してくれる(relative mutation rate: u1, u2, u3...un)
仮に3つのローカスS1, S2, S3 のうちS1の進化速度がU1であると分かっていたと仮定すると、
U2 = U1 * u2 / u1
U3 = U1 * u3 / u1
T= t / ((U1 * L1) * (U2 * L2) * (U3 * L3))^(1/3)
となる。
HISTOGRAM GROUP 2: MARGINAL DISTRIBUTION VALUES AND HISTOGRAMS OF POPULATION SIZE AND MIGRATION PARAMETERS
の項目においては、q (mutation scaled effective population size)およびm (mutation scaled migration rate)のヒストグラムが表示される。
Neを有効集団サイズ、Mを真の移入率とすると、
4Ne=q/(U1/u1/) * L
M=ml(U1/u1)
となり、それぞれの値が算出可能。
なお、塩基置換速度が不明な場合も、
NeM=qm/4
となり、一世代あたりの移入個体数であるNeMが算出可能。
HISTOGRAM GROUP 3: MARGINAL DISTRIBUTIONS OF TMRCA VALUES
の項目においては、TMRCA（もっとも最近の共通祖先の年代）が表示される。TMRCAはローカスごとにg0_tmrca、g1_tmrcaといった形で表示される。ローカス1のTMRCA（年）をTg1とすると、
Tg1=g1_tmrca/(U1/u1)/l
となる。
IMaモデルによって得られた上記のパラメータは、次のASCII Curves - Approximate Posterior Densitiesの項目で視覚的に確認できる。
Outfileの最後の表示されるASCII Plots of Parameter Trendsでは、出力の時系列変化を確認することができる。もし、解析が進むごとに値が一定の方向性をもって変化しているならば、十分な長さの解析が行われていない。たとえば、t0 trendが右肩上がりに変化しているとすれば、より長い解析を行い、より多くのburninを設定することで、t0の出力がより高い値に変化する可能性がある。この項目で、すべての値が何の時系列変化も示していないことを確認する。
（別記）infileの書式
ヘッダー
1-5行目がヘッダー部分です。
データはspaceで区切られます。
-------------------------------------------------------------------------------------------------------------------------
（解析のタイトル）
（集団数）
（集団1の名前）（space）（集団2の名前）・・・
（集団間の系統関係; Newickで記述）
（ローカス数）
-------------------------------------------------------------------------------------------------------------------------
ローカス情報
6行目以降に各ローカスの情報を記述します
DNA/RNAシーケンス
進化モデルがHKY（H）およびInfinite site（I）から選べます。
遺伝様式はAutosomal（=1）、X-linked（=0.75）およびmitochondrial（=0.25）。
配列全体の塩基置換速度を指定することができます（指定しなくてもよい）。
-------------------------------------------------------------------------------------------------------------------------
（ローカスの名前）（space）（集団１の個体数）（space）（集団2の個体数）・・・（space）前の行の続き；（塩基数）（space）（進化モデル）（space）（遺伝様式）（塩基置換速度）
（個体１の名前：10文字）（ローカス１の配列データ）
（個体2の名前：10文字）（ローカス１の配列データ）
・・・
-------------------------------------------------------------------------------------------------------------------------
Microsatellite
Stepwise mutation model（S）のほか、joint SSM（J）やHapSTR（IS）といったモデルも選択できます（筆者はまだやったことがない）。
-------------------------------------------------------------------------------------------------------------------------
（ローカスの名前）（space）（集団１の個体数）（space）（集団2の個体数）・・・（space）前の行の続き；（ローカス数）（space）（進化モデル）（ローカス数）（space）（遺伝様式）
（個体１の名前：10文字）（ローカス１の反復数）（space）（ローカス２の反復数）・・・
（個体2の名前：10文字）（ローカス１の反復数）（space）（ローカス２の反復数）・・・
・・・
-------------------------------------------------------------------------------------------------------------------------
実際の例
-------------------------------------------------------------------------------------------------------------------------
Analysis_1
2
Pop1 Pop2
(0,1):2
3
Locus1 2 2 16 H 0.25
Pop1-1 TAACGCACTTTTAACT
Pop1-2 TAACGCTCTTTTAACT
Pop2-1 CAACGCACTTTTATCA
Pop2-2 CAACGCTCTTTTATCT
Locus2 2 2 10 I 0.25
Pop1-1 CTCTTTTAAC
Pop1-2 CTCTTTTAAC
Pop2-1 CTCTATTAAT
Pop2-2 CTCTTTAAAT
Locus3 4 4 3 S3 1
Pop1-1 118 72 14
Pop1-1 114 72 15
Pop1-2 112 74 15
Pop1-2 112 74 15
Pop2-1 114 76 15
Pop2-1 114 78 16
Pop2-2 114 76 16
Pop2-2 114 78 16
-------------------------------------------------------------------------------------------------------------------------