Answers.Body: <p>I was going to write some clever code counting <em>k-mer</em>s, but this problem can easily be solved using the handy sequentially-operating feature of python's str.replace() method.</p>
  <p>This code doesn't work if you need to look at the motif occurring on the other strand. Just the forward strand.</p>
  <pre class="lang-py prettyprint-override"><code>#!/usr/bin/env python3
  
  mystring = &quot;CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCC\
  TCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT\
  CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT&quot;
  
  motif = &quot;CCCTCT&quot;
  
  # This replaces the motif with 1s
  newstring = mystring.replace(motif, &quot;1&quot;)
  # results in 1111111111111111111111111111111111111CCT11
  # silly, right?
  
  # This whole block just parses that string and prints everything
  #  the way that OP wanted.
  done = False
  space=&quot;&quot;
  accumulator = 0
  while not done:
      if len(newstring) == 0:
          done = True
          break
      else:
          if newstring[0] == &quot;1&quot;:
              accumulator += 1
          else:
              if accumulator != 0:
                  print(&quot;{}({}){} &quot;.format(space, motif, accumulator), end = &quot;&quot;)
                  accumulator = 0
              print(newstring[0], end = &quot;&quot;)
              space = &quot; &quot;
      newstring = newstring[1::]
  
  # This is our check after exiting the while loop.
  #  This only fires if we stepped off cleanly at the
  #  end of a motif.
  if accumulator &gt; 0:
      print(&quot;{}({}){} &quot;.format(space, motif, accumulator), end = &quot;&quot;)
  print()
  </code></pre>
  <p>Prints, <code>(CCCTCT)37 CCT (CCCTCT)2</code>, matching OP's spec.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-19T19:59:08.183Z
Answers.Text: That looks good, but what if I have thousand of reads so I can use fasta/bam files as input. This looks fine when I have a read, am I right?
Answers.UserDisplayName: Amaru
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-21T00:56:07.603Z
Answers.Text: Not sure what you mean by "this looks fine when I have a read", but this code can be put into a method and worked into your framework you already have.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-11-07T03:31:30.897Z
Answers.LastActivityDate: 2021-11-07T03:31:30.897Z
Answers.UserDisplayName: conchoecia
Body: <p>I have a lot of fasta files, each one with thousand of reads containing the hexameric motif &quot;CCCTCT&quot;. The hexameric motif is highly continuous in most cases but interruptions may occur. I need to count the hexameric motif keeping the read ID and identify their structure (CCCTCT)n. These interruptions are real and do not come from Sequencing error</p>
  <p>Example Sequence:(interruption in bold)</p>
  <blockquote>
  <p>ACATTTTTTTTCCACATCTGATGTGGAAAAAAAAAAAAATGAAATAGCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT<strong>CCT</strong>CCCTCTCCCTCTCCCCACGGTCTCCCTCTCATGCGGAGCCAAAGCTGGACTGTACTGCT</p>
  </blockquote>
  <p>Strategy 1: Count the CCCTCT motif</p>
  <pre><code>grep -on &quot;CCCTCT&quot; fasta.fa | cut -f 1 -d &quot;:&quot; | sort | uniq -c
  </code></pre>
  <blockquote>
  <p>Output: <code>23 CCCTCT</code></p>
  </blockquote>
  <p>As you can see we have one interruption (<code>CCT</code>) in this sequence which the codes doesn't detect.</p>
  <p>Strategy 2: Identify reads with the same structure and count repeats in each read</p>
  <pre class="lang-py prettyprint-override"><code>from Bio import SeqIO
  from collections import defaultdict
  
  dedup_records = defaultdict(list)
  for record in SeqIO.parse(&quot;fasta.fa&quot;, &quot;fasta&quot;):
      # Use the sequence as the key and then have a list of id's as the value
      dedup_records[str(record.seq)].append(record.id)
  with open(&quot;Output.fasta&quot;, 'w') as output:
      for seq, ids in sorted(dedup_records.items(), key=lambda t: len(t[1]), reverse=True):
          # Join the ids and write them out as the fasta
          output.write(&quot;&gt;{}_counts{}\n&quot;.format(ids[0], len(ids)))
          output.write(seq + &quot;\n&quot;)
  </code></pre>
  <p>Output: list ranking of reads with identical structures but does not show the structure</p>
  <pre><code>&gt;ID/ccs_counts2
  CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT.....
  </code></pre>
  <p>Can you please give me an idea of how to obtain the following output? example of a sequence of 23 hexamers with interruptions.</p>
  <pre><code>(CCCTCT)20 [INTERRUPTION] (CCCTCT)3
  </code></pre>
  <p>I have used expansion hunter, repeatFinder and repeatAnalysis-tool by PacBio but they were not particularly sensitive for this case. Any help would be appreciated.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-05T17:10:33.813Z
