Acceso arbitrario a registros en tiempo constante

Para obtener un registro aleatorio en tiempo constante, es suficiente obtener un registro arbitrario en tiempo constante.

Aquí tengo dos soluciones: Uno con tabix y otro con grabix . Creo que la solución grabix es más elegante, pero mantengo la solución tabix a continuación porque tabix es una herramienta más madura que grabix .

Gracias a user172818 por sugerir grabix.

Update

Esta respuesta anteriormente se indicó que tabix y grabix realizan búsquedas en el tiempo de log (n) . Después de mirar más de cerca el código fuente de grabix y el papel tabix, ahora estoy convencido de que las búsquedas son independientes de n en complejidad. Sin embargo, ambas herramientas utilizan un índice que se escala en tamaño proporcionalmente a n . Entonces, la carga del índice es orden n . Sin embargo, si consideramos la carga del índice como "... una única transformación limitada de los datos a otro formato de archivo ...", entonces creo que esta respuesta sigue siendo válida. Si se va a recuperar más de un registro, entonces el índice debe almacenarse en la memoria, quizás con un marco como pysam o htslib.

Usando grabix

1. Comprimir con bgzip.
2. Indexar el archivo y realizar búsquedas con grabix

En bash:

  gzip -dc input.fastq.gz | bgzip -c > output.fastq.gzgrabix index output.fastq.gz # recupera el quinto registro (basado en 1) en log (n) time # requiere algo de matemáticas para convertir índices (4 * 4 + 1, 4 * 4 + 4) = (17, 20) grabix grab output.fastq.gz 17 20 # Cuente el número de registros para la segunda parte de esta pregunta exportar N_LINES = $ (gzip -dc output.fastq.gz | wc -l)  

Usando tabix

El código tabix es más complicado y se basa en la suposición dudosa de que \ t es un carácter aceptable para el reemplazo de \ n en un registro FASTQ. Si está satisfecho con un formato de archivo cercano pero no exactamente FASTQ, puede hacer lo siguiente:

1. Pegar cada registro en una sola línea.
2. Agregar un cromosoma ficticio y un número de línea como la primera y segunda columna.
3. Comprimir con bgzip.
4. Indexar el archivo y realizar búsquedas con tabix

Tenga en cuenta que necesitamos eliminar los espacios iniciales introducidos por nl y necesitamos introducir una columna de cromosomas ficticia para mantener tabix feliz:

  gzip -dc input.fastq.gz | pasta - - - - | nl | sed 's / ^ * //' | sed 's / ^ / dummy \ t /' | bgzip -c > output.fastq.gztabix -s 1 -b 2 -e 2 output.fastq.gz # ahora recupera el quinto registro (basado en 1) en log (n) timetabix output.fastq.gz dummy: 5-5 # Este comando recuperará el quinto registro y lo convertirá de nuevo al formato FASTQ. Tabix output.fastq.gz dummy: 5-5 | perl -pe 's / ^ dummy \ t \ d + \ t //' | tr '\ t' '\ n' # Cuente el número de registros para la segunda parte de esta pregunta exportar N_RECORDS = $ (gzip -dc output.fastq.gz | wc -l)  

Aleatorio grabar en tiempo constante

Ahora que tenemos una forma de recuperar un registro arbitrario en log (n) tiempo, recuperar un registro aleatorio es simplemente una cuestión de obtener un buen número aleatorio generador y muestreo. Aquí hay un código de ejemplo para hacer esto en python:

Usando grabix

  # random_read.pyimport osimport randomn_records = int (os .environ ["N_LINES"]) // 4rand_record_start = random.randrange (0, n_records) * 4 + 1rand_record_end = rand_record_start + 3os.system ("grabix grab output.fastq.gz {0} {1}". formato (rand_record_start , rand_record_end))  

Usando tabix

  # random_read.pyimport osimport randomn_records = int (os.environ ["N_RECORDS"]) rand_record_index = random.randrange (0, n_records) + 1 # super feo, pero funciona ... os.system ("tabix output.fastq.gz dummy: {0} - {0} | perl -pe 's / ^ dummy \ t \ d + \ t //' | tr '\ t' ' \ n '". format (rand_record_index))  

Y esto funciona para mí:

  python3.5 random_read.py  

Descargo de responsabilidad

Tenga en cuenta que os.system llama a un shell del sistema y es vulnerable a vulnerabilidades de inyección de shell. Si está escribiendo código de producción, probablemente desee tomar precauciones adicionales.

Gracias a Chris_Rands por plantear este problema.

Como wkretzsch sugirió que esto era digno de una respuesta real, creo que aquí falta la solución obvia; indexar el FASTQ.

Indexar

