Alignement de séquences biologiques
Algorithme de Needleman et Wunsch en Rust
Article mis en ligne le 14 août 2021
dernière modification le 30 août 2021

par Laurent Bloch

Cet article est un nouveau chapitre de mes exercices de programmation en Rust, précédé de Premiers programmes en Rust, de Programmation Rust, suite, de Codage de séquences biologiques avec Rust et de Configurer des séquences biologiques pour les aligner. Il arrive au cœur du sujet, l’algorithme de Needleman et Wunsch pour l’alignement de séquences, pour lequel j’ai déjà publié les solutions en Scheme et en Python (pour ce dernier, désolé, il faudra acheter le livre, sans d’ailleurs pour autant vous priver de celui de Scheme). La présentation de cet algorithme doit beaucoup à William Saurin.

 Programmation dynamique

L’algorithme de Needleman et Wunsch calcule un alignement optimal de deux chaînes de caractères en fonction du coût attribué aux trous (pénalité de gap, par exemple lorsqu’une chaîne est plus courte que l’autre) et aux mutations (lorsque l’on aligne deux caractères différents ; notre implémentation, plutôt que d’infliger une pénalité aux mutations, donnera une prime aux concordances (matches)). C’est un algorithme de programmation dynamique, ce qui consiste à conserver en mémoire (ici sous la forme d’un tableau) les résultats intermédiaires partiels pour ne pas avoir à les recalculer de multiples fois ; dans le cas de Needleman et Wunsch le gain est spectaculaire, puisqu’il permet de passer d’une complexité en O(kn) pour comparer des chaînes de longueur n à O(mn) pour deux chaînes de longueurs m et n. Les applications courantes d’alignement considèrent des séquences de plusieurs centaines de nucléotides ou d’acides aminés, en sachant que les banques de séquences nucléiques contiennent des génomes bactériens entiers, soit plusieurs millions de nucléotides, même s’il n’y a sans doute pas grand sens à aligner des génomes entiers.

Pour la description de l’algorithme je renvoie le lecteur aux articles cités dans le chapeau de celui-ci, j’évoquerai surtout ici les particularités de son écriture en Rust. L’algorithme direct exposé ici calcule un score de similitude en fonction des paramètres précisés ci-dessous, il restera encore à traiter dans un article suivant le problème du retour sur trace (backtracking).

 Acquisition et codage des données de test

Depuis le précédent article j’ai amélioré le programme d’acquisition et de mise en forme des données, je renvoie à cet article pour la nouvelle version. On notera la modification du format du tuple qui représente une séquence : la ligne d’identification obtenue de la première ligne du format FASTA est désormais de type String, donc codée en UTF8 pour laisser la porte ouverte à d’autres écritures que celle de l’anglais. J’ai testé avec les sinogrammes « simplifiés » et avec l’alphabet arabe. Le type d’une séquence est désormais (String, Vec<u8>).

Pour pouvoir les saisir dans la ligne de commande, j’ai introduit dans la structure Config les paramètres match_bonus (prime de concordance si les deux caractères à aligner sont les mêmes) et gap_penalty (pénalité de trou si l’on aligne un caractère avec un espace, ce paramètre doit avoir en principe une valeur négative) ; ces paramètres seront transmis à la fonction build_matrix du module build_sequences_matrix, objet du présent article ; ils sont de type f32 (virgule flottante simple précision).

La lecture de caractères Ascii représentés conformément au type u8 est inconfortable, je les convertis désormais en UTF8 pour l’affichage, sans avoir besoin de recourir à la fonction :

  1. pub const unsafe fn from_utf8_unchecked(v: &[u8]) -> &str

dont l’attribut unsafe priverait le programme de la garantie d’intégrité de la mémoire, mais en procédant ainsi, par exemple :

  1. ident = str::from_utf8(fasta_sequence.split_at(index).0).unwrap().to_string();
  2.  
  3. // ou, si seq2 est de type &Vec<u8> :
  4.  
  5. print!("{} ", char::from(seq2[i-1]));

Télécharger

