fastaファイルに含まれる配列をアライメントして、距離行列を得るプログラム(python)

こんにちは。早雲です。

この記事はすこし専門的な内容になっており、『生物学に特化した短編小説とエッセイ』とは別のマガジンに組み込まれています。

遺伝学とプログラミングに興味がある方はぜひ試してみてください。

・fastaファイルに含まれる配列をアライメントして、距離行列を得る

標題の通り、配列のアライメントをとって距離行列を得るコードを書きましたので、共有します。

すでにアライメントされたfastaファイルでpython(biopython)で系統樹を描くコードの説明は日本語でも結構あったのですが、アライメントの実施と距離行列を得るコードを説明しているサイトが少なかったので苦労しました。

参考にしたのは以下のサイトです。

Bio.Phylo.TreeConstruction.DistanceCalculator

BioPythonで系統樹を書いてNewick出力 - ピエール瀧になりたい

・アライメントツール：muscleのインストール（これが地味に苦労した…）

まずはアライメントツールのmuscleを入れます。biopythonのアライメントはあくまでツールを呼び出し実行するというものなので、ツールをローカル（自分のパソコン）にダウンロードしなければなりません。

参考にしたのは次のサイトです。

井上 潤：MUSCLE

muscle installation

muscleのダウンロードは以下のページから行いました。

MUSCLE download

僕はwindowsユーザなので、windows用のファイルをダウンロードします。
muscle3.8.31_i86win32.exe(横にWindows Intel i86 32/64と書いてあります)
をクリックするとダウンロードが開始されます。

これは自分のパソコンのDownloadのフォルダに保存されます（もし場所がわからないときは、ウィンドウズのタスクバーの検索用ボックスに「ダウンロード」といれてみてください）

そして、Downloadのフォルダからmuscle3.8.31_i86win32.exeを見つけだし、任意のフォルダにコピーします。そのとき、名前をmuscle.exeに変更してください。

で、そのままじゃこのツールをpythonから呼び出しずらいです。pythonから呼び出すときはsubprocessという関数を使うのですが、パスをいちいち指定しなければいけないのです。なのでパスを通します。

パスの通し方は以下のサイトを参考にしました。

Windowsの環境パスを通す(path)

windows 10では以下の手順になります。

スタート→コントロールパネルを検索→システムとセキュリティ→システム→システムの詳細設定→システムの環境変数→pathを選択→編集→新規

で、この新規をおしたら、入力ボックスに先ほどmuscle.exeを保存したディレクトリ（フォルダ）のパスを入力してください（エクスプローラの上の　PC　>　デスクトップ　>　…　と書いてある部分をクリックするとパスが表示されますのでコピペしてください）

そして、OKを押すとパスが保存されます。その後、パソコンを再起動をしてください。

これでmuscleをインストールしました。

・そしてpythonコード

なかなか苦労した割にはコードは少ないです。

まず必要なモジュールをインポート

import sys
import subprocess
from Bio.Align.Applications import MuscleCommandline
from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio import AlignIO

そしてインプットするfastaファイルを指定します。また、アウトプット先も決めておきましょう。

input_path = 'C:/Users/hoge/Desktop/test.fasta'
out_path = 'C:/Users/hoge/Desktop/test.tsv'

そして、muscleを実行します。

muscle_cline = MuscleCommandline(input=input_path)
child = subprocess.Popen(str(muscle_cline),
                        stdout=subprocess.PIPE,
                         stderr=subprocess.PIPE,
                         universal_newlines=True,
                         shell=(sys.platform!="win32"))

ここでstr(muscle_cline)の中身は'muscle -in C:/Users/hoge/Desktop/test.fasta'という文字列です。

つまりsubprocessでコマンドプロンプト(windowsなんでね）の命令としてmuscle - in …　を実行してくれている、という事です。

これに気づくのに大分骨を折りましたが…

次に、アライメントされたデータをファスタ形式で読み込み、alignという変数に格納します。

align = AlignIO.read(child.stdout, "fasta")

そして、距離行列を算出します。

calculator = DistanceCalculator('identity')
dm = calculator.get_distance(align)

出力します。

with open(out_path,"w") as f:
   f.write(str(dm))

上記のコードをまとめたのがこちらです。

import sys
import subprocess
from Bio.Align.Applications import MuscleCommandline
from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio import AlignIO
input_path = 'C:/Users/hoge/Desktop/test.fasta'
out_path = 'C:/Users/hoge/Desktop/test.tsv'
muscle_cline = MuscleCommandline(input=input_path)
child = subprocess.Popen(str(muscle_cline),
                        stdout=subprocess.PIPE,
                         stderr=subprocess.PIPE,
                         universal_newlines=True,
                         shell=(sys.platform!="win32"))
align = AlignIO.read(child.stdout, "fasta")
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(align)
with open(out_path,"w") as f:
   f.write(str(dm))

いや、それにしてもほんと、こんなマニアックな記事をいったいどれだけの人が見るのやら、という気がしています（だってbiopythonでやらなくてもいいことだし…）。ですが、参考になれば幸いです。

ではまた。