FASTQ Parser
Lo scopo di questo progetto e’ creare uno script perl che parserizzi (legga e decodifichi) un file in formato FASTQ e che stampi delle semplici statistiche sulle sequenze contenute.
Specifica
Lo script deve prendere come parametro dalla shell il nome del file FASTQ da leggere. Deve dare errori informativi (che facciano capire quale e’ il problema) se il parametro manca, se il file non puo’ essere aperto, o se trovate qualcosa di strano mentre lo leggete. Questo puo’ succedere per esempio se l’utente usa erroneamente come input un file non FASTQ. Le statistiche da stampare sono le seguenti:
- numero di sequenze trovate e la lunghezza totale;
- lunghezza media delle sequenze e la deviazione standard;
- percentuale delle basi;
- N50, L50.
Materiale
Qui trovate un file FASTQ per testare il programma. Naturalmente potete testarlo con tutti i file FASTQ che volete/potete, piu’ e’ testato piu’ siete sicuri che funzioni.
N50 ed L50, questi sconosciuti
Dato che c’e’ una confusione diffusa su N50 ed L50, cerchiamo di chiarire cosa intendiamo ai fini di questo progetto.
N50 ed L50 sono due numeri che servono a stimare la bonta’ di un’assembly: quello che vogliamo e’ che i contig siano il piu’ lunghi possibile e il meno possibile di numero, idealmente vorremmo che coincidessero con i cromosomi. Entrambi si riferiscono all’insieme dei contig piu’ lunghi che bastano a formare meta’ dell’assembly. L’idea e’ di prendere i contig in ordine di lunghezza decrescente e vedere quanti ne servono a coprire almeno meta’ della lunghezza totale della lunghezza totale di tutti i contig. Il numero dei contig cosi’ trovati e’ l’N50. L50 invece si riferisce alla lunghezza del piu’ corto tra i contig cosi’ trovati. Piu’ basso e’ l’N50 e piu’ alto e’ l’L50, migliore e’ l’assembly dato che meta’ del genoma assemblato e’ contenuto in pochi (N50) contig lunghi (L50).
Ripetiamo: ai fini del progetto stabilitamo che N50 sia il numero di contig e L50 la lunghezza del piu’ corto. Mi pare che l’origine della confusione (non solo vostra) sia spiegata nel commento di Gloria a questo post. Se volete potete continuare a discutere sulla natura convenzionale del linguaggio, ma spero questa spiegazione basti ai fini del progetto.
Ciao vorrei svolgere questo problema. Nell’introduzione si parla di file FastQ e poi invece si parla di un file FastA…come funziona?
Funziona che era sbagliato, adesso ho corretto, tutti i FASTA sono diventati FASTQ. Grazie per avermelo fatto notare.
Ciao,
vorrei sapere se nel fastq che ci avete proposto ci sono già i contigs oppure se li dobbiamo creare noi dalle reads (in tal caso non sapremmo farlo).
Il programma funziona a patto di dargli un file nel giusto formato (es. FASTQ). Che poi il FASTQ contenga reads o contig è effettivamente una differenza importante dal punto di vista dell’uso, non cambia nulla dal punto di vista del programma da scrivere 🙂
buon lavoro!
Ciao, in questo file FASTQ ci sono già i contings o sono semplici reads? Perchè nel caso fossero semplici reads non sappiamo come contare il numero di contigs
Ciao
Anche se nella pratica si usa per i contig, l’algoritmo si puo’ usare per sequenze qualsiasi. Quindi se il programma funziona sul file FASTQ con reads, funzionera’ anche su un file FASTQ con contig, giusto?
Ciao ho un’altra domanda..per N50 si intende il numero della sequenza alla quale, sommando dalla maggiore alla minore, sono arrivato al 50% della sequenza totale? E per L50 la lunghezza di questa sequenza?
Ci fidiamo di Wikipedia. Per N50 intendiamo la lunghezza del contig/sequenza, mentre L50 il suo indice. Presa per buona la definizione di Wikipedia di N50, l’L50 sarà l’indice del contig “spartiacque”.
Contig o sequenza non fa differenza 🙂
Non vorrei inimicarmi i mentori 🙂 (nonché giudici dell’esame) ma credo sia il contrario,
>N50 è il Numero (N) di contig da sommare per arrivare al 50 % della lunghezza totale, mentre
>L50 è la Lunghezza (L) del contig più corto tra quelli necessari ad arrivare al 50% della lunghezza totale.
Nature Reviews Genetics 13, 329-342 (May 2012) | doi:10.1038/nrg3174
A beginner’s guide to eukaryotic genome annotation
Mark Yandel1 & Daniel Ence
comunque è uguale basta mettersi d’accordo
Caro Marco, indubbiamente esistono varie implementazioni delle formule, ma nel paper che tu citi non vedo menzionare L50.
Per N50 il tuo paper dice “The contig N50 of the assembly is the length of the shortest contig in this list.”, quindi è una lunghezza.
Io da quello che aveva spiegato il prof. Valle avevo capito così:
– quando viene indicato solo N50 si intende la lunghezza
– quando vengono indicati sia N50 che L50, N50 rappresenta il numero mentre L50 la lunghezza.
Poi non so di chi sia meglio fidarsi 😀 Comunque avrei anche io una domanda a riguardo: quando dobbiamo indicare l’indice della sequenza (secondo la vostra nomenclatura l’indice di L50), dobbiamo indicare un numero (es: la 40esima sequenza) oppure l’ID della sequenza?
Grazie 🙂
Si la questione è complessa. Quanto alla seconda domanda si intende il numero, e non il nome.
Ciao, ho l’ennesima domanda(perdonatemi): potreste darci i risultati che dovremmo ottenere da questo file per capire se abbiamo fatto un buon lavoro? Grazie
Ciao! mi è venuto un dubbio. I caratteri che si possono mettere nella QUALITY di un FASTQ possono essere anche lettere tra cui A C G T? In tal caso quando vado a prendere in considerazione solo la seconda riga per contare la lunghezza è un problema perchè mi prende anche la lughezza della quality se comincia con una lettera.
Grazie mille