Edited at

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

More than 1 year has passed since last update.

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

本家

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] #1番目の引数には上記のall_seq.faなどfastaファイルを指定する
query_in = sys.argv[2] #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