Argomenti Trattati:

[1] Biosequenze
[2] Parsing

Parte 2 - PARSING
Autore: Davide FARO'

PREFAZIONE

Spesso in Bioinformatica si ha la necessità di estrarre particolari informazioni da files depositati in banca dati (entry). In base alla tipologia di informazione trattata abbiamo a che fare con un tipo di file piuttosto che con un altro, per esempio: le sequenze di DNA, sono strutturate in files con estensioni (.fasta) o (.gbk), le coordinate atomiche di una determinata proteina stanno in un file con estensione (.pdb) etc. Ognuno di questi files ha una sua precisa struttura con cui l’informazione è organizzata al suo interno. L'estrazione di specifiche informazioni da un file viene generalmente detto parsing ed il parser è un programma che esegue questo compito. Fare il parsing di opportune informazioni, dai files contenuti nelle più comuni banche dati, è relativamente semplice grazie a Biopython. La libreria mette infatti a disposizione dei programmatori il modulo Bio.SeqIO contenente tutte le classi necessarie allo scopo. Il modulo Bio.SeqIO è disponibile dalla versione Biopython 1.4 e successive.

PARSING DI FILE FASTA

Useremo come primo esempio, un file FASTA generico di nome "ls_orchid.fasta". Se aperto con qualsiasi editor di testo questo file si presenta così:

>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG
AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG
...

I puntini di sospensione indicano che ciò che è visualizzato non è completo. Il file completo ha infatti 94 record e comincia con il carattere ">", tipico del formato FASTA. Estraiamo le informazioni contenute in questo file utilizzando il modulo Bio.SeqIO.

from Bio import SeqIO
handle = open("ls_orchid.fasta")
for seq_record in SeqIO.parse(handle, "fasta") :
   print seq_record.id
   print seq_record.seq
   print len(seq_record.seq)
handle.close()

L’esecuzione del codice precedente legge il file FASTA, e genera il seguente output:

gi|2765658|emb|Z78533.1|CIZ78533
Seq(’CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA...’, SingleLetterAlphabet())
740
...
gi|2765564|emb|Z78439.1|PBZ78439
Seq(’CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGTTTACTTTGGTC...’, SingleLetterAlphabet())
592

Il "ciclo for" estrae l’ id di ogni record, la sequenza contenuta in ogni record, e la lunghezza di quest’ultima, fino alla fine del file. Come si nota, la seconda riga di ogni output è un oggetto sequenza, visto nella prima parte di questo tutorial, ma l’oggetto alphabet non è specificato, perchè il file FASTA non specifica un alfabeto, così SeqIO ne ha uno di default: SingleLetterAlphabet().

PARSING DI FILE GENBANK

from Bio import SeqIO
handle = open("ls_orchid.gbk")
for seq_record in SeqIO.parse(handle, "genbank") :
   print seq_record.id
   print seq_record.seq
   print len(seq_record.seq)
handle.close()

L’ output:

Z78533.1
Seq(’CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA...’, IUPACAmbiguousDNA())
740
...
Z78439.1
Seq(’CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGTTTACTTTGGTC...’, IUPACAmbiguousDNA())
592

Abbiamo estratto le stesse informazioni estratte dal file FASTA e cioè l’ id (una stringa molto più breve rispetto a quella del file FASTA), la sequenza e la relativa lunghezza, nell’oggetto sequenza questa volta è presente l’oggetto alphabet: IUPACAmbiguousDNA().

ITERAZIONE RIGA DOPO RIGA

Negli esempi precedenti abbiamo utilizzato un ciclo for ed estratto pezzi specifici di file, ma possiamo anche iterare rigo dopo rigo dentro un file scorrendo in avanti il puntatore. La classe Bio.SeqIO utilizza un puntatore che punta ad un record, con il metodo next() possiamo spostarci dal record corrente a quello successivo, vediamo un esempio:

from Bio import SeqIO
handle = open("ls_orchid.fasta")
record_iterator = SeqIO.parse(handle, "fasta")
first_record = record_iterator.next()
print first_record.id
print first_record.description
handle.close()

Una delle strutture dati con cui tutti i programmatori python sono abituati a lavorare è senza dubbio la lista, i dati che stanno dentro uno dei file visti prima possono essere inseriti in una lista e trattati di conseguenza. La funzione list() di python mette le righe di un file (oggetti SeqRecord), in una lista, vediamo come.

from Bio import SeqIO
handle = open("ls_orchid.gbk")
records = list(SeqIO.parse(handle, "genbank"))
handle.close()
print "Found %i records" % len(records)
print "The last record"
last_record = records[-1] #sintassi di python per referenziare l’ultimo elemento della lista
print last_record.id
print last_record.seq
print len(last_record.seq)
print "The first record"
first_record = records[0] #le liste indicizzano partendo da zero
print first_record.id
print first_record.seq
print len(first_record.seq)

L’output generato.

Found 94 records
The last record
Z78439.1
Seq(’CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGTTTACTTTGGTC...’, IUPACAmbiguousDNA())
592
The first record
Z78533.1
Seq(’CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA...’, IUPACAmbiguousDNA())
740

Un’altra struttura dati potente di python è il Dizionario, possiamo inserire i dati contenuti nel nostro file dentro un dizionario e costruirci così un piccolo database in memoria, che possiamo interrogare usando delle chiavi in maniera molto semplice ed efficiente. Utilizzando la funzione SeqIO.to_dict() otteniamo un dizionario di oggetti SeqRecord in memoria, per impostazione predefinita ogni record sarà identificato dall’attributo "id" dell’oggetto SeqRecord, vediamo un esempio con il file GENBANK.

