如何选择fasta文件中标题的最后一个字符？

Question

我有一个这样的 fasta 文件：

>XP1987651-apple1
ACCTTCCAAGTAG
>XP1235689-lemon2
TTGGAGTCCTGAG
>XP1254115-pear1
ATGCCGTAGTCAA

我想创建一个文件选择以'1'结尾的标题，例如：

>XP1987651-apple1
ACCTTCCAAGTAG
>XP1254115-pear1
ATGCCGTAGTCAA

到目前为止，我创建了这个：

fasta = open('x.fasta')
output = open('x1.fasta', 'w')
seq = ''

for line in fasta:
    if line[0] == '>' and seq == '':
        header = line
    elif line[0] != '>':
        seq = seq + line


for n in header:
    n = header[-1]
    if '1' in n:
        output.write(header + seq)
    header= line
    seq = ''


if "1" in header:
    output.write(header + seq)

output.close()

但是，它不会在创建的新文件中产生任何输出。你能找出错误吗？ 谢谢

Answer 1

一种选择是将整个文件读入一个字符串，然后将re.findall 与以下正则表达式模式一起使用：

>[A-Z0-9]+-\w+1\r?\n[ACGT]+

示例脚本：

fasta = open('x.fasta')
text = fasta.read()
matches = re.findall(r'>[A-Z0-9]+-\w+1\r?\n[ACGT]+', text)
print(matches)

对于您在上面提供的示例数据，将打印：

['>XP1987651-apple1\nACCTTCCAAGTAG', '>XP1254115-pear1\nATGCCGTAGTCAA']

Answer 2

您可以首先获取由'>' 分隔的个人记录列表，然后使用换行符.split('\n', 1) 的单个拆分提取标题和正文

records = [
    line.split('\n', 1)
    for line in fasta.read().split('>')[1:]
]

然后你可以简单地过滤掉不以1结尾的记录

for header, body in records:
    if header.endswith('1'):
        output.write('>' + header + '\n')
        output.write(body)

Answer 3

当您看到匹配的标题行时，您可以非常简单地设置一个标志。

with open('x.fasta') as fasta, open('x1.fasta', 'w') as output:
    for line in fasta:
        if line.startswith('>'):
            select = line.endswith('1\n')
        if select:
            output.write(line)

这避免了将整个文件读入内存；您一次只检查一行。

也许注意到line 将在行尾包含换行符。我选择简单地保留它；有时，如果您使用line = line.rstrip('\n') 修剪它并在必要时将其添加回输出中，事情会变得更容易。