Articles

alignement par paires de séquences nucléotidiques à l’aide de correspondances exactes maximales

octobre 15, 2020 by admin

approche

dans notre algorithme proposé, la première étape vers l’alignement des séquences consiste à extraire les MEMs entre les séquences en les comparant directement. La Figure 3a est un exemple qui compare une cible et une séquence de requête où CTC et AAA sont deux MEMs identifiés par la comparaison. Chaque groupe de symboles identiques continus dans la comparaison donne un MEM même s’il n’est composé que d’un seul symbole correspondant., Afin d’extraire tous les MEMs entre les séquences, la séquence de requête doit être décalée vers la droite et vers la gauche un symbole à la fois (voir Fig. 3b). Après chaque quart de travail, l’étape de comparaison doit être répétée pour identifier de nouveaux MEMs. Par exemple, la troisième ligne de la Fig. 3b représente le cas où la séquence de requête est décalée vers le symbole de droite et est comparée à la séquence cible. Le résultat de la comparaison identifie AAAAGC comme un nouveau MEM. Tous les autres MEMs extraits par les opérations shift et compare sont également mis en évidence dans la Fig. 3b., Trois des MEMs (Mx,My et Mz) sont mis en évidence avec des couleurs différentes à utiliser pour une explication ultérieure.

dans le modèle de notation affine-gap, le score d’alignement tel qu’il est calculé à L’aide D’Eq., 1 Où Nm est le nombre de correspondances recevant chacune un score de correspondance de Rm,Nx est le nombre de discordances recevant chacune une pénalité de non-correspondance De Px,No est le nombre d’ouvertures d’écart recevant chacune une pénalité d’ouverture D’écart de Po et Ng est la longueur totale de tous les écarts, chaque écart recevant une pénalité Il y aurait une ouverture d’écart pour chaque groupe d’écart continu. Par exemple, s’il y a deux espaces dans l’alignement, où la longueur du premier espace est de trois et la longueur du deuxième espace est de quatre, alors il y a deux ouvertures d’espace (No=2) et la longueur totale de l’espace est de sept (Ng=3+4=7).,

$$ {}COMME = (N_{m} \times R_{m}) – ((N_{x} \times P_{x}) + (N_{o} \times P_{o}) + (N_{g} \times P_{g})) $$

(1)

etant Donné la liste de tous les MEMs, l’alignement peut être calculée à l’aide partielle des alignements. Par exemple, considérez MEMs Mx, My et Mz dans la Fig. 3b. les alignements partiels effectués en prenant différentes combinaisons de Mx, My et Mz ainsi que le nombre de correspondances, de discordances et de lacunes, ainsi que les scores d’alignement résultants sont illustrés à la Fig. 4. L’alignement qui inclut uniquement Mx et Mz donne le score d’alignement le plus élevé., Notez que My et Mz se chevauchent et lorsque les deux sont considérés dans le même alignement, le chevauchement est exclu de Mz. Considérant tous les MEMs dans la Fig. 3b entraîne beaucoup plus de combinaisons où aucune d’entre elles n’atteint un score plus élevé.

l’Examen de toutes les combinaisons possibles de MEMs serait exhaustive., Dans la section « algorithme D’alignement », nous décrivons un nouvel algorithme de programmation dynamique DP-MEM qui trouve efficacement la meilleure combinaison sans considérer tous les cas. DP-MEM doit savoir quelles parties des séquences correspondent mais pas les symboles réels dans les séquences. L’entrée dans DP-MEM est le positionnement des MEMs dans la cible et dans les séquences de requête obtenues lors du processus d’extraction mem décrit dans la section « extraction MEM”., La façon dont les MEMs sont représentés avec leurs positions et la façon dont le nombre de correspondances, de discordances et d’écarts est calculé lorsque les MEMs sont combinés dans un alignement sont expliquées dans le reste de cette section. La Figure 5 est un autre exemple d’alignement avec six MEMs (M1 à M6) qui forment l’alignement entre la séquence cible T et la séquence de requête Q. Pour plus de simplicité, il n’y a pas de chevauchement entre MEMs dans cet exemple. Chaque MEM Mi est représenté comme un triplet de nombres entiers: les positions de départ en T et en Q (STi et SQI respectivement) et sa longueur (Li)., Les positions de fin en T et en Q peuvent être calculées plus tard (Φ2E de L’algorithme 2). Le tableau 1 répertorie la longueur et le positionnement de M1 à M6 en T et en Q.

