<p>Read the sequences with Biopython, <a href="https://edinburgh-genome-foundry.github.io/DnaChisel/ref/biotools.html?highlight=reverse#dnachisel.biotools.reverse_translate" rel="nofollow noreferrer">reverse translate</a> the sequences and the pattern peptides, construct pattern regexes, search (<code>find_matches()</code>) for the <a href="https://edinburgh-genome-foundry.github.io/DnaChisel/ref/core_classes.html?highlight=sequencepattern#dnachisel.SequencePattern" rel="nofollow noreferrer">patterns</a> for each sequence in a loop using <a href="https://edinburgh-genome-foundry.github.io/DnaChisel/index.html" rel="nofollow noreferrer">DNA Chisel</a>.</p>
<p><em>Disclaimer: I'm the maintainer of DNA Chisel</em></p>

CC BY-SA 4.0
2022-11-16T22:49:00.737Z
2022-11-16T22:49:00.737Z
Peter
<p>First of all, you can't do this this way. Fasta sequences can, and indeed yours do, span multiple lines so in order to find a pattern, you need to join those lines into one to make sure you don't miss any that are found across more than one line.</p>
<p>Next, as a general rule, it is a bad idea to hardcode output file names. That means the user cannot control the output, and that any existing files will be overwritten, so just print to standard output and redirect to a file:</p>
<pre><code>foo.py &gt; out_file
</code></pre>
<p>You also don't want to hardcode the name of the input file since that just makes the script harder to use. Instead, read the input file name as an argument. With all this in mind, here is a working version of your script:</p>
<pre class="lang-py prettyprint-override"><code>#!/usr/bin/env python3
import re
import sys

def print_matches(seq_name, matches):
     if matches:
          for match in matches:
               print(name,match.start(), match.end(),match.group(),sep='\t' )
      

pattern = re.compile(r&quot;[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]&quot;)
seq = &quot;&quot;

try:
     with open(sys.argv[1],'r') as infile:
          for line in infile:
               line = line.strip(&quot;\n&quot;)
               if line.startswith('&gt;'):
                    name=line.lstrip('&gt;')
                    if seq != &quot;&quot;:
                         matches = re.finditer(pattern, seq)
                         print_matches(name, matches)
                         seq = &quot;&quot;
               else:
                    seq += line
          else:
               matches = re.finditer(pattern, seq)
               print_matches(name, matches)
                    
except IOError as e:
     print(&quot;Failed to open file %s with error %s&quot; % (sys.argv[1], e))
     sys.exit(1)
</code></pre>
<p>Running that on your example gives:</p>
<pre><code>$ foo.py sequence.fasta 
P2  44  53  LPNNTASWF
P2  122 131 YGANKDGII
P2  209 218 MAGNGGDAA
P2  219 228 ALLLLDRLN
P2  304 313 AQFAPSASA
P2  391 400 VTLLPAADL
P2  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
</code></pre>
<p>The output is 4 tab-separated columns per match with the sequence name, the match's start position, its end position and the actual sequence matched.</p>

CC BY-SA 4.0
2022-11-14T12:40:17.207Z
2022-11-14T12:40:17.207Z
terdon
<p>I am having trouble trying to search for motifs in a multi fasta file. I have used two techniques one gives me the name of the sequence where the motif is found but doesn't give me the the motif and its position. The other method does not return anything. Please see the code below, any assistance will be highly appreciated.</p>
<p>This is the protein sequence:
(Please note I have only put in the first 3 sequences, the actual file I am working with has 100 sequences)</p>
<pre><code>&gt;P1
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQXXTXKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

&gt;P2
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

&gt;P3
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*
</code></pre>
<p>The first one is:</p>
<pre><code>infile=open(&quot;sequence.fasta&quot;,'r')
out=open(&quot;Result.csv&quot;,'w')
pattern=re.compile(r&quot;[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]&quot;)
for line in infile:
  line = line.strip(&quot;\n&quot;)
  if line.startswith('&gt;'):
     name=line
else:
   s=re.finditer(pattern,line)
   print('%s:%s' %(name,s))
   out.write('%s:\t%s\n' %(name,s))
