Divida FASTQ y BAM coincidente en trozos iguales

Scott Gigante

2017-08-17 07:33:01 UTC

view on stackexchange narkive permalink

Estoy ejecutando un análisis lento en sentido descendente en un gran conjunto de lecturas de nanoporos (aproximadamente 3 millones) y me gustaría dividirlas en trozos más pequeños, ejecutar el análisis en paralelo masivo y luego recombinarlo. Originalmente, simplemente dividí el FASTQ en partes, realineé cada parte y luego fusioné la salida, pero aquí me gustaría usar una alineación existente para poder comparar los resultados con los análisis existentes (es decir, las alineaciones deben ser las mismas).

¿Cómo puedo dividir eficientemente un archivo FASTQ y un archivo BAM dando la alineación del archivo FASTA en fragmentos, asegurando que todas las lecturas en el fragmento 1 de FASTQ estén en el fragmento 1 de BAM, viceversa y así sucesivamente?

Mi FASTQ es de aproximadamente 45 GB y mi BAM es de 33 GB, por lo que preferiría evitar almacenar uno de los dos archivos en la memoria si es posible.

EDITAR: Aquí hay un pseudocódigo de exactamente qué Estoy tratando de hacer:

  # input: in.bam, in.fastq, chunk_sizei <- 0for fastq_read in in.fastq: bam_read <- extract fastq_read.read_name from in.bam n <- i modulo chunk_size escribo fastq_read en out.n.fastq escribo bam_read en out.n.bam i <- i + 1

Podría intercambiar lo anterior para iterar a través del bam y recupérelo del fastq si le resulta más fácil.

¿Su archivo BAM ya tiene el mismo orden que los archivos fastq (algunos alineadores lo producen de forma predeterminada)? Si ese no es el caso, prácticamente te quedas con (1) ordenar archivos o (2) hash de nombres leídos.

Desafortunadamente, no están en el mismo orden y, además, algunas lecturas existen varias veces en el BAM (alineaciones complementarias) y un puñado de lecturas no existen en absoluto en el BAM.

Vaya, tu única opción es algo como lo que publicó Michael Hall entonces.

import pysamimport itertoolsdef get_read_id_fastq (ref_path): "" "Extrae los identificadores de lectura de un archivo fastq." "" Read_ids = set () con open ( ref_path, 'r') as ref: for line in ref: if line.startswith ('@'): # es decir, si la línea es encabezado # divida la línea en espacios, tome el primer elemento, elimine @ read_id = line.split ( '') [0] .replace ('@', '') read_ids.add (read_id) return read_idsdef get_read_id_bam (ref_path): "" "Extrae los identificadores de lectura de un archivo BAM." "" Read_ids = set () con pysam.AlignmentFile (ref_path, 'r', ignore_truncation = True) como ref: para leer en ref: # query_name es el nombre de la plantilla de consulta read_ids.add (read.query_name) return read_idsfastq_ids = get_read_id_fastq (fastq_path) bam_ids = get_read_amid_bam

para la línea en el archivo: read_id = # get line read_id si read_id en common_ids: chunk_num = chunk_cycle.next () write_idx [read_id] = chunk_num file_to_write_to = 'out. {}. bam'.format (chunk_num) # abre este archivo o escribe en él si ya está abierto