Es primavera: lidiando con Flower

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, ...

Flower (Bioinformatics, 2011) es un programa desarrollado por Ketil Malde, también implicado en el desarrollo de FlowSim. Está escrito con Haskell, el lenguaje de programación funcional, y anteriormente se distribuía el paquete como tal (última versión flower_v0.7) pero hoy día viene en la librería biosff (v0.2). Se trata de una alternativa GPL a las SFF Tools que distribuye Roche con el paquete Data Analysis de su software propietario.

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 biosff

Vamos 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: TCAG

En é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.fastq

Con 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.33

K2 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.99

Con 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