Biologia Molecolare Computazionale · 3.Algoritmi euristici (BLAST) 4.Allineamento multiplo....

Post on 18-Feb-2019

217 views 3 download

Transcript of Biologia Molecolare Computazionale · 3.Algoritmi euristici (BLAST) 4.Allineamento multiplo....

Biologia Molecolare Computazionale

Paolo Provero - paolo.provero@unito.it

2008-2009

Argomenti

I Allineamento di sequenzeI Ricostruzione di alberi filogeneticiI Gene prediction

Allineamento

Allineamento di sequenze1. Scoring2. Algoritmi esatti (programmazione dinamica)3. Algoritmi euristici (BLAST)4. Allineamento multiplo

Allineamento di sequenze

ProblemaDate due sequenze (nucleotidiche, aminoacidiche o altro)determinare se sono sufficientemente simili da farci ritenereche siano derivate da un progenitore comune attraversoprocessi di mutazione.

Allineamento

Un esempio di allineamento delle sequenzeCGGGTATCCAA e CCCTAGGTCCCA è

C G G G T A - - T C C A AC C C - T A G G T C C C A

I Il simbolo “-” (indel) rappresenta un’inserzione o unadelezione avvenuta durante la storia evolutiva che haportato alle due sequenze

I Una sequenza di L indel consecutivi si dice gap dilunghezza L: questo allineamento ha un gap di lunghezza1 e uno di lunghezza 2

Tipi di allineamento

I Globale: si allineano le sequenze intereI Locale: si allineano sottosequenze delle sequenze di

partenzaI Gapped: sono permessi indelsI Ungapped: non sono permessi indelI Pairwise: di due sequenzeI Multiple: di più di 2 sequenze

Score di un allineamento

I Due possibili allineamenti di CTGTA e CGTA

1. C T G T AC - G T A

2. C T G T A -- C - G T A

I Il primo è migliore in quanto richiede un numero minore dicambiamenti (sostituzioni e inserzioni/delezioni)

I E’ necessario definire uno score per confrontarequantitativamente gli allineamenti

Scoring per allineamenti di sequenze nucleotidiche

I Per il momento adottiamo uno score semplice:I Ogni match conta +1I Ogni mismatch conta -1I Ogni indel conta -2

I Esempio 1:C T G T AC - G T A+1 -2 +1 +1 +1 +2

I Esempio 2:C T G T A -- C - G T A

-2 -1 -2 -1 -1 -2 -9

Algoritmi di allineamento

I Gli algoritmi di allineamento risolvono il seguenteproblema: dato un sistema di scoring e due sequenze,trovare l’allineamento (gli allineamenti) con lo score piùalto tra tutti quelli possibili

I Gli algoritmi esatti individuano certamente il miglioreallineamento, ma richiedono molte risorse di calcolo

1. Needleman-Wunsch: global gapped alignments2. Smith-Waterman: local gapped alignments

I Gli algoritmi euristici non garantiscono che l’allineamentotrovato sia il migliore possibile, ma sono molto più veloci

1. BLAST (gapped or ungapped local alignments)

Perché abbiamo bisogno di algoritmi?

I In linea di principio sarebbe possibile trovare l’allineamentoottimale di due sequenze semplicemente:

1. Scrivere tutti gli allineamenti possibili2. Calcolare il relativo score3. Scegliere il migliore

I Tuttavia questo è possibile soltanto per sequenze moltocorte, in quanto il numero di allineamenti possibili crescemolto rapidamente al crescere della lunghezza dellasequenza

Perchè abbiamo bisogno di algoritmi

I Supponiamo di avere due sequenze di n = 1000 residuiI Il numero di allineamenti possibili è

N =

(2nn

)∼ 22n√πn

=22000√

1000π∼ 2 · 10600

Programmazione dinamica

I Gli algoritmi di programmazione dinamica suddividono ilproblema iniziale (trovare l’allineamento ottimale di duesequenze) in sottoproblemi più semplici e costruisconoprogressivamente la soluzione del problema globaleusando le soluzioni dei sottoproblemi

I L’algoritmo di Needleman-Wunsch calcola lo scoredell’allineamento ottimale di due sequenze in base agliscore degli allineamenti ottimali delle loro sottosequenze,progressivamente più lunghe

Algoritmo di Needleman-Wunsch

I Sequenze:x = X1X2 . . .Xm

y = Y1Y2 . . .Yn

I L’algoritmo costruisce la matrice B tale che Bij è lo scoredell’allineamento ottimale tra le prime i lettere di x e leprime j lettere di y

I I valori di Bij vengono calcolati progressivamente fino aottenere Bmn, che è lo score dell’allineamento ottimale trax e y

Algoritmo di Needleman-Wunsch

I L’algoritmo di Needleman-Wunsch si basa sul fattoseguente:

se si conoscono Bi−1,j , Bi,j−1 e Bi−1,j−1 è possibilecalcolare Bij

I Una volta ottenuto lo score dell’allineamento ottimale iltraceback consente di ottenere l’allineamento o gliallineamenti corrispondenti

Esercizio

Trovare un allineamento globale ottimale tra le sequenzeCATTGAATCT

Soluzione: score -2Allineamento (non l’unico):

C - A T - TG A A T C T-1 -2 +1 +1 -2 +1 -2

Algoritmo di Smith-Waterman

I L’algoritmo di Smith-Waterman trova l’allineamento localeottimale tra due sequenze, ovvero l’allineamento di scorepiù alto tra tutti quelli possibili tra una qualsiasisottosequenza di x e una qualsiasi sottosequenza di y

I L’algoritmo procede in modo simile a Needleman-Wunscheccetto che:

1. Ogni volta che Bij è negativo, si scrive 02. Il traceback non parte da Bmn ma dallo score più alto

trovato sulla tabella

Esercizio

Trovare un allineamento locale ottimale tra le sequenzeCATTGAATCT

Soluzione: score 2Allineamento:

A TA T+1 +1 +2

Gap penalties

I Il tipo di gap penalty che abbiamo usato si chiama lineare:la penalty (score negativo) associata a un gap èproporzionale alla lunghezza del gap:

δ(L) = −d · L

I Un altro tipo di gap penalty usato spesso si chiama affine:

δ(L) = −d − (L− 1) · e

con e < d . In questo modo aprire un nuovo gap è piùpenalizzato che allungarne uno esistente

Gap penalty affine: esempio

I Esempio: d = 2, e = 1A C C C TA - - C T+1 -3 +1 +1 0

A C C C TA - C - T+1 -2 +1 -2 +1 -1

Allineamento di sequenze aminoacidiche

I Nel caso di allineamento di sequenze di aminoacidi, unsistema di scoring come quello usato per le sequenzenucleotidiche non è adeguato, in quanto alcune sostituzionisono più probabili di altre

I Quindi i mismatch avranno score diversi a seconda degliaminoacidi coinvolti (matrici di sostituzione)

I Gli algoritmi di allineamento funzionano nello stesso modo,eccetto che lo score per un match/mismatch varia aseconda degli aminoacidi coinvolti

I Le matrici di sostituzione più usate sono BLOSUM e PAM

Esercizio

Usare le matrici di sostituzione BLOSUM62 e una linear gappenalty uguale a -5 per trovare un allineamento ottimaleglobale tra le sequenzeHAGHEAE

Soluzione: score 5Allineamento:

H - A GH E A E8 -5 4 -2 5

BLAST

Algoritmo euristico per allineareI Query sequence vs target sequenceI Query sequence vs sequence database

Versione originale: ungapped alignments.BLAST2: gapped alignments.

BLAST steps (simplified)

1. Generare una lista di tutte le parole di k (p.es. k = 3 perproteine) lettere nella queryPQGEFG→ (PQG, QGE, GEF, EFG)

2. Per ognuna di queste generare una lista di tutte le possibiliparole di 3 lettere che allineate danno un punteggiomaggiore di una soglia T (p.es. T = 13 con BLOSUM62)PQG→ PQG (18) PEG (15) . . .

3. Cercare nella sequenza target queste parole ed allinearlealla parola della query

4. Estendere in modo massimale gli allineamenti nelle duedirezioni (ottenendo gli High Scoring Pairs)

5. Se possibile unire gli HSP per ottenere allineamenti piùlunghi.

BLAST

Grafi

I Un grafo è un insieme di nodi e di archiI Ogni arco unisce due nodiI Un cammino tra il nodo n e il nodo m è una successione di

archi consecutivi che porta da n a mI Un circuito è un cammino tra un nodo e se stessoI Un grafo è connesso se per ogni coppia di nodi (n,m)

esiste un cammino tra n e m

Alberi

I Un grafo connesso e privo di circuiti si dice alberoI Dati due nodi di un albero, esiste uno e un solo cammino

