# Binary classification of DNA sequence motifs

[I apologize in advance for the length of this question. but it seems unavoidable.]

I'm trying to do binary classification, using Bayesian methods, and I'm having some trouble figuring what methods I need to do so. Bottom line - I need some kind of decision rule. I've tentatively settled on a version of the Bayes factor, outlined below, but am open to other suggestions. I know little about Bayesian statistics or decision theory, and this project has become to some extent an exercise in learning some Bayesian theory. My formal statistical education has been rather frequentist in nature. I also know a little about machine learning, but that doesn't seem relevant here.

Here is a description of the background. The data is a set of sequences of letters of equal length, say $n$. The letters are from some alphabet with $r$ letters. The specific application is for DNA sequences i.e. strings composed from the four letters A, C, G, T. So here, $r=4$. The aim is to use the elements of such a set of DNA sequences (called motifs) to predict which other sequences belong to the family.

Personal/historical note: this was started some years ago as a project in motif finding, and was abandoned. I've recently picked it up again. The approach outlined here is outside the mainstream of motif finding, as far as I can tell, because motif finding is mostly concerned with finding new motif families, not extending existing ones. This approach could be used to find new members of existing families, but I don't see how it could be used to find new motif families.

Let $\mathbf{X}=(X_1,\dots, X_n)$ be a random vector. We model these sequences as a set of i.i.d. random vectors, $\{\mathbf{X^{(i)}}=(X_1^{(i)},\dots, X_n^{(i)})$, $i\in \mathcal{I}\}$, where $\mathbf{X^{(i)}}\sim\mathbf{X}$ $\forall i\in \mathcal{I}$.

We define a family of probability distributions for X. We first try to find the member of this family which is the best fit for the data, and then use it for prediction.

Our assumption is that $\mathbf{X}$ can be divided into disjoint subsets, such that variables in each subset are dependent among themselves, and the subsets are mutually independent.

It suffices to define the joint probability distribution on one of these subsets. Let $\mathbf{S}=(X_{j_1},\dots, X_{j_t})$ be a sub-vector of $\mathbf{X}$ for some subset $\{j_1,\dots, j_t\} \in \{1,2\dots, n\}$. So, $\mathbf{S}$ takes $k=r^t$ possible values. Then the joint probability distribution is defined as a categorical probability distribution on $\mathbf{S}$.

\begin{align}\label{prob} P(S=(s_{1}, s_{2}, \dots, s_{t})) = p_{s_{1}s_{2}\dots s_{t}} \end{align}

So the distribution of $\mathbf{S}$ is defined by $k$ probabilities.

Observe that the random variables $\{X_{j_1},\dots, X_{j_t}\}$ are not independent for general values of this probabilities.

Let us restrict our attention to the sequence of subsets of $\{\mathbf{X^{(i)}}, i\in \mathcal{I}\}$, corresponding to $\mathbf{S}$, which we call $\mathbf{Y} = \{\mathbf{S}^{(i)}, i\in \mathcal{I}\}$. Observe that this is a sequence of i.i.d. categorical variables. Then the likelihood of $\mathbf{Y}$ is \begin{align*} P(\mathbf{Y}|\mathbf{p}, \mathcal{M}) = \prod_{i=1}^k p_i^{f_i} \end{align*} where $(f_1, \dots, f_k)$ are the frequencies of the $(s_{1}, s_{2}, \dots, s_{t})$ in Equation~\ref{prob} for some fixed ordering, and $\mathcal{M}$ is the chosen model, namely the division of $\mathbf{X}$ into subsets. The corresponding uninformative prior for $\mathbf{p}$ is \begin{align*} \pi(\mathbf{p}) = (k-1)! I(\textstyle\sum_i p_i = 1) \end{align*} This is a special case of the Dirichlet distribution

\begin{align*} \pi(\mathbf{p}) = \frac{\Gamma(\sum_i{\alpha_i})}{\prod_i\Gamma(\alpha_i)} \prod_{i\in\mathcal{I}} p_i^{\alpha_i -1} I ( \textstyle\sum_i p_i = 1 ) \end{align*} with $\alpha_i=1$ $\forall i$.

Now, one can use Bayesian model selection to select a suitable model, by searching for the model that maximize $P(\mathbf{Y}|\mathcal{M})$. Let us assume we have found a suitable candidate for this model, namely $\mathcal{M}_1$. Call the corresponding uncorrelated model $\mathcal{M}_0$.

The idea is to try to discriminate between $\mathcal{H}_1=\mathcal{M}_1$ and $\mathcal{H}_0=\mathcal{M}_0$, to try to use this model for prediction. To test this, we use k-fold cross-validation. We divide a data set of related sequences (motifs), into $K$ equal parts. We use $K-1$ of these parts as a training set, and the remaining part as a testing set.

