将 FASTA 文件中的序列提取到多个文件中，文件基于单独文件中的 header_IDs

Question

我正在寻找一种 python 解决方案，根据与列表的匹配，将 multiple 序列从 FASTA 文件提取到 multiple 文件中标题 ID 在单独的文件中。

这是Extract sequences from a FASTA file based on entries in a separate file 和https://www.biostars.org/p/2822/ 上发布的问题的稍微复杂的版本，它只为所有匹配输出一个文件。

我是 python 新手，正在尝试找到一种方法：

• 获取包含将在 fasta 标头中的字符串的文件
  
• 将与字符串匹配的所有记录写入单独的 fasta 文件
  

header_ID_strings 文件如下所示：

CAP357_2030_09WPI、CAP357_2040_11WPI、CAP357_2050_13WPI等……

我的 fasta 文件示例如下所示：

>CAP357_2030_009wpi_v1v3_1_056_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGG
>CAP357_2040_011wpi_v1v3_1_008_00006_001.1
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
>CAP357_2040_011wpi_v1v3_1_030_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
>CAP357_2040_011wpi_v1v3_1_004_00001_000.2
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
>CAP357_2050_013wpi_v1v3_1_047_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT

预期输出

文件 1：CAP357_2030_009wpi_v1v3.fasta
>CAP357_2030_009wpi_v1v3_1_056_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGG

文件2：CAP357_2040_011wpi_v1v3.fasta
>CAP357_2040_011wpi_v1v3_1_008_00006_001.1
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
>CAP357_2040_011wpi_v1v3_1_030_00002_000.4
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
>CAP357_2040_011wpi_v1v3_1_004_00001_000.2
GTAAAATTAACCCCACTCTGTGTCACTCTAAATTGTACAACTGCAAAGGGT
等等……

此代码来自上面的链接，但我想要：
* 匹配写入单独的输出文件
* 如果可能，我不必单独指定每个输出文件（我将有多达 30 个输出文件）

#!/usr/bin/env python

import sys
from Bio import SeqIO

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file = sys.argv[3]

wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(id_file))
print "Found %i unique identifiers in %s" % (len(wanted), id_file)

index = SeqIO.index(input_file, "fasta")
records = (index[r] for r in wanted)
count = SeqIO.write(records, output_file, "fasta")
assert count == len(wanted)

print "Saved %i records from %s to %s" % (count, input_file, output_file)

到目前为止，这是我想出的（下面的脚本），但不知道如何手动指定所有输出文件和变量（我这里只包括了三个）

from Bio import SeqIO
import pandas as pd
import sys

input_file = sys.argv[1]
id_file = sys.argv[2]
output_file2020 = sys.argv[3]
output_file2030 = sys.argv[4]
output_file2040 = sys.argv[5]

colnames = ["2020", "2030", "2040"]
headerlist = pd.read_csv(id_file, names = colnames, header = None)
infile = list(SeqIO.parse(input_file, "fasta"))
2020_seq = tuple(headerlist.2020)
2030_seq = tuple(headerlist.2030)
2040_seq = tuple(headerlist.2040)

count2020 = 0
count2030 = 0
count2040 = 0
for record in infile:
    if record.id in 2020_seq:
        SeqIO.write([record], output_file2020, "fasta")
        countSU += 1
    elif record.id in PI_seq:
        SeqIO.write([record], output_file2030, "fasta")
        countPI += 1
    elif record.id in REC_seq:
        SeqIO.write([record], output_file2040, "fasta")
        countREC += 1
    else:
        print("no matches found")

print("number of SU is", count2020)
print("number of PI is", count2030)
print("number of REC is", count2040)

Answer 1

为了完整起见，这里是“最终”脚本：

#!/usr/bin/env python
# a script to extract fasta records from a fasta file to multiple separate fasta files based on a particular ID (time point) in a particular field, for a given delimiter
# to run, navigate to file location with command prompt and enter: python split_fasta_by_collections.py infile.fasta
from Bio import SeqIO
import os
import sys

records = SeqIO.parse(sys.argv[1], "fasta")
collected = {}
for record in records:
    descr = record.description.split("_")[1] # "_" sets the delimeter, "1" sets the field where counting starts at 0 for the first field
    try:
        collected[descr].append(record)
    except KeyError:
        collected[descr] = [record ,]

file_name = "outfile%s.fasta" 
file_path = os.getcwd() #sets the output file path to your current working directory

for (descr, records) in collected.items():  
    with open(os.path.join(file_path, file_name % descr), 'w') as f:
        SeqIO.write(records, f, 'fasta')

Answer 2

几个简短的建议：

如果您的所有标题都遵循相同的模式，那么您可以提取唯一元素：

record.description.split("_")[1] 

（从“CAP357_2040_011wpi_v1v3_1_008_00006_001.1”得到“2040”）

如果您使用 dict，您可以组装记录集合：

collected = {}
for record in records:
    descr = record.description.split("_")[1]
    try:
        collected[descr].append(record)
    except KeyError:
        collected[descr] = [record ,]

然后您可以将每个集合写入一个新文件：

file_name = "outfile%s" 
for (descr, records) in collected.items():   # iteritems in python2
    with open(os.path.join(file_path, file_name % descr), 'w') as f:
        SeqIO.write(records, f, 'fasta')