# Probability of finding a particular sequence of base pairs

Thinking about probability always makes me realize how bad I am at counting...

Consider a sequence of $n$ base letters $A,\; T, \; C, \text{ and } G$, each equally likely to appear. What is the probability that this sequence contains a particular sequence of base pairs of interest of length $r\leq n$?

There are $4^n$ different (equally likely) sequences possible. Start with the sequence of interest at the beginning of the full sequence; $4^{n-r}$ sequences like this are possible. We can start our sequence of interest in $n+1 -r$ different locations. Hence, my answer is $(n+1-r)/4^r$.

This probability is increasing in $n$, which makes sense to me. But this probability exceeds 1 when $n>4^r +r-1$. But that can't be. The probability should approach 1 in the limit (seems to me), but not exceed it.

I assume that I'm double counting something. What am I missing? Thanks.

(FYI, not homework, just a toy example in preparation for exams. A question posed by my molecular biologist friend.)

Let's contemplate a small version of this problem with $n=5$. What is the chance a sequence of five letters will contain the target $\ldots A C G T\ldots$? This is easy: $4^{-4}$ of all sequences begin with this string, another $4^{-4}$ end with it, and no sequence both begins and ends with this string. Therefore the chance is $2 \times 4^{-4}$.

On the other hand, what is the chance of $\ldots A A A A \ldots$? Once again, $4^{-4}$ of the sequences begin with this string, the same proportion end with this string, and $4^{-5}$ of all sequences do both. Therefore, by the Principle of Inclusion-Exclusion, the answer is $2 \times 4^{-4} - 4^{-5}$.

In general, the answer depends on the structure of the substring. To be more specific, when you're scanning a string (from left to right, say) for $ACGT$, you ignore all characters until you see that initial $A$. After that, there are three possibilities: the next character is a match for $C$, the next one is a non-match for $C$ but is not an $A$ (so you are back in the wait-for-an-$A$ state), or the next one is a non-match yet it's an $A$, placing you into the just-saw-an-$A$ state. In contrast, consider a search for $ACTACG$. Suppose you have seen the prefix $ACTAC$. The next character will match if it is $G$. When it's a non-match, (i) a $C$ puts you into the initial wait-for-an $A$ state, (ii) an $A$ has you watching out for a $C$, and (iii) a $T$ means you have already seen $\ldots ACT$ and you're already halfway to a match (and looking for the second $A$). The relevant "structure" evidently consists of patterns of substrings in the target that match the prefix of the target. That's why the chances depend on the target string.

The FSA diagrams I advocate in a reply at Time taken to hit a pattern of heads and tails in a series of coin-tosses can help with understanding this phenomenon.

A crude approximation would be $1-(1-1/4^r)^{n-r+1}$. You take the probability that your sequence do not occur at a particular location, put it to the power of the number of locations (falsely assuming independence), which is $n-r+1$ not $n-r$, and this is an approximation of its not occurring so you then need to subtract this from $1$.

A precise calculation will depend on the precise pattern you are looking for. $AAAAA$ is more likely to not occur than $ATCGT$.

• Maybe it's just me, but $1-(1-(1/4)^r)^{n-(r-1)}$ seems a bit clearer in terms of understanding how the equation was constructed.
– user136692
Nov 2, 2016 at 12:06
• @JoeRocc - I suspect this is personal. If you read from page $300$ through to page $400$ of a book, have you read $400-300+1=101$ pages or $400-(300-1)=101$ pages? Nov 2, 2016 at 13:05
• No worries, I was only going by my intuition of the problem. If we intuitively derive an equation to be $(a-(b-(c-1+d)))$, then when trying to explain it to someone I think it's best to leave it as that rather than to simplify it to $a-b+c-1+d$ (though this certainly may turn out to be more intuitive upon consideration). Your intuition may have been different in any case :)
– user136692
Nov 2, 2016 at 13:51

It is possible to obtain the exact probability of a particular subsequence by using a Markov chain representation of the problem. The specifics of how to construct the chain depend on the particular subsequence of interest, but I will give a couple of examples of how to do this.

Exact probability via Markov chain: Consider a discrete sequence of outcomes of $$A,T,C,G$$ where the outcomes in the sequence are exchangeable, and suppose we are interested in some substring of length $$k$$. For any given value of $$n$$, let $$\mathscr{W}$$ be the event that the substring of interest occurs, and let $$\mathscr{H}_a$$ be the event that the last $$a$$ outcomes are the first $$a < k$$ characters in the substring of interest (but no more than this). We use these events to give the following partition of $$k+1$$ possible states of interest:

$$\begin{matrix} \text{State 0} & & & \bar{\mathscr{W}} \cap \mathscr{H_0}, \text{ } \text{ } \text{ } \\[6pt] \text{State 1} & & & \bar{\mathscr{W}} \cap \mathscr{H_1}, \text{ } \text{ } \text{ } \\[6pt] \text{State 2} & & & \bar{\mathscr{W}} \cap \mathscr{H_2}, \text{ } \text{ } \text{ } \\[6pt] \text{State 3} & & & \bar{\mathscr{W}} \cap \mathscr{H_3}, \text{ } \text{ } \text{ } \\[6pt] \vdots & & & \vdots \\[6pt] \text{State }k-1 & & & \bar{\mathscr{W}} \cap \mathscr{H_{k-1}}, \\[6pt] \text{State }k & & & \mathscr{W}. \quad \quad \quad \text{ } \text{ } \\[6pt] \\[6pt] \end{matrix}$$