Afin de faciliter les tests, voici deux petites séquences factices, identiques à celles des livres de Scheme et de Python, et deux vraies séquences de gènes cousins de deux espèces d’orchidées :

  1. > Première séquence
  2. GAATTCAGTTA
  3.  
  4. > Seconde séquence
  5. GGATCGA
  6.  
  7. >gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
  8. CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG
  9. AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG
  10. CCGCCTCGGGAGCGTCCATGGCGGGTTTGAACCTCTAGCCCGGCGCAGTTTGGGCGCCAAGCCATATGAA
  11. AGCATCACCGGCGAATGGCATTGTCTTCCCCAAAACCCGGAGCGGCGGCGTGCTGTCGCGTGCCCAATGA
  12. ATTTTGATGACTCTCGCAAACGGGAATCTTGGCTCTTTGCATCGGATGGAAGGACGCAGCGAAATGCGAT
  13. AAGTGGTGTGAATTGCAAGATCCCGTGAACCATCGAGTCTTTTGAACGCAAGTTGCGCCCGAGGCCATCA
  14. GGCTAAGGGCACGCCTGCTTGGGCGTCGCGCTTCGTCTCTCTCCTGCCAATGCTTGCCCGGCATACAGCC
  15. AGGCCGGCGTGGTGCGGATGTGAAAGATTGGCCCCTTGTGCCTAGGTGCGGCGGGTCCAAGAGCTGGTGT
  16. TTTGATGGCCCGGAACCCGGCAAGAGGTGGACGGATGCTGGCAGCAGCTGCCGTGCGAATCCCCCATGTT
  17. GTCGTGCTTGTCGGACAGGCAGGAGAACCCTTCCGAACCCCAATGGAGGGCGGTTGACCGCCATTCGGAT
  18. GTGACCCCAGGTCAGGCGGGGGCACCCGCTGAGTTTACGC
  19.  
  20. >gi|2765657|emb|Z78532.1|CCZ78532 C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA
  21. CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAACAGAATATATGATCGAGTG
  22. AATCTGGAGGACCTGTGGTAACTCAGCTCGTCGTGGCACTGCTTTTGTCGTGACCCTGCTTTGTTGTTGG
  23. GCCTCCTCAAGAGCTTTCATGGCAGGTTTGAACTTTAGTACGGTGCAGTTTGCGCCAAGTCATATAAAGC
  24. ATCACTGATGAATGACATTATTGTCAGAAAAAATCAGAGGGGCAGTATGCTACTGAGCATGCCAGTGAAT
  25. TTTTATGACTCTCGCAACGGATATCTTGGCTCTAACATCGATGAAGAACGCAGCTAAATGCGATAAGTGG
  26. TGTGAATTGCAGAATCCCGTGAACCATCGAGTCTTTGAACGCAAGTTGCGCTCGAGGCCATCAGGCTAAG
  27. GGCACGCCTGCCTGGGCGTCGTGTGTTGCGTCTCTCCTACCAATGCTTGCTTGGCATATCGCTAAGCTGG
  28. CATTATACGGATGTGAATGATTGGCCCCTTGTGCCTAGGTGCGGTGGGTCTAAGGATTGTTGCTTTGATG
  29. GGTAGGAATGTGGCACGAGGTGGAGAATGCTAACAGTCATAAGGCTGCTATTTGAATCCCCCATGTTGTT
  30. GTATTTTTTCGAACCTACACAAGAACCTAATTGAACCCCAATGGAGCTAAAATAACCATTGGGCAGTTGA
  31. TTTCCATTCAGATGCGACCCCAGGTCAGGCGGGGCCACCCGCTGAGTTGAGGC

Télécharger

 gi|2765658| donne l’identifiant GenBank de la séquence ;
 emb|Z78533.1| donne son identifiant EMBL (European Molecular Biology Laboratory).

 Discipline de typage

La plupart des langages modernes appliquent un typage des données plus ou moins strict, mais tous ne tirent pas les conséquences ultimes de cette pétition de principe.

Rust et son ancêtre Ada sont parmi les plus exigeants sous ce rapport. Un typage exigeant impose plus de discipline au programmeur, mais permet la production de programmes plus sûrs, ce qui réduit d’autant les coûts de maintenance (c’était l’exigence primordiale du cahier des charges rédigé par le Department of Defense (DoD) américain qui a donné naissance à Ada). Il faut savoir qu’au cours de son cycle de vie un logiciel peut engendrer des coûts de maintenance bien supérieurs aux coûts de développement initiaux, sans parler des dégâts provoqués par des programmes critiques faux, dont l’exemple le plus spectaculaire mais pas forcément le plus onéreux fut la destruction du premier exemplaire de la fusée Ariane V.

Un typage plus rigoureux des données, en contrôlant mieux notamment les conversions entre types, numériques ou autres, permet aussi de produire du code plus efficace et plus compact.

Ainsi, il est bon d’observer que l’égalité n’est pas définie pour les nombres en virgule flottante, mais il est mieux d’en tirer toutes les conséquences, ce que fait Rust. Que le lecteur se rassure, on peut quand même comparer des flottants entre eux, mais il faut le dire, ce qui implique que l’on a conscience qu’il s’agit d’un ordre partiel (essentiellement en raison de la valeur spéciale NaN, Not a Number, qui peut être le résultat d’une opération). On écrit donc par exemple (j’ai mis un certain temps à trouver la solution) :

  1. let tmp = f32::max(v2, v3);

Les indices de vecteurs et de chaînes sont de type usize (taille d’agrégat, non signée), on ne peut les multiplier par un flottant sans le dire :

  1. the_mat.set(0, j, gap_penalty * j as f32) ;

On remarque au passage la notation de style par objets pour appliquer à la case (0,j) de la matrice the_mat la méthode set qui modifie la valeur de son contenu. Obtenir cette valeur pour un calcul nécessite de la déballer :

  1. the_mat.get(i,j-1).unwrap()

Ces exigences semblent bien fatigantes aux utilisateurs de langages plus laxistes, mais elles sont le prix d’un logiciel plus sûr, plus robuste et plus efficace.

 Calcul de la matrice des scores

Conformément à l’algorithme décrit ici le calcul des scores de similitude utilise une matrice définie comme suit, si la longueur de la première séquence est l_seq1 et celle de la seconde l_seq2 ; la première séquence sera disposée horizontalement, d’où le nombre de colonnes égal à l_seq1+1, et le seconde verticalement, d’où le nombre de lignes égal à l_seq2+1 :

  1. use simple_matrix::Matrix;
  2. ...
  3. let mut the_mat: Matrix::<f32> = Matrix::new(l_seq2+1, l_seq1+1);

Télécharger

La valeur de chaque case de la matrice, d’indice (i, j), qui correspond à la position i-1 dans la seconde séquence et à la position (j-1) de la première séquence, se calcule ainsi :

  1. for j in 1..col {
  2.             the_mat.set(0, j, gap_penalty * j as f32) ;
  3. }
  4. let mut score: f32 = 0.0;
  5. for i in 1..lin {
  6.             the_mat.set(i, 0, gap_penalty * i as f32) ;
  7.     for j in 1..col {
  8.         if seq1[j-1] == seq2[i-1] {
  9.             score = match_bonus} else {
  10.             score = 0.0}
  11.         the_mat.set(i, j, max3(the_mat.get(i-1,j-1).unwrap()
  12.                                + score,
  13.                                the_mat.get(i-1,j).unwrap()
  14.                                + gap_penalty,
  15.                                the_mat.get(i,j-1).unwrap()
  16.                                + gap_penalty));
  17.     }
  18. }

Télécharger

Ce que l’on peut représenter par le schéma suivant :

