Articles

Parvis justering av nukleotide-sekvenser ved hjelp av maksimal nøyaktig matcher

Tilnærming

I våre forslag til algoritme, det første skrittet mot å samkjøre sekvenser er å trekke ut MEMs mellom sekvenser ved direkte sammenligning av dem. Figur 3a er et eksempel som sammenligner et mål og en spørring sekvens der CTC og AAA er to MEMs identifisert av sammenligningen. Hver gruppe av kontinuerlig identiske symboler i sammenligningen, resultere i et MEM selv om den består av bare en enkelt matchende symboler., For å trekke ut alle MEMs mellom sekvenser, spørring-sekvensen må bli flyttet helt til høyre og til venstre et symbol på en gang (se Fig. 3b). Etter hvert skift, sammenligningen trinn må gjentas for å identifisere nye MEMs. For eksempel, den tredje linje i Fig. 3b representerer de tilfelle hvor spørring sekvensen er forskjøvet til høyre ett symbol og er sammenlignet med målsekvensen. Resultatet av sammenligningen identifiserer AAAAGC som en ny MEM. Alle andre MEMs hentet ut av skift og sammenlign operasjoner er også fremhevet i Fig. 3b., Tre av MEMs – (Mx,My og Mz) er uthevet med forskjellige farger som skal brukes for senere forklaring.

Fig. 3

MEM utvinning ved hjelp av skift og sammenlign operasjoner. en Identifisere MEMs ved direkte sammenligning av sekvenser. b spørringen er forskjøvet til venstre helt til siste symbolet i søket sekvens er justert i forhold til det første symbolet i målsekvensen. Deretter spørringen sekvensen er forskjøvet til høyre til det første symbolet i søket sekvens er justert i forhold til det siste symbolet av målsekvensen., Etter hvert skift den overlappende delen av spørringen, og målet sekvenser er i forhold til å identifisere nye MEMs. Tre av MEMs – (Mx,My og Mz) er uthevet med forskjellige farger som skal brukes for senere forklaring

I den affine-gap scoring modell, justering score SOM er beregnet ved hjelp av Eq., 1 hvor Nm er antall kamper hver motta en match score av Rm,Nx er antallet samsvarer ikke hver motta en mismatch straff av Px,Ingen er antall gap åpninger hvert mottak av en sprekk åpne straff av Po og Ng er total lengde av alle åpninger, hver gap motta et gap extension straff av Pg. Det ville være et gap åpner for hver gruppe av kontinuerlig gap. For eksempel, hvis det er to hull i justeringen, der lengden av den første gapet er tre og lengden av den andre gapet er fire, så det er to gap åpninger (Nei=2) og den totale lengden av gapet er syv (Ng=3+4=7).,

$$ {}SOM = (N_{m} \ganger R_{m}) – ((N_{x} \ganger P_{x}) + (N_{o} \ganger P_{o}) + (N_{g} \ganger P_{g})) $$
(1)

Gitt liste av alle MEMs, justeringen kan beregnes ved hjelp av delvis justeringer. For eksempel, tenk MEMs Mx, My og Mz i Fig. 3b. Delvis justeringer gjort ved å ta ulike kombinasjoner av Mx, My og Mz sammen med antall kamper, og samsvarer ikke hull, så vel som den resulterende justering score er vist i Fig. 4. Justeringen som bare omfatter Mx og Mz resultater i høyeste justering score., Vær oppmerksom på at Min og Mz overlapper hverandre, og når begge er vurdert i samme justering overlapper hverandre, er utelukket fra Mz. Tatt i betraktning alle MEMs i Fig. 3b resultater i mange kombinasjoner der ingen av dem oppnår en høyere score.

Fig. 4

Alle mulige kombinasjoner av MEMs i justeringen

