# What is the expected number of occurrences of a 9-mer in 500 random DNA strings, each of length 1000?

I have started learning bioinformatics. There are some matter of finding expected value. But I think I am very weak in calculating such types of things.

As expected value is related to statistics, its explanation is skipped in bioinformatics. So, I am posting it here.

### Question:

Suppose, I have 500 strings, each having length 1000.

Now, I have to calculate the expected number of occurrences of a sub-string having length exactly 9.

Notice that, the string contains only four letters A, T, G, C with same probability (each 0.25).

Another thing to be noted: Overlapping strings should be counted.

### My Approach:

The probability of existing a 9-length sub-string among all 9-length sub-strings = $$(0.25)^9$$

The number of occurrences of a 9-length sub-string in a string having length 1000 = $$(1000-9+1) * (0.25)^9$$

If the number of such string becomes 500, then the number of occurrences would be = $$500 * (1000-9+1) * (0.25)^9$$

But I did wrong somewhere, may be in assumption or in calculation.

Could you please guide me to get the actual solution?

### Source:

This problem is a part of Bioinformatics course track in Coursera.

### Accuracy:

Allowable error = 0.0001

As allowable error is 0.0001, the given calculation serves the purpose and gives a good approximation. It was my bad that I entered less digits there and got that wrong.

The answer is: $$1.8920898$$

However, This answer gives an approximation about the probability. But when it is converted to expected value by multiplying, it becomes a little bit bad and does not serve the purpose. It gives answer: $$1.8885179$$. According to whuber ♦'s calculation in comment, it came $$1.895678$$ which also does not serve the purpose.

• It depends on the details of each string. Possibly you are thinking of your 500 strings as being some kind of random sample from a large set, but if so, exactly what is that large set? If it's the set of all strings with the uniform distribution then the answer still depends on the length-9 probe string. (That surprises some people.) This question has been discussed a few times on this site (e.g., stats.stackexchange.com/questions/26988), so do some searching.
– whuber
Aug 19, 2016 at 18:15
• @whuber primarily, we should only think about this problem without considering that large dataset. Aug 19, 2016 at 18:19
• There is no way to answer it without creating a model for the probability distribution of those 500 strings. For instance, in one model, A never appears. The expected count of AAAAAAAAA therefore is zero. In another model, all strings start and end with 100 A's. The expected count of AAAAAAAAA therefore is at least 184.
– whuber
Aug 19, 2016 at 18:22
• @whuber I got this problem from here. I did something wrong with precision. According to them, my solution is correct. They allow error up to 0.0001. So, may be you are telling about something which deals with the expected number less than that. However, I am actually asking about expected number, not only probability. That means it is not exact duplicate of the question you have provided. Thanks. Aug 19, 2016 at 18:28
• Henry's answer gives 1.895678. You have reported 500 times the probability, but the probability is not the expectation. Unfortunately, because your link is inaccessible, we don't know how they arrived at the solution or what assumptions they made to do it--we only know they had to make some simplifying assumptions--so it isn't even clear (a) what the "correct" answer would be or (b) that they are reporting the correct answer!
– whuber
Aug 19, 2016 at 18:55

Finding the expectation is not so difficult.

The total number of substrings is

$$N = n_\text{strings} \cdot (1+l_\text{string}-l_\text{substring}) = 500 \cdot 992 = 496000$$

For each single substring the expectation of a match is (assuming identical and independent distribution for the occurrence of the letters in the strings)

$$E[\text{single substring match}] = \left( \frac{1}{4} \right)^9$$

And the expectation for all substrings is simply a multiple of this

$$E[\text{number of matches}] = N \cdot E[\text{single substring match}] = 496000 \cdot \left( \frac{1}{4} \right)^9 = \frac{7750}{4096} \approx 1.89208984375$$

These problems tend to become more difficult when you try to find more than just the expectation.

The expectation is not so difficult compared to finding the distribution. In the case of the distribution you have to take into account potential correlations between adjacent substrings. This occurs in problems like:

What is the probability for an N-char string to appear in an M-length random string?

Probability of a similar sub-sequence of length X in two sequences of length Y and Z

Probability of finding a particular sequence of base pairs

That was about a different question. The probability of finding at least one match. That is not the same as the expectation.

To find the expectation you have to sum the probability to get at least one match, the probability to get at least two matches, the probability to get at least 3 matches, etc.

See also the here the properties of the expectated value and how you can express it in terms of the CDF.

This is only a partial answer because it doesn't meet the accuracy requirements, however I'll explain where the answer comes from. As has been pointed out by whuber♦ above in the comments and in this post, the specific content of the strings influence the true probability of the event described in the OP.

Instead of analytically working out these many details, one might try a simulated approach. The following is an implementation in Python of an algorithm that simply generates k-mers of length 9 and 500 sequences of length 1000, counts occurrences of the k-mer in the other sequences, and puts the resulting count into a list. There are improvements that can be made to this code to improve the time and space requirements, but I think this simpler example communicates the process adequately.

from random import choice
import numpy as np

def randseq(n):
return ''.join([choice(['A', 'T', 'C', 'G']) for i in range(n)])

counts = []
i = 0
while len(counts) < 5 or 3 * np.std(counts) / np.sqrt(i) > 0.0001:
i += 1
kmer = randseq(9)
kcount = 0
for j in range(500):
seq = randseq(1000)
if kmer in seq:
kcount += seq.count(kmer)
counts.append(kcount)
print(i, kcount, np.mean(counts), 3 * np.std(counts) / np.sqrt(i))


Using the above, I found that a sample of $$n=14199$$ iterations of the 'experiment' yielded a mean of $$1.8833016409606311 \pm 0.03492470125969519$$ (3 times standard error). I stopped the calculations at this point because I realized that the mean would converge very slowly. The standard error of the mean is $$s_{\bar{x}} = \frac{1}{\sqrt{n}} \sqrt{\frac{(\sum x_i - \bar{x})^2}{n}} = \frac{s_x}{\sqrt{n}}$$

where $$n$$ is the sample size, $$x_i$$ is an instance of the random variable representing the number of occurrences of the k-mer in the other sequences, and $$\bar{x}$$ is the estimate of the mean of this random variable.

Since the mean and standard deviation are consistent estimators, we need only have a sufficiently large $$n$$ to reach an arbitrary precision in $$\bar{x}$$ as calculated in terms of $$s_{\bar{x}}$$. Thinking in terms of confidence intervals, perhaps we'd be satisfied with $$3 s_{\bar{x}} = 10^{-4}$$ being a way of declaring whether we've met the accuracy requirement, giving us:

$$\frac{s_x}{\sqrt{n}} = \frac{10^{-4}}{3}$$

Given that my data had $$s_x = 1.387203975627196$$ at $$n=14199$$, we'll use that value for a ball-park estimate. Substituting for $$s_x$$, and solving for $$n$$, I would expect that the required sample size using this approach would be $$n=1731901383$$. My approach here doesn't work to answer the question in a reasonable amount of time, but would eventually yield the correct answer.

### Update

I re-ran the simulations and obtained a closer $$3 \sigma$$ estimate of $$\bar{x} = 1.8897422883888795 \pm 0.0043822154484445775$$ where $$n=886650$$.