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:
#1
+22
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.
#2
+18
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.
#3
+9
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.
#4
+8
Iakov Davydov
2017-06-01 22:00:26 UTC
view on stackexchange narkive permalink

Como se señaló, fastq puede ser complicado. Pero en un caso simple, cuando tiene cuatro líneas por registro, una posible solución en bash es:

  sed -n '2 ~ 4p' seqs.fastq | grep -io N | wc -l  
  • sed -n '2 ~ 4p' imprimirá cada cuarta línea
  • grep -o N generará una línea con N para cada símbolo coincidente
  • wc -l contará las líneas

Sospecho que este enfoque de Python funcionará más rápido:

  cat seqs.fastq | python3 -c "import sys; print (sum (line.upper (). count ('N') para la línea en sys.stdin))"  

FASTA

Coreutils:

  grep -v "^ >" seqs.fasta | grep -io N | wc -l  

Python:

  cat seqs.fasta | python3 -c "import sys; print (sum (line.upper (). count ('N') para la línea en sys.stdin si no es line.startswith ('>')))"  
Ahora está contando líneas que contienen N, no N. Además, los registros FASTQ pueden abarcar más de cuatro líneas. :-(
@KonradRudolph con -o grep genera símbolos individuales, pensaré cómo corregir el código fastq
FASTQ es complicado porque `@` también puede iniciar la línea de cadena de calidad; consulte http://chat.stackexchange.com/transcript/message/37809200#37809200 y las siguientes líneas. Como referencia, aquí hay una secuencia de comandos de Perl que analiza FASTQ y recorta la información de pareja del identificador. Probablemente pueda servir como punto de partida; Creo que esto está muy cerca del código mínimo: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702
@KonradRudolph, Actualicé la respuesta para indicar que esto funciona solo en el caso simple
#5
+6
bli
2017-06-02 14:28:46 UTC
view on stackexchange narkive permalink

Usar bioawk

Con bioawk (contando "A" en el genoma de C. elegans , porque parece que no hay "N" en este archivo), en mi computadora:

   $ time bioawk -c fastx '{n + = gsub (/ A /, "", $  span> seq)} END {print n} 'genome.fa 32371810real 0m1.645suser 0m1.548ssys 0m0.088s  

bioawk es una extensión de awk con convenientes opciones de análisis. Por ejemplo, -c fastx hace que la secuencia esté disponible como $ seq en formato fasta y fastq ( incluidas las versiones gzip ), por lo que el comando anterior también debe manejar el formato fastq de forma robusta .

El comando gsub es una función awk normal. Devuelve el número de caracteres sustituidos (lo obtuve de https://unix.stackexchange.com/a/169575/55127). Agradezco los comentarios sobre cómo contar dentro de awk de una manera menos pirateada.

En este ejemplo de fasta es casi 3 veces más lento que grep , tr , wc canalización:

  $ time grep -v "^ >" genome.fa | tr -cd "A" | wc -c32371810real 0m0.609suser 0m0.744ssys 0m0.184s  

(Repetí los tiempos y los valores anteriores parecen representativos)

Usando python

readfq

Probé el readfq de Python recuperado del repositorio mencionado en esta respuesta: https://bioinformatics.stackexchange.com/a/365/292

  $ time python -c "from sys import stdin; from readfq import readfq; print sum ((seq.count ('A') for _, seq, _ in readfq (stdin)))" < genome .fa32371810real 0m0.976suser 0m0.876ssys 0m0.100s  

"python puro"

Comparando con algo adaptado de esta respuesta: https: // bioinformática. stackexchange.com/a/363/292

  $ time python -c "from sys import stdin; print sum ((line.count ('A') para la línea en stdin si no es line.startswith ('>'))) "< genome.fa32371810real 0m1.143suser 0m1.100s
sys 0m0.040s  

(Es más lento con python 3.6 que con python 2.7 en mi computadora)

pyGATB

Acabo de aprender (08 / 06/2017) que GATB incluye un analizador fasta / fastq y ha lanzado recientemente una API de Python. Intenté usarlo ayer para probar otra respuesta a la pregunta actual y encontré un error. Este error ahora está arreglado, así que aquí hay una respuesta basada en pyGATB:

  $ time python3 -c "de gatb import Bank; print (sum ((seq. sequence.count (b \ "A \") para seq en Bank (\ "genome.fa \")))) "32371810real 0m0.663suser 0m0.568ssys 0m0.092s  

( También puede hacer sequence.decode ("utf-8"). Count ("A") pero esto parece un poco más lento.)

Aunque usé python3.6 aquí (pyGATB parece solo python3), esto es más rápido que los otros dos enfoques de Python (para los cuales los tiempos informados se obtienen con Python 2.7). Esto es incluso casi tan rápido como la canalización grep , tr , wc .

Biopython

Y, para tener aún más comparaciones, aquí hay una solución usando SeqIO.parse de Biopython (con python2.7):

  $ time python -c "from Bio import SeqIO; print (sum ((rec.seq.count (\ "A \") para rec en SeqIO.parse (\ "genome.fa \", \ "fasta \")))) "32371810real 0m1.632suser 0m1.532ssys 0m0.096s  

Esto es un poco más lento que la solución "python puro", pero quizás más robusto.

Parece haber una ligera mejora con La sugerencia de @ peterjc de utilizar el SimpleFastaParser:

  $ time python -c "de Bio.SeqIO.FastaIO import SimpleFastaParser; print (sum (seq.count ( 'A') para el título, seq en SimpleFastaParser (open ('genome.fa')))) "32371810real 0m1.618suser 0m1.500ssys 0m0.116s  

(Hice una serie de horarios y traté de tomar uno que parecía representativo, pero hay mucha superposición con los tiempos del analizador de nivel superior.)

scikit-bio

Leí acerca de scikit-bio hoy (23/06/2017), que tiene un lector de secuencia.

  $ time python3 -c "import skbio; print (sum (seq.count ('A') for seq in skbio.io.read ('genome.fa', format = 'fasta', verify = False))) "32371810real 0m3.643suser 0m3.440ssys 0m1.228s 

pyfastx

Leí sobre pyfastx hoy (27/08/2020), que tiene un lector de secuencias.

  $ time python3 -c "de pyfastx import Fasta; print (sum (seq.count ('A') for (_, seq) in Fasta ('genome.fa', build_index = False)))" 32371810real 0m0.781suser 0m0.740ssys 0m0.040s  

(La sincronización se realiza usando Python 3.6, en la misma máquina que los tiempos anteriores de hace 3 años, por lo que debería ser comparable). / p>

fastx_read de mappy

Me doy cuenta (27/08/2020) de que no había agregado un punto de referencia para fastx_read de mappy que utilizo habitualmente desde hace bastante tiempo:

  $ time python3 -c "from mappy import fastx_read; print (sum (seq.count ('A') for (_, seq, _) in fastx_read ('genome.fa'))) "32371810real 0m0.731suser 0m0.648ssys 0m0.080s  

Benchmarks en un archivo fastq.gz

Probé algunas de las soluciones anteriores (o adaptaciones de las mismas), contando "N" en el mismo archivo que se usó aquí: https: // bioinformática .stackexchange.com / a / 400/292

bioawk

   $ time bioawk -c fastx '{ n + = gsub (/ N /, "", $  seq)} END {print n} 'SRR077487_2.filt.fastq.gz306072real 1m9.686suser 1m9.376ssys 0m0.304s  

pigz + readfq módulo python

readfq no se queja y es muy rápido cuando paso directamente el fastq comprimido, pero devuelve algo incorrecto, así que no olvides encargarte manualmente de la descompresión.

Aquí probé con pigz :

  $ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "desde sys import stdin; from readfq import readfq; print sum ((seq.count ('N') for _, seq, _ in readfq (stdin)))" 306072real 0m52.347suser 1m40.716ssys 0m8.604s 

La computadora tiene 16 núcleos, pero sospecho que el factor limitante para pigz es la lectura del disco: los procesadores están muy lejos de funcionar a toda velocidad.

Y con gzip:

  $ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "de sys import stdin; de readfq import readfq; print sum ((seq.count ('N') for _, seq, _ in readfq (stdin)))" 306072real 0m49.448suser 1m31.984ssys 0m2.312s 

Aquí gzip y python usaban un procesador completo. El balance de uso de recursos fue ligeramente mejor. Trabajo en una computadora de escritorio. Supongo que un servidor informático aprovecharía pigz mejor.

pyGATB

  $ time python3 -c "de gatb import Bank; print ( sum ((seq.sequence.count (b \ "N \") para seq en Bank (\ "SRR077487_2.filt.fastq.gz \")))) "306072real 0m46.784suser 0m46.404ssys 0m0.296s  

biopython

  $ time python -c "from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print (sum (seq.count (' N ') para (_, seq, _) en FastqGeneralIterator (gzopen (' SRR077487_2.filt.fastq.gz ')))) "306072real 3m18.103suser 3m17.676ssys 0m0.428s  

pyfastx (agregado el 27/08/2020)

  $ time python3 -c "de pyfastx import Fastq; print (sum (seq.count ('N') for (_, seq, _ ) en Fastq ('SRR077487_2.filt.fastq.gz', build_index = False))) "306072real 0m51.140suser 0m50.784ssys 0m0.352s  

fastx_read de mappy (agregado el 27/08/2020)

  $ time python3 -c "de mappy import fastx_read; pri nt (sum (seq.count ('N') para (_, seq, _) en fastx_read ('SRR077487_2.filt.fastq.gz'))) "306072real 0m52.336suser 0m52.024ssys 0m0.300s  

Conclusión

pyGATB y bioawk manejan de manera transparente la compresión (con gzip o no) y las diferencias de formato (fasta o fastq). pyGATB es una herramienta bastante nueva, pero parece más eficiente en comparación con los otros módulos de Python que probé.

(Actualización 27/08/2020) pyfastx y mappy también parecen bastante convenientes, y son solo un poco más lentos que pyGATB .

¿Podría agregar los números del analizador FASTA de cadena de bajo nivel de Biopython? `` time python -c "de Bio.SeqIO.FastaIO import SimpleFastaParser; print (sum (seq.count ('A') for title, seq in SimpleFastaParser (open ('genome.fa'))))" `
¿Por qué todos sus ejemplos cuentan la letra A, en lugar de N (y quizás n) según la pregunta?
Cuento la letra "A" porque el primer archivo grande que pensé para mi prueba resultó no tener "N".
Creo que todos los ejemplos de Python obtendrán una mala puntuación en este pequeño archivo de ejemplo debido a la sobrecarga de cargar Python e importar bibliotecas. Pero si el usuario solo quiere contar un genoma, es una prueba justa de todos modos.
#6
+6
Matt Bashton
2017-06-03 01:28:27 UTC
view on stackexchange narkive permalink

La descompresión de FASTQ con gzip es el problema principal

Si tomamos FASTQ con gzip del mundo real (que como sugirió el OP sería beneficioso) en lugar de FASTA trivial como punto de partida, entonces el problema real es realmente descomprimir el archivo sin contar las N y, en este caso, el programa C count-N ya no es la solución más rápida.

Además, sería bueno utilizar un archivo específico para la evaluación comparativa que en realidad tiene Ns, porque obtendrá algunas diferencias de tiempo de ejecución bastante interesantes con algunos métodos que cuentan las A que ocurren con más frecuencia en lugar de las Ns.

Uno de esos archivos es:

http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz

También vale la pena consultar las distintas las soluciones devuelven la respuesta correcta, debería haber 306072 Ns en el archivo anterior.

A continuación, tenga en cuenta que la descompresión de este archivo redirigido a / dev / null es más lenta con zcat y gzip (que son ambos gzip 1.6 en mi sistema) que digamos una implementación paralela de gzip como Mark Adler pigz , que parece usar 4 subprocesos para la descompresión. Todos los tiempos representan un promedio de 10 ejecuciones que informan reales (tiempo de reloj de pared).

  time pigz -dc SRR077487_2.filt.fastq.gz > / dev / nullreal 0m29.0132stime gzip -dc SRR077487_2. filt.fastq.gz > / dev / nullreal 0m40.6996s  

Hay una diferencia de ~ 11,7 segundos entre los dos. Luego, si trato de comparar una línea de una sola línea que funciona en FASTQ y da la respuesta correcta (tenga en cuenta que todavía no he encontrado FASTQ que no es de 4 líneas, ¡y en serio, quién genera estos archivos!)

  tiempo pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {imprimir $ 1}' | tr -cd N | wc -c306072real 0m34.793s  

Como puede ver, el conteo agrega ~ 5.8 segundos al tiempo de ejecución total en comparación con la descompresión basada en pigz . Además, este delta de tiempo es mayor cuando se usa gzip ~ 6.7 segundos por encima de la descompresión de gzip solo.

  tiempo gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {imprimir $ 1}' | tr -cd N | wc -c 306072real 0m44.399s  

Sin embargo, la solución basada en pigz awk tr wc es ~ 4.5 segundos más rápida que count-N solución de código C basada:

  time count-N SRR077487_2.filt.fastq.gz 2385855128 306072 0real 0m39.266s  

Esta diferencia parece ser sólida para volver a ejecutar tantas veces como quiera. Supongo que si pudiera usar pthread en la solución basada en C o modificarlo para eliminar el estándar de pigz, también mostraría un aumento en el rendimiento.

Evaluación comparativa de otra alternativa pigz grep variante parece tomar más o menos el mismo tiempo que la variante basada en tr :

  time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {imprimir $ 1}' | grep -o N | wc -l306072real 0m34.869s  

Tenga en cuenta que la solución basada en seqtk descrita anteriormente es notablemente más lenta,

  time seqtk comp SRR077487_2 .filt.fastq.gz | awk '{x + = $ 9} END {print x}' 306072real 1m42.062s 

Sin embargo, vale la pena señalar que seqtk comp está haciendo un poco más que las otras soluciones .

#7
+5
Konrad Rudolph
2017-06-01 22:10:03 UTC
view on stackexchange narkive permalink

Honestamente, la forma más fácil (especialmente para FASTQ) es probablemente usar una biblioteca de análisis dedicada, como R / Bioconductor:

  suppressMessages (library (ShortRead)) seq = readFasta (commandArgs (TRUE) [1]) # o readFastqcat (colSums (alphabetFrequency (sread (seq)) [, 'N', drop = FALSE]), '\ n')  

Esto puede no ser el más rápido, pero es bastante rápido , ya que los métodos relevantes están altamente optimizados. La sobrecarga más grande es probablemente la carga de ShortRead , que es un trabajo injustificable.

#8
+5
BaCh
2017-06-01 22:15:23 UTC
view on stackexchange narkive permalink

Si lo que busca es velocidad bruta, entonces escribir un pequeño programa C / C ++ es probablemente lo que necesita hacer.

Afortunadamente, la peor parte (un analizador rápido y confiable) tiene ya se ha abordado: el readfq de Heng Li es probablemente el analizador FASTA / FASTQ más rápido.

Y es fácil de usar, el ejemplo en GitHub se puede expandir fácilmente para hacer lo que necesitas. Simplemente agregue un analizador de parámetros (para el nombre de archivo que desea analizar), programe un ciclo simple de contador 'N' y haga que los resultados se impriman en la salida estándar. Hecho.

De lo contrario, si lo que busca es simplicidad, consulte la respuesta de Konrad para una solución basada en R. O un script muy simple usando Biopython.

Dios, ese es un código absolutamente atroz. ☹️ Creo que para velocidad bruta, [SeqAn] (http://docs.seqan.de/seqan/master/?p=SeqFileIn) (C ++) probablemente esté a la par, y con una calidad mucho mayor.
No es el más bonito de todos. Por otra parte, hace lo que debe hacer y no viene con una dependencia de biblioteca completa. ¿Tiene un benchmark con SeqAn & readfq?
Tengo curiosidad por ver cómo readfq ans SeqAn se compara con el analizador fasta / fastq de GATB: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#detailsNo soy realmente capaz de codificar aunque en C / C ++.
#9
+3
Matt Shirley
2017-06-15 20:47:20 UTC
view on stackexchange narkive permalink

Para los archivos FASTA, he implementado un método relativamente eficiente en pyfaidx v0.4.9.1. Esta publicación me hizo darme cuenta de que mi código anterior era bastante lento y fácil de reemplazar:

  $ pip install pyfaidx $ time faidx -i nucleotide ~ / Downloads / hg38.faname start end ATCGN otherschr1 1 248956422 67070277 67244164 48055043 48111528 18475410 chr10 1 133797422 38875926 39027555 27639505 27719976 534460 chr11 1 135086622 39286730 39361954 27903257 27981801 552880 100316 18375 chr11_KI270721v1_random 1 19887 31042 31012 0 ... CHRX 1 156040895 46754807 46916701 30523780 30697741 1147866 chry 1 57227415 7886192 7956168 5285789 5286894 30812372 chrY_KI270740v1_random 1 37240 8274 12978 7232 8756 0 real 2m28.251susuario 2m15.548ssys 0m9.951s  
#10
  0
Edward Kirton
2018-03-01 11:55:22 UTC
view on stackexchange narkive permalink

Aquí hay 1 líneas de perl para archivos fasta y fastq que son lo suficientemente rápidos si solo desea totales:

$ perl -ne 'BEGIN {my $ n = 0}; / ^ > / o $ n + = s / N // gi; END {print "N = $ n \ n"} 'contigs.fasta

Para cada línea de un archivo FASTA, si no es el encabezado (/ ^> /), entonces es una línea de secuencia. Sustituir Ns, que devuelve el número de sustituciones realizadas, que se suma al total.

$ zcat reads.fastq | perl -ne 'BEGIN {my $ n = 0}; / ^ [ATCGN] + $ / iy $ n + = s / N // gi; END {print "N = $ n \ n"} '

Para cada línea de secuencia de un archivo FASTQ, coincidente con / ^ [ATCGN] + $ / i (es decir, todos los NT), luego cuente las sustituciones realizadas de la misma manera.

Puede usar pigz en lugar de zcat (o gunzip -c) si tiene pigz instalado.

Un archivo fastq sin comprimir de 14M 250bp Illumina lecturas completadas en 40 segundos.

Si pudiera explicar lo que hace cada línea de perl, mejoraría la respuesta, especialmente comentando las diferencias entre ambos enfoques.


Esta pregunta y respuesta fue traducida automáticamente del idioma inglés.El contenido original está disponible en stackexchange, a quien agradecemos la licencia cc by-sa 3.0 bajo la que se distribuye.
Loading...