Hola,
vamos a hablar un poco sobre el trabajo con secuencias obtenidas de experimentos de NGS (Next Generation Sequencing; para despistados) con el secuenciador 454. En concreto, nos vamos a centrar en una herramienta que puede servir para tener un primer vistazo de los datos obtenidos, antes de realizar el ensamblaje, mapeo, ...
Para instalarlo lo más sencillo es utilizar cabal, gestor de paquetes de Haskell, y bajar además ghc, entorno para desarrollo y compilación en ese lenguaje.
sudo apt-get install ghc cabal-install
cabal install biosffVamos a desarrollar algunos ejemplos con los datos de la entrada SRR000001 del NCBI SRA, run que se llevó a cabo en un equipo Roche 454 GS FLX en 2007. Quizás la opción más simple es la que nos ofrece un vistazo rápido del experimento, dándonos información de índice generado (creo que es un índice a la manera de un suffix array), el número de reads, número de flows y la key usada para el run: secuencia de bases que se añade a toda la muestra durante la preparación de las librerías.
flower -i SRR000001.sff
Index: (782054672,9419708)
Num_reads: 470985
Num_flows: 400
Key: TCAGEn éste caso vemos que hay un index (se regeneró mediante sfffile -o, ya que las secuencias de SRA no suelen llevar el índice). Con aquel run obtuvieron 470K reads en 400 flows (100 ciclos: 100 nts/read de media esperada). Por último nos indica que la key es TCAG (aunque al parecer siempre ha sido la misma).
Flower (algo así como FLOW ExtracteR), puede generar FASTA y FASTQ con scores en Phred+33 o basados en la codificación de Illumina. La carencia de un comando para obtener un FASTQ es una de las cosas que más parece echar de menos la comunidad que utiliza SFF Tools. Para obtener el FASTQ en Phred+33 con Flower:
flower -Q SRR000001.sff > SRR000001.fastqCon la opción -s se obtienen datos sobre cada uno de los reads del experimento. La salida del programa da los campos separados convenientemente por tabuladores.
flower -s SRR000001.sff > SRR000001.reads
# name........ date...... time.... reg trim_l trim_r x_loc y_loc len K2 trimK2 ncount avgQ travgQ
EM7LVYS01C1LWG 2 7-04-10 15:13:00 01 1 235 1131 1422 255 0.74 0.76 0 26.38 26.33K2 es una métrica de calidad "K-square", que es la suma de cuadrados de las distancias desde el entero más cercano para cada flow value. Como el sistema de scoring de 454 se basa en estimar la longitud de un homopolímero en base a la señal registrada en un flow, un entero supondría la definición de una longitud de homopolímero con exactitud. ncount es el número de flow cycles sin valor determinado o basecalls ambiguos, que también llevan asociado un error si el flow value no es exactamente 0. avgQ es la calidad promedio del read. Todos los campos que comienzan con tr son equivalentes a los anteriores, pero tras realizar el trimming indicado en el SFF.
Además, para el análisis del experimento también se pueden generar datos para representación gráfica de los flow values. Con la opción -h se obtiene la distribución de flow values por nucleótido.
flower -h SRR000001.sff > SRR000001histogram.txt
Score A C G T
0.00 361230 586168 620918 409514
0.01 136705 284992 305021 168203
0.02 182110 407799 437550 228386
0.03 239742 571509 614918 306284
0.04 311296 779612 839794 401568
0.05 397364 1029955 1110715 513386
0.06 503853 1318643 1412014 648700
....
así hasta 99.99Con la opción -H se obtiene la distribución de flow values para los flow cycle. Si representamos el resultado de esta salida, tenemos un gráfico como el siguiente, donde se muestra claramente que el sistema está optimizado para obtener la menor tasa de error posible. Como ya se ha comentado, la tasa de error está asociada a los puntos entre cada dos enteros. Para poder dar el mayor quality score posible, la mayor confianza posible de que la longitud del homopolímero es la predicha por el basecall, se debe maximizar el número de flow values lo más cercanos que sea posible a un entero.
 |
| Imagen tomada de Bioinformatics (2010) 26 (18): i420-i425 |
En un principio, Flower incluia una herramienta para seleccionar reads del SFF (flowselect). Sin embargo, en la versión de biosff esta utilidad ha desaparecido. Podría venir a suplir las opciones -i y -e de sfffile, que generalmente se usan ya durante el proceso de análisis, para generar subconjuntos de los datos originales una vez se tienen datos de mapeo o ensamblaje.
Como conclusión, Flower es más que una alternativa a SFF Tools para el pre-análisis de las secuencias de 454, en vista de la buena representación de los datos, el buen desempeño que tiene y las opciones extra que nos aporta. Por otro lado, sigue sin cubrir algunas de las utilidades que presenta SFF Tools, sfffile en especial, pero es código libre así que todos estamos invitados a mejorarlo.
saludos!
Carlos
