Alignement multiple

Prérequis

Cette section prolonge l'alignement global de la section 2.3 et suppose connus les scores additifs, les pénalités de gap et les matrices de substitution présentées en 2.5. Elle utilise aussi la lecture d'un fichier FASTA vue en 1.2. Les séquences seront notées $s_1, \dots, s_k$, de longueurs $n_1, \dots, n_k$, et les positions restent 1-based dans le texte.

Introduction

Un alignement multiple est le point de départ de presque toute analyse comparative : domaines conservés, sites actifs, amorces, sélection ou phylogénie. Deux protéines homologues peuvent partager un domaine tout en ayant accumulé insertions et délétions ; un ARN ribosomique conserve des régions structurées séparées par des boucles variables. Il faut donc décider quelles positions se correspondent à travers toute la famille, pas seulement entre deux membres choisis au hasard.

Le problème computationnel change de nature lorsque l'on passe de deux à $k$ séquences. L'alignement par paires se résout exactement en temps quadratique. L'alignement multiple exact explore un espace de dimension $k$ : son coût croît exponentiellement avec le nombre de séquences, même pour un critère additif simple. Les outils courants utilisent donc des heuristiques progressives, guidées par un arbre de similarité puis éventuellement raffinées. Elles sont rapides, mais ne garantissent pas l'optimum et propagent les erreurs précoces.

La difficulté ne s'arrête pas à l'alignement brut. Les colonnes riches en gaps, les extrémités mal alignées et les symboles ambigus influencent l'arbre final ; un trimming mal documenté peut retirer des sites informatifs ou laisser des artefacts. Nous suivons donc une chaîne complète : score, complexité, stratégie, alignement, nettoyage et vérification. Un jeu synthétique de quatre taxons illustre les commandes sans prétendre reproduire une analyse biologique.

Théorie

Alignement par paires et alignement multiple

Un alignement par paires de deux séquences $s$ et $t$ est une suite de colonnes qui apparie un caractère de $s$ à un caractère de $t$, ou un caractère à un gap. Un alignement multiple de $k$ séquences est une matrice à $k$ lignes et $L$ colonnes : chaque ligne est une séquence dans laquelle des gaps ont été insérés, toutes les lignes ont la même longueur, et aucune colonne n'est entièrement composée de gaps. Les caractères de chaque séquence conservent leur ordre.

La différence essentielle est la cohérence. Si l'on aligne séparément toutes les paires d'une famille, on obtient $\binom{k}{2}$ alignements qui peuvent être incompatibles : une paire peut préférer un gap à une position, une autre paire le préférer ailleurs, et aucune matrice commune ne peut satisfaire les deux. Un alignement multiple impose une seule histoire de colonnes. En retirant les lignes autres que $i$ et $j$, puis les colonnes où ces deux lignes ont un gap, on obtient l'alignement par paires induit $A_{ij}$. C'est cet objet, et non un alignement par paires calculé isolément, qui entre dans le score.

Un alignement multiple n'est donc pas la concaténation d'alignements deux à deux. Les désaccords entre outils sont fréquents : deux heuristiques peuvent produire des matrices différentes, toutes deux plausibles, sans qu'aucune ne soit optimale au sens d'un critère global.

Le score somme-de-paires

Le critère le plus classique est le score somme-de-paires (sum-of-pairs, SP). Pour un alignement $A$ de $k$ séquences, on additionne les scores de tous les alignements par paires induits :

Si le score d'une paire est additif par colonne, le calcul se fait colonne par colonne. Pour une colonne contenant des résidus et des gaps, chaque paire de lignes contribue par un score de substitution si les deux résidus sont présents, par une pénalité de gap si exactement une des deux lignes porte un gap, et par zéro si les deux portent un gap. La dernière règle évite de compter deux fois un événement d'insertion déjà représenté ailleurs ; elle correspond à la suppression de la colonne dans l'alignement induit.

Le score SP peut être pondéré pour éviter qu'un groupe de séquences très proches domine le résultat. Il reste un compromis : selon les paramètres, deux substitutions rares peuvent coûter moins qu'un long gap. Un score élevé ne prouve donc pas que l'alignement est biologiquement correct ; il indique seulement sa cohérence avec le modèle choisi.