Tableau 1 début et de fin de la position de MEMs dans la Fig., 5

Tableau 2 Calcul du nombre de discordances et d’écarts entre les MEMs dans la Fig. 5

en cas d’incompatibilité et d’écart entre Mi et Mj, tous les écarts sont considérés comme continus pour réduire la pénalité ouverte d’écart (une seule pénalité ouverte d’écart est appliquée pour un écart continu). Ainsi, pour tous les MEMs adjacents qui ont des écarts entre eux, une seule pénalité ouverte d’écart est appliquée., Le placement des inadéquations et le seul écart continu n’est pas important, car cela n’affecterait pas le score d’alignement. Nous supposons que la pénalité de non-concordance est constante (c’est habituel pour les séquences D’ADN).

extraction MEM

Il existe des méthodes pour extraire des Correspondances exactes maximales entre de longues séquences telles qu’un génome entier. Cependant, ces méthodes sont basées sur le prétraitement et l’indexation d’une ou des deux séquences, ce qui est une opération fastidieuse. Par exemple, dans DNA read aligner, le génome de référence est indexé une fois, et le même index est utilisé chaque fois qu’une nouvelle lecture est alignée., Nous recherchons un algorithme rapide pour identifier les MEMs entre des séquences relativement courtes qui changent pour chaque alignement. Une méthode de force brute pour ce problème (fichier supplémentaire 1: Section II) est lente et inefficace(avec la complexité de O (n3)). Nous proposons une méthode parallèle rapide au niveau des bits pour accélérer le processus d’extraction MEM. Notre méthode D’extraction MEM est basée sur les opérations de décalage et de comparaison illustrées à la Fig. 3b. la première étape consiste à représenter des séquences avec des vecteurs de bits, où A, C, T et G sont codés comme 00, 01, 10 et 11, respectivement (fichier supplémentaire 1: Section III)., La Figure 6 illustre un exemple de paire de séquences, ainsi que les représentations bit-vectorielles correspondantes. Dans un ordinateur de base, le mot machine est généralement de 64 bits pouvant accueillir 32 symboles nucléotidiques. Comme une séquence est généralement supérieure à 32 Symboles, le vecteur binaire correspondant est stocké dans plusieurs mots machine. Chaque opération sur des vecteurs de bits de séquences de taille n symboles agit sur $\lceil \ frac {n} {32} \rceil $ mots machine.

avec la représentation bit-vectorielle des séquences, décaler une séquence d’un symbole revient à décaler le bit-vectoriel de deux bits, et comparer les séquences peut être fait avec l’instruction XOR (32 symboles à la fois). Dans la sortie XOR (X), 00 signifie que les symboles sont mis en correspondance, et une séquence de 00s affiche un MEM., Un ensemble d’opérations shift et bitwise comme indiqué dans L’algorithme 1 calcule X et ensuite le vecteur bit edge (E) dans lequel le début et la fin de chaque MEM sont mis en évidence avec des Bits set (bits avec une valeur de un). La Figure 6 montre les vecteurs de bits X et E avec des MEMs en surbrillance. Le positionnement des MEMs dans les séquences est ensuite calculé à partir du vecteur bit edge (fichier supplémentaire 1: Section IV).,

algorithme D’alignement

dans la section « approche”, nous montrons qu’en considérant différentes combinaisons de MEMs et en calculant le score d’alignement pour l’alignement correspondant, on peut identifier la combinaison de MEMs qui aboutit au score d’alignement maximum. Cependant, l’examen de toutes les combinaisons possibles de MEMs est une solution naïve. Une façon plus systématique de trouver l’alignement efficacement est d’utiliser la programmation dynamique.

