finding motifs in fasta file with multiple sequences

Question

I am able to find matching motifs on a single sequence, as well as their position as shown below :

demo = "ABCLXLYLZLABC"
motif="L.L.L"
for element in re.finditer(motif,demo):
    print(element.start(),element.end,element.group())

This works and gives me the desired output which is: 3 8 LXLYL

Now I want to know how can I then do this same process but on a multifasta file. As done in this example script from https://www.biostars.org/p/14305/ .

The protein sequences on the file are :

> P1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL
>P2
 MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
 LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
 DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
 ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
 ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
 ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
 EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
 PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR
>P3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

The script that was used to find the motif is below, the problem now is that it does not show the position where the motif is found. Could you assist in how one could do so.

import re
infile=open("seq.fasta",'r')
out=open("results.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P)")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:      
   s = re.findall(pattern,line)
   print '%s:%s' %(name,s)
   out.write('%s:\t%s\n' %(name,s))

Their script seems to have worked because the desired output was

>p1: 
PGCP

>p1: 
PHCP,
PKCP

The problem with their script is that, it does not provide the position where the sequence was found. Could you please assist on how that script can be modified? I want to try it but I want my output to have (1) Name Of Sequence, (2) The Motif Found , (3) The Position

Answer 1

This is easy just replace findall with finditer in the above code. The index will then appear where the match first appears and where it stops. However, you will to parse the finditer object which in the code is s

Then replace

print '%s:%s' %(name,s)

with

for match_obj in s:
    # print each re.Match object
    print(match_obj)
    
    # extract each matching number
    print(match_obj.group())

That should work. When it does you simply pass the match_obj and match_obj.group() to a file by replacing s with match_obj, match_obj.group()

Thus replace

  out.write('%s:\t%s\n' %(name,s))

with

  out.write('%s:\t%s\t%s\n' %(name,match_obj,match_obj.group()))

You might need to do an append BTW.

thus replace

out=open("results.csv",'w')

with

import os
if os.path.exists("results.csv"):
    os.unlink ("results.csv")
out=open("results.csv",'a')