はじめに
皆様は、アラインメントの残基番号って、どうやって数えていますか？
例えば、Seq1 の n 番目の残基が Seq2 の何番目にアラインされるか知りたい場合です。
Clustal Omega
このように折り返した端っこにしか残基番号が書いてなかったり……。
ESPript
アラインメントの作図によく使われる ESPript も一番上の配列の番号しか表示されません。でも、Seq1 の D49 にアラインされる Seq2 の D の残基番号を知りたい、みたいなケースがよくあります。
私は前後の配列を手打ちして検索し、Word の文字数カウントで数えたりしていました。
が、注目したい残基が複数 × 複数の配列 になると作業が超大変になります。
アラインメントファイルとは
そもそも、アラインメントのファイルの中身はどうなっているのでしょうか。
この配列のアラインメントを MEGA11 で作成し、
Data > Export Alignment > FASTA format
を選択して FASTA 形式で保存します。
メモ帳で開いて見てみると、通常の FASTA と同じく、> で始まり改行で終わるヘッダー行があり、ギャップが - で表されています。つまり、配列部分を等幅フォントで並べて表示してやれば勝手にアラインメントに見える、という仕組みです。
折り返さないで表示するとこんな感じ。
配列の - 以外の文字をカウントすれば残基番号を数えることができそうです。
すべての残基番号を表示し、csv ファイル（コンマ区切りの表ファイル）として出力してみることにしました。
完成図
Python にやらせる
プログラミング超初心者の時に作ったのですが、for 文、カウント、文字連結といった教科書の最初の方に出てくるような知識だけでできるので初心者の練習にピッタリだなと、初心者ながらに思いました（笑
import re
正規表現のモジュールを使用します。
fas = 'C:/Users/xxx/yyy/test.fas'#FASTAで保存したアラインメントのファイル
output = []
with open(fas, "r") as f:#ファイルを開く
fa = f.readlines()#リストとして読み込み
for i in fa:
m = re.match('>', i)#文頭の > を検索
if m:#マッチしたら
output.append(i.replace('\n', ''))#その行から改行を削除しリストoutputに格納
else:
num = ['number']
letters = ['']
count = 0# 変数 count に 0 を代入
for j in i:# for 文で一つひとつの文字を…
letters.append(j)# letters に格納
m2 = re.match('[A-Z]', j)#大文字アルファベットか判定
if m2:#マッチしたら
count += 1# count に 1 足して
num.append(str(count))# その値をリスト num に格納
else:
num.append('-')#ハイフンをリスト num に格納
output.append(','.join(num))# リストの要素を , で連結し、outputに格納
output.append(','.join(letters))# 文字のリストも同じく
保存します。
out_file = fas.replace('.fas', '.csv')#入力ファイルの拡張子を .csv にして出力ファイル名とする。
print(out_file)
with open(out_file, 'w') as r:
for d in output:
r.write("%s\n" % d)
Excel で開く
出力したファイルを Excel で開き、セルを全選択して
ホーム > 書式 > 列の幅の自動調整
で見やすくなります。
表示 > ウィンドゥ枠の固定 > 先頭列の固定
で、ヘッダーを常に表示しておくことができます。
まとめ
アラインメントの残基番号をすべて表示した csv ファイルを作成することができました。
残基番号の確認に使うだけでなく、Excel で開いて直接セルに色を付けたり、コメントを書き込むという使い方もできるので、個人的に重宝しています。
