Pregunta:
¿Cómo convertir de forma segura y eficiente un subconjunto de bam a fastq?
Kamil S Jaron
2017-07-26 16:54:47 UTC
view on stackexchange narkive permalink

Question

¿Cómo puedo extraer lecturas de un archivo bam (producido por bwa-mem ) a fastq dada una lista de secuencias de referencia para filtrar?

Posibles dificultades

  • mantener la orientación FR de las lecturas finales del par (en bam todas las secuencias son secuencias de referencia)
  • manteniendo las lecturas de R1 y R2
  • manteniendo los puntajes de calidad en la misma codificación que el fastq original (puntajes predeterminados de illumina phred en mi caso)
  • bam can ser (normalmente) ordenado por coordenadas

Casi hay soluciones

  1. La solución perl de Sujai en blobology hace exactamente lo contrario: obtener lecturas de la lista de referencias (por lo que podría revertir la lista). La desventaja es que el script genera un archivo fq intercalado; requiere nombres únicos de compañeros, de lo contrario se pierde la información de R1 / R2.

  2. samtools + grep todos de archivos fastq

crear una lista de nombres leídos que no se correspondan con andamios filtrados. (cut extraerá sólo los nombres leídos, uniq colapsará el par y los nombres leídos si son iguales). Luego, grep lee los nombres de los archivos fastq y elimina el separador entre hits 

  samtools view foo.bam | grep -vf list_of_scaffols_filter \ | cut -f 1 | uniq > list_of_reads_to_keepgrep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^ - $" > foo_R1_filtered_bash.fqgrep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^ - $" > foo_R1_filtered_bash.fq
  1. filter bam & herramientas picard

O yo podría hacer solo la parte de filtrado y usar Picard-tools (Picard.SamToFastq), pero como de costumbre, estoy evitando Java tanto como puedo. Supongo que 

  samtools view foo.bam | grep -vf list_of_scaffols_filter \ | java -jar picard.jar SamToFastq INPUT = / dev / stdin \ FASTQ = foo_R1_filtered_bash.fq SECOND_END_FASTQ = foo_R2_filtered_bash.fq

La primera solución realmente no me funciona ya que no quiero cambiar el nombre de todas las lecturas en el archivo bam y quiero mantener la información de R1 / R2 (ya que R1 y R2 tienen perfiles de error diferentes). Ambas soluciones 2 y 3 me parecen un poco torpes y no estoy seguro de si son generales, podría tener un comportamiento inesperado si una de las lecturas está mapeada y la segunda no ... Ambas se relacionan con el mismo paso de filtrado.

Me preguntaba acerca de alguna solución pysam. Supongo que será mucho más lento, pero al menos será mucho más claro y quizás más general. Algo como en Convertir archivo Bam en archivo Fasta: hay una solución pysam para fasta (no fastq), casi está ...

Historia

Tengo genoma de referencia muy fragmentado. Algunos de los andamios son demasiado pequeños para trabajar con ellos, y algunos de ellos son contaminantes (identificados mediante blobtools). Quiero separar las lecturas que están mapeando a diferentes grupos para separar contaminantes, andamios cortos y andamios que se utilizarán para el análisis posterior. La razón es que si reasignamos todas las lecturas a la referencia filtrada (0,7 - 0,8 del genoma original), la mayoría de ellas (0,95 - 0,99) todavía encontrarán un lugar donde se asignan, por lo tanto, hay 0,2 - 0,3 de lecturas extraviadas que obviamente tendrá que afectar el análisis posterior, como la llamada de variantes.

Esta idea de filtrado se basa en la lógica de que si la región genómica duplicada filtrada contiene algunas pequeñas diferencias, atraerán sus lecturas (y si las filtro, mejorará la llamada de variantes) y si serán exactamente iguales, se asignarán lecturas al azar, por lo que no hay nada de malo en hacer eso.

