4
$\begingroup$

Say we have two random sequences of length $y=3300$ and length $z=26$, where the sequences are made out of a combination of $k=4$ four letters ACGT.

  • the letter in each position of the sequence is independent from the other letters
  • the letter in each position in the sequence and can have equal probability for A, C, G or T. That is, every possible sequence has equal probability.

(Note:this relates to RNA sequences where the nucleotides are completely random iid and with equal probability for each of the four nucleotides)

What is the probability that both these two sequences contain a same subsequence of length $x=19$? Or what is an approximation of this probability?

A related problem is: A fair die is rolled 1,000 times. What is the probability of rolling the same number 5 times in a row?

but here we are not looking for the probability of rolling a specific sequence (like the same number five times in a row) but instead for any possible subsequence that occurs in another random sequence of some length.

$\endgroup$
0

1 Answer 1

0
$\begingroup$

An approximation.

  • Step 1: Compute the expectation value for the number of similar sub sequences.

    There are $y-x+1$ and $z-x+1$ subsequences. The probability for a match between two subsequences is $(1/4)^x$, then the expectation value for the number of matches $n$ is $$E[n] = (y-x+1)(z-x+1)(1/4)^x$$

  • Step 2 (a dead end, wrong step that will be adjusted in step 3): If we are naive, then we consider the distribution of the number of matches to be Poisson distributed. Then we put the expectation value of the step 1 approximately equal to the probability of having at least 1 or more sub sequence matches by using the Poisson distribution.

    $$P(n \geq 1) \approx 1-e^{-E[n]}$$

    But this is false. See for instance the linked question A fair die is rolled 1,000 times. What is the probability of rolling the same number 5 times in a row? where the probability of 5 times at least the same number in a row is not $1- \exp \left(-(1000-5+1)\left(\frac{1}{6^4}\right) \right)\approx 0.536$ and instead it is approximately $0.474$

    The reason that it is false, is because when we already have one matching subsequence, then the probability of another matching subsequence becomes more likely. The distribution of the number of matches $n$ is not a Poisson distribution.

  • Step 3 (adjustment for error in step 2): we need to figure out an estimate for the relative probabilities of 1 match, 2 matches, etc. and use this along with the expectation.

    Say that we already have one matching subsequence of length $x$, then another match is found when either on the left or right of the subsequences we can find another matching letter.

    If we ignore the situations where the subsequences are found at the boundary (if we assume that we can extend the subsequences indefinitely), say we have a 19 length sub-sequence like CTACGTGCCCGCCGAGGAG in both two sequences, as below where the * are wild cards and can be with 1/4-th probability any letter A, C, G or T.

    sequence 1: *******CTACGTGCCCGCCGAGGAG******
    sequence 2: *******CTACGTGCCCGCCGAGGAG******
    

    There is a 1/4 probability that the match extends with another pair on the right. There is a 1/4 probability that the match extends with another pair on the left. So the distribution of additional matches, $S$, will be geometric distributed, or more precisely the sum of two geometric distributions (on both sides). Then the expectation value can be expressed approximately as:

    $$E[n_{true}] \approx E[n_{naive}] \cdot (1+2\cdot E[S])$$

    where $E[S] = 1/(k-1)$ is the expectation value of the geometric distributed number of extensions $S$ (and $k$ is the number of options).

    Then we can use the approximation

    $$P(n \geq 1) \approx 1-e^{-E[n]/(1+2\cdot E[S])}$$

  • Step 4 (adjustment for error in step 3)

    When we apply the last expression to the dice roll situation, then we get the estimate $$1- \exp \left(-(1000-5+1)\left(\frac{1}{6^4}\right) /1.4 \right)\approx 0.422$$ and we get problems of over-compensating. Interrestingly when we do $$1- \exp \left(-(1000-5+1)\left(\frac{1}{6^4}\right) /1.2 \right)\approx 0.473$$ then we get very close to the right answer.

    We can see that the correction from step 3 is counting the cases multiple times. For instance, some subsequence in position $i$ is extended to the right and some subsequence in position $i+1$ is extended to the left are the same cases.

    So when we do the extension only to one single side, and use only a factor 1 instead of 2, then we get a better approximation because we are not counting cases multiple times. So the final result is: $$\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/4)^x \\ E[S] &=& \frac{1}{k-1} \end{array}$$ Where $x$ is the length of the sub-sequence, $y$ and $z$ are 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, or $k=4$ for RNA/DNA-sequences), $n$ is the number of sub-sequence matches, $S$ is the length by which the sequence of length $x$ is extended (finding more matches).

Conclusion

$$P \approx 1- e^\frac{3275 \cdot 8 \cdot 0.25^{19}}{1+1/3} \approx 7.14 \cdot 10^{-8}$$

Comparison with simulation

Below we use $x = 9$, $y=200$, $z=180$, $k=4$ and compare the formula with a simulation of one million pairs of sequences.

The differences are small with 0.09015627 for the estimate from the formula and 0.090055 for the estimate from the simulation.

sim = function(n1=100,n2=100, l = 9){   
  ### create two sequences   
  s1 = sample(c("a","c","g","t"), n1, replace=TRUE)
  s2 = sample(c("a","c","g","t"), n2, replace=TRUE)   
  
  ### perform a loop that searches with grepl for every possible substring from s1 in the string made from s2    
  count = 0   
  string2 = paste0(s2, collapse = "")   
  for (i in 1:(n1-l+1)) {
    substring = paste0(s1[c(1:l)+i-1], collapse = "")
    if (grepl(substring, string2)) {
      count = 1
    }   
  }   
  return (count)
}

x = 9 
y = 180 
z = 200 
k = 4

### this returns 0.09015627
1-exp(-(z-x+1)*(y-x+1)*(1/k)^x*(k-1)/k)

### this returns 0.090055 
set.seed(1)
nt = 10^6 
sum(replicate(nt,sim(y,z,x)))/nt
$\endgroup$
5
  • $\begingroup$ Since you're apparently interested in approximations, please so indicate in the question. $\endgroup$
    – whuber
    Commented Feb 25, 2022 at 13:39
  • 1
    $\begingroup$ @whuber I was initially not looking for an approximation (I thought it could be solved with some method as in the linked question). But, the current solution seems like an elegant approximation, and I don't believe that improvements can be made except for very long direct computations. I am curious to see whether I am wrong about that, do you have ideas whether there is something in between? $\endgroup$ Commented Feb 25, 2022 at 14:11
  • $\begingroup$ I think you're completely correct about that. It is clear that in the particular case you posit, where the alphabet is small, one string is long, and the other very much shorter, then highly accurate asymptotic approximations can easily be derived. The interest, then, seemed to be in the possibility you contemplated applications that didn't fit those specific circumstances and where asymptotic results would not work very well. $\endgroup$
    – whuber
    Commented Feb 25, 2022 at 14:17
  • 4
    $\begingroup$ @whuber the contemplated application is MSH3 Homology and Potential Recombination Link to SARS-CoV-2 Furin Cleavage Site which, I believe, got the computation in figure 2 completely wrong. It is the wrong computation for the wrong problem. Wrong computation because it should be like in this question. Wrong problem because it computes a probability that is not relating to the problem (like this dice rolling problem). So they computed $p=1.23*10^{-11}$ instead of $p=0.0035$ $\endgroup$ Commented Feb 26, 2022 at 12:08
  • $\begingroup$ In that case, a reasonably tight lower bound would serve you well. $\endgroup$
    – whuber
    Commented Feb 26, 2022 at 13:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.