使用条件匹配一行中的多个模式

Question

我有一个像这样的fasta文件：myfasta.fasta

>1_CDS
AAAAATTTCTGGGCCCCGGGGG
AAATTATTA
>2_CDS
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>3_CDS
TTTGGGAATTAAACCCT
>4_CDS
TTTGGGAATTAAACCCT
>5_rRNA
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>6_tRNA
TTAAAAATTTCTGGGCCCCGGGAAAAAA

我有一个代码，我想用它来根据它们的ID来分离序列，这些ID具有匹配的模式，如'CDS'，'tRNA'等。在下面的代码中，我试图使用startswith并且匹配模式在线，这不是'似乎工作。有人可以请帮助我如何在python中查找两个条件。

代码：python mycode.py myfasta.fasta

#!/usr/bin/env python
import sys
import os
myfasta = sys.argv[1]
fasta = open(myfasta)

for line in fasta:
    if line.startswith('>') and 'CDS' in line:
        print(line)
    else:
        print(line)

预期输出（如果我使用CDS）：

>1_CDS
AAAAATTTCTGGGCCCCGGGGG
AAATTATTA
>2_CDS
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>3_CDS
TTTGGGAATTAAACCCT
>4_CDS
TTTGGGAATTAAACCCT

Answer 1

这是一个适合您的代码。如果一行有CDS，它会打印该行和下一行。 strip()在打印行时删除结束字符。

#!/usr/bin/env python
import sys
import os
myfasta = sys.argv[1]

flag = False
with open(myfasta) as fasta:
    for line in fasta:
        if line.startswith('>') and 'CDS' in line:
            flag = True
        elif line.startswith('>'):
            flag = False
        if flag:
            print(line.strip())

编辑：您可以删除elif部分，如下代码：

#!/usr/bin/env python
import sys
import os
myfasta = sys.argv[1]

flag = False
with open(myfasta) as fasta:
    for line in fasta:
        if line.startswith('>'):
            flag = 'CDS' in line
        if flag:
            print(line.strip())

Answer 2

Maanijou的答案很好。

另外，请考虑使用迭代器替代。

EDIT: Updated the code based on your comments

#!/usr/bin/env python
import sys
import os
myfasta = sys.argv[1]
fasta = open(myfasta, "r+")

file_contents = iter(fasta)

try:
    print_flag = True
    while True:
        line = file_contents.next()
        if line.startswith('>'):
            if "CDS" in line:
                print (line.strip())
                print_flag = True
            else:
                print_flag = False
        else:
            if print_flag:
                print (line.strip())

except StopIteration:
    print ("Done")
    fasta.close()

Explanation

file_contents = iter(fasta)将可迭代文件对象转换为迭代器，您可以在其上简单地继续调用next()，直到您用完所有内容为止

为什么我不建议调用readlines，因为其他一些答案是有时fasta文件可能很大并且调用readlines消耗大量内存。

如果一行满足你的搜索需求你只需打印它和下一行，如果不是你只是阅读下一行而什么都不做，

Explanation for Update

1. 由于文件模式，你得到属性错误，我无法在本地重现它，但我认为用正确的模式打开文件应修复它
2. 你现在说可能有超过1个基因组序列为CDS更新代码打印文件中的1个CDS头的所有基因组序列

我用修改过的fasta文件测试了它

>1_CDS
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGCG
>2_CDS
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>3_CDS
TTTGGGAATTAAACCCT
>4_CDS
TTTGGGAATTAAACCCT
>5_rRNA
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>6_tRNA
TTAAAAATTTCTGGGCCCCGGGAAAAAA

而这个输出

python fasta.py fasta.fasta
>1_CDS
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGGG
AAAAATTTCTGGGCCCCGGGCG
>2_CDS
TTAAAAATTTCTGGGCCCCGGGAAAAAA
>3_CDS
TTTGGGAATTAAACCCT
>4_CDS
TTTGGGAATTAAACCCT
Done

Answer 3

这是你想要的吗？

#!/usr/bin/env python
import sys
import os
from collections import defaultdict

myfasta = sys.argv[1]
with open(myfasta) as fasta:
    data = fasta.read().splitlines()

pattern_data = defaultdict(list)
index = 0
while index < len(data):
    if data[index].startswith('>'):
        start = data[index].index('_') + 1
        key = data[index][start:]
        pattern_data[key].append(data[index + 1])
    index += 2

此时，您可以随意使用已排序的数据执行任何操作。

以上假设您解析的整个文件遵循上面显示的确切格式：1行以“>”开头，id是后面的单行。如果您有多行，则代码需要稍作修改。

编辑：我刚刚阅读了fasta文件。我现在知道它们实际上可能在识别后具有长于一行的序列。因此，需要修改上述代码以考虑多行序列。更通用的方法如下：

#!/usr/bin/env python
import sys
import os
from collections import defaultdict

myfasta = sys.argv[1]
with open(myfasta) as fasta:
    data = fasta.read().splitlines()

id_line_indices = [index for index, line in enumerate(data) if line.startswith('>')]
id_line_indices.append(len(data))
pattern_buckets = defaultdict(list)

i = 0
while i < len(id_line_indices) - 1:
    start = data[id_line_indices[i]].index('_') + 1
    key = data[id_line_indices[i]][start:]

    sequence = [data[index] for index in range(id_line_indices[i] + 1, id_line_indices[i + 1])]
    sequence = ''.join(sequence)

    pattern_buckets[key].append(sequence)
    i += 1

对于上述数据集，这仍然实现了相同的结果。例如，

print(pattern_buckets['CDS'])
print(pattern_buckets['rRNA'])

会得到你：

['AAAAATTTCTGGGCCCCGGGGG', 'TTAAAAATTTCTGGGCCCCGGGAAAAAA', 'TTTGGGAATTAAACCCT', 'TTTGGGAATTAAACCCT']
['TTAAAAATTTCTGGGCCCCGGGAAAAAA']