Comments.Text: I've also mulled on this sort of puzzle in the past (when dealing with microsatellites) and I feel like there should be some elegant approaches that could fit.  How strict do you want the definition of a repeat region and interruption between repeats?  For example, more than one perfect tandem repeat triggers a match, and less than one repeat-length mismatch before it resumes counts as an interruption (rather than the end of the whole repeat)?  Or even simpler, just however many perfect matches with any-length interruptions between?
Comments.UserDisplayName: Jesse
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-11-08T18:16:59.027Z
Comments.Text: Thanks for your thoughts. I have to be either flexible to detect some changes (if there is some) and strict to be able to report those sequences following the patterns. I guess I just figured it out.
Comments.UserDisplayName: Amaru
Tags: <fasta><python><phylogenetics><structural-variation>
Title: Counting hexamers in fasta sequence and identify its structure (and interruptions)
ThreadText: # Counting hexamers in fasta sequence and identify its structure (and interruptions)
  
  (Asked by: username_0 on 2021-11-02)
  
  I have a lot of fasta files, each one with thousand of reads containing the hexameric motif "CCCTCT". The hexameric motif is highly continuous in most cases but interruptions may occur. I need to count the hexameric motif keeping the read ID and identify their structure (CCCTCT)n. These interruptions are real and do not come from Sequencing error
  
  Example Sequence:(interruption in bold)
  
  > ACATTTTTTTTCCACATCTGATGTGGAAAAAAAAAAAAATGAAATAGCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT**CCT**CCCTCTCCCTCTCCCCACGGTCTCCCTCTCATGCGGAGCCAAAGCTGGACTGTACTGCT
  
  Strategy 1: Count the CCCTCT motif
  
  ```
  grep -on "CCCTCT" fasta.fa | cut -f 1 -d ":" | sort | uniq -c
  
  ```
  
  > Output: `23 CCCTCT`
  
  As you can see we have one interruption (`CCT`) in this sequence which the codes doesn't detect.
  
  Strategy 2: Identify reads with the same structure and count repeats in each read
  
  ```py
  from Bio import SeqIO
  from collections import defaultdict
  
  dedup_records = defaultdict(list)
  for record in SeqIO.parse("fasta.fa", "fasta"):
      # Use the sequence as the key and then have a list of id's as the value
      dedup_records[str(record.seq)].append(record.id)
  with open("Output.fasta", 'w') as output:
      for seq, ids in sorted(dedup_records.items(), key=lambda t: len(t[1]), reverse=True):
          # Join the ids and write them out as the fasta
          output.write(">{}_counts{}\n".format(ids[0], len(ids)))
          output.write(seq + "\n")
  
  ```
  
  Output: list ranking of reads with identical structures but does not show the structure
  
  ```
  >ID/ccs_counts2
  CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT.....
  
  ```
  
  Can you please give me an idea of how to obtain the following output? example of a sequence of 23 hexamers with interruptions.
  
  ```
  (CCCTCT)20 [INTERRUPTION] (CCCTCT)3
  
  ```
  
  I have used expansion hunter, repeatFinder and repeatAnalysis-tool by PacBio but they were not particularly sensitive for this case. Any help would be appreciated.
  
  **Tags:** fasta, python, phylogenetics, structural-variation
  
  > **Comments:**
  > **username_1**: I've also mulled on this sort of puzzle in the past (when dealing with microsatellites) and I feel like there should be some elegant approaches that could fit.  How strict do you want the definition of a repeat region and interruption between repeats?  For example, more than one perfect tandem repeat triggers a match, and less than one repeat-length mismatch before it resumes counts as an interruption (rather than the end of the whole repeat)?  Or even simpler, just however many perfect matches with any-length interruptions between?  
  > **username_0**: Thanks for your thoughts. I have to be either flexible to detect some changes (if there is some) and strict to be able to report those sequences following the patterns. I guess I just figured it out.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2021-11-07. Score: 1
  
  I was going to write some clever code counting *k-mer*s, but this problem can easily be solved using the handy sequentially-operating feature of python's str.replace() method.
  
  This code doesn't work if you need to look at the motif occurring on the other strand. Just the forward strand.
  
  ```py
  #!/usr/bin/env python3
  
  mystring = "CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCC\
  TCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT\
  CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT"
  
  motif = "CCCTCT"
  
  # This replaces the motif with 1s
  newstring = mystring.replace(motif, "1")
  # results in 1111111111111111111111111111111111111CCT11
  # silly, right?
  
  # This whole block just parses that string and prints everything
  #  the way that OP wanted.
  done = False
  space=""
  accumulator = 0
  while not done:
      if len(newstring) == 0:
          done = True
          break
      else:
          if newstring[0] == "1":
              accumulator += 1
          else:
              if accumulator != 0:
                  print("{}({}){} ".format(space, motif, accumulator), end = "")
                  accumulator = 0
              print(newstring[0], end = "")
              space = " "
      newstring = newstring[1::]
  
  # This is our check after exiting the while loop.
  #  This only fires if we stepped off cleanly at the
  #  end of a motif.
  if accumulator > 0:
      print("{}({}){} ".format(space, motif, accumulator), end = "")
  print()
  
  ```
  
  Prints, `(CCCTCT)37 CCT (CCCTCT)2`, matching OP's spec.
  
  > **Comments:**
  > **username_0**: That looks good, but what if I have thousand of reads so I can use fasta/bam files as input. This looks fine when I have a read, am I right?  
  > **username_2**: Not sure what you mean by "this looks fine when I have a read", but this code can be put into a method and worked into your framework you already have.  
  
  ---
  
  