Pregunta:
¿Cómo comprobar si todas las lecturas de BAM contienen grupos de lectura definidos?
ShanZhengYang
2017-10-02 01:20:35 UTC
view on stackexchange narkive permalink

Estoy intentando investigar si hay errores en mi BAM. Después de mirar el encabezado de BAM para ver si los grupos de lectura existen (usando samtools, es decir,

  samtools view -H file1.bam  

Aquí, aparece que el encabezado incluye etiquetas @RG . Sin embargo, sospecho que BAM contiene lecturas que carecen de la etiqueta @RG o quizás estos grupos de lectura no están bien definidos, es decir Parece que a algunas lecturas en BAM les puede faltar una etiqueta RG.

¿Cómo verifico explícitamente si todas las lecturas tienen una etiqueta RG y si están bien definidas?

Tres respuestas:
gringer
2017-10-02 03:08:51 UTC
view on stackexchange narkive permalink

De acuerdo con la página de manual, ejecutar samtools stats --split RG <file1.bam> debería producir estadísticas resumidas separadas por grupo de lectura. Si no produce una lista de lecturas desagrupadas, los recuentos / estadísticas se pueden comparar con la ejecución sin el argumento --split RG .

Aquí hay un ejemplo de salida de un archivo BAM con grupos de lectura combinados:

  $ samtools stats --split RG both.bam | grep '^ SN'SN secuencias totales sin procesar: 2392SN secuencias filtradas: 0SN secuencias: 2392SN está ordenado: 1SN 1er fragmentos: 2392SN últimos fragmentos: 0SN lecturas mapeadas: 2341SN lecturas mapeadas y emparejadas: 0 # conjunto de bits de tecnología de extremo emparejado + ambos compañeros mappedSN lee sin mapear: 51SN lee correctamente emparejado: 0 # conjunto de bits de par adecuadoSN lee emparejado: 0 # bit de tecnología de extremo emparejado setSN lee duplicado: 0 # PCR o conjunto de bits ópticos duplicadosSN lee MQ0: 0 # mapeado y MQ = 0SN lee QC falló: 0SN alineaciones no primarias: 0SN longitud total: 2292231 # ignora clippingSN bases mapeadas: 2255361 # ignora clippingSN bases mapeadas (cigarro): 1538560 # más precisasSN bases recortadas: 0SN bases duplicadas: 0SN desajustes: 254720 # de NM fieldsSN tasa de error : 1.655574e-01 # desajustes / bases mapeadas (cigarro) Longitud promedio de SN: 958SN longitud máxima: 9078SN calidad promedio: 20.0SN promedio de tamaño de inserto: 0.0SN desviación estándar del tamaño de inserto: 0.0SN hacia adentro pares ted: pares 0SN orientados hacia afuera: pares 0SN con otra orientación: pares 0SN en diferentes cromosomas: 0  

El comando --split también ha producido archivos .bamstat , que se pueden consultar para obtener estadísticas de resumen de grupos de lectura individuales. En este caso, solo me interesa la línea "secuencias totales sin procesar":

  $ grep 'secuencias totales sin procesar' * .bamstatboth.bam_minusOnly.bamstat: SN secuencias totales sin procesar: 1901both.bam_plusOnly .bamstat: SN total de secuencias sin procesar: 491  

Hay 1901 secuencias del grupo "minusOnly" y 491 secuencias del grupo "plusOnly", y sé por las estadísticas anteriores que hay 2392 secuencias en total. 1901 + 491 = 2392, entonces sé que no hay lecturas desagrupadas en este archivo.

No estoy seguro de cómo se extrapola de esto si cada una de las lecturas tiene una etiqueta RG. ¿Podría extrapolar un poco cómo deduciría de estas estadísticas?
Agregué información adicional para demostrar cómo se puede hacer esto.
Pierre
2017-10-02 12:07:44 UTC
view on stackexchange narkive permalink

usando samjdk e invocando la función getReadGroup ()

getReadGroup () devuelve el SAMReadGroupRecord del SAMFileHeader para este SAMRecord, o nulo si 1) este registro no tiene etiqueta RG, o 2) el encabezado no contiene el grupo de lectura con el ID dado. o 3) este registro no tiene SAMFileHeader

  java -jar dist / samjdk.jar -e 'return record.getReadGroup () == null;' input.bam  

la salida será un archivo sam / bam con RG faltante / incorrecto.

Bioathlete
2017-10-03 00:22:57 UTC
view on stackexchange narkive permalink

Una solución completa para validar un archivo SAM / BAM sería ValidateSamFile de picard. Esto encontrará muchos otros problemas con los archivos BAM.

http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile



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