3 STATE MODELS
L’idea di combinare le previsioni in una struttura genica completa è che i vincoli ‘grammaticali’ possono escludere alcuni assemblaggi di esoni sbagliati. La struttura grammaticale del problema è stata sottolineata da David Searls (Searls, 1992; Dong e Searls, 1994), che ha anche proposto di utilizzare i metodi delle grammatiche formali dell’informatica e della linguistica. La programmazione dinamica può spesso essere descritta convenientemente da una sorta di automa a stati finiti (Searls e Murphy, 1995; Durbin et al., 1998). Un modello potrebbe avere uno stato per l’inizio della traduzione (S), uno per i siti donatori (D), uno per i siti accettatori (A) e uno per la fine della traduzione (T). Ogni volta che viene fatta una transizione da uno stato all’altro, viene aggiunto un punteggio (o una penalità). Per la transizione dallo stato donatore allo stato accettore, il punteggio dell’introne viene aggiunto al punteggio totale, e così via. Nella Figura 11.2 è mostrato il diagramma di stato per il semplice algoritmo di programmazione dinamica di cui sopra. Per ogni variabile nell’algoritmo c’è uno stato corrispondente con lo stesso nome, e anche uno stato iniziale e uno finale sono necessari.
Figura 11.2. Un automa a stati finiti corrispondente all’algoritmo DP semplice.
Il vantaggio di una tale formulazione è che la programmazione dinamica per trovare il punteggio massimo (o la pena minima) è di tipo più generale, e quindi aggiungere nuovi stati o nuove transizioni è facile. Per esempio, disegnare il diagramma di stato per un algoritmo di programmazione dinamica più generale che permette qualsiasi numero di geni e anche di geni parziali è semplice (Figura 11.3), mentre è complicato da scrivere. Allo stesso modo il diagramma di stato per l’algoritmo frame-aware abbozzato sopra è mostrato nella Figura 11.4.
Figura 11.3. Il modello della Figura 11.2 con l’aggiunta di transizioni che permettono la predizione di qualsiasi numero di geni e geni parziali dove la sequenza inizia o finisce nel mezzo di un esone o introne.
Figura 11.4. Un modello che assicura la coerenza della cornice in tutto un gene. Come nelle due figure precedenti, le linee tratteggiate corrispondono alle regioni intergeniche, quelle tratteggiate agli introni e quelle piene alle regioni codificanti (esoni).
Se i punteggi utilizzati sono probabilità log o probabilità log, allora un automa a stati finiti è essenzialmente un modello di Markov nascosto (HMM), e questi sono stati introdotti di recente nel gene finding da diversi gruppi. L’unica differenza fondamentale dagli schemi di programmazione dinamica discussi in precedenza è che questi modelli sono completamente probabilistici, il che ha alcuni vantaggi. Uno dei vantaggi è che il problema della ponderazione è più semplice.
VEIL (Henderson et al., 1997) è un’applicazione di una HMM al problema del gene finding. In questo modello tutti i sensori sono HMM. Il modulo degli esoni è essenzialmente una catena di Markov disomogenea del primo ordine, descritta sopra. Questo è l’ordine naturale per l’implementazione in una HMM, perché allora ciascuna delle probabilità condizionali della catena di Markov disomogenea corrisponde alla probabilità di una transizione da uno stato al successivo nella HMM. Non è possibile evitare i codoni di stop nel frame di lettura quando si usa un modello del primo ordine, ma in VEIL alcuni stati in più sono aggiunti in modo intelligente, il che rende la probabilità di un codone di stop pari a zero. I sensori per i siti di splice sono fatti in modo simile. I singoli moduli sono poi combinati essenzialmente come nella Figura 11.2 (cioè la coerenza del frame non è applicata). Il modello combinato è un grande HMM, e tutte le transizioni hanno probabilità associate. Queste probabilità possono essere stimate da un insieme di dati di allenamento con un metodo di massima verosimiglianza. Per combinare i modelli questo si riduce essenzialmente a contare le occorrenze dei diversi tipi di transizioni nel set di dati. Pertanto, la ponderazione implicita dei singoli sensori non è davvero un problema.
Anche se il modo in cui viene trovata la struttura genetica ottimale è simile nello spirito alla programmazione dinamica di cui sopra, sembra molto diverso in pratica. Questo perché la programmazione dinamica è fatta a livello dei singoli stati in tutti i sottomodelli; ci sono più di 200 stati in VEIL. Poiché il modello è completamente probabilistico, si può calcolare la probabilità di qualsiasi sequenza di stati per una data sequenza di DNA. Questa sequenza di stati (chiamata percorso) determina l’assegnazione di esoni e introni. Se il percorso passa attraverso il modello degli esoni, quella parte della sequenza è etichettata come esone; se passa attraverso il modello degli introni è etichettata come introne, e così via. L’algoritmo di programmazione dinamica, chiamato algoritmo Viterbi, trova il percorso più probabile attraverso il modello per una data sequenza, e da questo si ricava la struttura genetica prevista (vedi Rabiner (1989) per un’introduzione generale alle HMM).
Questo modello probabilistico ha il vantaggio di risolvere il problema della ponderazione dei singoli sensori. La stima di massima verosimiglianza dei parametri può essere dimostrata essere ottimale se ci sono sufficienti dati di allenamento, e se la natura statistica dei geni può essere descritta da un tale modello. Una parte debole di VEIL è il modello degli esoni del primo ordine, che probabilmente non è in grado di catturare le statistiche delle regioni codificanti, e la maggior parte degli altri metodi usa modelli del quarto o quinto ordine.
Un gene finder basato su HMM chiamato HMMgene è attualmente in sviluppo. Il metodo di base è lo stesso di VEIL, ma include diverse estensioni alla metodologia HMM standard, che sono descritte da Krogh (1997). Una delle più importanti è che le regioni di codifica sono modellate da una catena di Markov disomogenea del quarto ordine invece di una catena del primo ordine. Questo è fatto da un’estensione quasi banale del formalismo HMM standard, che permette una catena di Markov di qualsiasi ordine in uno stato del modello, mentre l’HMM standard ha una semplice distribuzione di probabilità incondizionata sulle quattro basi (corrispondente allo 0° ordine). Il modello è frame-aware e può prevedere qualsiasi numero di geni e geni parziali, quindi la struttura generale del modello è come in Figura 11.4 con transizioni aggiunte per permettere l’inizio e la fine negli introni, come in Figura 11.3.
Come già detto, il metodo di stima della massima verosimiglianza funziona bene se la struttura del modello può descrivere le vere statistiche dei geni. Questa è un’assunzione molto idealizzata, e quindi HMMgene usa un altro metodo per stimare i parametri chiamato massima verosimiglianza condizionata (Juang e Rabiner, 1991; Krogh, 1994). In senso lato, la massima verosimiglianza massimizza la probabilità delle sequenze di DNA nel set di allenamento, mentre la massima verosimiglianza condizionata massimizza la probabilità delle strutture geniche di queste sequenze, che, dopo tutto, è ciò che ci interessa. Questo tipo di ottimizzazione è concettualmente simile a quella usata in GeneParser, dove viene ottimizzata anche la precisione della predizione. HMMgene utilizza anche un algoritmo di programmazione dinamica diverso dall’algoritmo di Viterbi per la predizione della struttura del gene. Tutti questi metodi hanno contribuito alle alte prestazioni di HMMgene.
Genie è un altro esempio di modello di stato probabilistico che viene chiamato HMM generalizzato (Kulp et al., 1996; Reese et al., 1997). La Figura 11.4 è in effetti la struttura di stato di Genie, e sia questa figura che la Figura 11.2 sono essenzialmente copiate da Kulp et al. (1996). In Genie, i sensori di segnale (siti di giunzione) e di contenuto (potenziale di codifica) sono reti neurali, e l’output di queste reti è interpretato come probabilità. Questa interpretazione richiede la stima di ulteriori parametri di probabilità che funzionano come pesi sui sensori. Così, anche se è formulato come un modello probabilistico, il problema della ponderazione appare ancora sotto mentite spoglie. L’algoritmo di previsione è quasi identico all’algoritmo di programmazione dinamica dell’ultima sezione. Una versione di Genie include anche le somiglianze del database come parte del sensore di esoni (Kulp et al., 1997).
Ci sono due vantaggi principali delle HMM generalizzate rispetto alle HMM standard. In primo luogo, i singoli sensori possono essere di qualsiasi tipo, come le reti neurali, mentre in una HMM standard sono limitati dalla struttura HMM. In secondo luogo, la distribuzione della lunghezza (per esempio, delle regioni di codifica) può essere presa in considerazione esplicitamente, mentre la distribuzione naturale della lunghezza per una HMM è una distribuzione geometrica, che decade esponenzialmente con la lunghezza. Tuttavia, è possibile avere una modellazione della lunghezza abbastanza avanzata in una HMM se vengono utilizzati diversi stati (Durbin et al., 1998). Il vantaggio di un sistema come HMMgene, d’altra parte, è che si tratta di un unico modello integrato, che può essere ottimizzato tutto in una volta per la massima precisione di predizione.
Un altro gene finder basato su un HMM generalizzato è GENSCAN (Burge e Karlin, 1997). Le principali differenze tra la struttura di stato di GENSCAN e quella di Genie o HMMgene è che GENSCAN modella la sequenza in entrambe le direzioni simultaneamente. In molti gene finder, come quelli descritti sopra, i geni sono prima predetti su un filamento e poi sull’altro. Modellare entrambi i filamenti simultaneamente è stato fatto con molto successo in GeneMark, e un metodo simile è implementato in GENSCAN. Un vantaggio (e forse il principale) è che questa costruzione evita le previsioni di geni sovrapposti sui due filamenti, che presumibilmente sono molto rare nel genoma umano. GENSCAN modella qualsiasi numero di geni e geni parziali come HMMgene. I sensori in GENSCAN sono simili a quelli usati in HMMgene. Per esempio, il sensore di codifica è una catena di Markov disomogenea del quinto ordine. I sensori di segnale sono essenzialmente matrici di peso dipendenti dalla posizione, e quindi sono anche molto simili a quelli di HMMgene, ma ci sono caratteristiche più avanzate nei modelli dei siti di giuntura. GENSCAN modella anche i promotori e il 5′ e 3′ UTR.