基于另一个文件中的共享字段，从一个文件中绘制多个序列答案

【问题标题】：Drawing multiple sequences from 1 file, based on shared fields in another file基于另一个文件中的共享字段，从一个文件中绘制多个序列
【发布时间】：2019-01-01 02:20:44
【问题描述】：

我正在尝试运行一个 python 脚本来从一个单独的文件 (merged.fas) 中绘制序列，相对于一个列表 (gene_fams_eggnog.txt) 我作为另一个程序的输出。

代码如下：

from Bio import SeqIO
import os, sys, re
from collections import defaultdict
sequences = "merged.fas"
all_seqs = SeqIO.index(sequences, "fasta")
gene_fams = defaultdict(list)

gene_fams_file = open("gene_fams_eggnog.txt")
for line in gene_fams_file:
     fields = re.split("\t", line.rstrip())
     gene_fams[fields[0]].append[fields[1]]

for fam in gene_fams.keys():
    output_filename = str(fam) + ".fasta"
    outh = open(output_filename, "w")
    for id in gene_fams[fam]:
        if id in all_seqs:
            outh.write(">" + all_seqs[id].description + "\n" + str(all_seqs[id].seq) + "\n")
        else:
            print "Uh oh! Sequence with ID " + str(id) + " is not in the all_seqs file!"
            quit()
    outh.close()

列表如下所示：

1 酿酒酵母_DAA09367.1
1 bieneu_EED42827.1
1 Asp_XP_749186.1
1 Mag_XP_003717339.1
2 Mag_XP_003716586.1
2 Mag_XP_003709453.1
3 Asp_XP_749329.1

字段0表示基于序列之间相似性的分组。该脚本旨在从 merge.fas 中获取与字段 1 中的代码相对应的所有序列，并将它们写入基于字段 0 的文件中。

因此，在我显示的列表部分的情况下，所有在字段 0 中具有 1 的序列（Saccharomycescerevisiae_DAA09367.1、bieneu_EED42827.1、Asp_XP_749186.1、Mag_XP_003717339.1）都将被写入名为 1.fasta 的文件。这应该从 2.fasta 开始——不管有多少组。

所以这已经奏效了，但是它不包括组中的所有序列，它只会包括最后一个被列为该组的一部分的序列。使用上面的示例，我只有一个文件 (1.fasta) 和一个序列 (Mag_XP_003717339.1)，而不是全部四个。

感谢您的任何帮助，谢谢， JT

【问题讨论】：

不要用普通字典定义gene_fams，而是from collections import defaultdict; gene_fams = defaultdict(list)
感谢 Chris 的建议，虽然它似乎仍然只有一个序列。
还应该将outh.close() 移回缩进级别
啊抱歉，我正在运行的实际脚本中没有它，这只是一个错误。不过，我会更改它，以避免进一步混淆。

标签： python bioinformatics biopython fasta

【解决方案1】：

虽然我没有发现您抱怨的问题的原因，但我很惊讶您的代码运行时出现了这个错误：

gene_fams[fields[0]].append[fields[1]]

即append[...] 而不是 append(...)。但也许这也是“我正在运行的实际脚本中不存在”。我在下面重写了你的脚本，它对我来说很好。一个问题是您使用了 Python 内置变量名称 id。您会看到我为了避免此类错误而采取了极端措施：

from Bio import SeqIO
from collections import defaultdict

SEQUENCE_FILE_NAME = "merged.fas"
FAMILY_FILE_NAME = "gene_families_eggnog.txt"

all_sequences = SeqIO.index(SEQUENCE_FILE_NAME, "fasta")
gene_families = defaultdict(list)

with open(FAMILY_FILE_NAME) as gene_families_file:
    for line in gene_families_file:
        family_id, gene_id = line.rstrip().split()
        gene_families[family_id].append(gene_id)

for family_id, gene_ids in gene_families.items():
    output_filename = family_id + ".fasta"

    with open(output_filename, "w") as output:
        for gene_id in gene_ids:
            assert gene_id in all_sequences, "Sequence {} is not in {}!".format(gene_id, SEQUENCE_FILE_NAME)

            output.write(all_sequences[gene_id].format("fasta"))

【讨论】：