\documentclass{article}
\begin{document}


\section{Fonction de score} \label{sec:score}

Soit $A$ un alphabet compl\'et\'e par la lettre "-" qui repr\'esente
un ``gap'' dans une s\'equence, par exemple pour de l'ADN:
\[ A=\{ "A", "C", "G", "T", "-" \} \]
On se donne une fonction score $\sigma (\alpha ,\beta )$ \`a
deux variables $\alpha $ et $\beta $ qui sont des lettres de
l'alphabet $A$. Cette fonction indique le score qu'on souhaite
attribuer \`a l'alignement du caract\`ere $\alpha $ d'une chaine
avec le caract\`ere $\beta $ \`a la position correspondante
de l'autre chaine.

En C, on peut utiliser deux instructions \verb|switch| imbrique\'es
pour impl\'ementer la fonction score (par exemple pour l'ADN),
ou (pour les prot\'eines par exemple) utiliser une matrice
\`a deux param\`etres (les index peuvent etre par exemple le code ASCII
de la lettre). On affecte un score positif lorsque les 2 caract\`eres
sont \'egaux et un score n\'egatif sinon.

\'Ecrire un exemple de fonction \verb|int score(char a,char b)|.

\section{Algorithme d'alignement global.} \label{sec:global}
On se donne deux chaines de caract\`eres \`a aligner le mieux possible.
On autorise l'insertion de gaps. Par exemple pour aligner
\verb|ACTTGTGA| et \verb|CTGTGAAA| on pourrait utiliser:\\
\verb|ACTTGTG--A|\\
\verb|-CT-GTGAAA|\\
Si on se donne un score de +2 pour une correspondance exacte et -1 sinon,
on obtiendrait ici un score total de:\\
-1+2+2-1+2+2+2-1-1+2

Notations: si $S$ est une chaines de caract\`eres, on note $|S|$
la longueur de $S$ et $S[i]$ le $i$-\`eme caract\`ere de la chaine $S$.

Hypoth\`ese: on suppose que le score de l'alignement de deux gaps "-" est
n\'egatif. Ceci assure qu'on n'ins\`ere un nombre fini de "-" pour
obtenir le score maximum, le probleme de maximisation est alors une
recherche de maximum dans un nombre fini de possibilit\'es et admet
donc une solution (pas forc\'ement unique).

Si $S$ et $T$ deux chaines de longueur $n$ et $m$, on note $V(n,m)$
le score maximum d'alignement possible pour $S$ et $T$. On va
calculer $V(n,m)$ par un proc\'ed\'e de r\'ecurrence \`a deux
indices qui ressemble dans son concept au calcul des coefficients
du triangle de Pascal.

On va donc calculer tous les $V(i,j)$ pour $i\leq n$ et
$j \leq m$. Supposons qu'on souhaite calculer $V(i,j)$ avec
$i$ et $j$ non nuls. On regarde quelle peut \^etre la derni\`ere
paire de caract\`eres align\'es. Vu l'hypoth\`ese faire sur la fonction
score, il y a 3 possibilit\'es:
\begin{itemize}
\item on aligne $S[i]$ avec $T[j]$, si c'est le cas:
\[ V(i,j)=V(i-1,j-1)+\sigma (S[i],T[j]) \]
\item on aligne $S[i]$ avec "-", si c'est le cas:
\[ V(i,j)=V(i-1,j)+\sigma (S[i],-) \]
\item on aligne "-" avec $T[j]$, si c'est le cas:
\[ V(i,j)=V(i,j-1)+\sigma (-,S[j]) \]
\end{itemize}
On en d\'eduit donc que pour $i>0$ et $j>0$:
\begin{equation} \label{eq:dynamic}
 V(i,j)=\mbox{max} ( V(i-1,j-1)+\sigma (S[i],T[j]),
V(i-1,j)+\sigma (S[i],-),V(i,j-1)+\sigma (-,T[j]) ) 
\end{equation}

Il suffit ensuite d'initialiser $V(i,0)$ et $V(0,j)$ pour $i\leq n$
et $j\leq m$ ce qui est relativement simple puisque $V(i,0)$ ou
$V(j,0)$ est l'alignement de $i$ lettres avec "-":
\[ V(i,0)=\sum _{k=1}^i \sigma (S[k],-), \quad
V(0,j)=\sum _{k=1}^j \sigma (-,T[k]) \]

Par exemple si on se donne S=ACTG et T=AACC, on obtient (avec la fonction
score valant 2 et -1):\\
\begin{tabular}{rr|r|r|r|r|r|}\\
  & j & 0  & 1  & 2  & 3  & 4 \\
i &   &    & A  & A  & C  & C \\ \hline
0 &   & 0  & -1 & -2 & -3 & -4 \\ \hline
1 & A & -1 & 2  & 1  & 0  & -1 \\ \hline
2 & C & -2 & 1  & 1  & 3  & 2  \\ \hline
3 & T & -3 & 0  & 0  & 2  & 2  \\ \hline
4 & G & -4 & -1 & -1 & 1  & 1  \\ \hline
\end{tabular}\\
On remplit la base ($i=0$ et $j=0$) puis par exemple ligne par
ligne en appliquant la formule (\ref{eq:dynamic}).

L'alignement optimal a donc pour score 1.

%En C, on utilisera un tableau d'entiers a deux indices:\\
%\verb|int V[n+1][m+1];|\\
%pour les scores.
En C++, on utilisera un tableau d'entiers \`a deux indices, qu'on peut
cr\'eer de la mani\`ere suivante:\\
\verb|vector<int> ligne(m+1,0);|\\
\verb|vector< vector<int> > V(n+1,ligne);|\\
On initialise ensuite \verb|V[0][j]| et \verb|V[j][0]| (faire une boucle
\`a chaque fois).
On ex\'ecute alors une double boucle (l'une imbriqu\'ee dans
l'autre), par exemple on fait varier \verb|i| de 1 \`a $n$
puis \`a \verb|i| fix\'e, \verb|j| de 1 \`a $m$ et on applique
la formule de r\'ecurrence (\ref{eq:dynamic})

Am\'elioration: on remarque que la formule (\ref{eq:dynamic}) ne
fait intervenir que la ligne pr\'ec\'edente ($i-1$) 
et la ligne actuelle ($i$), on peut donc se contenter de stocker
seulement les 2 lignes \verb|Vprecedent[m+1]| et \verb|Vactuel[m+1]|

On peut ensuite chercher \`a retracer les alignements de score 1.
Pour cela on part de $V(n,m)$ et on retrace dans (\ref{eq:dynamic})
quel(s) terme(s) ont donn\'e lieu au max. Cela revient \`a regarder
dans les cases au-dessus, \`a gauche, ou en diagonale \`a gauche au-dessus
quel(s) est (sont) la (les) case(s) telles que
la diff\'erence entre les 2 scores correspond bien
au score de l'alignement des 2 lettres (ou gap) qu'on rajoute. Par exemple
ici, on part de $(4,4)$, on peut remonter en $(3,3)$ ou $(3,4)$
puis on se retouve en $(2,3)$ dans les deux cas. Ensuite on ne
peut passer qu'en $(1,2)$ m\^eme si $(2,2)$ a la m\^eme valeur que
$(1,2)$ car le score de $(2,3)$ est plus grand que les pr\'ec\'edents
il y a donc eu alignement de 2 lettres. Depuis $(1,2)$ on a \`a
nouveau 2 possibilit\'es $(1,1)$ et $(0,1)$ et on termine en $(0,0)$
ce qui donne 4 possibilit\'es d'alignement de score 1:\\
\verb|A-CTG  -ACTG  A-CTG  -ACTG|\\
\verb|AACC-  AACC-  AAC-C  AAC-C|

Pour programmer la recherche des chemins, on note dans un tableau au fur 
et \`a mesure du calcul de \verb|V[i][j]|
les directions qui ont donn\'e lieu au \verb|max| de la formule. 
On peut pour cela utiliser un type structur\'e, par exemple:
\begin{verbatim}
struct element {
  int valeur;
  bool en_haut;
  bool a_gauche;
  bool diag;
};

vector<element> ligne(m+1);
vector< vector<element> > V(n+1,ligne);
\end{verbatim}
Si le score \verb|V[i][j]| est obtenu en utilisant la diagonale, on
\'ecrira par exemple:\\
\verb|V[i][j].diag=true;|\\
On cr\'ee ensuite une fonction r\'ecursive \verb|calcule_chemin| dont le r\^ole
est d'ajouter \`a une liste d'alignements les alignements dont on connait
la fin, plus pr\'ecis\'ement qui se 
terminent \`a la position (\verb|pos_a|,\verb|pos_b|)
 par les chaines de caract\`eres \verb|fin_a| et \verb|fin_b|:
\begin{verbatim}
void calcule_chemin(const string & a,const string & b, 
		    const vector< vector<element> > & v,
		    string fin_a, string fin_b, int pos_a, int pos_b,
		    vector<string> & alignements){
  ...
}
\end{verbatim}
Dans \verb|calcule_chemin|, on teste d'abord si \verb|pos_a| ou \verb|pos_b| 
est nul. Par exemple si \verb|pos_a| est nul,
l'alignement est obtenu en rajoutant \verb|pos_b| carcat\`eres \verb|-|
\`a \verb|fin_a| et en rajoutant le d\'ebut (les \verb|pos_b| premiers
caract\`eres) de \verb|b| \`a \verb|fin_b|.

Si \verb|pos_a| et \verb|pos_b| sont non nuls, on teste les trois 
possibilit\'es (par en haut, par la gauche, par la diagonale)
et on appelle r\'ecursivement la fonction d'affichage d'alignement
par exemple pour un d\'eplacement vertical:
\begin{verbatim}
  if (v[pos_a][pos_b].en_haut)
    calcule_chemin(a,b,v,a[pos_a-1]+fin_a,'-'+fin_b,pos_a-1,pos_b,alignements);
\end{verbatim}

\section{Alignement local.} \label{sec:local}
On n'impose plus un alignement de toutes les lettres de chaque chaine
mais seulement d'un morceau de chaque chaine. La m\'ethode de recherche
de l'alignement local optimal est la m\^eme \`a ceci pr\`es que:
\begin{itemize}
\item quand on calcule le max dans l'\'equation (\ref{eq:dynamic}),
on doit rajouter 0, ce qui correspond \`a rejeter un alignement
de valeur n\'egative en pr\'ef\'erant un alignement local avec 0
lettres
\item l'initialisation se fait avec $V(i,0)=V(0,j)=0$
\item le score de l'alignement local doit \^etre recherc\'e dans tout
le tableau $V(i,j)$ et pas seulement dans la derni\`ere ligne.
\end{itemize}


\section{Simulation d'alignements} \label{sec:simu}
\'Ecrire un programme qui simule un grand nombre d'alignements 
de chaines de m\^eme longueur prises au hasard et observer
la r\'epartition des scores. La compr\'ehension de cette r\'epartition
est essentielle pour \'eviter de donner un sens faux \`a des
r\'esultats d'alignements obtenus en bilologie mol\'eculaire.

\section{Le programme de s\'equen\c{c}age BLAST.} \label{sec:blast}
Les s\'equences se pr\'esentent sous forme de fichier texte au
format \verb|FASTA|: la premi\`ere ligne d\'ecrit la s\'equence
consid\'er\'ee (par exemple nom de la prot\'eine correspondante),
les lignes suivantes sont des suites de lettres choisies dans
l'alphabet ad\'equat (par exemple A, C, G, T pour l'ADN, 20
lettres possibles pour les prot\'eines, cf. le tutoriel).
Sur \verb|pcm2| le r\'epertoire \verb|/var/db/blast/data/| contient
des bases de donn\'ees d'ADN et de prot\'eines au format \verb|FASTA|
ce qui peut servir d'exemple.

Le programme \verb|blast| peut \^etre ex\'ecut\'e en mode distant
depuis un navigateur Internet, (cherchez BLAST dans les signets 
Bioinformatique, on y trouve \'egalement un tutoriel (en anglais)).\\
Avantages: l'interface est plus conviviale et il n'y a
pas de probl\`eme d'installation
ni de mise-\`a-jour du logiciel et des bases de donn\'ees.\\
Inconv\'enients: il n'y a aucun contr\^ole de la puissance de calcul 
allou\'ee et d\'epend du bon fonctionnement du r\'eseau, 
ne peut \^etre utilis\'e \`a l'int\'erieur d'un autre programme.

On a donc aussi install\'e \verb|blast| sur le serveur local
\verb|pcm2| avec un nombre limit\'e de bases de donn\'ees de s\'equences
d'ADN et de prot\'eines. Ce qui permet \'egalement de s'exercer \`a
l'usage du logiciel en \'evitant d'encombrer ces serveurs Internet.\\
En mode local, il existe deux commandes~: \verb|blastall| et \verb|blastpgp|.
\begin{enumerate}
\item Exemples de lignes de commandes pour tester contre une s\'equence
de la base \verb|ecoli.nt| (comparaison d'ADN):\\
\verb|blastall -p blastn -d /var/db/blast/data/ecoli.nt -i test.txt|\\
Cette commande compare la s\'equence du fichier \verb|test.txt| avec
la base de donn\'ee \verb|ecoli.nt| et affiche le r\'esultat de la
comparaison. Utilisez Shift-PageUp et Shift-PageDown pour parcourir
les r\'esultats. En utilisant l'option \verb|-o test.out| ou la redirection
\verb|> test.out|, on peut envoyer le r\'esultat dans le fichier 
\verb|test.out|, puis on peut le lire avec un \'editeur de texte, par
exemple~:\\
\verb|emacs test.out &|
\item Options importantes de \verb|blastall|:\\
\verb|-p ...|: \verb|blastp| ou \verb|blastn| (prot\'eine ou nucl\'eotide), 
\verb|blastx| (nucl\'eotide contre prot\'eine avec toutes les possibilit\'es
de lecture), \verb|tblastn| (prot\'eine par rapport \`a nucl\'eotide,
toutes possibilit\'es de lecture) ou \verb|tblastx| (nucl\'eotide contre
nucl\'eotide toutes possibilit\'es de lecture).\footnote{blastn and 
blastp offer fully gapped alignments.  
blastx and tblastn have 'in-frame' gapped alignments and use sum 
statistics to link alignments from different frames. tblastx provides 
only ungapped alignments.}
\verb|-d ...|: nom de la base de s\'equences par d\'efaut \verb|nr|.
Plusieurs noms de base de donn\'ee peuvent \^etre donn\'es.\\
\verb|-i ...|: nom du fichier contenant la s\'equence \`a tester\\
\verb|-o ...|: nom du fichier contenant le r\'esultat de l'analyse.\\
Taper \verb|blastall -| pour avoir la liste compl\`ete des options.
\item \verb|blastpgp|:\\
Il s'agit d'une version sp\'eciale de \verb|blastp| qui recherche
des profils, ce qui peut par exemple servir \`a rechercher des s\'equences 
homologues\footnote{The PSI-BLAST program uses the information from any 
significant alignments returned to construct a position-specific score matrix, 
which replaces the query sequence for the next round of database searching. 
PSI-BLAST may be iterated until no new significant alignments are found.}
\end{enumerate}


\end{document}