<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>

CC BY-SA 4.0
2021-11-19T19:59:08.183Z
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?
Amaru
CC BY-SA 4.0
2021-11-21T00:56:07.603Z
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.
conchoecia
CC BY-SA 4.0
2021-11-07T03:31:30.897Z
2021-11-07T03:31:30.897Z
conchoecia
<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>

CC BY-SA 4.0
2021-11-05T17:10:33.813Z
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?
Jesse
CC BY-SA 4.0
2021-11-08T18:16:59.027Z
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.
Amaru
<fasta><python><phylogenetics><structural-variation>
Counting hexamers in fasta sequence and identify its structure (and interruptions)
# 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)

ACATTTTTTTTCCACATCTGATGTGGAAAAAAAAAAAAATGAAATAGCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCTCCCCACGGTCTCCCTCTCATGCGGAGCCAAAGCTGGACTGTACTGCT

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
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-mers, 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.
#!/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.  

---