Then we apply a decision rule to every element in the testing set to determine whether the element is a better fit for $\mathcal{M}_1$ or $\mathcal{M}_0$.

The obvious approach is to use the Bayes factor, namely

\begin{align}\label{bayesfact} K = \frac{P(D|\mathcal{M}_1)} {P(D|\mathcal{M}_0)} = \frac{\int P(D|\theta, \mathcal{M}_1)P(\theta|\mathcal{M}_1)\, d\theta} {\int P(D|\theta, \mathcal{M}_0)P(\theta|\mathcal{M}_0)\, d\theta} \end{align} where $D$ is an individual element in the testing set.

However, it turns out that with an uninformative prior, the values of this are the same for all sequences, since, with an uninformative prior, there is not enough information to distinguish between the sequences.

After some confusion, it became clear that one needs an informative prior here, otherwise one cannot do anything. It then occurred to me that one can use the posterior distribution for the parameters instead of the prior, namely $P(\theta|\mathcal{M}_1, D_0)$ in the place of $P(\theta|\mathcal{M}_1)$, in Equation~\ref{bayesfact}. Here $D_0$ corresponds to some other data, that is used to get an informative distribution for the parameters. In this case, the obvious choice is the training set. I wasn't sure if that was actually done or not; however after some searching, I came across this variation described as pseudo-Bayes factors in Ch 5 of Robert's Bayesian Choice''. However, from his description (I have not looked into this otherwise) the main motivation for doing this is that bane of Bayesians, improper priors. I.e. if one conditions on a training sample. the posterior corresponding to an improper prior may be proper. I don't know if anyone is using this in the context I have described.

So, let $K$ be a pseudo-Bayes factor. Then \begin{align*} K = \frac{P(D|\mathcal{M}_1, D_0)} {P(D|\mathcal{M}_0, D_0)} = \frac{\int P(D|\theta, \mathcal{M}_1, D_0)P(\theta|\mathcal{M}_1, D_0)\, d\theta} {\int P(D|\theta, \mathcal{M}_0, D_0)P(\theta|\mathcal{M}_0, D_0)\, d\theta} \end{align*} where $D_0$ is the training set, and $D$ is an individual element in the testing set as before.

This is the same as the ratio of the the post-predictive distributions corresponding to $\mathcal{M}_1$ and $\mathcal{M}_0$ respectively. As usual for Bayes factors, greater than 1 accepts the alternative hypothesis, less than 1 accepts the null.

What I did for testing was permute the testing set, which is an array, along the columns. This destroyed any correlation structure that might be present. Call $D^{\prime}$ the testing set, and $D^{\prime}_p$ the permuted set. Then I took all the sequences from $D^{\prime}$ and $D^{\prime}_p$ (which are two sets of equal size) and calculated the log of the pseudo-Bayes factor for them. One would expect the sequences in $D^{\prime}_p$ to fall to the left of 0 and $D^{\prime}$ to the right of 0. Preliminary testing gave something close to that. I was wondering if this was a reasonable testing framework, and how to summarize the results. The motivation for the permutation was to get a dataset which had a similar (empirical) marginal distribution to the real dataset, but also which had had any possible correlation structure removed from it. The idea was to construct a set that would be as hard to distinguish from the original dataset as possible. It would nice to have a more formal criteria than this, though.

So, here are some questions.

1. Is it reasonable to test on $D^{\prime}$ and $D^{\prime}_p$? Can one give a formal reason for doing this?

2. I thought I would plot the PBF for one instance of $D^{\prime}$ and $D^{\prime}_p$. What would be a good way to graphically represent this? Some kind of bar plot?

3. I thought of summarizing the overall classification results by using a classification table, i.e. a 4-way table with the counts that

• correctly accept $\mathcal{H}_0$
• correctly accept $\mathcal{H}_1$
• incorrectly accept $\mathcal{H}_0$
• incorrectly accept $\mathcal{H}_1$

Does this look reasonable? It is Ok to add total classification counts across all training sets? That seems like a good summary approach.

4. Does this approach using a pseudo-Bayes factor makes sense, and are there any good alternatives to doing it this way?

5. One is expected to show how this method performs in comparison with other motif finding approaches. Do you have specific suggestions about which methods to compare it to?

6. There is something called the 0.632+ Bootstrap, which is supposed to be an improvement on regular cross-validation. Is this something worth considering?

I'm also a little fuzzy about the justification of the Bayes factor, but that can wait for another question.

If you think the notation in this question could be improved, please comment.