How to calculate the probability of finding two proteins that share a 5 amino acid long motif from a proteome of around 1067 proteins that have an average length of 65 residues. The probability of a specific sequence that is completely the same is 1067/20^65
$\begingroup$
$\endgroup$
0
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
3
Its just significant at 2%
0.055 . 65 . 1067 = 0.021 (its actually 0.02 - see note)
Where critical probability is 0.05, i.e. 5%
There are 20 aa so 1/20 probability of random occurrence. The oligopeptide is 5 amino acids long, hence to the power 5. Obviously there are 1067 proteins and the probability is per protein, thus if these were independent its a straight multiplication because that is how independent events are measured.
Note The 65 is really 60 because it's a sliding window moving 1 amino acid at a time and the proteins are only 65 in length. The actual value is 60 rather than 65 because at the 3' terminal of the sequence the sliding window can't capture the remaining 4 amino acids, because once the oligopeptide < 5 the calculation is no longer valid. I used 65 so you could see where the value is originating.
Obviously the proteins in question are not independent so you need to state the assumption, i.e. assuming independence.
-
1$\begingroup$ There are 22 amino acids, not 20. Almost all animals (the only exceptions I know are a few insect species) have 21, because of selenocysteine, the 21st amino acid. I know less about pyrrolysine, the 22nd, but I believe it is found in some bacteria and archaea only. Given how few occurreneces of this amino acid most species have, I doubt it will change the analysis much, but there are 21 amino acids in most non-plant eukaryotes. $\endgroup$– terdon ♦Commented Jan 2 at 10:33
-
$\begingroup$ Thanks @terdon, now I know. Selenocysteine is interesting because of the importance of cysteine disulphide bonds. I assume selenocysteine doesn't form them. Anyway, the probability of random occurrence is 1/21 or 1/22 depending. $\endgroup$– M__ ♦Commented Jan 2 at 14:37
-
2$\begingroup$ I... don't remember. Wow. I did my PhD on these things but it's been a while and I always focused on detecting them rather than on their biochemistry, but still! Anyway, Sec has selenium where Cys has sulfur, so presumably it cannot make disulphide bonds by definition, but I do seem to recall that the selenium allows for very similar bonding. We tended to treat Sec as a more reactive version of Cys, but my ignorance of biochemistry is deep so I don't really know. $\endgroup$– terdon ♦Commented Jan 2 at 16:12
$\begingroup$
$\endgroup$
1
This question has two aspects.
1 birthday problem
You state that
The probability of a specific sequence that is completely the same is 1067/20^65
but the probability is a bit larger. You can solve this as the birthday problem. If you have $k = 1067$ random draws out of $n=20^{65}$ proteins, then the probability of at least one double case is approximately $$p \approx \frac{0.5k^2}{d} = \frac{569244.5}{20^{65}} $$
2 Markov chain problem
For two random sequences to contain a similar subsequence the probability is approximately
$$\begin{array}{rcl} P(n \geq 1) &\approx& 1-e^{-E[n]/(1+1\cdot E[S])} \\ E[n] &=& (y-x+1)(z-x+1)(1/k)^x \\ E[S] &=& \frac{1}{k-1} \end{array}$$
Where $x$ is the length of the sub-sequence, $y$ and $z$ are the lengths of the random sequences among which we search for a matching subsequence of length $x$, $k$ is the number of different types of letters out of which the sequences can be composed (e.g. $k=6$ for six-sided dice rolls, $k=4$ for RNA/DNA-sequences, $k=20$ for your proteins), $n$ is the number of sub-sequence matches, $S$ is the length by which the sequence of length $x$ is extended (finding more matches).
With $k=20$, $x=5$ and $y=z=65$ this gives
$$P(n \geq 1) \approx 0.0011$$
This is described in this question on cross validated Probability of a similar sub-sequence of length X in two sequences of length Y and Z and this question on biostars Probability of finding a common sub-sequence of at least 'k' nucleotides between two sequences
Combining the two
For $m=1067$ proteins there are approximately $n(n-1)/2$ pairs and each of those have a $0.0011$ probability of having a match. This give a total probability of approximately
$$P_{total} = 1-(1-P_{individual})^{1067 \cdot 1066/2} = 1-e^{\log(1-P_{individual})*1067\cdot 1067/2} \approx 1-e^{-628} $$
comparing the formula's with simulation
We can verify the computation with a formula by testing it with a simulation for a smaller number of proteins. Say we have 30 proteins, then the probability of at least one match is 0.3815459.
... to be continued, my phone battery is running out.
-
$\begingroup$ Of course, the above is only true for a model where each protein is a complete random draw of a sequence of amino acids. When certain patterns of aminoacids have different probabilities, then matches may be more likely. $\endgroup$ Commented Apr 10 at 22:06