How to score how densely bases are distributed along a given oligo?

Question

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?

Answer 1

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

Answer 2

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 3

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