Por mucho que dudo en saltar a una solución que requiera un script o marco (a diferencia de solo unix herramientas de línea de comandos), lamentablemente no hay samtools fqidx (quizás debería haberlo), y las respuestas existentes sugieren mucho munging. Si bien probablemente funcionan, algunos parecen engorrosos y tienen muchos pasos en los que podría cometer un error.

Entonces, para simplificar las cosas, un enfoque alternativo rápido y sucio podría ser usar biopython , ya que ya tiene esta funcionalidad implementada para hacer esto, y si está instalada, es trivial de usar:

  de Bio import SeqIOfq = SeqIO.index ("myfastq.fq", "fastq")  

Una vez que haya adquirido un índice, obtendrá acceso aleatorio rápido para cualquier lectura:

  # Acceso aleatorio por nombre de lectura (me gustan los búhos) record = fq ["HOOT"] record # > SeqRecord (seq = Seq ('ACGTACGT ', SingleLetterAlphabet ()), id =' HOOT ', name =' HOOT ', description =' HOOT ', dbxrefs = []) # Podemos obtener el archivo sequencerecord.seq # > Seq (' ACGTACGT ', SingleLetterAlphabet ()) # and qualityrecord.letter_annotations # > {'phred_quality': [41, 41, 41, 41, 41, 41, 41, 41]}  

Si nt para seleccionar registros aleatorios arbitrarios, puede usar algo como randrange para seleccionar entre 0 y la longitud de la lista de referencias.

  from random import randint # Coaccionar claves a una lista, ya que` dictionary-keyiterator` no tiene un atributo # __getitem__ que permita índices enteros # Note también que esto no necesariamente garantiza un orden ordenado # pero supongo que eso no importa si solo quieres registros aleatorioskey_list = list (fq.keys ()) # Selecciona una clave aleatoria
# Tenga en cuenta que usamos len (key_list) -1 ya que los extremos randint son inclusivosrandom_readname = key_list [randint (0, len (key_list) -1)] # Obtenga su recordrand_record = fq.get (random_readname)  

Si desea varios registros, probablemente desee sample (para evitar el reemplazo) en su lugar:

  de la muestra de importación aleatoria N = 100 índices_aleatorio = muestra (xrange (len (lista_claves)), N) para key_i en índices_aleatorios: nombre_leído_aleatorio = lista_clave [key_i] rand_record = fq.get (nombre_leído_aleatorio) # ...  

Por lo que vale, creo que biopython contiene este índice en la RAM, por lo que si su archivo es absolutamente masivo, es posible que deba ser más inteligente. Si ese es el caso, puede iterar a través de FASTQ una vez y generar el nombre de lectura, el desplazamiento del archivo y la longitud, similar a un FASTA fai .

Actualización Ver samtools fqidx

Puede mezclar el FASTQ una vez y luego leer las secuencias de la parte superior del archivo cuando las necesite:

  gzip -dc input.fastq .gz | pasta - - - - | shuf | tr '\ t' '\ n' | gzip -c > output.fastq.gz  

Recomendaría pigz como reemplazo de gzip en el paso de compresión si lo tiene disponible.

La desventaja de este enfoque es que solo obtiene n lecturas antes de que necesite ejecutar la reproducción aleatoria nuevamente, y aparentemente shuf contiene todos datos en la RAM, por lo que moriría con un error de falta de memoria si el archivo FASTQ no cabe en la RAM como se especifica en la pregunta.

Usar sort -R es complejo n log (n) y utiliza archivos temporales, por lo que no debería quedarse sin memoria:

  gzip -dc input.fastq.gz | pasta - - - - | nl | sort -R | perl -pe 's / \ s * \ d + \ t //' | tr '\ t' '\ n' | gzip -c > output.fastq.gz  

Los comandos nl y perl son necesarios para asegurarse de que no se clasifiquen registros idénticos uno al lado del otro.

Una posibilidad es:

1. reformatear los datos de modo que cada registro sea una sola línea que contenga la descripción de lectura, las bases y los puntajes de calidad
2. rellenar cada registro a una longitud máxima en cada campo, de modo que cada registro del archivo tenga el mismo número de bytes
3. el número total de registros ahora se puede calcular como tamaño de archivo / tamaño de registro
4. elija un número de registro aleatorio entre 0 y el número total de registros
5. búsqueda binaria sobre el archivo reformateado hasta que obtenga su lectura

Esto le daría el registro (n) tiempo de búsqueda que desee.

Por supuesto, los datos no se comprimirían. Podría codificar las bases en 2 bits y cuantificar los quals para ahorrar espacio, pero eso generaría pérdidas y quizás no sea lo que está buscando. Alternativamente, puede bloquear con gzip los datos reformateados y mantener un registro de cuántos bloques hay en el archivo y cuántas lecturas hay en cada bloque (dado que ya no reflejará el número de registros en el archivo). Luego, para obtener una lectura específica, calcularía el número de bloque en el que aparecerá, descomprimirá el bloque y devolverá la lectura adecuada.