Complexité : pourquoi l'exact est hors de portée

Pour $k$ séquences de longueur maximale $L$, une programmation dynamique exacte sur le score SP doit décrire l'alignement de $k$ préfixes. Le nombre d'états est de l'ordre de $L^k$, et chaque état compare jusqu'à $2^k - 1$ transitions. Le temps est donc $O(2^k L^k)$ et la mémoire $O(L^k)$ pour la matrice complète. Pour $k=3$ et $L=100$, un million d'états reste accessible ; pour $k=10$ et $L=100$, $10^{20}$ états ne le sont plus.

Le problème d'optimisation du score SP est NP-difficile lorsque $k$ fait partie de l'entrée. Les outils pratiques abandonnent donc l'optimalité globale. Ils calculent d'abord des distances entre paires, ce qui coûte $O(k^2 L^2)$ avec un alignement naïf, puis construisent un arbre guide et fusionnent les séquences ou les profils. Le coût des fusions dépend de l'implémentation, mais il reste polynomial. Le prix de cette réduction est l'absence de garantie : un mauvais choix précoce ne sera pas nécessairement corrigé.

Heuristiques progressives et arbre guide

Une heuristique progressive procède en trois temps. Elle estime une distance entre chaque paire, souvent à partir d'un alignement global corrigé pour les substitutions multiples, construit un arbre guide par regroupement hiérarchique (Neighbor-Joining ou clustering rapide), puis fusionne les feuilles des plus proches aux plus éloignées. Chaque fusion aligne deux profils, c'est-à-dire deux groupes déjà alignés.

Un profil contient, pour chaque colonne, la fréquence des résidus et l'information de gap. Aligner deux profils revient à aligner deux matrices de profils avec une programmation dynamique dont le score intègre ces fréquences. Une fois qu'un gap a été introduit dans un profil, il est généralement conservé lors des fusions suivantes : c'est la règle du « once a gap, always a gap ». Elle stabilise l'alignement, mais elle transforme une erreur précoce en erreur définitive.

Les raffinements itératifs tentent de corriger ce défaut. Ils retirent un groupe, réaligne les séquences restantes, puis réinsèrent le groupe et conservent la version qui améliore le score. MAFFT et Clustal Omega proposent de telles itérations. Elles améliorent souvent la qualité, au prix du temps de calcul. L'arbre guide doit être vu comme un ordre de fusion, pas comme l'arbre phylogénétique final : il est estimé sur des distances approximatives et sert uniquement à l'alignement.

Clustal Omega : mBed et HHalign

Clustal Omega, décrit par Sievers et al. en 2011, combine deux idées pour passer à l'échelle. La construction de l'arbre guide repose sur mBed, une méthode qui estime les distances à partir d'un petit ensemble de séquences de référence et produit un arbre en temps quasi linéaire en nombre de séquences. L'alignement progressif n'utilise pas un simple score de profil, mais des profils HMM comparés par HHalign, l'algorithme de la suite HH-suite. Cette combinaison permet de traiter des jeux de plusieurs milliers de séquences avec une qualité souvent proche de méthodes plus coûteuses.

Les options utiles sont peu nombreuses : --outfmt=clu ou fasta, --iterations pour un raffinement, --guidetree-out pour archiver l'arbre guide et --threads pour le parallélisme. La version 1.2.4 est celle de ce cours. Les limites restent celles des méthodes progressives : l'arbre guide est approximatif, les premiers gaps sont difficiles à défaire et les régions divergentes peuvent être mal alignées. Inspectez l'alignement, pas seulement son score.

MAFFT : précision ou vitesse

MAFFT propose plusieurs stratégies. Le mode L-INS-i (--localpair --maxiterate 1000) réaligne les séquences par paires locales et raffine le résultat par itérations ; il est recommandé pour un jeu modéré, de quelques dizaines à quelques centaines de séquences divergentes, lorsque la précision prime. Le mode FFT-NS-2 (--retree 2 --maxiterate 0) est une méthode progressive rapide fondée sur la transformée de Fourier rapide, adaptée aux grands jeux et aux séquences proches. Entre les deux, G-INS-i privilégie un alignement global et E-INS-i gère les domaines séparés par de longues insertions.

