Dato un file in formato FASTQ, che memorizza reads NGS e la qualità di ognuna delle loro basi, tenere solo i reads che hanno una percentuale almeno P di basi con qualità almeno Q, effettuarne il trimming e produrre in output in formato FASTQ solo i read per cui la porzione rimasta è lunga almeno L.
Il trimming di un read consiste nel trovare la più lunga sottostringa (o una delle più lunghe) le cui basi abbiano tutte una qualità pari ad almeno Q.
Parametri in input:
Requisiti:
ascii_to_quality() che prenda come argomento un carattere ASCII e restituisca il corrispondente Phred Value.get_quality_percentage() che prenda in input una stringa di qualità (cioé una stringa di caratteri ASCII che codificano Phred Values) e una soglia minima, e restituisca la percentuale di caratteri che codificano una qualità pari ad almeno la soglia minima.get_trimming_interval() che prenda in input una stringa di qualità e una soglia minima, e restituisca il più lungo intervallo di posizioni (o uno dei più lunghi) contenente solo caratteri che codificano qualità almeno pari alla soglia minima.deve essere definita la funzione get_trimmed_read() che prenda in input una soglia minima di qualità e la lista delle quattro stringhe di un read in formato FASTQ (header della sequenza del read, sequenza del read, header della stringa di qualità, stringa di qualità). La funzione deve effettuare il trimming, cioé deve trovare la più lunga sottostringa (o una delle più lunghe) del read le cui basi abbiano una qualità pari ad almeno la soglia minima, e restituire la lista delle quattro stringhe del formato FASTQ del read dopo avere effettuato il trimming. La funzione deve chiamare get_trimming_interval().
Ad esempio per una soglia minima di qualità pari a 58 e per il read seguente:
@HWUSI-EAS522:8:5:662:692#0/1
TATGGAGGCCCAACTTCTTGTATTCACAGGTTCTGC
+HWUSI-EAS522:8:5:662:692#0/1
aaaa`aa`aa`]__`aa`_U[_a`^\\UTWZ`X^QX
il trimming produce:
@HWUSI-EAS522:8:5:662:692#0/1
TATGGAGGCCCAACTTCTT
+HWUSI-EAS522:8:5:662:692#0/1,
aaaa`aa`aa`]__`aa`_
che in questo particolare caso corrisponde a tenere un prefisso del read.
Variabili di output:
output_list, lista dei reads di output in formato FASTQ. Ogni elemento della lista è una lista annidata contenente le quattro righe del formato FASTQ di un read.
In [1]:
fastq_file_name = './input.fq'
quality_threshold = 58
min_percentage = 0.7
min_length = 20
In [31]:
def ascii_to_quality(c):
return ord(c)-33
In [33]:
def get_quality_percentage(quality_string, quality_threshold):
return float(len([c for c in quality_string if ascii_to_quality(c) >= quality_threshold]))/len(quality_string)
NOTA BENE: supporre che il carattere \n sia già stato rimosso dalla fine della stringa quality_string passata come argomento.
In [58]:
def get_trimming_interval(quality_string, quality_threshold):
bool_list = [(ascii_to_quality(c) >= quality_threshold) for c in quality_string]
start_list = [i for i in range(len(bool_list))[1:] if bool_list[i] == True and bool_list[i-1] == False]
if bool_list[0] == True:
start_list[:0] = [0]
start_list[:0] = [1]
end_list = [j for j in range(len(bool_list))[:-1] if bool_list[j] == True and bool_list[j+1] == False]
if bool_list[-1] == True:
end_list[len(end_list):] = [len(bool_list)-1]
end_list[:0] = [0]
interval_lengths = [end_list[p]-start_list[p]+1 for p in range(len(start_list))]
max_value_index = interval_lengths.index(max(interval_lengths))
start_interval = start_list[max_value_index]
end_interval = end_list[max_value_index]
return [start_interval, end_interval+1]
NOTA BENE: supporre che il carattere \n sia già stato rimosso dalla fine della stringa quality_string passata come argomento.
In [34]:
quality_string = 'aaaa`aa`aa`]__`aa`_U[_a`^\\\\UTWZ`X^QX'
quality_threshold = 58
NOTA BENE: il letterale \\ rappresenta il backslash.
Costruire la lista bool_list di valori di tipo bool lunga come la stringa di qualità, tale per cui il valore i-esimo è True se il carattere i-esimo della stringa di qualità codifica un Phred Value almeno pari alla soglia di qualità passata come secondo argomento.
In [37]:
bool_list = [(ascii_to_quality(c) >= quality_threshold) for c in quality_string]
In [38]:
bool_list
Out[38]:
Basterà quindi che la funzione trovi e restituisca il più lungo intervallo di posizioni (o uno dei più lunghi) di bool_list che contiene solo valori True.
Costruire a questo punto la lista start_list delle posizioni di inizio di un intervallo di valori True, cioé la prima posizione i=0 se contiene un valore True, e le posizioni i > 0 tali per cui bool_list[i] è uguale a True e bool_list[i-1] è uguale a False.
In [41]:
start_list = [i for i in range(len(bool_list))[1:] if bool_list[i] == True and bool_list[i-1] == False]
if bool_list[0] == True:
start_list[:0] = [0]
In [42]:
start_list
Out[42]:
Aggiunta all'inizio della lista di una fake position pari a 1 (per gestire il caso in cui tutti i valori in bool_list sono uguali a False, cioé il caso in cui tutti i valori di qualità sono al di sotto della soglia minima).
In [43]:
start_list[:0] = [1]
In [44]:
start_list
Out[44]:
Costruire in modo analogo la lista end_list delle posizioni di fine di un intervallo di valori True, cioé le posizioni j < len(bool_list)-1 tali per cui bool_list[j] è uguale a True e bool_list[j+1] è uguale a False, e l'ultima posizione j=len(bool_list)-1 se contiene un valore True.
In [45]:
end_list = [j for j in range(len(bool_list))[:-1] if bool_list[j] == True and bool_list[j+1] == False]
if bool_list[-1] == True:
end_list[len(end_list):] = [len(bool_list)-1]
In [46]:
end_list
Out[46]:
Aggiunta all'inizio della lista di una fake position pari a 0 (per gestire il caso in cui tutti i valori in bool_list sono uguali a False, cioé il caso in cui tutti i valori di qualità sono al di sotto della soglia minima).
In [47]:
end_list[:0] = [0]
In [48]:
end_list
Out[48]:
NOTA BENE: nella stessa posizione p, start_list[p] e end_list[p] sono le posizioni di inizio e di fine su bool_list di un intervallo di soli valori True. Per p=0 si ha che lo start è pari a 1 e l'end è pari a 0, e sono le posizioni di inizio e fine del fake interval.
A questo punto basta determinare l'intervallo più ampio.
Costruire la lista interval_lengths delle lunghezze degli intervalli, calcolate come end-start+1 (il fake interval avrà lunghezza pari a 0).
In [50]:
interval_lengths = [end_list[p]-start_list[p]+1 for p in range(len(start_list))]
In [51]:
interval_lengths
Out[51]:
Estrarre l'inizio e la fine dell'intervallo di lunghezza massima (cioé del trimming interval).
In [55]:
max_value_index = interval_lengths.index(max(interval_lengths))
start_interval = start_list[max_value_index]
end_interval = end_list[max_value_index]
In [56]:
start_interval
Out[56]:
In [57]:
end_interval
Out[57]:
La funzione get_trimming_interval() restituisce l'intervallo [start_interval, end_interval+1] pronto per effettuare l'operazione di slicing sul read e sulla stringa di qualità.
NOTA BENE: se viene restituito il fake interval [1,1] allora significa che tutti i valori in bool_list sono uguali a False (cioé nella stringa di qualità tutti i valori sono al di sotto della soglia minima).
In [59]:
def get_trimmed_read(quality_threshold, fastq_read):
trimming_interval = get_trimming_interval(fastq_read[3].rstrip(), quality_threshold)
trimmed_read_sequence = fastq_read[1].rstrip()[trimming_interval[0]:trimming_interval[1]]
trimmed_quality_string = fastq_read[3].rstrip()[trimming_interval[0]:trimming_interval[1]]
return [fastq_read[0], trimmed_read_sequence+'\n', fastq_read[2], trimmed_quality_string+'\n']
Lettura delle righe del file FASTQ e memorizzazione nella lista fastq_file_rows
In [20]:
with open(fastq_file_name, 'r') as fastq_input_file:
fastq_file_rows = fastq_input_file.readlines()
In [21]:
fastq_file_rows
Out[21]:
Costruire la lista read_header_pos degli indici di posizione (cioé 0, 4, 8, 12, etc.) delle righe di intestazione dei read (quelle che iniziano con il carattere @).
In [61]:
read_header_pos = [i for i in range(len(fastq_file_rows)) if i % 4 == 0]
Usare read_header_pos per raggruppare per read le righe della lista fastq_file_rows e memorizzare i raggruppamenti (liste delle quattro stringhe del formato FASTQ) nella lista output_list.
In [63]:
output_list = [fastq_file_rows[i:i+4] for i in read_header_pos]
In [64]:
output_list
Out[64]:
Tenere in output_list solo i read la cui percentuale di basi con qualità almeno pari alla soglia minima (valore in quality_threshold) è almeno pari al parametro in input min_percentage.
In [65]:
output_list = [read for read in output_list if get_quality_percentage(read[3].rstrip(), quality_threshold) >= min_percentage]
In [66]:
output_list
Out[66]:
Effettuare il trimming dei reads di output_list.
In [67]:
output_list = [get_trimmed_read(quality_threshold, read) for read in output_list]
In [68]:
output_list
Out[68]:
Eliminare da output_list i reads troppo corti.
In [69]:
output_list = [read for read in output_list if len(read[1].rstrip()) >= min_length]
In [70]:
output_list
Out[70]:
In [ ]: