使用 Biopython 从 ID 列表中提取 fasta 文件答案

【问题标题】：Extract fasta files from ID list with Biopython使用 Biopython 从 ID 列表中提取 fasta 文件
【发布时间】：2021-12-28 07:30:43
【问题描述】：

我正在使用 Biopython 在 fasta 文件中查找与包含所选 ID 的 .txt 文件中的 ID 匹配的序列。在 fasta 文件中手动搜索 ID 名称时，我确实得到了命中，但以下脚本没有找到/提取任何序列：

#!/usr/bin/env python3
from Bio import SeqIO
wanted_ids = "transcript.orthogroup7.txt"
input_filename = "hq_isoseq_transcripts.fasta"
output_filename = "wanted_hq_isoseq_transcripts.fasta"
count = 0
total = 0
output_handle = open(output_filename, "w")
for record in SeqIO.parse(input_filename, "fasta"):
    total = total + 1
    if record.id in wanted_ids:
        count = count + 1
        SeqIO.write(record, output_handle, "fasta")
output_handle.close()
print(str(count) + " records selected out of " + str(total))

这是 .fasta 文件的摘录

>20190807-P1_PBI19-RO01_8pM_HQ_transcript/0
GGTACGTATAGGAATTATAAGGGCTCATTGAGACAGTCTGAGCAAAAAAGTCTCGTTCAG
GCTCGCTAAAAGCCTTTGGTGACGCATCTGTCCGTGAGCGAGTCACATTTTGTGAGCGAG
TCACGAATCTCGGCGCCGGCGAATCATGAATCTCCGTAACAGTGAGTCGCGAATCTCCTA
AGCCGTGTGTTTTGCGCACGTCACTTCATTCCTCGGTACAGGGCTTCCTCTTAACAGGAA ... and so on
>20190807-P1_PBI19-RO01_8pM_HQ_transcript/1
GACTATATAGGAATTATAAGGGCTCATTGAGACAGTCTGAGCAAAAAAGTCTCGTTCAGG
CTCGCTAAAAGCCTTTGGTGACGCATCTGTCCGTGAGCGAGTCACATTTTGTGAGCGAGT
CACGAATCTCGGCGCCGGCGAATCATGAATCTCCGTAACAGTGAGTCGCGAATCTCCTAA
GCCGTGTGTTTTTGCGCACGTCACTTCATTCCTCGGTACAGGGCTTCCTCTCAGCAGGAA ... and so on

还有我的 ID .txt 文件

>20190807-P1_PBI19-RO01_8pM_HQ_transcript/0
>20190807-P1_PBI19-RO01_8pM_HQ_transcript/1

我在 conda 环境中执行脚本并收到以下意外消息

(biopython_transcripts) $ python3 script.py
0 records selected out of 35629

将“计数”和“总数”更改为例如 1000 时，我得到

(biopython_transcripts) $ python3 script.py
1000 records selected out of 35629

所以脚本运行但没有提取任何序列，尽管我确信序列存在并且 .txt 和 .fasta 文件中的名称完全匹配。

我对 Biopython 以外的其他解决方案持开放态度，但非常感谢帮助理解为什么这不起作用以及可能的解决方案。

【问题讨论】：

这能回答你的问题吗？ check if a line exists in a file python
您正在检查 id 是否包含在 filename 中，而它们显然不是。

标签： python bioinformatics biopython fasta

【解决方案1】：

您需要使用所需的 id 解析文件并构建要保留的 id 列表。目前，您的代码会根据所需 id 的 文件名 检查 id（即，仅当 id 是 transcript.orthogroup7.txt 的子字符串时才匹配）

这是对您的代码的修改：

from Bio import SeqIO
wanted_filename = "transcript.orthogroup7.txt"           ## changed name
input_filename = "hq_isoseq_transcripts.fasta"
output_filename = "wanted_hq_isoseq_transcripts.fasta"
count = 0
total = 0
output_handle = open(output_filename, "w")

with open(wanted_filename) as f:                         ## added lines to
    wanted_ids = [l[1:-1] for l in f]                    ## build list of ids
## now wanted_ids = ['20190807-P1_PBI19-RO01_8pM_HQ_transcript/0', '20190807-P1_PBI19-RO01_8pM_HQ_transcript/1']

for record in SeqIO.parse(input_filename, "fasta"):
    total = total + 1
    if record.id in wanted_ids:
        count = count + 1
        SeqIO.write(record, output_handle, "fasta")
output_handle.close()
print(str(count) + " records selected out of " + str(total))

【讨论】：