Le mode change le résultat et le temps de calcul. --auto laisse MAFFT choisir selon la taille du jeu : documentez la stratégie retenue, car elle fait partie du résultat. La version 7.526 est celle du cours. Un alignement MAFFT n'est pas « meilleur » dans l'absolu, seulement pour un compromis donné entre temps, mémoire et précision. Comparer deux modes sur un sous-jeu est souvent plus instructif qu'un débat théorique.

Trimming, colonnes conservées et ambiguïté

Toutes les colonnes ne méritent pas d'entrer dans l'analyse. La fraction de gaps d'une colonne est le nombre de gaps divisé par le nombre de séquences ; l'occupation est son complément. Une colonne conservée n'a aucun gap et tous ses caractères sont identiques. Une colonne variable contient au moins deux résidus distincts. Une colonne informative au sens de la parcimonie contient au moins deux états distincts présents chacun au moins deux fois : elle peut soutenir un regroupement, contrairement à une colonne où l'état minoritaire n'apparaît qu'une fois.

trimAl automatise le nettoyage. L'option -gt fixe un seuil de gaps, mais son sens est inversé par rapport à l'intuition : -gt 0.8 conserve les colonnes dont la fraction de gaps est au plus $1 - 0.8 = 0.2$. L'option -automated1 choisit une heuristique automatique adaptée à la phylogénie ; elle ne se réduit pas au seul seuil de gaps et peut être plus agressive. La sortie doit être comparée à l'alignement non trimmé, et la commande archivée.

Les symboles d'ambiguïté compliquent le tableau. N signifie « base inconnue », R « purine », Y « pyrimidine » ; ce ne sont pas des gaps. Une colonne N,N,N a une occupation de 1 et un score de conservation de 1, alors qu'elle n'apporte aucune information. De même, un gap - signale une absence dans l'alignement, pas nécessairement une délétion biologique : il peut provenir d'une région mal séquencée ou d'un artefact d'alignement. Avant de calculer des distances, il faut décider comment traiter gaps et ambiguïtés, et l'écrire.

Exemple travaillé à la main

Prenons trois séquences alignées sur quatre colonnes :

s1  ACGT
s2  AC-T
s3  ATGT

Avec un score d'appariement $+1$, une substitution $-1$ et un gap linéaire $-1$, calculons le score SP. Colonne 1 : A, A, A, trois appariements identiques, contribution $3$. Colonne 2 : C, C, T, contributions $+1 -1 -1 = -1$. Colonne 3 : G, gap, G, contributions $-1 +1 -1 = -1$. Colonne 4 : T, T, T, contribution $3$. Le total vaut $3 - 1 - 1 + 3 = 4$.

Vérifions par les paires induites. $A_{12}$ vaut ACGT contre AC-T, soit $1+1-1+1=2$. $A_{13}$ vaut ACGT contre ATGT, soit $1-1+1+1=2$. $A_{23}$ vaut AC-T contre ATGT, soit $1-1-1+1=0$. La somme $2+2+0=4$ confirme le calcul colonne par colonne. L'occupation de la troisième colonne vaut $2/3$, sa conservation vaut $1$ car les deux résidus observés sont des G. Le trimming à 30 % de gaps retire cette colonne, alors qu'un seuil de 34 % la conserve. Ce petit exemple montre qu'un score élevé de conservation peut coexister avec un gap.

Implémentation

Les commandes ci-dessous supposent les outils installés et exécutées depuis la racine du dépôt. Elles écrivent leurs sorties dans le répertoire courant ; dans un projet réel, créez un dossier de résultats daté et conservez les fichiers intermédiaires.

clustalo -i datasets/phylogeny/demo.fasta -o demo.clustal \
  --outfmt=clu --force --verbose
