按列表文件顺序提取fasta序列

Question

我需要从 "goodProteins.fasta" 文件（第一个输入）中提取一些 fasta 序列，id 列表文件存在于单独的文件夹（第二个输入）中。

fasta序列文件格式为：

>1_12256
FSKVJLKDFJFDAKJQWERTYU......
>1_12257
SKJFHKDAJHLQWERTYGFDFHU......
>1_12258
QWERTYUHKDJKDJOKK......
>1_12259
DJHFDSQWERTYUHKDJKDJOKK......
>1_12260
ADKKHDFHJQWERTYUHKDJKDJOKK......

其中一个 id 文件的格式是：

1_12258
1_12256
1_12257

我正在使用以下脚本：

from Bio import SeqIO
import glob

def process(wanted_file, result_file):
    fasta_file = "goodProteins.fasta" # First input (Fasta sequence)

    wanted = set()
    with open(wanted_file) as f:
        for line in f:
            line = line.strip()
            if line != "":
                wanted.add(line)

    fasta_sequences = SeqIO.parse(open(fasta_file),'fasta')
    with open(result_file, "w") as f:
        for seq in fasta_sequences:
            if seq.id in wanted:
                SeqIO.write([seq], f, "fasta")

listFilesArr = glob.glob("My_folder\*txt") # takes all .txt files as
                                           # Second input in My_folder
for wanted_file in listFilesArr:
    result_file = wanted_file[0:-4] + ".fasta"
    process(wanted_file, result_file)

它应该根据id文件中的信息和顺序列表提取fasta序列，所需的输出将是：

>1_12258
QWERTYUHKDJKDJOKK......
>1_12256
FSKVJLKDFJFDAKJQWERTYU......
>1_12257
SKJFHKDAJHLQWERTYGFDFHU......

但我得到：

>1_12256
FSKVJLKDFJFDAKJQWERTYU......
>1_12257
SKJFHKDAJHLQWERTYGFDFHU......
>1_12258
QWERTYUHKDJKDJOKK......

也就是说，在我的最终输出中，我得到 header 根据它们的较低值排序，但我希望它们的顺序与列表文件中描述的顺序完全相同。我不知道该怎么做...请帮忙。

Answer 1

我认为排序问题的根本原因是因为 wanted 是一个未排序的 set。由于您希望 wanted_files 中的序列 ID 确定顺序，因此您需要将它们存储在保留顺序的其他内容中，例如列表。

或者，您可以只在读取时处理 wanted_file 的每一行。这种方法的一个问题是它可能需要您多次通读 "goodProteins.fasta" 文件 — 如果 wanted_file 的内容未按排序顺序，可能每行一次。

为避免这种情况，可以使用 SeqIO.to_dict() 函数将整个文件读入内存驻留字典，其键为序列 ID，然后对每个 wanted_file 重复使用。你说这个文件是 50-60 MB，但这对于今天的大多数硬件来说并不算太大。

无论如何，这是尝试执行此操作的代码。为了避免全局变量，有一个 Process class 读取 "goodProteins.fasta" 文件并在创建它的实例时将其转换为字典。实例是可调用和可重用的，这意味着相同的过程对象可以与每个 wanted_file 一起使用，而无需重复读取序列文件。

请注意，代码未经测试，因为我的系统上没有安装数据文件或 Bio 模块 — 但希望它足够接近以提供帮助。

from Bio import SeqIO
import glob

class Process(object):
    def __init__(self, fasta_file_name):
        # read entire fasta file into memory as a dictionary indexed by ID
        with open(fasta_file_name, "rU") as fasta_file:
            self.fasta_sequences = SeqIO.to_dict(
                                        SeqIO.parse(fasta_file, 'fasta'))

    def __call__(self, wanted_file_name, results_file_name):
        with open(wanted_file_name, "rU") as wanted, \
             open(results_file_name, "w") as results:
            for seq_id in (line.strip() for line in wanted):
                if seq_id:
                    SeqIO.write(self.fasta_sequences[seq_id], results, "fasta")

process = Process("goodProteins.fasta")  # create process object

# process each wanted file using it
for wanted_file_name in glob.glob(r"My_folder\*.txt"):
    results_file_name = wanted_file_name[:-4] + ".fasta"
    process(wanted_file_name, results_file_name)

按列表文件顺序提取fasta序列

Extracting fasta sequences in list files order

python

file-io

biopython