Teoria e ricerca di patterns in biosequenze

 5
Il nostro intento è estrarre solo patterns diagnostici da famiglie costituite sperimentalmente, cioè 
formate da stringhe che in laboratorio abbiano presentato caratteristiche biologiche comuni; questa 
scelta è motivata dal fatto che se i patterns diagnostici       ( fin tanto che sono tali) si presentano in 
stringhe recentemente decodificate possono essere usati per inserire le stesse nelle famiglie che 
abbiano quel pattern associato. A volte però si sceglie di inserire una nuova stringa in una famiglia  
perché presenta le caratteristiche salienti della famiglia stessa, e può accadere che con tale 
inserimento i patterns diagnostici per quella famiglia non si adattino alla nuova stringa, cioè 
perdano la loro “diagnosticità”: in tal caso bisogna calcolare un nuovo insieme di patterns 
diagnostici per la famiglia allargata. Questo procedimento è senza dubbio più corretto del 
precedente in quanto inserire una stringa in una famiglia solo perché in essa occorre un pattern 
proprio di quella famiglia è un processo puramente computazionale che non sempre si accorda con 
gli esperimenti di laboratorio. Queste due situazioni - cioè l’usare vecchie informazioni per 
catalogare nuove sequenze e usare le nuove sequenze per modificare le vecchie informazioni - sono 
i due aspetti di un metodo chiamato “machine learning”.       
Il motivo per cui vengono usati i patterns piuttosto che altri strumenti per la caratterizzazione di 
gruppi di biosequenze è di natura strettamente biologica. In natura le specie si evolvono per 
diversificazione mediante la manifestazione di mutazioni del DNA e l’eredità di tali caratteri da 
parte delle generazioni successive.  
Questi caratteri si trasmettono dalle regioni codificanti del DNA alle proteine, quindi si può 
rappresentare la storia evolutiva di una famiglia di proteine con un albero le cui foglie sono le 
proteine della famiglia attualmente esistenti attraverso vari organismi  e costituiscono quindi una 
famiglia di proteine. La radice dell’albero è la proteina primordiale (detta proteina madre) che ha 
dato origine a tutte le proteine dell’albero, mentre i nodi intermedi contengono tutte le proteine 
modificate. I rami uscenti da un nodo rappresentano eventi di diversificazione in cui le mutazioni 
hanno dato origine a nuove proteine; ciò che risulta più difficile è trovare nelle sequenze 
discendenti dei residui della stringa madre, e qui entra in gioco un altro fattore biologico: le proteine 
possono essere decomposte in due regioni, quelle vitali per le operazioni della proteina e quelle che 
hanno un ruolo di supporto. Le seconde sono molto soggette alle mutazioni  perenni, mentre le zone 
vitali operano una sorta di “pulizia” delle modifiche nelle generazioni discendenti, in modo tale da 
rimanere sempre simili alla sequenza madre. 
Grazie a queste regioni in cui si preserva la struttura madre il modo più adatto per modellizzare quel 
che si conserva è dato dai patterns, che esprimono bene sia agli amminoacidi conservati che le zone 
con mutazioni (dette regioni jolly). 
 6
I problemi che intendiamo risolvere fanno parte del così detto “pattern matching”, estensione dello 
“string matching” detto anche “corrispondenza tra stringhe”, che è la ricerca di stringhe in un testo. 
È utile allora vedere quali siano i principali algoritmi di risoluzione di quest’ultimo problema per 
comprendere come essi abbiano influenzato gli algoritmi di pattern matching. 
Nel primo capitolo viene data la formalizzazione del problema dello string matching e la soluzione 
detta “banale”, cioè quella che usa tutti i possibili confronti tra il testo e la stringa da ricercare. 
Naturalmente tale processo dà un costo esponenziale, quindi l’intento è migliorare l’efficienza di 
questo algoritmo con opportuni accorgimenti o con algoritmi di carattere completamente diverso. A 
quest’ultima categoria appartiene l’algoritmo proposto da Rabin e Karp: esso converte il testo e la 
stringa di input in stringhe numeriche e tramite il confronto dei loro valori determina tutte le 
possibili occorrenze della stringa. Sebbene tale soluzione nel caso peggiore dia lo stesso costo 
dell’algoritmo banale, è preferibile a quest’ultimo nel caso medio in quanto ne dimezza i costi. Una 
seconda e più efficiente soluzione è quella proposta da Knuth, Morris e Pratt: essa ha come base 
l’algoritmo banale ma a differenza di esso opera dei confronti all’interno della stringa per saltare il 
confronto con alcune posizioni del testo. Nonostante questo algoritmo sia sensibilmente migliore 
del precedente non costituisce la migliore soluzione conosciuta che invece è data dall’algoritmo di 
Boyer e Moore. Anche questo opera dei salti di posizione nel testo, ma il metodo è molto più 
semplice perché più vicino all’algoritmo banale: esso utilizza due euristiche, quella del carattere 
discordante e quella del buon suffisso. Entrambe vengono applicate quando si trova nel testo una 
serie di caratteri uguali alla stringa seguita da un carattere discordante come per l’algoritmo di 
Knuth, Morris e Pratt, ma si sviluppano autonomamente e viene scelta l’euristica che dà il salto 
minore.  
Il secondo capitolo si occupa del problema del pattern matching: esso consiste nella ricerca di 
parole ripetute all’interno di un testo, ma a differenza della corrispondenza tra stringhe non si 
cercano ripetizioni esatte di tali parole, bensì sottostringhe che abbiano un certo numero di errori 
rispetto ad esse. Data la natura genetica del problema il testo di input è un insieme di sequenze 
biologicamente affini per le quali viene specificato anche il tipo e la quantità di errori ammessi. 
L’algoritmo che risolve tale problema viene raggiunto mediante una serie di algoritmi intermedi in 
ognuno dei quali si risolve un problema simile a quello del pattern matching; l’algoritmo “base” è 
quello di Karp, Miller e Rosenberg che cerca tutte le parole ripetute almeno q volte (per q fissato) 
nelle N stringhe di input, dove la lunghezza di tali parole è un numero fissato oppure la più grande 
possibile. Questo algoritmo esegue una serie di passi in ognuno dei quali raddoppia la lunghezza 
delle parole trovate e si ferma quando raggiunge la lunghezza massima oppure quando non è più 
possibile estendere le parole trovate nei passi precedenti.  
 7
Dal momento che in genetica in certe situazioni alcuni caratteri dell’alfabeto (che sia quello degli 
amminoacidi o quello dei nucleotidi) hanno lo stesso significato, si crea una sorta di relazione di 
similarità tra tali caratteri: infatti se in una stringa compare ad esempio il carattere A che è simile al 
carattere B, allora da essa possiamo estrarre due parole che si differenziano solamente per tali 
caratteri. Questo ci porta a modificare l’algoritmo di Karp, Miller e Rosenberg in modo tale da 
ammettere caratteri simili nell’estrazione di parole dalle stringhe di input. In questo modo la 
quantità di parole estratte da un insieme di biosequenze aumenta proporzionalmente alla quantità di 
caratteri simili. Il passo successivo è la modifica dell’algoritmo di K-M-R in modo tale da 
ammettere errori; diciamo che tra due stringhe c’è un errore se si differenziano per un carattere in 
più, in meno o semplicemente diverso e non simile. Una miglioria all’algoritmo di K-M-R è data 
dall’algoritmo con alberi suffissi che utilizza una forma non compatta di alberi, i suffix tries. Con 
tali alberi K-M-R subisce l’ultima modifica necessaria per avere l’algoritmo più efficiente possibile 
e sfrutta il fatto che nelle stringhe di input  si presentano molte copie esatte di ogni pattern.  
Il capitolo 3 affronta lo stesso argomento ma modifica leggermente la struttura dei patterns che 
diventano qui espressioni contenenti al loro interno le posizioni e il tipo di errori ammissibili. 
Vengono usati cioè i patterns in notazione PROSITE, di cui vedremo una definizione più precisa. I 
patterns vengono prima estratti dai cammini di una struttura denominata “ grafo pattern” ( di cui 
vedremo la costruzione) e poi modificati dalle operazioni di “generalizzazione”  che hanno lo scopo 
di allargare l’insieme dei patterns trovato con tutte le soluzioni possibili: questo perché il grafo 
pattern è una struttura talmente compatta che fa estrarre solamente i patterns “minimali”, cioè quelli 
che costituiscono un sistema di generatori per l’insieme delle soluzioni.  
 8
 
Capitolo 1 
 
 
 
 
 
 
Algoritmi di string matching 
 
 
Il problema dello string matching è il seguente: dato in input un testo T e una stringa P trovare tutte 
le occorrenze di P all’interno del testo. Se indichiamo con m la lunghezza del testo e con n quella 
della stringa, le occorrenze sono allora un insieme di numeri compresi tra 1 e m che indicano le 
posizioni in cui compaiono le copie di P. Formalmente possiamo formulare il problema in questo 
modo:  
 
Def: Dato un testo T[1,…,n] di n caratteri e una stringa P[1,…,m] di lunghezza m, il problema dello 
string matching consiste nel trovare tutte le posizioni s di T tali che: 
 
                          T[s+1,…s+m]=P[1,…m]  dove 0≤s≤n-m. 
 
Ogni posizione s soddisfacente tali condizioni è detta “spostamento valido”. 
 
Il problema dello string matching consiste quindi nel ricercare in un testo tutti gli spostamenti validi 
rispetto una data stringa. 
Un primo algoritmo per risolvere tale problema è quello banale, che confronta tutte le posizioni di T 
con quelle di P, quindi ha un costo O((n-m+1)m) nel caso peggiore (che si verifica quando occorre 
confrontare tutti i caratteri di P con m-1 caratteri di T per ognuna delle n posizioni di T). 
In genere però non occorre eseguire tutti i confronti possibili, infatti si può confrontare ogni 
carattere da sinistra verso destra fino a quando non si trova una discordanza, e dato che raramente 
tali discordanze si presentano all’ultimo carattere, la maggior parte delle volte non si eseguono tutti 
i confronti, ma mediamente la metà di essi. 
 9
Esempio: Consideriamo su Σ={A,B,C} il testo T=AABCABBCAC e la stringa P=CAC. 
In tal caso si confronta solo il primo carattere di P con le prime tre posizioni di T in quanto sono 
tutte e tre diverse da esso, si confrontano poi la quarta e la quinta posizione di T rispettivamente con 
i primi due caratteri di P e poi di nuovo la quinta con il primo carattere di P. 
Da questo passaggio si nota come un miglioramento dell’algoritmo possa consistere nel non dover 
eseguire dei confronti successivi sulla base dei precedenti, ad esempio in quest’ultimo passaggio si 
potrebbe evitare di riconfrontare la quinta posizione con la stringa dato che si è trovato al passo 
precedente una discordanza.  
Gli algoritmi che vedremo in questo capitolo sfruttano proprio questo “recupero” delle informazioni 
avute da confronti precedenti per saltare confronti successivi, e ognuno di essi avrà una modalità 
diversa per eseguire tali salti.         
Tali algoritmi vengono presentati in ordine crescente di efficienza e il costo di ognuno sarà 
confrontato con quello dell’algoritmo ingenuo e dell’algoritmo precedente.  
Il primo algoritmo che presenta una miglioria rispetto a quello banale è quello di Rabin-Karp: esso 
decrementa il costo medio (che si verifica quando m=⎣n/2⎦ ) portandolo a O(n
2
) e non porta alcuna 
miglioria al costo della situazione peggiore. Tale algoritmo opera una conversione delle stringhe in 
numeri in base d per un dato d e poiché si potrebbero avere stringhe con valori molto elevati tutti i 
calcoli vengono eseguiti modulo un appropriato numero q scelto in base alla più lunga parola 
rappresentabile dalla macchina su cui l’algoritmo viene eseguito. Un secondo e più efficiente 
algoritmo è quello di Knuth-Morris-Pratt: esso migliora sia il caso medio che quello pessimo 
sfruttando le informazioni sui confronti precedenti unitamente ai caratteri della parola da ricercare 
nel testo. La soluzione più efficiente che si conosca per lo string-matching è data dall’algoritmo di 
Boyer-Moore che è molto simile all’algoritmo banale ma lo migliora sensibilmente con uno studio 
dei caratteri interni alla parola da ricercare. 
 10
 
Algoritmo di Rabin-Karp 
 
In questo algoritmo il testo (che è la stringa di input) T viene precedentemente codificato in cifre in 
base d dove d=|Σ| e Σ è l’alfabeto dei caratteri del testo. 
Assumiamo inizialmente d=10 e vediamo poi come estendere l’algoritmo a d generico. 
Indichiamo con p il valore decimale della stringa P e con t
s
 il valore decimale di T[s+1,…,s+m]. E’ 
chiaro che  s è uno spostamento valido se e solo se t
s
=p. 
Mostriamo che tutte le sottostringhe di T di lunghezza m (cioè tutti i  t
s 
) possono essere calcolate in 
tempo O(m); è ovvio che t
0 
si calcola in tempo O(m), ma per i successivi t
i  
si procede 
ricorsivamente usando la formula: 
 
t
s+1
= 10(t
s
-10
m-1
T[s+1])+T[s+m+1] 
 
infatti la cifra t
i+1
 si calcola da t
i 
rimuovendo la prima cifra con 10
m-1
T[i+1] (detta cifra 
maggiormente significativa)  e aggiungendo la nuova cifra a destra con T[i+m+1] (detta cifra meno 
significativa). 
Ogni t
i 
 viene determinato in tempo costante, quindi il calcolo di tutti i restanti n-m  t