mafft --localpair --maxiterate 1000 datasets/phylogeny/demo.fasta > demo.linsi.fasta
mafft --retree 2 --maxiterate 0 datasets/phylogeny/demo.fasta > demo.fftns2.fasta
trimal -in demo.clustal -out demo.gt50.clustal -clustal -gt 0.5
trimal -in demo.clustal -out demo.auto.clustal -clustal -automated1

Clustal Omega écrit un alignement CLUSTAL, MAFFT un FASTA aligné et trimAl un alignement au format demandé. Vérifiez toujours le format réellement produit : --outfmt=clu et -clustal ne sont pas les valeurs par défaut de tous les outils, et une extension .aln ne garantit rien.

Le code Python suivant montre une version compacte du parser et des statistiques de colonnes. Il est volontairement autonome : il ne lit pas de fichier et ne dépend que de la bibliothèque standard.

Aucune sortie. Exécutez le code pour voir le résultat.

La sortie vérifiée le 2026-11-01 avec Python 3.13 est :

{'seq1': 'ACGT', 'seq2': 'AC-T', 'seq3': 'ATGT'}
[1.0, 1.0, 0.67, 1.0]
{'seq1': 'ACT', 'seq2': 'ACT', 'seq3': 'ATT'}
[1.0, 0.67, 1.0, 1.0]

Trois points méritent attention. D'abord, les lignes de consensus commencent par un espace et sont ignorées ; un parser qui se contente de découper chaque ligne en deux champs les interpréterait comme des séquences. Ensuite, la longueur $L$ doit être vérifiée après concaténation des blocs, pas bloc par bloc : un fichier peut avoir des blocs de tailles différentes mais des séquences finales de même longueur. Enfin, trim conserve les colonnes dont la fraction de gaps est inférieure ou égale au seuil. Le seuil 0.3 retire la troisième colonne, dont la fraction vaut $1/3 \approx 0.33$.

La version des exercices ajoute les contrôles de type, le refus des longueurs différentes et la reconnaissance des symboles d'ambiguïté, ainsi que les gaps . de certains formats : ces détails évitent des résultats silencieusement faux.

Choisir et documenter les paramètres

Le choix des paramètres dépend de la question. Pour une phylogénie de gènes codants, on aligne les acides aminés ou les codons, on vérifie le cadre de lecture, puis on retire les colonnes ambiguës. Pour rechercher un domaine, un alignement local ou un profil HMM peut être plus pertinent qu'un alignement global. Pour un jeu de milliers de séquences, la vitesse de FFT-NS-2 ou de Clustal Omega peut primer sur un gain de précision marginal.

La reproductibilité repose sur une trace écrite : nom et version de l'outil, commande exacte, format d'entrée et de sortie, seuil de trimming, matrice et pénalités de gap, date et somme de contrôle des fichiers. Si l'outil expose une graine aléatoire, fixez-la. Conservez l'alignement non trimmé à côté de l'alignement filtré, et notez la longueur et le nombre de colonnes conservées. Enfin, comparez au moins deux stratégies sur un sous-jeu : un désaccord localisé est une information, pas un détail à masquer.

Étude de cas

Exemple synthétique. Le fichier datasets/phylogeny/demo.fasta contient quatre séquences de 16 nucléotides créées pour ce cours. Elles ne représentent aucune espèce réelle et ne permettent aucune conclusion évolutive ; leur intérêt est de rendre colonnes et gaps faciles à inspecter.

La question est double : les outils produisent-ils le même alignement sur un jeu sans indel, et le trimming change-t-il quelque chose ? Nous procédons en quatre étapes : inspecter le FASTA, aligner avec Clustal Omega, comparer avec les deux modes MAFFT, puis appliquer trimAl et calculer les statistiques.

clustalo -i datasets/phylogeny/demo.fasta -o demo.clustal --outfmt=clu --force
mafft --localpair --maxiterate 1000 datasets/phylogeny/demo.fasta > demo.linsi.fasta
mafft --retree 2 --maxiterate 0 datasets/phylogeny/demo.fasta > demo.fftns2.fasta
trimal -in demo.clustal -out demo.gt50.clustal -clustal -gt 0.5