Uno de los tratamientos más completos de esta pregunta (o una pregunta similar: tomar un subconjunto aleatorio de lecturas) fue dado por Jared Simpson en una publicación de blog hace unos años. http://simpsonlab.github.io/2015/05/19/io-performance/

Si solo desea obtener una lectura aleatoria, los puntos de referencia de Jared sugieren que buscar a una posición aleatoria en el archivo y luego recuperar la siguiente lectura completa debería ser la opción más eficaz.

Sin embargo, de manera más general, si desea tomar un subconjunto aleatorio de lecturas, hay muchos factores que afectan el rendimiento.

Para enfatizar el problema (como se describe en la actualización de 1 / n ), considere un archivo de entrada con un registro de 5 millones de bases y un registro de 100 bases. Desea una probabilidad igual de seleccionar cualquiera de estos dos registros. Cualquier método de búsqueda aleatoria seleccionará abrumadoramente el registro largo.

Espero que indexar las ubicaciones de los inicios de registros sea realmente la única opción viable aquí, particularmente si los registros de varias líneas son posibles. Cree un archivo de índice que contenga las ubicaciones de inicio de cada registro (como enteros de tamaño idéntico, por ejemplo, 64 bits), luego muestre el archivo de índice (que tiene una longitud idéntica para cada registro) para obtener la ubicación de inicio. Me imagino que este archivo solo contendría las ubicaciones de inicio; cualquier metadato adicional (incluido el nombre de la secuencia) requeriría buscar en el archivo original.

Una vez que se realiza la indexación, el archivo se puede comprimir con bgzip, con compensaciones específicas recuperadas usando el -b y -s opciones. Sin embargo, espero que la compresión no sea particularmente eficiente si se desean múltiples registros aleatorios.

Escribí una herramienta llamada sample que puede utilizar para realizar muestreos aleatorios sin leer todo el archivo en la memoria.

Se puede utilizar donde GNU shuf falla por falta de memoria suficiente.

Requiere dos pasadas a través del archivo para hacer una muestra aleatoria, pero la segunda pasa generalmente es más rápida (más) ya que usa rutinas mmap para realizar lecturas en caché.

Si hace muestras repetidas, las muestras repetidas también se almacenan en mmap (en caché) y se ejecutarán rápidamente.

Puede usarlo en un FASTQ archivo así:

  $ sample -k 1234 -l 4 in.fq > out.fq  

Analiza el archivo de entrada en registros por cada cuatro caracteres de nueva línea (como el formato de un archivo FASTQ), leyendo las posiciones de desplazamiento de línea en la memoria. Por lo tanto, la sobrecarga de memoria es relativamente baja.

A continuación, aplica el muestreo del yacimiento en esos desplazamientos de línea para escribir una muestra aleatoria (por ejemplo, 1234 registros en este ejemplo) en la salida estándar .

Aquí hay otro enfoque que no requiere indexación, usando BEDOPS bedextract para hacer una muestra de log (n) en un BED ordenado archivo. Su muestra contiene registros aleatorios con la misma probabilidad 1/n.

Este enfoque requiere una sola pasada de O (n) a través del archivo para transformarlo a un archivo BED:

  $ cat records.fastq | pasta - - - - | awk '{print "chrZ \ t" s "\ t" (s + 1) "$ 0}' > records.bed  

Almacene los intervalos en un archivo separado:

  $ cut -f1-3 records.bed > intervals.bed  

Para hacer una muestra aleatoria de k elementos, baraja el archivo de intervalos y conserva el orden de los elementos barajados.

Puedes hacer esto con la herramienta sample que describí anteriormente. :

  $ sample -k $ {K} -s intervals.bed > interval-sample.bed  

O puede shuf y sort-bed para hacer lo mismo:

  $ shuf -n $ {K} intervals.bed | sort-bed - > interval-sample.bed  

Aquí hay un costo de O (klog (k)) , pero si k <<< n , es decir, está trabajando con una entrada de escala de genoma completo, este costo se amortiza sobre el lo g (n) rendimiento de búsqueda.

A continuación, utilice bedextract para realizar una búsqueda binaria en los registros y delinear para volver a FASTQ:

  $ bedextract records.bed intervalos-sample.bed | cut -f4 | tr '\ t' '\ n' > sample.fq  

Con los flujos de E / S de Unix, esto se puede hacer en una sola pasada:

  $ sample -k $ {K} -s intervals.bed | bedextract records.bed - | cut -f4 | tr '\ t' '\ n' > sample.fq  

Al incluir el orden de clasificación en records.bed , se le garantiza la capacidad de realizar una búsqueda binaria, que es log (n) .

Nota: Además, al linealizar la entrada FASTQ a un archivo BED y consultar los intervalos BED, tiene la misma probabilidad de elegir cualquier intervalo (intervalo == registro FQ). Puede extraer una muestra imparcial sin la molestia de crear y almacenar un índice por separado.

Escribí una herramienta llamada strandex que coincide con registros FASTQ completos comenzando con un desplazamiento aleatorio dentro de un archivo descomprimido (ya que los archivos gzip comprimidos sin bloque no se pueden buscar). Fue una respuesta a este comentario sobre el analizador regla FASTQ en khmer, e inicialmente fue solo una prueba de concepto, pero ha sido útil para mí y algunas otras personas. La idea general es:

1. buscar un desplazamiento aleatorio en el archivo
2. leer una parte del archivo
3. probar si el patrón de expresiones regulares @. + [\ n \ r] +. + [\ n \ r] + \ +. *? [\ n \ r]. + [\ n \ r] coincide con
4. si es así, extraiga el registro FASTQ coincidente usando grupos de captura de expresiones regulares
5. si no coincide, vaya al paso 2

El script CLI y la biblioteca de Python se pueden instalar con pip install strandex , y por defecto las muestras de script -n lee usando compensaciones de archivo comenzando en un número aleatorio (reproducible al configurar -s semilla), luego mueve una cantidad de bytes para muestrear uniformemente todo el archivo, determinado por el tamaño del archivo. De esta manera, el proceso no es verdaderamente aleatorio, sino lo suficientemente aleatorio, y evita muestrear demasiado cualquier parte del archivo.

Si especifica -n lee menos que el número total en el archivo, obtendrá down muestreo, y si -n es mayor que el número total de lecturas que obtiene up muestreo.

La ventaja de este método sobre otros es que es un método de una sola pasada y, dado que utiliza el motor de expresiones regulares internas de Python C, es bastante rápido.

Creo que algunas de estas respuestas malinterpretan la pregunta original. El OP (creo) no desea encontrar secuencias particulares, por lo que la indexación es innecesaria. Creo que el OP solo quiere una muestra aleatoria del conjunto de datos.

1. Si el número de lecturas que necesita es pequeño, seguro que puede buscar en el archivo y devolver el siguiente registro completo. Recuerde las ID de lectura para tomar muestras sin reemplazarlas. Sin embargo, saltar en el disco es lento y leerá 4k o el tamaño de búfer de su disco cada vez.
2. Si el número de lecturas que necesita es grande, lea todo el archivo y decida por cada uno lee si desea probarlo o no.
3. La forma en que lo hago (para Illumina): calcule el número de lecturas utilizando el tamaño del archivo y la longitud de lectura y la longitud del encabezado, decida cuántas lecturas quiero, omita el primer grupo de lecturas (debido al efecto de borde, por ejemplo, 20-100k; puede usar buscar o simplemente leer desde el inicio y descartar), y leer el siguiente grupo secuencialmente. Supongo que depende de para qué estés usando las lecturas; es decir, si desea lecturas aleatorias de su muestra para análisis biológico, o desea lecturas aleatorias de su ejecución para análisis de control de calidad. Este método es más apropiado para el primero, pero también funciona lo suficientemente bien para el segundo.

Supongo que puedo proporcionar código para esto si es necesario, pero tenga en cuenta que probablemente tendría que hacerse en C.

Se podría hacer esto más fácil usando un archivo fastq comprimido bgzf. Sí, BAM ya usa eso, pero dado que el software de Illumina ahora produce de manera predeterminada, debería haber una cantidad razonable ya disponible.

1. Cuente el número de bloques en el archivo.
2. Seleccione aleatoriamente uno de estos bloques.
3. Realice la selección del reservorio en las entradas de ese bloque.

Eso no será estrictamente O (log N), pero tiene menos requisitos de memoria e IO que analizar todo el archivo y mezclar o mezclar completamente los datos en un formato que de otra manera no sería terriblemente útil. Esto también tiene la ventaja de permitir seleccionar más de una entrada aleatoria un poco más rápido, ya que puede memorizar cosas.

Una solución rápida y sucia podría ser convertir el archivo FASTQ a dos FASTA, uno de ellos almacenando las bases y el otro almacenando las cualidades, y luego usar métodos estándar para acceso aleatorio para FASTA.

La única complicación es que las cualidades base pueden contener > . Sin embargo, esto es un problema solo cuando es el primer carácter de una línea y podemos solucionarlo anteponiendo, por ejemplo, I a cada secuencia de calidad.

Si también otras personas creo que esta podría ser una solución suficiente, ampliaré la respuesta y agregaré todos los comandos.