<p>You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As.  You might also need to scale by the sequence length, if you have different length sequences that you want to compare.</p>

<p>Edit: also, there might be a correction to make if one or both ends of the sequence are an A</p>

CC BY-SA 4.0
2019-02-02T02:35:42.147Z
Interestingly idea. Can you think of any edge cases that would cause problems for this scheme?
Tom
CC BY-SA 4.0
2019-02-04T13:51:54.357Z
Zero As would throw a spanner in the works.
Jay Moore
CC BY-SA 4.0
2019-02-01T10:39:11.890Z
2019-02-01T10:49:40.930Z
2019-02-01T10:49:40.930Z
Jay Moore
<p>The Markov chain approach is very interesting. You could deploy a maximum likelihood approach where instead of scoring the genetic distance within a homologue, the algorithm looks at its physically linked neighbour. The genetic distance of the transition point i.e. AAAA then CCCC is deducted from the final genetic distance so that AAAAACCCCCC would be zero. It would be fairly easy to script and use PAML or equivalent where the default phylogeny would be a star tree.  </p>

<p>The challenge is to establish a Deep Learning approach. Its not impossible, you just assign pre-determined classes of oligos which fall into the definition you want and use these as your training set.</p>

CC BY-SA 4.0
2019-02-03T08:19:00.933Z
2019-02-03T08:19:00.933Z
M__
<p>To define the base density score, or clustering score, of a sequence <span class="math-container">$S$</span> we can think of each position in the sequence as being part of a <a href="https://en.wikipedia.org/wiki/Markov_chain" rel="nofollow noreferrer">Markov process</a>. For the sake of symmetry around a single base we consider odd-order Markov models from 3-onward where we ask what is the probability that we see the character at <span class="math-container">$S_i$</span> in the <span class="math-container">$j$</span> bases from <span class="math-container">$[S_{i-j}, S_{i+j}]$</span>. For each column we have <code>length(S) -1</code>, <span class="math-container">$|S|-1$</span>, possible probabilities to calculate per window. </p>

<p>We can't look left or right of the edges of the sequences so the probabilities are skewed, but we can get a rough idea of what is going on.</p>

<ul>
<li>For each base position we calculate a clustering score of <span class="math-container">$ \frac{1}{w}\Sigma_{j=1}^{j=w} -log_2 \frac{m(S_{i-j}, S_{i+j})-1 + \frac{1}{|S|}}{|S_{\{i-j, i+j\}}| -1 }$</span> where:

<ul>
<li><span class="math-container">$w$</span> is the window size we can look at. The max is <span class="math-container">$|S|-1$</span> and the min is <span class="math-container">$1$</span>. This just means that we can calculate a score in the context of the whole sequence or we maybe only care within some smaller window, like 5 bp flanking either side of the base in question.</li>
<li>The function <span class="math-container">$m(S_{i-j}, S_{i+j})$</span> is just the number of counts of matching bases in the window <span class="math-container">$j$</span> on either side of the base <span class="math-container">$S_i$</span>. For example in the sequence <code>gatAgga</code>, if <span class="math-container">$S_i$</span> is the uppercase <code>A</code> and <span class="math-container">$j$</span> is 3, the function <span class="math-container">$m(S_{i-3}, S_{i+3})$</span> returns 3, since there is the <code>A</code> base at <span class="math-container">$S_i$</span> and the two additional <code>a</code> bases at <span class="math-container">$i-2$</span> and at <span class="math-container">$i+3$</span>. We have to subtract one from this number because we don't want to factor in the matching base at <span class="math-container">$S_i$</span>.</li>
<li>We then add <span class="math-container">$\frac{1}{|S|}$</span> to that number. This is a pseudocount that we need to add to prevent from taking the log of zero. This number is pretty arbitrary but I chose <span class="math-container">$\frac{1}{|S|}$</span> so that non-clusting regions' scores aren't too high.</li>
<li>We then divide that numerator by the length of <span class="math-container">$S$</span> from the little window that we just looked at minus one, <span class="math-container">$|S_{\{i-j, i+j\}}| -1 $</span>. For example in the sequence <code>gatAgga</code> from above we only want to consider the bases <code>gat gga</code>. This is 6 bases.</li>
<li>We then take the negative log to convert the proportional scale to a log scale.</li>
<li>The summation just says that we do this for every symmetrical window size around the base in question then sum them up. For example we would calculate the ratio for <code>tAg</code>, <code>atAgg</code>, and <code>gatAgga</code>.</li>
<li>The <span class="math-container">$\frac{1}{w}$</span> just gives the average log clustering score of all the windows that we looked at. One nice thing about this is that the impact of each increasing window size has less effect on the ratio, so there is a sort of convergence that happens.</li>
</ul></li>
</ul>