å Undersøke alle mulige kombinasjoner av MEMs ville være uttømmende., I «Alignment algoritmen» – delen beskriver vi en roman dynamisk programmering algoritme DP-MEM som effektivt finner den beste kombinasjonen uten å vurdere alle tilfeller. DP-MEM behov for å vite hvilke deler av sekvenser kampen, men ikke selve symbolene i sekvenser. Innspill til DP-MEM er lokalisering av MEMs i mål og i spørringen sekvenser som er innhentet i løpet av MEM utvinning prosessen beskrevet i «MEM utvinning» – delen., Hvordan MEMs er representert med sine posisjoner, og hvordan antall kamper, og samsvarer ikke hull beregnes når MEMs er kombinert i en justering er forklart nærmere i resten av dette avsnittet. Figur 5 er et annet eksempel er i samsvar med seks MEMs – (M1 til M6) som danner justeringen mellom målsekvens T og spørring sekvens Q. For enkelhets det er ingen overlapp mellom MEMs i dette eksemplet. Hver MEM Mi er representert som en triplett av heltall numre: start posisjon i T og Q (STi og SQi henholdsvis) og lengde (Li)., Avslutningen posisjoner i T og Q kan bli beregnet senere (Φ2E av Algoritmen 2). Tabell 1 viser lengde og plassering av M1 til M6 i T og Q.

Fig. 5

Et eksempel justering med uthevet MEM

Tabell 1 Start-og slutt posisjon av MEMs i Fig., 5
Tabell 2 Computing antallet samsvarer ikke og avvik mellom MEMs i Fig. 5

I tilfelle det er både samsvarer ikke og avvik mellom Mi og Mj, alle hullene anses som kontinuerlig for å redusere gapet åpne straff (bare ett gap åpne straff er søkt om en fortsetter gap). Derfor, for alle tilstøtende MEMs som har mellomrom mellom dem, bare ett gap åpne straff er brukt., Plassering av og samsvarer ikke den eneste kontinuerlig gapet er ikke viktig, som det ikke ville påvirke justering score. Vi antar at misforholdet straff er konstant (dette er vanlig for DNA-sekvenser).

MEM utvinning

Det finnes metoder for å trekke ut maksimal eksakte samsvar mellom lengre sekvenser for eksempel et helt genom. Men disse metodene er basert på forbehandling og indeksering av en eller begge sekvenser som er en tidkrevende operasjon. For eksempel, i DNA les aligner, referanse-genom er indeksert gang, og den samme indeksen er brukt hver gang en ny lese er justert., Vi er ute etter en rask algoritme for å identifisere MEMs mellom relativt korte sekvenser som endres for hver justering. Et brute force-metoden for dette problemet (Ekstra fil 1: Del II) er trege og ineffektive (med kompleksitet O(n3)). Vi foreslår en rask bit-nivå parallell metode for å øke hastigheten på MEM utvinning prosessen. Våre MEM utvinning metode er basert på skift og sammenlign operasjoner vist i Fig. 3b. Det første trinnet er å representere sekvenser med bit-vektorer, der A, C, T og G, er kodet som 00, 01, 10 og 11, henholdsvis (Ekstra fil 1: Del III)., Figur 6 viser et eksempel sekvensen par, sammen med tilsvarende bit-vektor representasjoner. I en vare datamaskinen, maskinen ord er vanligvis 64 bits som kan romme 32 nukleotid-symboler. Siden en sekvens er vanligvis større enn 32-symboler, tilsvarende bit-vector er lagret i flere maskinen ord. Hver operasjon på bit-vektorer av sekvenser av størrelse n-symboler fungerer på \(\lceil \frac {n}{32} \rceil \) maskinen ord.

Fig., 6

gjengivelse av sekvenser med bit-vektorer. XOR-utgang (X) med uthevet MEMs. Kantene bit-vektor (E) identifiserer starten og slutten av hver MEMs

Med bit-vektorer representasjon av sekvenser, skiftende en sekvens av ett symbolet er det samme som å flytte bit-vektor av to biter, og sammenligne sekvenser kan gjøres med XOR undervisning (32-symboler på en gang). I XOR-utgang (X), 00 betyr at symboler er matchet, og en sekvens av 00s viser et MEM., Et sett av skift og bitvis operasjoner som vist i Algoritme 1 regner ut X og senere kanten bit-vektor (E) som i starten og slutten av hver MEM er merket med sett bits bits (med en verdi på én). Figur 6 viser X og E bit-vektorer med uthevet MEMs. Plasseringen av MEMs i sekvenser er deretter beregnet fra kanten bit-vektor (Ekstra fil 1: Kapittel IV).,

Justering algoritme

