4

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.

whuber
  • 322,774
Enamul Hassan
  • 143
  • 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., http://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
  • @whuber Henry's answer gives output 1.8885179 which does not pass. My solution is 1.8920898 which passes. – Enamul Hassan Aug 19 '16 at 18:42
  • 1
    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 '16 at 18:55
  • @whuber So, how did you get 1.895678 from Henry's answer? However, your calculated answer also did not pass. – Enamul Hassan Aug 19 '16 at 19:09

2 Answers2

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.

Sextus Empiricus
  • 77,915
0

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

Galen
  • 8,442