i  
costa O(n-m). 
Il costo totale dei t
i 
 è dunque O(n). 
Per p il discorso è analogo, lo si calcola da P con la regola di Horner: 
 
p =P[m]+10(P[m-1]+10(P[m-2]+…+10(P[2]+10P[1])…)) 
 
e quindi viene calcolato in un tempo O(m). 
Complessivamente la costruzione di p e degli t
i
 è O(n+m). Una volta ottenute queste stringhe si 
eseguono i confronti di p con ciascuno dei t
i 
∀ i=0,…n-m, confronti che non vengono fatti cifra a 
cifra, ma che vengono fatti considerando il valore decimale di p e t
i
.  
Il costo totale della procedure è allora O(n+m), ma quando  p e t
i 
 sono molto grandi sorgono dei 
problemi causati da cifre che potrebbero non essere conformi alla parola della macchina. In tali 
situazioni si assumono i numeri modulo q con q in modo tale che 10q sia rappresentabile dentro una 
parola della macchina.  
Possiamo allora generalizzare al caso non decimale assumendo un alfabeto {0,1,…d-1} e 
scegliendo q in modo tale che dq sia rappresentabile in una parola della macchina.  
 11
Tutti i numeri saranno allora in base d e modulo q, e tale calcolo è eseguibile in tempo O(n+m).   
I valori p e t
i
 si calcolano con: 
t
s+1
=(d(t
s
-T[s+1]h)+T[s+m+1])     mod q 
 
dove h≡d
m-1
    mod q è il valore della cifra “1” nella posizione più significativa di una porzione di 
testo di m cifre. 
L’unico problema è che dobbiamo cercare tutte le posizioni s tali che t
s
=p, ma in questo modo  
troviamo solo quando t
s
≡p   mod q, che non è la stessa cosa. L’unico modo in cui sarebbe utile tale 
congruenza sarebbe di usarla per trovare le posizioni s tali che  t
s
Τp   mod q così da escludere delle 
posizioni di T. Occorrerà quindi fare un confronto diretto tra le restanti posizioni di T e la stringa P.  
Il programma ha una parte introduttiva in cui h viene inizializzato al valore della cifra più 
significativa di una porzione di m cifre e p e t
0
  a 0. Viene eseguito un ciclo for con cui vengono 
calcolati p e t
0
 in base d e modulo q e un altro ciclo con cui vengono analizzati tutti i t
s
 ( a partite da 
s=0) per cercare le posizioni valide e dentro il quale viene generato il successivo  t
s
.  
Il caso peggiore si presenta quando tutte le posizioni sono potenzialmente valide, cioè t
s
≡p   mod q 
per ogni s e quindi occorre confrontare sempre m caratteri del testo con gli m della stringa, 
situazione che dà un costo uguale a quello dell’algoritmo banale, cioè un costo quadratico, ma R-K 
si comporta molto meglio nel caso medio, quando cioè sono poche le posizioni potenzialmente 
valide e quindi il costo è O(n+m) sommato alla quantità di volte in cui occorre confrontare le 
posizioni potenzialmente valide. 
 
Esempio: Consideriamo il testo e la stringa in base 10 cioè Σ={0,1,2,3,4,5,6,6,7,8,9} e sia P=1234 
la stringa da ricercare nel testo T=478292123478. 
In tale situazione la lunghezza della stringa è m=4 e assumiamo q=11. Dal momento che la stringa P 
dà p≡2   mod 11, cerchiamo nel testo tutte le posizioni s tali che t
s
≡2   mod11 dove le t
s
 sono 
assunte di lunghezza m=4.  
Per s=1 si ha t
1
=4782≡ 8   mod 11 che non va bene; il successivo valore t
2
 viene calcolato in tempo 
costante con la formula sopra che elimina da t
1
 la prima cifra e inserisce alla sua estremità destra la 
cifra del testo nella posizione s+m+1. Otteniamo così t
2
=7829≡8 mod 11, che scartiamo; lo stesso 
ragionamento porta a t
3
=8292≡2349 mod 11, t
4
=2921≡6 mod11, t
5
=9212≡5 mod11, t
6
=2123≡0 
mod11,t
7
=1234≡2 mod11, quindi per s=7 abbiamo un candidato, ma ciò come sappiamo non basta, 
occorre verificare posizione per posizione se è effettivamente valida.  
 12
In questo caso lo è, ma ciò non esclude che ce ne possano essere altre, allora continuiamo a 
confrontare i restanti t
s
 con p: t
