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

Question

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

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

我有一个代码，我想用它来根据它们的I来分离序列，这些I具有“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的回答很好。

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

编辑:根据您的注释更新代码

#!/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()

解释

file_contents = iter(fasta)将可迭代文件对象转换为迭代器，您只需继续调用next()，直到用完要读取的内容为止。

我不建议像其他一些答案那样调用readlines的原因是，有时候fasta文件可能很大，调用readlines会消耗大量内存。

如果一条线满足你的搜索要求，你只需打印它和下一行，如果不是，你只需读下一行，什么也不做，

对更新的解释

1. 由于文件模式导致属性错误，我无法在本地复制它，但我认为使用正确的模式打开文件应该修复它。
2. 您现在已经说过，CDS可能有超过一个基因组序列，更新了代码，以打印文件中的一个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']