Answers.Body: <p>I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):</p>
  
  <pre><code>#!/usr/bin/env python3
  
  import sys
  
  sequence = sys.stdin.readlines()
  d = {}
  temp_genename = None
  temp_sequence = None
  # The lines returned by readlines() contain end-of-line characters,
  # this is why strip() is needed
  motif = sequence[0].strip()
  
  def motif_count(motif, temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
          if start &lt;0:
              return count
          count += 1
  
  if sequence[1][0] != '&gt;':
  
      print("ERROR")
  
      exit(1)
  
  for line in sequence[1:]:
      if line[0] == '&gt;':
          # Record the previous one in the dictionary
          if temp_genename is not None:
              d[temp_genename] = temp_sequence
          temp_genename = line.split(' ')[0].strip()
          temp_genename = temp_genename.replace('&gt;', '')
          temp_sequence = ""
  
      else:
  
          temp_sequence += line.strip()
  
  # This should be done each time a new sequence is found,
  # not only at the end
  d[temp_genename] = temp_sequence
  for key in d:
      # We need to use the result
      # motif_count(motif, key)
      # We want to count the motif in the sequence, not in its name
      count = motif_count(motif, d[key])
      print(f"{key}: {count}")
  </code></pre>
  
  <p>Testing it:</p>
  
  <pre><code><span class="math-container">$ cat tmp.fa 
  AT
  &gt;seq1 Comment......
  AGGTATA
  TGGCGCC
  &gt;seq2 Comment.....
  GGCCGGCGC
  GCCGGCGGG
  $</span> cat tmp.fa | ./count_motif.py
  seq1: 2
  seq2: 0
  </code></pre>
  
  <p>Here is another way to do, using <code>valmap</code> to make a dict of the motif counts from the dict of sequences, and with a few other programming style changes:</p>
  
  <pre><code>#!/usr/bin/env python3
  
  import sys
  from cytoolz import valmap
  
  motif = sys.stdin.readline().strip()
  
  
  # We can make this a function of only one parameter
  # because motif is fixed
  def motif_count(temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
          if start &lt;0:
              return count
          count += 1
  
  
  seq_dict = {}
  temp_genename = None
  temp_sequence = None
  
  for line in sys.stdin:
      if line[0] == '&gt;':
          # Record the previous one in the dictionary
          if temp_genename is not None:
              seq_dict[temp_genename] = temp_sequence
          temp_genename, *_ = line[1:].strip().split()
          temp_sequence = ""
      else:
          temp_sequence += line.strip()
  
  msg = f"""File format might be wrong"
  motif: {motif}
  temp_genename: {temp_genename}
  temp_genename: {temp_sequence}"""
  
  assert temp_sequence is not None and len(temp_sequence), msg
  
  # Record the last sequence
  seq_dict[temp_genename] = temp_sequence
  
  # Make a new dict where the values are the motif counts
  # valmap needs the function motif_count to take only one parameter
  counts_dict = valmap(motif_count, seq_dict)
  # Print the results, iterating over dict items
  for (name, count) in counts_dict.items():
      print(f"{name}: {count}")
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-05-09T21:18:06.467Z
Answers.Text: The lack of "exit (1)" indentation was an important bug
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-05-09T14:35:39.677Z
Answers.LastActivityDate: 2019-05-09T14:52:03.103Z
Answers.LastEditDate: 2019-05-09T14:52:03.103Z
Answers.UserDisplayName: bli
Answers.Body: <p>This is a good task for <a href="https://biopython.org/" rel="nofollow noreferrer">Biopython</a>:</p>
  
  <pre><code>from Bio import SeqIO
  
  with open('input_file') as f:
      motif = next(f).rstrip().lower()
      for r in SeqIO.parse(f, 'fasta'):
          print(r.id, ':', r.seq.lower().count_overlap(motif))
  </code></pre>
  
  <p>Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-05-08T20:22:41.167Z
Answers.Text: The OP seems to want to read the input from stdin. Using a file is a far superior approach, of course, but maybe stdin is a requirement.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-05-08T19:46:08.877Z
Answers.LastActivityDate: 2019-05-08T19:46:08.877Z
Answers.UserDisplayName: Chris_Rands
Answers.Body: <p>I like @Chris_Rands tidy use of OO, objects are really important per se.</p>
  
  <p>There appears a paradox with your assignment that you might want to refer back to your supervisor/tutors.</p>
  
  <p>If "AT" is the motif, this cannot overlap with another "AT" motif, therefore the request for "overlapping motifs" makes this part of the code superfluous.</p>
  
  <p>It would be better expressed if the motif was "ATA" for example. Thus if the sequence was ATATATA the motif is present 3 times, but only twice if the motif was contiguous. </p>
  
  <p>More generally the code essentially searches for AT microsatellite motifs and formally performed using e.g. tandem repeat finder.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-05-09T04:09:23.063Z
Answers.LastActivityDate: 2019-05-09T04:15:24.647Z
Answers.LastEditDate: 2019-05-09T04:15:24.647Z
Answers.UserDisplayName: M__
Body: <p>I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with >, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:</p>
  
  <pre><code> AT
   &gt;seq1 Comment......
   AGGTATA
   TGGCGCC
   &gt;seq2 Comment.....
   GGCCGGCGC
   GCCGGCGGG
  </code></pre>
  
  <p>The output should be:</p>
  
  <pre><code>seq1: 2
  seq2: 0
  </code></pre>
  
  <p>I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif_count does not produce any result when called on my d values.</p>
  
  <p>Here is my code:</p>
  
  <pre><code> #!/usr/bin/env python3
  
   import sys
  
   sequence = sys.stdin.readlines()
   d = {}
   temp_genename = None
   temp_sequence = None
   motif = sequence[0]
  
  def motif_count(motif, temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
      if start &lt;0:
          return count
      count += 1
  
  if sequence[1][0] != '&gt;':
  
      print("ERROR")
  
  exit(1)
  
  for line in sequence[1:]:
      if line[0] == '&gt;':
         temp_genename = line.split(' ')[0].strip()
         tmp_genename = temp_genename.replace('&gt;', '')
         temp_sequence = ""
  
  else:
  
        temp_sequence += line.strip()
  
  d[tmp_genename] = temp_sequence
  for key in d:
      motif_count(motif, key)
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-05-08T19:06:59.737Z
Comments.Text: Hello Dari, this looks a bit like a homework. It's ok to ask, but it's nice to disclaim it.
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-05-09T13:10:12.900Z
Comments.Text: Maybe it's just an error when copying your actual code in your question, but in `motif_count`, the `while` loop does not contain the test that can trigger the `return`, nor the incrementation of `count`. `exit(1)` should also probably be re-indented. I tried to fix this. Hopefully this corresponds to your actual code now.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-05-09T13:16:51.600Z
Comments.Text: Your code does nothing with the values returned by `motif_count`. To create a new dictionary with values modified by a function, the `valmap` function from the `cytoolz` module could be useful.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-05-09T13:19:28.060Z
Comments.Text: @bli please don't correct code in questions! We can't know if an error is because the OP didn't copy/paste correctly or if that's the actual code being used, so when you fix errors in the question, it could end up causing everyone to waste their time trying to find the error in code that is now correct.
Comments.UserDisplayName: terdon
Tags: <python><phylogenetics><homework>
Title: Calculate the occurrence of motif in sequences (including overlaps) using Python
ThreadText: # Calculate the occurrence of motif in sequences (including overlaps) using Python
  
  (Asked by: username_0 on 2019-05-08)
  
  I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with >, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:
  
  ```
   AT
   >seq1 Comment......
   AGGTATA
   TGGCGCC
   >seq2 Comment.....
   GGCCGGCGC
   GCCGGCGGG
  
  ```
  
  The output should be:
  
  ```
  seq1: 2
  seq2: 0
  
  ```
  
  I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif\_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif\_count does not produce any result when called on my d values.
  
  Here is my code:
  
  ```
   #!/usr/bin/env python3
  
   import sys
  
   sequence = sys.stdin.readlines()
   d = {}
   temp_genename = None
   temp_sequence = None
   motif = sequence[0]
  
  def motif_count(motif, temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
      if start <0:
          return count
      count += 1
  
  if sequence[1][0] != '>':
  
      print("ERROR")
  
  exit(1)
  
  for line in sequence[1:]:
      if line[0] == '>':
         temp_genename = line.split(' ')[0].strip()
         tmp_genename = temp_genename.replace('>', '')
         temp_sequence = ""
  
  else:
  
        temp_sequence += line.strip()
  
  d[tmp_genename] = temp_sequence
  for key in d:
      motif_count(motif, key)
  
  ```
  
  **Tags:** python, phylogenetics, homework
  
  > **Comments:**
  > **username_1**: Hello username_0, this looks a bit like a homework. It's ok to ask, but it's nice to disclaim it.  
  > **username_2**: Maybe it's just an error when copying your actual code in your question, but in `motif_count`, the `while` loop does not contain the test that can trigger the `return`, nor the incrementation of `count`. `exit(1)` should also probably be re-indented. I tried to fix this. Hopefully this corresponds to your actual code now.  
  > **username_2**: Your code does nothing with the values returned by `motif_count`. To create a new dictionary with values modified by a function, the `valmap` function from the `cytoolz` module could be useful.  
  > **username_3**: @username_2 please don't correct code in questions! We can't know if an error is because the OP didn't copy/paste correctly or if that's the actual code being used, so when you fix errors in the question, it could end up causing everyone to waste their time trying to find the error in code that is now correct.  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_2 on 2019-05-09. Score: 1 | Accepted answer
  
  I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):
  
  ```
  #!/usr/bin/env python3
  
  import sys
  
  sequence = sys.stdin.readlines()
  d = {}
  temp_genename = None
  temp_sequence = None
  # The lines returned by readlines() contain end-of-line characters,
  # this is why strip() is needed
  motif = sequence[0].strip()
  
  def motif_count(motif, temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
          if start <0:
              return count
          count += 1
  
  if sequence[1][0] != '>':
  
      print("ERROR")
  
      exit(1)
  
  for line in sequence[1:]:
      if line[0] == '>':
          # Record the previous one in the dictionary
          if temp_genename is not None:
              d[temp_genename] = temp_sequence
          temp_genename = line.split(' ')[0].strip()
          temp_genename = temp_genename.replace('>', '')
          temp_sequence = ""
  
      else:
  
          temp_sequence += line.strip()
  
  # This should be done each time a new sequence is found,
  # not only at the end
  d[temp_genename] = temp_sequence
  for key in d:
      # We need to use the result
      # motif_count(motif, key)
      # We want to count the motif in the sequence, not in its name
      count = motif_count(motif, d[key])
      print(f"{key}: {count}")
  
  ```
  
  Testing it:
  
  ```
  $ cat tmp.fa 
  AT
  >seq1 Comment......
  AGGTATA
  TGGCGCC
  >seq2 Comment.....
  GGCCGGCGC
  GCCGGCGGG
  $ cat tmp.fa | ./count_motif.py
  seq1: 2
  seq2: 0
  
  ```
  
  Here is another way to do, using `valmap` to make a dict of the motif counts from the dict of sequences, and with a few other programming style changes:
  
  ```
  #!/usr/bin/env python3
  
  import sys
  from cytoolz import valmap
  
  motif = sys.stdin.readline().strip()
  
  # We can make this a function of only one parameter
  # because motif is fixed
  def motif_count(temp_sequence):
      count = 0
      start = -1
      while True:
          start = temp_sequence.find(motif, start +1)
          if start <0:
              return count
          count += 1
  
  seq_dict = {}
  temp_genename = None
  temp_sequence = None
  
  for line in sys.stdin:
      if line[0] == '>':
          # Record the previous one in the dictionary
          if temp_genename is not None:
              seq_dict[temp_genename] = temp_sequence
          temp_genename, *_ = line[1:].strip().split()
          temp_sequence = ""
      else:
          temp_sequence += line.strip()
  
  msg = f"""File format might be wrong"
  motif: {motif}
  temp_genename: {temp_genename}
  temp_genename: {temp_sequence}"""
  
  assert temp_sequence is not None and len(temp_sequence), msg
  
  # Record the last sequence
  seq_dict[temp_genename] = temp_sequence
  
  # Make a new dict where the values are the motif counts
  # valmap needs the function motif_count to take only one parameter
  counts_dict = valmap(motif_count, seq_dict)
  # Print the results, iterating over dict items
  for (name, count) in counts_dict.items():
      print(f"{name}: {count}")
  
  ```
  
  > **Comments:**
  > **username_4**: The lack of "exit (1)" indentation was an important bug  
  
  ---
  
  ### Answer by username_5 on 2019-05-08. Score: 3
  
  This is a good task for [Biopython](https://biopython.org/):
  
  ```
  from Bio import SeqIO
  
  with open('input_file') as f:
      motif = next(f).rstrip().lower()
      for r in SeqIO.parse(f, 'fasta'):
          print(r.id, ':', r.seq.lower().count_overlap(motif))
  
  ```
  
  Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.
  
  > **Comments:**
  > **username_3**: The OP seems to want to read the input from stdin. Using a file is a far superior approach, of course, but maybe stdin is a requirement.  
  
  ---
  
  ### Answer by username_4 on 2019-05-09. Score: 0
  
  I like @Chris\_Rands tidy use of OO, objects are really important per se.
  
  There appears a paradox with your assignment that you might want to refer back to your supervisor/tutors.
  
  If "AT" is the motif, this cannot overlap with another "AT" motif, therefore the request for "overlapping motifs" makes this part of the code superfluous.
  
  It would be better expressed if the motif was "ATA" for example. Thus if the sequence was ATATATA the motif is present 3 times, but only twice if the motif was contiguous.
  
  More generally the code essentially searches for AT microsatellite motifs and formally performed using e.g. tandem repeat finder.
  
  
  ---
  
  