Un aproccio Dividi et impera

A dispetto di quanto detto nella sezione precedente esiste un modo molto elegante che permette non solo di calcolare il costo di un allineamento con complessità nello spazio lineare, ma anche di ricostruire un allineamento con tale costo in tempo quadratico ($O(n\cdot m)$). Per poter ottenere questo risultato dobbiamo appliccare un approccio basato sul dividi et impera (in inglese divide and conquer).

Prima di procedere dobbiamo notare che sebbene noi abbiamo utilizzato la convenzione di cercare un allineamento ottimo con il massimo punteggio a partire dall'inizio delle due sequenze fino al termine, è comunque equivalente allineare le due sequenze a partire dal fondo. In pratica, se consideriamo, che le nostre due sequenze partano da 0 e siano lunghe rispettivamente $M$ e $N$ (gli indici da $0...M-1$ e $0..N-1$), allora le relazioni di ricorrenza precedenti per il calcolo di un allineamento con punteggio ottimale divengono per l'inizializzazione:

$$\begin{eqnarray*} C(0,0)=& 0 & \\ C(0,j)=& g*j &\quad \forall j=1,N-1 \\ C(i,0)=& g*i &\quad \forall i=1,M-1 \\ \end{eqnarray*}$$

e per la ricorrenza (con $\forall i=1...M \quad j=1...N$)

$$\begin{eqnarray*} C(i,j) =& max\{ C(i-1,j-1) + M(seq1(i-1),seq2(j-1)), \\ & C(i-1,j)+g, C(i,j-1)+g \} \\ \end{eqnarray*}$$

Considerando quanto abbiamo viasto fino ad ora, possiamo immaginare che sia equivalente ricavare la medesima soluzione partendo dal termine della sequenza (aggiungendo in fondo il termine di inizializzazione), e avremo l'inizializzazione

$$\begin{eqnarray*} C(M,N)=& 0 & \\ C(M,j)=& g*j &\quad \forall j=N,1 \\ C(i,N)=& g*i &\quad \forall i=M-1,1 \\ \end{eqnarray*}$$

e per la ricorrenza (con $\forall i=M-1...0 \quad j=N-1...0$)

$$\begin{eqnarray*} C(i,j) =& max\{ C(i+1,j+1) + M(seq1(i),seq2(j)), \\ & C(i+1,j)+g, C(i,j+1)+g \} \\ \end{eqnarray*}$$

Notiamo che quando consideriamo l'allineamento dall'inizio delle due sequenze, la posizione $C(i,j)$ viene calcolata considerando nel caso di match, i residui che occuperanno la posizione $i,j$, ossia (siccome le sequenze sono più corte di uno rispetto alla matrice) $M(seq1(i-1),seq2(j-1))$. Mentre quando si allinea dal termine delle due sequenze, la posizione che dovremmo considerare ler il match sarà $i+1,j+1$, che per le nostre sequenze è $M(seq1(i),seq2(j))$.

Come esempio dell'equivalenza degli algoritmo, vediamo in Figura , lo stesso costo e il medesimo all'allineamento per le due sequenze 'ALS; e 'AS'.

Traducendo il tutto in codice Python, in Figura , vediamo l'analogo di linear_cost che però calcola il costo di un allineamento tra due sequenze, occupando uno spazio lineare, a partire dal termine e arrivando all'inizio (come effettuato manualmente in Figura ).

Figura: Esempio di calcolo di allinamento partendo dall'inizio e dalla fine. Come al solito '0' per inizio, 'r' per riga, 'c' per colonna e 'd' per diagonale. Le maiucole sono state utilizzate per la traiettoria ottimale.

Figura: Calcolo del costo dell'allinamento di due sequenze con occupazione di spazio lineare, partendo dal fondo delle sequenze

Adesso che siamo convinti che sia equivalente allineare dall'inizio o dalla fine le due sequenze, possiamo passare a vedere come utilizzare un approccio dividi et impera.

Se noi consideriamo la linea che taglia a metà la matrice dello score rispetto alla prima sequenza, abbiamo che esiste un punto, in cui un allineamento con punteggio ottimo attraversa quella riga. Ricordandoci che siamo in grado di calcolare il costo di un allineamento sia partendo dall'inizio delle sequenze, sia dalla fine, possiamo valutare qual'è il punto in cui il nostro allineamento ottimo attraversa la frontiera, semplicemente calcolando la posizione per cui il punteggio è massimo tra i vettori che ci ritornano calcolando gli allineamenti parziali. Più precisamante, se abbiamo che le due sequenze sono s1=s1[0:m] e s1=s2[0:n] e la matrice dei punteggi è C possiamo verificare che la colonna (jmax) in cui il nostro allineamento attraversa la riga i-esima (C[i][jmax]), la possiamo calcolare come

$$jmax= \max_{j\in [0,N]}S^{top}[j]+S_{bottom}[j]$$ (8.5)

