Articles

alinierea Pairwise a secvențelor nucleotidice folosind potriviri exacte maxime

octombrie 15, 2020 by admin

abordare

în algoritmul nostru propus, primul pas spre alinierea secvențelor este extragerea MEMs între secvențe prin compararea directă a acestora. Figura 3A este un exemplu care compară o țintă și o secvență de interogare în care CTC și AAA sunt două MEMs identificate prin comparație. Fiecare grup de simboluri identice continue în comparație, duce la o MEM, chiar dacă acesta este compus dintr-un singur simbol de potrivire., Pentru a extrage toate MEMs între secvențe, secvența de interogare trebuie să fie deplasată până la dreapta și la stânga câte un simbol la un moment dat (vezi Fig. 3b). După fiecare schimbare, pasul de comparație trebuie repetat pentru a identifica noi MEMs. De exemplu, a treia linie din Fig. 3B reprezintă cazul în care secvența de interogare este deplasată spre simbolul din dreapta și este comparată cu secvența țintă. Rezultatul comparației identifică AAAAGC ca un nou MEM. Toate celelalte MEMs extrase prin operații shift și compare sunt, de asemenea, evidențiate în Fig. 3b., Trei dintre MEMs (Mx,My și Mz) sunt evidențiate cu culori diferite pentru a fi utilizate pentru explicații ulterioare.

În afine-diferența de notare model, scorul de aliniere CA este calculată folosind Ecuația., 1 Unde Nm este numărul de meciuri fiecare primind un scor de meci de Rm, Nx este numărul de neconcordanțe fiecare primind o penalizare de nepotrivire de Px, nu este numărul de deschideri de goluri fiecare primind o penalizare deschisă de goluri de Po și Ng este lungimea totală a tuturor golurilor, fiecare decalaj primind o penalizare de prelungire a decalajului de Pg. Ar exista o deschidere decalaj pentru fiecare grup de decalaj continuu. De exemplu, dacă există două goluri în aliniere, unde lungimea primului decalaj este de trei și lungimea celui de-al doilea decalaj este de patru, atunci există două deschideri de decalaj (nu=2) și lungimea totală a decalajului este de șapte (Ng=3+4=7).,

$ $ {} AS = (n_{m} \times R_{m}) – ((N_{x} \times p_{x}) + (n_{o} \times p_{o}) + (n_{g} \times p_{g})) $ $

(1)

având în vedere lista tuturor MEMs, alinierea poate fi calculată folosind aliniamente parțiale. De exemplu, luați în considerare MEMs Mx, My și Mz în Fig. 3B. aliniamentele parțiale realizate prin luarea diferitelor combinații de Mx, My și Mz împreună cu numărul de potriviri, nepotriviri și goluri, precum și scorurile de aliniere rezultate sunt prezentate în Fig. 4. Alinierea care include doar MX și MZ are ca rezultat cel mai mare scor de aliniere., Rețineți că, My și Mz se suprapun reciproc și atunci când ambele sunt considerate în aceeași aliniere suprapunerea este exclusă din Mz. Având în vedere toate MEMs din Fig. Rezultatele 3B în mai multe combinații în cazul în care nici unul dintre ei atinge un scor mai mare.

Examinarea tuturor combinație posibilă de MEMs-ar fi exhaustivă., În secțiunea” algoritmul de aliniere ” descriem un nou algoritm dinamic de programare DP-MEM care găsește eficient cea mai bună combinație fără a lua în considerare toate cazurile. DP-MEM trebuie să știe ce părți ale secvențelor se potrivesc, dar nu simbolurile reale din secvențe. Intrarea în DP-MEM este poziționarea MEMs în țintă și în secvențele de interogare care sunt obținute în timpul procesului de extracție MEM descris în secțiunea „extracție MEM”., Modul în care MEMs sunt reprezentate cu pozițiile lor și modul în care numărul de potriviri, nepotriviri și lacune sunt calculate atunci când MEMs sunt combinate într-o aliniere sunt explicate în restul acestei secțiuni. Figura 5 este un alt exemplu de aliniere cu șase MEMs (M1 la M6) care formează alinierea între secvența țintă T și secvența de interogare Q. pentru simplitate, nu există nicio suprapunere între MEMs în acest exemplu. Fiecare MEM Mi este reprezentat ca un triplet de numere întregi: pozițiile de pornire în T și în Q (STi și respectiv SQi) și lungimea sa (Li)., Pozițiile finale în T și în Q pot fi calculate mai târziu (Φ2E din algoritmul 2). Tabelul 1 enumeră lungimea și poziționarea M1 până la M6 în T și în Q.

