diff --git a/BINF2025_TP4.ipynb b/BINF2025_TP4.ipynb
index be4f132..7bc60da 100644
--- a/BINF2025_TP4.ipynb
+++ b/BINF2025_TP4.ipynb
@@ -1,347 +1,404 @@
{
- "nbformat": 4,
- "nbformat_minor": 0,
- "metadata": {
- "colab": {
- "provenance": [],
- "authorship_tag": "ABX9TyOkaYHMV5kAI2JXjqDKbiVa",
- "include_colab_link": true
- },
- "kernelspec": {
- "name": "python3",
- "display_name": "Python 3"
- },
- "language_info": {
- "name": "python"
- }
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "colab_type": "text",
+ "id": "view-in-github"
+ },
+ "source": [
+ "![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"
+ ]
},
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {
- "id": "view-in-github",
- "colab_type": "text"
- },
- "source": [
- ""
- ]
- },
- {
- "cell_type": "markdown",
- "source": [
- "# BINF TP4 - Algorithmes d'alignement par HMM"
- ],
- "metadata": {
- "id": "lSPnN2vrOzhO"
- }
- },
- {
- "cell_type": "markdown",
- "source": [
- "Dans ce TP nous allons implémenter un algorithme d'alignement par chaine de Markov cachée.\n",
- "\n",
- "## Notation\n",
- "\n",
- "Soit un alignement $A$, sa profondeur $prof(A) = P$ est le nombre de séquences alignées, sa largeur $larg(A) = N$ est la taille (identique) de chacune de ses séquences et $A_{(k,l)}$ est le ième caractère de sa kème séquence avec $0 \\geq k < P$ et $0 \\geq l < N$. On note $P$:Acide Aminé $\\rightarrow$ Entier la fonction qui associe une valeur entière unique dans [0,20] à chaque acide aminé et donnée par :\n",
- "\n",
- "\n"
- ],
- "metadata": {
- "id": "c6NL0YIMPDS3"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "P = {\"A\": 0, \"C\": 1, \"D\": 2, \"E\": 3, \"F\": 4, \"G\": 5, \"H\": 6, \"I\": 7,\n",
- " \"K\": 8, \"L\": 9, \"M\": 10, \"N\": 11, \"P\": 12, \"Q\": 13, \"R\": 14, \"S\": 15,\n",
- " \"T\": 16, \"V\": 17, \"W\": 18, \"Y\": 19}"
- ],
- "metadata": {
- "id": "A8LhiAAeP9u-"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "source": [
- "En plus, voici quelques fonctions de support pour charger des séquences d'un fichier fasta et sauvegarder et charger une chaîne de Markov :"
- ],
- "metadata": {
- "id": "qEi7qNzIQ0fV"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "def read_fasta(fname):\n",
- " #Read a fasta file with multiple entries.\n",
- " #Returns a list of sequences.\n",
- " res = []\n",
- " with open(fname, 'r') as f:\n",
- " tmp = None\n",
- " for line in f:\n",
- " if line.startswith(\">\"):\n",
- " if tmp:\n",
- " res.append(tmp)\n",
- " tmp = \"\"\n",
- " else:\n",
- " tmp += line.rstrip()\n",
- " if tmp:\n",
- " res.append(tmp)\n",
- " return res\n",
- "\n",
- "def save_HMM(T, e_M, e_I, f):\n",
- " #Save the HMM specified by (T, e_M, e_I) into the opened file f\n",
- " f.write(\"{}\\n\".format(T.shape[0]))\n",
- " for i in range(T.shape[0]):\n",
- " f.write(\",\".join([\"{:.3f}\".format(T[i,j]) for j in range(9)]) + \"\\n\")\n",
- " for i in range(e_M.shape[0]):\n",
- " f.write(\",\".join([\"{:.3f}\".format(e_M[i,j]) for j in range(20)]) + \"\\n\")\n",
- " for i in range(e_I.shape[0]):\n",
- " f.write(\",\".join([\"{:.3f}\".format(e_I[i,j]) for j in range(20)]) + \"\\n\")\n",
- "\n",
- "def load_HMM(f):\n",
- " #Load the HMM stored in the opened file f into matrices\n",
- " #Returns (T, e_M, e_I)\n",
- " line = f.readline()\n",
- " N = int(line.rstrip(\"\\n\"))\n",
- " T = np.zeros((N, 9))\n",
- " e_M = np.zeros((N, 20))\n",
- " e_I = np.zeros((N, 20))\n",
- "\n",
- " for i in range(N):\n",
- " T[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
- " for i in range(N):\n",
- " e_M[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
- " for i in range(N):\n",
- " e_I[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
- "\n",
- " return (T, e_M, e_I)"
- ],
- "metadata": {
- "id": "30qlDWTmQ47E"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "source": [
- "## Exemples:\n",
- "Vous trouverez des exemples de fichiers d'entrées et de sorties [ici](https://www.dropbox.com/scl/fo/llpj27nr17qkzrcc2w43f/ADdTZVnoqQYK1Cgrplggdao?rlkey=83glsunr3k3fgxhp4w4nidgxj&st=wqo11d1h&dl=0).\n",
- "\n",
- "Voici un extrait :\n",
- ""
- ],
- "metadata": {
- "id": "OkQ9ze0RXbKk"
- }
- },
- {
- "cell_type": "markdown",
- "source": [
- "# Modèle Plan 9\n",
- "\n",
- "On va utiliser le modèle de HMM Plan 9 vu en cours. Ce modèle est structuré selon les colonnes de l'alignement et possède 3 types d'états différents : $M$ représentant les correspondances et substitutions, $D$ représentant les délétions et $I$ représentant les insertions. Chaque colonne contient ces trois états sauf la première qui n'a pas d'état $D$. Les états $M$ et $I$ émettent des caractères et ont donc des vecteurs d'émissions associés, les états $D$ son silencieux.\n",
- "\n",
- "\n",
- "\n",
- "Une HMM représentant un alignement contenant $N$ colonnes est représenté par trois matrices :\n",
- "\n",
- "1. **La matrice de transition** $T$ de taille $(N+1) \\times 9$ où chaque ligne $T_u$ contient les $9$ probabilités de transition des états de la colonne $u$ $(M_u,I_u,D_u)$ vers leurs successeurs. L'ordre des transitions est : $M_u \\rightarrow M_{(u+1)}, M_u \\rightarrow D_{(u+1)}, M_u \\rightarrow I_u, D_u \\rightarrow M_{(u+1)}, D_u \\rightarrow D_{(u+1)}, D_u \\rightarrow I_u, I_u \\rightarrow M_{(u+1)}, I_u \\rightarrow D_{(u+1)}, I_u \\rightarrow I_u$.\n",
- "\n",
- " La première colonne $u=0$ est spéciale, c'est l'état de départ du modèle, elle ne contient pas d'état $D$ donc les probabilités associées sont fixées à $0$. De la même manière, la dernière colonne ne possède que des transitions vers l'état $M$, les valeurs des autres transitions sont fixées à $0$.\n",
- "\n",
- "\n",
- "\n",
- "2. **Les matrices d'émissions** $e^M$ et $e^I$ de tailles $(N+1) \\times 20$ où $e_{u,P(c)}^M$ et $e_{u,P(c)}^I$ contiennent respectivement la probabilité d'émission du caractère $c$ pour les états $M$ et $I$ de la colonne $u$. L'état $M$ de la première colonne n'émet pas de caractère donc $e_{0,*}^M=0$."
- ],
- "metadata": {
- "id": "CEkLeT2-M4YO"
- }
- },
- {
- "cell_type": "markdown",
- "source": [
- "# Exercice 1 - Construction de la chaîne\n",
- "\n",
- "Construire une HMM revient à remplir les matrices $T$, $e^M$ et $e^I$ à partir d'un alignement $A$. Pour commencer, on va déterminer la taille de la chaîne. Pour cela, on introduit le paramètre $α \\in [0,1]$ qu'on va utiliser pour déterminer les colonnes à marquer dans $A$. Une colonne marquée correspond à une colonne de la chaîne.\n",
- "\n",
- "---\n",
- "\n",
- "Q1. Ecrire la fonction\n",
- "\n",
- "mark_columns(A: list, alpha: float) -> list\n",
- "\n",
- "qui prend en entrée un alignement (une liste de séquences protéiques de même taille) et une valeur de $\\alpha$ et retourne la liste des indices des colonnes marquées. On marque une colonne quand elle contient une proportion de caractères \"-\" (de trou) $< \\alpha$."
- ],
- "metadata": {
- "id": "hEobsQ4yP_un"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "#Votre code ici"
- ],
- "metadata": {
- "id": "ikf5pcUtRYvV"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "source": [
- "---\n",
- "Maintenant qu'on a $N$ colonnes marquées, on va construire la chaîne en tant que tel en suivant les etapes suivantes :\n",
- "\n",
- "1. Construire la matrice $T$, de tailles $(N+1)\\times 9$ et les matrices $e^M$ et $e^I$ de tailles $(N+1) \\times 20$. Initialiser les valeurs de ces matrices à $1$ sauf les cas particuliers mentionnés ci-dessus.\n",
- "\n",
- "2. Pour remplir les matrices, on va parcourir chaque séquence de l’alignement et commencer par déterminer à quelle colonne $k$ et état $π_k$ est associée chaque caractère, puis mettre à jours les matrices en utilisant la procédure suivant:\n",
- "\n",
- "\n",
- "3. Normaliser chaque élément de $T$ de telle sorte à ce que la somme des probabilités de transition pour chaque état de chaque colonne soit égale à $1$. Normaliser $e^M$ et $e^I$ de telle sorte à ce que la somme des probabilités d'émission des caractères pour chaque colonne soit égale à $1$.\n",
- "\n",
- "Remarque : si vous rencontrez un caractère X dans les séquences alors ignorez le pour le remplissage de $e^M$ et $e^I$.\n",
- "\n",
- "---\n",
- "\n",
- "Q2. Ecrire la fonction\n",
- "\n",
- "build_HMM(A: list, mark: list) -> (T, eM, eI)\n",
- "\n",
- "qui retourne les matrices $T$, $eM$ et $eI$ representant la chaine construite a partir de l'alignement $A$ et de la liste de colonnes marques $mark$."
- ],
- "metadata": {
- "id": "WaTWmtVRRaj6"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "#Votre code ici"
- ],
- "metadata": {
- "id": "Y3TXPqNhTwNo"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "source": [
- "# Exercice 2 : Alignement d'une séquence sur le modèle\n",
- "\n",
- "Maintenant qu'on a construit le modèle, on peut l'utiliser pour aligner des séquences sur celui-ci. Pour cela, on va utiliser l'algorithme de Viterbi qui est un algorithme de programmation dynamique (structure similaire à Levenshtein) qui va calculer la séquence d'états cachés $\\Pi$ de probabilité maximale étant donnée la séquence $s$. A partir de ces états, on pourra calculer le score de l'alignement ainsi que la superposition correspondante.\n",
- "\n",
- "\n",
- "Pour une HMM de taille $N$ et une séquence à aligner de taille $M$, on va remplir une matrice de score $V$ et de retour $B$ de taille $(3N+1) \\times (M+1)$ où $3N$ est le nombre d'état de la HMM, + 1 pour l'état spécial END et $M+1$ car la première colonne correspond à la séquence vide (comme pour l'alignement).\n",
- "\n",
- "* Bords : La première colonne et les premières et secondes lignes sont fixées à $-\\infty$ sauf la case haut-gauche qui vaut $0$ : $v_{(i,0)}= -\\infty$ pour $i=1…(3N+1)$, $v_{(0,j)}= -\\infty$ pour $j=1…(M+1)$, $v_{(1,j)}= - \\infty$ pour $j=0…(M+1)$ et $v_{0,0}=0$. Notez que $log(1)=0$, donc en donnant cette valeur on indique une probabilité $1$ de démarrer dans l'état de départ. Pour $B$, la première colonne et les premières et deuxièmes lignes sont initialisées à la valeur vide : $b_{(i,0)}=\\emptyset$ pour $i=0…(3N+1)$ et $b_{(0,j)}=b_{(i,j)}=\\emptyset$ pour $j=0…(M+1)$.\n",
- "\n",
- "\n",
- "\n",
- "* Remplissage de V : on va utiliser les équations suivantes pour remplir la case à la position $2 \\leq i < 3N$ (indice de la HMM) et $0< j < M$ (indice de la séquence) :\n",
- "\n",
- "1. Si on est sur un état $M$ ($i$ multiple de $3$) et $i>2$, alors\n",
- "\n",
- "$$\n",
- "\\begin{equation*}\n",
- "v_{(i,j)}=log(e_{(i//3,P(s_{(j-1)}))}^M + \\epsilon) + max\n",
- "\\begin{cases}v_{(i-3,j-1)}+log(t_{((i-3)//3,M \\rightarrow M)})\\\\v_{(i-2,j-1)}+log(t_{((i-2)//3,D \\rightarrow M)})\\\\ v_{(i-1,j-1)}+log(t_{((i-1)//3,I \\rightarrow M))}) \\end{cases},\n",
- "\\end{equation*}\n",
- "$$\n",
- "\n",
- " Où $\\epsilon > 0$ est un paramètre de l'algorithme et $a//b$ est la division entière de $a$ par $b$.\n",
- "\n",
- " 2. Si on est sur un état $D$ ($i-1$ multiple de $3$) et $i > 2$, alors\n",
- "\n",
- "$$\n",
- "\\begin{equation*}\n",
- "v_{(i,j)} = max \\begin{cases} v_{(i-4,j)} + log(t_{((i-4)//3,M \\rightarrow D)})\\\\ v_{(i-3,j)} + log(t_{((i-3)//3,D \\rightarrow D)})\\\\ v_{(i-2,j)} + log(t_{((i-2)//3,I \\rightarrow D))}) \\end{cases}\n",
- "\\end{equation*}\n",
- "$$\n",
- "\n",
- " 3. Si on est sur un état $I$ ($i - 2$ multiple de $3$) alors\n",
- "\n",
- "$$\n",
- "\\begin{equation*}\n",
- "v_{(i,j)} = log(e_{(i//3,P(s_{(j-1)}))}^I + \\epsilon) + max \\begin{cases} v_{(i-2,j-1)} + log(t_{((i-2)//3,M \\rightarrow I)})\\\\ v_{(i-1,j-1)} + log(t_{((i-1)//3,D \\rightarrow I)})\\\\ v_{(i,j-1)} + log(t_{(i//3,I \\rightarrow I)})\\end{cases}\n",
- "\\end{equation*}\n",
- "$$\n",
- "\n",
- "4. La dernière ligne $i=3N$ est un cas particulier (c'est l'état de fin), on la remplie comme un état $M$ qui n'émet aucun caractère (donc on néglige la première partie de l'équation avec le log) et sans changer de position dans la séquence (utiliser j au lieu de j-1) :\n",
- "\n",
- "$$\n",
- "\\begin{equation*}\n",
- "v_{(3N,j)}= max \\begin{cases} v_{(i-3,j)} + log(t_{((i-3)//3,M \\rightarrow M)})\\\\ v_{(i-2,j)} + log(t_{((i-2)//3,D \\rightarrow M)})\\\\ v_{(i-1,j)} + log(t_{((i-1)//3,I \\rightarrow M)}) \\end{cases}\n",
- "\\end{equation*}\n",
- "$$\n",
- "\n",
- "En cas d'égalité dans les scores, suit l'ordre : $M > D > I$.\n",
- "\n",
- "Le score d'alignement est donnée dans la case base-droite de $V$.\n",
- "\n",
- "* Remplissage de $B$ : en même temps que l'on remplit la matrice $V$ on va aussi remplir la matrice de retour $B$ tel que $b_{(i,j)} = (k,l)$ les coordonnées de la case utilisée pour remplir $v_{(i,j)}$.\n",
- "\n",
- "---\n",
- "Q1. Ecrire la fonction\n",
- "\n",
- "viterbi_fwd(T: array, eM: array, eI: array, s: str) -> (array, array)\n",
- "\n",
- "qui retourne les matrices $V$ et $B$ a partir de la chaine (T, eM, eI) et de la sequence $s$."
- ],
- "metadata": {
- "id": "Zr8PUqENTxcd"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "#Votre code ici"
- ],
- "metadata": {
- "id": "2AxydaWBbPe9"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "markdown",
- "source": [
- "\n",
- "**Etape de retour** : L'étape de retour permet de reconstituer la séquence des états cachés $\\Pi$ à partir de laquelle on obtient la superposition des séquences.\n",
- "\n",
- "On va démarrer de la case bas-droite de $B$ et se déplacer en suivant les coordonnées des cases de $B$ jusqu'à atteindre une valeur nulle $\\emptyset$. Pour chaque nouvel état visité sauf le premier, on ajoute son type au début de la séquence $\\Pi$ et on met à jour la superposition :\n",
- "\n",
- "* Si c'est un état $M$ ou $I$, on écrit $s_{(j-1)}$ (la première colonne de $V$ et $B$ correspond à la séquence vide).\n",
- "* Si c'est un état $D$ on écrit « - ».\n",
- "\n",
- "---\n",
- "Q2. Ecrire la fonction\n",
- "\n",
- "viterbi_bwd(V: array, B: array, s: str) -> (float, str, str)\n",
- "\n",
- "qui retourne les matrices le score d'alignement, la séquence d'états et la séquence $s$ alignée etant donnée les matrices $V$ et $B$.\n",
- "\n"
- ],
- "metadata": {
- "id": "F-3Vn80UbQsQ"
- }
- },
- {
- "cell_type": "code",
- "source": [
- "#Votre code ici"
- ],
- "metadata": {
- "id": "48vtEZxoAkeE"
- },
- "execution_count": null,
- "outputs": []
- }
- ]
-}
\ No newline at end of file
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "lSPnN2vrOzhO"
+ },
+ "source": [
+ "# BINF TP4 - Algorithmes d'alignement par HMM"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "c6NL0YIMPDS3"
+ },
+ "source": [
+ "Dans ce TP nous allons implémenter un algorithme d'alignement par chaine de Markov cachée.\n",
+ "\n",
+ "## Notation\n",
+ "\n",
+ "Soit un alignement $A$, sa profondeur $prof(A) = P$ est le nombre de séquences alignées, sa largeur $larg(A) = N$ est la taille (identique) de chacune de ses séquences et $A_{(k,l)}$ est le ième caractère de sa kème séquence avec $0 \\geq k < P$ et $0 \\geq l < N$. On note $P$:Acide Aminé $\\rightarrow$ Entier la fonction qui associe une valeur entière unique dans [0,20] à chaque acide aminé et donnée par :\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "id": "A8LhiAAeP9u-"
+ },
+ "outputs": [],
+ "source": [
+ "P = {\"A\": 0, \"C\": 1, \"D\": 2, \"E\": 3, \"F\": 4, \"G\": 5, \"H\": 6, \"I\": 7,\n",
+ " \"K\": 8, \"L\": 9, \"M\": 10, \"N\": 11, \"P\": 12, \"Q\": 13, \"R\": 14, \"S\": 15,\n",
+ " \"T\": 16, \"V\": 17, \"W\": 18, \"Y\": 19}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "qEi7qNzIQ0fV"
+ },
+ "source": [
+ "En plus, voici quelques fonctions de support pour charger des séquences d'un fichier fasta et sauvegarder et charger une chaîne de Markov :"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "id": "30qlDWTmQ47E"
+ },
+ "outputs": [],
+ "source": [
+ "def read_fasta(fname):\n",
+ " #Read a fasta file with multiple entries.\n",
+ " #Returns a list of sequences.\n",
+ " res = []\n",
+ " with open(fname, 'r') as f:\n",
+ " tmp = None\n",
+ " for line in f:\n",
+ " if line.startswith(\">\"):\n",
+ " if tmp:\n",
+ " res.append(tmp)\n",
+ " tmp = \"\"\n",
+ " else:\n",
+ " tmp += line.rstrip()\n",
+ " if tmp:\n",
+ " res.append(tmp)\n",
+ " return res\n",
+ "\n",
+ "def save_HMM(T, e_M, e_I, f):\n",
+ " #Save the HMM specified by (T, e_M, e_I) into the opened file f\n",
+ " f.write(\"{}\\n\".format(T.shape[0]))\n",
+ " for i in range(T.shape[0]):\n",
+ " f.write(\",\".join([\"{:.3f}\".format(T[i,j]) for j in range(9)]) + \"\\n\")\n",
+ " for i in range(e_M.shape[0]):\n",
+ " f.write(\",\".join([\"{:.3f}\".format(e_M[i,j]) for j in range(20)]) + \"\\n\")\n",
+ " for i in range(e_I.shape[0]):\n",
+ " f.write(\",\".join([\"{:.3f}\".format(e_I[i,j]) for j in range(20)]) + \"\\n\")\n",
+ "\n",
+ "def load_HMM(f):\n",
+ " #Load the HMM stored in the opened file f into matrices\n",
+ " #Returns (T, e_M, e_I)\n",
+ " line = f.readline()\n",
+ " N = int(line.rstrip(\"\\n\"))\n",
+ " T = np.zeros((N, 9))\n",
+ " e_M = np.zeros((N, 20))\n",
+ " e_I = np.zeros((N, 20))\n",
+ "\n",
+ " for i in range(N):\n",
+ " T[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
+ " for i in range(N):\n",
+ " e_M[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
+ " for i in range(N):\n",
+ " e_I[i,:] = np.array([float(e) for e in f.readline().rstrip().split(\",\")])\n",
+ "\n",
+ " return (T, e_M, e_I)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "OkQ9ze0RXbKk"
+ },
+ "source": [
+ "## Exemples:\n",
+ "Vous trouverez des exemples de fichiers d'entrées et de sorties [ici](https://www.dropbox.com/scl/fo/llpj27nr17qkzrcc2w43f/ADdTZVnoqQYK1Cgrplggdao?rlkey=83glsunr3k3fgxhp4w4nidgxj&st=wqo11d1h&dl=0).\n",
+ "\n",
+ "Voici un extrait :\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "CEkLeT2-M4YO"
+ },
+ "source": [
+ "# Modèle Plan 9\n",
+ "\n",
+ "On va utiliser le modèle de HMM Plan 9 vu en cours. Ce modèle est structuré selon les colonnes de l'alignement et possède 3 types d'états différents : $M$ représentant les correspondances et substitutions, $D$ représentant les délétions et $I$ représentant les insertions. Chaque colonne contient ces trois états sauf la première qui n'a pas d'état $D$. Les états $M$ et $I$ émettent des caractères et ont donc des vecteurs d'émissions associés, les états $D$ son silencieux.\n",
+ "\n",
+ "\n",
+ "\n",
+ "Une HMM représentant un alignement contenant $N$ colonnes est représenté par trois matrices :\n",
+ "\n",
+ "1. **La matrice de transition** $T$ de taille $(N+1) \\times 9$ où chaque ligne $T_u$ contient les $9$ probabilités de transition des états de la colonne $u$ $(M_u,I_u,D_u)$ vers leurs successeurs. L'ordre des transitions est : $M_u \\rightarrow M_{(u+1)}, M_u \\rightarrow D_{(u+1)}, M_u \\rightarrow I_u, D_u \\rightarrow M_{(u+1)}, D_u \\rightarrow D_{(u+1)}, D_u \\rightarrow I_u, I_u \\rightarrow M_{(u+1)}, I_u \\rightarrow D_{(u+1)}, I_u \\rightarrow I_u$.\n",
+ "\n",
+ " La première colonne $u=0$ est spéciale, c'est l'état de départ du modèle, elle ne contient pas d'état $D$ donc les probabilités associées sont fixées à $0$. De la même manière, la dernière colonne ne possède que des transitions vers l'état $M$, les valeurs des autres transitions sont fixées à $0$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "2. **Les matrices d'émissions** $e^M$ et $e^I$ de tailles $(N+1) \\times 20$ où $e_{u,P(c)}^M$ et $e_{u,P(c)}^I$ contiennent respectivement la probabilité d'émission du caractère $c$ pour les états $M$ et $I$ de la colonne $u$. L'état $M$ de la première colonne n'émet pas de caractère donc $e_{0,*}^M=0$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "hEobsQ4yP_un"
+ },
+ "source": [
+ "# Exercice 1 - Construction de la chaîne\n",
+ "\n",
+ "Construire une HMM revient à remplir les matrices $T$, $e^M$ et $e^I$ à partir d'un alignement $A$. Pour commencer, on va déterminer la taille de la chaîne. Pour cela, on introduit le paramètre $α \\in [0,1]$ qu'on va utiliser pour déterminer les colonnes à marquer dans $A$. Une colonne marquée correspond à une colonne de la chaîne.\n",
+ "\n",
+ "---\n",
+ "\n",
+ "Q1. Ecrire la fonction\n",
+ "\n",
+ "mark_columns(A: list, alpha: float) -> list\n",
+ "\n",
+ "qui prend en entrée un alignement (une liste de séquences protéiques de même taille) et une valeur de $\\alpha$ et retourne la liste des indices des colonnes marquées. On marque une colonne quand elle contient une proportion de caractères \"-\" (de trou) $< \\alpha$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "ikf5pcUtRYvV"
+ },
+ "outputs": [],
+ "source": [
+ "#Votre code ici\n",
+ "def mark_columns(A, alpha):\n",
+ " if A == []:\n",
+ " return []\n",
+ " res=[]\n",
+ " for i in range(len(A[0])):\n",
+ " counter = 0\n",
+ " for seq in A:\n",
+ " if seq[i]== '-':\n",
+ " counter=+1\n",
+ " if counter / len(A) < alpha:\n",
+ " res.append(i)\n",
+ " return res"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "WaTWmtVRRaj6"
+ },
+ "source": [
+ "---\n",
+ "Maintenant qu'on a $N$ colonnes marquées, on va construire la chaîne en tant que tel en suivant les etapes suivantes :\n",
+ "\n",
+ "1. Construire la matrice $T$, de tailles $(N+1)\\times 9$ et les matrices $e^M$ et $e^I$ de tailles $(N+1) \\times 20$. Initialiser les valeurs de ces matrices à $1$ sauf les cas particuliers mentionnés ci-dessus.\n",
+ "\n",
+ "2. Pour remplir les matrices, on va parcourir chaque séquence de l’alignement et commencer par déterminer à quelle colonne $k$ et état $π_k$ est associée chaque caractère, puis mettre à jours les matrices en utilisant la procédure suivant:\n",
+ "\n",
+ "\n",
+ "3. Normaliser chaque élément de $T$ de telle sorte à ce que la somme des probabilités de transition pour chaque état de chaque colonne soit égale à $1$. Normaliser $e^M$ et $e^I$ de telle sorte à ce que la somme des probabilités d'émission des caractères pour chaque colonne soit égale à $1$.\n",
+ "\n",
+ "Remarque : si vous rencontrez un caractère X dans les séquences alors ignorez le pour le remplissage de $e^M$ et $e^I$.\n",
+ "\n",
+ "---\n",
+ "\n",
+ "Q2. Ecrire la fonction\n",
+ "\n",
+ "build_HMM(A: list, mark: list) -> (T, eM, eI)\n",
+ "\n",
+ "qui retourne les matrices $T$, $eM$ et $eI$ representant la chaine construite a partir de l'alignement $A$ et de la liste de colonnes marques $mark$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "Y3TXPqNhTwNo"
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "def build_HMM(A, mark):\n",
+ "\"\"\"(A: list, mark: list) -> (T, eM, eI)\"\"\"\n",
+ "#init matrices\n",
+ " T = np.full((1 + len(mark),9),1)\n",
+ " eM = np.full((1+len(mark),20),1)\n",
+ " eI = np.full((1+len(mark),20),1)\n",
+ " eM[:,0] = 0\n",
+ " eI[:,0] = 0\n",
+ " \n",
+ " \n",
+ " \n",
+ " for k in range(1, len(A[0]+1)): #1, prof de A \n",
+ " pi = [0]\n",
+ " for l in range(1, len[A]+1): #1, larg de A\n",
+ " if l in mark:\n",
+ " if A[k][l] == '-':\n",
+ " pi.append('D')\n",
+ " else:\n",
+ " pi.append('M')\n",
+ " else:\n",
+ " if A[k][l] == '-':\n",
+ " pi.append(0)\n",
+ " else:\n",
+ " pi.append('I')\n",
+ " \n",
+ " l0= 1\n",
+ " while pi[l0] == 0:\n",
+ " l0+=1\n",
+ " \n",
+ " if A[k][l0] == 'I':\n",
+ " eI[0][P(A[k][l0])]\n",
+ " \n",
+ " u=0\n",
+ " for l in range()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "Zr8PUqENTxcd"
+ },
+ "source": [
+ "# Exercice 2 : Alignement d'une séquence sur le modèle\n",
+ "\n",
+ "Maintenant qu'on a construit le modèle, on peut l'utiliser pour aligner des séquences sur celui-ci. Pour cela, on va utiliser l'algorithme de Viterbi qui est un algorithme de programmation dynamique (structure similaire à Levenshtein) qui va calculer la séquence d'états cachés $\\Pi$ de probabilité maximale étant donnée la séquence $s$. A partir de ces états, on pourra calculer le score de l'alignement ainsi que la superposition correspondante.\n",
+ "\n",
+ "\n",
+ "Pour une HMM de taille $N$ et une séquence à aligner de taille $M$, on va remplir une matrice de score $V$ et de retour $B$ de taille $(3N+1) \\times (M+1)$ où $3N$ est le nombre d'état de la HMM, + 1 pour l'état spécial END et $M+1$ car la première colonne correspond à la séquence vide (comme pour l'alignement).\n",
+ "\n",
+ "* Bords : La première colonne et les premières et secondes lignes sont fixées à $-\\infty$ sauf la case haut-gauche qui vaut $0$ : $v_{(i,0)}= -\\infty$ pour $i=1…(3N+1)$, $v_{(0,j)}= -\\infty$ pour $j=1…(M+1)$, $v_{(1,j)}= - \\infty$ pour $j=0…(M+1)$ et $v_{0,0}=0$. Notez que $log(1)=0$, donc en donnant cette valeur on indique une probabilité $1$ de démarrer dans l'état de départ. Pour $B$, la première colonne et les premières et deuxièmes lignes sont initialisées à la valeur vide : $b_{(i,0)}=\\emptyset$ pour $i=0…(3N+1)$ et $b_{(0,j)}=b_{(i,j)}=\\emptyset$ pour $j=0…(M+1)$.\n",
+ "\n",
+ "\n",
+ "\n",
+ "* Remplissage de V : on va utiliser les équations suivantes pour remplir la case à la position $2 \\leq i < 3N$ (indice de la HMM) et $0< j < M$ (indice de la séquence) :\n",
+ "\n",
+ "1. Si on est sur un état $M$ ($i$ multiple de $3$) et $i>2$, alors\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation*}\n",
+ "v_{(i,j)}=log(e_{(i//3,P(s_{(j-1)}))}^M + \\epsilon) + max\n",
+ "\\begin{cases}v_{(i-3,j-1)}+log(t_{((i-3)//3,M \\rightarrow M)})\\\\v_{(i-2,j-1)}+log(t_{((i-2)//3,D \\rightarrow M)})\\\\ v_{(i-1,j-1)}+log(t_{((i-1)//3,I \\rightarrow M))}) \\end{cases},\n",
+ "\\end{equation*}\n",
+ "$$\n",
+ "\n",
+ " Où $\\epsilon > 0$ est un paramètre de l'algorithme et $a//b$ est la division entière de $a$ par $b$.\n",
+ "\n",
+ " 2. Si on est sur un état $D$ ($i-1$ multiple de $3$) et $i > 2$, alors\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation*}\n",
+ "v_{(i,j)} = max \\begin{cases} v_{(i-4,j)} + log(t_{((i-4)//3,M \\rightarrow D)})\\\\ v_{(i-3,j)} + log(t_{((i-3)//3,D \\rightarrow D)})\\\\ v_{(i-2,j)} + log(t_{((i-2)//3,I \\rightarrow D))}) \\end{cases}\n",
+ "\\end{equation*}\n",
+ "$$\n",
+ "\n",
+ " 3. Si on est sur un état $I$ ($i - 2$ multiple de $3$) alors\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation*}\n",
+ "v_{(i,j)} = log(e_{(i//3,P(s_{(j-1)}))}^I + \\epsilon) + max \\begin{cases} v_{(i-2,j-1)} + log(t_{((i-2)//3,M \\rightarrow I)})\\\\ v_{(i-1,j-1)} + log(t_{((i-1)//3,D \\rightarrow I)})\\\\ v_{(i,j-1)} + log(t_{(i//3,I \\rightarrow I)})\\end{cases}\n",
+ "\\end{equation*}\n",
+ "$$\n",
+ "\n",
+ "4. La dernière ligne $i=3N$ est un cas particulier (c'est l'état de fin), on la remplie comme un état $M$ qui n'émet aucun caractère (donc on néglige la première partie de l'équation avec le log) et sans changer de position dans la séquence (utiliser j au lieu de j-1) :\n",
+ "\n",
+ "$$\n",
+ "\\begin{equation*}\n",
+ "v_{(3N,j)}= max \\begin{cases} v_{(i-3,j)} + log(t_{((i-3)//3,M \\rightarrow M)})\\\\ v_{(i-2,j)} + log(t_{((i-2)//3,D \\rightarrow M)})\\\\ v_{(i-1,j)} + log(t_{((i-1)//3,I \\rightarrow M)}) \\end{cases}\n",
+ "\\end{equation*}\n",
+ "$$\n",
+ "\n",
+ "En cas d'égalité dans les scores, suit l'ordre : $M > D > I$.\n",
+ "\n",
+ "Le score d'alignement est donnée dans la case base-droite de $V$.\n",
+ "\n",
+ "* Remplissage de $B$ : en même temps que l'on remplit la matrice $V$ on va aussi remplir la matrice de retour $B$ tel que $b_{(i,j)} = (k,l)$ les coordonnées de la case utilisée pour remplir $v_{(i,j)}$.\n",
+ "\n",
+ "---\n",
+ "Q1. Ecrire la fonction\n",
+ "\n",
+ "viterbi_fwd(T: array, eM: array, eI: array, s: str) -> (array, array)\n",
+ "\n",
+ "qui retourne les matrices $V$ et $B$ a partir de la chaine (T, eM, eI) et de la sequence $s$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "2AxydaWBbPe9"
+ },
+ "outputs": [],
+ "source": [
+ "#Votre code ici"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "F-3Vn80UbQsQ"
+ },
+ "source": [
+ "\n",
+ "**Etape de retour** : L'étape de retour permet de reconstituer la séquence des états cachés $\\Pi$ à partir de laquelle on obtient la superposition des séquences.\n",
+ "\n",
+ "On va démarrer de la case bas-droite de $B$ et se déplacer en suivant les coordonnées des cases de $B$ jusqu'à atteindre une valeur nulle $\\emptyset$. Pour chaque nouvel état visité sauf le premier, on ajoute son type au début de la séquence $\\Pi$ et on met à jour la superposition :\n",
+ "\n",
+ "* Si c'est un état $M$ ou $I$, on écrit $s_{(j-1)}$ (la première colonne de $V$ et $B$ correspond à la séquence vide).\n",
+ "* Si c'est un état $D$ on écrit « - ».\n",
+ "\n",
+ "---\n",
+ "Q2. Ecrire la fonction\n",
+ "\n",
+ "viterbi_bwd(V: array, B: array, s: str) -> (float, str, str)\n",
+ "\n",
+ "qui retourne les matrices le score d'alignement, la séquence d'états et la séquence $s$ alignée etant donnée les matrices $V$ et $B$.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "48vtEZxoAkeE"
+ },
+ "outputs": [],
+ "source": [
+ "#Votre code ici"
+ ]
+ }
+ ],
+ "metadata": {
+ "colab": {
+ "authorship_tag": "ABX9TyOkaYHMV5kAI2JXjqDKbiVa",
+ "include_colab_link": true,
+ "provenance": []
+ },
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.11.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}