la programmation dynamique est la méthode d’approche de la solution à un problème en définissant et en résolvant des sous-problèmes plus petits., Les Solutions pour les sous-problèmes sont utilisées pour résoudre un problème plus important à chaque étape. Le processus est répété jusqu’à ce que tous les sous-problèmes soient résolus. Finalement, la solution à l’un des sous-problèmes serait la solution au problème initial. Lorsque tous les sous-problèmes sont résolus, un processus de retour en arrière identifie une série de solutions qui contribuent à la solution finale. En programmation dynamique, il devrait y avoir un ordre des données d’entrée le long duquel la procédure de récursivité se déroule.

nous trions tous les MEMs en fonction de la position de leur extrémité dans la séquence de requête (EQ)., Les MEMs qui se terminent dans la même position sont ordonnés de manière arbitraire. Le sous-problème jth consiste à trouver l’alignement des sous-séquences de T et Q qui se terminent par jth MEM MJ (T et Q respectivement). Nous montrerons que cet ordre de MEM est suffisant pour supporter la récursivité correcte.

dans la liste triée des MEMs, EQi=EQj indique que L’un des Mi ou Mj chevauche complètement l’autre MEM dans la séquence de requête. Puisque dans Φ2B de L’algorithme 2, la région de chevauchement est exclue, Mi et Mj ne peuvent pas être dans le même alignement., Ainsi, ith et jth sous-problèmes sont résolus indépendamment les uns des autres et l’ordre de i et j dans la liste triée pourrait être arbitraire. Si EQk>EQj (k> j dans la liste triée), Mk ne peut pas faire partie de l’alignement qui se termine par Mj. Ainsi, le jème sous-problème peut être résolu indépendamment de la solution au kème sous-problème. Notez qu’il est également possible de trier les MEMs en fonction de leur position de fin dans la séquence cible (ET) en utilisant une justification similaire.

notre algorithme de programmation dynamique proposé (DP-MEM) est élaboré dans L’algorithme 2., Pour L’exemple MEMs extrait à la Fig. 3b, la table de programmation dynamique et la valeur intermédiaire calculées dans l’algorithme sont représentées sur les figures. 7 et 8 respectivement. L’entrée de DP-MEM est la liste des MEMs où chaque MEM (Mj) est un triplet d’entiers . La deuxième entrée n est le nombre de MEMs dans la liste. La sortie s est le score d’alignement pour les séquences. L’algorithme imprime les indices de tous les MEMs qui forment l’alignement où le premier et le dernier nombre imprimé sont les indices des MEMs les plus à droite et les plus à gauche de l’alignement respectivement., Toutes les étapes de L’algorithme 2 sont commentées comme suit:

Φ1: notation de chaque MEM pour tous ses symboles correspondants., Notez qu’il existe des symboles correspondant à Lj dans Mj. Sj représente le score d’alignement le plus élevé pour l’alignement se terminant à Mj. L’initialisation de Sj dans cette étape est similaire au calcul du score d’alignement partiel lorsque seul Mj est inclus dans l’alignement. W est utilisé pour le retour en arrière. La valeur de -1 indique que le SJ actuel est obtenu en considérant MJ seul dans l’alignement.
Φ2: calcul de Sj pour chaque MEM (Mj)., Pour calculer Sj, pour chaque MEM Mi où Mi apparaît avant Mj dans la liste, l’algorithme ajoute Mj à l’alignement se terminant à Mi (extension des alignements précédemment trouvés) et recherche l’extension qui maximise Sj (résolution d’un sous-problème plus grand en utilisant des sous-problèmes précédemment résolus).
Φ2A: ignorer l’extension lorsque ce n’est pas possible. Si eti > ETj alors Mi contient une partie de la séquence cible qui est au-delà de l’alignement se terminant à Mj et l’extension n’est pas possible. Si EQi=EQj ou ETi=ETj ou SQI≥SQJ ou STi≥STj, alors l’un des MEMs chevauche complètement l’autre MEM., Dans ce cas, Mi et Mj ne peuvent pas être alignés ensemble.
Φ2B: Calcul de la longueur du chevauchement entre Mi et Mj. Si ${MO}_{i}^{j}$ est inférieur ou égal à zéro, alors aucun chevauchement n’existe.
Φ2C: conserver une copie de Mj avant d’exclure le chevauchement.
Φ2D: si le chevauchement existe, à l’exclusion du chevauchement de Mj
Φ2E: calcul de la position finale de Mj en T et Q.
Φ2F: calcul de la distance (Nombre de symboles) entre Mi et MJ en T et Q.
Φ2G: calcul du nombre,
Φ2H: calcul de la pénalité pour les discordances et les écarts entre Mi et Mj ($p_{i}^{j}$). S’il existe un écart, une seule pénalité ouverte est soustraite.
Φ2I: calcul du score d’alignement $\left (s_{i}^{j}\right)$ lorsque Mj est ajouté à L’alignement se terminant à Mi. Le score de tous les symboles correspondants dans Mj (LJ×Rm) est ajouté au score d’alignement pour L’alignement se terminant à Mi (Si). Ensuite, la pénalité pour les écarts et les discordances entre Mi et MJ$\left (p_{i}^{j}\right)$ est soustraite.,
Φ2J: si l’extension de Mj à L’alignement se terminant par Mi donne un score $\left (s_{i}^{j}\right)$ supérieur au score actuel pour Mj (Sj) alors le nouveau score est stocké dans Sj. W est également défini sur i pour garder une trace du Mi qui maximise le score pour Mj.
Φ2K: restauration de la valeur de Mj avant l’exclusion afin que Mj puisse être utilisé dans d’autres extensions d’alignement.
Φ3: recherche du MEM avec le Sj le plus élevé. Ce MEM est le dernier MEM dans l’alignement (Me)., Le score le plus élevé (Se) est renvoyé sous la forme S qui est le score d’alignement le plus élevé pour les séquences données. L’index du MEM qui maximise Sj est stocké dans e pour commencer à revenir en arrière de moi.
Φ4: dans l’alignement, le MEM précédent immédiat pour moi est celui qui maximise le score d’alignement pour moi. L’index d’un tel mem est stocké dans W. par conséquent, l’itération de f←W Visite l’index de tous les MEMs dans l’alignement. Lorsque W est égal à -1, Mf est le premier MEM dans l’alignement, et l’itération est arrêté.,

Dans notre algorithme, nous ne pénaliserait pas les incohérences et les lacunes avant la première MEM et après la dernière MEM dans l’alignement. Il en résulte un algorithme d’alignement local. En considérant ces pénalités, l’algorithme génère un alignement global (fichier supplémentaire 1: Section V).

l’équation à calculer $p_{I}^{j}$ dans Φ2H de L’algorithme 2 suppose qu’il n’y a pas de symbole correspondant entre T et Q dans la zone entre Mi et Mj (tous les symboles sont comptés comme des incompatibilités ou des lacunes)., Bien que cette hypothèse ne soit pas vraie pour tous les Mi, elle est toujours vraie pour le Mi qui conduit au maximum $s_{i}^{j}$ ce qui annule l’effet de l’hypothèse étant incorrecte pour les autres Mi. À titre de preuve, supposons qu’il existe un symbole correspondant dans la zone entre Mi et Mj. Le symbole correspondant serait un MEM (Mk). Mk est déjà étendu à L’alignement se terminant à Mi. Ainsi, lors de l’extension de Mj à Mk, un score plus élevé est obtenu par rapport à l’extension de MJ à Mi.

le chaînage des graines colinéaires comme indiqué dans ont été largement utilisés dans l’alignement de grandes séquences, c’est-à-dire l’alignement génome à génome., Il a également été utilisé pour identifier les régions candidates pour une lecture étant donné un ensemble de MEMs dans BWA. Les algorithmes de chaînage avec scoring sont similaires à L’algorithme de programmation dynamique que nous avons proposé (DP-MEM). Cependant, il existe des différences qui rendent DP-MEM approprié pour l’alignement par paires de courtes séquences. DP-MEM prend en compte le fait que tous les MEMs dans une certaine taille d’écart sont fournis dans l’entrée et optimise le nombre d’itérations dans l’algorithme. DP-MEM met également en œuvre une approche heuristique pour compenser l’effet des MEMs courts retirés de la liste d’entrée résultant des écarts entre MEMs.,

