<p>Its just significant at 2%</p>
<blockquote>
<p>0.05<sup>5</sup> . 65 . 1067 = 0.021 (<em>its actually 0.02 - see note</em>)</p>
</blockquote>
<p>Where critical probability is 0.05, i.e. 5%</p>
<p>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.</p>
<p><strong>Note</strong> 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 &lt; 5 the calculation is no longer valid. I used 65 so you could see where the value is originating.</p>
<p>Obviously the proteins in question are not independent so you need to state the assumption, i.e. assuming independence.</p>

CC BY-SA 4.0
2024-01-02T10:33:38.313Z
There are 22 amino acids, not 20. Almost all animals (the only exceptions I know are [a few insect species](https://pubmed.ncbi.nlm.nih.gov/18698431/)) have 21, because of [selenocysteine](https://en.wikipedia.org/wiki/Selenocysteine), the 21st amino acid. I know less about [pyrrolysine](https://en.wikipedia.org/wiki/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.
terdon
CC BY-SA 4.0
2024-01-02T14:37:43.217Z
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.
M__
CC BY-SA 4.0
2024-01-02T16:12:00.090Z
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 di*sulphide* 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.
terdon
CC BY-SA 4.0
2024-01-02T04:10:50.577Z
2024-01-02T04:15:51.257Z
2024-01-02T04:15:51.257Z
M__
<p>This question has two aspects.</p>
<h3>1 birthday problem</h3>
<p>You state that</p>
<blockquote>
<p>The probability of a specific sequence that is completely the same is 1067/20^65</p>
</blockquote>
<p>but the probability is a bit larger. You can solve this as the <a href="https://en.m.wikipedia.org/wiki/Birthday_problem" rel="nofollow noreferrer">birthday problem</a>. If you have <span class="math-container">$k = 1067$</span> random draws out of <span class="math-container">$n=20^{65}$</span> proteins, then the probability of at least one double case is approximately <span class="math-container">$$p \approx \frac{0.5k^2}{d} = \frac{569244.5}{20^{65}} $$</span></p>
<h3>2 Markov chain problem</h3>
<p>For two random sequences to contain a similar subsequence the probability is approximately</p>
<blockquote>
<p><span class="math-container">$$\begin{array}{rcl}
P(n \geq 1) &amp;\approx&amp; 1-e^{-E[n]/(1+1\cdot E[S])} \\
E[n] &amp;=&amp; (y-x+1)(z-x+1)(1/k)^x \\
E[S] &amp;=&amp; \frac{1}{k-1}
\end{array}$$</span></p>
<p>Where <span class="math-container">$x$</span> is the length of the sub-sequence, <span class="math-container">$y$</span> and <span class="math-container">$z$</span> are the lengths of the random sequences among which we search for a matching subsequence of length <span class="math-container">$x$</span>, <span class="math-container">$k$</span> is the number of different types of letters out of which the sequences can be composed (e.g. <span class="math-container">$k=6$</span> for six-sided dice rolls, <span class="math-container">$k=4$</span> for RNA/DNA-sequences, <span class="math-container">$k=20$</span> for your proteins), <span class="math-container">$n$</span> is the number of sub-sequence matches, <span class="math-container">$S$</span> is the length by which the sequence of length <span class="math-container">$x$</span> is extended (finding more matches).</p>
</blockquote>
<p>With <span class="math-container">$k=20$</span>, <span class="math-container">$x=5$</span> and <span class="math-container">$y=z=65$</span> this gives</p>
<p><span class="math-container">$$P(n \geq 1) \approx 0.0011$$</span></p>
<p>This is described in this question on cross validated <a href="https://stats.stackexchange.com/questions/565676/">Probability of a similar sub-sequence of length X in two sequences of length Y and Z</a> and this question on biostars <a href="https://www.biostars.org/p/249573/" rel="nofollow noreferrer">Probability of finding a common sub-sequence of at least 'k' nucleotides between two sequences</a></p>
<h3>Combining the two</h3>
<p>For <span class="math-container">$m=1067$</span> proteins there are approximately <span class="math-container">$n(n-1)/2$</span> pairs and each of those have a <span class="math-container">$0.0011$</span> probability of having a match. This give a total probability of approximately</p>
<p><span class="math-container">$$P_{total} = 1-(1-P_{individual})^{1067 \cdot 1066/2} = 1-e^{\log(1-P_{individual})*1067\cdot 1067/2} \approx 1-e^{-628} $$</span></p>
<h3>comparing the formula's with simulation</h3>
<p>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.</p>
<p>... to be continued, my phone battery is running out.</p>

CC BY-SA 4.0
2024-04-10T22:06:58.427Z
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.
Sextus Empiricus
CC BY-SA 4.0
2024-04-10T22:00:46.000Z
2024-04-10T22:28:16.777Z
2024-04-10T22:28:16.777Z
Sextus Empiricus
<p>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</p>

<phylogenetics><proteins><r><sequence-analysis><homework>
probability of finding a 5 amino acids in a row within a proteome
# probability of finding a 5 amino acids in a row within a proteome

(Asked by: username_0 on 2024-01-02)

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

**Tags:** phylogenetics, proteins, r, sequence-analysis, homework

---

## 2 answers

### Answer by username_1 on 2024-01-02. Score: 2 | Accepted answer

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.

> **Comments:**
> **username_2**: There are 22 amino acids, not 20. Almost all animals (the only exceptions I know are [a few insect species](https://pubmed.ncbi.nlm.nih.gov/18698431/)) have 21, because of [selenocysteine](https://en.wikipedia.org/wiki/Selenocysteine), the 21st amino acid. I know less about [pyrrolysine](https://en.wikipedia.org/wiki/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.  
> **username_1**: 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.  
> **username_2**: 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 di*sulphide* 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.  

---

### Answer by username_3 on 2024-04-10. Score: 0

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](https://en.m.wikipedia.org/wiki/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](https://stats.stackexchange.com/questions/565676/) and this question on biostars [Probability of finding a common sub-sequence of at least 'k' nucleotides between two sequences](https://www.biostars.org/p/249573/)

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

> **Comments:**
> **username_3**: 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.  

---