I «Tilnærming» – delen, viser vi at ved å vurdere ulike kombinasjoner av MEMs og computing justeringen score for tilsvarende justering, kan man identifisere kombinasjon av MEMs som resulterer i maksimal justering score. Imidlertid, undersøke alle mulige kombinasjoner av MEMs er en naiv løsning. En mer systematisk måte å finne justeringen er effektivt å bruke dynamisk programmering.

Dynamisk programmering er metoden nærmer seg løsning på et problem ved å definere og løse mindre subproblems., Løsninger for subproblems brukes til å løse et større problem på hvert trinn. Prosessen gjentas til alle subproblems er løst. Til slutt, løsningen til en av de subproblems ville være løsningen på det første problemet. Når alle subproblems er løst på en tilbakesporing prosessen identifiserer en rekke løsninger, som bidrar til den endelige løsningen. I dynamisk programmering, bør det være en bestilling av input data sammen som recursion prosedyren fortsetter.

Vi sortere alle MEMs i henhold til den posisjonen til sin ende i søket sekvens (EQ)., MEMs som ender i samme posisjon er ordnet på en vilkårlig måte. Den jth subproblem er å finne plasseringen av undersekvenser av T og Q som slutten på jth MEM Mj (T og Q henholdsvis). Vi vil vise at dette bestilling av MEM er tilstrekkelig til å støtte den rette recursion.

I den sorterte listen av MEMs, EQi=EQj indikerer at ett av Mi eller Mj fullt overlapper den andre MEM i søket rekkefølge. Siden i Φ2B av Algoritmen 2 overlapper regionen er utelukket, Mi og Mj-kan ikke være i samme stilling., Dermed ed og jth subproblems løses uavhengig av hverandre og rekkefølgen av i og j i den sorterte listen kan være vilkårlig. Hvis EQk>EQj (k>j i sortert liste), Mk kunne ikke være en del av den justering som ender i Mj. Dermed jth subproblems kan løses uavhengig fra løsning til kth subproblem. Vær oppmerksom på at det er også mulig å sortere MEMs basert på deres posisjon i forhold målsekvensen (ET) ved hjelp av en lignende begrunnelse.

Våre forslag dynamisk programmering-algoritmen (DP-MEM) er utarbeidet i Algoritme 2., For eksempel MEMs ut i Fig. 3b, dynamisk programmering bord og middels verdi beregnet i algoritmen er vist i Fig. 7 og 8, henholdsvis. Innspill til DP-MEM er listen over MEMs hvor hver MEM (Mj) er en triplett av heltall . Den andre input n er antall MEMs i listen. Produksjonen S er justeringen score for sekvenser. Algoritmen skriver ut indekser for alle MEMs som danner justering hvor den første og den siste trykte tall er indekser av høyre og venstre MEMs i justeringen henholdsvis., Alle trinnene i Algoritmen 2 er kommentert i det følgende:

Fig. 7

Dynamisk programmering tabellen brukes i Algoritmen 2 til å behandle hentet MEMs i Fig. 3b. Celle i og j representerer verdien av \(S_{i}^{j}\). Tomme celler er ikke evaluert i Φ2. Evaluering av celler med kryss hoppet i Φ2A. Første Verdi av Sj er beregnet i Φ1. Endelige verdien av Sj og dens kilde (hva maksimerer Sj) er markert for hver rad. Den høyeste Sj (S13) er justeringen score., M13 er den siste MEM i alignment og MEM før det er MW=M3. Siden W=-1,M3 er den første MEM i justeringen. Poengsystemet for denne justeringen er Rm=2,Px=3,Po=4 og Pe=1

Fig. 8

Mellomliggende verdier for å beregne \(S_{i}^{j}\) i Fig. 7. Vær oppmerksom på at Sij i dette tallet refererer til \(S_{i}^{j}\)

  • Φ1: Scoring hver MEM for alle sine matchende symboler., Merk at det er Lj matchende symboler i Mj. Sj representerer den høyeste justering score for justering slutter på Mj. Initialiserer Sj i dette trinnet er lik computing delvis justering score når bare Mj er inkludert i justeringen. W er brukt for å spore tilbake. Verdien -1 betyr at den gjeldende Sj er oppnådd ved å vurdere Mj alene i justeringen.

  • Φ2: Databehandling Sj for hver MEM (Mj)., For å beregne Sj, for hver MEM Mi der Mi vises før Mj i listen, algoritmen legger Mj til justering slutter på Mi (som strekker seg tidligere funnet justeringer) og ser for filtypen som maksimerer Sj (løse en større subproblem ved hjelp av tidligere løst subproblems).

  • Φ2A: Hopp extension når det ikke er mulig. Hvis ETi>ETj deretter Mi inneholder en del av mål-sekvens som er utenfor justering slutter på Mj og utvidelse er ikke mulig. Hvis EQi=EQj eller ETi=ETj eller SQi≥SQj eller STi≥STj deretter ett av MEMs fullt overlapper den andre MEM., I dette tilfellet, Mi og Mj-kan ikke være i en justering sammen.

  • Φ2B: Databehandling lengden av overlapp mellom Mi og Mj. Hvis \({MO}_{i}^{j}\) er mindre enn eller lik null, så ingen overlapp eksisterer.

  • Φ2C: å Holde en kopi av Mj før eksklusive overlapper hverandre.

  • Φ2D: Hvis overlapper eksisterer, unntatt overlapp fra Mj –

    – >

  • Φ2E: Databehandling slutter posisjon Mj i T og Q.

  • Φ2F: beregning av avstand (antall symboler) mellom Mi og Mj i T og Q.

  • Φ2G: Computing antallet samsvarer ikke og avvik mellom Mi og Mj.,

  • Φ2H: beregning av straff for den samsvarer ikke og avvik mellom Mi og Mj (\(P_{i}^{j}\)). Hvis gapet, bare en gap åpne straff er trukket fra.

  • Φ2I: Databehandling justering score \(\left (S_{i}^{j}\right)\) når Mj er lagt til justering slutter på Mi. Poengsummen for alle matchende symboler i Mj (Lj×Rm) er lagt til justering score for justering slutter på Mi (Si). Så straffen for hull og samsvarer ikke mellom Mi og Mj\(\left (P_{i}^{j}\right)\) er trukket fra.,

  • Φ2J: Hvis utvide Mj til justering slutter i Mi resulterer i en poengsum som er \(\left (S_{i}^{j}\right)\) høyere enn nåværende resultat for Mj (Sj) og deretter den nye poengsum er lagret i Sj. Også W er satt til jeg til å holde styr på Mi som maksimerer score for Mj.

  • Φ2K: Gjenopprette verdien av Mj før utelukkelse slik at Mj kan brukes i andre justering utvidelser.

  • Φ3: på Jakt etter den MEM med den høyeste Sj. Dette MEM er den siste MEM i alignment (Meg)., Den høyeste scoren (Se) er tilbake som S, som er den høyeste justering score for den gitte sekvenser. Indeksen for MEM som maksimerer Sj er lagret i e for å begynne å spore tilbake fra Meg.

  • Φ4: I justeringen, umiddelbar tidligere MEM for Meg er den som maksimerer justeringen score for Meg. Indeksen for slike MEM er lagret i W. Som et resultat, iterasjonen av f←W besøk indeksen for alle MEMs i justeringen. Når B er lik -1, Mf er den første MEM i justering og iterasjonen er stoppet.,

I vår algoritme, kan vi ikke straffe samsvarer ikke og sprekker før den første MEM og etter siste MEM i justeringen. Dette resulterer i en lokal justering algoritme. Ved å vurdere disse straffer algoritmen genererer en global justering (Ekstra fil 1: Kapittel V).

ligningen for å beregne \(P_{i}^{j}\) i Φ2H av Algoritmen 2 forutsetter at det er ingen passende symbol mellom T og Q i området mellom Mi og Mj (alle symboler er regnet som samsvarer ikke eller hull)., Selv om denne forutsetningen er ikke sant for alle Mi, det er alltid sann for Mi som fører til maksimal \(S_{i}^{j}\), som overrules effekten av antagelsen blir feil for andre Mi. Som et bevis, antar det er et passende symbol i området mellom Mi og Mj. Matchende symboler ville være et MEM (Mk). Mk er allerede utvidet til justering slutter på Mi. Dermed, når utvide Mj å Mk et høyere resultat er oppnådd i forhold til å utvide Mj til Mi.

Kjeding colinear frø som diskutert i har vært mye brukt i justeringen av store sekvenser, dvs. genom-til-genom justering., Det har også blitt brukt til å identifisere kandidaten regioner for en leser gitt et sett av MEMs i BWA. Sammenkjeding av algoritmer med scoring er like dynamisk programmering algoritmen vi har foreslått (DP-MEM). Det er imidlertid forskjeller som gjør DP-MEM egnet for parvis justering av korte sekvenser. DP-MEM tar hensyn til at alle MEMs innen et visst gap størrelse er gitt i inn-og optimerer antall iterasjon i algoritmen. DP-MEM også implementerer en heuristisk tilnærming for å kompensere for effekten av kort MEMs fjernet fra input-listen som følge gap mellom MEMs.,

Optimalisering

Gitt sekvenser av lengde n, algoritmen for å trekke ut MEMs (gitt i «MEM utvinning» – delen) krever at 2(n−1) skift og 2n−1 sammenligne operasjoner på bit-vektorer (hver akt på \(\lceil \frac {n}{32} \rceil \) maskinen ord) som resulterer i en algoritme med kompleksitet O(n2) å produserer kanten bit-vektorer for gitt par av sekvenser. Kompleksiteten i funksjon som regner ut posisjonering av MEMs fra kanten bit-vektor og sorterer dem basert på EQ er ennå ikke lagt. Videre, hvis m MEMs er pakket ut, Φ2 av Algoritmen 2 (DP-MEM) har kompleksitet O(m2)., Tatt i betraktning kompleksiteten av MEM utvinning og DP-MEM -, MEM-Align er mye lavere enn tilgjengelig justering algoritmer. For å fremskynde prosessen, kan vi presentere flere optimaliseringer for MEM-Align som oppnår raskere kjøring av går på bekostning av nøyaktighet. På den annen side, vi introdusere metoder for å forbedre nøyaktigheten med en minimal ytelse tap.

Banded justering

Banded alignment er en kjent heuristisk metode for å øke hastigheten på justeringen. Denne teknikken begrenser mønster av hull i alignment (Ekstra fil 1: Kapittel VI)., Følgelig, hvis justeringen mellom to sekvenser ikke følger dette mønsteret, algoritmen vil ikke identifisere justering. I tradisjonell dynamisk programmering, justeringen er innhentet etter beregning av verdien av alle celler i tabellen. Men, med banded justering optimalisering, bare celler på diameter og nær diagonal er evaluert. gl er bredden av bandet i banded justering der cellene lenger enn gl til diameteren er ikke evaluert. Mørkere celler i Fig. 1 vis det tilfelle hvor gl=1.,

i Motsetning til tradisjonelle dynamisk programmering tilnærming, MEM-Align ikke har en tilsvarende tabell for å bruke banded justering. Men, vi fant ut at vi kan simulere den samme effekten ved å begrense antall av skiftarbeid i MEM utvinning prosessen. For eksempel, hvis vi shift spørringen sekvens opp til gl til høyre og til venstre, har vi nå kommet i samsvar med bandet av gl. Banded-justering redusere kompleksiteten av MEM utvinning fra O(n2) til O(n.(2gl+1)) hvor gl er en liten og fast verdi. Dermed kompleksiteten av MEM utvinning er O(n) når banded justering er brukt., Også, med sa banded justeringen, er det sannsynlig at færre MEMs er trukket ut som hastigheter opp den påfølgende algoritmisk trinn.

Kort MEM fjerning

Vi har observert at de fleste hentet MEMs er korte og er et resultat av tilfeldig matchende symboler. For å øke hastigheten MEM-Align, MEMs kortere enn sl er filtrert ut under MEM utvinning prosessen. Dette reduserer antall MEMs å bli hentet og behandlet (senere påskynde algoritme). Filtrering kort MEM er gjort ved å utvide Algoritme 1 med et sett av skift og bitvis operasjoner (Ekstra fil 1: Kapittel VII).,

På den annen side, det er i sjeldne tilfeller der kort MEMs er en del av justeringen, for eksempel, en matchende symboler omgitt av samsvarer ikke. Uten å ha alle MEMs i kilde-listen, DP-MEM er ikke i stand til å finne den samme justeringen som det finner når alle MEMs finnes i kilde-listen. For å kompensere for tapte kort MEMs-i inngang, vil vi endre Φ2H av DP-MEM til å vurdere muligheten av å ha korte kamper mellom MEMs (Ekstra fil 1: Kapittel VIII).

