通过文件头中的ID号从fasta文件中提取序列

Question

我有一个带有多个序列的fasta文件，其标题看起来像这样：

>1016BSA34080.1
MTHSVRIITVTVNFLQHRFFIDYMSEIGLLDGEIEQMVSALQEQVHIVARARTLPEMKNLERDTHVIVKT
LKKQLTAFHSEVKKIADSTQRSRYEGKHQTYEAKVKDLEKELRTQIDPPPKSVSEKHMEDLMGEGGPDGS
GFKTTDQVLRAGIRIQNDA

>1038BSA81955.1
MQQQQARRRMEEPTAAAATASSTTSFAAQPLLSRSVAPQAASSPQASARLAESAGFRSAAVFGSAQAAVG
GRGRGGFGAPPGRGGFGAPPAAGFGAAPAFGAPPTLQAFSAAPAPGGFGAPPAPQGFGAPRAAGFGAPPA
PQAFSAVAPASSTAIPLDVTTYLGDTFGSAPTRGPP

标题开头的4位数字是序列的唯一ID。

[您能帮我写一个Python脚本来提取4位ID的序列吗（在一个文本文件中，每行一个ID）？

我尝试修改此脚本（我在此网站上找到：Extract sequences from a FASTA file based on entries in a separate file）以适应我的目的（徒然）：

f2 = open('accessionids.txt','r')
f1 = open('fasta.txt','r')
f3 = open('fasta_parsed.txt','w')

AI_DICT = {}
for line in f2:
    AI_DICT[line[:-1]] = 1

skip = 0
for line in f1:
    if line[0] == '>':
        _splitline = line.split('|')
        accessorIDWithArrow = _splitline[0]
        accessorID = accessorIDWithArrow[1:-1]
        # print accessorID
        if accessorID in AI_DICT:
            f3.write(line)
            skip = 0
        else:
            skip = 1
    else:
        if not skip:
            f3.write(line)

f1.close()
f2.close()
f3.close()

我是Python的新手，我们将不胜感激！谢谢-Divya

Answer 1

accessionids.txt仅包含四位数的代码吗？

如果是这样，将accessorID更改为：

accessorID = accessorIDWithArrow[1:5]

一些使它变得更Pythonic的方法是：

对于AI_DICT，使用集合而不是字典，使用strip()而不是切片来删除换行符，并使用生成器表达式来构建集合

AI_SET = set((line.strip() for line in f2))

True使用False和skip，而不是0和1。

我将这样重做主循环：

in_accession_ids = False
for line in f1:
    if line[0] == '>':
        _splitline = line.split('|')
        accessorIDWithArrow = _splitline[0]
        accessorID = accessorIDWithArrow[1:5]
        # print accessorID
        in_accession_ids = accessorID in AI_SET
    if in_accession_ids:
        f3.write(line)

我认为这种方式的逻辑更加明显。同样，从原始文件的skip = 0或我的文件中的in_accession_ids=True开始，这意味着您将在查找第一个序列标头之前先打印所有内容。那可能是你想要的，那可能不是-我以为我不在重写中。

您可能最终想研究一下Biopython集合-对于此特定任务而言，它有些过分，但总体而言还不错。许多用于读取FASTA文件和相关格式的工具，等等。

http://biopython.org/wiki/Biopython

Answer 2

使用Biopython，您可以这样做（需要安装biopyhton）：

from Bio import SeqIO

f1 = "fasta.fa"
f2 = "accessionids.txt"
f3 = "selected_seqs.fa"
selected_seqs = list()

with open(f2, "r") as seq_ids:
    accessionids = [line.rstrip("\n") for line in seq_ids]

for seq_record in SeqIO.parse(f1, "fasta")
    header = seq_record.name # (or .id or so)
    for accession_id in accessionids:
        if accession_id == header[0:4]:
            selected_seqs.append(seq_record)


SeqIO.write(selected_seqs, f3, "fasta")

这将遍历您的序列记录（fasta文件），并为每个条目检查是否存在与accessionids文件中的id相匹配的内容。

注意：

• seq_record可以具有不同的标签，请检查您的标识符位于哪个标签中。
• 如果您的ID不在开头（并且对于单个序列头来说是唯一的，则只需使用if accession_id in header:
• 有关SeqIO的更多详细信息，请参见biopython tutorial第5章。>>