Tabelul 1 poziția de Pornire și se termină de MEMs în Fig., 5

Tabelul 2 Calcul numărul de neconcordanțe și lacune între MEMs în Fig. 5

în cazul În care există neconcordanțe și lacune între Mi și Mj, toate golurile sunt considerate continuă să reducă decalajul de pedeapsă deschis (doar un singur gol deschis de penalizare se aplică pentru o continuă gap). Astfel, pentru toate MEMs adiacente care au lacune între ele, se aplică o singură penalizare deschisă., Plasarea neconcordanțelor și singurul decalaj continuu nu este importantă, deoarece nu ar afecta scorul de aliniere. Presupunem că pedeapsa de nepotrivire este constantă (aceasta este obișnuită pentru secvențele ADN).

extracția MEM

există metode de extragere a potrivirilor exacte maxime între secvențe lungi, cum ar fi un întreg genom. Cu toate acestea, aceste metode se bazează pe preprocesarea și indexarea uneia sau a ambelor secvențe, care este o operație consumatoare de timp. De exemplu, în ADN read aligner, genomul de referință este indexat o dată și același indice este utilizat de fiecare dată când o nouă citire este aliniată., Căutăm un algoritm rapid pentru identificarea MEMs între secvențe relativ scurte care se schimbă pentru fiecare aliniere. O metodă de forță brută pentru această problemă (fișier suplimentar 1: Secțiunea II) este lentă și ineficientă(cu complexitatea O (n3)). Propunem o metodă paralelă rapidă la nivel de biți pentru a accelera procesul de extracție MEM. Metoda noastră de extracție MEM se bazează pe operațiunile de schimbare și comparare prezentate în Fig. 3b. primul pas este reprezentarea secvențelor cu vectori de biți, unde A, C, T și G sunt codificate ca 00, 01, 10 și, respectiv, 11 (fișier suplimentar 1: Secțiunea III)., Figura 6 ilustrează o pereche de secvențe de exemplu, împreună cu reprezentări bit-vectoriale corespunzătoare. Într-un computer de mărfuri, cuvântul mașină este de obicei 64 biți care pot găzdui 32 de simboluri nucleotide. Deoarece o secvență este de obicei mai mare de 32 de simboluri, vectorul de biți corespunzător este stocat în mai multe cuvinte de mașină. Fiecare operație pe vectorii de biți ai secvențelor de simboluri de mărime n acționează asupra cuvintelor mașinii $\lceil \frac {N}{32} \rceil$.

Cu pic-vectori de reprezentare a secvențelor, trecerea de la o secvență cu un simbol este la fel ca trecerea de la bit-vector de doi biți, și compararea secvențelor se poate face cu instrucțiunea XOR (32 de simboluri la un moment dat). În ieșirea XOR (X), 00 înseamnă că simbolurile sunt potrivite, iar o secvență de 00s arată un MEM., Un set de operații shift și bitwise așa cum se arată în algoritmul 1 calculează X și ulterior marginea bit-vector (E) în care începutul și sfârșitul fiecărui MEM sunt evidențiate cu biți set (biți cu o valoare de unul). Figura 6 prezintă X și e bit-vectori cu MEMs evidențiate. Poziționarea MEMs în secvențe sunt apoi calculate de la marginea bit-vector (fișier suplimentar 1: Secțiunea IV).,

Aliniere algoritm

În „Abordare” secțiune, vom arăta că prin luarea în considerare diferite combinații de MEMs și de calcul scorul de aliniere pentru aliniere corespunzătoare, se poate identifica o combinație de MEMs că rezultatele în maxim scorul de aliniere. Cu toate acestea, examinarea tuturor combinațiilor posibile de MEMs este o soluție naivă. Un mod mai sistematic de a găsi alinierea eficient este de a utiliza programarea dinamică.programarea dinamică este metoda de abordare a soluției la o problemă prin definirea și rezolvarea subproblemelor mai mici., Soluțiile pentru subprobleme sunt folosite pentru a rezolva o problemă mai mare la fiecare pas. Procesul se repetă până când toate subproblemele sunt rezolvate. În cele din urmă, soluția la unul dintre subprobleme ar fi soluția la problema inițială. Când toate subproblemele sunt rezolvate, un proces de backtracking identifică o serie de soluții care contribuie la soluția finală. În programarea dinamică, ar trebui să existe o ordonare a datelor de intrare de-a lungul cărora se desfășoară procedura de recursivitate.

sortăm toate MEMs în funcție de poziția sfârșitului lor în secvența de interogare (EQ)., MEMs care se termină în aceeași poziție sunt ordonate într-un mod arbitrar. Subproblem jth este de a găsi alinierea subsecțiuni de T și Q care se termină la jth MEM Mj (T și Q respectiv). Vom arăta că această ordonare a MEM este suficientă pentru a susține recursivitatea corectă.

în lista sortată de MEMs, EQi=EQj indică faptul că unul dintre Mi sau Mj se suprapune complet celuilalt MEM în secvența de interogare. Deoarece în Φ2B din algoritmul 2 regiunea de suprapunere este exclusă, Mi și Mj nu pot fi în aceeași aliniere., Astfel, subproblemele ith și jth sunt rezolvate independent una de cealaltă, iar ordinea i și j din lista sortată ar putea fi arbitrară. Dacă EQk>EQj (k>j în lista sortată), Mk-ar putea să nu fie o parte a procesului de aliniere care se termină în Mj. Astfel, subproblemele jth pot fi rezolvate independent de la soluție la subproblema kth. Rețineți că este de asemenea posibil să sortați MEMs pe baza poziției lor finale în secvența țintă (ET) folosind o justificare similară.algoritmul nostru de programare dinamică (DP-MEM) propus este elaborat în algoritmul 2., Pentru exemplul MEMs extras în Fig. 3B, tabelul de programare dinamică și valoarea intermediară calculată în algoritm sunt prezentate în Fig. 7 și respectiv 8. Intrarea în DP-MEM este lista MEMs unde fiecare MEM (Mj) este un triplet de numere întregi . A doua intrare n este numărul de MEMs din listă. Ieșirea S este scorul de aliniere pentru secvențe. Algoritmul imprimă indicii tuturor MEMs care formează alinierea unde primul și ultimul număr imprimat sunt indicii MEMs din dreapta și din stânga din aliniere., Toate etapele Algoritmului 2 sunt comentate în următoarele:

Φ1: Notare fiecare MEM pentru toate simboluri identice., Rețineți că există simboluri de potrivire Lj în MJ. Sj reprezintă cel mai mare scor de aliniere pentru alinierea care se termină la Mj. Inițializarea Sj în această etapă este similară calculării scorului de aliniere parțială atunci când numai Mj este inclus în aliniere. W este folosit pentru backtracking. Valoarea -1 indică faptul că SJ curent este obținut prin luarea în considerare Mj singur în alinierea.Φ2: calcul Sj pentru fiecare MEM (Mj)., Pentru a calcula Sj, pentru fiecare MEM Mi unde Mi apare înainte de Mj în listă, algoritmul adaugă MJ la alinierea care se termină la Mi (extinderea aliniamentelor găsite anterior) și caută extensia care maximizează Sj (rezolvarea unei subprobleme mai mari folosind subprobleme rezolvate anterior).Φ2A: săriți extensia atunci când nu este posibil. Dacă ETi> ETj, atunci Mi conține o parte din secvența țintă care este dincolo de alinierea care se termină la Mj și extensia nu este posibilă. Dacă EQi=EQj sau ETi=ETj sau SQI≥SQJ sau STi≥STj, atunci unul dintre MEMs se suprapune complet celuilalt mem., În acest caz, Mi și Mj nu pot fi într-o aliniere împreună.Φ2B: calculând lungimea suprapunerii dintre Mi și Mj. Dacă ${Mo}_{i}^{j}$ este mai mic sau egal cu zero, atunci nu există nicio suprapunere.Φ2C: păstrarea unei copii a Mj înainte de a exclude suprapunerea.
Φ2D: Dacă există suprapunere, cu excepția suprapunere de la Mj
Φ2E: Calcul poziția termină de Mj în T și Î.
Φ2F: Calcul distanța (numărul de simboluri) între Mi și Mj în T și Î.
Φ2G: Calcul numărul de neconcordanțe și lacune între Mi și Mj.,
Φ2H: calcularea pedepsei pentru nepotrivirile și decalajele dintre Mi și Mj ($P_{i}^{j}$). Dacă există un decalaj,se scade o singură penalizare deschisă.
Φ2I: Calcul scorul de aliniere $\left (S_{i}^{j}\right)$ când Mj este adăugat la alinierea se termină la Km. Scorul pentru toate simbolurile de potrivire din Mj (LJ×Rm) se adaugă la scorul de aliniere pentru alinierea care se termină la Mi (Si). Apoi se scade pedeapsa pentru golurile și nepotrivirile dintre Mi și Mj$\left (P_{i}^{j}\right)$.,Φ2J: dacă extinderea Mj la alinierea care se termină în rezultatele Mi într-un scor $\left (s_{i}^{j}\right)$ mai mare decât scorul curent pentru Mj (Sj), atunci noul scor este stocat în SJ. De asemenea, W este setat la i pentru a urmări Mi care maximizează scorul pentru Mj.Φ2K: restaurarea valoarea Mj înainte de excludere, astfel încât Mj poate fi utilizat în alte extensii de aliniere.Φ3: Privind pentru MEM cu cea mai mare Sj. Acest MEM este ultimul MEM din aliniere (Me)., Cel mai mare scor (Se) este returnat ca S, Care este cel mai mare scor de aliniere pentru secvențele date. Indicele MEM care maximizează Sj este stocat în e pentru a începe backtracking de la mine.Φ4: în aliniere, MEM-ul anterior imediat pentru mine este cel care maximizează scorul de aliniere pentru mine. Indicele unui astfel de MEM este stocat în W. ca urmare, iterația lui F←W vizitează indicele tuturor MEMs din aliniere. Când W este egal cu -1, Mf este primul MEM în aliniere și iterația este oprită.,

în algoritmul nostru, nu penalizăm neconcordanțele și lacunele înainte de prima MEM și după ultima MEM din aliniere. Acest lucru are ca rezultat un algoritm de aliniere locală. Luând în considerare aceste penalizări algoritmul generează o aliniere globală (fișier suplimentar 1: Secțiunea V).

ecuația de calcul $P_{i}^{j}$ în Φ2H din algoritmul 2 presupune că nu există un simbol de potrivire între T și Q în zona dintre Mi și Mj (toate simbolurile sunt considerate nepotriviri sau goluri)., Deși această ipoteză nu este adevărată pentru toate Mi, este întotdeauna adevărată pentru Mi care duce la maxim $s_{i}^{j}$ care anulează efectul presupunerii fiind incorect pentru alte Mi. Ca dovadă, să presupunem că există un simbol de potrivire în zona dintre Mi și Mj. Simbolul de potrivire ar fi un MEM (Mk). Mk este deja extins la alinierea care se termină la Mi. Astfel, atunci când se extinde Mj la Mk, se obține un scor mai mare în comparație cu extinderea Mj la Mi.înlănțuirea semințelor colineare așa cum s-a discutat în au fost utilizate pe scară largă în alinierea secvențelor mari, adică alinierea genomului la genom., Acesta a fost, de asemenea, utilizat pentru a identifica regiunile candidate pentru o citire dat un set de MEMs în BWA. Algoritmii de înlănțuire cu punctaj sunt similari cu algoritmul de programare dinamic pe care l-am propus (DP-MEM). Cu toate acestea, există diferențe care fac DP-MEM potrivit pentru alinierea perechilor de secvențe scurte. DP-MEM ia în considerare faptul că toate MEMs într-o anumită dimensiune decalaj sunt furnizate în intrare și optimizează numărul de iterație în algoritmul. DP-MEM implementează, de asemenea, o abordare euristică pentru a compensa efectul MEMs scurte eliminate din lista de intrare care rezultă lacune între MEMs.,

Optimizare

Având în vedere secvențe de lungime n, algoritmul pentru a extrage MEMs (prevăzute în „MEM extracție” secțiune) necesită 2(n−1) shift și 2n−1 compara operații pe bit-vectori (fiecare act pe $\lceil \frac {n}{32} \rceil $ masina de cuvinte) care duce la un algoritm cu complexitatea O(n2) pentru a produce marginea pic-vectori pentru o anumită pereche de secvențe. Complexitatea funcției care calculează poziționarea MEMs de la marginea bit-vector și Sortează-le pe baza EQ este încă să fie adăugate. Mai mult, dacă m MEMs sunt extrase, Φ2 din algoritmul 2(DP-MEM) are complexitatea lui O (m2)., Având în vedere complexitatea extracției MEM și DP-MEM, mem-Align este mult mai lent decât algoritmii de aliniere disponibili. Pentru a accelera procesul, vă prezentăm mai multe optimizări pentru MEM-Align, care realizează o rulare mai rapidă prin sacrificarea preciziei. Pe de altă parte, introducem metode de îmbunătățire a preciziei cu o pierdere minimă de performanță.alinierea în bandă este o metodă euristică cunoscută pentru a accelera procesul de aliniere. Această tehnică limitează modelul golurilor din aliniere (fișier suplimentar 1: secțiunea VI)., În consecință, dacă alinierea dintre două secvențe nu urmează acest model, algoritmul nu va identifica alinierea. În programarea dinamică tradițională, alinierea se obține după calcularea valorii tuturor celulelor din tabel. Cu toate acestea, cu optimizarea alinierii în bandă, sunt evaluate numai celulele de pe diametru și aproape de diagonală. gl este lățimea benzii în alinierea banded unde celulele mai departe decât GL la diametrul nu sunt evaluate. Celule mai întunecate în Fig. 1 Arată cazul în care gl=1.,spre deosebire de abordarea tradițională de programare dinamică, MEM-Align nu are un tabel similar pentru a aplica alinierea banded. Cu toate acestea, am constatat că putem simula același efect prin limitarea numărului de operații de schimbare în procesul de extracție MEM. De exemplu, dacă schimbăm secvența de interogare până la gl spre dreapta și spre stânga, obținem alinierea bandată cu banda gl. Alinierea în bandă reduce complexitatea extracției MEM de la O(n2) la o(n.(2GL+1)) unde gl este o valoare mică și fixă. Astfel, complexitatea extracției MEM este O (n) atunci când se aplică alinierea în bandă., De asemenea, cu alinierea bandată menționată, este probabil ca mai puține MEMs să fie extrase, ceea ce accelerează pașii algoritmici ulteriori.

îndepărtarea MEM scurt

am observat că majoritatea MEMs extrase sunt scurte și sunt rezultatul de potrivire aleatoriu simboluri. Pentru a accelera alinierea MEM, MEMs mai scurte decât sl sunt filtrate în timpul procesului de extracție MEM. Acest lucru reduce numărul de MEMs care trebuie extrase și prelucrate (accelerarea ulterioară a algoritmului). Filtrarea MEM scurt se face prin extinderea algoritmului 1 cu un set de operații shift și bitwise (fișier suplimentar 1: secțiunea VII).,pe de altă parte, există cazuri rare în care MEMs scurte fac parte din aliniere; de exemplu, un simbol de potrivire înconjurat de nepotriviri. Fără a avea toate MEMs în lista de intrare, DP-MEM nu este capabil să găsească aceeași aliniere ca se găsește atunci când toate MEMs există în lista de intrare. Pentru a compensa MEMs scurte pierdute în intrare, modificăm Φ2H din DP-MEM pentru a lua în considerare posibilitatea de a avea meciuri scurte între MEMs (fișier suplimentar 1: secțiunea VIII).

pot exista cazuri mai dificile în care, în aliniere, există mai multe MEMs scurte între două MEMs (vezi Fig. 9)., Singura modalitate de a identifica corect scorul pentru zona dintre Mi și Mj în Φ2H este aplicarea unei alinieri globale la această regiune. Cu toate acestea, Φ2H este o operație frecventă și ar trebui să rămână rapidă. În consecință, am decis să depășim parțial problema prin limitarea cazurilor posibile (o metodă euristică).

Dacă există lacune în zona cuprinsă între Mi și Mj, presupunem că există doar un continuu decalaj fie la stânga sau la capătul din dreapta al zonei. Apoi, doar două alinieri sunt posibile pentru zonă., Numărul de simboluri de potrivire este contorizat pentru ambele cazuri într-o manieră secvențială, iar cel care are ca rezultat potriviri maxime este luat ca număr de potriviri între Mi și Mj (fișier suplimentar 1: Secțiunea IX). Comparația secvențială este o operație costisitoare și concepem o metodă pentru a evita comparația secvențială atunci când este posibil (fișier suplimentar 1: secțiunea X).orice alt caz care nu se potrivește presupunerii de mai sus are ca rezultat o aliniere cu un scor mai mic., Cu toate acestea, având în vedere rata scăzută a lacunelor și a neconcordanțelor, posibilitatea de a avea mai multe lacune și neconcordanțe într-o zonă mică este scăzută.în Φ2 din DP-MEM, MJ extinde toate aliniamentele care se termină în {M1 … MJ-1} (dacă este posibil). Cu toate acestea, pentru fiecare Mj acolo este un mic subset Ωj⊆{M1,…Mj−1} astfel încât prin extinderea Mj pentru toate aliniamente se încheie într-o Mi∈Ωj alinierea care se termină în Mj este găsit (Eq. 2). Cu alte cuvinte, ar fi mai puține $s_{i}^{j}$ care trebuie evaluate. Definiția setului Ωj și dovada Eq. 2 sunt furnizate în fișierul suplimentar 1: secțiunea XI., Definiția Ωj este afectată atunci când se aplică optimizarea de eliminare MEM scurt (fișier suplimentar 1: secțiunea XII).

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

(2)

Hibrid de aliniere

Pentru a menține acuratețea algoritmului, ne-am decis să utilizăm o metodă hibridă, care este o combinație de MEM-Alinierea și Smith-Waterman algoritm. Definim trei cazuri în care alinierea MEM poate fi inexactă., Dacă alinierea unei perechi de secvențe cade într-unul din aceste cazuri, folosim algoritmul Smith-Waterman pentru a alinia secvențele. Aceste cazuri sunt:

când secvențele sunt repetitive, iar numărul de MEMs extrase depășește pragul TM. Am constatat că MEM-Align este probabil să producă aliniere inexactă la alinierea secvențelor repetitive. O valoare TM adecvată scade șansa de a raporta alinierea inexactă cu o creștere neglijabilă a timpului mediu de procesare.,
când scorul de aliniere calculat pentru alinierea generată de mem-Align este mai mic decât un prag TS. Acest caz apare mai ales atunci când există un decalaj în alinierea care nu pot fi identificate din cauza alinierii banded.
atunci când nu MEM mai mult decât sl există pentru a fi extras (un caz rar)., Dacă SL este setat la o valoare ridicată și similitudinea dintre secvențe este scăzută,

deși trimiterea perechilor de secvențe la un algoritm extern are ca rezultat un calcul suplimentar, numărul secvențelor trimise algoritmului extern rămâne mic dacă valorile corespunzătoare sunt alese pentru TM și TS.când distanța dintre Mi și Mj este mare, nu este probabil să aibă mi și Mj ca MEMs adiacente în aliniere., Prin urmare, algoritmul ignoră extensia dacă distanța dintre Mi și Mj este mai lungă decât un prag TD (reducând în continuare numărul de $s_{i}^{j}$ care trebuie evaluat). Această optimizare îmbunătățește ușor performanța cu un efect secundar neglijabil asupra preciziei.

Be Able