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

ゲノム規模のタンパク質情報をもとに種の系統作成をするパイプラインです。

本家
http://www.iqtree.org

mafft
https://mafft.cbrc.jp/alignment/software/

trimal
http://trimal.cgenomics.org

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