che li unisceI Un albero si dice binario senza radice se ogni nodo è

connesso a uno o tre archiI Un albero si dice binario con radice se ogni nodo è

connesso a uno o tre archi, tranne un unico nodo, laradice, che è connesso a due archi

I I nodi connessi a un solo arco si dicono foglie

Ricostruzione di alberi filogenetici

I Problema: dato un insieme di specie attualmente esistentiricostruire l’albero filogenetico che riproduce il percorsoevolutivo attraverso il quale esse sono derivate da unantenato comune

I L’albero filogenetico è rappresentato da un albero binario,con o senza radice

I La lunghezza degli archi è proporzionale al tempotrascorso tra gli eventi di speciazione

I Le specie esistenti sono rappresentate dalle foglieI La radice, se specificata, rappresenta l’antenato comuneI Se l’albero filogenetico è senza radice, non è specificata la

direzione del tempo evolutivo tra nodi interni

Ricostruzione di alberi filogenetici

I Distance methodsI Parsimony methodsI Statistical methods

Distanze

I Una distanza è una regola per associare a ogni coppia(x , y) di punti di un insieme S un numero d(x , y) ≥ 0 taleche

1. d(x , y) = 0 se e solo se x = y2. d(x , y) = d(y , x)3. d(x , y) + d(y , z) ≥ d(x , z)

I Per esempio la distanza geometrica tra punti di un pianosoddisfa questi assiomi

Distanza evolutiva

I La distanza evolutiva tra due specie esistenti si puòdefinire come il doppio del tempo trascorso dal loroantenato comune più recente

I E’ facile convincersi che questa definizione soddisfa gliassiomi che definiscono una distanza

I Se si considera un albero filogenetico con radice, ladistanza evolutiva tra le specie x e y è proporzionale allalunghezza del cammino tra x e y

Distanza ultrametrica

I La distanza evolutiva ha la seguente proprietà:dati tre punti x , y , z e le tre distanze tra di essi, due diqueste sono uguali tra loro e maggiori della terza

I Una distanza che soddisfi queste proprietà si diceultrametrica

Metodi basati sulle distanze

I Vedremo due metodi basati sulle distanze:1. UPGMA: ricostruisce alberi filogenetici con radice2. Neighbor joining: ricostruisce alberi filogenetici senza

radiceI L’algoritmo UPGMA (Unweighted Pair Group Method using

Arithmetic averages) permette di risolvere il problemaseguente:dato un insieme di specie esistenti e le loro distanzeevolutive, ricostruire l’albero filogenetico con radice che nerappresenta l’evoluzione da un antenato comune

UPGMA

1. Individuare le due specie più vicine tra le N esistenti2. Definire un nuovo nodo come loro antenato comune,

ponendolo all’altezza uguale a metà della distanza tra idue nodi

3. Eliminare le due specie dalla lista e sostituirle conl’antenato comune

4. Ripetere con la nuova lista di N − 1 specie5. Ripetere fino a che la lista contiene una sola specie (la

radice = antenato comune)6. Si può dimostrare che il metodo riproduce l’unico albero

filogenetico con radice che riproduce le distanze evolutive

Distanza tra gruppi

I Per calcolare la distanza tra nodi interni si usa la formula

d(A,B) =1

nAnB

∑x∈A;y∈B

d(x , y)

Esempio

I Le distanze evolutive tra 4 specie sono date nella tabellaseguente (Mya)

Human Chimp Bonobo GorillaHuman 0 12 12 14Chimp 12 0 4 14

Bonobo 12 4 0 14Gorilla 14 14 14 0

I Verificare la proprietà ultrametricaI Ricostruire l’albero filogenetico con radice usando UPGMA

UPGMA in pratica

I In pratica le distanze evolutive in genere non sono noteI Si usano allora distanze surrogate, chesi suppone siano

approssimativamente proporzionali alle distanze evolutiveI Le distanze surrogate in genere non sono ultrametricheI La conseguenza è che le distanze calcolate sull’albero

così prodotto non sono uguali alle distanze in input

Esempio

I Considerare la matrice di distanza tra 3 specie:A B C

A 0 6 8B 6 0 4C 8 4 0

I Mostrare che si tratta di una matrice di distanzaI Mostrare che questa distanza non è ultrametricaI Costruire l’albero filogenetico con radice usando UPGMA e

calcolare le distanze sull’alberoI Mostrare che queste distanze sono ultrametriche e che

non coincidono con quelle in input

Neighbors

I Dato un albero filogenetico senza radice, due foglie sidicono neighbors se il cammino che le unisce passa perun solo nodo interno

I Notare che non è necessariamente vero che il neighbor diuna specie è la specie più vicina nel senso della distanzacalcolata sull’albero

Tree-derived distances

I Data una matrice di distanze il metodo del neighbor joiningcostruisce l’albero unrooted che riproduce le distanze date,se questo esiste

I Una distanza per la quale esista tale albero si dice“tree-derived” o “additive”

I Una distanza ultrametrica è tree-derived, ma l’opposto nonè necessariamente vero

δ(x , y)

I Supponiamo di avere una matrice di distanze d(x , y) traspecie

I Definiamo la quantità

δ(x , y) = (N − 4)d(x , y)−∑

n 6=x ,y

(d(x ,n) + d(y ,n))

I Se x e y sono tali che la δ(x , y) è minima, allora x e ysono neighbors

Neighbor joining

I Calcolare la matrice delle δI Unire la coppia di specie con il delta più piccolo attraverso

un nodo interno r1 con le distanze

d(x , r1) =d(x ,n)− d(y ,n) + d(x , y)

2

d(y , r1) =d(y ,n)− d(x ,n) + d(x , y)

2dove n è un’altra foglia qualsiasi

I Sostituire x e y con r1. La distanza tra r1 e le altre foglie èdata da

d(r1,n) =d(x ,n) + d(y ,n)− d(x , y)

2

I Ripetere fino a che non si ricostruisce l’intero albero

Esempio

I Calcolare le distanze trale foglie dell’alberorappresentato in figura

I Applicare l’algoritmo"neighbor joining" alledistanze calcolate emostrare che l’alberoottenuto riproduce quellodi partenza

x y

z

w

1

1

1

4

5

Rooting

I E’ possibile localizzare la radice di un albero ottenuto conneighbor-joining aggiungendo un outgroup (specie che sisa essere più distante dalle altre di quanto queste siano traloro)

I La radice dell’albero originario coincide con il punto da cuisi diarma l’outgroup.

Esempio

I Costruire l’albero unrooted di 3 specie dalle seguentidistanze:

A B CA 0 3 4B 3 0 5C 4 5 0

I Localizzare la radice usando l’outgroup D, con distanze

d(D,A) = 12 d(D,B) = 11 d(D,C) = 14

Maximum Parsimony

I Trovare l’albero che spiega la discendenza delle specieesistenti da un antenato comune con il minimo numero dimutazioni.

I Per semplicità ci limiteremo a considerare sostituzioniI Esempio: le specie esistenti sono rappresentate dalle

seguenti sequenze:AAGAAAGGAAGA

Costo di un albero

AAA

AAA AGA

AAG AAA GGA AGA

1

0101

Cost = 3

AAA

AAA AGA

AAG AAA GGA AGA

01

0101

Cost = 3

AAA

AAA AAA

AAG AGA AAA GGA

0

0

2011

Cost = 4

I Elencare i possibili alberiI Calcolare il costo per ogni alberoI Scegliere l’albero di costo minore

Alberi possibili

Gli alberi si distinguono perI TopologiaI Assegnazione delle specie esistenti alle foglieI Assegnazione di sequenze ai nodi intermedi

Topologie per 4 specie

Assegnazioni per 4 specie

A B C D A C B D A D B C

B A C D B C A D B D A C

C A B D C B A D C D B A

D A C B D B A C D C A B

A B C D

A C B D

A D C D

Numero di alberi

I Il numero di assegnazioni (sommato su tutte le topologie) è

(2N − 3)!

2N−2(N − 2)!

I Esempio: per N = 4

(2N − 3)!

2N−2(N − 2)!=

5!

4 · 2!= 15

I Questo numero diventa rapidamente molto grande: perN = 10:

17!

288!= 34.459.425

I Con metodi euristici si può usare parsimony anche percentinaia di specie

Allineamenti multipli: SP scoresIl modo più comune di assegnare uno score a un allineamentomultiplo è un Sum of Pairs (SP) score:

I Lo score totale è la somma degli score delle singolecolonne:

S =∑

i

Si

I Lo score di una colonna è la somma degli score di tutte lecoppie di simboli (usando una matrice di sostituzioneappropriata)

Si =∑k<l

s(mki ,m

li )

I Gaps (linear gap penalty):

s(−,a) = −d

s(−,−) = 0

Esempio

match = 1, mismatch = -1, d = -2

