Pythonのfileinputでフィルタとかファイルの書き換えとか
Pythonでファイルの読み込みや書き込みをする際にはopen関数を使うのが基本ですが、下に挙げるような少し複雑な操作をしたい場合には頭を抱えることになります。
- 標準入力からも読み込みたい、ただし指定したときにはファイルからも読み込んでほしい。
- ファイルの読み書きを同時に行いたい。
特に2.が個人的に深刻で、良い方法が見つからず毎回無駄な中間ファイルを作っては、消すほどではないし残しとくかぁみたいな思考になること度々でした。今回上の問題をまとめて解決してくれる便利なモジュールfileinputを見つけたので備忘録を兼ねたまとめです。
標準入力から読み込む
以下公式のドキュメントからのコピペですが、fileinputからファイルを読み込む基本形は次の通りです。
import fileinput for line in fileinput.input(): process(line)
要するにfileinputのinput関数にコマンドライン引数として与えたファイルが順次渡されて開かれていきますよ、ということらしいです。開かれた後は1行ずつ処理することができます。複数のファイルを渡す場合でもループ処理のコードが1回で済むので便利ですね。
ファイル名を省略したり-にしたりすると標準入力から読み込む挙動のため、フィルタプログラムのようなモノを簡単に作成することができます。例えばファイルまたは標準入力を1行ずつ全て足して出力するというプログラムは次のように書けます。
上の例と少し違いますが、with文と一緒に使えるのでそちらで書いています。まぁ、この辺りはopenを使う場合と同じですね。
#!/usr/bin/python3 import fileinput total = 0 with fileinput.input() as f: for line in f: line = line.rstrip() total += int(line) print(total)
上のスクリプトを仮にsum.pyとすると次のように動作します。
$ echo -e "1\n2\n3\n4\n5" | ./sum.py 15
こういった処理はawkでやるのが王道な気がしますがPythonで書けるのは便利ですね。
ファイルの書き換え
fileinputのinput関数にはinplace引数と呼ばれるものがあり、これをTrueにすることでファイルの読み書きを同時に行うことができます。
公式ドキュメントの説明によると、inplace=Trueを指定した場合には入力ファイルがバックアップファイルに一時的に移動され(デフォルトではバックファイルの名前は入力ファイル名+.bak)、標準出力が入力ファイルに設定されます。また、バックアップファイルは標準出力に設定されたファイルが閉じられると自動で削除されます。
これだけ聞くと、新しくファイルを作ってそこに書き込んだ後、元の入力ファイルと同じ名前にリネームするような単純な処理と大差ない気もするのですが、コードの行数が減るのは良いですね。
例えばファイル中のアルファベットを全て小文字に書き換えるようなプログラムは次のようになります。
#!/usr/bin/python3 import fileinput with fileinput.input(inplace=True) as f: for line in f: line = line.rstrip() line = line.lower() print(line)
このプログラムを仮にto_lower.pyとすると挙動は次の通りです。
$ cat upper.txt AAA BBB $ ./lower.py upper.txt $ cat upper.txt aaa bbb
行を削除したい場合には単純にスキップすれば良いようです。>で始まる行を削除するコードを考えてみました。
#!/usr/bin/python3 import fileinput with fileinput.input(inplace=True) as f: for line in f: if line.startswith(">"): continue print(line, end="")
挙動は次の通りです。一応スクリプト名はdelete.pyです。
$ cat test.fas >A aaa >B bbb $ ./delete.py test.fas $ cat test.fas aaa bbb
惜しむらくは置換したり、削除する部分をプログラム中で指定してしまっている点ですね。argparseと組み合わせるとsedもどきが作れそうな気がするので、機会があれば挑戦してみたい所です。
何はともあれ最後までお読みいただきありがとうございました。
参考リンク
全自動でシングルコピー遺伝子の連結系統樹推定(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()
Pythonでディレクトリを一時的に移動して元のディレクトリに戻る
Pythonスクリプト内で外部プログラムを呼び出して解析を実行することが多々あるのですが、プログラムの中には出力先を指定するオプションがなくカレントディレクトリに結果を出力することを想定しているものがあります。
条件が一つだけの場合はそれで良いのですが、わざわざPythonを使って別のプログラムを呼び出そうとしているときは、複数の条件を一度に実行したいという場合が多いです。
そのためこのような場合は、条件ごとに解析ディレクトリを作りPythonスクリプト内でそこに移動し、プログラムを実行するという処理が必要になります。面倒ですね
そこで今回はPythonを使ってディレクトリを移動する方法と、最後に実行時のディレクトリに帰ってくる方法を備忘録を兼ねてまとめました。
カレントディレクトリを調べる
実行時のカレントディレクトリを取得し最後に帰ってこられるようにします。osモジュールのgetcwd関数を使うと簡単に調べられます。
>>> import os >>> # /home/userにいるとする ... os.getcwd() '/home/user'
あとはこれを適当な変数に記憶させるだけです。
ディレクトリを移動する
Pythonスクリプト内でディレクトリを移動する場合は、osモジュールのchdir関数を使うと簡単です。
... 上の続き >>> os.chdir("..") >>> os.getcwd() '/home'
コード
以上をまとめて今回作りたかった外部プログラム呼び出しのために一時的にディレクトリを移動し、実行後始めのディレクトリに戻ってくるスクリプトは以下の通りです。
#!/usr/bin/python3 import os import subprocess def main(indir): cwd = os.getcwd() print(cwd) os.chdir(indir) cmd = "pwd" subprocess.call(cmd) os.chdir(cwd) subprocess.call(cmd) if __name__ == "__main__": dest = "/bin" main(dest)
このスクリプトを上と同じく"/home/user"で実行すると出力結果は以下の通りになります。
/home/user /bin /home/user
期待通りディレクトリを移動することができました。後は任意のプログラムを実行するだけです。
最後までお読みいただきありがとうございました。
Guakeがexitでフリーズするので解決策と代替アプリを探した件
キーボードからさっと簡単に呼び出せるドロップダウン型のターミナルGuakeを愛用していたのですが、Ubuntu 18.04 LTSにしてから、exitで端末を閉じようとするとフリーズするようになってしまいました。
まぁ、Guakeの特性上基本的に端末は開きっぱなしでOKなので、気をつければよいだけなのですが、つい他のターミナルエミュレータと同じくCtrl+Dで終了しようとしてしまい、フリーズしてストレスを募らせるという精神衛生上良くない状況に陥っていました。
というわけで、Guakeのフリーズを解消する方法を探したのでそのまとめです。あと良い機会なので、ついでにGuakeっぽい使い方ができる他のターミナルエミュレータも調べてみました。
書いておいてなんですが、最新版では既にこのバグは解消されているようで、そのうちリポジトリにもアップデートが来るような気もします。
解決手順
同様のバグ報告が既に上がっており、アッサリ解決法が見つかりました。こちらです。要するに「Githubにコードの修正箇所が載っているからそれと同じようにしたら直ったぜ」とのことです。。私はこの方法で直りましたが、一応システム側のファイルを変更することになるのでバックアップを取るなど自己責任でお願いします。
問題のコードは、/usr/lib/python3/dist-packages/guake/guake_app.pyにあります
log.debug("Terminal exited: %s", term) if libutempter is not None: libutempter.utempter_remove_record(term.get_pty()) self.delete_tab(self.notebook.page_num(widget), kill=False, prompt=False)
となっている箇所の、最終行のインデントに問題があるらしいので、その部分のインデントを一段階下げ
log.debug("Terminal exited: %s", term) if libutempter is not None: libutempter.utempter_remove_record(term.get_pty()) self.delete_tab(self.notebook.page_num(widget), kill=False, prompt=False)
のようにします。Githubでは1435行目となっていましたが、私の環境では該当箇所は1420行目にありました。
変更を保存した後Guakeを再起動させると、exitを打ってもフリーズすることなく画面が閉じるはずです。
代替アプリ
直ったので良いのですが、一行だけの変更で直るバグが結構長い間リポジトリで放置されているのがどことなく不安なので、類似のターミナルエミュレータも探してみました。
色々あったような気もしますが、この記事を書いている現在覚えているのは1つだけです(笑)。
それがtilixです。画面分割等の豊富な機能を揃えており、かつキーボードから簡単に呼び出せるようなので、今後困ったらこれを使おうと思います。
使い方はこちらの方の記事を参考にさせていただきました。
公式のリポジトリにあるのでUbuntuならインストールはaptで一発です。
sudo apt install tilix
tilixにはQuakeモードというものがあり、こちらのモードではドロップダウン型のターミナルのように使えます。
Quakeモードで起動するコマンドは次のとおりです。
tilix --quake
とはいえ呼び出すたびにこのコマンドを打つ必要があり、このコマンドを打つためのターミナルはどうするんだという鶏卵問題のような事態が発生するので、便利に使うためにはGnomeのキーボードショートカットに上のコマンドを登録するなどの必要があるかと思います。
私はタイル型ウィンドウマネージャのi3を使っていますので、configファイルに次の行を追記してF1でいつでも呼び出せるように設定しました。
bindsym F1 exec --no-startup-id "tilix --quake"
使用感は、Guakeより若干重いかなという程度、画面分割を使いこなせればこちらも重宝しそうです。
最後までお読みいただきありがとうございました。
Pythonでディレクトリ下のファイルを再帰的にリネームして通し番号を振る
pythonでディレクトリ(フォルダ)下のファイルを一括で変更して通し番号を振るスクリプトを作ったので備忘録を兼ねてまとめます。
Windowsでも動くような気はしますが、テストした環境はUbuntu 18.04 LTSです。
目標
以下のような構成のフォルダがあるとします。
$ tree dir1
dir1/
├── a.txt
└── b
├── b.txt
├── c
│ └── c.txt
└── d
└── d.txt
これを以下のようにリネームするのが今回作りたかったスクリプトです。
$ python rename.py dir1 file $ tree dir1 dir1/ ├── b │ ├── c │ │ └── file_3.txt │ ├── d │ │ └── file_4.txt │ └── file_2.txt └── file_1.txt
要するに、ディレクトリの構造は変えずにファイルだけリネームし、与えられた接頭辞を付けて、通し番号を振りたかったんです。
スクリプト
完成したスクリプトは以下のとおりです。コメントで簡単な解説を入れました。最初は画像ファイルを対象とするつもりだったので変数名がそれに関連したものになってしまっていますが、ファイルの種類にかかわらず動作します。
ただ、globの再帰的なマッチを使用していますので、Python 3.5以降でないと正しく動作しません。
ファイル名の変更はファイルの上書き等によりデータのロストが起こりやすいものですので、お試しになる場合はバックアップを取るなどして自己責任でお願いします。
# rename.py import glob import os import sys indir = sys.argv[1] prefix = sys.argv[2] # ディレクトリ下のファイル名を再帰的に取得する path = os.path.join(indir, "**/*.*") all_imgs = glob.glob(path, recursive=True) all_imgs.sort() # ファイル数が何桁か調べる digit = len(str(len(all_imgs))) for index, img in enumerate(all_imgs, 1): # 元のディレクトリ構造とファイル名、拡張子を取得する ori_dir, ori_name = os.path.split(img) filename, suffix = os.path.splitext(ori_name) # 自分自身はリネームしない if ori_name == "rename.py": continue # 通し番号の桁を揃える num = str(index).zfill(digit) # 新しいファイル名の生成、ディレクトリと拡張子は元のまま new_name = "{0}_{1}{2}".format(prefix, num, suffix) new_path = os.path.join(ori_dir, new_name) # 変更後と同じ名前のファイルが存在する場合はスキップ if os.path.exists(new_path): sys.stderr.write( "{0} already exists\n".format(new_path)) continue # リネーム os.rename(img, new_path) print("{0} -> {1}".format(ori_name, new_name))
Newick形式で系統樹の樹形が同じか確かめる
複数のNewick形式の系統樹から樹形が同じものを抽出する必要に駆られたので備忘録を兼ねたまとめです。
当たり前ですが、
((A,B),(C,D)); ((C,D),(A,B));
の2つのNewick形式は同じ樹形を表しており、
((A,B),(C,D)) ((A,C),(B,D));
の2つのNewick形式は異なる樹形を表しています。
問題は、文字列としては異なる2つのNewick形式で樹形が同一かを判定しなければならない点にあります。枝長を無視して判定したいような場合には、さらに処理がややこしくなります。この辺りで頭が痛くなってきました。
という訳で既存のコードを探します。
調べた結果、Rのapeパッケージにシンプルに樹形比較をしてくれる関数があるようなのでこれを利用します。
インストール方法はごく普通です。
install.packages("ape")
使い方もシンプルです。同じくapeのread.tree関数でNewick形式を読み込んで、all.equal.phylo関数で樹形が同じかを判定します。今回は文字列をそのまま読み込ませていますが、ファイルのパスを渡すこともできます。
> t1 <- read.tree(text='((A,B),(C,D));') > t2 <- read.tree(text='((C,D),(A,B));') > all.equal.phylo(t1, t2) [1] TRUE
> t1 <- read.tree(text='((A,B),(C,D));') > t3 <- read.tree(text='((A,C),(B,D));') > all.equal.phylo(t1, t3) [1] FALSE
枝長を無視して樹形比較を行いたい場合は、use.edge.length引数をFALSEにします。デフォルトではTRUEになっています。
以上をまとめて今回作った、Newick形式の系統樹ファイルを2つ受け取って、樹形が同じかを表示するプログラムを下に載せました。
#!/usr/bin/Rscript
library('ape')
args <- commandArgs(trailingOnly = TRUE)
nwk1 <- args[1]
nwk2 <- args[2]
tree1 <- read.tree(nwk1)
tree2 <- read.tree(nwk2)
is.same.topology <- all.equal.phylo(tree1, tree2, use.edge.length=FALSE)
if (is.same.topology == TRUE) {
cat("identical\n")
} else {
cat("different\n")
}
本当はPythonで作りたかったんですけどね。Rも勉強しないといけないってことですかね。
UbuntuにNeovimをインストールする
この度VimからNeovimに乗り換えたので、備忘録を兼ねてそのインストール手順をまとめます。環境はUbuntu 18.04 LTSです。
正直な所VimとNeovimとの違いは分かっていませんし、プラグインの導入等一通りの設定が終わった今でも、Vimとの使用感の違いはないのですが、何というか新しそうなので使っています(笑)。
ただ、GUIがカッコイイOnivimが開発されているのは魅力だなーと思います。見た目は大事ですからね。
さて、インストールは簡単だったんですが自動補完プラグインのdeopleteを使う関係で二点躓きました。
neovimのバージョンが古い
Ubuntu 18.04ではリポジトリにneovimが追加されているんですが、deopleteが要求するバージョンよりも古いためそのままインストールしてしまうと、エラーが出てしまいます。ということで、ppaを登録してインストールすることにしました。
こちらの記事を参考にneovimのppaを追加しました。
コマンドは次の通りです。
sudo add-apt-repository ppa:neovim-ppa/stable sudo apt-get update sudo apt-get install neovim
pipでneovimではなくpynvimをインストールする
deopleteを使いたい場合はpipでneovimパッケージもインストールしましょう、というのが定石ぽかったんですがエラーが出てしまいました。結構困ってたんですが、原因は単純でした。今はneovimではなくpynvimに名前が変わっているようです。
公式の説明もそう変わっていますし、やっぱり面倒でも最新の情報を追わないとダメですね。参考にさせていただいたのはこちらの記事です。
コマンドは次の通りです。
sudo -H pip3 install pynvim
設定でも色々苦労した記憶はあるのですが、とりあえずここまでにしておきます。UpdateRemotePluginsとCheckHealthすると、案外何とかなったりします。
ただ、終わってしまえばdeopleteは設定が本当に楽で設定ファイルが前と比べかなりスッキリしたので、そこは素直に乗り換えてよかったなと思っています。
まあ、VIｍ8ならdeopleteも使えますけど…。
仕方ないからって訳じゃないんですが、せっかく設定ファイルが分かれているのでVimはプラグイン無し、Neovimはプラグインてんこ盛り、と使い分けています。
最後までお読みいただきありがとうございました。
