\documentclass[11pt]{article}
\usepackage[margin=1.2in,a4paper]{geometry} % try also "letterpaper" instead of "a4paper"
\usepackage{enumitem} % this package is used to generate the custom chapter prefix "[Chapter number].[Exercise Number]"
\usepackage{tikz} % this package is used for drawings
\usetikzlibrary{automata,arrows,shapes,decorations.pathreplacing} % we give some instructions to tikz
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{url}

%%%%%%%%%% macros %%%%%%%%%%
  \newcommand{\Fig}[1]{Figure~\ref{#1}}
  \newcommand{\Thm}[1]{Theorem~{#1}}
  \newcommand{\Sec}[1]{Section~{#1}}   
  \newcommand{\Prob}[1]{Problem~{#1}}  
  \newcommand{\Exe}[1]{Exercise~{#1}}  
  \newcommand{\Algo}[1]{Algorithm~{#1}}
  \newcommand{\ltdots}{..}  
  \newcommand{\wt}{\::\:}
  \newcommand{\st}{\:|\:}
  \newcommand{\Eq}[1]{Equation~({#1})} 
  \newcommand{\Chap}[1]{Chapter~{#1}}
  \newcommand{\Ex}[1]{Example~{#1}}     
  \DeclareMathOperator*{\argmax}{argmax}  
%%%%%%%%%% macros %%%%%%%%%%

\begin{document}

%%%%%%%%%% header %%%%%%%%%%
\noindent {\sf\LARGE GENOME-SCALE ALGORITHM DESIGN} \\
\noindent {\sf by Veli M{\" a}kinen, Djamal Belazzougui, Fabio Cunial and Alexandru I. Tomescu}\\
\noindent {\sf Cambridge University Press, 2nd edition, 2023}\\
\noindent {\sf www.genome-scale.info}\\
\noindent\makebox[\linewidth]{\rule{\paperwidth}{0.4pt}}
%%%%%%%%%% header %%%%%%%%%%

\section*{Exercises for Chapter 7. Hidden Markov models (HMMs)}

\begin{enumerate}[label=7.\arabic*]

\item Argue that $O(|H|)$ working space suffices in \Prob{7.3} to compute the final numerical value (not the actual solution path) with the Viterbi algorithm, and in \Prob{7.4} with the forward and backward algorithms.
\item
\label{exe:spaceefficientviterbi}
To traceback the solution of Viterbi algorithm, the naive solution requires $O(n|H|)$ working space. Show that this can be improved to $O(\sqrt{n}|H|)$ by sampling every $\sqrt{n}$th position in $S$, without asymptotically affecting the running time of the traceback.
\item The standard algorithm for the multiplication of two matrices of $m$ rows by $m$ columns takes time $O(m^3)$.
More sophisticated approaches achieve time $O(m^\omega)$ with $\omega<3$ (the current record is $\omega\leq 2.38$
and $\omega=2$ is the best one can hope for).
Suppose that we have to compute many instances of the forward and backward algorithms on different sequences, but on the same HMM. Show how to use fast matrix multiplication to improve the running time of such a computation.
\item Given an HMM and a sequence $S = s_1\cdots s_n$, derive an $O(n|H|)$ time algorithm to compute, for all $h \in H$,
\begin{equation}
\argmax_i \{\mathbb{P}(P, s_1\cdots s_n) \st P = p_0p_1\cdots p_np_{n+1} \in \mathcal{P}(n) \text{ and } p_i = h\}.
\end{equation}
\item Multiplication is the source of numerical errors in HMM algorithms, when the numbers become too large to fit into a computer word. Show how the Viterbi algorithm can be implemented using a sum of logarithms to avoid these numerical problems.
\label{exe:viterbi-numerical-errors}
\item For the forward and backward algorithms the sum of logarithms conversion is not enough for numerical stability. Browse the literature to find a solution for this problem.
\label{exe:BF-numerical-errors}
\item The flexibility of choosing the states, transitions, emissions, and their probabilities makes
HMMs a powerful modeling device. So far we have used a \emph{zeroth-order Markov model} for emission probabilities
(probabilities only depended on the state, not on the sequence context). We could just as well use 
\emph{first-order Markov chains} or, more generally, \emph{$k$-th order Markov chains}, in which the probability depends on the state and on the last $k$ symbols preceding the current one: $\mathbb{P}(s_i \mid s_{i-k} \cdots s_{i-1})=\mathbb{P}(s_i \mid s_{1} \cdots s_{i-1})$.

Notice that the states of the HMM are independent, in the sense that each state can choose a Markov chain of a different order from that of the Markov chain for its emission probabilities. In addition to the use of different order Markov chains, we could adjust how many symbols are emitted in each state. Use these considerations to design a realistic HMM for \emph{eukaryote gene prediction}. Try to take into account intron/exon boundary di-nucleotides, codon adaptation, and other known features of eukaryote genes. Consider also how you can train the HMM.
\item 
\emph{Profile HMMs} are an extension of HMMs to the problem of aligning a sequence with an existing multiple alignment (profile).
Consider for example a multiple alignment of a protein family:
\begin{center}
\begin{tabular}{c}
\texttt{AVLSLSKTTNNVSPA} \\
\texttt{AV-SLSK-TANVSPA} \\
\texttt{A-LSLSK-TANV-PA} \\
\texttt{A-LSSSK-TNNV-PA} \\
\texttt{AS-SSSK-TNNV-PA} \\
\texttt{AVLSLSKTTANV-PA}
\end{tabular}
\end{center}
We considered the problem of aligning a sequence $A$ against a profile in the context of progressive multiple alignment in \Sec{6.6.4}, and
the idea was to consider the multiple alignment as a sequence of columns and apply normal pair-wise alignment with proper extensions of substitution and indel scores.
Consider $A=\mathtt{AVTLSLSTAANVSPA}$ aligned to our example profile above, for example, as follows:
\begin{center}
\begin{tabular}{c}
\texttt{AVTLSLS--TAANVSPA} \\
\\[-.3cm]
\texttt{AV-LSLSKTTN-NVSPA} \\
\texttt{AV--SLSK-TA-NVSPA} \\
\texttt{A--LSLSK-TA-NV-PA} \\
\texttt{A--LSSSK-TN-NV-PA} \\
\texttt{AS--SSSK-TN-NV-PA} \\
\texttt{AV-LSLSKTTA-NV-PA}
\end{tabular}
\end{center}
Here we have added two gaps to the sequence and two gap columns to the profile following the ``once a gap, always a gap'' principle.

Profile HMMs are created
using \emph{inhomogeneous} Markov chains, such that each of the columns will form separate match, insertion, and deletion states, and transitions go from left to right, as illustrated in \Fig{fig:profileHMM}. Match and deletion states emit the columns of the profile, so they do not contain self-loops. Insertion states emit symbols from the input sequence, so they contain self-loops to allow any number of symbols emitted between states that emit also columns of the profile.

\begin{figure}
\begin{center}
\begin{tikzpicture}[xscale=2.1,yscale=1.7]
\node[circle,fill=none,draw=black,minimum size=12mm] (start) at (0,-1) {\texttt{start}};
\node[diamond,fill=none,draw=black,minimum size=12mm] (I0) at (0,0) {$I^{0}$};
\node[circle,fill=none,draw=black,minimum size=10mm] (D1) at (1,1) {$D^{1}$};
\node[diamond,fill=none,draw=black,minimum size=12mm] (I1) at (1,0) {$I^{1}$};
\node[rectangle,fill=none,draw=black,minimum size=9mm] (M1) at (1,-1) {$M^{1}$};
\node[circle,fill=none,draw=black,minimum size=10mm] (D2) at (2,1) {$D^{2}$};
\node[diamond,fill=none,draw=black,minimum size=12mm] (I2) at (2,0) {$I^{2}$};
\node[rectangle,fill=none,draw=black,minimum size=9mm] (M2) at (2,-1) {$M^{2}$};
\node[circle,fill=none,draw=none,minimum size=12mm] (cdotsplus) at (3,1) {$\cdots$};
\node[circle,fill=none,draw=none,minimum size=12mm] (cdots) at (3,0) {$\cdots$};
\node[circle,fill=none,draw=none,minimum size=12mm] (cdotsminus) at (3,-1) {$\cdots$};
\node[circle,fill=none,draw=black,minimum size=10mm] (Dm) at (4,1) {$D^{m}$};
\node[diamond,fill=none,draw=black,minimum size=12mm] (Im) at (4,0) {$I^{m}$};
\node[rectangle,fill=none,draw=black,minimum size=9mm] (Mm) at (4,-1) {$M^{m}$};
\node[circle,fill=none,draw=black,minimum size=12mm] (end) at (5,-1) {\texttt{end}};
\draw[->,out=105,in=160,loop,looseness=2] (I0) to (I0);
\draw[->,out=105,in=160,loop,looseness=2] (I1) to (I1);
\draw[->,out=105,in=160,loop,looseness=2] (I2) to (I2);
\draw[->,out=105,in=160,loop,looseness=2] (Im) to (Im);
\draw[->] (start) to (I0);
\draw[->] (start) to (D1);
\draw[->] (start) to (M1);
\draw[->] (I0) to (D1);
\draw[->] (I0) to (M1);
\draw[->] (D1) to (I1);
\draw[->] (D1) to (D2);
\draw[->] (D1) to (M2);
\draw[->] (M1) to (I1);
\draw[->] (M1) to (D2);
\draw[->] (M1) to (M2);
\draw[->] (I1) to (D2);
\draw[->] (I1) to (M2);
\draw[dashed] (I2) to (cdotsplus);
\draw[dashed] (I2) to (cdotsminus);
\draw[dashed] (D2) to (cdotsplus);
\draw[dashed] (D2) to (cdotsminus);
\draw[dashed] (M2) to (cdotsplus);
\draw[dashed] (M2) to (cdotsminus);
\draw[->] (D2) to (I2);
\draw[->] (M2) to (I2);
\draw[->,dashed] (cdots) to (Dm);
\draw[->,dashed] (cdots) to (Mm);
\draw[->,dashed] (cdotsplus) to (Dm);
\draw[->,dashed] (cdotsplus) to (Mm);
\draw[->,dashed] (cdotsminus) to (Dm);
\draw[->,dashed] (cdotsminus) to (Mm);
\draw[->] (Dm) to (end);
\draw[->] (Dm) to (Im);
\draw[->] (Im) to (end);
\draw[->] (Mm) to (end);
\draw[->] (Mm) to (Im);
\end{tikzpicture}
\caption{Profile HMM illustration without showing the transition and emission probabilities.\label{fig:profileHMM}}
\end{center}
\end{figure}

Since the resulting HMM is reading only one sequence, the Viterbi, forward, and backward algorithms are almost identical to the ones we have studied so far.
The only difference is that deletion states are \emph{silent} with respect to the input string, since they do not emit any symbol.
\begin{itemize}
\item[a)] Modify the Viterbi recurrences to handle both emitting and silent states.
\item[b)] Derive the Viterbi recurrences specific to profile HMMs.
\end{itemize}

\item Derive a local alignment version of profile HMMs.
\item 
\emph{Pair HMMs} are a variant of HMMs emitting two sequences, such that a path through the HMM can be interpreted as an alignment of the input sequences. Such pair HMMs have a \emph{match} state emitting a symbol from both sequences simultaneously, and symmetric \emph{insertion} and \emph{deletion} states to emit only from one input sequence.
\begin{itemize}
\item[a)] Fix a definition for pair HMMs and derive the corresponding Viterbi, forward, and backward recurrences.
\emph{Hint.} The result should look very similar to Gotoh's algorithm for global alignment with affine gap costs from \Sec{6.4.4}. 
\item[b)] Apply a derivation for pair HMMs similar to the ones we used in obtaining relation ({7.13}), in order to define the probability of $a_i$ aligning to $b_j$ over all alignments of $A=a_1 \cdots a_m$ and $B=b_1 \cdots b_n$.
\item[c)] Let $p_{ij}$ denote the probability derived above to align $a_i$ to $b_j$. We say that the \emph{most robust alignment}
of $A$ and $B$ is the alignment maximizing the sum of values $p_{ij}$ over $i,j$ such that the $a_i \to b_j$ substitution is part of the alignment. Derive a dynamic programming algorithm to compute this most robust alignment.
\end{itemize}

\section*{Additional exercises not in the book}

\item 

Assume a set of DNA sequences with coding/non-coding labeling. 
\emph{Implement} a program that trains the emission/transition probabilities for the coding/non-coding HMM considered at the book, given the training data.

\item 

\emph{Implement} the Viterbi algorithm in the special case of the coding/non-coding HMM.
Implement also the tracing back of the optimal path and deduce the labels for some example sequence.
If you did the previous assignment, train the HMM with some other example sequence. 

\item 

Familiarize yourself with the HMM package of Biopython: \url{http://biopython.org/DIST/docs/api/Bio.HMM-module.html}. Explain how the provided functions match the description in the book.


\item 

Use Biopython to run the Viterbi algorithm on the coding/non-coding HMM.




\end{enumerate}

\end{document}  