<p>Implemented in python:</p>

<pre class="lang-python prettyprint-override"><code>from math import log as mlog
s = "GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA"
assert len(s) &gt; 0
pseudocount = 1/len(s)
#maxdepth = 5
maxdepth = len(s)
score = []
for i in range(len(s)): # iterate positions
    temp = []
    for j in range(1, maxdepth): # iterate depth
        rmin = max(i-j, 0)
        rmax = i+j+1
        num_matches_min_one = sum([1 for k in s[rmin:rmax] if k == s[i]])-1
        len_s_min_one = len(s[rmin:rmax]) -1
        sim = abs(-1 * mlog((num_matches_min_one/len_s_min_one) + pseudocount, 2))
        #print("i={}, j={}, rmin={}, rmax={}, sim={}".format(i, j, rmin, rmax, sim))
        temp.append(sim)
    print("{0:d}\t{1}\t".format(i, s[i]), end='\t')
    score.append(sum(temp)/len(temp))
    print("{0:.4f}".format(score[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x = [i+1 for i in range(len(s))]
plt.bar(x, score)
plt.xticks(x, [k for k in s])
ax.set_yscale("log", nonposy='clip')
plt.xlabel('sequence')
plt.ylabel('-log2(cluster-ness)')
plt.title("cluster-ness score")
plt.savefig('output.png', dpi=150)
</code></pre>

<p>This calculates the matrix, prints out a tab-delimited file of the position, base, and log clusting score, then saves a plot of the results.</p>

<p>For the test sequence <code>GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA</code> and using the whole sequence as context, I got the following plot:</p>

<p><a href="https://i.sstatic.net/C3SdP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/C3SdP.png" alt="enter image description here"></a></p>

<p>As you can see, the log clustering score is low around the things that are cluster-y and high around things that aren't cluster-y.</p>

<p><strong>Note:</strong> OP requested that the score be scaled from 0 to 1. I instead used the log scale to make visualization easier.</p>

<p><strong>Warning 1:</strong> If you plan on doing this for large sequences, it is really important to limit the <code>maxdepth</code>, or <span class="math-container">$w$</span>, parameter. The time complexity for this algorithm is <span class="math-container">$O(|S|*w)$</span>, so the number of calculations can blow up to <span class="math-container">$O(|S|^2)$</span>.</p>

<p><strong>Warning 2:</strong> It is also important to limit the <span class="math-container">$w$</span> parameter because the edge effects propagate through the entire sequence when <span class="math-container">$w=|S|$</span>. Look at these examples built off of OP's initial suggestions of <code>AAAAACCCCC</code> and <code>ACACACACAC</code>. The clustering score never stabilizes when <span class="math-container">$w=|S|$</span>, however it does stabilize away from the edges when <span class="math-container">$w=5$</span>.</p>

<p>For 50 characters of <code>ACACAC...</code>:</p>

<p><a href="https://i.sstatic.net/t8E0M.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/t8E0M.png" alt="enter image description here"></a></p>

<p>For 50 characters of <code>[25xA][25xC]</code>:</p>

<p><a href="https://i.sstatic.net/bPPvA.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/bPPvA.png" alt="enter image description here"></a></p>

CC BY-SA 4.0
2019-02-03T05:24:14.000Z
2019-02-03T15:51:15.817Z
2019-02-03T15:51:15.817Z
conchoecia
<p>Suppose I have a short chain string of oligos, and I want to rate them from 0 to 1 based on how clustered the "A"s are in respect to the chain. For example:</p>

<pre><code>"AAAAACCCCC" #should score as 1 for being as clustered as it possibly could
"ACACACACAC" #should be scored as 0 for being as thinly distributed as it could
</code></pre>

<p>Current idea is to strip both ends of A's, then remove first and last characters, then </p>

<pre><code>gaps = ''.join(str(x) for x in deck).split("A")
listofCs = [x for x in gaps if x] 
</code></pre>

<p>to join it all together as a list of how big the gaps are (i.e. the "C"s), and simply count the average size </p>

<pre><code>average = sum(map(len,listofCs))/len(listofCs)
</code></pre>

<p>But I feel like there should be a much more elegant solution to this. It feels messy. How can I improve the scoring?</p>

CC BY-SA 4.0
2019-01-31T11:24:02.217Z
I would look at how much consecutive the different elements are in the sequence. But how would you rate a sequence like `AACCAACACCAACC` ? Is it for DNA or protein sequences (how much elements do you expect, 4, [5 with an N], 21) ?
llrs
<ngs><algorithms>
How to score how densely bases are distributed along a given oligo?
# How to score how densely bases are distributed along a given oligo?

(Asked by: username_0 on 2019-01-30)

Suppose I have a short chain string of oligos, and I want to rate them from 0 to 1 based on how clustered the "A"s are in respect to the chain. For example:
"AAAAACCCCC" #should score as 1 for being as clustered as it possibly could
"ACACACACAC" #should be scored as 0 for being as thinly distributed as it could

Current idea is to strip both ends of A's, then remove first and last characters, then 
gaps = ''.join(str(x) for x in deck).split("A")
listofCs = [x for x in gaps if x] 

to join it all together as a list of how big the gaps are (i.e. the "C"s), and simply count the average size 
average = sum(map(len,listofCs))/len(listofCs)

But I feel like there should be a much more elegant solution to this. It feels messy. How can I improve the scoring?

**Tags:** ngs, algorithms
> **Comments:**
> **username_1**: I would look at how much consecutive the different elements are in the sequence. But how would you rate a sequence like `AACCAACACCAACC` ? Is it for DNA or protein sequences (how much elements do you expect, 4, [5 with an N], 21) ?  

---

## 3 answers

### Answer by username_2 on 2019-02-01. Score: 2

You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As. You might also need to scale by the sequence length, if you have different length sequences that you want to compare.

Edit: also, there might be a correction to make if one or both ends of the sequence are an A

> **Comments:**
> **username_0**: Interestingly idea. Can you think of any edge cases that would cause problems for this scheme?  
> **username_2**: Zero As would throw a spanner in the works.  

---

### Answer by username_3 on 2019-02-03. Score: 0

The Markov chain approach is very interesting. You could deploy a maximum likelihood approach where instead of scoring the genetic distance within a homologue, the algorithm looks at its physically linked neighbour. The genetic distance of the transition point i.e. AAAA then CCCC is deducted from the final genetic distance so that AAAAACCCCCC would be zero. It would be fairly easy to script and use PAML or equivalent where the default phylogeny would be a star tree.

The challenge is to establish a Deep Learning approach. Its not impossible, you just assign pre-determined classes of oligos which fall into the definition you want and use these as your training set.


---

### Answer by username_4 on 2019-02-03. Score: 0

To define the base density score, or clustering score, of a sequence $S$ we can think of each position in the sequence as being part of a Markov process. For the sake of symmetry around a single base we consider odd-order Markov models from 3-onward where we ask what is the probability that we see the character at $S_i$ in the $j$ bases from $[S_{i-j}, S_{i+j}]$. For each column we have length(S) -1, $|S|-1$, possible probabilities to calculate per window. 
We can't look left or right of the edges of the sequences so the probabilities are skewed, but we can get a rough idea of what is going on.

For each base position we calculate a clustering score of $ \frac{1}{w}\Sigma_{j=1}^{j=w} -log_2 \frac{m(S_{i-j}, S_{i+j})-1 + \frac{1}{|S|}}{|S_{\{i-j, i+j\}}| -1 }$ where:


$w$ is the window size we can look at. The max is $|S|-1$ and the min is $1$. This just means that we can calculate a score in the context of the whole sequence or we maybe only care within some smaller window, like 5 bp flanking either side of the base in question.
The function $m(S_{i-j}, S_{i+j})$ is just the number of counts of matching bases in the window $j$ on either side of the base $S_i$. For example in the sequence gatAgga, if $S_i$ is the uppercase A and $j$ is 3, the function $m(S_{i-3}, S_{i+3})$ returns 3, since there is the A base at $S_i$ and the two additional a bases at $i-2$ and at $i+3$. We have to subtract one from this number because we don't want to factor in the matching base at $S_i$.
We then add $\frac{1}{|S|}$ to that number. This is a pseudocount that we need to add to prevent from taking the log of zero. This number is pretty arbitrary but I chose $\frac{1}{|S|}$ so that non-clusting regions' scores aren't too high.
We then divide that numerator by the length of $S$ from the little window that we just looked at minus one, $|S_{\{i-j, i+j\}}| -1 $. For example in the sequence gatAgga from above we only want to consider the bases gat gga. This is 6 bases.
We then take the negative log to convert the proportional scale to a log scale.
The summation just says that we do this for every symmetrical window size around the base in question then sum them up. For example we would calculate the ratio for tAg, atAgg, and gatAgga.
The $\frac{1}{w}$ just gives the average log clustering score of all the windows that we looked at. One nice thing about this is that the impact of each increasing window size has less effect on the ratio, so there is a sort of convergence that happens.


Implemented in python:
from math import log as mlog
s = "GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA"
assert len(s) > 0
pseudocount = 1/len(s)
#maxdepth = 5
maxdepth = len(s)
score = []
for i in range(len(s)): # iterate positions
    temp = []
    for j in range(1, maxdepth): # iterate depth
        rmin = max(i-j, 0)
        rmax = i+j+1
        num_matches_min_one = sum([1 for k in s[rmin:rmax] if k == s[i]])-1
        len_s_min_one = len(s[rmin:rmax]) -1
        sim = abs(-1 * mlog((num_matches_min_one/len_s_min_one) + pseudocount, 2))
        #print("i={}, j={}, rmin={}, rmax={}, sim={}".format(i, j, rmin, rmax, sim))
        temp.append(sim)
    print("{0:d}\t{1}\t".format(i, s[i]), end='\t')
    score.append(sum(temp)/len(temp))
    print("{0:.4f}".format(score[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x = [i+1 for i in range(len(s))]
plt.bar(x, score)
plt.xticks(x, [k for k in s])
ax.set_yscale("log", nonposy='clip')
plt.xlabel('sequence')
plt.ylabel('-log2(cluster-ness)')
plt.title("cluster-ness score")
plt.savefig('output.png', dpi=150)

This calculates the matrix, prints out a tab-delimited file of the position, base, and log clusting score, then saves a plot of the results.
For the test sequence GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA and using the whole sequence as context, I got the following plot:

As you can see, the log clustering score is low around the things that are cluster-y and high around things that aren't cluster-y.
Note: OP requested that the score be scaled from 0 to 1. I instead used the log scale to make visualization easier.
Warning 1: If you plan on doing this for large sequences, it is really important to limit the maxdepth, or $w$, parameter. The time complexity for this algorithm is $O(|S|*w)$, so the number of calculations can blow up to $O(|S|^2)$.
Warning 2: It is also important to limit the $w$ parameter because the edge effects propagate through the entire sequence when $w=|S|$. Look at these examples built off of OP's initial suggestions of AAAAACCCCC and ACACACACAC. The clustering score never stabilizes when $w=|S|$, however it does stabilize away from the edges when $w=5$.
For 50 characters of ACACAC...:

For 50 characters of [25xA][25xC]:


---

