Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-11-16T22:49:00.737Z
Answers.LastActivityDate: 2022-11-16T22:49:00.737Z
Answers.UserDisplayName: Peter
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-11-14T12:40:17.207Z
Answers.LastActivityDate: 2022-11-14T12:40:17.207Z
Answers.UserDisplayName: terdon
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-14T10:13:39.147Z
Comments.Text: @terdon I have made the amendments and put in more detail
Comments.UserDisplayName: thole
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-14T10:22:35.050Z
Comments.Text: 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!
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-11-20T10:18:42.913Z
Comments.Text: Cross-post https://www.biostars.org/p/9545061/#9545129
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-30T00:08:47.043Z
Comments.Text: can we use it with a file in stockholm format?
Comments.UserDisplayName: Babacar Dione
Tags: <proteins><python><motifs><multi-fasta>
Title: motif searching in a multi-fasta file
ThreadText: # 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:<callable_iterator object at 0x000001D0611B7AF0
  
  ```
  
  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**: @username_1 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:
  
  ```py
  #!/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*
  
  
  ---
  
  