SNP連結配列の最尤系統樹を推定するときのTips
2021-12-30
大貫が引っかかったところを解決するためのTipsの紹介です。
VCF→phylip（or fasta）の変換
SNPデータを使う場合，フィルタリングなどを行った後のデータはVCFになっていることが多いと思います。
RAxML-NGやIQ-TREEにインプットする際の一つの方法は，VCFからSNPデータのみを取り出して連結配列とし，phylip形式に変換することです（fastaでもいいですが，後述のraxml_ascbiasのインプットがphylipなのでこちらをおすすめします）。
なお，MIG-Seqのように，invariant sitesも含めた全配列をインプットしても現実的な計算時間で系統樹を推定できる場合には，別の方法をとった方がいいと思います。機会があれば追記します。
ASC関連のトラブルシューティング
SNP連結配列の最尤系統樹推定では，塩基置換モデルにascertainment bias correction（ASC）を含めて解析する必要があります。SNP連結配列にはinvariant sitesがないためです（詳しいやり方はYF_bioさんのQiita，詳しい理由はLeaché et al. (2015)を読んでください）。 ところが，ASCを含めた塩基置換モデルで解析しようとすると，RAxML-NGやIQ-TREEが以下のようなエラーを吐くことがあります。
code: RAxML-NGのエラー
ERROR: You enabled ascertainment bias correction for partition noname, but it contains （数字） invariant sites.
This is not allowed! Please either remove invariant sites from MSA or disable ascertainment bias correction.
code: IQ-TREEのエラー
ERROR: Invalid use of +ASC because of （数字） invariant sites in the alignment
もしくは
Skipped since +ASC is not applicable
ぼくは「いわゆるinvariant sitesは除いたはずなのにおかしいな？」となりました。
これはinvariant sitesの定義の齟齬（もしくは勘違い？）が原因でした。ASCが仮定するinvariant sitesにはいくつかタイプがあり，それら全てを除去しないと動かないそうです。IQ-TREEのFAQにあった解説がわかりやすかったので，和訳しておきます。 ASC（ascertainment bias correction）モデルを使用した場合、エラーメッセージが表示されることがあります。
ERROR: Invaid use of +ASC because of ... invariant sites in the alignment.
またはモデルテストを行う際に
Skipped since +ASC is not applicable
これは、アライメントにモデルの数学的条件に違反するinvariant sites（columns）が含まれているためです。invariant sitesには次のようなものがあります。
・Constant sites：すべての配列で一文字の状態を含む。例えば、DNAアライメントにおいて、すべての配列が特定の部位でA（アデニン）を示す場合。
・Partially constant sites：単一の文字、ギャップ、または未知の文字を含む。例えば、ある部位で、ある配列はG（Guanine）を示し、ある配列は-（gap）を示し、他の配列はNを示す。
・Ambiguously constant sites。例えば、ある配列はC（シトシン）、ある配列はY（CまたはTの意味）、ある配列は-（ギャップ）を示す。
これらの部位は、+ASCモデルを適用する前に、すべてアライメントから削除する必要があります。
つまり，ASCを含めたモデルを適用するには，Constant sites，Partially constant sites，Ambiguously constant sitesを全て除去する必要があるということです。
これを自動でやってくれるのが，raxml_ascbiasです。これはpythonスクリプトで，動かすには，python3を入れた上，numpyとpandasとbiopythonを入れておく必要があります（このブログも導入がそこそこ詳しく書いてあります。アンチcondaの人なら，おそらくそんなに支障はないと思います）。raxml_ascbiasの具体的な使い方はGitHubに書いてあるので各自調べてください。 勝手にアミノ酸配列と判定されちゃうトラブルシューティング（IQ-TREE）
SNP連結配列には複合塩基が多数含まれるケースがしばしばあります。
IQ-TREEでは，与えられたアライメント中に複合塩基が多すぎてATGCの割合が小さいと，「アライメントはアミノ酸配列である」と判断し，コマンドに書かれた塩基置換モデルを受けつけなかったり，ModelFinderを使った場合には自動でアミノ酸置換モデルを当てはめてしまったり（悪質！）することがあります。