<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Esercizio 4
EMBL (http
Step1:  NOTA BENE
Step2:  Importazione del modulo re per utilizzare le espressioni regolari.
Step3:  Lettura del file del codice genetico in una lista di righe
Step4:  Costruzione del dizionario del codice genetico
chiave
Step5:  Costruire il dizionario con la funzione dict().
Step6:  Lettura del file EMBL in un'unica stringa
Step7:  Estrazione dell'identificatore univoco e dell'organismo relativo all'entry.
Estrarre dal record ID
Step8:  Estrazione dell'organismo.
Step9:  Estrazione della sequenza nucleotidica
Ottenere la lista dei record della sequenza nucleotidica escludendo da ognuno di essi solo l'intero finale.
Step10:  Estrarre dalla lista appena costruita, la lista contenente i pezzi di lunghezza 10 della sequenza nucleotidica in modo tale che l'elemento i-esimo sia una lista annidata contenente i sei pezzi di sequenza relativi all'i-esimo record della lista appena costruita.
NOTA BENE
Step11:  Concatenare i pezzi della lista in un'unica stringa.
Step12:  Estrazione della sequenza della proteina
Estrarre il prefisso della proteina contenuto nel record
Step13:  Ottenere la lista degli altri record di proteina
Step14:  Aggiungere in testa alla lista il prefisso trovato prima e concatenare tutti gli elementi della lista per ottenere la sequenza della proteina in un'unica stringa.
Step15:  Coding sequence (CDS) in formato FASTA
Estrarre dal record
FT   CDS             139..4287
lo start e l'end (1-based) della CDS.
Step16:  Ottenere la sequenza della CDS.
Step17:  Produrre la sequenza della CDS in formato FASTA con il seguente header
Step18:  Distribuzione di frequenza dei codoni della CDS
Estrarre la lista dei codoni della CDS.
Step19:  In alternativa (by Andrea Vignati), usare la funzione range() per creare la lista delle posizioni di inizio 0, 3, 6, 9, 12, ... dei codoni e usare una list comprehension che estragga i codoni sulla base di tali posizioni.
Step20:  Costruire la lista di tuple (codone, frequenza) elencate per frequenze decrescenti.
Step21:  Distribuzione di frequenza degli amminoacidi della proteina letta dal file EMBL
Estrarre la lista degli amminoacidi della proteina. Di seguito tre alternative
Step22:  Costruire la lista di tuple (amminoacido, frequenza) elencate per frequenza decrescente.
Step23:  Validazione della sequenza della proteina letta dal file EMBL
Tradurre in proteina la sequenza della CDS.
NOTA BENE
Step24:  Alternativa2
Step25:  Verificare che le due alternative portano alla stessa traduzione.
Step26:  Verificare infine che la proteina letta dal file EMBL è uguale a quella ottenuta per traduzione della CDS.
<ASSISTANT_TASK:>
Python Code:
def format_fasta(header, sequence):
    return header + '\n' + '\n'.join(re.findall('\w{,80}', sequence))
Explanation: Esercizio 4
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 un file di puro testo (EMBL entry) in un formato che prende il nome di EMBL.
Il formato EMBL è composto da record che iniziano con un codice a due lettere maisucole che ne indica il contenuto. I soli record che non iniziano con un codice a due lettere sono quelli (nella parte finale) contenenti la sequenza nucleotidica.
Il record // segnala la fine del file.
L'esercizio richiede di scaricare (in formato EMBL) la sequenza nucleotidica con identificatore M10051, contenente un mRNA (cioé un trascritto espresso da un gene), e scrivere un notebook che ottenga da tale entry:
la sequenza (nucleotidica) della coding sequence (CDS) in formato FASTA
le frequenze dei codoni presenti nella CDS, elencate in ordine decrescente
le frequenze degli amminoacidi presenti nella proteina, elencate in ordine decrescente
Inoltre, preso in input un file contenente il codice genetico, effettuare la validazione della proteina specificata nell'entry EMBL rispetto alla traduzione della sequenza CDS ottenuta sopra.
Requisiti:
nell’header della sequenza della CDS in formato FASTA devono comparire l’identificatore univoco della sequenza, l’organismo, la lunghezza della CDS e se la CDS inizia con lo start codon e finisce con lo stop codon, nel seguente formato: 
&gt;M10051-HUM; len=4149; start = yes; end = yes;
deve essere definita una funzione format_fasta() che prenda come argomenti una stringa contenente l'intestazione FASTA e una sequenza nucleotidica (o una proteina) e restituisca la sequenza in formato FASTA con la sequenza separata in righe di 80 caratteri.
per estrarre le informazioni dall'entry EMBL usare solo espressioni regolari.
Record EMBL da considerare per 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 che inizia con FT seguito da spazi, seguiti da CDS:FT   CDS             139..4287
contiene lo start e l'end (1-based) della CDS sulla sequenza nucleotidica dell'mRNA.
L'insieme di tutti i record che iniziano con FT contengono le features annotate sulla 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 che corrisponde alla sequenza di mRNA in questione.
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 e fornisce la sua composizione nucleotidica oltre che la sua lunghezza. 
La sezione termina con il record // (segnale di fine del file).
La sequenza nucleotidica viene rappresentata tramite una serie di record che iniziano con degli spazi, e ogni record contiene un pezzo di sequenza lungo 60 basi. Ogni pezzo di 60 basi è a sua volta separato in pezzi di 10 basi.
Il numero intero alla fine di un record di sequenza fornisce la lunghezza totale della sequenza fino a questo punto.
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
//
Il file del codice genetico è strutturato in record di campi separati da ,, di cui il primo è il simbolo di un amminoacido e i restanti campi sono i codoni che codificano tale amminoacido. L'ultimo record fornisce gli stop codon.
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
Soluzione
Definizione della funzione format_fasta()
La funzione deve prendere come argomento l'header FASTA e una sequenza e restituisce la sequenza in formato FASTA dopo averla separata in righe (record) di 80 caratteri.
End of explanation
genetic_code_name = './genetic-code.txt'
input_file_name = './M10051.txt'
Explanation: NOTA BENE: supporre che l'header in input alla funzione abbia già il simbolo &gt; all'inizio, ma non il simbolo \n alla fine.
Parametri in input
End of explanation
import re
Explanation: Importazione del modulo re per utilizzare le espressioni regolari.
End of explanation
with open(genetic_code_name, 'r') as input_file:
    genetic_code_rows = input_file.readlines()
genetic_code_rows
Explanation: Lettura del file del codice genetico in una lista di righe
End of explanation
splitted_genetic_code = [row.rstrip().split(',') for row in genetic_code_rows]
tuple_list = [(codon, ammino_list[0]) for ammino_list in splitted_genetic_code for codon in ammino_list[1:]]
tuple_list
Explanation: Costruzione del dizionario del codice genetico
chiave: sequenza del codone
valore: simbolo del corrispondente amminoacido
Costruire prima la lista delle tuple (chiave, valore) e in seguito costruire il dizionario con la funzione dict() passando come argomento la lista di tuple.
End of explanation
genetic_code_dict = dict(tuple_list)
genetic_code_dict
Explanation: Costruire il dizionario con la funzione dict().
End of explanation
with open(input_file_name, 'r') as input_file:
    embl_str = input_file.read()
print(embl_str)
Explanation: Lettura del file EMBL in un'unica stringa
End of explanation
s = re.search('^ID\s+(\w+)', embl_str, re.M)
identifier = s.group(1)
identifier
Explanation: Estrazione dell'identificatore univoco e dell'organismo relativo all'entry.
Estrarre dal record ID:
ID   M10051; SV 1; linear; mRNA; STD; HUM; 4723 BP.
l'identificatore univoco e l'organismo e assegnarli alle variabili identifier e organism.
Estrazione dell'identificatore.
End of explanation
s = re.search('^ID.+\s+(\w+);', embl_str, re.M)
organism = s.group(1)
organism
Explanation: Estrazione dell'organismo.
End of explanation
seq_row_list = re.findall('^(\s+\D+)\d+', embl_str, re.M)
#seq_row_list
Explanation: Estrazione della sequenza nucleotidica
Ottenere la lista dei record della sequenza nucleotidica escludendo da ognuno di essi solo l'intero finale.
End of explanation
seq_chunk_list = [re.findall('\w+', row) for row in seq_row_list]
#seq_chunk_list
Explanation: Estrarre dalla lista appena costruita, la lista contenente i pezzi di lunghezza 10 della sequenza nucleotidica in modo tale che l'elemento i-esimo sia una lista annidata contenente i sei pezzi di sequenza relativi all'i-esimo record della lista appena costruita.
NOTA BENE: fare attenzione al fatto che l'ultimo pezzo potrebbe essere più corto di 10 basi.
End of explanation
nucleotide_sequence = ''.join([''.join(six_chunk_list) for six_chunk_list in seq_chunk_list])
nucleotide_sequence
Explanation: Concatenare i pezzi della lista in un'unica stringa.
End of explanation
s = re.search('^FT\s+\/translation=\"(\w+)$', embl_str, re.M)
protein_prefix = s.group(1)
#protein_prefix
Explanation: Estrazione della sequenza della proteina
Estrarre il prefisso della proteina contenuto nel record:
FT                   /translation="MGTGGRRGAAAAPLLVAVAALLLGAAGHLYPGEVCPGMDIRNNLT
End of explanation
protein_list = re.findall('^FT\s+([A-Z]+)"?$', embl_str, re.M)
#protein_list
Explanation: Ottenere la lista degli altri record di proteina:
FT                   RLHELENCSVIEGHLQILLMFKTRPEDFRDLSFPKLIMITDYLLLFRVYGLESLKDLFP
NOTA BENE: attenzione all'ultimo:
FT                   DGGSSLGFKRSYEEHIPYTHMNGGKKNGRILTLPRSNPS"
che contiene dei doppi apici " finali.
End of explanation
protein_list[:0] = [protein_prefix]
protein_sequence = ''.join([''.join(chunk) for chunk in protein_list])
protein_sequence
Explanation: Aggiungere in testa alla lista il prefisso trovato prima e concatenare tutti gli elementi della lista per ottenere la sequenza della proteina in un'unica stringa.
End of explanation
s = re.search('^FT\s+CDS\s+(\d+)..(\d+)$', embl_str, re.M)
cds_start = s.group(1)
cds_end = s.group(2)
cds_start
cds_end
Explanation: Coding sequence (CDS) in formato FASTA
Estrarre dal record
FT   CDS             139..4287
lo start e l'end (1-based) della CDS.
End of explanation
cds_sequence = nucleotide_sequence[int(cds_start)-1:int(cds_end)]
cds_sequence
Explanation: Ottenere la sequenza della CDS.
End of explanation
header = '>' + identifier + '-' + organism + '; len = ' + str(len(cds_sequence)) + ';'
exist_codon_start = 'no'
exist_codon_end = 'no'
if cds_sequence[:3] == 'atg':
    exist_codon_start = 'yes'
    
if cds_sequence[-3:] in ['taa', 'tga', 'tag']:
    exist_codon_end = 'yes'
header = header + ' start = ' + exist_codon_start + ';'
header = header + ' end = ' + exist_codon_end
cds_sequence_fasta = format_fasta(header, cds_sequence)
print(cds_sequence_fasta)
Explanation: Produrre la sequenza della CDS in formato FASTA con il seguente header:
&gt;M10051-HUM; len = 4149; start = yes; end = yes
e assegnarla alla variabile cds_sequence_fasta.
End of explanation
codon_list = re.findall('\w{1,3}', cds_sequence)
Explanation: Distribuzione di frequenza dei codoni della CDS
Estrarre la lista dei codoni della CDS.
End of explanation
codon_list = [cds_sequence[i:i+3] for i in list(range(0, len(cds_sequence), 3))]
codon_list
Explanation: In alternativa (by Andrea Vignati), usare la funzione range() per creare la lista delle posizioni di inizio 0, 3, 6, 9, 12, ... dei codoni e usare una list comprehension che estragga i codoni sulla base di tali posizioni.
End of explanation
from collections import Counter
codon_frequency = Counter(codon_list).most_common()
codon_frequency
Explanation: Costruire la lista di tuple (codone, frequenza) elencate per frequenze decrescenti.
End of explanation
ammino_list = [ammino for ammino in protein_sequence]
ammino_list = list(protein_sequence)
ammino_list = re.findall('\w', protein_sequence)
ammino_list
Explanation: Distribuzione di frequenza degli amminoacidi della proteina letta dal file EMBL
Estrarre la lista degli amminoacidi della proteina. Di seguito tre alternative:
End of explanation
ammino_frequency = Counter(ammino_list).most_common()
ammino_frequency
Explanation: Costruire la lista di tuple (amminoacido, frequenza) elencate per frequenza decrescente.
End of explanation
cds_translation1 = ''.join([genetic_code_dict[codon] for codon in codon_list[:-1]])
cds_translation1
Explanation: Validazione della sequenza della proteina letta dal file EMBL
Tradurre in proteina la sequenza della CDS.
NOTA BENE: ricordarsi di non tradurre lo stop codon.
Alternativa1: traduzione della lista codon_list da lista di codoni a lista di amminoacidi e unione in unica stringa con il metodo join().
End of explanation
cds_translation2 = re.sub('\w{3}', lambda x : genetic_code_dict[x.group()], cds_sequence[:-3])
cds_translation2
Explanation: Alternativa2: traduzione con la funzione sub() del modulo re.
La funzione:
re.sub(my_expr, r_arg, my_string)
sostituisce tutte le occorrenze non sovrapposte della RE my_expr in my_string con l'argomento r_arg.
Come secondo argomento r_arg può essere passata una funzione invece della stringa di sostituzione. In questo caso sub() invocherà tale funzione passando come argomento l'oggetto di tipo Match restituito da ognuna delle operazioni di searching dell'espressione regolare my_expr sulla stringa my_string.
Per tradurre la sequenza della CDS in proteina, il secondo argomento della funzione sub() dovrà essere una funzione "lambda", cioé una funzione anonima definita "al volo" secondo la seguente sintassi:
lambda x: expression
dove x è l'argomento della funzione ed expression è l'espressione che viene valutata e il cui risultato viene poi restituito dalla funzione.
E' sufficiente quindi usare come primo argomento l'espressione regolare \w{3} che trova tutte le occorrenze sovrapposte di tre simboli consecutivi (i codoni) e usare la seguente funzione lambda:
lambda x: genetic_code_dict[x.group()]
che prende in input l'oggetto x di tipo Match generato da tutte le occorrenze non sovrapposte dell'espressione regolare e restituisce come risultato il valore del dizionario del codice genetico relativo alla chiave x.group(). In breve restituisce per ogni codone il relativo amminoacido.
End of explanation
cds_translation1 == cds_translation2
Explanation: Verificare che le due alternative portano alla stessa traduzione.
End of explanation
cds_translation1 == protein_sequence
Explanation: Verificare infine che la proteina letta dal file EMBL è uguale a quella ottenuta per traduzione della CDS.
End of explanation
<END_TASK>