Articles

parvis inriktning av nukleotidsekvenser med maximal exakt matchning

oktober 15, 2020 by admin

tillvägagångssätt

i vår föreslagna algoritm är det första steget mot att anpassa sekvenser att extrahera MEMs mellan sekvenser genom att direkt jämföra dem. Figur 3a är ett exempel som jämför ett mål och en frågesekvens där CTC och AAA är två MEMs som identifieras genom jämförelsen. Varje grupp av kontinuerliga identiska symboler i jämförelsen resulterar i en MEM även om den består av endast en enda matchande symbol., För att extrahera alla MEMs mellan sekvenserna måste frågesekvensen skiftas hela vägen till höger och till vänster en symbol åt gången (se Fig. 3b). Efter varje skift måste jämförelsesteget upprepas för att identifiera nya MEMs. Till exempel, den tredje raden i Fig. 3B representerar fallet där frågesekvensen skiftas till höger en symbol och jämförs med målsekvensen. Resultatet av jämförelsen identifierar AAAAGC som en ny MEM. Alla andra MEMs som extraheras genom skift-och jämförelseoperationer markeras också i Fig. 3b., Tre av MEMs (Mx,My och MZ) är markerade med olika färger som ska användas för senare förklaring.

i affine-gap-poängmodellen beräknas justeringspoängen enligt Eq., 1 där Nm är antalet matcher varje som får en match poäng av Rm, Nx är antalet missmatchningar varje emot en obalans straff Px, No är antalet gap öppningar varje emot ett gap Öppet straff av Po och Ng är total längd av alla luckor, varje gap som får en lucka förlängning straff av Pg. Det skulle finnas en lucka öppning för varje grupp av kontinuerlig lucka. Om det till exempel finns två luckor i inriktningen, där längden på det första gapet är tre och längden på det andra gapet är fyra, så finns det två gapöppningar (no=2) och den totala längden på gapet är sju (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)

Med tanke på listan över alla MEMs kan justeringen beräknas med hjälp av partiella anpassningar. Tänk till exempel MEMs Mx, My och MZ i Fig. 3B. de partiella inriktningar som görs genom att ta olika kombinationer av Mx, My och MZ tillsammans med antalet matcher, mismatch och luckor, samt de resulterande inriktnings poäng visas i Fig. 4. Justeringen som bara innehåller Mx och MZ resulterar i högsta justeringspoäng., Observera att MY och Mz överlappar varandra och när båda betraktas i samma inriktning är överlappningen utesluten från Mz. Med tanke på alla MEMs i Fig. 3B resulterar i många fler kombinationer där ingen av dem uppnår en högre poäng.

att undersöka alla möjliga kombinationer av MEMs skulle vara uttömmande., I avsnittet ”Inriktningsalgoritm” beskriver vi en ny dynamisk programmeringsalgoritm DP-MEM som effektivt hittar den bästa kombinationen utan att överväga alla fall. DP-MEM behöver veta vilka delar av sekvenserna som matchar men inte de faktiska symbolerna i sekvenserna. Inmatningen till DP-MEM är placeringen av MEMs i målet och i frågesekvenserna som erhålls under mem-extraktionsprocessen som beskrivs i avsnittet ”mem-extraktion”., Hur MEMs representeras med sina positioner och hur antalet matcher, missmatchningar och luckor beräknas när MEMs kombineras i en anpassning förklaras i resten av detta avsnitt. Figur 5 är ett annat exempel justering med sex MEMs (M1 till M6) som bildar anpassningen mellan målsekvens T och frågesekvens Q. för enkelhetens skull finns det ingen överlappning mellan MEMs i detta exempel. Varje mem Mi representeras som en triplett av heltal: startpositionerna i T och i Q (STi och SQi respektive) och dess längd (Li)., Slutpositionen i T och Q kan beräknas senare (Φ2E av Algoritm 2). Tabell 1 visar längden och placeringen av M1 till M6 i T och i Q.

Tabell 1 Start-och slutposition för MEMs i Fig., 5

Tabell 2 Beräkning antal felmatchningar och luckor mellan MEMs i Fig. 5

om det finns både missmatchningar och luckor mellan Mi och MJ anses alla luckor vara kontinuerliga för att minska gapet Öppet straff (endast ett gap Öppet straff tillämpas för en fortsätter gap). Således tillämpas endast ett gap Öppet straff för alla intilliggande MEMs som har luckor mellan dem., Placeringen av missmatchningar och det enda kontinuerliga gapet är inte viktigt, eftersom det inte skulle påverka inriktningsresultatet. Vi antar att mismatch straff är konstant (detta är vanligt för DNA-sekvenser).

mem extraction

det finns metoder för att extrahera maximala exakta matchningar mellan långa sekvenser som ett helt genom. Dessa metoder är dock baserade på förbehandling och indexering av en eller båda sekvenser som är en tidskrävande operation. Till exempel, i DNA läs aligner, referensgenomet indexeras en gång, och samma index används varje gång en ny läsning är inriktad., Vi letar efter en snabb algoritm för att identifiera MEMs mellan relativt korta sekvenser som ändras för varje inriktning. En brute force-metod för detta problem (ytterligare fil 1: avsnitt II) är långsam och ineffektiv(med komplexiteten hos O (n3)). Vi föreslår en snabb bit-nivå parallell metod för att påskynda mem extraktionsprocessen. Vår mem-extraktionsmetod är baserad på skiftet och jämför operationer som visas i Fig. 3b. det första steget är att representera sekvenser med bitvektorer, där A, C, T och G kodas som 00, 01, 10 respektive 11 (Ytterligare fil 1: avsnitt III)., Figur 6 illustrerar ett exempelsekvens par, tillsammans med motsvarande bit-vektor representationer. I en råvarudator är maskinordet vanligtvis 64 bitar som rymmer 32 nukleotidsymboler. Eftersom en sekvens vanligtvis är större än 32 symboler lagras motsvarande bitvektor i flera maskinord. Varje operation på bitvektorer av sekvenser av storlek N-symboler verkar på $\lceil \frac {n}{32} \ rceil$ maskinord.

med bitvektorer representation av sekvenser, flytta en sekvens med en symbol är densamma som att flytta bitvektorn med två bitar, och jämföra sekvenser kan göras med xor instruktion (32 symboler åt gången). I xor-utgången (X) betyder 00 Att symboler matchas, och en sekvens av 00s visar en MEM., En uppsättning skift-och bitwise-operationer som visas i algoritm 1 beräknar X och därefter edge bit-vector (E) där start och slutet av varje MEM markeras med inställda bitar (bitar med ett värde av en). Figur 6 visar X-och e-bitvektorerna med markerade MEMs. Placeringen av MEMs i sekvenser beräknas sedan från kantbitvektorn (ytterligare fil 1: avsnitt IV).,

Justeringsalgoritm

i ”Approach” – avsnittet visar vi att genom att överväga olika kombinationer av MEMs och beräkna justeringsresultatet för motsvarande justering kan man identifiera kombinationen av MEMs som resulterar i den maximala justeringsresultatet. Att undersöka alla möjliga kombinationer av MEMs är dock en naiv lösning. Ett mer systematiskt sätt att hitta inriktningen effektivt är att använda dynamisk programmering.

dynamisk programmering är metoden att närma sig lösningen till ett problem genom att definiera och lösa mindre underproblem., Lösningar för underproblem används för att lösa ett större problem vid varje steg. Processen upprepas tills alla underproblem är lösta. Så småningom skulle lösningen på ett av underproblemen vara lösningen på det ursprungliga problemet. När alla underproblem löses en backtracking process identifierar en rad lösningar som bidrar till den slutliga lösningen. Vid dynamisk programmering bör det finnas en beställning av indata längs vilka rekursionsproceduren fortsätter.

vi sorterar alla MEMs enligt positionen för deras slut i query sequence (EQ)., MEMs som slutar i samma position beställs på ett godtyckligt sätt. Underproblemet är att hitta anpassningen av efterföljanden av T och Q som slutar vid JTH MEM MJ (t respektive Q). Vi kommer att visa att denna beställning av MEM är tillräcklig för att stödja rätt rekursion.

i den sorterade listan över MEMs indikerar EQi=EQJ att en av Mi eller mj helt överlappar den andra MEM i frågesekvensen. Eftersom i Φ2B av algoritm 2 är överlappningsområdet uteslutet, kan Mi och mj inte vara i samma inriktning., Således löses ith-och jth-delproblemen oberoende av varandra och ordningen för i och j i den sorterade listan kunde vara godtycklig. Om EQK>EQJ (k>j i den sorterade listan) kunde Mk inte vara en del av justeringen som slutar i mj. Därmed jth problem kan lösas oberoende av lösningen på kth subproblem. Observera att det också är möjligt att sortera MEMs baserat på deras slutposition i målsekvensen (ET) med hjälp av en liknande motivering.

vår föreslagna dynamiska programmeringsalgoritm (DP-MEM) utarbetas i algoritm 2., Till exempel MEMs extraheras i Fig. 3B, den dynamiska programmeringstabellen och mellanliggande värde som beräknas i algoritmen visas i fikon. 7 respektive 8. Inmatningen till DP-MEM är listan över MEMs där varje MEM (Mj) är en triplett av heltal . Den andra ingången n är antalet MEMs i listan. Utgången S är justeringsresultatet för sekvenserna. Algoritmen skriver ut indexen för alla MEMs som bildar inriktningen där de första och sista tryckta siffrorna är indexen för de högsta och de vänstra Memsna i justeringen respektive., Alla steg i algoritm 2 kommenteras i följande:

Φ1: Scoring varje MEM för alla dess matchande symboler., Observera att det finns lj matchande symboler i MJ. Sj representerar den högsta justeringspoängen för justeringen som slutar vid mj. Initiera Sj i detta steg liknar att beräkna partiella inriktningen poäng när endast Mj ingår i inriktningen. W används för backtracking. Värdet på -1 indikerar att den nuvarande Sj erhålls genom att överväga MJ ensam i anpassningen.
Φ2: Computing Sj för varje MEM (Mj)., För att beräkna Sj, för varje MEM Mi där Mi visas före Mj i listan, lägger algoritmen till MJ till justeringen slutar vid Mi (förlängning tidigare funna anpassningar) och letar efter förlängningen som maximerar Sj (lösa ett större subproblem med tidigare lösta subproblem).
Φ2A: hoppa över förlängning när det inte är möjligt. Om eti>ETj innehåller Mi en del av målsekvensen som ligger utanför inriktningen som slutar vid MJ och förlängningen är inte möjlig. Om EQi=EQj eller ETi=ETj eller SQi≥SQj eller STi≥STj då en av MEMs helt överlappar andra MEM., I detta fall kan Mi och mj inte vara i en anpassning tillsammans.
Φ2B: beräkna längden på överlappningen mellan Mi och mj. Om\ ({MO}_{I}^{j}\) är mindre än eller lika med noll, finns ingen överlappning.
Φ2C: att behålla en kopia av Mj innan överlappning utesluts.
Φ2D: om överlappning föreligger, exklusive överlappning från MJ
Φ2E: Computing slutposition av MJ i T och Q.
Φ2F: Computing avståndet (antal symboler) mellan Mi och Mj i T och Q.
Φ2G: Computing antal missmatcher och luckor mellan Mi och Mj.,
Φ2H: beräkna straffet för missmatchningar och luckor mellan Mi och MJ ($P_{i}^{j}$). Om gapet finns, subtraheras endast ett gap Öppet straff.
Φ2I: Computing alignment score $\left (s_{i}^{j}\right)$ när MJ läggs till justeringen slutar vid Mi. Poängen för alla matchande symboler i MJ (lj×Rm) läggs till justeringsresultatet för justeringsändelsen vid Mi (Si). Därefter subtraheras straffet för luckorna och missmatchningarna mellan Mi och MJ$\left (p_{i}^{j}\right)$.,
Φ2J: om du utvidgar Mj till justeringen som slutar med Mi-resultat till en poäng $\left (s_{i}^{j}\right)$ högre än nuvarande poäng för MJ (sj) lagras den nya poängen i sj. Även W är inställd på att jag ska hålla reda på Mi som maximerar poängen för MJ.
Φ2K: återställa värdet för MJ före uteslutning så att Mj kan användas i andra inriktningsförlängningar.
Φ3: Letar efter MEM med högsta sj. Detta MEM är den sista MEM i anpassningen (mig)., Den högsta poängen (Se) returneras som S vilket är den högsta justeringspoängen för de givna sekvenserna. Indexet för MEM som maximerar Sj lagras i e för att börja backtracking från mig.
Φ4: i justeringen är den omedelbara föregående MEM för mig den som maximerar justeringspoängen för mig. Indexet för sådan MEM lagras i W. som ett resultat besöker iterationen av F W W indexet för alla MEMs i inriktningen. När W är lika med -1 är Mf den första MEM i inriktningen och iterationen stoppas.,

i vår algoritm straffar vi inte missmatchningar och luckor före den första MEM och efter den sista MEM i justeringen. Detta resulterar i en lokal inriktningsalgoritm. Genom att överväga dessa påföljder genererar algoritmen en global anpassning (ytterligare fil 1: avsnitt V).

ekvationen för att beräkna $P_{i}^{j}$ i Φ2H av algoritm 2 förutsätter att det inte finns någon matchande symbol mellan T och Q i området mellan Mi och MJ (alla symboler räknas som felaktigheter eller luckor)., Även om detta antagande inte är sant för alla Mi, är det alltid sant för Mi som leder till maximal $s_{i}^{j}$ som åsidosätter effekten av antagandet är felaktig för andra Mi. Som ett bevis, anta att det finns en matchande symbol i området mellan Mi och mj. Matchande symbol skulle vara en MEM (Mk). Mk förlängs redan till inriktningen som slutar vid Mi. Således uppnås en högre poäng vid förlängning av MJ till Mk jämfört med att förlänga MJ till Mi.

kedja colinear frön som diskuterats i har använts i stor utsträckning i anpassningen av stora sekvenser, dvs genom-till-genom justering., Det har också använts för att identifiera kandidatregioner för en läsning ges en uppsättning MEMs i BWA. Kedjningsalgoritmer med poäng liknar den dynamiska programmeringsalgoritmen vi föreslog (DP-MEM). Det finns dock skillnader som gör DP-MEM lämplig för parvis anpassning av korta sekvenser. DP-MEM tar hänsyn till att alla MEMs inom en viss gapstorlek tillhandahålls i ingången och optimerar antalet iterationer i algoritmen. DP-MEM genomför också ett heuristiskt tillvägagångssätt för att kompensera för effekten av korta MEMs som tagits bort från inmatningslistan, vilket resulterar i luckor mellan MEMs.,

optimering

givna sekvenser av längd n, algoritmen för att extrahera MEMs (tillhandahålls i avsnittet ”mem extraction”) kräver 2(n−1) shift och 2n−1 Jämför operationer på bitvektorer (varje handling på $\lceil \frac {n}{32} \rceil $ maskinord) som resulterar i en algoritm med komplexitet av O(n2) för att producera kantbitvektorer för det givna paret av sekvenser. Komplexiteten i funktionen som beräknar positionering av MEMs från edge bit-vektorn och sorterar dem baserat på EQ är ännu inte tillagd. Vidare, om m MEMs utvinns, Φ2 av Algoritm 2 (DP-MEM) har komplexitet O(m2)., Med tanke på komplexiteten i mem-extraktion och DP-MEM är mem-Align mycket långsammare än tillgängliga inriktningsalgoritmer. För att påskynda processen presenterar vi flera optimeringar för MEM-Align som uppnår snabbare runtime genom att offra noggrannhet. Å andra sidan introducerar vi metoder för att förbättra noggrannheten med minimal prestandaförlust.

banded alignment

Banded alignment är en känd heuristisk metod för att påskynda anpassningsprocessen. Denna teknik begränsar mönstret för luckorna i justeringen (ytterligare fil 1: avsnitt VI)., Följaktligen, om anpassningen mellan två sekvenser inte följer detta mönster, identifierar algoritmen inte inriktningen. I traditionell dynamisk programmering erhålls inriktningen efter beräkning av värdet av alla celler i tabellen. Med den bandade justeringsoptimeringen utvärderas dock endast celler på diametern och nära diagonalen. gl är bandets bredd i bandad inriktning där celler längre än gl till diametern inte utvärderas. Mörkare celler i Fig. 1 Visa fallet där gl=1.,

Till skillnad från traditionell dynamisk programmeringsmetod har MEM-Align inte en liknande tabell för att tillämpa banded alignment. Vi fann dock att vi kan simulera samma effekt genom att begränsa antalet skiftoperationer i mem-utvinningsprocessen. Till exempel, om vi flyttar frågesekvensen upp till gl till höger och till vänster, uppnår vi banded inriktning med bandet av gl. Banded-justering minskar komplexiteten hos mem-extraktion från O(N2) till O(n.(2gl+1)) där gl är ett litet och fast värde. Således är komplexiteten hos mem-extraktionen O (n) när bandad inriktning appliceras., Med den nämnda banded inriktningen är det också troligt att färre MEMs extraheras vilket påskyndar de efterföljande algoritmiska stegen.

kort mem-borttagning

Vi observerade att majoriteten av extraherade MEMs är korta och är resultatet av slumpmässigt matchande symboler. För att påskynda mem-Align filtreras MEMs kortare än sl ut under mem-extraktionsprocessen. Detta minskar antalet MEMs som ska extraheras och bearbetas (därefter påskyndar algoritmen). Filtrering kort MEM görs genom att förlänga algoritm 1 med en uppsättning skift och bitwise operationer (ytterligare fil 1: avsnitt VII).,