8
= 2347≡4 mod 11, t
9
=3478≡2 mod 11 che dà nella nona posizione 
un possibile candidato, ma confrontando tale posizione con la stringa si osserva che non è valida, e 
quindi l’unica occorrenza della stringa è nella settima posizione. 
 
Algoritmo di Knuth-Morris-Pratt 
   
La principale causa dell’inefficienza dell’algoritmo ingenuo è la perdita di informazioni avute dai 
precedenti confronti tra la stringa P e il testo T .  
L’algoritmo di Knuth-Morris-Pratt sfrutta proprio questo principio per migliorare il tempo di 
esecuzione dell’algoritmo banale.  
Supponiamo che in un confronto tra P e T i primi q caratteri siano coincidenti ma che i restanti     
m-q+1 non lo siano; può accadere che alcune posizioni di T (vedremo che esse sono k) comprese 
nei q caratteri combacianti siano valide, e quindi sia possibile saltare alcune prime posizioni dei q 
caratteri combacianti senza doverle confrontare. 
E’ possibile determinare questa situazione semplicemente eseguendo alcuni confronti con i caratteri 
interni alla stringa; in dettaglio, se  P[1,…,q]=T[s+1,…,s+q] per un certo q (cioè dalla posizione s 
di T  i primi q caratteri di P combaciano per q caratteri in T) bisogna stabilire se esiste una 
posizione s
1
 di T tale che: 
 
P[1,…,k]=T[s
1
+1,... s
1
+k]  con k tale che : s
1
+k=s+q  
cioè se esiste una posizione s
1 
dopo s che sia un potenziale candidato basandosi solo sulla 
conoscenza di  T[s+1,…,s+q]. 
Tutto questo ovviamente è utile se s
1
 è determinabile senza effettuare tutti i confronti tra s e s
1
. Ne 
segue anche che la situazione migliore si presenta quando s
1
=s+q, quando cioè nessuna delle q-1 
posizioni di T è un  possibile spostamento valido. 
 L’unico problema nella richiesta precedente è che il problema ha due incognite, s
1
 e k, ma vedremo 
come k possa essere ottenuto indipendentemente da s
1
.   
 13
E’ proprio k ad essere calcolabile mediante un confronto interno dei caratteri di P: infatti dal 
momento che T[s
1
+1,... s
1
+k] deve essere parte della porzione conosciuta di T, deve essere un 
suffisso di P[1,…,q]:=P
q
, quindi quel che bisogna determinare è semplicemente il più grande k 
minore di q  tale che  P
k
 sia suffisso di P
q
. In tal modo il successivo spostamento potenzialmente 
valido è s
1
=s+(q-k).  
Il valore k viene calcolato in base al numero di caratteri coincidenti tra testo e stringa, quindi viene 
naturale indicarlo con una funzione così definita: 
 
π(q)=max{k<q : P
k
 sia suffisso di P
q
} 
 
detta “funzione prefisso”. 
Informalmente
 
π(q) è la lunghezza del più lungo prefisso di P che sia suffisso di P
q
. 
In realtà una volta calcolato π(q) si esegue un processo iterativo con cui si determinano tutti i 
suffissi del suffisso precedente, cioè π(q), π(π(q)), π(π(π(q))),…, ottenendo così l’insieme π*(q). 
Dato che la funzione prefisso è calcolabile in un tempo O(m) e che l’algoritmo K-M-P ha un costo 
O(m+n), il costo totale dell’algoritmo banale con la procedura K-M-P è O(n+m). 
 
Esempio: consideriamo la stringa P =abcababcababa nel testo: 
T 
abaabcababcabbabaabacbaabc
1234567891011121314151617181920212223242526
 
 
allora la posizione s=3 dà come valida la stringa P
10
=abcababcab, dunque q=10, e si calcola 
π(10)=2 in quanto il più lungo suffisso di P
10
 che sia anche prefisso è ab che compare nella nona 
posizione di P
10
. Si può allora calcolare l’insieme dei valori di π∗(10) con la tabella della funzione 
π: 
 
π(i) 
i 012 34567 8910111213
P[i]a bca babca ba b a
-100-10200-102021
  
 
con tale tabella si trova che π(10)=2, e iterando tale risultato si trovano tutti gli altri k =π
i
(π(10)) che 
danno tutte le posizioni potenzialmente valide: π(2)=0, cioè π*(10)={2,0} ottenendo così tutte le 
prossime posizioni da analizzare applicando la formula s’=s+q-π(10), cioè l’insieme {3+10-2,3+10-
0}={11,13}.