Pregunta:
¿Cómo concatenar "por cromosoma" -VCF?
ShanZhengYang
2018-01-15 03:38:51 UTC
view on stackexchange narkive permalink

Tengo varios VCF que son VCF que solo contienen información por cromosoma. Es decir, hay un VCF del cromosoma 1 (con solo chr1), un VCF del cromosoma 2 (con solo chr2), etc.

Verifiqué que estos VCF fueran válidos a través de VCFtools , es decir, 

  $ vcf-validator chr1.vcf

que funciona --- estos son VCF válidos que me dieron.

Ahora, me gustaría combinar estos VCF en un VCF.

Intenté ingenuamente la siguiente operación cat : 

  $ cat chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf > total_chroms.vcf

Sin embargo, esto no funciona correctamente 

  $ vcf-validator total_chroms.vcfLa etiqueta de encabezado 'contig' no está presente para CHROM = chr1. (No es obligatorio pero muy recomendable). No se pudo analizar la línea, número incorrecto de columnas: [## fileformat = VCFv4.2 \ n] en /path/vcftools-0.1.14/perl/Vcf.pm línea 172, <__ANONIO__> línea 191016. Vcf :: throw ('Vcf4_2 = HASH (0x1ae7208)', 'No se pudo analizar la línea, número incorrecto de columnas: [## filefor ...') llamado en /path/vcftools-0.1.14/perl/ Vcf.pm línea 335 Vcf :: next_data_array ('Vcf4_2 = HASH (0x1ae7208)') llamado en /path/vcftools-0.1.14/perl/Vcf.pm línea 3457 Vcf4_1 :: next_data_array ('Vcf4_2 = HASH7 (0x1a)' ) llamado en /path/vcftools-0.1.14/perl/Vcf.pm línea 2574 VcfReader :: run_validation ('Vcf4_2 = HASH (0x1ae7208)') llamado en /path/vcftools-0.1.14//bin/vcf-validator línea 60 main :: do_validation ('HASH (0x16ada68)') llamado en /path/vcftools-0.1.14//bin/vcf-validator línea 14 $

¿Qué herramientas están disponibles para fusionar estos VCF juntos en un VCF total?

Tres respuestas:
Bioathlete
2018-01-15 05:25:11 UTC
view on stackexchange narkive permalink

Recomendaría bcftools concat. No puede simplemente cat juntos porque cada archivo tiene una sección de encabezado. El comando bcftools se encargará de todo eso por usted. Cada archivo vcf debe ordenarse antes de llamar a concat

bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf

Es un poco extraño. `vcf-validator chr1.vcf` no muestra problemas. Creo que la salida de `bcftools concat -o chroms12.vcf chr1.vcf chr2.vcf` muestra un error:` FIXME: nombre de secuencia chr1 en chr1.vcf`
@ShanZhengYang los errores parecen indicar que tiene `chr1` en el campo CHROM para algunas filas, pero no tiene la línea correspondiente en el encabezado. ¿Puede comprobar si eso es cierto en el archivo `chr1.vcf`? ¿Quizás al validador de alguna manera le falta este error?
@juod El encabezado VCF me parece correcto. `#CHROM POS ID REF ALT QUAL FILTER INFO`
@juod "pero no hay una línea correspondiente en el encabezado" Quizás no entiendo esto. Además del código pegado arriba, no hay otra información "cromática" en el encabezado
De acuerdo con la especificación VCF https://samtools.github.io/hts-specs/VCFv4.1.pdf, debe haber líneas de encabezado para cada cromosoma, esto es del primer ejemplo en los documentos vinculados - `## contig = `
@Bioathlete Ya veo. Creo que lo entiendo ahora --- este problema se ha discutido antes: https://github.com/samtools/bcftools/issues/326 Sé que estos VCF son de hg38. ¿Podría simplemente agregar la línea `## contig = `, etc.a los VCF en cuestión? Solo tomaría las longitudes de los cromosomas de hg38
@ShanZhengYang prueba eso, debería funcionar. Ni siquiera creo que sea necesario proporcionar la longitud, pero podrías experimentar y hacérnoslo saber.
@juod ¡Lo anterior funcionó! ¡Gracias por la ayuda!
terdon
2018-01-15 21:13:06 UTC
view on stackexchange narkive permalink

La mejor herramienta para el trabajo es probablemente bcftools como sugerido por Bioathlete, pero también puede hacerlo manualmente. Solo necesita recopilar todas las líneas de encabezado de todos los archivos vcf, eliminar los duplicados y luego imprimir todos los encabezados + los datos reales en el nuevo: 

  grep '^ ##' chr * vcf | ordenar | uniq > all.vcfgrep -m1 '^ # CHR' chr1.vcf >> all.vcf ## Obtiene la línea de encabezado chr grep -v '^ #' chr * vcf >> all.vcf
Parece que obtengo un encabezado VCF mal formado con los comandos anteriores
@ShanZhengYang malformado ¿cómo? ¿Podría editar su pregunta y mostrarnos un ejemplo de sus archivos?
Entonces, veo que los comandos `grep` anteriores dan el mismo comportamiento con todos mis VCF estándar. Obtendré algo como esto en `all.vcf`:` chr10.vcf: ## contig = chr10.vcf: ## fileDate = 20121011chr10.vcf: ## fileformat = VCFv4. 2chr10.vcf: ## reference = path / reference.fachr10.vcf: ## source = foobarchr11.vcf: ## contig = chr11.vcf: ## fileDate = 20121011chr11.vcf: # # fileformat = VCFv4.2chr11.vcf: ## referencia = ruta / referencia.fachr11.vcf: ## fuente = foobar .... `
Pierre
2018-01-16 00:30:09 UTC
view on stackexchange narkive permalink

use picard GatherVcf : http://broadinstitute.github.io/picard/command-line-overview.html#GatherVcfs

Recopila varios archivos VCF de una operación de dispersión en un solo archivo VCF. Los archivos de entrada deben proporcionarse en orden genómico y no deben tener eventos en posiciones superpuestas.

El error que obtengo aquí es `Excepción en el hilo" principal "picard.PicardException: para indexar la entrada VCF resultante, los VCF deben contener ## líneas de contig. Es cierto, mis VCF "por cromosoma" no contienen una línea `## contig` en el encabezado ... pero no estoy seguro de que eso importe
Editar arriba: Aparentemente, este problema se ha discutido antes: github.com/samtools/bcftools/issues/326 Pregunta: ¿cómo podría agregar estas líneas de contig si no conozco esta información? Sé que estos VCF son de hg38. ¿Podría simplemente agregar la línea `## contig = `, etc.a los VCF en cuestión? Solo tomaría las longitudes de los cromosomas de hg38
"pero no estoy seguro de que eso importe" comprueba el orden de los contigs usando el dict. use UpdateVcfSequenceDictionary https://broadinstitute.github.io/picard/command-line-overview.html
"comprueba el orden de los contigs usando el dict. use UpdateVcfSequenceDictionary". Ah, ahora lo entiendo. ¡Gracias por la ayuda!


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