# 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

• 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 '16 at 18:15
• @whuber primarily, we should only think about this problem without considering that large dataset. – Enamul Hassan Aug 19 '16 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 '16 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. – Enamul Hassan Aug 19 '16 at 18:28
• A solution that shows how to find the full probability distribution automatically includes everything you need to find expectations. Sometimes there are shortcuts to finding expectations, without needing the full distribution, but a close reading of the duplicate should reveal that cannot be the case here. The fact that you need only an approximate answer, though, opens up genuinely new possibilities. (But please see Henry's answer to the duplicate first: it appears to solve your problem.) Please include that information in an edit to your post. – whuber Aug 19 '16 at 18:30

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.

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 still isn't clear what question this answers! Please edit your post to show all the assumptions that must be made to obtain a definite answer. – whuber Aug 19 '16 at 19:51