dove abbiamo indicato con $S^{top}$ e $S_{bottom}$, i vettori ottenuti allineando con linear_cost le sequenze s1[0:i] e s2, e con back_linear_cost (dalla fine all'inizio) le sequenze s1[i:m] e s2. Per verivicare questa affermazione, possiamo utilizzare il precedente algoritmo global_al, per calcolare la matrice di backtrace, e utilizzando l'equazione verificare che le colonne indicate di volta in volta dai valori di $jmax$ corrispondono alle colonne che si ottengono per ciascuna riga della matrice calcolata con global_al. Da cui utilizzando l'interprete Python

>>> s1,s2='ABCDX','ACLD'
>>> B # backtrace matrix calulated as usual
[['0', 'R', 'R', 'R', 'R'], 
['C', 'D', 'R', 'R', 'R'], 
['C', 'C', 'D', 'D', 'D'], 
['C', 'C', 'D', 'R', 'R'], 
['C', 'C', 'C', 'D', 'D'], 
['C', 'C', 'C', 'D', 'C']]
>>> # best path from the bottom-right corner: C D (R D) C D
...
>>> def find_jmax(sTop,sBot): 
...    jmax=0
...    smax=sTop[0]+sBot[0]
...    for j in range(1,len(sTop)):
...        if sTop[j]+sBot[j]>= smax :
...            smax=sTop[j]+sBot[j]
...            jmax=j
...    return jmax
... 
>>> len1=len(s1)
>>> for i in range(len1,0,-1): # from the last to the first row
...     (sTop,b)=linear_cost(s1[0:i],s2) # first half
...     (sBot,br)=back_linear_cost(s1[i:],s2) # second half
...     print i,s1[0:i],s1[i:],sTop,sBot
...     print "backtrace",B[i][find_jmax(sTop,sBot)]
... 
5 ABCDX  [-5.0, -3.0, -1.0, -1.0, 0.0] [-4.0, -3.0, -2.0, -1.0, 0.0]
backtrace C
4 ABCD X [-4.0, -2.0, 0.0, 0.0, 1.0] [-4.0, -3.0, -2.0, -1.0, -1.0]
backtrace D
3 ABC DX [-3.0, -1.0, 1.0, 0.0, -1.0] [-3.0, -2.0, -1.0, 0.0, -2.0]
backtrace R
2 AB CDX [-2.0, 0.0, 0.0, -1.0, -2.0] [-1.0, 0.0, -1.0, -1.0, -3.0]
backtrace C
1 A BCDX [-1.0, 1.0, 0.0, -1.0, -2.0] [-1.0, -1.0, -2.0, -2.0, -4.0]
backtrace D

L'aproccio è quindi quello di partire dal calcolo della colonna ($jmax$) in cui l'allineamento attraversa la riga centrale (mid) della matrice, dopo di che il problema è scomposto in due sottoproblemi più piccoli che riguardano la valutazione del passaggio dell'allineamento ottimale nel quarto superiore sinistro e nel quarto inferiore destro. Questa procedura viene iterata dimezzando ogni volta l'area della matrice che dobbiamo calcolare. La divisione in sottoproblemi piccoli riduce drasticamente la complessità computazionale nel tempo passando da una complessità $O(M^2N)$ a $O(MN)$, come mostrato in Figura .

La tecnica del divide and conquer, può quindi essere schematizzata come

• se il problema é ancora grande:
  • risolvi il problema come somma di k (2 per noi) sottoproblemi di dimensione N/k
• se il problema è diventato atomico:
  • risolvi l'istanza singola e ritorna il risultato al chiamante

L'agoritmo che implementa questo approccio basato sul divide and conquer è riportato in Figura . È particolarmente elegante il fatto che pur utilizzando soltanto dei vettori lunghi quanto la seconda sequenza si possa ottenere il costo ed anche un allineamento ottimo con un tempo che rimane ancora $O(MN)$.

Per meglio capire il funzionamento di quest'approccio ricorsivo in Figura riportiamo le varie fasi dell'algoritmo utilizzando le due sequenze 'ABCDX' e 'ACLD'. Ciascuna chiamata è etichettata da 0,T,B in dipendenza del fatto che siamo al livello iniziale, si considera la parte sopra della sequenza (top), o quella finale (bottom). Inotre la funzione ricorsiva, a seconda del caso che ci si trovi nel caso base o in quello generale, ritornerà immediatamente il risultato richiesto, oppure genererà a sua volta due chiamate ricorsive. Per poter meglio visualizzare le chiamate annidate le riportiamo indentate (Figura ), per cui per esempio 0.T.T.B identifica il fatto che a partire dall'inizio ('0') abbiamo diviso il problema in sottoproblemi di dimensione $1/2^3$ e che la chiamata attuale (0.T.T.B) è stata ottenuta prendendo la: metà sotto a destra ('B') della metà sopra a sinistra ('T') della metà sopra a sinistra ('T') della matrice di score.

Figura: Esempio della procedura di taglio delle soluzioni

Dalla figura soprastante, possiamo notare che calcoliamo la matrice 1 volta

poi al secondo giro la metà della matrice, al terzo firo 1/4 ecc.

Allora le operazioni che compiamo sono $NM\cdot (1+1/2+1/4+1/8+..)$

ossia se $T=MN\sum_{i=0}^N(1/2^i)=NM S_N$ per valutare $S_N$

possiamo moltiplicarlo per $(1-1/2)$ così si ha

$(1-1/2)S_N=(1-1/2)+(1/2-1/4)+..=1-1/2^{N+1}$ il che implica

$S_N=\frac{1-1/2^{N+1}}{1-1/2}=2(1-1/2^{N+1})\le 2$

osssi la complessitá nel tempo è

$C_{tempo}=O(NM\cdot S_N)=O(2NM)=O(NM)$

in pratica raddoppiamo soltanto il tempo mantenendo

comunque una complessità spaziale lineare $O(N)$

rispetto al caso di una complessità spaziale $O(NM)$

Figura: Algoritmo per l'allineamento basato sulla tecnica Dividi et impera

Figura: Esempio del funzionamento di rec_al. I livelli di indentazione corrispondono ai livelli di chiamate della funzione. T=parte sopra della sequenza, B=parte sotto.