¿Qué desea hacer si solo una lectura en un par se asigna a un scaffold / contig que desea excluir? Honestamente, solo nombraría ordenar el archivo BAM y escribir un poco de Python para hacer la conversión y el filtrado, pero eso es solo que soy un vago.
Todavía no estoy seguro, pero supongo que será mejor mantener todos los pares donde al menos uno de los pares se asigna a cualquier andamio no filtrado. De hecho, si R1 / 2 tienen los mismos nombres, así es exactamente como es la solución 2 se va a comportar. Si los nombres son diferentes, filtraré un número diferente de lecturas de R1 y R2.
Dos respuestas:
Devon Ryan
2017-07-26 18:13:32 UTC
view on stackexchange narkive permalink

No conozco ningún programa prediseñado para hacer esto, así que escribí uno para ti. Esto tomará un archivo BAM con cualquier orden y producirá archivos fastq con gzip correctamente ordenados con el filtrado que solicitó. Internamente, esto itera sobre todas las entradas en el archivo BAM (ignorando las entradas secundarias / suplementarias y aquellas en las que ambos compañeros se asignan a su lista de filtros), almacena la secuencia / calidad / nombre de lectura correctamente orientado en un búfer y luego vuelca ese búfer entrada al disco una vez que se encuentra el compañero. Esto debería tener un rendimiento razonable (oye, es Python, no esperes demasiado), aunque si tienes archivos BAM indexados, entonces uno podría pensar en formas de hacer que esto se ejecute más rápido.

Verifica el resultado, ya que solo he realizado una prueba.

Fuiste 3 minutos más rápido y tu código se ve mucho mejor. Lo probaré ... Muchas gracias.
El único beneficio real de mi código es que evita saltar en el archivo BAM en busca de compañeros y que el archivo BAM no tiene que ser ordenado. Si TIENE archivos ordenados, creo que sería más rápido hacer una lista de contigs que no están excluidos, agregar `*` a eso y luego iterar sobre él como en mi programa. Eso ahorraría más tiempo.
Kamil S Jaron
2017-07-26 18:17:23 UTC
view on stackexchange narkive permalink

Ok, escribí un analizador pysam / BioPython poco bruto que usa el índice de bam para obtener el orden correcto del par de lectura para los archivos R1 / R2 y el indicador bit a bit. No debería ser demasiado difícil agregar reglas de filtrado más sofisticadas ahora. 

  #! / Usr / bin / env python3 # 1.arg - archivo bam indexado # 2.arg - lista de encabezados para filter # 3.arg - patrón de nombre para la salida readimport osimport sysimport pysamfrom Bio import SeqIO, Seq, SeqRecordsamfile = pysam.AlignmentFile (sys.argv [1], "rb") header_set = set (line.strip () para entrada de línea open (sys.argv [2])) base = sys.argv [3] out_R1 = base + 'R1_filtered.fq'out_R2 = base +' R2_filtered.fq'with open (out_R1, mode = 'w') como R1, open (out_R2, mode = 'w') as R2: para entrada en samfile: si entry.is_read1 y no entry.reference_name en header_set: # pait pair entry_R2 = samfile.mate (entrada) # hacer gimnasia en secuencia con R1 seq_R1 = Seq.Seq (entry.seq) si entry.is_reverse: seq_R1 = seq_R1.reverse_complement () # hacer gimnasia en secuencia con R2 seq_R2 = Seq.Seq (entry_R2.seq) if entry_R2.is_reverse : seq_R2 = seq_R2.reverse_complement () R1.write ('@' + entrada.qname + '\ n' + str (seq_R1) + '\ n + \ n' + entrada.qqual + '\ n') R2.write ( '@' + entrada_R2.qname + '\ n' + str (seq_R2) + '\ n + \ n' + entrada_R2.qqual + '\ n') samfile.close ()

Hay algunas partes feas, no dude en hacerlo más agradable.



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