EMBL (http://www.ebi.ac.uk/cgi-bin/sva/sva.pl/) è una banca di sequenze nucleotidiche sviluppata da EMBL-EBI (European Bioinformatics Institute, European Molecular Biology Laboratory), in cui ogni sequenza nucleotidica viene memorizzata, con altre informazioni, in file di testo (entry EMBL) in un formato che prende il nome di formato EMBL.
Il formato EMBL è composto da record che iniziano con un codice a due lettere maisucole che specifica il contenuto del record. I soli record che non iniziano con il codice a due lettere sono quelli contenenti la sequenza nucleotidica.
Dato un file in formato EMBL, contenente la sequenza nucleotidica (sequenza di basi) di un mRNA (trascritto espresso da un gene), produrre:
Inoltre, preso in input un file contenente il codice genetico, effettuare la validazione della proteina contenuta nel file EMBL rispetto alla traduzione della sequenza della CDS tramite il codice genetico.
Parametri in input:
Il file del codice genetico è strutturato in record di campi separati da ,, in cui il primo campo è il simbolo di un amminoacido e gli altri campi sono i codoni che codificano quell'amminoacido.
A,gct,gcc,gca,gcg
C,tgt,tgc
D,gat,gac
E,gaa,gag
F,ttt,ttc
G,ggt,ggc,gga,ggg
H,cat,cac
I,att,atc,ata
K,aaa,aag
L,tta,ttg,ctt,ctc,cta,ctg
M,atg
N,aat,aac
P,cct,ccc,cca,ccg
Q,caa,cag
R,cgt,cgc,cga,cgg,aga,agg
S,tct,tcc,tca,tcg,agt,agc
T,act,acc,aca,acg
V,gtt,gtc,gta,gtg
W,tgg
Y,tat,tac
s,tga,taa,tag
Dove trovare le informazioni che servono per risolvere l'esercizio:
Il record che inizia con ID
ID M10051; SV 1; linear; mRNA; STD; HUM; 4723 BP.contiene l'identificatore univoco della sequenza (M10051) e l'organismo (HUM). Il fatto che il file si riferisca alla sequenza nucleotidica di un gene è indicato dalla presenza della parola mRNA.
Il record
FT CDS 139..4287contiene lo start e l'end (1-based) della coding sequence (CDS) sulla sequenza nucleotidica.
L'insieme dei record che iniziano con FT sono quelli che contengono le features della sequenza nucleotidica. In particolare tutti i record della sezione:
FT /translation="MGTGGRRGAAAAPLLVAVAALLLGAAGHLYPGEVCPGMDIRNNLT
FT RLHELENCSVIEGHLQILLMFKTRPEDFRDLSFPKLIMITDYLLLFRVYGLESLKDLFP
FT NLTVIRGSRLFFNYALVIFEMVHLKELGLYNLMNITRGSVRIEKNNELCYLATIDWSRI
FT LDSVEDNHIVLNKDDNEECGDICPGTAKGKTNCPATVINGQFVERCWTHSHCQKVCPTI
FT [...]
FT DGGSSLGFKRSYEEHIPYTHMNGGKKNGRILTLPRSNPS"contengono la sequenza della proteina espressa dal gene.
Il record che inizia con SQ:
SQ Sequence 4723 BP; 1068 A; 1298 C; 1311 G; 1046 T; 0 other;introduce la sezione della sequenza nucleotidica che termina con il record // (file del file). Ogni record contenente la sequenza nucleotidica inizia con una serie di spazi iniziali, e contiene un chunk di sequenza lungo 60 basi. L'intero alla fine del record fornisce la lunghezza totale dei chunks fino a tale record. Ogni chunk in un record viene poi separato in chunks più piccoli di 10 basi.
SQ Sequence 4723 BP; 1068 A; 1298 C; 1311 G; 1046 T; 0 other;
ggggggctgc gcggccgggt cggtgcgcac acgagaagga cgcgcggccc ccagcgctct 60
tgggggccgc ctcggagcat gacccccgcg ggccagcgcc gcgcgcctga tccgaggaga 120
ccccgcgctc ccgcagccat gggcaccggg ggccggcggg gggcggcggc cgcgccgctg 180
ctggtggcgg tggccgcgct gctactgggc gccgcgggcc acctgtaccc cggagaggtg 240
tgtcccggca tggatatccg gaacaacctc actaggttgc atgagctgga gaattgctct 300
gtcatcgaag gacacttgca gatactcttg atgttcaaaa cgaggcccga agatttccga 360
gacctcagtt tccccaaact catcatgatc actgattact tgctgctctt ccgggtctat 420
gggctcgaga gcctgaagga cctgttcccc aacctcacgg tcatccgggg atcacgactg 480
[...]
tttttcgttc cccccacccg cccccagcag atggaaagaa agcacctgtt tttacaaatt 4620
cttttttttt tttttttttt tttttttttg ctggtgtctg agcttcagta taaaagacaa 4680
aacttcctgt ttgtggaaca aaatttcgaa agaaaaaacc aaa 4723
//
NOTA BENE:
t con simboli u.Requisiti:
nell’header della sequenza della CDS in formato FASTA devono comparire l’identificatore univoco della sequenza, l’organismo a cui si riferisce e la lunghezza, nel seguente formato:
>M10051-HUM; len=4149le sequenze in formato FASTA devono essere prodotte in righe di 80 caratteri
non effettuare la sostituzione da t a u.
deve essere definita una funzione format_fasta() che prenda come argomenti un'intestazione FASTA, una sequenza nucleotidica/proteica, e restituisca la sequenza in formato FASTA con la sequenza separata in righe di 80 caratteri.
usare solo espressioni regolari per estrarre le informazioni (tranne per estrarre le informazioni dal file del codice genetico)
Variabili di output:
cds_sequence_fasta: sequenza nucleotidica in formato FASTAcodon_frequency: lista di tuple (codon, frequency) elencate per frequenza decrescenteammino_frequency: lista di tuple (amminoacid, frequency) elencate per frequenza decrescente
In [1]:
def format_fasta(header, sequence):
return header + '\n' + '\n'.join(re.findall('\w{1,80}', sequence))
NOTA BENE: supporre che l'header in input alla funzione non abbia il simbolo newline \n alla fine ma che abbia il simbolo > all'inizio.
In [2]:
genetic_code_name = './genetic-code.txt'
input_file_name = './M10051.txt'
Importazione del modulo re per utilizzare le espressioni regolari (RE).
In [3]:
import re
In [4]:
with open(genetic_code_name, 'r') as genetic_file:
file_str_list = genetic_file.readlines()
In [5]:
file_str_list
Out[5]:
Ottenere da file_str_list il dizionario genetic_code_dict in cui le chiavi sono i codoni e i valori sono il corrispondente amminoacido.
Costruire prima la lista di tuple (chiave, valore) e usare poi la funzione dict() per costruire il dizionario.
In [6]:
key_value_tuple_list = [(codon, row.rstrip().split(',')[0]) for row in file_str_list for codon in row.rstrip().split(',')[1:]]
In [7]:
key_value_tuple_list
Out[7]:
Costruire il dizionario con la funzione dict().
In [8]:
genetic_code_dict = dict(key_value_tuple_list)
In [9]:
genetic_code_dict
Out[9]:
In [10]:
with open(input_file_name,'r') as input_file:
file_str = input_file.read()
In [11]:
file_str
Out[11]:
Estrarre dal record ID l'identificatore univoco e l'organismo nelle variabili identifier e organism.
ID M10051; SV 1; linear; mRNA; STD; HUM; 4723 BP.
Estrazione dell'organismo.
In [12]:
s = re.search('([\w\s]+;){5}\s+(\w+);', file_str, re.M)
organism = s.group(2)
In [13]:
organism
Out[13]:
Estrazione dell'identificatore.
In [14]:
s = re.search('^ID\s+(\w+);', file_str, re.M)
identifier = s.group(1)
In [15]:
identifier
Out[15]:
Estrarre nella lista seq_row_list i record della sequenza nucleotidica escludendo solo l'intero finale (mantenendo gli spazi iniziali e gli spazi prima dell'intero).
tgggggccgc ctcggagcat gacccccgcg ggccagcgcc gcgcgcctga tccgaggaga 120
In [16]:
seq_row_list = re.findall('^\W{2}(\D+)\d+', file_str, re.M)
In [17]:
seq_row_list
Out[17]:
Estrarre da seq_row_list la lista seq_chunk_list contenente i chunks di (al più) lunghezza 10 della sequenza nucleotidica.
NOTA BENE: l'elemento seq_chunk_list[i] è una lista annidata e contiene i sei chunks relativi all'i-esimo record di seq_row_list.
In [18]:
seq_chunk_list = [re.findall('\w+', row) for row in seq_row_list]
In [19]:
seq_chunk_list
Out[19]:
Concatenare i chunks della lista seq_chunk_list per ottenere la sequenza nucleotidica nella variabile nucleotide_sequence.
In [20]:
nucleotide_sequence = ''.join([''.join(list_six_chunks) for list_six_chunks in seq_chunk_list])
In [21]:
nucleotide_sequence
Out[21]:
Estrarre nella variabile protein_prefix il prefisso della proteina contenuto nel record contenente la parola /translation:
FT /translation="MGTGGRRGAAAAPLLVAVAALLLGAAGHLYPGEVCPGMDIRNNLT
In [22]:
s = re.search('^FT\s+\/translation=\"(\w+)$', file_str, re.M)
protein_prefix = s.group(1)
In [23]:
protein_prefix
Out[23]:
Estrarre nella lista protein_row_list gli altri record (compreso l'ultimo) che contengono la sequenza della proteina.
FT RLHELENCSVIEGHLQILLMFKTRPEDFRDLSFPKLIMITDYLLLFRVYGLESLKDLFP
Attenzione all'ultimo:
FT DGGSSLGFKRSYEEHIPYTHMNGGKKNGRILTLPRSNPS"
In [24]:
protein_row_list = re.findall('^FT\s+([ACDEFGHIKLMNPQRSTVWY]+)\"?$', file_str, re.M)
In [25]:
protein_row_list
Out[25]:
Aggiungere in testa alla lista il prefisso trovato prima e concatenare nella variabile protein_sequence tutti i blocchi della lista protein_row_list per ottenere la sequenza della proteina.
In [26]:
protein_row_list[:0] = protein_prefix
protein_sequence = ''.join([''.join(chunk) for chunk in protein_row_list])
In [27]:
protein_sequence
Out[27]:
Estrarre dal record
FT CDS 139..4287
lo start e l'end della CDS sulla sequenza nucleotidica dell'mRNA.
In [28]:
s = re.search('^FT\s+CDS\s+(\d+)\.\.(\d+)$', file_str, re.M)
cds_start = s.group(1)
cds_end = s.group(2)
In [29]:
cds_start
Out[29]:
In [30]:
cds_end
Out[30]:
Ottenere la coding sequence (CDS)
In [31]:
cds_sequence = nucleotide_sequence[int(cds_start)-1:int(cds_end)]
In [32]:
cds_sequence
Out[32]:
Produrre nella variabile cds_sequence_fasta la sequenza della CDS in formato FASTA con il seguente header:
>M10051-HUM; len = 4149
In [33]:
header = '>' + identifier + '-' + organism + '; len = ' + str(len(cds_sequence))
cds_sequence_fasta = format_fasta(header, cds_sequence)
In [34]:
print(cds_sequence_fasta)
Ottenere in codon_list la lista dei codoni (stop codon escluso) della coding sequence.
In [35]:
codon_list = re.findall('\w{3}', cds_sequence)[:-1]
In [36]:
codon_list
Out[36]:
Costruire la lista codon_frequency di tuple (codon, frequency) elencate per frequenza decrescente.
In [37]:
from collections import Counter
codon_counter = Counter(codon_list)
codon_frequency = codon_counter.most_common()
In [38]:
codon_frequency
Out[38]:
Ottenere in ammino_list la lista degli amminoacidi della proteina.
In [39]:
ammino_list = re.findall('\w', protein_sequence)
In [40]:
ammino_list
Out[40]:
Costruire la lista ammino_frequency di tuple (amminoacid, frequency) elencate per frequenza decrescente.
In [41]:
from collections import Counter
ammino_counter = Counter(ammino_list)
ammino_frequency = ammino_counter.most_common()
In [42]:
ammino_frequency
Out[42]:
Ottenere nella variabile cds_translation la traduzione della CDS tramite il codice genetico.
In [43]:
cds_translation = ''.join([genetic_code_dict[codon] for codon in codon_list])
In [44]:
cds_translation
Out[44]:
Verificare che la proteina estratta dal file in input sia uguale a quella ottenuta dalla traduzione della coding sequence dell'mRNA.
In [45]:
cds_translation == protein_sequence
Out[45]: