From silico.biotoul.fr

Il s'agit de programmer en perl l'algorithme de Needleman-Wunsch exploitant la programmation dynamique décrit dans Alignement de deux séquences d'acides nucléiques.

Etant données 2 séquences nucléiques ou protéiques, une matrice contenant le score du meilleur alignement est progressivement remplie.

Vous commencerez par programmer la méthode pour un alignement global sans brèche en utilisant un score d'homologie :

• identité = +1
• substitution = -1
• indel = -2

Ce score sera représenté par une matrice de substitution assez triviale pour les séquences nucléiques :

   -  A  T  C  G
-  1 -2 -2 -2 -2
A -2  1 -1 -1 -1
T -2 -1  1 -1 -1
C -2 -1 -1  1 -1
G -2 -1 -1 -1  1

Cela se traduit en perl de la manière suivante :

# Default substitution matrix (homology)
#                         -   A   T   C   G
my $homologyMatrix = [ [  9, -2, -2, -2, -2 ], # -
                       [ -2,  1, -1, -1, -1 ], # A
                       [ -2, -1,  1, -1, -1 ], # T or U
                       [ -2, -1, -1,  1, -1 ], # C
                       [ -2, -1, -1, -1,  1 ]  # G
];

Le fait d'utiliser une telle matrice rendra le programme générique en terme du type de séquences alignées, qu'elles soient nucléiques ou protéiques. Pour les séquences protéiques, l'utilisateur pourra passer quelle matrice utiliser en paramètres, par exemple PAM120 ou BLOSUM62, mais vous commencerez par des séquences nucléiques.

Une fois la matrice de programmation dynamique remplie, faites-la afficher en indiquant le(s) chemin(s) correspondant au(x) meilleur(s) alignement(s). Par exemple, pour les séquences TTGTG et ACTTGCAG :

           A      C      T      T      G      C      A      G 
    0   - -2   - -4     -6     -8    -10    -12    -14    -16 
                      \                                       
 T -2     -1     -3     -3     -5     -7     -9    -11    -13 
                             \                                
 T -4     -3     -2     -2     -2     -4     -6     -8    -10 
                                    \                         
 G -6     -5     -4     -3     -3     -1   - -3     -5     -7 
                                           \      \           
 T -8     -7     -6     -3     -2     -3     -2   - -4     -6 
                                                         \    
 G-10     -9     -8     -5     -4     -1     -3     -3     -3 

Faites ensuite afficher le(s) meilleur(s) alignement(s), par exemple :

--TTG-TG
ACTTGCAG

--TTGT-G
ACTTGCAG

Pour passer les paramètres à votre programme, utilisez le module Getopt (cf. la documentation CPAN) :

#!/usr/bin/perl -w
 
use strict;
use Data::Dumper;
use Getopt::Long;
 
my %options;
 
GetOptions(\%options,
	'help', 'sequences=s@'
);
 
usage() if $options{help};
 
sub usage {
print STDERR << "EOF";
 
  Dynamic programming sequence alignment.
 
 usage: $0 [--help] [--sequences first_sequence --sequences second_sequence]
 
EOF
	exit;
}
 
# default sequences used for development
my $s1 = $options{sequences}[0] ||   'TTGTG';
my $s2 = $options{sequences}[1] || 'ACTTGCAG';

Algorithme 1 : alignement global

DYNAMIC-PROGRAMMING(x,m,y,n)
 D(0,0) <- 0
 for i = 1..m
   do D(i,0) <- D(i-1,0) + Del(x[i])
 for j = 1..n
   do D(0,j) <- D(0,j-1) + Ins(y[j])
   for i = 1..m
     do D(i,j) <- max{D(i-1,j-1)+Sub(x[i],y[j]),
                      D(i-1,j) + Del(x[i]),
                      D(i,j-1) + Ins(y[j])}
 return D(m,n)

Alignement avec brèches

Fonction de coût d’une brèche : br(k) = g + h × (k − 1) avec k sa longueur, g le coût de son ouverture et h le coût de sa prolongation. Les relations suivantes permettent de définir un algorithme mettant en oeuvre l’alignement avec brèches (s’inspirer de l’algorithme DYNAMIC-PROGRAMMING), notamment au moyen des matrice I et S en plus de la matrice D déjà utilisée dans DYNAMIC-PROGRAMMING :

Pour i = 0..m − 1 et j = 0..n − 1 :
 S[i, j] = min{S[i − 1, j] + h, D[i − 1, j] + g},
 I[i, j] = min{I[i, j − 1] + h, D[i, j − 1] + g},
 D[i, j] = min{D[i − 1, j − 1] + Sub(x[i], y[j]), S[i, j], I[i, j]} pour la récurrence, et :
 S[−1, −1] = S[i, −1] = S[−1, j] = ∞,
 I[−1, −1] = I[i, −1] = I[−1, j] = ∞ pour l’initialisation de S et T , et enfin pour D :
 D[−1, −1] = 0,
 D[0, −1] = g,
 D[−1, 0 = g,
 D[i, −1] = D[i − 1, −1] + h,
 D[−1, j] = D[−1, j − 1] + h.