La sortie illustrative de Clustal Omega est la suivante. Les quatre séquences ont la même longueur et aucun gap ; l'alignement se réduit donc à une superposition.

CLUSTAL O(1.2.4) multiple sequence alignment
 
taxon_A          ACGTACGTACGTACGT
taxon_B          ACGTACGTACGTACGA
taxon_C          ACGTTCGTACGTACGT
taxon_D          TTGTACGTACGTACGT
                 "  ** ********** "

Les guillemets de la dernière ligne sont ajoutés pour rendre les espaces visibles ; ils ne font pas partie du fichier. La ligne de consensus signale les colonnes identiques. Les positions 1, 2, 5 et 16 sont variables : A/A/A/T, C/C/C/T, A/A/T/A et T/A/T/T. Les douze autres colonnes sont conservées. Aucune colonne ne contient de gap.

Le rapport de alignment_report donne : 4 séquences, 16 colonnes, 12 conservées, 4 variables, 0 colonne vide, une fraction de gaps nulle, une occupation moyenne de 1,0 et une conservation moyenne de $15/16 = 0.9375$. Les quatre colonnes variables ne sont pas informatives : l'état minoritaire n'y apparaît qu'une fois. Le rapport émet donc un avertissement, même si l'alignement est valide.

Le trimming illustre un point souvent mal compris. Avec -gt 0.5, trimAl conserve les colonnes dont la fraction de gaps est au plus 0,5 ; sans gap, aucune colonne n'est retirée et la longueur reste 16 (sortie illustrative). Ce seuil ne peut pas corriger un manque de signal évolutif. -automated1 applique une autre heuristique : comparez sa sortie à l'alignement non trimmé et notez la longueur obtenue.

Les séquences ne diffèrent que par une à trois positions. Les aligneurs produisent donc le même alignement non gapé (sortie illustrative) : la comparaison porte sur les statistiques, pas sur la supériorité d'un outil. Pour départager Clustal Omega et MAFFT, il faudrait un jeu réel plus long, avec insertions et régions divergentes. Ce mini-jeu vérifie seulement que le pipeline, le parser et les métriques sont cohérents avant les données biologiques.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-5/5.1-alignement-multiple -q

python -m pytest exercises/module-5/5.1-alignement-multiple -q

python -m pytest exercises/module-5/5.1-alignement-multiple/projet -q

Résumé

• Un alignement multiple est une matrice cohérente ; il ne se réduit pas à une collection d'alignements par paires.
• Le score somme-de-paires additionne les scores des alignements induits pour toutes les paires de séquences.
• L'optimisation exacte du score SP coûte $O(2^k L^k)$ et devient rapidement hors de portée ; le problème est NP-difficile.
• Les heuristiques progressives utilisent un arbre guide et des profils ; une erreur précoce est difficile à corriger.
• Clustal Omega combine des arbres guides mBed et des profils HMM via HHalign ; MAFFT propose L-INS-i pour la précision et FFT-NS-2 pour la vitesse.
• trimAl filtre les colonnes selon les gaps ou une heuristique ; -gt 0.8 autorise au plus 20 % de gaps.
• Une colonne conservée n'est pas forcément informative, et un symbole d'ambiguïté n'est pas un gap.
• La reproductibilité exige la version, la commande, le seuil et la conservation de l'alignement non trimmé.

Références

• Sievers, F. et al. (2011). Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Molecular Systems Biology, 7, 539.
• Katoh, K. et Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular Biology and Evolution, 30(4), 772–780.
• Capella-Gutiérrez, S., Silla-Martínez, J. M. et Gabaldón, T. (2009). trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics, 25(15), 1972–1973.
• Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research, 32(5), 1792–1797.
• Wang, L. et Jiang, T. (1994). On the complexity of multiple sequence alignment. Journal of Computational Biology, 1(4), 337–348.
• Documentation Clustal Omega, EMBL-EBI. https://www.ebi.ac.uk/Tools/msa/clustalo/
• Documentation MAFFT. https://mafft.cbrc.jp/alignment/software/
• Documentation trimAl. https://trimal.cgenomics.org/