全自動でシングルコピー遺伝子の連結系統樹推定(OrthoFinder, MAFFT, trimAL, IQ-TREE)
オーソログ推定プログラムのOrthoFinderですが、便利なことにオーソログ推定の結果予測されたシングルコピー遺伝子の情報を個別に抽出してディレクトリにまとめてくれます。
シングルコピー遺伝子の系統樹は種の系統関係を反映している可能性が高いと考えられるため、シングルコピー遺伝子の配列を連結して系統樹を推定すればそれなりに信頼のできる種の系統樹を簡単に得ることができます。
そんなわけで、系統関係を知りたい生物種の全遺伝子配列を入力としてOrthoFinderを実行し、その結果を基にシングルコピー遺伝子の連結系統樹推定まで全自動で実行するプログラムを作成したのでそのまとめです。
既に同様のパイプラインをQiitaで公開されている方がおりましたが、OrthoFinderの実行まで全自動化したかった、アライメントと非保存領域のトリミングをマルチスレッド実行したかった、これら2点の理由と後学のために自分でも作ることにしました。
コードは一番下に載せました。コピペで多分動くと思います。
プログラム実行書式
コマンドの書式は以下の通りです。解析対象生物ごとに全遺伝子情報が記載されたFASTAまたはGenBankファイルを用意しディレクトリにまとめておけば、シングルコピー遺伝子の連結系統樹推定まで実行可能です。なお実行には、Python 3.4以降が必須です。
$ python reconstruct_species_tree.py -i [入力ディレクトリ] -o [出力ディレクトリ]
実行後出力ディレクトリ以下に下記のように結果がまとめられます
$ tree [出力ディレクトリ] [出力ディレクトリ] ├── 01_fasta ├── 02_orthofinder ├── 03_single_copy ├── 04_alignment ├── 05_triming ├── 06_mltree └── species_tree.nwk
species_tree.nwkが推定したシングルコピー遺伝子の連結系統樹です。
依存関係
外部ライブラリとしてpandas、biopythonに依存しています。pipでインストール可能です。
$ pip install biopython pandas
また、OrthoFinder、MAFFT、trimAl、IQ-TREEが実行可能になっている必要があります。IQ-TREEはバージョン1.6以降をご使用ください。
プログラム概要
プログラムは6つのステップで実行されます。
- OrthoFinderの入力データ準備
- OrthoFinderによるオーソログ推定
- シングルコピー遺伝子配列の抽出及び連結系統樹推定のためのFASTAファイルのヘッダー書き換え
- MAFFTによるマルチプルアライメント
- trimALによるアライメント困難領域の削除
- IQ-TREEによる最尤連結系統樹推定
1. OrthoFinderの入力データ準備
OrthoFinderの実行には解析対象生物の全遺伝子配列が記載されたFASTAファイルが必要になりますので準備します。
FASTA形式でファイルを持っている場合は何もする必要がありませんが、GenBank形式で全遺伝子情報を入手している場合もあるためプログラム中でCDSの翻訳配列を抜き出すことで対応することにしました。
ファイル形式の判別及びCDS配列の抽出はBiopythonを使いました。
2. OrthoFinderによるオーソログ推定
OrthoFinderでシングルコピー遺伝子を推定します。配列だけ入手できれば良いので、-ogオプションを使って解析を途中で止めています。
内部で実行されるコマンドは次のとおりです。
$ orthofinder -f [in] -a [threads] -og
3. シングルコピー遺伝子配列の抽出及び連結系統樹推定のためのFASTAファイルのヘッダー書き換え
この後IQ-TREE中で各シングルコピー遺伝子についてモデル推定を行い連結系統樹推定を実行しますが、この際FASTAファイルのヘッダーが同一の配列が同じOTU由来として認識されるため、個々の生物種の遺伝子についてヘッダーを同一にしておく必要があります。OrthoFinderでは入力ファイルのファイル名が生物種名として認識されるのでそれで統一することにしました。
OrthoFinderの結果として各オーソロググループのFASTAファイルの他に、各オーソロググループについて配列名と生物種名をまとめたTSVファイル(Orthogroups.tsv)とシングルコピー遺伝子のグループ番号をまとめたリスト(Orthogroups_SingleCopyOrthologues.txt)が出力されるため、これを利用してシングルコピー遺伝子配列ファイルの抽出とヘッダーの書き換えを行い、改めて別ディレクトリに配列ファイルを出力しています。
tsvファイルの処理にpandasを利用しています。
4. MAFFTによるマルチプルアライメント
MAFFTでマルチプルアライメントを行います。MAFFTは並列実行にも対応していますが、再現性が無い点が少し気になったので、シングルコアで実行しています。ただ、それだけだと遅くなってしまうためmultiprocessingを利用して複数コアで同時に各配列のアライメントを行うことにしました。
内部で実行されるコマンドは次の通りです。
$ mafft --maxiterate 1000 --localpair --anysymbol [in] > [out]
5. trimAlによる非保存領域の削除
アライメントが不完全な領域を系統樹推定に用いることを防ぐためこれを取り除きます。トリミングにはtirmAlを用いました。オプションが多いプログラムですが、とりあえず最尤系統樹推定に適しているらしい-automated1オプションで実行しています。
またMAFFTと同様にtrimAlの実行はシングルスレッドですがmultiprocessingモジュールを利用して同時に複数ファイルについてトリミングを行っています。
内部で実行されるコマンドは次の通りです。
$ trimal -automated1 -in [in] -out [out]
6. IQ-TREEによる最尤連結系統樹推定
トリミング済みのシングルコピー遺伝子配列を基に、連結系統樹をIQ-TREEで推定します。IQ-TREEでは設定ファイルでパーティションを指定することにより、配列ごとに最適置換モデル推定を行うことができるためこれを利用しています。
必要なパーティションファイルをプログラム中で自動生成して読み込ませています。
内部で実行されるコマンドの書式は次の通りです。
$ iqtree -sp [partition] -bb 1000 -pre [prefix] -nt [threads]
コード
#!/usr/bin/python3 """\ reconstruct_species_tree.py -i <dir> -o <dir> [-t <int>] [-f] [-h] [OPTIONS]: --indir |-i <dir> Path to the directory containing all gene sequence files (FASTA or GenBank format) for each organism --outdir |-o <dir> Path to the result directory --threads|-t <int> Number of threads to use [default: 'max'] --fast |-f Use --auto option in mafft to speed up alignment --help |-h show this help message and exit """ import argparse import multiprocessing import os from pathlib import Path import shutil import subprocess from Bio import SeqIO import pandas as pd class ArgumentParser(argparse.ArgumentParser): def format_help(self): return __doc__ parser = ArgumentParser() parser.add_argument( "-i", "--indir", required=True, type=str, metavar="<dir>", help="Path to the directory containing all gene sequence files" "(FASTA or GenBank format) for each organism" ) parser.add_argument( "-o", "--outdir", required=True, type=str, metavar="<dir>", help="Path to the result directory" ) parser.add_argument( "-t", "--threads", default=len(os.sched_getaffinity(0)), type=int, metavar="", help="Number of threads to use [default:all available threads]" ) parser.add_argument( "-f", "--fast", default=False, action="store_true", help="Use --auto option in mafft to speed up alignment" ) options = parser.parse_args() class SpeciesTreeReconstructor: def __init__(self, indir, outdir, num_of_threads, fast=False): self.threads = num_of_threads self.fast = fast self.seq_dir = Path(indir) self.out = Path(outdir) self.fasdir = self.out / "fasta" self.orthofinder_out = self.out / "orthofinder" self.single_copy_dir = self.out / "single_copy_orthologues" self.mafft_out = self.out / "alignment" self.trimal_out = self.out / "triming" self.iqtree_out = self.out / "mltree" def orthofinder_input(self): seqs = self.seq_dir.glob("*.*") for seq in seqs: seq = str(seq) record = SeqIO.parse(seq, "fasta") is_fasta = any(record) if is_fasta: shutil.copy(seq, str(self.fasdir)) else: self.extarct_proteins(seq) def reconstruct_species_tree(self): self.out.mkdir() self.fasdir.mkdir() self.orthofinder_out.mkdir() self.single_copy_dir.mkdir() self.mafft_out.mkdir() self.trimal_out.mkdir() self.iqtree_out.mkdir() print("***** Ortholog inference by OrthoFinder *****") self.orthofinder_input() self.orthofinder(self.fasdir) self.move_orthofinder_result() print("***** Extract single copy gene sequences") self.replace_headers_by_species_name() print("***** Multiple alignment by mafft *****") self.mafft() print("***** Trimming non-conserved sequence region by trimAL *****") self.trimal() print("***** Reconstruct ML tree by IQ-TREE *****") self.iqtree() def extarct_proteins(self, gbk): gbk = Path(gbk) fasdir = Path(self.fasdir) prefix = gbk.stem outfas = fasdir / "{0}.faa".format(prefix) with open(str(outfas), "w", encoding="UTF-8") as out: for record in SeqIO.parse(str(gbk), "genbank"): for feature in record.features: if feature.type != "CDS": continue try: protein_id = feature.qualifiers["protein_id"][0] anno = protein_id translation = feature.qualifiers["translation"][0] print(">{0}\n{1}".format(anno, translation), file=out) except KeyError: continue def orthofinder(self, fasdir): fasdir = str(fasdir) cmd = ["orthofinder", "-f", fasdir, "-a", str(self.threads), "-og"] subprocess.run(cmd) def move_orthofinder_result(self): tmp = self.fasdir.glob("**/Results_*") results = next(tmp) dest = str(self.orthofinder_out) for result_dir in results.glob("*"): shutil.move(str(result_dir), dest) empty_result = self.fasdir / "OrthoFinder" shutil.rmtree(str(empty_result)) def replace_headers_by_species_name(self): result = self.orthofinder_out table = next(result.glob("**/Orthogroups.tsv")) id_file = next(result.glob("**/Orthogroups_SingleCopyOrthologues.txt")) tmp = pd.read_table(id_file, header=None) single_copy_list = list(tmp[0]) df = pd.read_table(str(table)) unnamed = df.columns[0] new_name = "group" df = df.rename(columns={unnamed: new_name}) singles = df[df[new_name].isin(single_copy_list)] seq_id2spe = {} for species, column in singles.iteritems(): tmp = {seq_id: species for seq_id in column} seq_id2spe.update(tmp) tmp = result.glob("**/Single_Copy_Orthologue_Sequences") single_copy = next(tmp) single_copy_seqs = single_copy.glob("*.*") for single_copy_seq in single_copy_seqs: file_name = single_copy_seq.name dest = self.single_copy_dir / file_name with open(str(dest), "w", encoding="UTF-8") as out: for record in SeqIO.parse(str(single_copy_seq), "fasta"): seq_id = record.description species = seq_id2spe[seq_id] species = species.replace(" ", "_") seq = record.seq new_record = ">{0}\n{1}".format(species, seq) print(new_record, file=out) @staticmethod def exec_cmd(cmd): return subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL) def mafft(self): indir = self.single_copy_dir seqs = indir.glob("*.*") if self.fast: args_tuples = [("mafft", "--anysymbol", "--auto", str(seq)) for seq in seqs] pos = 3 else: args_tuples = [("mafft", "--maxiterate", "1000", "--localpair", "--anysymbol", str(seq)) for seq in seqs] pos = 5 with multiprocessing.Pool(self.threads) as pool: for proc in pool.imap_unordered(self.exec_cmd, args_tuples): args = proc.args result = proc.stdout.decode("utf8") seq = Path(args[pos]) group_id = seq.stem outfile = self.mafft_out / "{0}.aln".format(group_id) with open(str(outfile), "w", encoding="UTF-8") as out: out.write(result) def trimal(self): indir = self.mafft_out seqs = indir.glob("*.*") args_tuples = [] for seq in seqs: group_id = seq.stem infile = str(seq) outfile = str(self.trimal_out / "{0}.trim".format(group_id)) args = ("trimal", "-automated1", "-in", infile, "-out", outfile) args_tuples.append(args) with multiprocessing.Pool(self.threads) as pool: for proc in pool.imap_unordered(self.exec_cmd, args_tuples): args = proc.args result = args[5] size = os.path.getsize(result) if not size > 0: seq = args[3] dest = result shutil.copy(seq, dest) def iqtree(self): indir = self.trimal_out seqs = indir.glob("*.*") config = "partitions.nex" config_path = self.iqtree_out / config with open(str(config_path), "w", encoding="UTF-8") as out: begin = "#nexus\nbegin sets;" print(begin, file=out) for index, seq in enumerate(seqs, 1): index = str(index) seq = seq.resolve() seq = str(seq) partition = "\tcharset part{0} = {1}: ;".format(index, seq) print(partition, file=out) end = "end;" print(end, file=out) prefix = "tree" cmd = ["iqtree", "-spp", config, "-bb", "1000", "-pre", prefix, "-nt", str(self.threads)] subprocess.run(cmd, cwd=str(self.iqtree_out)) tree_file = str(self.iqtree_out / "{0}.treefile".format(prefix)) new = str(self.out / "species_tree.nwk") shutil.copy(tree_file, new) if __name__ == "__main__": tree_indir = options.indir tree_outdir = options.outdir threads = options.threads mafft_fast = options.fast species_tree_reconstructor = SpeciesTreeReconstructor(tree_indir, tree_outdir, threads, mafft_fast) species_tree_reconstructor.reconstruct_species_tree()