Det kan være flere vanskelige saker der i justeringen, flere kort MEMs eksisterer mellom to MEMs (se Fig. 9)., Den eneste måten å identifisere score for området mellom Mi og Mj i Φ2H er å bruke en global justering til denne regionen. Imidlertid, Φ2H er en hyppig betjening og bør være rask. Derfor bestemte vi oss for å delvis løse problemet ved å begrense mulige tilfeller (en heuristisk metode).

Fig., 9

Et eksempel som viser flere kort MEM i det lille området mellom Mi og Mj i en justering

Hvis det er hull i området mellom Mi og Mj, vi antar det er bare en kontinuerlig gap enten til venstre slutt eller til den høyre enden av området. Så, bare to tilnærmingene er mulig for området., Antall matchende symboler regnes for begge tilfeller i en sekvensiell måte, og som resulterer i maksimal kamper er tatt som antall kamper mellom Mi og Mj (Ekstra fil 1: Kapittel IX). De fortløpende sammenligning er en kostbar operasjon, og vi utarbeide en metode for å unngå sekvensielle sammenligningen når det er mulig (Ekstra fil 1: Kapittel X).

Alle andre tilfeller som ikke passer den ovennevnte forutsetning om resultater i tråd med en lavere score., Men, tatt i betraktning den lave prisen av hull og samsvarer ikke, muligheten for å ha flere hull og samsvarer ikke i et lite område er lav.

Effektiv alignment extension

I Φ2 av DP-MEM -, Mj strekker seg alle linjer som ender i {M1…Mj−1} (hvis mulig). Men for hver Mj det er en mindre undergruppe Ωj⊆{M1…Mj−1} slik at ved å utvide Mj til alle justeringer ender i en Mi∈Ωj justeringen som ender i Mj er funnet (Eq. 2). Med andre ord, det ville bli færre \(S_{i}^{j}\) skal vurderes. Definisjon av set-Ωj og bevis av Eq. 2 er gitt i annen fil 1: Kapittel XI., Definisjonen av Ωj påvirkes når MEM-kort fjerning optimalisering er aktivert (Ekstra fil 1: Kapittel XII).

$$ \maks\limits_{M_{i} \i \Omega_{j}}{S_{i}^{j}} = \maks\limits_{1 \leq jeg \leq j-1}{S_{i}^{j}} $$
(2)

Hybrid justering

for Å opprettholde nøyaktigheten av algoritmen, bestemte vi oss for å utnytte en hybrid metode som er en kombinasjon av MEM-Rett og Smith-Waterman-algoritmen. Vi definerer tre tilfeller der MEM-Align kan være unøyaktig., Hvis justeringen av et par sekvenser faller ned i en av disse sakene, bruker vi Smith-Waterman algoritme for å justere sekvenser. Disse sakene er:

  • Når sekvenser er repeterende, og antall trukket ut MEMs overstiger terskelen TM. Vi fant MEM-Align er egnet til å gi unøyaktige justering når du retter repeterende sekvenser. En passende TM verdi reduserer sjansen for rapportering unøyaktig tråd med en ubetydelig økning i den gjennomsnittlige saksbehandlingstiden.,

  • Når den beregnede justering score for justering generert av MEM-Align er lavere enn en terskel TS. Dette er tilfellet for det meste oppstår når det er et gap i den justering som ikke kan identifiseres på grunn av banded justering.

  • Når ingen MEM lenger enn sl eksisterer for å bli trukket ut (et sjeldent tilfelle)., Hvis sl er satt til en høy verdi, og likheten mellom sekvenser er lavt,

Selv om du sender sekvens par til en ekstern algoritmen resulterer i ekstra beregninger, antall sekvenser som sendes til eksterne algoritmen er fortsatt lite om aktuelle verdiene er valgt for TM og TS.

Hoppe over fjernt MEMs

Når avstanden mellom Mi og Mj er stor, er det ikke sannsynlig å ha Mi og Mj som tilstøtende MEMs i justeringen., Derfor, algoritmen hopper utvidelsen hvis avstanden mellom Mi og Mj er lengre enn en terskel TD (ytterligere å redusere antallet \(S_{i}^{j}\) skal vurderes). Denne optimaliseringen litt forbedrer ytelsen med en ubetydelig effekt på nøyaktighet.