ほぼ中立ブログ

少しだけ百合に偏った雑記ブログ

BiopythonでFASTAファイルの重複をチェック

【スポンサーリンク】

プログラムで自動で配列を集めてきたり、トリミングなどの処理を行ったりした場合に、予期せず全く同じ配列ファイルが違う名前で紛れてしまうことがあります。私のやり方がまずいのかもしれませんが。

系統樹作成などにそうしたファイルを使ってしまうと解析に悪影響が出るため、プログラムで、指定したディレクトリ内のFASTAファイルの内容が他のファイルと重複していないか調べることにしました。

表現がややこしくて申し訳ないのですがヘッダ行の内容は無視し、あるFASTAファイルに記載されている配列の全セットが、他のファイルに記載されている配列セットと完全に一致したときにファイル名を出力するようにしました

完成したコードは以下の通りです。

#!/usr/bin/env python3

from collections import defaultdict
from glob import glob
import os
import sys
from Bio import SeqIO


def check_duplicate(fas_dir):
    path = os.path.join(fas_dir, "*.*")
    fastas = glob(path)
    dup_comb = set()
    for fas1 in fastas:
        seq_dict1 = fasta_to_dict(fas1)
        for fas2 in fastas:
            if fas1 == fas2:
                continue
            comb = frozenset({fas1, fas2})
            if comb in dup_comb:
                continue
            seq_dict2 = fasta_to_dict(fas2)
            if seq_dict1 == seq_dict2:
                dup_comb.add(comb)
                print("{0}\t{1}".format(fas1, fas2))


def fasta_to_dict(fas):
    '''FASTAファイルを入力とし、配列をキー、その個数を値とする辞書を生成する'''
    seq2count = defaultdict(int)
    for record in SeqIO.parse(fas, "fasta"):
        seq = str(record.seq)
        seq2count[seq] += 1
    return seq2count


if __name__ == '__main__':
    indir = sys.argv[1]
    check_duplicate(indir)

一応想定した動作をしたのですが、2重ループでファイルを開いている箇所が重いらしく、ファイル数を増やすと結構時間がかかってしまうため、何か改善策があれば修正したいと思います。

読みづらい文章になってしまいましたが、お読みいただきありがとうございました。

バイオインフォマティクス関連書籍

Dr.Bonoの生命科学データ解析 [ 坊農 秀雅 ]

価格:3,240円
(2018/10/17 23:30時点)
感想(0件)