Найти все узоры в файле мультифаста, включая перекрывающиеся мотивы. ⇐ Python

[Anonymous]

У меня есть мультифастовый файл, он выглядит так:

Код: Выделить всё

>NP_001002156.1
MKTAVDRRKLDLLYSRYKDPQDENKIGVDGIQQFCDDLMLDPASVSVLIVAWKFRAATQCEFSRQEFLDG
MTDLGCDSPEKLKSLLPRLEQELKDSGKFRDFYRFTFSFAKSPGQKCLDLEMAVAYWNLILSGRFKFLGL
WNTFLLEHHKKSIPKDTWNLLLDFGNMIADDMSNYAEEGAWPVLIDDFVEFARPIVTAENLQTL
>NP_957070.2
MAKDAGLKETNGEIKLFINQSPGKAAGVLQLLTVHPASITTVKQILPKTLTVTGAHVLPHMVVSTPQRPT
IPVLLTSPHTPTAQTQQESSPWSSGHCRRADKSGKGLRHFSMKVCEKVQKKVVTSYNEVADELVQEFSSA
DHSSISPNDAVSSCHVYDQKNIRRRVYDALNVLMAMNIISKDKKEIKWIGFPTNSAQECEDLKAERQRRQ
ERIKQKQSQLQELIVQQIAFKNLVQRNREVEQQSKRSPSANTIIQLPFIIINTSKKTIIDCSISNDKFEY
LFNFDSMFEIHDDVEVLKRLGLALGLESGRCSAEQMKIATSLVSKALQPYVTEMAQGSVNQPMDFSHVAA
ERRASSSTSSRVETPTSLMEEDEEDEEEDYEEEDD
>NP_123456.1
MALLLLLGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
...

Хотя существует отличный скрипт Python для обработки поиска мотивов в мультифаст-файле (https://www.biostars.org/p/14305/), если использовался шаблон «[KHR]{3}»,
он вернет только список мотивов и множество пустых результатов:

Код: Выделить всё

>NP_001002156.1
:['RRK']
>NP_001002156.1
:[]
>NP_001002156.1
:['HHK']
>NP_957070.2
:[]
>NP_957070.2
:['RRR']
...

и в той же последовательности просочился некоторый мотив (HKK).

Здесь я нашел еще один скрипт Python:

Код: Выделить всё

#coding:utf-8
import re
pattern = "[KHR]{3}"
with open('seq.fasta') as fh:
fh.readline()
seq = ""
for line in fh:
seq += line.strip()
rgx = re.compile(pattern)
result = rgx.search(seq)
patternfound = result.group()
span = result.span()
leftpos = span[0]-10
if leftpos < 0:
leftpos = 0
print(seq[leftpos:span[0]].lower() + patternfound + seq[span[1]:span[1]+10].lower())

он возвращает первый совпадающий мотив, найденный в контексте (10 аминокислот вперед после совпадающего мотива,
и назад 10 перед совпадающим мотивом) только для одной фастовой (первой) последовательности, для первой фастовой
последовательности NP_001002156.1 с использованием скрипта, возвращаемый результат:

Код: Выделить всё

mktavdRRKldllysrykd

но у него нет заголовка файла ">NP_001002156.1" и все остальные 2 мотива в контексте были опущены:

Код: Выделить всё

glwntfllehHHKksipkdtwnl
lwntfllehhHKKsipkdtwnll

Здесь я хочу, чтобы желаемый скрипт возвращал совпадающий мотив с его позицией в контексте каждой
последовательности fasta в файле multifasta, и он представлял бы результаты следующим образом:

Код: Выделить всё

>NP_001002156.1_matchnumber_1_(7~9)
mktavdrRRKldllysrykd
>NP_001002156.1_matchnumber_2_(148~150)
glwntfllehHHKksipkdtwnl
>NP_001002156.1_matchnumber_3_(149~151)
lwntfllehhHKKsipkdtwnll
>NP_957070.2_matchnumber_1_(163~165)
chvydqknirRRRvydalnvlma
>NP_123456.1
no match found

Примечание:Положение совпадающего шаблона не является положением контекста.

Кто-нибудь может мне помочь? Заранее спасибо.

Подробнее здесь: https://stackoverflow.com/questions/541 ... ing-motifs