Algoritmi per la concatenazione di match per il mapping del cDNA Flavia Rustici Miguel Gutierrez Michele Orbini Michelucci Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 2 Introduzione Genomica comparativa branca della biologia che si occupa del confronto di sequenze genomiche Confronto tra genomi di organismi in stretta relazione Confronto di regioni che conservano la sintenia (regioni in cui geni ortologhi occorrono nello stesso ordine) Stabilire le regioni codificanti e non codificanti di sequenze (struttura esoni-introni) 3 Introduzione (2) cDNA DNA complementare ottenuto dall’mRNA nella fase di trascrizione inversa Contiene solo esoni cDNA mapping Obiettivo: trovare il gene (e riconoscere la struttura esoniintroni) nel genoma da cui il cDNA è stato estratto Permette ulteriori analisi degli elementi regolatori del gene (promotori, fattori di trascrizione, binding sites, etc.) 4 Introduzione (3) Approccio software (anchor-based): 1. Calcolo dei frammenti (segmenti uguali in sequenze diverse) 2. Calcolo di una catena globale di score massimo di frammenti (ancore) 3. Allineamento delle regioni comprese tra le ancore 5 Sommario 1. Introduzione 2. Terminologia e definizioni 3. Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 6 Terminologia e definizioni Siano P e T due sequenze tali che P: pattern (cDNA) T: testo (DNA) Sia M=(m1,m2, …,mn)una lista di match Definizioni Un match (frammento) mi tra due sequenze P e T è una tripla (pi,ti,li) tale che P[ pi pi li 1] T [ti ti li 1] T (DNA) = … G G A A A T A G T T A … mi P (cDNA) = … A C C A T A G T C G C C … Esempio match mi tra 2 sequenze 7 Terminologia e definizioni(2) Due match mi e mj hanno overlap sul pattern quando pi p j pi li Due match mi e mj hanno overlap sul testo quando ti t j ti li T (DNA) = … G G A C A C A G T ….….. G T A C C C G C C mi pj mj P (cDNA) = … G G A C A C A G T A C C C G C C … overlap sulppattern pEsempio i i+li 8 Terminologia e definizioni(3) Possono esistere casi in cui due match mi e mj hanno overlap sia su pattern che su testo text_overlap(mi, mj) = max { 0, ti + li - tj}=2 T (DNA) = … G G A C A C A G T G T G T A C C C G C C mi mj P (cDNA) = … G G A C A C A G T G T A C C C G C C … pattern_overlap(mi, mj) = max { 0, pi + li - pj}=4 overlap(mi, mj) = max {pattern_overlap(mi,mj), text_overlap(mi,mj)}=4 9 Sommario 2. Introduzione Terminologia e definizioni 3. Match Chaining Problem su 2 sequenze 1. 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 10 Match Chaining Problem Definizione del problema: Input M = {m1, m2, … , mn} una lista di match ordinati in base alla loro posizione sul pattern pi i j pi p j i j ti t j Output Catena di match C { mc1 , mc2 ,, mc f } : CM i,j: mci ,mc j C, i j (pci pc j ) (tci tc j ) C deve massimizzare una funzione di score 11 Match Chaining Problem (2) score(C) misura della “bontà” di una catena di frammenti C Esempi score( C ) f , dove f è il numero di match in C f score( C ) lci , somma delle lunghezze dei match in C i 1 Algoritmo di base per la concatenazione 1. For i from 1 to n find prev(mi) O( n g( n )) compute score(mi) 2. Select the match with the largest score() and the construct the actual chain by tracing prev() entries. 12 Range Maximum Query (RMQ) Problema la cui soluzione sarà utile per il calcolo del match predecessore ( prev(m) ) Problema: Dato un array A di n interi trovare l’indice di massimo valore nel subarray A[i … j] per ogni i e j. Soluzioni: Tempo O(n) Array statico Esempio: i j A … 4 2 8 1 3 5 … RMQ( A, i, j) = 8 13 RMQ Dinamico Definizione del problema: Input insieme S di elementi, inizialmente vuoto ogni elemento I S ha associati 2 valori a(I) b(I) Output RMQ(S, p, q) elemento I di S tale che I S’ & b(I)= max { b(I’) | I’ S’} dove S’ = {I’ | I’ S & p ≤ a(I’) ≤ q } S 14 RMQ Dinamico: implementazione Necessaria una struttura dati che supporti le seguenti operazioni: Ricerca Aggiornamento (inserzione/cancellazione) AVL TREE Caratteristiche generali: Albero binario di ricerca bilanciato in altezza Ricerca e Aggiornamento O(log2 n), con n = n° nodi correnti 6 3 1 8 4 9 15 Nozioni sugli alberi Un albero è definito da un insieme di nodi V collegati tra loro da un insieme di archi E tale che: esiste un nodo speciale detto radice (root) tutti i nodi tranne la radice hanno un padre ogni nodo ha un insieme di nodi figli ad eccezione di un sottoinsieme di nodi detti foglia (leaf) che non hanno figli la relazione padre-figlio è rappresentata da un arco che collega il nodo padre al nodo figlio 16 Albero(V,E) rappresentazione grafica a b e f c g d i h m n l o p q radice nodo interno r s foglia Ogni nodo è radice di un sottoalbero 17 Nozioni sugli alberi(2) Un albero è detto binario quando ogni nodo ha al più due figli Dato un nodo v indichiamo con sx(v): il figlio sinistro di v dx(v): il figlio destro di v sinistra(v): il sottoalbero radicato nel figlio sinistro di v destra(v): il sottoalbero radicato nel figlio destro di v Esempio: a sx(a) = b dx(a) = c b sinistra(a) c destra(a) 18 Nozioni sugli alberi(3) Un albero binario è detto di ricerca quando ogni suo nodo v rispetta le seguenti proprietà: v' sin istra( v ) v' v v' ' destra( v ) v v' ' 6 3 1 8 4 9 Albero binario di ricerca ordinato su valori interi 19 Nozioni sugli alberi(4) Un albero binario è detto bilanciato in altezza quando, in ogni nodo, le altezze del sottoalbero destro e sinistro differiscono di al più 1 l’altezza di un albero binario bilanciato in altezza è O(log2 n) a hdx = 1 c b hsx = 3 d b a d e hsx - hdx = 2 f Albero binario non bilanciato f e c Albero binario bilanciato 20 Nozioni sugli alberi(5) Un albero AVL è un albero binario di ricerca bilanciato in altezza 6 3 1 8 4 9 Albero AVL L’albero AVL è la struttura dati che utilizzeremo per risolvere il problema dell’RMQ dinamico Ricerca e Aggiornamento O(log n), con n = n° nodi correnti 21 RMQ Dinamico: implementazione(2) Ogni nodo v dell’AVL-tree rappresenta un oggetto IS, ed ha associato: a(v), chiave (viene usata per l’ordinamento dei nodi nell’AVL) b(v), rappresenta lo score MAX_POINTER, puntatore al nodo con il massimo valore di b nel sottoalbero radicato in v 6 a(v) 3 6 3 1 4 8 1 3 b(v) 3 MAX_POINTER(v) 4 7 9 5 22 RMQ Dinamico: esempio Aggiornamento: insert(11,8) Fase 1: ricerca posizione b(11) = 8 6 b(MAX_POINTER(6)) b(MAX_POINTER(8)) = 57 b(MAX_POINTER(9)) 579 <non 8 aggiornare ha figli: aggiungo il nodo come figlio dx MAX_POINTER(8) MAX_POINTER(9) MAX_POINTER(6) 3 1 4 11 3 8 8 1 4 7 11 3 8 9 11 5 8 11 8 23 RMQ Dinamico: esempio Aggiornamento: insert(11,8) Fase 2: rotazione Dopo l’inserimento il nodo 8 è sbilanciato, c’è bisogno di una ROTAZIONE. 6 3 1 4 3 8 1 4 7 3 9 h( sx(8) ) - h( dx(8) ) = -2 5 11 8 24 RMQ Dinamico: esempio Aggiornamento: insert(11,8) Fase 2: rotazione 6 8 è diventato nodo foglia, si deve aggiornare MAX_POINTER(8). 3 1 4 3 9 1 4 7 8 3 5 11 8 Rotazione e aggiornamento dei MAX_POINTER aumentano la complessità dell’inserzione? 25 RMQ Dinamico: aggiornamento Complessità di un aggiornamento con rotazione v1 v2 v2 Rotazione v1 T1 T3 T2 T3 T1 T2 La rotazione richiede tempo O(1) in quanto i puntatori da aggiornare sono 3. Anche l’aggiornamento del MAX_POINTER richiede tempo O(1), in quanto riguarda solo v1 e v2: basta confrontare i loro MAX_POINTER con quelli delle radici dei sottoalberi T1, T2 e T3. La rotazione e l’aggiornamento dei MAX_POINTER non aumentano l’ordine di complessità: O(log n), dovuto a fase di ricerca della posizione in cui eseguire l’inserzione/cancellazione. 26 RMQ Dinamico: algoritmo rmq(S,p,d) Fase 1: Ricerca del nodo vleft tale che a(vleft) = min { a(I) | p ≤ a(I) , IS } Ricerca del nodo vright tale che a(vright) = max { a(I) | a(I) ≤ q , IS } Ricerca del nodo vlca, lowest common ancestor di vleft e vright Fase 2: Controllando i MAX_POINTER dei nodi lungo il cammino da vleft a vright troviamo il nodo v con il massimo b tale che p ≤ a(v) ≤ q 27 RMQ Dinamico: esempio p q Ricerca: rmq(S, 4, 9) Fase 1: ricerca vleft Pleft = [ 6 , 3 , 4 ] ricerca vright Pright= [ 6 , 9 ] ricerca vlca (lowest common ancestor) 6< > 94 → dx sx p=4 93 = < 94 ,→ trovato dx vright 4 = 4 , trovato vleft 6 3 1 4 q=9 3 9 1 4 7 8 3 5 11 8 28 RMQ Dinamico: esempio Fase2: checkSX visita ricorsivamente tutti i nodi nel cammino P1 = [ sx(vlca), …,vleft] per trovare il massimo b nel range [p, a(vlca)) checkDX è speculare per cammino P2 = [ dx(vlca), …,vright] nel range (a(vlca), q] checkSX(v,P1,m):: begin if (v == null) then return m; if (isLeaf(v)) then return maxb(v,m); if (dx(v)P1) then Vlca return checkSX(dx(v),P1,m); fuori range else //dx(v)P1 MAX_POINTER(v1) = MAX_POINTER(v3) checkSX(dx(v1)) V1 fuori range m1 = maxb(MAX_POINTER(dx(v)),v,m); V2 return checkSX(sx(v),P1,m1); end Vright V3 fuori range Vleft V4 O(log n) max_tmp V5 max_tmp = maxb { v3, MAX_POINTER(dx(v3)), max_tmp} 29 RMQ Dinamico: complessità rmq(S, p, q):: begin Pleft = search_vleft(S, p); Pright = search_vright(S, q); Plca = compute_vlca(Pleft, Pright); P1 = Pleft – Plca; O(log n) P2 = Pright – Plca; O(log n) maxleft = checkSX(vlca, P1, null); maxright = checkDX(vlca, P2, null); return maxb(maxleft, maxright, vlca); O(1) end 30 Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 2. Algoritmo base Range Maximum Query (RMQ) 3. Algoritmo MCCM 1. 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 31 MCCM (Shibuya e Kurochkin 2003) Input: M = {m1, m2, … , mn} una lista di match ordinati in base alla loro posizione sul pattern pi Output: catena di match C C={m …,,mm } : {m mc2c ,, c1c, ,m c cf} CM i,j: mci ,mc j C, i j 1 2 f (pci pc j ) (tci tc j ) C deve massimizzare una funzione di score Caratteristiche: Si applica a due sequenze genomiche (un genoma ed una sequenza di cDNA) Basato su Dynamic Range Maximum Query (RMQ) Permette: Mapping efficiente e accurato Limitazione della distanza tra frammenti Overlap di frammenti 32 MCCM: distanza massima tra match È necessario limitare la distanza nel testo tra due match consecutivi per evitare situazioni come questa MAX_LEN DNA prev(mi) mi cDNA Chiamiamo MAX_LEN la distanza massima ammissibile tra due match nella catena C. 33 MCCM con MAX_LEN f score (C ) lci i 1 1. 2. R = struttura dati RMQ, inizialmente vuota; per ogni i da 1 a nn a. per ogni match mj tale che j<i, pj+lj≤pi e mmjjR R insert(R, mm tj+l -1, -1, score(m j,j, j)); insert(R, tjj+l score(m j j)); b. c. 3. 4. prev(mi)= rmq(R, rmq(R, tit - iMAX_LEN -1, t-1, - MAX_LEN i-1);ti-1); se (prev(mi) = null) allora score(mi)=li; altrimenti score(mi)= score(prev(mi))+li; m = <trova il match con il massimo score>; C = <catena costruita tracciando i prev a partire da m> tempo O(n log(n)) Non permette overlap! 34 MCCM con overlap È utile permettere gli overlap perchè: sono molto comuni nel mapping del cDNA, occorrendo spesso alle estremità degli esoni del cDNA aumentano la quantità di sequenza ricoperta dalla catena T (DNA) = … G G A C A C A G T G T G T A C C C G C C mi mj P (cDNA) = … G G A C A C A G T G T A C C C G C C … 35 MCCM con overlap: idea L’algoritmo considera 4 tipi di match candidati per prev(mi): A. B. C. D. match che non si sovrappongono con mi match mj che hanno pattern_overlap(mj,mi)>text_overlap(mj,mi)≥0 match mj che hanno pattern_overlap(mj,mi)=0 e text_overlap(mj,mi) >0 match mj che hanno text_overlap(mj,mi)>pattern_overlap(mj,mi)>0 S2 (cDNA) mi b2 b1 d a c S1 (genoma) 36 MCCM con overlap: idea (2) 1. Per ogni match mi : a) b) 2. per ciascun insieme A, B, C e D si cerca il miglior match candidato ad essere prev(mi) prev(mi) sarà il “migliore” tra i 4 candidati Si determina il match m con il massimo score e si costruisce la catena a partire da m utilizzando l’informazione prev() f f 1 i 1 i 1 score( C ) lci overlap _ length( mci , mci1 ) 37 MCCM con overlap Ricerca del match candidato predecessore per mi di tipo A: struttura dati RMQ dinamica ( R1 ) in R1 sono inseriti i match mj tali che j<ie pattern_overlap(mj, mi) = 0 rmq(R1, ti MAX_LEN 1, ti 1) MAX_LEN : mi : mj (j < i) R1 T P Non ci sono match in R1 che abbiano overlap su pattern con mi 38 MCCM con overlap Ricerca del match candidato predecessore per mi di tipo B: struttura dati RMQ dinamica ( R2 ) in R2 sono inseriti i match mj tali che ji e pattern_overlap(mj, mi) 0 rmq(R2, ti pi MAX_LEN 1, ti pi 1) : mi : mj (j < i) R2 T P In R2 ci sono solo match che hanno overlap su pattern con mi Con l’RMQ cerchiamo il match che massimizza lo score tra quelli che hanno pattern_overlap(mj, mi) text_overlap(mj, mi) 39 MCCM con overlap Ricerca del match candidato predecessore per mi di tipo C: struttura dati coda ( Q ) implementata con albero binario di ricerca, che supporta inserimenti, cancellazioni e ricerche (non c’è bisogno di rmq) in Q sono inseriti i match mj tali che ji e pattern_overlap(mj, mi) = 0 find_best(Q,ti) : mi : mj (j < i) Q T P Non ci sono match in Q che abbiano overlap su pattern con mi. In Q andiamo a cercare l’elemento che massimizza lo score di mi tra quelli che hanno text_overlap(mj, mi) 0 41 MCCM con overlap Ricerca del match candidato predecessore per mi di tipo D: struttura dati RMQ dinamica ( R3 ) in R3 sono inseriti i match mj tali che ji e pattern_overlap(mj, mi) 0 rmq(R3, ti pi 1, ti 1) : mi : mj (j < i) R3 T P In R3 ci sono solo match che hanno overlap su pattern con mi Con l’RMQ cerchiamo il match che massimizza lo score tra quelli che hanno text_overlap(mj, mi) pattern_overlap(mj, mi) NB: non si garantisce che tj ti 42 MCCM con overlap Come si calcola il migliore tra i 4 match candidati? si calcola lo score che mi avrebbe usando come predecessore di mi ognuno dei 4 candidati si sceglie quello che massimizza score(mi) score( mi ) score( m ) li overlap( m , mi ) m prev( mi ) 44 MCCM con overlap R1, R2, R3 = strutture dati RMQ, inizialmente vuote; Q = coda, inizialmente vuota; per ogni i da 1 a n in ordine 1. 2. 3. per ogni match mj tale che j<i, pj+lj≤pi e mjR1 a) m(a) = rmq(R1,ti - MAX_LEN - 1, ti - 1); m(b) = rmq(R2,ti - pi - MAX_LEN - 1, ti - pi - 1); m(c) = find_best(Q,ti); m(d) = rmq(R3,ti - pi, ti- 1); se ((m(d) null)AND(t(m(d))>ti )) allora m(d)=nil; se ((m(a) = null)AND(m(b) = null) AND (m(c) = null) AND(m(d) = null)) allora b) c) d) e) f) g) h) i) j) k) 4. 5. insert(R1, mj, tj+lj-1, score(mj)); delete(R2, mj) insert_queue((Q, mj, tj, tj+lj-1, score(mj)-tj-lj ); delete(R3, mj) tempo O(n log(n)) prev(mi)= nil; score(mi)=li; altrimenti prev(mi)= max_score(mi, m(a), m(b), m(c), m(d)); //max_score: restituisce il match che massimizza la funzione di score per mi score(mi) = score(prev(mi) + li – overlap(prev(mi), mi); insert(R2, mi, ti-pi, score(mi)- pi - li); insert(R3, mi, ti-pi, key(score(mi)- ti - lj, ti - pi); m = <match con il massimo score>; C = <catena costruita tracciando i prev a partire da m> 45 Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 47 Mapping su sequenze multiple I genomi di specie in stretta relazione oppure i genomi di ceppi differenti di una stessa specie, condividono una grande quantità di sequenza lo stesso avviene per una sequenza di cDNA con il genoma da cui proviene Estendiamo MCCM per mappare una sequenza di cDNA su più sequenze genomiche per: identificare geni comuni a più genomi individuare le regioni di sintenia 48 Terminologia e definizioni Definizioni Un match (frammento) m su k sequenze (k-1 sequenze genomiche e una di cDNA) può essere rappresentato come un iper-rettangolo in uno spazio k-dimensionale Rk, ossia una k-tupla ([p1 … q1], …, [pk … qk]) tale che S1[p1 … q1] = …= Sk[pk … qk] S2 S3 q1 p1 S1 p2 q2 p3 q3 49 Terminologia e definizioni Possiamo definire un match m anche dai suoi angoli estremi, ossia da una coppia di k-tuple (beg(m),end(m)) tale che: beg(m) = (beg(m).x1, … , beg(m).xk) = (p1, … , pk) end(m) = (end(m).x1, … , end(m).xk) = (q1, … , qk) Relazione << tra frammenti: m’<<m (m’ precede m) sse i 1..k 1 ) beg( m' ).xi beg( m ).xi 2 ) end ( m' ).xi end ( m ).xi m’ e m si dicono colineari se m’<<m oppure m<<m’ ossia, appaiono nello stesso ordine in tutte le sequenze 50 Problema più generale Obiettivo: trovare una catena C di match non-overlaping frammenti tale che l'ammontare di sequenza ricoperta dai frammenti della catena sia massimizzato. Input M = {m1, m2, … , mn}, un insieme di n match su k sequenze Output trovare un sottoinsieme CC={f , fc2, ,… {1m m,c2fl,}, mc f } di match colineari tale 1 che: elementi di C non si sovrappongono la quantità di sequenza ricoperta dai frammenti è massimizzato. f score( C ) mci .length i 1 51 Higher-dimensional Chaining Gli algoritmi di approssimazione risolvono il problema della concatenazione dei match per più di 2 genomi Problema: complessità quadratica nel numero n di match Può essere un grande svantaggio per grandi valori di n Alternativa: Algoritmo di Zhang costruisce una catena ottima usando la divisione dello spazio basandosi su kd-alberi 52 Terminologia e definizioni KD-Albero: albero K-Dimensionale Un Albero kd (k-dimensionale) è una struttura dati per il partizionamento dello spazio, che organizza i punti su uno spazio euclídeo a k dimensioni Utilizza soltanto piani perpendicolari a uno degli assi del sistema di coordinate 53 K-Dimensional Tree Esempio di suddivisione di spazio per la creazione di kd-tree, con k = 3 La prima divisione (rosso) taglia la cella principale (bianco) in 2 sottocelle le sottocelle sono divise a loro volta (verde) in due sottocelle Queste 4 sottocelle sono ancora divise (blu) a loro volta in 2 sottocelle ciscuna Poiché non ho più divisioni da fare, le 8 celle finali sono chiamate celle foglie Le sfere gialle rappresentano i nodi dell’albero 54 Higher-dimensional Chaining Gli algoritmi di approssimazione risolvono il problema della concatenazione dei match per più di 2 genomi Problema: compessità quadratica nel numero n di match Può essere un grande svantaggio per grandi valori di n Alternativa: Algoritmo di Zhang costruisce una catena ottima usando la divisione dello spazio basandosi su kd-alberi Problema: l’analisi rigorosa della complessità in tempo di questo algoritmo è difficile perché la costruzione della catena è integrata nella struttura dell’kd-tree 55 Higher-dimensional Chaining (2) Alternativa: Algoritmo di Myers e Miller Basato sul paradigma line-sweep Ricerca ortogonale su Range Tree Tempo O(n logk n) Adesso presenteremo un algoritmo che migliora la complessità dell’algoritmo di Myers e Miller 56 Higher-dimensional Chaining (3) Ricerca ortogonale su range tree Costruzione di range tree 57 Costruzione di un range tree Dato un insieme S di n punti d-dimensionali, si costruisce il suo range tree in questo modo: Se d = 1: il range tree di S è un albero binario di ricerca di altezza minima o un array ordinato sugli alementi di S Se d > 1: il range tree di S è un albero T binario di ricerca di altezza minima con n foglie, la cui i-esima foglia più a sinistra contiene l’elemento in S con l’i-esimo valore pù piccolo della coordinata x1 associamo ad ogni nodo interno v di T un sottoinsieme canonico Cv S che contiene punti memorizzati nelle foglie del sottoalbero radicato in v Per ogni v, lv è la più piccola coordinata x1 di ogni punto di Cv hv è la più grande coordinata x1 di ogni punto di Cv C*v {(p.x 2 ..., p.x d ) R d -1 : (P.x 1 , p.x 2 ..., p.x d ) C v } Ogni nodo interno v contiene lv, hv e un range tree a d-1 dimensioni costruito su C*v Per ogni dimensione fissata d la struttura dati è costruita in 58 O(n logd-1 n) tempo e spazio Higher-dimensional Chaining (3) Ricerca ortogonale su range tree Costruzione di range tree Range Query 59 Range Query (RQ) Dato un insieme S di punti in uno spazio a d dimensioni Rd, una Range Query restituisce tutti i punti di S che si trovano nel iper-rettangolo R (p, q) = [l1... H1] × [L2... H2] × ... × [ld ... hd] Se d = 1: la query può essere risolta in tempo O(log n), dovuto alla ricerca binaria sull’albero Se d > 1: si visita il range tree a partire dalla radice 60 Visita del range tree per ogni nodo visitato v: se v è una foglia: restituiamo il punto corrispondente se si trova nell’iperrettangolo R(p, q) se v è un nodo interno: [l v . .. h v ] [l1... h1 ] {} non si fa niente 1 [l v . .. h v ] [l1... h1 ] si ricerca ricorsivamente nel range tree di dimensione d-1 memorizzato in v con l’iper-rettagolo [l2... h2] × ... × [ld ... hd] altrimenti, visitiamo ricorsivamente entrambi i figli di v Tempo O (logd n + z), con z il numero di punti nell’iperrettangolo R(p, q) 61 Range Query (RQ) Dato un insieme S di punti in uno spazio a d dimensioni Rd, una Range Query restituisce tutti i punti di S che si trovano nel iper-rettangolo R (p, q) = [l1... H1] × [L2... H2] × ... × [ld ... hd] Se d = 1: la query può essere risolta in tempo O(log n), dovuto alla ricerca binaria sull’albero Se d > 1: la query può essere risolta in tempo O (logd n + z), con z il numero di punti nell’iper-rettangolo R(p, q) visitando ricorsivamente il range tree a partire dalla radice 62 Higher-dimensional Chaining (3) Ricerca ortogonale su range tree Costruzione di range tree Range Query Possiamo migliorare l’efficienza se come ultima dimensione dell’albero se utiliziamo un albero di priorità anziché un range tree, ottenendo un algoritmo di concatenazione di complessità in tempo O(n logk-1 n) Possiamo migliorare ancora l’algoritmo mediante la tecnica di fractional cascading, che risparmia un fattore logaritmico nel rispondere alle Range Query 63 Fractional Cascading Vogliamo modificare la tecnica per RQ(0,q) per rispondere efficientemente al Range Maximum Query RMQ(0,q) Consideriamo il caso d = 2 il range tree è un x-albero binario di ricerca di y-array per ogni nodo v del x-albero indichiamo con v.L e v.R rispettivamente i suoi figli sinistro e destro l’y-array di v, indicato con Av, contiene tutti i punti in Cv ordinati in modo crescente in base alla loro coordinata y ogni elemento p Av contiene due puntatori Lptr punta all’elemento q1 più grande in Av.L tale che q1 p Rptr punta all’elemento q2 più grande in Av.R tale che q2 p se questi elementi non esistono i puntatori hanno valore NULL Esempio di questa struttura 64 Esempio Range tree con Fractional Cascading per i punti nella figura in alto a destra 65 Esempio (2) I nodi colorati sono quelli visitati dalla Range Query I nodi tratteggiati sono i nodi canonici I cerchi indicano i puntatori NULL Gli elementi cololati dell’y-array dei nodi canonici sono i punti all’interno del rettangolo di query 66 Fractional Cascading (2) La ricerca di tutti i punti nel rettangolo R(0,(h1,h2)) avviene in 2 passi: Primo passo: la ricerca binaria sul y-array del nodo radice del x-albero per localizzare il punto ph2 più a destra, tale che ph2.y [0... h2] Secondo passo: visitiamo ricorsivamente x-albero tenendo traccia dei puntatori per localizzare la foglia più a destra ph1 tale che ph1.x [0... h1] Durante la visita ricorsiva del x-albero, identifichiamo un insieme di nodi che chiameremo canonici, rispetto alla nostra query 67 Terminologia e definizioni Un insieme di nodi canonici è il più piccolo insieme P di nodi v1, ... , vl del x-albero tale che l l P : j 1 Av j j 1 C v j contiene tutti i punti p S tali che p.x [0...h1] 68 Fractional Cascading (3) Se nessun p P soddisfa p.y [0... h1] analiziamo i puntantori durante la visita del x-albero Se incontriamo un nodo canonico vj ,l’elemento ej puntato dell’ultimo puntatore visitato suddivide la lista AVj nel modo seguente: gli elementi a destra di ej non sono contenuti nell’iperrettangolo R (0, (h1, h2)) tutti gli altri elementi di AVj sono contenuti nell’iperrettangolo R (0, (h1, h2)) ej è detto split element. Il numero di nodi canonici è O(log n) Possiamo trovare tutti i nodi canonici e gli split element del loro y-array in tempo O(log n) 69 Fractional Cascading (4) Quindi, sia z il numero di punti nell’iperrettangolo R(0, q), utilizzando la tecnica di fractional cascading possiamo realizzare i Range Query in tempo: se d = 2 O(log n + z) se d > 2 O(logd-1 n + z) 70 RMQ su Range Tree Possiamo rispondere al RMQ migliorando gli yarray nel fractional cascading implementandoli con una coda di priorità costruita sull’ordinamento dei punti contenuti nell’ y-array il rank space dei punti nel y-array è l’intervallo [1 … N ], dove N è la dimensione dell’array la coda di priorità supporta: inserzioni, cancellazione, ricerca del predecessore e successore in tempo O(log log N) 71 RMQ su Range Tree (2) Per rispondere a RMQ(0,q): 1. 2. si cercano tutti i nodi canonici per ogni nodo canonico vj 1. si cerca lo split element in AVj e sia rj il suo indice all’interno del y-array 2. 3. tutti i punti in AVj a sinistra dello split elemento (esso stesso compreso) sono contenuti in R(0,q) i restanti elementi di AVj non appartengono a R(0,q) qj = predecessor(rj) è l’indice dell’elemento di massimo score in AVj l’elemento con il massimo score nel range(0,q) sarà l’elemento che ha il massimo score tra i massimi trovati su ogni nodo canonico 72 RMQ su Range Tree: complessità Rispondere a una RMQ a due dimensioni (d=2) costa in tempo O(log n log(log n)) Il numero di nodi canonici è O(log n) Ogni operazione sulla coda di priorità costa O(log log n) Se d ≥ 2 costa in tempo O(logd-1 n log(log n)) 73 Teorema Dati k>2 sequenze genomiche e n match si può trovare una catena ottima globale in : tempo O(n logk-2 n log (log n)) spazio O(n logk-2 n) Dimostrazione (cenno) Per ogni match eseguiamo una RMQ Poichè i punti sono già ordinati sulla prima dimensione l’RMQ necessita solo di d = k – 1 dimensioni Per d ≥ 2 la ricerca dell’RMQ costa spazio O(n logd-1 n) e tempo O(n logd-1 n log log n) Pertanto, trovare la catena ottima globale costa in tempo O(n logk-2 n log (log n)) e spazio O (n logk-2 n) 74 Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 1. 2. 3. 5. Limitazione della distanza tra i match Overlap Riduzione a grafo Algoritmo basato su geometria Risultati Sperimentali Riferimenti 75 Terminologia e definizioni Ad ogni match m è associato un valore di rarità r definito come il numero massimo di volte che la sottostringa Si[pi … pi+l-1] occorre nella stessa sequenza Si, per 1 i k S1 r=4 S2 S3 Dati due match m e m’ tali che i 1..k beg( m' ).xi beg( m ).xi diremo che hanno overlap sse end (m' ). xk end (m). xk & end (m' ). xk beg (m). xk S1 m’ S2 m S3 Non c’è ovelap tra m’ e m 76 Terminologia e definizioni Per ogni coppia di frammenti m’ e m su k sequenze, dove la kesima è la sequenza di cDNA, la quantità di overlap sul cDNA è calcolata come end ( m' ).xk beg( m ).xk 1 overlapk ( m' , m ) 0 altrimenti Relazione << tra frammenti: m’<<m (m’ precede m) sse i 1..k se beg( m ).xk end ( m' ).xk end ( m ).xk 1 ) beg( m' ).xi beg( m ).xi 2 ) end ( m' ).xi end ( m ).xi m’ e m si dicono colineari se m’<<m oppure m<<m’ ossia, appaiono nello stesso ordine in tutte le sequenze 77 Match chaining problem su k sequenze Input M = {m1, m2, … , mn} una lista di match su K ( 2) sequenze: S1, S2, ... , Sk-1 sequenze di DNA Sk sequenza di cDNA Output {m mc2 Catena di match colineari C C ={m …,m , cmf cf} , cioè c2 ,, c1c,1 ,m mc1 mmc2c2 << mc f , tale che la funzione m <<… c1 << cf c f 1 i 1 i 1 score C mi .length overlapk (mi , mi 1 ) sia massimizzata 78 Soluzione con grafo Il match chaining problem può essere ridotto al problema della “ricerca del cammino di costo massimo su un grafo orientato aciclico” Un grafo G(V, E) è definito da una coppia di insiemi V, insieme dei nodi (o vertici) E, insieme degli archi che collegano coppie di nodi: E (V, V) Un grafo G(V, E) si dice pesato quando ad ogni arco ei E è associato un valore w(ei) 79 G(V,E) rappresentazione grafica V = { a, b, c, d, e, f } E = { (a, b), (b, c), (b, f), (c, f), (e, f),} a w((a, b))=5 5 w((b, c))=3 w((b, f))=8 w((c, f))=0 w((e, f))=2 d 8 b 3 f 2 e 0 c non pesato Grafo pesato 80 Nozioni sui grafi Un grafo G(V, E) è orientato quando i suoi archi hanno un verso, cioè sono percorribili in un solo senso (vi, vj) E, è un arco percorribile dal nodo vi al nodo vj (vi, vj) (vj, vi) Un cammino di lunghezza l in G(V, E) è dato da una sequenza di vertici v0, v1, ... , vl tali che (vi, vi+1) E , i = 0, … , l – 1 Il costo di un cammino v0, v1, ... , vl su un grafo pesato G(V, E) è dato dalla somma dei pesi degli archi che lo compongono l 1 w(v0 , , vl ) w((vi , vi 1 )) i 0 81 Nozioni sui grafi(2) Un cammino è chiuso (ciclo) quando v0 = vl, cioè quando inizia e termina sullo stesso nodo. Un grafo aciclico è un grafo in cui non esistono cammini chiusi b f e c Grafo orientato ciclico b f e c Grafo orientato aciclico 82 G(V, E) rappresentazione grafica V = { a, b, c, d, e, f } E = { (a, b), (b, c), (b, f), (c, f), (e, f),} a w((a, b))=5 5 w((b, c))=3 w((b, f))=8 w((c, f))=0 w((e, f))=2 d 8 b 3 f 0 c 2 e esempio: w(b,c,f,e)= 3 + 0 + 2 = 5 Grafo pesato orientato aciclico 83 Riduzione a grafo Costruiamo un grafo G(V, E) tale che V rappresenta l’insieme dei match E rappresenta le relazioni di precedenza tra i match ogni match in input è rappresentato come un nodo nel grafo G sia v’ il nodo che rappresenta m’ sia v il nodo che rappresenta m esiste in E un arco orientato (v’, v) solo se m’ << m w((v’, v)) = m.length – overlapk(m’, m) si aggiungono in V i nodi 0 e t tali che il nodo 0 ha un arco uscente verso ogni altro nodo vi = mi del grafo w((0, vi)) = mi.length il nodo t ha un arco entrante proveniente da ogni altro nodo vi del grafo w((vi, t)) = 0 84 G(V, E) rappresentazione grafica V = { m1, m2, m3, m4, m5, m6,}0, 0 }t } ),}(0,m1), … , (0,m6), E = { (m2, m1), (m2, m6), (m2, m3), (m2, m5), (m3, m5), (m3, m6), (m6, m5), )} (m1,t), … , (m6,t)} 3 m2 6 m4 0 0 5 5 0 7 0 8 m1 m5 2 m6 t 0 5 3 4 3 0 4 3 0 m3 “ricerca del cammino di costo massimo da 0 a t sul grafo pesato orientato aciclico” 85 Risoluzione del problema sul grafo È un problema classico della teoria dei grafi risolto in tempo O(n2), con programmazione dinamica Inefficiente con un numero di frammenti n molto alto La riduzione costa comunque tempo O(n2) 86 Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. 5. Limitazione della distanza tra i match Overlap RMQ a dimensioni multiple Concatenzione di match overlapping 1. Riduzione a grafo 2. Algoritmo basato su geometria 3. Risultati Sperimentali Riferimenti 87 Soluzione geometrica Usando una soluzione basata sulla geometria è possibile ridurre la complessità in tempo Concetto di regione (k-dimensionale): S2 (cDNA) C(m) m A(m) = Rect (0,beg(m)) B(m) = ([beg(m).x1 … end(m).x1], [0 … beg(m).x2]) C(m) = ([0 … end(m).x1], [beg(m).x2 … end(m).x2]) AB(m) = A(m) B(m) A(m) B(m) S1 (genoma) 88 Idea ordinare i match in input m1, ..., mn e processarli in base ad un criterio di ordinamento. Per ogni match mi: cercare il miglior match predecessore tra quelli già processati, distinguendo in base al fatto che abbia o meno sovrapposizione sul pattern (cDNA) calcolare score(mi) marcare il match mi come già processato (attivo) 89 Considerazioni iniziali Consideriamo i match come coppie di punti a k dimensioni (beg(m),end(m)) con beg(m) = beg(m).x1, … , beg(m).xk end(m) = end(m).x1, … , end(m).xk D1, D2 sono strutture dati semi-dinamiche che supportano RMQ a 2k dimensioni con attivazione per manipolare i match (RMQD1,RMQD2) RMQ con attivazione: nel calcolo dell’RMQ si considerano solo i match attivi, ignorando quelli inattivi activate(D, m, p): procedura che attiva il match m (già inserito nella struttura dati D) con priorità p 90 Algoritmo Ordiniamo i punti di inizio dei match in base alla loro coordinata x1 Tutti i match sono caricati come inattivi sia in D1 che in D2 Processiamo i match come a simulare una linea che percorre i punti sulla coordinata x1 (line-sweep paradigm) 1. 2. 3. 4. Quando il match mi viene processato: gli viene associato prev(mi) viene calcolato score(mi) mi viene attivato in D1 e in D2 Prendiamo il match di massimo score e costruiamo la catena utilizzando l’informazione sul predecessore 91 Ricerca del predecessore Si sfrutta il concetto di regione per individuare i match già attivi che sono in relazione << con il match corrente mi, e tra questi distinguiamo quelli che hanno overlap sul pattern da quelli che non lo hanno: mj << mi sse beg(mj) A(mi) and end(mj) { A(mi) B(mi) C(mi) } mj ha overlapk con mi sse end(mj) C(mi) S2 (cDNA) m1 Esempio: mj è un match non (ancora) attivo non è vero che m1 << mi m2 << mi e non ha overlapk m3 << mi ed ha overlapk mi C(mi ) mj m3 m2 A(mi ) B(mi ) S1 (genoma) 92 Ricerca del predecessore(2) La ricerca dei match candidati viene fatta utilizzando RMQ a 2k dimensioni su due regioni disgiunte le prime k dimensioni devono assicurare che beg(m’) A(mi) le ultime k dimensioni devono assicurare, distinguendo i due casi, che end(m’) {AB(mi)} oppure end(m’) C(mi) m1 m2 = RMQD1(A(mi), {AB(mi)}) per la ricerca del match di massimo score che non ha sovrapposizione sul pattern (cDNA) con il match corrente mi end(m1).xk< beg(mi).xk = RMQD2(A(mi), C(mi)) per la ricerca del match di massimo score che ha sovrapposizione sul pattern (cDNA) con il match corrente mi beg(mi).xk≤ end(m2).xk ≤ end(mi).xk mi S2 (cDNA) C(mi ) m2 A(mi ) m1 B(mi ) S1 (genoma) 93 Ricerca del predecessore(3) Tra i 2 match candidati si sceglie quello che massimizza lo score di mi: score1 = m1.score + mi.length score2 = m2.score + mi.length – overlapk(m2, mi) if (score1 ≥ score2) then prev(mi)=m1 score(mi)=score1 else prev(mi)=m2 score(mi)=score2 mi S2 (cDNA) C(mi ) m2 A(mi ) m1 B(mi ) S1 (genoma) 94 Attivazione di un match Dopo aver stabilito il prev(mi) e score(mi), si deve attivare mi: in D1 il match mi viene attivato con priorità score(mi) in D2 il match mi viene attivato con priorità score(mi) – end(mi).xk Perchè usiamo priorità diverse? Quando si calcola RMQD1 si cercano frammenti senza overlap su cDNA con il match corrente mi, perciò quello che massimizza score(mi) sarà quello col più alto score score(mi) = score(m1) + mi.length Con RMQD2 si cerca il match migliore tra quelli che hanno overlap con il match corrente mi, perciò non è detto che il match di score massimo in D2 sia quello che massimizza lo score(mi) score(mi) = score(m2) + mi.length – overlapk(m2, m) = score(m2) + mi.length – ( end(m2).xk – beg(mi).xk ) = score(m2) – end(m2).xk + mi.length + beg(mi).xk 95 Ricerca della catena ottima Esempio m1 Input: M = {m1, m2, m3, m4, m5, m6, m7, m8} 5 m2 11 m3 20 m4 12 m5 9 m6 23 m7 16 m8 8 Output: C = {m1, m2, m4, m7, m6} Questa parte dell’algoritmo ha costo O(n) NB: l’ordine dei match in input non è quello con cui sono processati nella fase precedente 96 Complessità La complessità dell’algoritmo dipende dalla struttura dati utilizzata per implementare l’RMQ con attivazione in D1 e D2. Range tree: eseguire RMQ costa O(log2k-1 n log (log n)) tempo e O(log2k-1 n) spazio L’algoritmo esegue 2 RMQ per ogni match, quindi calcola per O(n) volte l’RMQ La complessità totale è O(n log2k-1 n log (log n)) tempo e O(n log2k-1 n) spazio 97 Abbassiamo la complessità È possibile ridurre la complessità dell’algoritmo riducendo le dimensioni dell’RMQ in D1 e D2 il primo range [0 … beg(mi).x1] può essere ignorato perché i match sono processati in ordine crescente secondo la coordinata x1 del punto di inizio tutti i match già attivi rispettano questo range il k-esimo range [0 … beg(mi).xk] può essere ignorato in D1: con il k-esimo range si vuole un match m1 tale che beg(m1).xk < beg(mi).xk, ma questa condizione è ridondante in quanto il 2k-esimo range controlla che end(m1).xk < beg(mi).xk in D2 si possono verificare 2 casi: beg(m2).xk < beg(mi).xk , OK beg(m2).xk ≥ beg(mi).xk, in questo caso m2 deve essere scartato. Questo non compromette l’ottimalità della catena poiché non esiste un altro match che termina in C(mi) che possa dare ad mi uno score migliore di uno che temina in AB(mi) La complessità scende a O(n log2k-3 n log (log n)) tempo e O(n log2k-3 n) spazio 98 Casi speciali L’algoritmo può essere specializzato per 1. usare in input match di tipo multi_MUM, riducendo la complessità 2. limitare la quantità di overlap ammissibile tra due match consecutivi nella catena di output, riducendo la complessità 99 Casi speciali: multi-MUM Maximal Multiple Unique Match: match con valore di rarità 1, cioè per ogni sequenza Si la sottostringa Si[pi … pi + l - 1] occorre in Si esattamente una sola volta Conseguenze: basta controllare le k dimensioni sulla fine dei match in quanto S2 (cDNA) anche il primo range [0 … end(m).x1] può essere ignorato, a causa del paradigma line-sweep S1 (genoma) end(m’) < end(m) beg(m’) < beg(m) ogni m’ con end(m’).x1 > end(m).x1 avrà beg(m’).x1 > beg(m).x1, quindi m’ non è stato ancora attivato D1 e D2 contengono punti a k - 1 dimensioni (end(m).x2,…,end(m).xk) La complessità si riduce a O(n logk-2 n log (log n)) 100 Casi speciali: limitazione sull’overlap Assunzioni: l: lunghezza minima di un match l -1: massimo overlap fra due match consecutivi nella catena Per trovare il predecessore m’ del match m, imponiamo che i 1k end (m' ).xi beg (m).xi l ma poiché l è la lunghezza minima di un match, sicuramente avremo che i 1k end (m' ).xi beg (m).xi l m S2 (cDNA) beg (m' ) A(m) C(m) l A(m) B(m) 101 l S1 (genoma) Casi speciali: limitazione sull’overlap(2) i 1k end (m' ).xi beg (m).xi l beg (m' ) A(m) S2 (cDNA) v(m) m l OK OK A (m) OK l : lunghezza minima di un match l S1 (genoma) Conseguenze: basta controllare le k dimensioni sulla fine dei match! se applichiamo il paradigma linesweep in base al punto virtuale v(m) = (beg(m).x1+l,…,beg(m).xk+l) si può eliminare anche la prima dimensione D1 e D2 contengono punti a k - 1 dimensioni (end(m).x2,…,end(m).xk ) La complessità si riduce a O(n logk-2 n log (log n)) 102 Sommario 1. 2. 3. Introduzione Terminologia e definizioni Match Chaining Problem su 2 sequenze 1. 2. 3. Algoritmo base Range Maximum Query (RMQ) Algoritmo MCCM 1. 2. 4. Match Chaining Problem su k sequenze 1. Terminologia e definizioni (2) 2. Concatenzione di match non overlapping 1. 3. RMQ a dimensioni multiple Concatenzione di match overlapping 2. Riduzione a grafo Algoritmo basato su geometria 3. Risultati Sperimentali 1. 5. Limitazione della distanza tra i match Overlap Riferimenti 103 Risultati sperimentali Per gli esperimenti sono stati utilizzati come sequenze: il cDNA del D.melanogaster (23 Mbp) cromosoma 2R del D.melanogaster (21 Mbp) cromosoma 2R del D.simulans (20 Mbp) Con il programma ramaco si sono generati frammenti di lunghezza minima 15. Sono stati confrontati 3 algoritmi su differenti valori di rarità massima dei match in input: soluzione basata su grafo (graph) soluzione geometrica (geom1) soluzione geometrica con limite di overlap (geom2) 104 Riferimenti 1. Match Chaining Algorithms for cDNA Mapping di T.Shibuya e I.Kurochkin, proceedings of 3rd Workshop on Algorithms in Bioinformatics (WABI), 462-475, 2003. 2. Chaining Algorithms and applications to comparative genomics di M.I.Abouelhoda e E.Ohlebusch, Journal of Discrete Algorithms 3(2-4), 321-341, 2005. 3. A chaining Algorithm for Mapping cDNA Sequences to Multiple Genomic Sequences di M.Abouelhoda, proceedings of 14th International Symposium on String Processing and Information Retrieval (SPIRE), 1-13, 2007. 105