从一个大的 fasta 文件中提取特定的 fasta 序列答案

【问题标题】：Extract specific fasta sequences from a big fasta file从一个大的 fasta 文件中提取特定的 fasta 序列
【发布时间】：2016-05-13 17:30:30
【问题描述】：

我想使用以下脚本从一个大的 fasta 文件中提取特定的 fasta 序列，但输出为空。

transcripts.txt 文件包含我要从assembly.fasta 导出到selected_transcripts.fasta 的列表转录ID（ID 和序列）。例如：

成绩单.txt：成绩单_00004|5601 Transcript_00005|5352

assembly.fasta:

>成绩单_00004|5601
GATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT
>成绩单_00004|5360
CGATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT

ID 前面有 > 符号：>Transcripts_00004|5601。

我必须读取assembly.fasta 文件，如果assembly.fasta 中的成绩单ID 与transcripts.txt 中写入的相同，我必须在selected_transcripts.fasta 中写入此成绩单ID 及其序列。所以，在上面的例子中，我只需要写第一个成绩单。

有什么建议吗？谢谢。

from Bio import SeqIO

my_list = [line.split(',') for line in open("/home/universita/transcripts.txt")]

fin = open('/home/universita/assembly.fasta', 'r')
fout = open('/home/universita/selected_transcripts.fasta', 'w')

for record in SeqIO.parse(fin,'fasta'):
    for item in my_list:
        if item == record.id:
            fout.write(">" + record.id + "\n")
            fout.write(record.seq + "\n")

fin.close()
fout.close()

【问题讨论】：

见biostars.org/p/68718
您能否edit 提出您的问题并包括transcripts.txt 的一部分以及assembly.fasta 的一部分，以便我们可以处理一些数据？
您在每个冒号之后拆分您的成绩单行，但它是空格分隔的。这是故意的吗？
您的示例 FASTA 文件的标题行中没有 >。
对不起，在 trascripts.txt 中，每个 ID 都是文件的一行。所以示例中的第一个成绩单在第一个泳道，第二个成绩单在第二个泳道，依此类推。

标签： python biopython

【解决方案1】：

根据您的示例，有几个小问题可以解释为什么您什么也得不到。您的 transcripts.txt 在一行中有多个条目，因此 my_list 将在 my_line[0] 中包含 first_line 的所有项目，在您的循环中，您逐行遍历 my_list，因此您的第一项将是

['Transcript_00004|5601', 'Transcript_00005|5352']

此外，如果assembly.fasta 在标题行中没有>，您将不会返回任何带有 ID 和序列的记录。下面的代码应该可以解决这些问题，假设您将> 添加到标题中并且split 函数现在使用的是空格而不是冒号。

from Bio import SeqIO

my_list = []
with open("transcripts.txt") as transcripts:
    for line in transcripts:
        my_list.extend(line.split(' '))

fin = open('assembly.fasta', 'r')
fout = open('selected_transcripts.fasta', 'w')

for record in SeqIO.parse(fin,'fasta'):
    for item in my_list:
        if item.strip() == record.id:
            fout.write(">" + record.id + "\n")
            fout.write(record.seq + "\n")


fin.close()
fout.close()

对成绩单的读取进行了更改，以便将所有 ID 分别附加到 my_list。此外，每个项目都被去除空白，以避免在与record.id 进行比较时字符串中出现换行符。

【讨论】：