å andra sidan finns det sällsynta fall där korta MEMs är en del av justeringen; till exempel en matchande symbol omgiven av missmatchningar. Utan att ha alla MEMs i inmatningslistan kan DP-MEM inte hitta samma justering som den hittar när alla MEMs finns i inmatningslistan. För att kompensera för förlorade korta MEMs i ingången modifierar vi Φ2H av DP-MEM för att överväga möjligheten att ha korta matcher mellan MEMs (ytterligare fil 1: avsnitt VIII).

det kan finnas svårare fall där det i justeringen finns flera korta MEMs mellan två MEMs (se Fig. 9)., Det enda sättet att korrekt identifiera betyg för området mellan Mi och Mj i Φ2H är att tillämpa en global anpassning till denna region. Φ2H är dock en frekvent operation och bör förbli snabb. Följaktligen bestämde vi oss för att delvis övervinna problemet genom att begränsa möjliga fall (en heuristisk metod).

om det finns luckor i området mellan Mi och MJ, antar vi att det bara finns ett kontinuerligt gap antingen till vänster eller till vänster.den högra änden av området. Då är det bara två anpassningar som är möjliga för området., Antalet matchande symboler räknas för båda fallen på ett sekventiellt sätt och det som resulterar i maximala matchningar tas som antalet matchningar mellan Mi och MJ (ytterligare fil 1: avsnitt IX). Den sekventiella jämförelsen är en dyr operation och vi utarbetar en metod för att undvika Sekventiell jämförelse när det är möjligt (ytterligare fil 1: Avsnitt X).

alla andra fall som inte passar ovanstående antagande resulterar i en justering med en lägre poäng., Med tanke på den låga graden av luckor och missmatchningar är dock möjligheten att ha flera luckor och missmatchningar i ett litet område låg.

efficient alignment extension

I Φ2 av DP-MEM utökar MJ alla anpassningar som slutar i {M1…MJ−1} (om möjligt). Men för varje Mj finns en mindre delmängd Ωj{M1…MJ−1} så att genom att förlänga Mj till alla anpassningar som slutar i en mi Ωj hittas justeringen som slutar i Mj (Eq. 2). Med andra ord skulle det finnas färre $s_{i}^{j}$ som ska utvärderas. Definition av uppsättningen Ωj och bevis på Eq. 2 finns i ytterligare fil 1: avsnitt XI., Definitionen av Ωj påverkas när kort mem borttagning optimering tillämpas (ytterligare fil 1: avsnitt XII).

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

(2)

Hybridinriktning

för att upprätthålla algoritmens noggrannhet bestämde vi oss för att använda en hybridinriktning.metod som är en kombination av mem-align och Smith-waterman algoritm. Vi definierar tre fall där mem-Align kan vara felaktiga., Om anpassningen av ett par sekvenser faller ner i ett av dessa fall använder vi Smith-Waterman-algoritmen för att anpassa sekvenser. Dessa fall är:

När sekvenserna är repetitiva och antalet extraherade MEMs överstiger tröskelvärdet TM. Vi hittade MEM-Align kommer sannolikt att producera felaktig inriktning när man anpassar repetitiva sekvenser. Ett lämpligt TM-värde minskar risken för att rapportera felaktig anpassning med en försumbar ökning av den genomsnittliga behandlingstiden.,
När den beräknade justeringspoängen för justeringen som genereras av MEM-Align är lägre än en tröskel TS. Detta fall uppstår oftast när det finns ett gap i anpassningen som inte kan identifieras på grund av bandad inriktning.
När ingen MEM längre än sl existerar för att extraheras (ett sällsynt fall)., Om sl är inställt på ett högt värde och likheten mellan sekvenserna är låg,

även om sändning av sekvenspar till en extern algoritm resulterar i ytterligare beräkning, är antalet sekvenser som skickas till den externa algoritmen fortfarande liten om lämpliga värden väljs för TM och TS.

hoppa över avlägsna MEMs

När avståndet mellan Mi och Mj är stort, är det inte troligt att ha Mi och MJ som intilliggande MEMs i justeringen., Därför hoppar algoritmen över förlängningen om avståndet mellan Mi och Mj är längre än en tröskel TD (ytterligare minskning av antalet $s_{i}^{j}$ som ska utvärderas). Denna optimering förbättrar prestandan något med en försumbar bieffekt på noggrannheten.

Be Able