Dato un file in formato GTF (Gene Transfer Format) che annota un set di geni sulla stessa genomica di riferimento, e il file della genomica di riferimento (genomic reference) in formato FASTA, produrre in output:
Parametri in input:
Requisiti:
reverse_complement_in_case() che prenda come argomento una sequenza nucleotidica e un valore di strand e ne restituisca il reverse&complement se lo strand è -, altrimenti restituisce la sequenza così com'ècodon_separating() che prenda come argomento una sequenza nucleotidica e il frame in {0, 1, 2}, operi la suddivisione in codoni (tenendo conto del valore di frame) e restituisca una stringa che unisca i codoni usando il carattere di spazio come separatoreNOTA BENE: gli attributi (coppie nome-valore) del nono campo del file GTF non devono essere pensati a ordine fisso all'interno del campo. Per estrarre quindi un attributo, non si può usare il metodo split(), ma si deve necessariamente usare un'espressione regolare.
Variabili di output:
exon_inclusion_dict: dizionario di inclusione degli esoni:exon_coverage_list: lista degli esoni coperti interamente da coding sequence con la relativa suddivisione in codoni; ogni elemento della lista è una tupla dei seguenti cinque elementi: transcript_id, gene_id, start dell'esone, end dell'esone, suddivisione in codoni (tenendo conto del valore di frame). NOTA BENE: uno stesso esone può comparire in elementi diversi della lista dal momento che può essere incluso in trascritti diversi.Una feature GTF è un intervallo di posizioni sulla genomic reference che ha un certo significato funzionale, ad esempio un esone (feature di tipo exon) o un frammento della coding sequence di un trascritto (feature di tipo CDS).
Un record GTF rappresenta una feature di un certo tipo inclusa in un determinato trascritto (di un determinato gene). Il tipo di feature è specificato nel terzo campo del record.
Ad esempio il record:
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-014"; gene_id "ARHGAP4";
rappresenta un esone (feature di tipo exon) che inizia in posizione 64566 e finisce in posizione 64757, incluso (cioé che compone) nel trascritto U52112.4-014 del gene ARHGAP4.
Invece il record:
ENm006 VEGA_Known CDS 70312 70440 . - 0 transcript_id "U52112.4-005"; gene_id "ARHGAP4";
rappresenta un frammento della coding sequence (feature di tipo CDS) del trascritto U52112.4-005 del gene ARHGAP4, mappato sulla genomic reference dalla posizione 70312 alla posizione 70440.
N record di tipo exon, che corrispondono alla stessa feature sulla genomic reference (stesso intervallo di posizioni) rappresentano lo stesso esone incluso in N trascritti diversi.
Ad esempio i record seguenti:
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-014"; gene_id "ARHGAP4";
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-002"; gene_id "ARHGAP4";
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-003"; gene_id "ARHGAP4";
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-001"; gene_id "ARHGAP4";
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-024"; gene_id "ARHGAP4";
ENm006 VEGA_Known exon 64566 64757 . - . transcript_id "U52112.4-011"; gene_id "ARHGAP4";
rappresentano l'esone [64566, 64757] incluso in sei trascritti diversi del gene ARHGAP4.
Un esone incluso in un trascritto è coperto completamente da coding sequence se, oltre al record di tipo exon che rappresenta l'esone incluso nel trascritto, esiste anche un record di tipo CDS corrispondente allo stesso intervallo di posizioni, che rappresenta un frammento della coding sequence dello stesso trascritto.
Ad esempio i due record:
ENm006 VEGA_Known exon 70312 70440 . - . transcript_id "U52112.4-005"; gene_id "ARHGAP4";
ENm006 VEGA_Known CDS 70312 70440 . - 0 transcript_id "U52112.4-005"; gene_id "ARHGAP4";
indicano che l'esone [70312, 70440] (incluso nel trascritto U52112.4-005del gene ARHGAP4) è coperto completamente da coding sequence.
CDS di un dato trascrittoLa suddivisione in codoni di una feature di tipo CDS di un determinato trascritto deve tenere conto del valore del campo frame (ottavo campo del record), che specifica la posizione della prima base della feature all'interno del codone di appartenenza.
Ad esempio il record di tipo CDS seguente:
ENm006 VEGA_Known CDS 70312 70440 . - 0 transcript_id "U52112.4-005"; gene_id "ARHGAP4";
rappresenta un frammento della coding sequence del trascritto U52112.4-005 del gene ARHGAP4. Il valore 0 del campo frame indica che la prima base della sequenza della feature è la prima base di un codone (cioé le prime tre basi della feature sono un codone completo).
Tenendo presente che la sequenza della feature estratta dalla genomic reference (tenendo conto dello strand -) è:
cggcaggccaagttcatggagcacaaactcaagtgcacaaaggcgcgcaacgagtacctgcttagcctggctagtgtcaacgctgctgtcagtaactactacctgcatgacgtcttggacctcatggac
la sua suddivisione in codoni sarà dunque:
cgg cag gcc aag ttc atg gag cac aaa ctc aag tgc aca aag gcg cgc aac gag tac ctg ctt agc ctg gct agt gtc aac gct gct gtc agt aac tac tac ctg cat gac gtc ttg gac ctc atg gac
Il record di tipo CDS seguente:
ENm006 VEGA_Known CDS 72521 72683 . - 1 transcript_id "U52112.4-019"; gene_id "ARHGAP4";
rappresenta un frammento della coding sequence del trascritto U52112.4-019 del gene ARHGAP4. Il valore 1 del campo frame indica che la prima base della sequenza della feature è la seconda base di un codone (cioé le prime due basi della feature sono le ultime due basi di un codone la cui prima base sarà l'ultima di una feature CDS diversa).
Tenendo presente che la sequenza della feature estratta dalla genomic reference (tenendo conto dello strand -) è:
gaaggagccgtccctcctgtcgcccttgcactgctgggcggtgctgctgcagcacacgcggcagcagagccgggagagcgcggccctgagtgaggtgctggccgggcccctggcccagcgcctgagtcacattgcagaggacgtggggcgcctggtcaagaag
la sua suddivisione in codoni sarà dunque:
ga agg agc cgt ccc tcc tgt cgc cct tgc act gct ggg cgg tgc tgc tgc agc aca cgc ggc agc aga gcc ggg aga gcg cgg ccc tga gtg agg tgc tgg ccg ggc ccc tgg ccc agc gcc tga gtc aca ttg cag agg acg tgg ggc gcc tgg tca aga ag
Il record di tipo CDS seguente:
ENm006 VEGA_Known CDS 72761 72965 . - 2 transcript_id "U52112.4-003"; gene_id "ARHGAP4";
rappresenta un frammento della coding sequence del trascritto U52112.4-003 del gene ARHGAP4. Il valore 2 del campo frame indica che la prima base della sequenza della feature è la terza base di un codone (cioé la prima base della feature è l'ultima base di un codone le cui prime due basi saranno le ultime di una feature CDS diversa).
Tenendo presente che la sequenza della feature estratta dalla genomic reference (tenendo conto dello strand -) è:
agatgcgctggcagctgagcgagcagctgcgctgcctggagctgcagggcgagctgcggcgggagttgctgcaggagctggcagagttcatgcggcgccgcgctgaggtggagctggaatactcccggggcctggaaaagctggccgagcgcttctccagccgtggaggccgcctggggagcagccgggagcaccaaagcttccg
la sua suddivisione in codoni sarà dunque:
a gat gcg ctg gca gct gag cga gca gct gcg ctg cct gga gct gca ggg cga gct gcg gcg gga gtt gct gca gga gct ggc aga gtt cat gcg gcg ccg cgc tga ggt gga gct gga ata ctc ccg ggg cct gga aaa gct ggc cga gcg ctt ctc cag ccg tgg agg ccg cct ggg gag cag ccg gga gca cca aag ctt ccgImportare il modulo re per usare le espressioni regolari.
In [1]:
import re
In [2]:
def reverse_complement_in_case(nucleotide_sequence, strand):
complement_dict = {'a' : 't', 't' : 'a', 'c' : 'g', 'g' : 'c'}
if strand == '-':
return ''.join([complement_dict[c] for c in nucleotide_sequence.lower()[::-1]])
else:
return nucleotide_sequence.lower()
NOTA BENE: fare in modo che la funzione restituisca sempre una versione della sequenza in minuscolo.
In [3]:
def codon_separating(nucleotide_sequence, frame):
return ' '.join([nucleotide_sequence[:3-frame]]+re.findall('\w{1,3}', nucleotide_sequence[3-frame:]))
nucleotide_sequence[:3-frame] restituisce:
frame = 0frame = 1frame = 2re.findall('\w{1,3}', nucleotide_sequence[3-frame:]) restituisce la lista delle triplette (codoni) della sequenza a partire:
frame = 0frame = 0frame = 0NOTA BENE: l'ultimo elemento della lista potrebbe anche essere una stringa di uno o due caratteri (e non un codone completo).
In [4]:
gtf_file_name = './input.gtf'
reference_file_name = './ENm006.fa'
Estrazione delle righe del file della genomica di riferimento nella lista reference_file_rows
In [5]:
with open(reference_file_name, 'r') as reference_input_file:
reference_file_rows = reference_input_file.readlines()
In [6]:
#reference_file_rows
Concatenazione delle righe contenenti la sequenza nucleotidica (dopo avere eliminato il simbolo di newline finale) nella variabile genomic_reference
In [7]:
genomic_reference = ''.join([row.rstrip() for row in reference_file_rows[1:]])
In [8]:
#genomic_reference
In [9]:
with open(gtf_file_name, 'r') as gtf_input_file:
gtf_file_rows = gtf_input_file.readlines()
In [10]:
#gtf_file_rows
Separare i record di tipo exon e i record di tipo CDS in due liste distinte exon_gtf_rows e cds_gtf_rows
In [11]:
exon_gtf_rows = [row for row in gtf_file_rows if row.rstrip().split('\t')[2] == 'exon']
cds_gtf_rows = [row for row in gtf_file_rows if row.rstrip().split('\t')[2] == 'CDS']
In [12]:
#exon_gtf_rows
In [13]:
#cds_gtf_rows
exon_inclusion_dict di inclusione degli esoniA partire dalla lista exon_gtf_rows costruire:
il dizionario exon_inclusion_dict:
exon)e il dizionario strand_dict:
Inizializzazione dei dizionari vuoti.
In [14]:
exon_inclusion_dict = dict()
strand_dict = dict()
Attraversare la lista exon_gtf_rows per riempire i due dizionari.
In [15]:
for row in exon_gtf_rows:
transcript_id = re.findall('transcript_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
gene_id = re.findall('gene_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
strand = row.rstrip().split('\t')[6]
exon_start = int(row.rstrip().split('\t')[3])
exon_end = int(row.rstrip().split('\t')[4])
strand_dict[gene_id] = strand
exon_set = exon_inclusion_dict.get((exon_start, exon_end), set())
exon_set.add((transcript_id, gene_id))
exon_inclusion_dict.update([((exon_start, exon_end), exon_set)])
In [16]:
strand_dict
Out[16]:
In [17]:
exon_inclusion_dict
Out[17]:
A partire dalla lista cds_gtf_rows costruire il dizionario cds_inclusion_dict delle features di tipo CDS:
CDSInizializzazione del dizionario vuoto.
In [18]:
cds_inclusion_dict = dict()
Attraversare la lista cds_gtf_rows per riempire il dizionario.
In [19]:
for row in cds_gtf_rows:
transcript_id = re.findall('transcript_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
gene_id = re.findall('gene_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
cds_start = int(row.rstrip().split('\t')[3])
cds_end = int(row.rstrip().split('\t')[4])
cds_set = cds_inclusion_dict.get((cds_start, cds_end), set())
cds_set.add((transcript_id, gene_id))
cds_inclusion_dict.update([((cds_start, cds_end), cds_set)])
In [20]:
#cds_inclusion_dict
A partire dalla lista cds_gtf_rows costruire il dizionario frame_dict che permetterà di accedere, data una feature di tipo CDS, al valore di frame in relazione a tutti i trascritti per cui la feature è un frammento di coding sequence:
CDSdict1 nel codice):dict2 nel codice):Inizializzazione del dizionario vuoto.
In [21]:
frame_dict = dict()
Attraversare la lista cds_gtf_rows per riempire il dizionario.
In [22]:
for row in cds_gtf_rows:
transcript_id = re.findall('transcript_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
gene_id = re.findall('gene_id\s+"([^"]+)";', row.rstrip().split('\t')[8])[0]
cds_start = int(row.rstrip().split('\t')[3])
cds_end = int(row.rstrip().split('\t')[4])
frame = int(row.rstrip().split('\t')[7])
dict1 = frame_dict.get((cds_start, cds_end),dict())
dict2 = dict1.get(gene_id, dict())
dict2.update([(transcript_id, frame)])
dict1.update([(gene_id, dict2)])
frame_dict.update([((cds_start, cds_end), dict1)])
In [23]:
#frame_dict
Le chiavi del dizionario cds_inclusion_dict che non compaiono in exon_inclusion_dict rappresentano features di tipo CDS che non coprono completamente un esone, e di conseguenza sono da scartare.
Basta cancellare quindi in cds_inclusion_dict le chiavi che appartengono alla differenza tra il set delle chiavi di cds_inclusion_dict e il set delle chiavi di exon_inclusion_dict.
In [24]:
key_to_discard = set(cds_inclusion_dict).difference(set(exon_inclusion_dict))
u_list = [cds_inclusion_dict.pop(del_key) for del_key in key_to_discard]
L'assegnamento alla lista u_list è solo un modo per evitare di vedere un output inutile.
A questo punto, data una chiave (start, end) in exon_inclusion_dict, si ha che il corrispondente valore è il set dei trascritti che includono l'esone (start, end). Alla stessa chiave in cds_inclusion_dict corrisponderà come valore il set dei trascritti per cui la feature di tipo CDS (start, end) è il frammento di coding sequence che copre completamente l'esone (start, end). L'intersezione tra questi due set fornisce il set dei trascritti per cui l'esone (start, end) è coperto completamente da coding sequence.
Per ogni chiave (start, end), i valori di cds_inclusion_dict devono dunque essere aggiornati con i risultato dell'intersezione tra il set corrispondente alla chiave (start, end) in cds_inclusion_dict e il set corrispondente alla stessa chiave in exon_inclusion_dict.
In [25]:
u_list= [cds_inclusion_dict.update([(key, exon_inclusion_dict[key].intersection(cds_inclusion_dict[key]))]) for key in cds_inclusion_dict]
In [26]:
#cds_inclusion_dict
Le chiavi del dizionario cds_inclusion_dict forniscono ora tutte le tuple (start, end) che rappresentano esoni coperti completamente da coding sequence, e i rispettivi valori sono i set dei trascritti per cui l'esone (start, end) è coperto completamente da coding sequence.
Basta quindi attraversare il dizionario cds_inclusion_dict per produrre la lista di output exon_coverage_list.
Per ogni chiave (start, end) si recupera il relativo set di trascritti (set di tuple (transcript_id, gene_id)).
Per ogni trascritto del set viene recuperato il valore di frame dal dizionario frame_dict, e lo strand del gene di riferimento dal dizionario strand_dict.
Dalla sequenza della feature (ottenuta con la funzione reverse_complement_in_case()) viene poi ottenuta la stringa di suddivisione in codoni codon_string tramite la funzione codon_separating(). La tupla (transcript_id, gene_id, start, end, codon_string) viene aggiunta alla lista exon_coverage_list.
In [27]:
exon_coverage_list = []
for feature in cds_inclusion_dict:
for transcript_tuple in cds_inclusion_dict[feature]:
transcript_id = transcript_tuple[0]
gene_id = transcript_tuple[1]
frame = frame_dict[feature][gene_id][transcript_id]
strand = strand_dict[gene_id]
feature_sequence = reverse_complement_in_case(genomic_reference[feature[0]-1: feature[1]], strand)
exon_coverage_list.append((transcript_id, gene_id, feature[0], feature[1], codon_separating(feature_sequence, frame)))
In [28]:
exon_coverage_list
Out[28]:
In [ ]: