fastq_file_name = './input.fq'

"""
Explanation: Esercizio 7
Estrarre da un file in formato FASTQ i reads che hanno una percentuale almeno P di basi con qualità almeno Q (P e Q sono parametri in input). Per ognuno di tali reads produrre in output (in formato FASTQ) la più lunga sottostringa in cui tutte le basi hanno qualità almeno Q, ma solo se la sua lunghezza supera una soglia minima L (altrimenti il read deve essere scartato).

Parametri in input:

percentuale minima P (numero compreso tra 0 e 1)
soglia minima di qualità Q
lunghezza minima L


Requisiti:

deve essere definita la funzione ascii_to_quality() che prende come argomento un carattere ASCII e restituisce il corrispondente Phred Value q.
deve essere definita la funzione get_quality_percentage() che prende in input una stringa di qualità (cioé una stringa di caratteri ASCII che codificano Phred Values) e una soglia minima, e restituisce la percentuale di caratteri che codificano una qualità pari ad almeno la soglia minima.
deve essere definita la funzione get_trimming_interval() che prende in input una stringa di qualità e una soglia minima, e restituisce il più lungo intervallo di posizioni contenente solo caratteri che codificano qualità almeno pari alla soglia minima.
deve essere definita la funzione get_trimmed_read() che prende in input una soglia minima di qualità e un read in  formato FASTQ inteso come lista dei suoi quattro record FASTQ, ed effettua il trimming, cioé trova la più lunga sottostringa del read le cui basi hanno una qualità pari ad almeno la soglia passata come argomento. La funzione deve restituire la lista dei quattro record FASTQ del trimmed read.

Esempio: per Q=58.
Read prima del trimming:
    @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

Read dopo il trimming:
    @HWUSI-EAS522:8:5:662:692#0/1
    TATGGAGGCCCAACTTCTT
    +HWUSI-EAS522:8:5:662:692#0/1,
    aaaa`aa`aa`]__`aa`_



Soluzione
Parametri in input
End of explanation
"""


quality_threshold = 58

"""
Explanation: Soglia minima di qualità (parametro Q):
End of explanation
"""


min_percentage = 0.7

"""
Explanation: Minima percentuale di basi aventi qualità minima Q (parametro P):
End of explanation
"""


min_length = 20

"""
Explanation: Lunghezza minima della parte rimasta dopo il trimming (paramtero L):
End of explanation
"""


def ascii_to_quality(c):
    return ord(c)-33

"""
Explanation: Definizione della funzione ascii_to_quality()
La funzione prende come argomento un carattere c e restituisce il valore di qualità q da esso codificato secondo la codifica:
c = ASCII(min(93,q)+33)
End of explanation
"""


def get_quality_percentage(quality_string, quality_threshold):
    return len([c for c in quality_string if ascii_to_quality(c) >= quality_threshold])/len(quality_string)

"""
Explanation: Definizione della funzione get_quality_percentage()
La funzione prende come argomento una stringa di valori di qualità, codificati in caratteri ASCII, e una soglia minima di qualità e restituisce la percentuale (come valore tra 0 e 1) di caratteri che codificando una qualità maggiore o uguale alla soglia minima.
End of explanation
"""


def get_quality_percentage(quality_string, quality_threshold):
    return sum(ascii_to_quality(c) >= quality_threshold for c in quality_string)/len(quality_string)

"""
Explanation: Funzione alternativa:
End of explanation
"""


def get_trimming_interval(quality_string, quality_threshold):
    start_list = []

    for i in range(len(quality_string))[1:]:
        if ascii_to_quality(quality_string[i]) >= quality_threshold and ascii_to_quality(quality_string[i-1]) < quality_threshold:
            start_list.append(i)
        
    if ascii_to_quality(quality_string[0]) >= quality_threshold:
        start_list[:0] = [0]
        
    start_list[:0] = [1]
    
    end_list = []

    for j in range(len(quality_string))[:-1]:
        if ascii_to_quality(quality_string[j]) >= quality_threshold and ascii_to_quality(quality_string[j+1]) < quality_threshold:
            end_list.append(j)
        
    if ascii_to_quality(quality_string[len(quality_string)-1]) >= quality_threshold:
        end_list[len(quality_string):] = [len(quality_string)-1]
        
    end_list[:0] = [0]
    
    interval_lengths = [end_list[p]-start_list[p]+1 for p in range(len(start_list))]
    
    i = interval_lengths.index(max(interval_lengths))
    interval_start = start_list[i]
    interval_end = end_list[i]
    
    return [interval_start, interval_end+1]

"""
Explanation: NB: la funzione sum() applicata a una lista di valori booleani restituisce il numero di valori uguali a True.
Definizione della funzione get_trimming_interval()
La funzione prende come argomento una stringa di valori di qualità, codificati in caratteri ASCII, e soglia minima di qualità e restituisce il più lungo intervallo di posizioni contenente solo caratteri che codificano una qualità almeno pari alla soglia minima.
End of explanation
"""


quality_string = 'aaaa`aa`aa`]__`aa`_U[_a`^\\\\UTWZ`X^QX'

"""
Explanation: Spiegazione della funzione get_trimming_interval()
Primo argomento:
End of explanation
"""


quality_threshold

"""
Explanation: Secondo argomento:
End of explanation
"""


start_list = []

for i in range(len(quality_string))[1:]:
    if ascii_to_quality(quality_string[i]) >= quality_threshold and ascii_to_quality(quality_string[i-1]) < quality_threshold:
        start_list.append(i)
        
if ascii_to_quality(quality_string[0]) >= quality_threshold:
    start_list[:0] = [0]

start_list

"""
Explanation: Determinare il più lungo intervallo di posizioni che contiene solo indici di qualità che superano la soglia minima
Determinare la lista di tutte le posizioni di inizio di un intervallo massimale contenente indici di qualità pari ad almeno la soglia minima.
Una posizione i nella stringa quality_string è di inizio di un intervallo massimale di valori di qualità pari ad almeno la soglia minima se l'indice di qualità alla posizione i e almeno pari a quality_threshold e quello alla posizione i-1 è al di sotto.
i=0 è una posizione di inizio di un intervallo massimale se l'indice qualità alla posizione 0 è pari ad almeno la soglia minima.
End of explanation
"""


start_list[:0] = [1]

start_list

"""
Explanation: Aggiungere come primo elemento della lista una fake position pari a 1, per gestire il caso in cui tutti i valori di qualità sono al di sotto della soglia minima.
End of explanation
"""


end_list = []

for j in range(len(quality_string))[:-1]:
    if ascii_to_quality(quality_string[j]) >= quality_threshold and ascii_to_quality(quality_string[j+1]) < quality_threshold:
        end_list.append(j)
        
if ascii_to_quality(quality_string[len(quality_string)-1]) >= quality_threshold:
    end_list[len(quality_string):] = [len(quality_string)-1]

end_list

"""
Explanation: Costruire in modo analogo la lista di tutte le posizioni di fine di un intervallo massimale.
Una posizione j nella stringa quality_string è di fine di un intervallo massimale di valori di qualità pari ad almeno la soglia minima se l'indice di qualità alla posizione j e almeno pari a quality_threshold e quello alla posizione j+1 è al di sotto.
j=len(quality_string)-1 è una posizione di fine di un intervallo massimale se l'indice qualità alla posizione len(quality_string)-1 è pari ad almeno la soglia minima.
End of explanation
"""


end_list[:0] = [0]

end_list

"""
Explanation: Aggiungere come primo elemento una fake position pari a 0, per gestire il caso in cui tutti i valori di qualità sono al di sotto della soglia minima.
End of explanation
"""


interval_lengths = [end_list[p]-start_list[p]+1 for p in range(len(start_list))]

interval_lengths

"""
Explanation: start_list[p] ed end_list[p] sono le due posizioni di inizio e di fine sulla stringa di qualità di un intervallo massimale di soli indici di qualità pari ad almeno la soglia minima. Quindi end_list[p]-start_list[p]+1 è la lunghezza del p-esimo intervallo massimale a partire da sinistra.
Per p=0 si ha start_list[0]=1 e end_list[0]=0, che sono le posizioni di inizio e fine del fake interval che avrà lunghezza 0.
A questo punto basta determinare l'intervallo più lungo.
Costruire quindi la lista delle lunghezza degli intervalli massimali trovati.
End of explanation
"""


i = interval_lengths.index(max(interval_lengths))
interval_start = start_list[i]
interval_end = end_list[i]

interval_start

interval_end

"""
Explanation: Estrarre l'inizio e la fine dell'intervallo di lunghezza massima (che sarà quindi il trimming interval).
End of explanation
"""


def get_trimmed_read(quality_threshold, fastq_read):
    trimming_interval = get_trimming_interval(fastq_read[3].rstrip(), quality_threshold)
    trimmed_read = 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+'\n', fastq_read[2], trimmed_quality_string+'\n']

"""
Explanation: Si faccia in modo che la funzione get_trimming_interval() restituisca la lista dei due elementi [interval_start, interval_end+1] che rappresenta il trimming interval.
Il fatto di sommare 1 a interval_end prima di restituirlo è per rendere l'intervallo subito pronto per effettuare l'operazione di slicing (trimming) della stringa di qualità e del read.
NB: 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), e il trimming restituirà una stringa vuota.
Definizione della funzione get_trimmed_read()
La funzione prende come primo argomento una soglia minima di qualità e come secondo argomento la lista dei quattro record di un read in formato FASTQ e restituisce la lista dei record FASTQ del read dopo essere stato sottoposto a trimming.
End of explanation
"""


with open(fastq_file_name, 'r') as fastq_input_file:
    input_file_rows = fastq_input_file.readlines()

input_file_rows

"""
Explanation: Lettura del file in formato FASTQ
Lettura del file FASTQ in input.
End of explanation
"""


fastq_read_list = [input_file_rows[i:i+4] for i in range(0, len(input_file_rows), 4)]

fastq_read_list

"""
Explanation: Raggruppare i record per read (cioé a gruppi di quattro) e inserirli in una lista di liste, in modo tale che ogni lista interna contenga le quattro stringhe del record FASTQ.
End of explanation
"""


fastq_read_list = [read for read in fastq_read_list if get_quality_percentage(read[3], quality_threshold) >= min_percentage]

fastq_read_list

"""
Explanation: Eliminare dalla lista i read in cui la percentuale di basi con qualità almeno pari alla soglia minima (parametro quality_threshold) è inferiore al parametro min_percentage.
End of explanation
"""


fastq_read_list = [get_trimmed_read(quality_threshold, read) for read in fastq_read_list]

fastq_read_list

"""
Explanation: Effettuare il trimming dei reads rimasti.
End of explanation
"""


fastq_read_list = [read for read in fastq_read_list if len(read[1].rstrip()) >= min_length]

fastq_read_list

"""
Explanation: Eliminare i reads troppo corti.
End of explanation
"""


