iq-treeとは

モデル選択も自動でしてくれる最尤系統樹作成ソフト（高速）です。
ゲノム時代の系統解析ソフトを目指したという触れ込み通り、
複数遺伝子を使った解析が非常にしやすくなっています。
各遺伝子を区別して、それぞれに最適なモデルを選択したあと、個別にパラメータ推定をして、その全データを使った系統樹を書くなどできます。
PAUP時代のヒトには懐かしい感じ。

参考サイト

http://biopython.org/wiki/SeqIO
http://www.stevekellylab.com/software/orthofinder

1,シングルコピー遺伝子を選出する

orthofinder の出力ファイル（Orthogroups_.csvなど）からsingle copy遺伝子だけを選抜し、以下の様なスペース区切りのファイルを作る。

1列目がオーソロガスグループ名(OG)で、以降が各OGに属する配列の名前になっている。

1	2	3	4
OG0000000	fasta_ID_sp_A_seq1	fasta_ID_sp_A_seq2	fasta_ID_sp_A_seq3
OG0000001	fasta_ID_sp_B_seq2	fasta_ID_sp_B_seq2	fasta_ID_sp_B_seq2
OG0000002	fasta_ID_sp_C_seq3	fasta_ID_sp_C_seq3	fasta_ID_sp_C_seq3

OG_list.txt

OG0000000 fasta_ID_sp_A_seq1 fasta_ID_sp_B_seq2 fasta_ID_sp_C_seq3
OG0000001 fasta_ID_sp_A_seq2 fasta_ID_sp_B_seq2 fasta_ID_sp_C_seq3 
OG0000002 fasta_ID_sp_A_seq3 fasta_ID_sp_B_seq2 fasta_ID_sp_C_seq3

全種の全タンパク質の配列情報が入っているmultiple fastaも準備する。
(各生物種毎に全CDS情報をfastaでとってきてcatで結合などする）

all_seq.fa

>fasta_ID_sp_A_seq1
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
>fasta_ID_sp_B_seq2
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
>fasta_ID_sp_C_seq3
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
........

2,OG毎に各配列情報を別ファイルにいれ分ける

以下のスクリプトを使って、データを分ける

各OG名がファイル名になったfastaが作成される。
各ファイルに、それぞれのOGに属する配列が入っている。

OG_separater.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys
from Bio import SeqIO

fasta_in = sys.argv[1]                            #１番目の引数には上記のall_seq.faなどfastaファイルを指定する
query_in = sys.argv[2]                            #２番目の引数には上記のOG_list.txtなどオーソログファイルを指定する

for q in open(query_in, "r"):                     #オーソログファイルを開いて1行づつ読み込む
    query = q.split()                             #スペース毎に切りとってリスト形式でqueryに保存する
    f = open(query[0], 'w')                       #最初の列(OG名）と同じ名前のファイルを作成する
    for record in SeqIO.parse(fasta_in, 'fasta'): #fastaファイルを開くSeqIOを使ってパースする(1項目づつ読み込む）
        id_part = record.id                       #fastaのID部分を読み込む
        desc_part = record.description            #fastaのdescription部分を読み込む
        seq = record.seq                          #fastanの配列部分を読み込む
        for i in range(len(query)):              #オーソログファイル中の各OGに含まれる配列数を数えて、その分繰り返す(python2の人はrange を x rangeにする)
            if desc_part == query[i] :            #オーソログファイルの配列descriptionとfastaの配列descriptionが一致したら、、、
                fasta_seq = '>' + desc_part + '\n' + seq + '\n'     #fasta形式に整えて
                print(fasta_seq)                  #標準出力にfastaを出力（進行状況把握用）
                f.write(str(fasta_seq))           #各OGファイルにfastaを出力
    f.close()

3,アライメントとトリミング

mafft と trimalを使って各OGファイルをアライメント＋トリミングする
trimalの設定はML法に適しているというautomated1 をつかった。

align.sh

#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
awk '{print($1)}' $1 | while read x; do #引数に前述のOG_list.txtなどのOGリストを指定する。
    mafft --auto $x > $x.maffted.fa
    trimal -in $x.maffted.fa -out $x.maffted.trimed.fa -htmlout $x.maffted.trimed.fa.html -automated1
done

htmlファイルにトリミングの様子が出るので、変なトリミングなどが無いか確認しておく。

4.OTU名を揃える

各fastaファイル中の配列名は、ファイル間で統一した種名に変える。
色々方法はあるが、例えば、
OG_list.txtと同じ順番で種名をタブ区切りで並べたファイルspecies.list.txtを作成し
以下のスクリプトで変更する。

change_OTU.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys
from Bio import SeqIO

query_in = sys.argv[1]                            #1番目の引数には上記のOG_list.txtなどオーソログファイルを指定する
species_in = sys.argv[2]                            #OG_list.txtと同じ順番で種名を記述したファイルを指定する

for sp in open(species_in, "r"):                     #種名ファイルを開く
    sp_list = sp.split()                            #スペース毎に切りとってリスト形式でsp_listに保存する

for q in open(query_in, "r"):                     #オーソログファイルを開いて1行づつ読み込む
    query = q.split()                             #スペース毎に切りとってリスト形式でqueryに保存する
    f = open(query[0]+".maffted.trimed.edit.fa", 'w')
    for record in SeqIO.parse(query[0]+".maffted.trimed.fa", 'fasta'): #fastaファイルを開くSeqIOを使ってパースする(1項目づつ読み込む）
        desc_part = record.description            #fastaのdescription部分を読み込む
        seq = record.seq                          #fastanの配列部分を読み込む
        for i in range(len(query)):              #オーソログファイル中の各OGに含まれる配列数を数えて、その分繰り返す
            if desc_part == query[i] :            #オーソログファイルの配列descriptionとfastaの配列descriptionが一致したら、、、
                fasta_seq = '>' + sp_list[i-1] + '\n' + seq + '\n'     #配列名を該当する種名に置き換えて、fasta形式に整えて
                print(fasta_seq)                  #標準出力にfastaを出力（進行状況把握用）
        f.write(str(fasta_seq))           #各OGファイルにfastaを出力
    f.close()

5,iq-treeのrun設定ファイルを作る

iq-treeはRunの設定をNexusファイル（懐かし！）で作成する。

以下のような感じ
part_の部分は任意の名前
=以下で各OGのファイルを指定することで、それぞれのOG（遺伝子）ごとに別個にモデル選択などをしてくれる。

run.nex

#nexus
begin sets;
    charset part1 = PsbL.maffted.trimed.fa: ;
    charset part2 = accD.maffted.trimed.fa: ;
    charset part3 = atpB.maffted.trimed.fa: ;
    charset part4 = atpE.maffted.trimed.fa: ;
    charset part5 = atpI.maffted.trimed.fa: ;
    charset part6 = chlB.maffted.trimed.fa: ;
    charset part7 = chlL.maffted.trimed.fa: ;
    charset part8 = chlN.maffted.trimed.fa: ;
    charset part9 = clpP.maffted.trimed.fa: ;
    charset part10 = cysA.maffted.trimed.fa: ;
    charset part11 = infA.maffted.trimed.fa: ;
    charset part12 = petA.maffted.trimed.fa: ;
end;

6,run実行

以下のコマンドで1000回のBSをしてくれる

iqtree -sp run.nex -nt AUTO -bb 1000

orthofinder + mafft + trimal + iq-treeでゲノム規模のデータから系統樹を作る