python - 从fasta文件中有选择地选择核苷酸序列？

Question

如果基因名称存储在文本文件中，如何使用biopython从fasta文件中剪切我感兴趣的基因？

#extract genes                
f1 = open('ortholog1.txt','r')
f2 = open('all.fasta','r')
f3 = open('ortholog1.fasta','w')

genes = [line.rstrip('\n') for line in f1.readlines()]


i=0
for seq_record in SeqIO.parse(f2, "fasta"):
    if genes[i] == seq_record.id:              
        print genes[i]
        f3.write('>'+genes[i])
        i=i+1
        if i==18:
            break
        f3.write('\n')
        f3.write(str(seq_record.seq))
        f3.write('\n')


f2.close()
f3.close()

我正在尝试上面的代码。但它有一些错误并且不通用，因为像 ortholog1.txt (包含基因名称)一样，还有 5 个类似的文件。每个文件中的基因数量也各不相同(并非像这里一样总是 18 个)。这里all.fasta是包含所有基因的文件。 ortholog1.fasta 必须包含剪断的核苷酸序列。

Answer 1

基本上，您可以让 Biopython 完成所有工作。

我猜测“ortholog1.txt”中的基因名称与 fasta 文件中的基因名称完全相同，并且每行有一个基因名称。如果没有，您需要根据需要调整它们以使它们对齐。

from Bio import SeqIO

with open('ortholog1.txt','r') as f:
    orthologs_txt = f.read()
orthologs = orthologs_txt.splitlines()

genes_to_keep = []
for record in SeqIO.parse(open('all.fasta','r'), 'fasta'):
    if record.description in orthologs:
        genes_to_keep.append(record)

with open('ortholog1.fasta','w') as f:
    SeqIO.write(genes_to_keep, f, 'fasta')

编辑:这是使输出基因保持与直系同源文件中相同顺序的一种方法:

from Bio import SeqIO

with open('all.fasta','r') as fasta_file:
    record_dict = SeqIO.to_dict(open(SeqIO.parse(fasta_file, 'fasta')

with open('ortholog1.txt','r') as text_file:
    orthologs_txt = text_file.read()

genes_to_keep = []     
for ortholog in orthologs_txt.splitlines():
    try:
        genes_to_keep.append( record_dict[ortholog] )
    except KeyError:
        pass

with open('ortholog1.fasta','w') as output_file:
    SeqIO.write(genes_to_keep, output_file, 'fasta')