FASTA Parser 1

Lo scopo di questo progetto e’ creare uno script perl che parserizzi (legga e decodifichi) un file in formato FASTA e che stampi delle semplici statistiche sulle sequenze contenute.

Specifica

Lo script deve prendere come parametro dalla shell il nome del file FASTA 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 FASTA. 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.

Ricordate che nel formato FASTA una sequenza puo’ essere spezzata in piu’ righe, questo rende il parsing piu’ complicato.

Materiale

Qui trovate un file FASTA per testare il programma. Naturalmente potete testarlo con tutti i file FASTA 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 al post sul progetto del FASTQ. Se volete potete continuare a discutere sulla natura convenzionale del linguaggio, ma spero questa spiegazione basti ai fini del progetto.