通过标题中的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]

使其更具Python化的一些方法是：

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

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

对skip使用True和False，而不是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可以有不同的标签，请检查您的标识符所在的标签。
• 如果您的标识不在开头(并且对于单个序列标头是唯一的)，则只需使用if accession_id in header:
• See biopython tutorial第5章即可了解有关SeqIO的更多详细信息。