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

Question

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.

Answer 1

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

However, This answer gives an approximation about the probability.

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.

Answer 2

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