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

2023/09/27
Salichos L, Rokas A 2013.
Inferring ancient divergences requires genes with strong phylogenetic signals. Nature 497: 327-+.
などによって複数配列をconcatenateしての系統樹作成は推奨されていません
ASTRAL等を使って、遺伝子ごとの系統樹を統合する形での系統樹作成をおすすめします。

2019/03/06 scriptを変更

2019/09/30
改良したスクリプトをgithubに上げました。
PhyHat; A pipeline for ML-tree analysis with multiple gene set
https://github.com/maedat/PhyHat

#prequal

trimalの代わりになりそうな、配列処理ソフトprequalが発表されました。

PREQUAL: a program for detecting errors in sets of unaligned homologous sequences. Whelan, Irisarri & Burki. In Press. Bioinformatics. https://doi.org/10.1093/bioinformatics/bty448

詳しくは上坂さんのページと本家サイトを参照

繰り返し領域などをマスクすることができます。
ゲノムスケールデータには、Nベースやモデリングの不備に由来する欠損領域や明らかにおかしい繰り返しなどがあります。
異常なポジションを削っていくtrimalでは、他の種ではちゃんと読めている配列情報を取り逃しているかもしれません。prequalはそこを、変な配列だけをXでマスクすることで、他の種の情報を利用可能にしています。

#パイプライン作った。

ケースバイケースで使い分けが必要と思いますが、評判も良いので、それを使ったパイプラインを作成しました。

また以前に作成した、３つのパートに分かれていた種系統樹作成スクリプトを、１つに統合してみました。
前回のはこちら

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

tree.py

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

#python >3.5

import argparse
from Bio import SeqIO
import subprocess

parser = argparse.ArgumentParser()

parser.add_argument('-f','--fasta',  help="File path to a fasta including all species sequence", required=True)
parser.add_argument('-d','--group', help="Ortholog grouping file", required=True)
parser.add_argument('-n','--name',help="Species name in order (Space-separated)", type=str, required=True)

args = parser.parse_args()

fasta_in = args.fasta
query_in = args.group
sp = args.name
sp_list = sp.split()

print(fasta_in)
print(query_in)
print(sp_list)

fnex = open("iqtree.run.nex", 'w')
fnex.write("#nexus \n begin sets;\n")

for q in open(query_in, "r"):
    query = q.split()
    f = open(query[0], 'w')
    print("<" + query[0]+ ">")
    for record in SeqIO.parse(fasta_in, 'fasta'):
        id_part = record.id
        desc_part = record.description
        seq = record.seq
        for i in range(len(query)):
            if desc_part == query[i] :
                fasta_seq = '>' + desc_part + '\n' + seq + '\n'
                print(desc_part)
                f.write(str(fasta_seq))
    f.close()
    print("<OTU rename>")
    fsp = open(query[0] + ".sp.fa", 'w')
    for record in SeqIO.parse(fasta_in, 'fasta'):
        id_part = record.id
        desc_part = record.description
        seq = record.seq
        for i in range(len(query)):
            if desc_part == query[i] :
                print(sp_list[i-1])
                fasta_seq = '>' + sp_list[i-1] + '\n' + seq.rstrip("*") + '\n'
                fsp.write(str(fasta_seq))
    fsp.close()    
    subprocess.run("prequal "+query[0] + ".sp.fa", shell=True)
    subprocess.run("mafft --auto "+query[0]+".sp.fa.filtered"+" > "+query[0]+".sp.fa.filtered.maffted.fa", shell=True)    
    fnex.write("charset "+query[0]+" = "+query[0]+".sp.fa.filtered.maffted.fa: ; \n ")
fnex.write("end;")
fnex.close()

subprocess.run("iqtree -sp iqtree.run.nex -nt AUTO -bb 1000", shell=True)

使用例

python3 tree.py -f all.fasta -d og.single.txt -n "Mmu Hsa Dme"