from Bio import SeqIO
handle = open("ls_orchid.gbk")
orchid_dict = SeqIO.to_dict(SeqIO.parse(handle, "genbank"))
handle.close()

La variabile orchid_dict è un dizionario, vediamo quali sono le chiavi.

>>> print orchid_dict.keys()
[’Z78484.1’, ’Z78464.1’, ’Z78455.1’, ’Z78442.1’, ’Z78532.1’, ’Z78453.1’, ..., ’Z78471.1’]

A questo punto possiamo utilizzare queste chiavi per accedere ad un oggetto sequenza nel dizionario.

>>> seq_record = orchid_dict["Z78475.1"]
>>> print seq_record.description
P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA
>>> print seq_record.seq
Seq(’CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACATAATAAT ...’, IUPACAmbiguousDNA())

In modo altrettanto semplice si può fare lo stesso con un file FASTA, vediamo un esempio.

from Bio import SeqIO
handle = open("ls_orchid.fasta")
orchid_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta"))
handle.close()
print orchid_dict.keys()

La funzione to_dict() restituisce un dizionario di oggetti sequenza, le chiavi, dopo aver invocato la funzione print, sono:

[’gi|2765596|emb|Z78471.1|PDZ78471’, ’gi|2765646|emb|Z78521.1|CCZ78521’, ...
..., ’gi|2765613|emb|Z78488.1|PTZ78488’, ’gi|2765583|emb|Z78458.1|PHZ78458’]

Nel caso di un file FASTA, le chiavi sono stringhe abbastanza lunghe e strutturate, se volessimo rendere le chiavi di ricerca meno complesse, si potrebbe utilizzare una parte della chiave originale, il numero di adesione. Possiamo farlo scrivendoci da soli una funzione che dal record del dizionario ci tiri fuori la chiave da noi scelta. Nell’esempio, la funzione get_accession() da noi definita, usa la funzione built-in di python, split() che, preso in input un carattere, spezza la stringa in pezzi usando il carattere "|", vediamo un esempio.

def get_accession(record) :
   """Da un SeqRecord, ritorna un numero di adesione come una stringa
   e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1" """
   parts = record.id.split("|")
   assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"
   return parts[3]

Usiamo la nostra funzione.

from Bio import SeqIO
handle = open("ls_orchid.fasta")
orchid_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta"),key_function=get_accession)
handle.close()
print orchid_dict.keys()

L’ output sarà:

>>> print orchid_dict.keys()
[’Z78484.1’, ’Z78464.1’, ’Z78455.1’, ’Z78442.1’, ’Z78532.1’, ’Z78453.1’, ..., ’Z78471.1’]

Come vedete adesso le chiavi non sono più delle stringhe lunghe e complesse ma stringhe brevi e coincise. Supponiamo ora di voler estrarre una informazione specifica da un record, per esempio il nome della specie a cui appartiene la sequenza di DNA contenuta nel file GENBANK, queste informazioni sono registrate nel primo record, vediamo.

from Bio import SeqIO
record_iterator = SeqIO.parse(open("ls_orchid.gbk"), "genbank")
first_record = record_iterator.next()
print first_record

L’ output:

ID: Z78533.1
Name: Z78533
Desription: C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA.
source=Cypripedium irapeanum
taxonomy=[’Eukaryota’, ’Viridiplantae’, ’Streptophyta’, ..., ’Cypripedium’]
keywords=[’5.8S ribosomal RNA’, ’5.8S rRNA gene’, ’internal transcribed spacer’, ’ITS1’, ’ITS2’]
references=[...]
data_file_division=PLN
date=30-NOV-2006
organism=Cypripedium irapeanum
gi=2765658
Seq(’CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA...’, IUPACAmbiguousDNA())

L’informazione che ci interessa si trova nelle annotazioni, accessibile come in un dizionario sotto le voci “source” e “organism” evidenziati in rosso per metterle in risalto, accediamo ad esse in questo modo.

print first_record.annotations["source"]

Oppure:

print first_record.annotations["organism"]

Adesso creiamo una lista dei nomi delle specie di ogni sequenza contenuta nel file e la stampiamo.

from Bio import SeqIO
handle = open("ls_orchid.gbk")
all_species = []
for seq_record in SeqIO.parse(handle, "genbank") :
all_species.append(seq_record.annotations["organism"])
handle.close()
print all_species

Dalla versione 2.0 di python è possibile ottenere lo stesso risultato con una sintassi più concisa.

from Bio import SeqIO
all_species == [seq_record.annotations["organism"] for seq_record in \
SeqIO.parse(open("ls_orchid.gbk"), "genbank")]
print all_species

L’ output che otterremo in entrambi i casi sarà:

[’C.irapeanum’, ’C.californicum’, ’C.fasciculatum’, ’C.margaritaceum’, ..., ’P.barbatum’]

Con i file GENBANK questa operazione è piuttosto semplice perché sono annotati in modo standard, ma se volessimo fare lo stesso con un file FASTA? Basta spezzare il record descrizione usando i caratteri di spazio bianco ed ottenere una struttura dati tipo array o lista, l’ informazione che ci serve sta nel campo uno, vediamo un esempio.

from Bio import SeqIO
handle = open("ls_orchid.fasta")
all_species = []
for seq_record in SeqIO.parse(handle, "fasta") :
all_species.append(seq_record.description.split()[1])
handle.close()
print all_species

In modo più conciso:

from Bio import SeqIO
all_species == [seq_record.description.split()[1] for seq_record in \
SeqIO.parse(open("ls_orchid.fasta"), "fasta")]
print all_species

In generale, estrarre informazioni dal record descrizione FASTA non è facile, si preferisce farlo da file in formato GENBANK o EMBL, che conservano queste info strutturate in modo standard.

[Torna su]