PAOLO CAMPIGOTTO - MATTEO PEDROTTI
Prefazione
L'idea di intraprendere uno studio coordinato del genoma umano sorse tra il 1985 e il 1987, nel corso di una serie di congressi scientifici. L'avviamento effettivo del Progetto genoma umano risale, invece, al 1990, grazie ai fondi messi a disposizione negli Stati Uniti da due istituzioni governative, i National Institutes of Health (NIH) e il Departement of Energy (DOE).
Il "Genoma Umano" contiene tutte le istruzioni biochimiche sotto forma di combinazioni delle quattro basi del DNA: adenina, timina, citosina, guanina (A,T,C,G,) per costruire e mantenere in vita un essere umano. Nel 1986 venne varata dalla comunità scientifica internazionale una delle tappe più ambiziose della ricerca biologica: la possibilità di identificare la sequenza nucleotidica del genoma (cioè l'insieme dei geni presenti nel corredo cromosomico di una specie).
Il genoma è la collezione completa di tutto il materiale genetico di un organismo vivente. Il genoma umano è composto da circa 50.000-100.000 geni, localizzati sulle 23 coppie di cromosomi presenti nel nucleo di ciascuna cellula umana. Un singolo cromosoma può arrivare a contenere più di 250 milioni di coppie di basi, mentre l'intero genoma umano, secondo alcune stime, è costituito da circa 3 miliardi di coppie di basi.
L'informatica svolge due servizi importanti per la mappatura del genoma umano:
- analisi dei risultati di laboratorio
- registrazione dei dati in database globali
Il progetto proposto appartiene al primo tipo di servizio. In particolare intende fornire al ricercatore le occorrenze di 8 basi in un dato tratto genomico.
I problemi che sorgono durante l'analisi sono legati ad imperfezioni nelle fasi di laboratorio: è infatti possibile che venga riconosciuta una base invece che un'altra.
A conoscenza di questi rischi, il progetto proponeva, non solo di separare nel tratto genomico le sequenze di 8 basi, ma anche di ipotizzare che una base all'interno della stringa sia stata rilevata in maniera erronea; questo significa che per ogni base della stringa viene calcolata l'occorrenza di 4 possibili stringhe.
In questo modo il ricercatore ha a disposizione la probabilità del numero di occorrenze di una stringa e non più un numero preciso che potrebbe rivelarsi errato.
Programma
Il linguggio utilizzato è stato il Pascal standard.
Si è pensato di implementare le funzionalità in modalità batch in modo che una volta avviato il programma esso produca automaticamente i file di testo contenenti tutte le analisi possibili (2 file standard e da 0 a 16 file jolly).
Il programma sovrascrive ad ogni esecuzione i vecchi file di output.
L'eseguibile è avviato tramite il comando
genoma < file_genoma> [< posizione_jolly>]
dove < file_genoma> è il file che viene letto e < posizione_jolly> e' la posizione del jolly.
Se la posizione_jolly e' 0 nessun file jolly verra' stampato, se e' 9 verranno stampati tutti gli 8 file jolly, infine un numero x compreso tra 1 e 8 rappresentera' la posizione del jolly voluta.
NOTE AGGIUNTIVE: visto l'utilizzo di file di grosse dimensioni, e vista la limitatezza del buffer utilizzabile in Pascal (255 caratteri) si è deciso di stampare a video i MegaByte letti fino a quel momento, in modo che l'utente possa rendersi conto del punto di avanzamento nel file.
I calcoli richiesti influenzano in maniera minima la velocità di lettura (si parla di milionesimi di secondo per le operazioni e di millesimi per l'accesso al file su disco rigido) e forniscono un servizio importante all'utente.
Descrizione delle strutture dati utilizzate
Il progetto e' stato realizzato sfruttando la strutturazione del linguaggio Pascal tramite la creazione di un "main program" che organizza le chiamate alle procedure/funzioni che implementano le varie "azioni" richieste nella traccia.
Per rappresentare i dati ricavati dai files di input si sono utilizzate le seguenti strutture:
- un albero di grandezza standard (poiche' la radice e' un nodo fittizio 8 livelli "effettivi" per 65536 nodi complessivi), chiamato "albero alfabeto", che rappresenta il "dominio" del nostro problema, ovvero tutte le parole di otto lettere che si possono creare a partire dall'insieme {A,C,G,T}.
Ciascun nodo rappresenta un carattere dell'insieme e contiene:
- 4 puntatori ("figlioA","figlioC","figlioG","figlioT") che lo collegano al carattere successivo della parola; questi puntatori sono ordinati da sinistra a destra in base al carattere che contengono (ordine alfabetico);
- un campo "frequenza" che contiene la frequenza con cui compare nel file di input la stringa che si ottiene eseguendo il percorso dalla radice al nodo in questione. Ovviamente questo campo e' incrementato ogni volta che si attraversa il nodo.
L'albero cosi' realizzato e'dotato di un primo nodo fittizio che serve unicamente a "legare" i 4 possibili caratteri iniziali
di una parola.
Chiaramente una parola di otto lettere si ottiene eseguendo un percorso completo radice-foglia e la sua frequenza si ricava leggendo il campo frequenza della radice.
Poiche' nel progetto si richiede di lavorare unicamente con parole di 8 lettere, sembrerebbe uno spreco tenere un campo frequenza anche per i nodi intermedi; in verita' questa nostra scelta deriva dalla volonta' di conferire maggior "flessibilita'" al programma, cosicche', se in futuro si volesse lavorare con parole di 1..7 lettere, basterebbe una semplice modifica al codice esistente.
E' interessante notare come i nodi non possiedano un campo contenente il carattere che rappresentano: esso e'dedotto dal nome del puntatore presente nel padre del nodo stesso.
- 2 alberi binari di grandezza variabile, detti "alberi frequenza", che memorizzano le parole trovate nel file di input in maniera ordinata in base alla frequenza. Il primo albero e' l'albero contenente le stringhe complete di 8 basi, mentre il secondo albero contiene le stringhe con una base jolly. L'albero jolly e' appositamente creato per una data posizione del jolly percorrendo in maniera opportuna l'"albero alfabeto", e pertanto viene distrutto da un'apposita funzione ogni qual volta se ne e'stampato il contenuto in un file, per un totale di 8 alberi-jolly (creati e distrutti)
Ogni nodo contiene un campo parola, un campo frequenza e i due puntatori ai figli.
Descrizione del codice
Ciascuna azione che si compie sulle strutture sopra descritte e' implementata da una apposita funzione/procedura.
Ecco le principali funzioni usate per lavorare nell'albero alfabeto:
- inizializza_alfabeto(var radice: albero_alfabeto; livello: integer);
Nel "main" viene richiamata passando la radice dell'albero alfabeto, che sara' un puntatore cui si e' gia' assegnata una porzione di memoria tramite la "new", e il valore zero al parametro livello. Infatti l'albero e' formato da 9 livelli (0..8), di cui il primo (livello 0) e' quello del nodo fittizio.
La procedura e' ricorsiva: dopo aver settato a zero il campo frequenza del nodo corrente,se quest'ultimo non e' una foglia,si allocano i puntatori ai 4 figli e si richiama su essi la procedura stessa, altrimenti si settano a NIL i 4 puntatori.
- procedure aggiungi_frequenza_alfabeto(var radice: albero_alfabeto; var parola: stringa; livello:integer);
Questa procedura aggiorna l'albero dell'alfabeto quando nel testo di input si e' trovata una nuova parola di 8 lettere.
Essa viene chiamata nel "main" passando in ingresso la radice dell'albero "alfabeto", il valore uno al parametro "livello" e al parametro "parola" la parola individuata nel testo che rappresenta il percorso che si dovra' effettuare nell'albero stesso.
Al parametro "livello" viene passato il valore 1 sia perche' il primo livello dell'albero (livello 0) e' fittizio sia perche' il questo parametro fungera' da indice della stringa parola che va da 1 a 8.
La procedura e' ricorsiva: se non ci si trova nella foglia si richiama la procedura per i figli ("figlioA","figlioC","figlioG","figlioT"), altrimenti si itera il valore del nodo.
- procedure leggi_file(nome_ingresso: string);
Vengono letti maxbuffer (255) caratteri dal file e li mette in un buffer.
Poi si analizza carattere per carattere il buffer per tutta la sua lunghezza.
Se il carattere è valido (ossia se è in "ACGTacgt") i caratteri della stringa vengono spostati a sinistra di una posizione e il carattere appena letto viene inserito in coda.
(per es: ACCGTGAA è la parola già in stringa; leggiamo dal buffer la lettera G; spostiamo a sinistra i caratteri della stringa e aggiungiamo in coda la G: otteniamo CCGTGAAG.)
A questo punto si inserisce la parola in albero_alfabeto.
NOTE: i caratteri in minuscolo vengono automaticamente trasformati in maiuscolo.
Durante la lettura viene calcolato l'ammontare di byte letti.
- procedure stampa_alfabeto(var radice: albero_alfabeto; var parola: stringa; livello:integer; var canale_uscita: text);
Stampa sul file canale_uscita tutte le parole con frequenza maggiore di zero. Scorre l'albero alfabetico radice fino ad arrivare ad una foglia; mentre lo scorre aggiunge nella parola la lettera caratteristica di ogni nodo. Al termine dello scorrimento controlla se la frequenza è maggiore di zero: se sì la parola viene scritta sul file. Lo scorrimento viene effettuato esaminando il ramo A poi il ramo C, G e infine T. In questo modo le parole sono lette in ordine alfabetico.
La procedura è ricorsiva, la posizione in cui inserire il carattere è fornita da livello e la condizione di uscita è data dal controllo sul figlio della radice.
Ecco le principali funzioni usate per lavorare nell'albero frequenza:
- procedure inizializza_frequenza(var radice: albero_frequenza);
Banale,vedasi direttamente il codice.
- procedure aggiungi_parola_frequenza(var radice: albero_frequenza; var parola: stringa; var frequenza: integer);
Questa procedura lavora nell'albero delle frequenze.Viene richiamata nel "main" passando in ingresso la radice dell'albero delle frequenze, la parola da inserire e la relativa frequenza. L'albero e' binario e ordinato: il figlio sinistro contiene una parola con frequenza minore e il figlio destro una parola con frequenza maggiore della stringa contenuta nel nodo stesso.
Pertanto nell'aggiungere una nuova parola all'albero si confronta la sua frequenza con la frequenza della parola contenuta nel nodo corrente e si decide se richiamare la procedura sul figlio destro o sinistro; la ricorsione termina non appena non esiste il figlio scelto: in questo caso lo si creera' ex-novo e vi si inseriranno gli opportuni dati.
Se la frequenza della parola che si vuole inserire e' uguale a quella della stringa contenuta nel nodo corrente la parola va inserita scendendo verso sinistra nel sottoalbero radicato nel figlio sinistro e fermandosi appena si trova una parola di frequenza minore.
Cio' e' dovuto alla volonta' di produrre, in fase di output, una stampa in ordine alfabetico delle parole che hanno la stessa frequenza.
- procedure costruisci_frequenza_alfabeto(var radice: albero_alfabeto; var radice_frequenza: albero_frequenza; var parola: stringa; livello:integer);
Questa procedura costruisce l'albero delle frequenze ricavando dall'albero alfabeto tutte le parole di otto lettere di frequenza non nulla.
Per cercare i dati necessari alla costruzione dell'albero delle frequenze essa esamina tutti i percorsi che portano a ciascuna foglia dell'albero alfabeto e, se quest'ultima non contiene un campo frequenza nullo, inserisce la frequenza letta e il percorso effettuato nell'albero delle frequenze chiamando la procedura "aggiungi_parola_frequenza".
- procedure stampa_frequenza(var radice: albero_frequenza; var canale_uscita: text);
Stampa sul file tutte le parole in ordine di frequenza.
L'albero frequenza viene scorso prima nel ramo destro (cioè verso frequenze maggiori).
Arrivato in fondo al ramo destro la parola viene scritta sul file, poi ricorsivamente passa ad esaminare il ramo sinistro.
Il ramo sinistro comprende, oltre che le frequenze minori anche quelle uguali in ordine alfabetico. In questo modo le parole vengono scritte sul file in ordine decrescente di frequenza e, a parità di frequenza, in ordine alfabetico.
Infine si descrivono le procedure/funzioni usate per elaborare i dati in modalita' "jolly":
- procedure converti_e_copia(var stringa1: percorso; var stringa2: stringa; pos: integer);
Copia il percorso dato come stringa1 in stringa2 a partire dalla posizione pos trasformando i numeri (0..3) in lettere (A..T).
- procedure inizializza_percorso;
Banale. Imposta il percorso_corrente ad AAAAAAA.
- procedure stampa_percorso;
Stampa il percorso_corrente dopo averlo trasformato in stringa.
- function fornisci_percorso (lunghezza: integer): boolean;
Prende il percorso corrente (globale) e dalla lunghezza specificata in poi lo incrementa (a
Per farlo considera il percorso come un numero in base 4 (0>1>2>3) e somma 1 all'ultimo carattere. La somma è un ciclo for che parte dall'ultimo carattere del percorso e continua fino alla lunghezza specificata. Il ciclo parte considerando che ci sia un riporto da sommare al numero in ultima posizione. Poi, nel caso che il numero sia 3 (sommando 1 diventa 0) il riporto rimane 1, altrimenti il riporto diventa 0. Il ciclo continua a sommare il riporto fino a lunghezza. Se c'è un riporto anche al primo numero (posizione lunghezza), cioè se si verifica una specie di overflow, la funzione restituisce false, altrimenti true.
In questo modo il percorso_corrente viene modificato in modo che punti al ramo successivo, cioè quello più a destra.
- function esegui_percorso (var radice: albero_alfabeto; pos: integer):integer;
Legge il percorso corrente (globale) e scorre l'albero_alfabeto (radice) seguendo il percorso indicato dalla posizione pos. Arrivato in fondo prende il valore della foglia e lo restituisce.
- procedure esamina_sottoalbero (var radice:nodo_alfabeto; var albero_frequenza_jolly: albero_frequenza; var posiz_jolly, livello:integer; var uscita_alfabeto: text; var parola: stringa);
Questa procedura lavora nei 4 sottoalberi radicati nel nodo corrente dell'albero alfabeto, rappresentato dal parametro "radice": chiamando la funzione fornisci_percorso prende un pattern di lunghezza 8-jolly, lo rintraccia nei quattro sottoalberi radicati in figlioA, figlioC, filglioG, figlioT, computa la somma delle frequenze estratte dalle foglie ed inserisce nell'albero frequenza jolly il risultato della computazione e la stringa contenuta nell'array parola.
Chiaramente e' trattato in maniera particolare il caso in cui il jolly si trovi in ottava posizione; in questo caso non si chiama la funzione fornisci_percorso perche' il pattern sarebbe lungo 8-posizione del jolly,cioe' zero!
Infatti i 4 sottoalberi radicati nei figli sono degeneri poiche' questi ultimi sono delle foglie. Ci si limitera' dunque a sommare le frequenze delle 4 foglie e, se la somma non e' nulla, a inserire nell'albero frequenza jolly i risultati ottenuti.
- procedure conteggia_in_albero_alfabetico (var radice:albero_alfabeto;var albero_frequenza_jolly:albero_frequenza; var parola: stringa; posiz_jolly,livello: integer; var uscita_alfabeto: text);
Questa procedura e' molto delicata: permette di elaborare i dati in modalita' "jolly".
E' chiamata ricorsivamente sui nodi di "albero alfabeto" obbedendo al seguente schema logico: se il livello del nodo che sto esaminando e' minore del livello precedente a quello del "jolly" allora chiama ricorsivamente la procedura sui 4 figli provvedendo a tener traccia nella stringa "parola" del carattere delrappresentato dal figlio scelto;
altrimenti, se mi trovo appena al di sopra del livello del "jolly" chiama la procedura esamina_sottoalbero.
Mentre conteggia stampa su uscita le stringhe trovate.
- procedure distruggi_albero_jolly (var radice: albero_frequenza);
Cancella ricorsivamente l'albero frequenza.
La procedura e' chiamata solo per distruggere l'albero jolly ma funziona per qualsiasi radice di albero_frequenza passata.
- procedure salva_jolly (i: integer);
Stampa in "alfabeto_jollyX.txt" in ordine alfabetico e in "frequenza_jollyX.txt" le stringhe con jolly in posizione specificata i.
Il file alfabeto_jollyX.txt viene stampato automaticamente chiamando conteggia_in_albero_alfabeto, mentre per stampare la frequenza viene chiamata la procedura stampa_frequenza, la stessa utilizzata per il file senza jolly.
Al termine della stampa l'albero creato da conteggia_in_albero_alfabeto viene distrutto con la procedura distruggi_albero_jolly e poi reinizializzato.
MAIN PROGRAM
1. Inizializzazione:
- stabilito il nome del file contenente il tratto genomico
Se il file non e' indicato come argomento viene stampata la schermata di errore e il programma viene fermato.
- stabilita la posizione del jolly
Se il secondo argomento non e' presente, il jolly diventa 9 ossia vengono creati tutti i file jolly.
Se il secondo argomento e' valido viene analizzato il primo carattere:
- se e' un numero valido jolly assume quel valore
- altrimenti jolly diventa 0.
- stabilita la grandezza del buffer di accesso al file (da 1 a 255)
maxbuffer:=255; {1..255}
- settato il livello di partenza degli alberi a zero
livello:=0;
- inizializzato l'albero che contiene le stringhe in ordine alfabetico
new (radice_alfabeto);
inizializza_alfabeto(radice_alfabeto, livello);
- inizializzato l'albero che contiene le stringhe in ordine di frequenza
new (radice_frequenza);
inizializza_frequenza(radice_frequenza);
- inizializzato l'albero che contiene le stringhe con una lettera jolly in ordine di frequenza
new (radice_jolly);
inizializza_frequenza(radice_jolly);
2. Caricamento:
- si legge il file di input e si mettono le stringhe di 8 lettere trovate nell'albero alfabetico
leggi_file(nome_file);
3. Elaborazione:
- basandosi sull'albero alfabetico viene costruito l'albero frequenza
costruisci_frequenza_alfabeto(radice_alfabeto, radice_frequenza, parola, 1);
4. Stampa file senza jolly:
- apre il file chiamato per default "alfabeto.txt"
open (uscita, 'alfabeto.txt', NEW);
rewrite(uscita);
- stampa sul file l'albero alfabetico
stampa_alfabeto(radice_alfabeto, parola, 1, uscita);
- chiude il file "alfabeto.txt"
close(uscita);
- apre il file chiamato per default "frequenza.txt"
open (uscita, 'frequenza.txt', NEW);
rewrite(uscita);
- stampa sul file l'albero ordinato per frequenza
stampa_frequenza(radice_frequenza^.dx, uscita);
- chiude il file "frequenza.txt"
close(uscita)
5. Elaborazione e stampa file jolly
- controlla il valore della variabile jolly.
Se jolly e' 0 salta tutto e termina il programma.
Se jolly e' 9 chiama un ciclo for per calcolare tutti i jolly, altrimenti
salva_jolly(jolly);
- ciclo compiuto otto volte, una per ogni posizione del jolly nella stringa di 8 caratteri
for i:=1 to 8 do
- viene chiamata la procedura salva_jolly passando la posizione del jolly da calcolare.
salva_jolly(i);