Algoritmi per la concatenazione di match per il mapping del cDNA

1. Algoritmi per la concatenazione di match per il mapping del cDNA Flavia Rustici Miguel Gutierrez Michele Orbini Michelucci

2. Sommario • Introduzione • Terminologia e definizioni • Match Chaining Problem su 2 sequenze • Algoritmo base • Range Maximum Query (RMQ) • Algoritmo MCCM • Limitazione della distanza tra i match • Overlap • Match Chaining Problem su k sequenze • Terminologia e definizioni (2) • Concatenzione di match non overlapping • RMQ a dimensioni multiple • Concatenzione di match overlapping • Riduzione a grafo • Algoritmo basato su geometria • Risultati Sperimentali • Riferimenti

3. 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)

4. 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 esoni-introni) nel genoma da cui il cDNA è stato estratto • Permette ulteriori analisi degli elementi regolatori del gene (promotori, fattori di trascrizione, binding sites, etc.)

5. Introduzione (3) • Approccio software (anchor-based): • Calcolo dei frammenti (segmenti uguali in sequenze diverse) • Calcolo di una catena globale di score massimo di frammenti (ancore) • Allineamento delle regioni comprese tra le ancore

6. Sommario • Introduzione • Terminologia e definizioni • Match Chaining Problem su 2 sequenze • Algoritmo base • Range Maximum Query (RMQ) • Algoritmo MCCM • Limitazione della distanza tra i match • Overlap • Match Chaining Problem su k sequenze • Terminologia e definizioni (2) • Concatenzione di match non overlapping • RMQ a dimensioni multiple • Concatenzione di match overlapping • Riduzione a grafo • Algoritmo basato su geometria • Risultati Sperimentali • Riferimenti

7. T (DNA) = … G G A A A T A G T T A … P (cDNA) = … A C C A T A G T C G C C … Esempio match mi tra 2 sequenze 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) mitra due sequenze P e Tè una tripla (pi,ti,li) tale che mi

8. T (DNA) = … G G A C A C A G T ….….. G T A C C C G C C mj pj mi P (cDNA) = … G G A C A C A G T A C C C G C C … pi pi+li Terminologia e definizioni(2) • Due match mi e mj hanno overlap sul pattern quando • Due match mie mj hanno overlap sul testo quando Esempio overlap sul pattern

9. T (DNA) = … G G A C A C A G T G T G T A C C C G C C mj mi P (cDNA) = … G G A C A C A G T G T A C C C G C C … Terminologia e definizioni(3) • Possono esistere casi in cui due match mie mj hanno overlap sia su pattern che su testo text_overlap(mi, mj) = max { 0, ti + li - tj}=2 pattern_overlap(mi, mj) = max { 0, pi + li - pj}=4 overlap(mi, mj) = max {pattern_overlap(mi,mj), text_overlap(mi,mj)}=4

10. Sommario • Introduzione • Terminologia e definizioni • Match Chaining Problem su 2 sequenze • Algoritmo base • Range Maximum Query (RMQ) • Algoritmo MCCM • Limitazione della distanza tra i match • Overlap • Match Chaining Problem su k sequenze • Terminologia e definizioni (2) • Concatenzione di match non overlapping • RMQ a dimensioni multiple • Concatenzione di match overlapping • Riduzione a grafo • Algoritmo basato su geometria • Risultati Sperimentali • Riferimenti

11. Match Chaining Problem Definizione del problema: • Input • M = {m1, m2, … , mn} una lista di match ordinati in base alla loro posizione sul pattern pi • Output • Catena di match : • C deve massimizzare una funzione di score

12. Match Chaining Problem (2) • score(C) • misura della “bontà” di una catena di frammenti C • Esempi • , dove fè il numero di match in C • , somma delle lunghezze dei match in C • Algoritmo di base per la concatenazione • For i from 1 to n • find prev(mi) • compute score(mi) • Select the match with the largest score() and the construct the actual chain by tracing prev() entries.

13. A … 4 2 8 1 3 5 … j i RMQ( A, i, j) = 8 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: • Array statico • Esempio: Tempo O(n)

14. 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

15. Caratteristiche generali: • Albero binario di ricerca bilanciato in altezza • Ricerca e Aggiornamento O(log2 n), con n = n° nodi correnti 6 3 8 1 4 9 RMQ Dinamico: implementazione • Necessaria una struttura dati che supporti le seguenti operazioni: • Ricerca • Aggiornamento (inserzione/cancellazione) AVL TREE

16. 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

17. a b c d i l f g h e radice p q m n o nodo interno r s foglia Albero(V,E) rappresentazione grafica • Ogni nodo è radice di un sottoalbero

