python - 从一个大的 fasta 文件中提取特定的 fasta 序列

Question

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

transcripts.txt 文件包含我要从 assembly.fasta 导出到 selected_transcripts 的列表转录本 ID(包括 ID 和序列)。斋戒。例如:

1. 成绩单.txt:

   Transcript_00004|5601Transcript_00005|5352

2. assembly.fasta:

   >Transcript_00004|5601GATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT>Transcript_00004|5360CGATCTGGCGCTGAGCTGGGTGCTGATCGACCCGGCGTCCGGCCGCTCCGTGAACGCCTCGAGTCGGCGCCCGGTGTGCGTTGACCGGAGATCGCGATCTGGGGAGACCGTCGTGCGGTT

The IDs are preceded by the > symbol: >Transcripts_00004|5601.

I have to read the assembly.fasta file, if the transcript ID in assembly.fasta is the same of that write in transcripts.txt, I have to write this transcript ID and its sequence in selected_transcripts.fasta. So, in the example above, I have to write only the first transcript.

Any suggestions?Thanks.

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

Answer 1

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

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

此外，如果 assembly.fasta 在标题行中没有 >，您将不会返回任何包含 ID 和序列的记录。下面的代码应该解决这些问题，假设您将 > 添加到 header 并且 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 进行比较时字符串中出现换行符。