根据单独文件中的条目从FASTA文件中提取序列

Question

我有两个文件。

文件1：一个FASTA格式的文件，里面是基因序列，格式像这个例子：

>PITG_00002 | Phytophthora infestans T30-4 conserved hypothetical protein (426 nt)
ATGCATCGCTCGGGTTCCGCACGGAAAGCCCAAGGTCTGGGATTACGGGGTGGTGGTCGG
TTACACTTGGAATAACCTCGCAAATTCAGAATCTCTACAGGCTACGTTCGCGGATGGAAC
>PITG_00003 | Phytophthora infestans T30-4 protein kinase (297 nt)
ATGACGGCTGGGGTCGGTACGCCCTACTGGATCGCACCGGAGATTCTTGAAGGCAAACGG
TACACTGAGCAAGCGGATATTTACTCGTTCGGAGTGGTTTTATCCGAGCTGGACACGTGC
AAGATGCCGTTCTCTGACGTCGTTACGGCAGAGGGAAAGAAACCCAAACCAGTTCAGATC
>PITG_00004 | Phytophthora infestans T30-4 protein kinase, putative (1969 nt)
ATGCGCGTGTCTGGTCTCCTTTCAATTCTTGCAGCCACTTTGACCACGGCCCAAGACTAC

文件2：一个简单的文本文件，里面只有基因的识别号，像这样。

PITG_00003
PITG_00005
PITG_00023

文件2中的每个条目在文件1中都有，但文件1中的每个条目不一定在文件2中。我需要从文件1中删除所有不在文件2中的条目。我觉得biopython模块里应该有可以帮我解决这个问题的东西，只是我不知道是什么。例如，我最开始想用SeqIO.parse函数从我的FASTA文件中提取出识别号，但这样做只是让我得到了两个识别号的文件。我不知道怎么选择性地提取出在另一个文件中的识别号。也许可以把文件2中的所有条目读到一个字典里，然后把这个条目和文件1中对应的条目关联起来，再用SeqIO.parse提取整个序列……但我真的不知道……任何人能给我的帮助我都非常感激！

Answer 1

试试这个：

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()

简单说一下这里发生了什么…… accessionids.txt 是你的 文件 2，而 fasta.txt 是你的 文件 1。显然，你需要在代码中把这些文件名替换成你实际使用的文件名。

首先，我们创建一个字典（有时也叫哈希表或关联数组），对于 文件 2 中的每个 Accession ID，我们都会创建一个条目，其中 键 是 Accession ID，值 设置为 1（在这个情况下，值其实并不重要）。

接下来，我们查看 文件 1，再次查看该文件中的每一行。如果文件中的某一行以 > 开头，那么我们就知道这行包含一个 Accession ID。我们把这一行按 | 切分，因为每个包含 Accession ID 的行都会在字符串中有一个 |。然后，取切分后的第一部分，按照 _splitline[0] 来指定。我们使用 accessorIDWithArrow[1:-1] 来去掉字符串的第一个和最后一个字符，这两个字符分别是前面的 > 符号和后面的空格。

到这里，accessorID 现在包含了我们从 文件 2 期待的格式的 Accession ID。

接下来，我们检查之前创建并填充的字典中是否有这个 Accession ID 作为键。如果有，我们就立即把包含这个 Accession ID 的行写入一个新文件 fasta_parsed.txt，并将 skip 这个“标记”变量设置/重置为 0。接下来的 else 语句中包含的 if not skip 部分将允许与我们找到的 Accession ID 相关的后续行被打印到 fasta_parsed.txt 文件中。

对于在字典中找不到的 文件 1 中的 Accession ID（也就是说不在 文件 2 中），我们就不把这一行写入 fasta_parsed.txt，并将 skip 标记设置为 0。因此，直到在 文件 1 中找到另一个在 文件 2 中存在的 Accession ID，所有后续行都会被跳过。