Algoritmo di Nussinov

L'algoritmo di Nussinov, come la maggior parte degli algoritmi che cercano di predire la struttura secondaria dell'rna, trova strutture secondarie palindrome. Per cui gli pseudoknots non sono predicibili con questi approcci. L'algoritmo che vedremo ora cerca di massimizzare il numero di coppie, o più in generale un punteggio assegnato a ciascuna possibile coppia base(i),base(j). La più semplice funzione è ovviamente quella di associare un valore di $+1$ se le due basi sono complementari e $0$ altrimenti. Se definiamo rna la nostra sequenza di basi e vogliamo valutare il punteggio del possibile accoppiamento tra le basi in posizione i e j, possiamo utilizzare la funzione definita in Figura , dove abbiamo anche introdotto tol, che serve ad evitare di calcolare legami tra basi la cui separazione lungo la sequenza (|i-j|>=tol) sia troppo corta (per poter fare un legame richiediamo almeno un loop lungo una base).

Figura: Funzione per il calcolo del punteggio come basi complementari

Il problema è quello di trovare un algoritmo efficiente per trovare tra le possibili soluzioni il valore di una di quelle ottimali. L'idea di base è sempre quella di utilizzare la soluzione dei sottoproblemi per ottenere la soluzione del problema generale.

Nel caso in esame dobbiamo valutare la possibilità dei vari accoppiamenti, per cui possiamo partire considerando i possibili accoppiamenti locali (vicini in sequenza) e poi aumentarne la separazione di sequenza tenendo conto del risultato ottenuto per le separazioni di sequenza inferiori. Ciò significa partire cercando i possibili accoppiamenti tra $i$ e $i+2$, per ogni $i \in [0,L-2[$ dove $L$ è la lunghezza della sequenza di rna, poi continuare la ricerca incrementando la separazine tra le due possibili basi (cioè $i$ e $i+s$ con $s=3,4,..L-1$. Se quindi visualizziamo una matrice quadrata che rappresenta le possibili coppie che si possono creare all'interno di una sequenza, essa andrà valutata per diagonali, dalla maggiore all'angolo della matrice, come riportato in Tabella . Ricordiamo infatti che la diagonale della matrice rappresenta la base con se stessa $base(i),base(i)$, mentre la diagonale con separazione di sequenza uguale ad $s$ è $base(i),base(i+s)$.

Tabella: Esempio del calcolo per separazioni di sequenza sempre maggiori

	A	C	U	G	G	T
A	$\searrow$	$\searrow$	$\searrow$	?	?	?
C		$\searrow$	$\searrow$	?	?	?
U			$\searrow$	$\searrow$	?	?
G				$\searrow$	$\searrow$	?
G					$\searrow$	$\searrow$
T						$\searrow$

Secondo l'algoritmo che andremo a considerare, ci sono solo quattro possibli scelte che di volta in volta dovremo compiere, valutando quale tra le quattro ci fornisce il risultato massimo. Dato il fatto che noi abbiamo già valutato tutte le soluzioni per separazioni in sequenza pari a $s-1$, se vogliamo valutare quelle per $s=j-i$ con $j>i$ avremo che la soluzione ottimale per la coppia $i,j$ sará da scegliere tra

(A)

aggiungere la coppia $i$,$j$ alla miglior struttura trovata per la sottosequenza che va dalla posizione $i+1$ a $j-1$,

(B)

aggiungere una posizione spaiata i nella miglior struttura trovata per la sottosequenza che va dalla posizione $i+1$ a $j$,

(C)

aggiungere una posizione spaiata j nella miglior struttura trovata per la sottosequenza che va dalla posizione $i$ a $j-1$,

(D)

biforcare la soluzione, cercando quella posizione $k$ (con $i<k<j$) per la quale abbiamo trovato le due sottostrutture ottimali tra $i, k$ e $k+1, j$, e prendere la combinazione delle due.

Possiamo vedere la rappresentazione di quanto detto in Figura . Il corrispondente algoritmo che dovremo implementare ricercherà il massimo tra le quattro soluzioni. Data la nostra matrice dei punteggi $S$, dobbiamo risolvere la ricorrenza

$$\begin{eqnarray*} S(i,j)&=\max\{ &(A), (B), (C),(D) \} \\ &=\max\{ &S(i+1,j-1)+... ...\ & &S(i+1,j),S(i,j-1), \\ & &\max_{i<k<k}(S(i,k)+S(k+1,j)) \} \end{eqnarray*}$$

con $j-i\ge 2$. Per poter utilizzare la ricorrenza, necessitiamo anche delle condizioni al contorno, in modo da avere sempre inizializzato i valori di $S(i,j)$, da cui si ha che

$$\begin{eqnarray*} S(i,i) = S(i-1,i)=S(i,i-1)=0 \quad \forall i \end{eqnarray*}$$

È da notare che il valore dell'ultima cella della prima colonna delle matrice $S(0,L-1)$ contiene il punteggio ottenuto. La nostra implementazione in Python dell'algoritmo di Nussinov è quindi riportata in Figura . Da tale figura è evidente che la complessità computazionale dell'algorimo di Nussinov è $O(n^2)$ nello spazio, e $O(n^3)$ nel tempo come mostrano i tre diversi loops.

Figura: I quattro possibili casi considerati dall'algoritmo di Nussinov. x rappresenta la coppia di partenza di valutazione mentre o è la coppia che si vuole valutare correntemente

Figura: Algoritmo di Nussinov per il calcolo della struttura secondaria dell'RNA

Una volta ottenuta la matrice dei punteggi $S(i,j)$ (smat in Figura ) di solito si vuole anche ottenere uno dei possibili accoppiamenti con tale punteggio. Il metodo utilizzato è quello di utilizzare una coda LIFO (Last In First Out), chiamata anche stack. Le operazioni che possiamo fare su questa coda sono quelle che possiamo fare su di una pila di libri di cui esaminiamo quello in cima. Possiamo quindi:

• aggiungere un libro in cima (PUSH)
• togliere un libro dalla cima (POP)
• controllare se ci sono rimasti libri nella coda (isempty)

Lo stack, ci serve per memorizzare quale percorso si è scelto rispetto alla matrice di score (appaiamento, lungo $i$,lungo $j$, biforcazione). L'algoritmo che implementa il calcolo del backtrace sulla matrice di score è presentato in Figura , ed un esempio del suo utilizzo è

           Score is  4.0
           AACUUCUUCAA
           ((.)).((.))

Che rappresenta la struttura

             C       C
           A = U   U = A
           A = U   U = A
                 C

Figura: Algoritmo per trovare una soluzione di appaiamento secondo il punteggio ottenuto dall'algoritmo di Nussinov