Since the sequence of outcomes is assumed to be exchangeable we have independent outcomes conditional on their respective probabilities $$\theta_A + \theta_T + \theta_C + \theta_G = 1$$. Your process of interest can be represented as a discrete-time Markov chains that begins in $$\text{State 0}$$ at $$n=0$$ and transitions according to a probability matrix that depends on the particular substring of interest. The transition matrix will always be a $$(k+1) \times (k+1)$$ matrix representing the probabilities of transition using the above states. If the substring of interest has not been reached then each transition can either bring you one step closer to the substring or it can set you back to a previous state that depends on the particular substring. Once the substring is reached, this is an absorbing state of the chain, representing the fact that the event of interest has occurred.

For example, if the substring of interest is $$AAAAAA$$ then the transition matrix is:

$$\mathbf{P} = \begin{bmatrix} 1-\theta_A & \theta_A & 0 & 0 & 0 & 0 & 0 \\[6pt] 1-\theta_A & 0 & \theta_A & 0 & 0 & 0 & 0 \\[6pt] 1-\theta_A & 0 & 0 & \theta_A & 0 & 0 & 0 \\[6pt] 1-\theta_A & 0 & 0 & 0 & \theta_A & 0 & 0 \\[6pt] 1-\theta_A & 0 & 0 & 0 & 0 & \theta_A & 0 \\[6pt] 1-\theta_A & 0 & 0 & 0 & 0 & 0 & \theta_A \\[6pt] 0 & 0 & 0 & 0 & 0 & 0 & 1. \\[6pt] \end{bmatrix}$$

Contrarily, if the substring of interest is $$ACTAGC$$ then the transition matrix is:

$$\mathbf{P} = \begin{bmatrix} 1-\theta_A & \theta_A & 0 & 0 & 0 & 0 \\[6pt] 1-\theta_A-\theta_C & \theta_A & \theta_C & 0 & 0 & 0 & 0 \\[6pt] 1-\theta_A-\theta_T & \theta_A & 0 & \theta_T & 0 & 0 & 0 \\[6pt] 1-\theta_A & 0 & 0 & 0 & \theta_A & 0 & 0 \\[6pt] 1-\theta_A-\theta_C-\theta_G & \theta_A & \theta_C & 0 & 0 & \theta_G & 0 \\[6pt] 1-\theta_A-\theta_C & \theta_A & 0 & 0 & 0 & 0 & \theta_C \\[6pt] 0 & 0 & 0 & 0 & 0 & 0 & 1. \\[6pt] \end{bmatrix}$$

As can be seen above, constructing the transition matrix requires attention to the particular substring. An incorrect outcome sets you back to a previous state in the string that depends on the particular substring of interest. Once the transition matrix is constructed, for a given value of $$n$$ the probability of having the substring in the chain is $$\mathbb{P}(\mathscr{W} | n) = \{ \mathbf{P}^n \}_{0,k}$$. (This probability is zero for all $$n.)

Programming this in R: You can program this as a function in R by creating a function that generates the transition matrix for the Markov chain and an array of its powers up to some desired number of trials. You can then read the appropriate transition probability for the value of $$n$$ that is of interest. Here is an example of some code to do this:

#Create function to give n-step transition matrix for n = 1...N
#We will use the example of the substring of interest "AAAAAA"

#a is the probability of A
#t is the probability of T
#c is the probability of C
#g is the probability of G
#N is the last value of n
PROB <- function(N,a,t,c,g) { TOT <- a+t+c+g;
a <- a/TOT;
t <- t/TOT;
c <- c/TOT;
g <- g/TOT;

P <- matrix(c(1-a, a, 0, 0, 0, 0, 0,
1-a, 0, a, 0, 0, 0, 0,
1-a, 0, 0, a, 0, 0, 0,
1-a, 0, 0, 0, a, 0, 0,
1-a, 0, 0, 0, 0, a, 0,
1-a, 0, 0, 0, 0, 0, a,
0, 0, 0, 0, 0, 0, 1),
nrow = 7, ncol = 7,
byrow = TRUE);
PPP <- array(0, dim = c(7,7,N));
PPP[,,1] <- P;
for (n in 2:N) { PPP[,,n] <- PPP[,,n-1] %*% P; }
PPP }

#Calculate probability for N = 100 for equiprobable outcomes
N <- 100;
a <- 1/4;
t <- 1/4;
c <- 1/4;
g <- 1/4;
PROB(N,a,t,c,g)[1,7,N];

 0.01732435


As you can see from this calculation, the probability of getting the substring $$AAAAAA$$ in $$n=100$$ tosses with equiprobable outcomes is $$0.01732435$$. This is just one example using a particular substring and a given number of trials, but it can be varied to obtain probabilities with respect to other substrings of interest.

• That's a really cool application of Markov chains. I made a script, which automatically constructs the transition matrices given a search string. I hope it is correct. It works at least for the two cases given in your post. Btw, the probability to find the BZ2 magic bit string 0x314159265359 in 1 TiB of random bits, seems to be ~3%. Nov 29, 2020 at 13:58
• @mxmlnkn: Since you found this example interesting, you might be interested in some of the broader statistical literature on pattern Markov chains and deterministic finite automata (DFAs) (see e.g., Nuel 2008). The function you are trying to construct here is the DFA of an input string.
– Ben
Nov 29, 2020 at 21:30

You are double counting the sequences that include several times your target subsequence, for example both at position A and at position B!=A. That's why your erroneous probability can exceed 1

• Very well done ! +1 Jan 12, 2017 at 21:33