0

我有一个包含多个序列的 fasta 文件，其标题如下所示：

>1016BSA34080.1
MTHSVRIITVTVNFLQHRFFIDYMSEIGLLDGEIEQMVSALQEQVHIVARARTLPEMKNLERDTHVIVKT
LKKQLTAFHSEVKKIADSTQRSRYEGKHQTYEAKVKDLEKELRTQIDPPPKSVSEKHMEDLMGEGGPDGS
GFKTTDQVLRAGIRIQNDA

>1038BSA81955.1
MQQQQARRRMEEPTAAAATASSTTSFAAQPLLSRSVAPQAASSPQASARLAESAGFRSAAVFGSAQAAVG
GRGRGGFGAPPGRGGFGAPPAAGFGAAPAFGAPPTLQAFSAAPAPGGFGAPPAPQGFGAPRAAGFGAPPA
PQAFSAVAPASSTAIPLDVTTYLGDTFGSAPTRGPP

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

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

我尝试修改此脚本（我在此网站上找到：根据单独文件中的条目从 FASTA 文件中提取序列）以适合我的目的（徒劳）：

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

4

2 回答 2

1

是否accessionids.txt只包含四位代码？

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

accessorID = accessorIDWithArrow[1:5]

一些使它更 Pythonic 的方法是：

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

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

使用TrueandFalse而不是 0 和 1 表示skip.

我将重做主循环：

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

于 2013-08-07T22:18:34.117 回答
1

使用 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 教程第 5 章。
于 2020-06-13T15:52:41.360 回答