</code></pre>
<p>This one above returns, the output below. It gives me the name of the sequence which is P100, but it does not give me (1) The motif found (2) the position of the motif. And the Results.csv file that it generates is blank.</p>
<pre><code>#output
&gt;P100:&lt;callable_iterator object at 0x000001D0611B7AF0
</code></pre>
<p>I then tried a different technique:</p>
<pre><code>infile=open(&quot;sequence.fasta&quot;,'r')
open=open(&quot;result.csv&quot;,'w')
pattern=re.compile(r&quot;(P[A-Z]{2}P&quot;)
for line in infile:
  line=line.strip(&quot;\n&quot;)
  if line.startswith('&gt;'):
     name=line
else:
   s = re.finditer(pattern,line)
   for match_obj in s:
       print(match_obj)
       print(match_obj.group())
  out.write('%s:\t%s\t%s\n' %(name, match_obj,match_obj.group()))
</code></pre>
<p>The one above does not return anything at all. the Results.csv file is also blank.
I'm still new in this python, and some of the techniques I used here I found on bioinformatics stackexchange.</p>
<p>The output I desire, is where I can get (1) Name of sequence, (2) Motif found and (3) position of motif (where it starts and ends) as shown in the example below in the context of every sequence in the multi fasta file</p>
<pre><code>&gt;P100
motif:
position:
</code></pre>
<p>Any format is fine as long as  I can get something similar to the format above. Could I also be assistet with   the indention as that cofuses me.</p>

CC BY-SA 4.0
2022-11-14T10:13:39.147Z
@terdon I have made the amendments and put in more detail
thole
CC BY-SA 4.0
2022-11-14T10:22:35.050Z
Once again: please don't make us debug your example code. You need to make sure that the code you give us actually executes so we don't waste time on syntax errors. The code you have given does not run. for one thing, it is missing the `import re`. For another, you are still using `open` as a variable name and then trying to write to `out`. All of these would give an error if you tried to execute the code. So please make sure you give us the actual, real code that you are executing!
terdon
CC BY-SA 4.0
2022-11-20T10:18:42.913Z
Cross-post https://www.biostars.org/p/9545061/#9545129
M__
CC BY-SA 4.0
2023-10-30T00:08:47.043Z
can we use it with a file in stockholm format?
Babacar Dione
<proteins><python><motifs><multi-fasta>
motif searching in a multi-fasta file
# motif searching in a multi-fasta file

(Asked by: username_0 on 2022-11-13)

I am having trouble trying to search for motifs in a multi fasta file. I have used two techniques one gives me the name of the sequence where the motif is found but doesn't give me the the motif and its position. The other method does not return anything. Please see the code below, any assistance will be highly appreciated.
This is the protein sequence:
(Please note I have only put in the first 3 sequences, the actual file I am working with has 100 sequences)
>P1
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQXXTXKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P2
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P3
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

The first one is:
infile=open("sequence.fasta",'r')
out=open("Result.csv",'w')
pattern=re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s=re.finditer(pattern,line)
   print('%s:%s' %(name,s))
   out.write('%s:\t%s\n' %(name,s))

This one above returns, the output below. It gives me the name of the sequence which is P100, but it does not give me (1) The motif found (2) the position of the motif. And the Results.csv file that it generates is blank.
#output
>P100:
I then tried a different technique:
infile=open("sequence.fasta",'r')
open=open("result.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.finditer(pattern,line)
   for match_obj in s:
       print(match_obj)
       print(match_obj.group())
  out.write('%s:\t%s\t%s\n' %(name, match_obj,match_obj.group()))

The one above does not return anything at all. the Results.csv file is also blank.
I'm still new in this python, and some of the techniques I used here I found on bioinformatics stackexchange.
The output I desire, is where I can get (1) Name of sequence, (2) Motif found and (3) position of motif (where it starts and ends) as shown in the example below in the context of every sequence in the multi fasta file
>P100
motif:
position:

Any format is fine as long as  I can get something similar to the format above. Could I also be assistet with   the indention as that cofuses me.

**Tags:** proteins, python, motifs, multi-fasta
> **Comments:**
> **username_0**: @terdon I have made the amendments and put in more detail  
> **username_1**: Once again: please don't make us debug your example code. You need to make sure that the code you give us actually executes so we don't waste time on syntax errors. The code you have given does not run. for one thing, it is missing the `import re`. For another, you are still using `open` as a variable name and then trying to write to `out`. All of these would give an error if you tried to execute the code. So please make sure you give us the actual, real code that you are executing!  
> **username_2**: Cross-post https://www.biostars.org/p/9545061/#9545129  
> **username_3**: can we use it with a file in stockholm format?  

---

## 2 answers

### Answer by username_1 on 2022-11-14. Score: 3 | Accepted answer

First of all, you can't do this this way. Fasta sequences can, and indeed yours do, span multiple lines so in order to find a pattern, you need to join those lines into one to make sure you don't miss any that are found across more than one line.
Next, as a general rule, it is a bad idea to hardcode output file names. That means the user cannot control the output, and that any existing files will be overwritten, so just print to standard output and redirect to a file:
foo.py > out_file

You also don't want to hardcode the name of the input file since that just makes the script harder to use. Instead, read the input file name as an argument. With all this in mind, here is a working version of your script:
#!/usr/bin/env python3
import re
import sys

def print_matches(seq_name, matches):
     if matches:
          for match in matches:
               print(name,match.start(), match.end(),match.group(),sep='\t' )
      

pattern = re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
seq = ""

try:
     with open(sys.argv[1],'r') as infile:
          for line in infile:
               line = line.strip("\n")
               if line.startswith('>'):
                    name=line.lstrip('>')
                    if seq != "":
                         matches = re.finditer(pattern, seq)
                         print_matches(name, matches)
                         seq = ""
               else:
                    seq += line
          else:
               matches = re.finditer(pattern, seq)
               print_matches(name, matches)
                    
except IOError as e:
     print("Failed to open file %s with error %s" % (sys.argv[1], e))
     sys.exit(1)

Running that on your example gives:
$ foo.py sequence.fasta 
P2  44  53  LPNNTASWF
P2  122 131 YGANKDGII
P2  209 218 MAGNGGDAA
P2  219 228 ALLLLDRLN
P2  304 313 AQFAPSASA
P2  391 400 VTLLPAADL
P2  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA

The output is 4 tab-separated columns per match with the sequence name, the match's start position, its end position and the actual sequence matched.


---

### Answer by username_4 on 2022-11-16. Score: 1

Read the sequences with Biopython, [reverse translate](https://edinburgh-genome-foundry.github.io/DnaChisel/ref/biotools.html?highlight=reverse#dnachisel.biotools.reverse_translate) the sequences and the pattern peptides, construct pattern regexes, search (`find_matches()`) for the [patterns](https://edinburgh-genome-foundry.github.io/DnaChisel/ref/core_classes.html?highlight=sequencepattern#dnachisel.SequencePattern) for each sequence in a loop using [DNA Chisel](https://edinburgh-genome-foundry.github.io/DnaChisel/index.html).

*Disclaimer: I'm the maintainer of DNA Chisel*


---

