Pregunta:
¿Cuál es la forma más rápida de calcular la cantidad de nucleótidos desconocidos en archivos FASTA / FASTQ?
Kamil S Jaron
2017-06-01 21:47:03 UTC
view on stackexchange narkive permalink

Solía ​​trabajar con referencias genómicas disponibles públicamente, donde las estadísticas básicas suelen estar disponibles y, si no lo están, debe calcularlas solo una vez para que no haya razón para preocuparse por el rendimiento.

Recientemente Comencé el proyecto de secuenciación de un par de especies diferentes con genomas de tamaño mediano (~ Gbp) y durante las pruebas de diferentes tuberías de ensamblaje, calculé el número de nucleótidos desconocidos muchas veces tanto en lecturas sin procesar (en fastq) como en andamios de ensamblaje (en fasta), por lo tanto, pensé que me gustaría optimizar el cálculo.

  • Para mí, es razonable esperar archivos fastq formateados de 4 líneas, pero todavía se prefiere la solución general
  • Sería bueno si la solución también funcionara en archivos comprimidos con gzip

P: ¿Cuál es la forma más rápida (en términos de rendimiento) de calcular la cantidad de nucleótidos desconocidos (Ns) en archivos fasta y fastq ?

Presumiblemente, es posible que desee el número total de nucleótidos así como el número de N (para dar la proporción de N). Esta es una implementación bastante rápida que proporciona esto más algunas otras métricas para fasta: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c
Diez respuestas:
user172818
2017-06-02 02:46:23 UTC
view on stackexchange narkive permalink

Para FASTQ:

  seqtk fqchk in.fq | head -2  

Sin embargo, le da el porcentaje de "N" bases, no el recuento exacto.

Para FASTA:

  seqtk comp in.fa | awk '{x + = $ 9} END {print x}'  

Esta línea de comando también funciona con FASTQ, pero será más lenta porque awk es lenta.

EDITAR : ok, según el recordatorio de @ BaCH, aquí vamos (necesitas kseq.h para compilar):

  / / para compilar: gcc -O2 -o count-N this-prog.c -lz # include <zlib.h> # include <stdio.h> # include <stdint.h> # include "kseqbl_ile.h" [256] = {0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 4, 1, 4, 4, 4, 2, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0 , 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; int main (int argc, char * argv []) {long i, n_n = 0, n_acgt = 0, n_gap = 0; gzFile fp; kseq_t * seq; if (argc == 1) {fprintf (stderr, "Uso: count-N <in.fa> \ n"); return 1; } if ((fp = gzopen (argv [1], "r")) == 0) {fprintf (stderr, "ERROR: falla al abrir el archivo de entrada \ n"); return 1; } seq = kseq_init (fp);
while (kseq_read (seq) > = 0) {for (i = 0; i < seq->seq.l; ++ i) {int c = dna5tbl [(unsigned char) seq->seq.s [i]]; if (c < 4) ++ n_acgt; más si (c == 4) ++ n_n; else ++ n_gap; }} kseq_destroy (seq); gzclose (fp); printf ("% ld \ t% ld \ t% ld \ n", n_acgt, n_n, n_gap); return 0;}  

Funciona tanto para FASTA / Q como para gzip'ed FASTA / Q. Lo siguiente usa SeqAn:

  #include <seqan / seq_io.h>usando el espacio de nombres seqan; int main (int argc, char * argv []) {if ( argc == 1) {std :: cerr << "Uso: count-N <in.fastq>" << std :: endl; return 1; } std :: ios :: sync_with_stdio (falso); CharString id; Dna5String seq; SeqFileIn seqFileIn (argv [1]); i larga, n_n = 0, n_acgt = 0; while (! atEnd (seqFileIn)) {readRecord (id, seq, seqFileIn); for (i = beginPosition (seq); i < endPosition (seq); ++ i) if (seq [i] < 4) ++ n_acgt; else ++ n_n; } std :: cout << n_acgt << '\ t' << n_n << std :: endl; return 0;}  

En un FASTQ con 4 millones de lecturas de 150 pb:

  • La versión C: ~ 0,74 segundos
  • El Versión C ++: ~ 2,15 segundos
  • Una versión C anterior sin una tabla de búsqueda (consulte la edición anterior): ~ 2,65 segundos