Optimisation

étant donné des séquences de longueur n, L’algorithme d’extraction de MEMs (fourni dans la section « extraction MEM”) nécessite 2(n−1) shift et 2n−1 compare des opérations sur des vecteurs de bits (chacun agit sur $\lceil \frac {n}{32} \rceil $ mots machine) qui aboutissent à un algorithme de complexité O(n2) pour produire des vecteurs de bits de bord pour la paire de séquences donnée. La complexité de la fonction qui calcule le positionnement des MEMs à partir du vecteur bit edge et les trie en fonction de L’égaliseur n’a pas encore été ajoutée. De plus, si m MEMs sont extraits, Φ2 de L’algorithme 2 (DP-MEM) a la complexité de O(m2)., Compte tenu de la complexité de l’extraction MEM et DP-mem, MEM-Align est beaucoup plus lent que les algorithmes d’alignement disponibles. Pour accélérer le processus, nous présentons plusieurs optimisations pour MEM-Align qui permet une exécution plus rapide en sacrifiant la précision. D’autre part, nous introduisons des méthodes pour améliorer la précision avec une perte de performance minimale.

alignement des bandes

l’alignement des bandes est une méthode heuristique connue pour accélérer le processus d’alignement. Cette technique limite le motif des lacunes dans l’alignement (fichier supplémentaire 1: Section VI)., Par conséquent, si l’alignement entre deux séquences ne suit pas ce modèle, l’algorithme n’identifiera pas l’alignement. En programmation dynamique traditionnelle, l’alignement est obtenu après avoir calculé la valeur de toutes les cellules du tableau. Cependant, avec l’optimisation de l’alignement des bandes, seules les cellules du diamètre et proches de la diagonale sont évaluées. gl est la largeur de la bande dans l’alignement des bandes où les cellules plus loin que gl au diamètre ne sont pas évaluées. Cellules plus sombres dans la Fig. 1 montre le cas où gl=1.,

contrairement à l’approche de programmation dynamique traditionnelle, MEM-Align n’a pas de table similaire pour appliquer l’alignement en bandes. Cependant, nous avons constaté que nous pouvons simuler le même effet en limitant le nombre d’opérations de décalage dans le MEM processus d’extraction. Par exemple, si nous décalons la séquence de requête vers gl vers la droite et vers la gauche, nous obtenons un alignement en bandes avec la bande de gl. L’alignement en bandes réduit la complexité de l’extraction MEM de O(n2) à O(N.(2GL+1)) où gl est une valeur petite et fixe. Ainsi, la complexité de l’extraction MEM est O(n) lorsque l’alignement par bandes est appliqué., De plus, avec ledit alignement en bandes, il est probable que moins de MEMs sont extraits ce qui accélère les étapes algorithmiques suivantes.

Short mem removal

Nous avons observé que la majorité des MEMs extraits sont courts et sont le résultat de symboles appariés aléatoirement. Pour accélérer MEM-Align, les MEMs plus courts que sl sont filtrés pendant le processus d’extraction MEM. Cela réduit le nombre de MEMs à extraire et à traiter (accélérant ensuite l’algorithme). Le filtrage de mem Court se fait en étendant L’algorithme 1 avec un ensemble d’opérations de décalage et de bits (fichier supplémentaire 1: Section VII).,

d’autre part, il existe de rares cas dans lesquels court MEMs font partie de l’alignement; par exemple, un symbole correspondant entouré par l’asymétrie. Sans avoir tous les MEMs dans la liste d’entrée, DP-MEM ne peut pas trouver le même alignement qu’il trouve quand tous les MEMs existent dans la liste d’entrée. Afin de compenser la perte de MEMs courts dans l’entrée, nous modifions Φ2H de DP-MEM pour envisager la possibilité d’avoir des correspondances courtes entre MEMs (fichier supplémentaire 1: Section VIII).

Il pourrait y avoir des cas plus difficiles où dans l’alignement, plusieurs MEMs courts existent entre deux MEMs (voir Fig. 9)., La seule façon d’identifier correctement le score pour la zone entre Mi et MJ dans Φ2H est d’appliquer un alignement global à cette région. Cependant, Φ2H est une opération fréquente et devrait rester rapide. Par conséquent, nous avons décidé de surmonter partiellement le problème en limitant les cas possibles (une méthode heuristique).

Si il y a des lacunes dans la zone entre Mi et Mj, nous supposons qu’il n’y est qu’un seul écart, soit à l’extrémité gauche ou à droite à la fin de la zone. Ensuite, seuls deux alignements sont possibles pour la zone., Le nombre de symboles correspondants est compté pour les deux cas de manière séquentielle et celui qui entraîne des correspondances maximales est considéré comme le nombre de correspondances entre Mi et Mj (fichier supplémentaire 1: Section IX). La comparaison séquentielle est une opération coûteuse et nous concevons une méthode pour éviter la comparaison séquentielle lorsque cela est possible (fichier supplémentaire 1: Section X).

tout autre cas qui ne correspond pas à l’hypothèse ci-dessus entraîne un alignement avec un score inférieur., Cependant, compte tenu du faible taux de lacunes et de inadéquations, la possibilité d’avoir plusieurs lacunes et inadéquations dans une petite zone est faible.

extension d’alignement efficace

dans Φ2 de DP-MEM, Mj étend tous les alignements qui se terminent par {M1…Mj−1} (si possible). Cependant, pour chaque Mj, il existe un sous-ensemble plus petit Ωj {{M1…Mj−1} tel qu’en étendant Mj à tous les alignements se terminant par un Mi∈Ωj, l’alignement qui se termine par Mj est trouvé (Eq. 2). En d’autres termes, il y aurait moins de $s_{i}^{j}$ à évaluer. Définition de L’ensemble Ωj et de la preuve D’Eq. 2 sont fournis dans le fichier supplémentaire 1: Section XI., La définition de Ωj est affectée lorsque l’optimisation de la suppression de MEM court est appliquée (fichier supplémentaire 1: Section XII).

$$ \max\limits_{M_{i} \in \Omega_{j}}{S_{i}^{j}} = \max\limits_{1 \leq i \leq j-1}{S_{i}^{j}} $$

(2)

Hybride alignement

afin De maintenir la précision de l’algorithme, nous avons décidé d’utiliser une méthode hybride qui est une combinaison de MEM-Aligner et de Smith-Waterman algorithme. Nous définissons trois cas dans lesquels MEM-Align peut être inexact., Si l’alignement d’une paire de séquences tombe dans l’un de ces cas, nous utilisons L’algorithme de Smith-Waterman pour aligner les séquences. Ces cas sont:

lorsque les séquences sont répétitives et que le nombre de MEMs extraits dépasse le seuil TM. Nous avons constaté que MEM-Align est susceptible de produire un alignement inexact lors de l’alignement de séquences répétitives. Une valeur TM appropriée diminue les chances de signaler un alignement inexact avec une augmentation négligeable du temps de traitement moyen.,
lorsque le score d’alignement calculé pour L’alignement généré par MEM-Align est inférieur à un seuil TS. Ce cas se produit principalement lorsqu’il y a un espace dans l’alignement qui ne peut pas être identifié en raison de l’alignement des bandes.
Lorsqu’il n’existe pas de MEM plus long que sl à extraire (un cas rare)., Si sl est défini sur une valeur élevée et que la similitude entre les séquences est faible,

bien que l’envoi de paires de séquences à un algorithme externe entraîne un calcul supplémentaire, le nombre de séquences envoyées à l’algorithme externe reste faible si des valeurs appropriées sont choisies pour TM et TS.

ignorer les MEMs distants

lorsque la distance entre Mi et Mj est grande, il est peu probable que Mi et Mj soient des MEMs adjacents dans l’alignement., Par conséquent, l’algorithme ignore l’extension si la distance entre Mi et Mj est plus longue qu’un seuil TD (réduisant encore le nombre de $s_{i}^{j}$ à évaluer). Cette optimisation améliore légèrement les performances avec un effet secondaire négligeable sur la précision.

Be Able