Voici le module de calcul de la matrice des scores dans son ensemble :

  1. // src/sequences_matrix/build_sequences_matrix.rs :
  2.  
  3. // This module was inspired by Vincent Esche's Seal crate,
  4. // but simplified and much more basic, without mmap and so on.
  5. // For pedagogic use.
  6.  
  7. pub mod build_sequences_matrix {
  8.  
  9.     use simple_matrix::Matrix;
  10.     use std::str;
  11.     use std::char;
  12.    
  13.     pub fn build_matrix(sequence1: (String, Vec<u8>),
  14.                         sequence2: (String, Vec<u8>),
  15.                         match_bonus: f32,
  16.                         gap_penalty: f32) {
  17. //                      -> Matrix<f32> {
  18.         print_seq(&sequence1);
  19.         print_seq(&sequence2);
  20.  
  21.         let l_seq1: usize = (sequence1.1).len();
  22.         let l_seq2: usize = (sequence2.1).len();
  23.  
  24.         println!("Longueur première séquence : {} ", l_seq1);
  25.         println!("Longueur seconde séquence : {} ", l_seq2);
  26.        
  27.         let mut the_mat: Matrix::<f32> = Matrix::new(l_seq2+1, l_seq1+1);
  28.  
  29.         init_matrix(&mut the_mat, l_seq2+1, l_seq1+1, 0.0);
  30.  
  31.         nw_matrix(&mut the_mat, l_seq2+1, l_seq1+1, match_bonus, gap_penalty, &sequence1.1, &sequence2.1);
  32.  
  33.         print_vector(&sequence1.1);
  34.            
  35.         print_matrix(&the_mat, &sequence2.1, l_seq2+1, l_seq1+1);
  36.     }
  37.  
  38.     pub fn nw_matrix(the_mat: &mut Matrix::<f32>,
  39.                      lin: usize,
  40.                      col: usize,
  41.                      match_bonus: f32,
  42.                      gap_penalty: f32,
  43.                      seq1: &Vec<u8>,
  44.                      seq2: &Vec<u8>) {
  45.         for j in 1..col {
  46.             the_mat.set(0, j, gap_penalty * j as f32) ;
  47.         }
  48.         let mut score: f32 = 0.0;
  49.         for i in 1..lin {
  50.             the_mat.set(i, 0, gap_penalty * i as f32) ;
  51.             for j in 1..col {
  52.                 if seq1[j-1] == seq2[i-1] {
  53.                     score = match_bonus} else {
  54.                     score = 0.0}
  55.                 the_mat.set(i, j, max3(the_mat.get(i-1,j-1).unwrap()
  56.                                        + score,
  57.                                        the_mat.get(i-1,j).unwrap()
  58.                                        + gap_penalty,
  59.                                        the_mat.get(i,j-1).unwrap()
  60.                                        + gap_penalty));
  61.             }
  62.         }
  63.     }
  64.  
  65.     pub fn max3(v1: f32, v2: f32, v3: f32) -> f32 {
  66.         let tmp = f32::max(v2, v3);
  67.         if v1 > tmp {
  68.             return v1 } else {
  69.             return tmp };
  70.     }
  71.    
  72.     pub fn init_matrix(the_mat: &mut Matrix::<f32>, lin: usize, col: usize, val: f32) {
  73.         for i in 0..lin {
  74.             for j in 0..col {
  75.                 the_mat.set(i, j, val) ;
  76.             }
  77.         }
  78.     }
  79.  
  80. // "print_seq" affiche une séquence selon différents formats.
  81.     pub fn print_seq(sequence: &(String, Vec<u8>)) {
  82.                 println!("Ident : {:?}", sequence.0);
  83.                 println!("Séquence : {:?}", sequence.1);
  84.                 let sequence_str = str::from_utf8(&sequence.1).unwrap().to_string();
  85.                 println!("Séquence : {}", &sequence_str);
  86.     }
  87.  
  88.     pub fn print_vector(the_vec: &Vec<u8>) {
  89.         let vec_str = str::from_utf8(the_vec).unwrap().to_string();
  90.         print!("{} ", "   ");
  91.         for c in vec_str.chars() {
  92.             print!("{} ", c);
  93.         }
  94.         print!("{}", "   \n");
  95.     }
  96.    
  97.     pub fn print_matrix(the_mat: &Matrix::<f32>, seq2: &Vec<u8>, lin: usize, col: usize) {
  98.         for i in 0..lin {
  99.             if i > 0 {print!("{} ", char::from(seq2[i-1]))} else
  100.             {print!("{} ", " ")};
  101.             for j in 0..col {
  102.                 print!("{} ", the_mat.get(i, j).unwrap());
  103.             }
  104.             print!("{}", "\n")
  105.         }
  106.     }
  107. }

Télécharger

Ce programme s’invoque ainsi, pour une prime de concordance égale à 1 et une pénalité de trou égale à 0,5 :

  1. cargo run Data/sequence3.fasta Data/sequence4.fasta 1 -0.5

Si nous mettons la pénalité de trou à 0, ainsi :

  1. cargo run Data/sequence3.fasta Data/sequence4.fasta 1 0

la matrice résultante sera :

G A A T T C A G T T A
0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 3
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 4 4 5 5 5 5
A 0 1 2 3 3 3 4 5 5 5 5 6

Le score final est donc 6 (en principe dans la dernière case en bas à droite).

 Et ensuite ?

Quel résultat ce programme nous donne-t-il ? Le score de l’alignement de deux séquences. Souvent ce résultat suffit : si la première séquence est celle que nous étudions, et que nous cherchons dans une banque celles qui lui ressemblent le plus, le calcul des scores nous permettra de sélectionner les 20 ou les 50 séquences les plus proches, selon les paramètres que nous aurons fournis. Souvent c’est ce que l’on veut.

Si maintenant nous voulons calculer les alignements correspondants et les représenter, par exemple ainsi :

G A A T T C A G T T A
¦ ¦ ¦ ¦ ¦ ¦
G G A _ T C _ G _ _ A

il faut remonter dans le tableau à partir de la dernière case en reconstituant le chemin parcouru, ce qui donne à chaque pas quels caractères ont été mis en correspondance, ainsi :

G A A T T C A G T T A
0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 3
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 4 4 5 5 5 5
A 0 1 2 3 3 3 4 5 5 5 5 6

La programmation de ce retour sur trace fera l’objet d’un prochain article.