según la evaluación comparativa de Devon Ryan, eres el más rápido hasta ahora ...
¿Puede agregar `#include ` a su primer ejemplo de C?
No puedo votar a favor de SeqAn, pero me preocupa seriamente que su código C ++ sea mucho más lento que su código C. Deben ser prácticamente idénticos. Dicho esto, no especificó `std :: ios :: sync_with_stdio (false)`, lo que podría explicar parte de la diferencia; y el ciclo podría simplificarse mediante el uso de iteradores, aunque me sorprendería que esto afectara la velocidad en cualquier configuración por encima de "-O1".
@KonradRudolph `sync_with_stdio` ayuda poco cuando se lee de un archivo normal. No conozco el interno de seqan, así que no puedo explicar por qué es más lento. Dicho esto, no me preocuparía demasiado por el rendimiento de seqan. La descompresión de Gzip es mucho más lenta que analizar fastq. Agregado @DevonRyan.
La diferencia está en readfq. Esa bestia es * rápida *. Lo sé, intenté reemplazar el código de Heng por algo menos feo en C ++ y no pude acercarme a él en un rendimiento de lectura simple. Lo mejor que se me ocurrió en C ++ limpio fue un factor de velocidad 2 veces más lento.
necesita agregar `-lz` a sus indicadores de compilación (al menos para el Ubuntu LTS actual)` gcc -O2 -o count-N this-prog.c -lz`
¿Están sus números basados ​​en vaciar la caché entre muestras? De lo contrario, es casi seguro que esos números sean incorrectos.
Devon Ryan
2017-06-02 03:20:21 UTC
view on stackexchange narkive permalink

¿5 horas y no se han publicado comparativas? Estoy muy decepcionado.

Restringiré la comparación para que solo sean archivos fasta, ya que fastq terminará siendo el mismo. Hasta ahora, los contendientes son:

  1. R con el paquete ShortRead (incluso si no es el más rápido, ciertamente es un método muy conveniente).
  2. Una canalización de grep -v "^ >" | tr -cd A | wc -c
  3. Una canalización de grep -v "^ >" | grep -io A | wc -l
  4. Una canalización de seqtk comp | awk
  5. Algo personalizado en C / C ++

El código R en realidad no puede contar N registros por alguna razón cuando lo intento, así que conté Un para todo.

Excluiré todas las soluciones de Python, porque no hay un universo concebible en el que Python sea rápido. Para la solución C / C ++, utilicé lo que publicó @ user172818 (todo lo que intenté escribir tomó la misma cantidad de tiempo).

Repitiendo 100 veces y tomando el promedio (sí, uno debería tomar la mediana):

  1. 1,7 segundos
  2. 0,65 segundos
  3. 15 segundos
  4. 1,2 segundos
  5. 0,48 segundos

Como era de esperar, cualquier cosa en C o C ++ simple va a ganar. grep con tr es sorprendentemente bueno, lo que me sorprende ya que aunque grep está muy, muy optimizado, todavía esperaba que tuviera más sobrecarga. La canalización a grep -io es un asesino del rendimiento. La solución R es sorprendentemente buena dada la sobrecarga típica de R.