A A C GA A C TA - C TA G C T+6 -7 +6 0 =5

SP scores: problema

A A C T G GA A C T G GA A C A G GA A C T G GA A C T G G

+10 +10 +10 +2 +10 +10 =52

I Per N sequenze una colonna di T vale N(N−1)2

I Un singolo mismatch costa N − 1I Costo relativo di un singolo mismatch: 2/N: decresce

all’aumentare di N

Over-counting of evolutionary events

Algoritmi per l’allineamento multiplo

I Esistono algoritmi di programmazione dinamica perl’allineamento multiplo, ma sono utilizzabili solo per Npiccolo

I Tra gli algoritmi euristici i più comuni sono gli algoritmi diallineamento progressivo:

I Allineare due sequenzeI Allineare una terza sequenza all’allineamento ottenutoI Continuare ad aggiungere una sequenza alla volta

all’allineamento

Algoritmo di Feng-Doolitle

1. Trovare tutti gli N(N − 1)/2 pairwise alignments delle Nsequenze

2. Costruire un albero “filogenetico” usando come distanzauna funzione dello score dei pairwise alignments

3. Allineare le due sequenze “evolutivamente” più vicine4. Allineare gli altri nodi dell’albero nell’ordine i cui sono stati

aggiunti all’albero

Per allineare due insiemi di sequenze, calcolare tutti i pairwisealignments e usare il migliore

Esempio

Usando lo scoring:match = 1; mismatch = -1; gap=-2allineare AAG, ATG, AA :

A A GA T GA A -+3 -1 -3 =-1

Applicazione: Positional Weight Matrices

I Siti di legame di TF:I Brevi sequenze di DNA (∼ 6− 15 bps)I Non completamente conservate

I Una descrizione possibile è la sequenza consensus,ricavabile da un allineamento multiplo di siti di legame noti:

A C C G G TA C T G G T- C T G G AA C C G G TA C Y G G W

IUPAC codes:I Y = pyrimidine (C or T)I W = T or A

Positional Weight Matrices

Una descrizione più accurata è data dalla PWM:

A C C G G TA C T G G T- C T G G AA C C G G TA C Y G G W

n A C G T1 3 0 0 02 0 4 0 03 0 2 0 24 0 0 4 05 0 0 4 06 1 0 0 3

Ricerca di TFBS con PWM

Le PWM possono essere usate per identificare potenzialibinding sites di TFs:

I Considerare le sequenze regolatrici dei geni di interesse(promotore, primi introni, ...)

I Fare “scorrere” la PWM sulla sequenza calcolando unoscore a ogni posizione

I Lo score esprime la somiglianza tra la sequenza e lamatrice

Log-likelihood ratio

Un sistema di scoring usato in pratica è il seguente: data unaPWM M e una sequnenza S:

S = log2P(S|M)

P(S|B)

doveI P(S|M) è la probabilità di generare S a partire da MI P(S|B) è la probabilità di generare S a partire dalle

(appropriate) frequenze nucleotidiche di background

Esempio

n A C G T1 3 0 0 02 0 4 0 03 0 2 0 24 0 0 4 05 0 0 4 06 1 0 0 3

S = ACTGGA

P(S|M) = 1 · 1 · 12· 1 · 1 · 1

4=

18

P(S|B) =

(14

)6

= 2−12 = 4096

P(S|M)

P(S|B)= 2−3 + 212 = 29

score = 9

(N.B. per evitare problemi con log(0)si aggiungono “pseudocounts” per lebasi che non compaiono mai inciascuna posizione)

Ricerca di TFBS

I Ogni volta che la sequenza supera un certo score minimodeciso a priori, la si identifica come un candidato TFBS

I Problema: in genere si trovano troppi candidati:I PWM di lunghezza 6I Accettiamo solo il massimo score possibile (1 sola

sequenza)I Su sequenze casuali, otteniamo un candidato ogni

46 = 4,000 basiI Sul genoma umano, ci aspettiamo

3 · 109

4 · 103 ∼ 7.5 · 105

candidati semplicemente per caso

Soluzioni

1. Spesso i TFBS funzionali sono ripetuti nella regioneregolatrice⇒ Selezionare soltanto le regioni regolatriciin cui il numero di candidati è significativamente piùalto di quanto ci si aspetta per caso

2. Spesso i TFBS funzionali sono evolutivamente conservati⇒ Selezionare soltanto i candidati TFBS che siritrovano in altre specie