#2020/08/16追記

以下に更新版の情報を置いてあります。
「orthofinder + prequal+ mafft + iq-treeでゲノム規模のデータから系統樹を作る」
https://qiita.com/MaedaTaro_Umiushi/items/d4004d3fb219a4f991e8

「orthofinder + mafft + trimal + iq-treeでゲノム規模のデータから系統樹を作る」
https://qiita.com/MaedaTaro_Umiushi/items/7c8cee435347eeee1cf5

#はじめに

ゲノム情報をもとに数百の遺伝子を使って系統解析をしようとしたのですが、ファイルの結合などのいいツールが見つからなかったので、Biopythonでそこいらの処理が楽になるようにしました。またConselをつかったAU検定でちょっとつまずいたので、そこについても記録しておきます。

初心者なので、間違いなどのコメント歓迎です。

#参考サイト
http://biopython.org/wiki/SeqIO
http://www.stevekellylab.com/software/orthofinder

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

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

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

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 xrange(len(query)):              #オーソログファイル中の各OGに含まれる配列数を数えて、その分繰り返す
            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ファイルをアライメント＋トリミングする

align.sh

#!/bin/sh
#$ -S /bin/bash
#$ -cwd
#$ -v PATH
mafft --maxiterate 1000 --localpair --thread 10 $1 > $1.maffted
trimal -automated1 -in $1.maffted -out $1.maffted.trimal.fa

#4,複数のOGファイルから、順番そのまま、fastaを横に結合していく
順番そのまま、fastaを横に結合していく（bashのpasteコマンドっぽく結合する）スクリプトを作成した

例

GeneA.fa

>sp1_A
AAAAAAAAAA
>sp2_A
BBBBBBBBB
>sp3_A
CCCCCCCCCC
>sp4_A
DDDDDDDDDD

GeneB.fa

>sp1_B
EEEEEEEEEE
>sp2_B
FFFFFFFFFF
>sp3_B
GGGGGGGGGG
>sp4_B
HHHHHHHHHH

out.fa

>sp1_A;sp1_B
AAAAAAAAAAEEEEEEEEEE
>sp2_A;sp2_B
BBBBBBBBBFFFFFFFFFF
>sp3_A;sp3_B
CCCCCCCCCCGGGGGGGGGG
>sp4_A;sp4_B
DDDDDDDDDDHHHHHHHHHH

無限に引数をとれるようになっているので以下の様にワイルドカードをつかってもOK

joint.sh

python fasta_merge_orderkeep.py *.maffted.trimal.fa

以下スクリプト

fasta_merge_orderkeep.py

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

#TARO MAEDA 2017 14 Jul. 

import sys
from Bio import SeqIO

num_file=len(sys.argv)-1  #入力ファイル数を数える
print("Detected file number")
print(num_file)

id_list = {}              #IDの辞書を宣言
seq_list = {}             #seqの辞書を宣言
counter_file = int("0")   #ファイルカウンターを初期化


print("File opening...") 

for i in xrange(num_file):                    #ファイル数分繰り返す（fastaファイルの中身を辞書に入力)
    counter_file += 1                         #ファイルカウンターをまわす
    counter_order = int("0")                  #ファイル中のfastaの順番をカウントするorderカウンターを初期化する
    for file1 in SeqIO.parse(sys.argv[counter_file], 'fasta'):  #引数N番目のファイルを読み込む
        counter_order +=1                     #orderカウンターをまわす
        print("now getting " +  file1.id)     #状況確認用（取り込み中のfasta idを出力する）
        id_list["ID_file" + str(counter_file) + "order" + str(counter_order)] = file1.id    #fastaIDを辞書に入力していく、orderカウンターとファイルカウンタでキーの値を指定していく。
        seq_list["Seq_file" + str(counter_file) + "order" + str(counter_order)] = file1.seq #fastaの配列情報を辞書に入力していく、orderカウンターとファイルカウンタでキーの値を指定していく。


max_count_file = counter_file    #最終的に読まれたすべてのファイル数を獲得する
max_count_order = counter_order  #最終的に読まれたすべての種数（1fastaファイル中のfastaデータ数）を獲得する
now_order = int("0")                    #今出力している種が何番目の種か（fastaファイル中の何番目のfastaデータか）を数えるカウンターを初期化する


print("File writing...") 


f = open("marged_out.fa", 'w')            #結果出力用ファイルを開く
for now_order in xrange(max_count_order): #種数分処理を繰り返す
    out_id  = ""                          #出力用変数を初期化する。
    out_seq  = ""                         #出力用変数を初期化する。
    for now_file in xrange(max_count_file): #全ファイル数分繰り返す
        out_id += id_list["ID_file" + str(now_file+1) + "order" + str(now_order+1)] +";" 
        out_seq += str(seq_list["Seq_file" + str(now_file+1) + "order" + str(now_order+1)]) 
    outfasta =">" + out_id + "\n" + out_seq + "\n"
    f.write(outfasta)
f.close()

print("Output file was made as marged_out.fa")