18. Esempio: sx(a) = b dx(a) = c a sinistra(a) destra(a) c b 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

19. 6 3 8 1 4 9 Nozioni sugli alberi(3) • Un albero binario è detto di ricerca quando ogni suo nodo v rispetta le seguenti proprietà: Albero binario di ricerca ordinato su valori interi

20. a b hdx = 1 c b a d hsx = 3 hsx - hdx = 2 d e f e c f 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) Albero binario bilanciato Albero binario non bilanciato

21. 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 6 3 8 Albero AVL 1 4 9 Nozioni sugli alberi(5) • Un albero AVL è un albero binariodi ricerca bilanciato in altezza

22. 6 6 1 4 9 a(v) 3 3 5 7 4 b(v) 3 8 1 3 MAX_POINTER(v) RMQ Dinamico: implementazione(2) • Ogni nodo v dell’AVL-tree rappresenta un oggetto IS, 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

23. 11 11 11 8 8 8 6 4 9 11 1 3 5 7 4 8 3 8 1 3 RMQ Dinamico: esempio • Aggiornamento: insert(11,8) • Fase 1: ricerca posizione b(11) = 8 b(MAX_POINTER(6)) = 7 b(MAX_POINTER(8)) = 5 b(MAX_POINTER(9)) = 5 5 < 8 aggiornare MAX_POINTER(9) 5 < 8 aggiornare MAX_POINTER(8) 7 < 8 aggiornare MAX_POINTER(6) 9 non ha figli: aggiungo il nodo come figlio dx

24. 8 6 1 4 9 11 3 3 4 7 5 8 3 1 RMQ Dinamico: esempio • Aggiornamento: insert(11,8) • Fase 2: rotazione Dopo l’inserimento il nodo 8 è sbilanciato, c’è bisogno di una ROTAZIONE. 8 3 h( sx(8) ) - h( dx(8) ) = -2 9 5 11 8

25. 8 6 4 1 11 9 3 3 5 7 4 8 3 1 RMQ Dinamico: esempio • Aggiornamento: insert(11,8) • Fase 2: rotazione 8 è diventato nodo foglia, si deve aggiornare MAX_POINTER(8). Rotazione e aggiornamento dei MAX_POINTER aumentano la complessità dell’inserzione?

26. La rotazione e l’aggiornamento dei MAX_POINTERnon aumentano l’ordine di complessità: • O(log n), dovuto a fase di ricerca della posizione in cui eseguire l’inserzione/cancellazione. RMQ Dinamico: aggiornamento • Complessità di un aggiornamento con rotazione • La rotazione richiede tempo O(1) in quanto i puntatori da aggiornare sono 3. • Anche l’aggiornamento del MAX_POINTERrichiede 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. v1 v2 v2 T1 Rotazione v1 T3 T2 T3 T1 T2

27. RMQ Dinamico: algoritmo • rmq(S,p,d) • Fase 1: • Ricerca del nodovlefttale chea(vleft) = min { a(I) | p ≤ a(I) , IS } • Ricerca del nodovright tale chea(vright) = max { a(I) | a(I) ≤ q , IS } • Ricerca del nodovlca, lowest common ancestor di vleft e vright • Fase 2: • Controllando i MAX_POINTER dei nodi lungo il cammino da vleftavrighttroviamo il nodo v con il massimo b tale che p≤ a(v) ≤ q

28. q p q = 9 p = 4 8 1 11 3 4 8 3 1 RMQ Dinamico: esempio • Ricerca: rmq(S, 4, 9) • Fase 1: • ricerca vleft • ricerca vright • ricerca vlca (lowest common ancestor) , 3 , 4 ] Pleft = [ 6 Pright= [ 6 , 9 ] 6 6 < 9 → dx 6 > 4 → sx 3 3 < 4 → dx 9 = 9 , trovato vright 9 5 4 4 = 4 , trovato vleft 7

29. Vlca checkSX(v,P1,m):: begin if (v == null) then return m; if (isLeaf(v)) then return maxb(v,m); if (dx(v)P1) then return checkSX(dx(v),P1,m); else //dx(v)P1 m1 = maxb(MAX_POINTER(dx(v)),v,m); return checkSX(sx(v),P1,m1); end V4 Vleft V3 V5 Vright V1 V2 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] fuori range MAX_POINTER(v1) = MAX_POINTER(v3)  checkSX(dx(v1)) fuori range fuori range O(log n) max_tmp max_tmp = maxb { v3, MAX_POINTER(dx(v3)), max_tmp}

30. rmq(S, p, q):: begin Pleft = search_vleft(S, p); Pright = search_vright(S, q); Plca = compute_vlca(Pleft, Pright); P1 = Pleft– Plca; P2 = Pright – Plca; maxleft = checkSX(vlca, P1, null); maxright = checkDX(vlca, P2, null); return maxb(maxleft, maxright, vlca); end O(1) RMQ Dinamico: complessità O(log n) O(log n)

31. Sommario • Introduzione • Terminologia e definizioni • Match Chaining Problem su 2 sequenze • Algoritmo base • Range Maximum Query (RMQ) • Algoritmo MCCM • Limitazione della distanza tra i match • Overlap • Match Chaining Problem su k sequenze • Terminologia e definizioni (2) • Concatenzione di match non overlapping • RMQ a dimensioni multiple • Concatenzione di match overlapping • Riduzione a grafo • Algoritmo basato su geometria • Risultati Sperimentali • Riferimenti

32. 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 = {mc1, mc2, … , mcf} : • 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

33. MAX_LEN MCCM: distanza massima tra match • È necessario limitare la distanza nel testo tra due match consecutivi per evitare situazioni come questa • Chiamiamo MAX_LEN la distanza massima ammissibile tra due match nella catena C. DNA prev(mi) mi cDNA

34. MCCM con MAX_LEN • R = struttura dati RMQ, inizialmente vuota; • per ogni i da 1 a n • per ogni match mjtale che j<i, pj+lj≤pie mjR insert(R, mj, tj+lj-1, score(mj)); • prev(mi)= rmq(R, ti- MAX_LEN -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> n mjR insert(R, mj, tj+lj-1, score(mj)); rmq(R, ti- MAX_LEN -1, ti-1); tempo O(n log(n)) Non permette overlap!

35. T (DNA) = … G G A C A C A G T G T G T A C C C G C C mj mi P (cDNA) = … G G A C A C A G T G T A C C C G C C … 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

36. S2 (cDNA) mi d a b1 b2 c S1 (genoma) MCCM con overlap: idea • L’algoritmo considera 4 tipi di match candidati per prev(mi): • 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

37. MCCM con overlap: idea (2) • Per ogni match mi : • 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()

38. MAX_LEN T P : mi : mj (j < i)  R1 MCCM con overlap • Ricerca del match candidato predecessore per mi di tipo A: • struttura dati RMQ dinamica ( R1) • in R1 sono inseriti i match mjtali che • j < i e • pattern_overlap(mj, mi) = 0 • rmq(R1, ti  MAX_LEN  1, ti  1) Non ci sono match in R1 che abbiano overlap su pattern con mi

39. T P : mi : mj (j < i)  R2 MCCM con overlap • Ricerca del match candidato predecessore per mi di tipo B: • struttura dati RMQ dinamica ( R2) • in R2 sono inseriti i match mjtali che • j  i e • pattern_overlap(mj, mi)  0 • rmq(R2, ti  pi  MAX_LEN  1, ti  pi  1) 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)

40. T P : mi : mj (j < i)  Q 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 mjtali che • j  i e • pattern_overlap(mj, mi) = 0 • find_best(Q,ti) Non ci sono match in Qche 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. T P : mi : mj (j < i)  R3 MCCM con overlap • Ricerca del match candidato predecessore per mi di tipo D: • struttura dati RMQ dinamica ( R3) • in R3 sono inseriti i match mjtali che • j  i e • pattern_overlap(mj, mi)  0 • rmq(R3, ti  pi  1, ti  1) 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)

43. R1, R2, R3= strutture dati RMQ, inizialmente vuote; Q = coda, inizialmente vuota; per ogni i da 1 a n in ordine per ogni match mjtale che j<i, pj+lj≤pie mjR1 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) 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 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> MCCM con overlap tempo O(n log(n))

44. Sommario • Introduzione • Terminologia e definizioni • Match Chaining Problem su 2 sequenze • Algoritmo base • Range Maximum Query (RMQ) • Algoritmo MCCM • Limitazione della distanza tra i match • Overlap • Match Chaining Problem su k sequenze • Terminologia e definizioni (2) • Concatenzione di match non overlapping • RMQ a dimensioni multiple • Concatenzione di match overlapping • Riduzione a grafo • Algoritmo basato su geometria • Risultati Sperimentali • Riferimenti

45. 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 Mapping su sequenze multiple

46. q1 p1 S1 q2 p2 S2 q3 p3 S3 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]

47. Terminologia e definizioni • Possiamo definire un matchm 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 • m’ e m si dicono colineari se m’<<m oppure m<<m’ • ossia, appaiono nello stesso ordine in tutte le sequenze

48. 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 C={f1, f2, … , fl} di match colineari tale che: • elementi di C non si sovrappongono • la quantità di sequenza ricoperta dai frammenti è massimizzato.

49. 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

50. 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