2
$\begingroup$

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

$\endgroup$

1 Answer 1

0
$\begingroup$

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')
$\endgroup$
1
  • $\begingroup$ Thank you so much for your suggestion, I am still new python user and will read up on biopython. With regards to the code will it now be: else:s=re.finditer (pattern,line) print'%s:%' %(name,s) ? Also by parse the finditer object what does that mean? I am struggling to link it to the original code. Could you perhaps try to edit your answer such that it looks like the original code, so that I can see where the changes will need to be made. $\endgroup$
    – thole
    Nov 13, 2022 at 6:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.