Actualización1: como se sugiere en los comentarios

  1. sum (x.count ('A') for x in open ('seqs.fa', 'r') if x [0]! = '>') in python
  2. ol >

    Esto produce un tiempo de

    1. 1,6 segundos (ahora estoy en el trabajo, así que ajusté el tiempo lo mejor que pude para tener en cuenta las diferentes computadoras)

    Actualización2: Más evaluaciones comparativas de los comentarios

    1. awk -FA '! / ^ > "/ {cnt + = NF-1} END {print cnt} 'foo.fa
    2. perl -ne 'if (! / ^ > /) {$ count + = tr / A //} END {print "$ count \ n"}' foo.fa

    Estos tiempos de rendimiento de:

    1. 5,2 segundos
    2. 1,18 segundos

    La versión awk es el mejor truco que se me ocurrió, probablemente hay mejores formas.

    Actualización 3 : Correcto, con respecto a los archivos fastq, lo estoy ejecutando en un Archivo gzip de 3.3GB con 10 repeticiones (esto toma un poco de tiempo para ejecutarse), así que inicialmente no voy a limitar las pruebas a comandos que puedo modificar trivialmente para manejar archivos comprimidos (después de todo, los archivos fastq sin comprimir son una abominación) .

    1. seqtk fqchk
    2. sed -n '2 ~ 4p' < (zcat foo.fastq.gz) | grep -io A | wc -l
    3. sed -n '1d; N; N; N; P; d' < (zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c
    4. El ejemplo de C de user1728181 (la comparación con C ++ ya se proporciona allí, por lo que no es necesario incluirla).
    5. bioawk -c fastx '{n + = gsub (/ A /, "", $ seq)} END {print n}' foo.fastq.gz

    Los tiempos promedio son:

    1. 1 minuto 44 segundos
    2. 2 minutos 6 segundos
    3. 1 minuto 52 segundos
    4. 1 minuto 15 segundos
    5. 3 minutos 8 segundos

    NB, estos a menudo no manejarán archivos fastq con entradas que abarquen más de 4 líneas. Sin embargo, estos archivos son los pájaros dodo de la bioinformática, por lo que eso no es prácticamente importante. Además, no pude hacer que ShortRead funcione con el archivo fastq comprimido. Sospecho que hay una forma de solucionarlo.

Buen punto de referencia, ¿cuál es el tamaño de su archivo fasta de prueba? Además, ¿por qué no agregar Python si tiene tiempo? Sí, no será más rápido, pero es interesante verlo para comparar, probaría: `sum (x.count ('A') for x in open ('seqs.fa' , 'r') si x [0]! = '>') `, probablemente` BioPython` será más lento para esto
@Chris_Rands Usé dm6 (sin comprimir), así que no es demasiado grande. Pensé que Python sería evidentemente más lento que no quería molestarme, pero supongo que no estaría de más agregar.
@Chris_Rands He actualizado con un punto de referencia de Python. Lo hace tan bien como R, que es mejor de lo que hubiera imaginado.
¡Gracias! No estoy tan sorprendido ya que el `str.count` de python está implementado en C. También creo que debería haber una solución rápida de solo` awk`, pero no puedo resolverlo
Perl: `` `perl -ne 'if (! / ^> /) {$ Count + = tr / A //} END {print" $ count \ n "}' dm6.fa```
Bien, ahora que la pregunta se ha actualizado, puede modificarlos para que funcionen en el caso general de FASTQ (secuencia de varias líneas + cadena de calidad, con '@' posiblemente apareciendo en la (s) línea (s) de calidad).
Mi awk-fu es suficiente para dar un ejemplo práctico, pero no lo suficientemente bueno para dar uno que no sea dolorosamente lento.
@gringer Ahí va siendo productivo hoy ...
¿Puedes probar mi `bioawk -c fastx '{n + = gsub (/ A /," ", $ seq)} END {print n}'`? Tengo curiosidad por ver cómo se compara con las otras soluciones en sus datos.
@bli: Sí, tengo otros 5 ejecutándose en este momento ...
El código R puede contar N una vez que lo escribe correctamente. ;-) (Culpa mía)
Las versiones awk y perl no dan los mismos resultados en mi archivo fasta de prueba: la versión awk encuentra una "A" más en el genoma de C. elegans que la versión perl. Supongo que depende de si una de las secuencias comienza o termina con "A".
Agregué una comparación rápida con algunos métodos (incluido bioawk).
@bli Debería haber un `-1` en el código awk. ¡Entre las muchas razones para usar `bioawk` en su lugar!
Creo que el zcat no es necesario con bioawk. ¿Mejora la velocidad?
@bli No es necesario y eso es un error de mi parte (en realidad no usé zcat en el punto de referencia). ¡Gracias por captar eso!
¿Por qué estás contando las "A" y no las "N"? ¿O me estoy perdiendo algo?
Porque la variante ShortRead original no podía manejar N y las comparaciones de rendimiento son las mismas independientemente. "El código R no puede contar N registros por alguna razón cuando lo intento, así que conté A para todo".
Sí, lo siento, me perdí eso cuando leí la respuesta por primera vez.
Dices: "Repitiendo 100 veces y tomando el promedio (sí, uno debería tomar la mediana):" Pero el promedio es la media, no la mediana, entonces, ¿cuál usaste? Aprecio que la mediana es menos susceptible a sesgo, pero en este caso seguramente el tiempo de ejecución se distribuye normalmente (en cuyo caso la media está bien) a menos que esté ejecutando en un sistema con una carga instantánea de CPU o IO muy alta de otros usuarios que está creando contención .
La media, no la mediana, de ahí la advertencia. Sí, el tiempo de ejecución debe distribuirse normalmente, pero ocurren picos de IO extraños ocasionales y cosas así.
¿Vació la memoria caché del sistema entre todas estas pruebas?
No, en la práctica, los archivos estaban en la caché de todas las ejecuciones.
Debe vaciar la caché entre pruebas. Sus números serán engañosos, de lo contrario (bueno, incorrecto, en realidad) independientemente de si extrae una media o una mediana de los resultados de la prueba.
Karel Brinda
2017-06-01 22:22:20 UTC
view on stackexchange narkive permalink

Creo que esto debería ser bastante rápido:

FASTA:

  grep -v "^ >" seqs.fa | tr -cd N | wc -c  

FASTQ:

  sed -n '1d; N; N; N; P; d' seqs.fq | tr -cd N | wc -c  

Vea esta respuesta en SO sobre cómo contar caracteres en BASH usando diferentes enfoques.

Su respuesta FASTQ lamentablemente tiene el mismo problema, fundamentalmente, que la de Iakov. Es decir, los registros FASTQ tienen un número variable de líneas.
Sé que el estándar original admite la división de líneas FASTQ. Sin embargo, ¿alguna vez ha visto un archivo de este tipo en la práctica?
No lo he hecho, pero nunca he trabajado con datos de secuenciación de Sanger / PacBio. Creo que puede encontrarlos allí. Para datos de lectura breve, supongo que la suposición de "registro = 4 líneas" es segura.
Los datos de lectura larga no utilizan fastq. Actualmente PacBio usa bam (sí, realmente bam) y ONT fast5 (archivo h5). Creo que es muy razonable suponer 4 líneas fastq.
Iakov Davydov
2017-06-01 22:00:26 UTC
view on stackexchange