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 :
dont l’attribut unsafe
priverait le programme de la garantie d’intégrité de la mémoire, mais en procédant ainsi, par exemple :
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 :
– 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) :
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 :
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 :
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
:
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 :
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 :
Ce programme s’invoque ainsi, pour une prime de concordance égale à 1 et une pénalité de trou égale à 0,5 :
Si nous mettons la pénalité de trou à 0, ainsi :
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.