Pregunta:
¿Existe una manera fácil de crear un resumen de un archivo VCF (v4.1) con variaciones estructurales?
Kamil S Jaron
2017-05-29 02:40:14 UTC
view on stackexchange narkive permalink

Tengo un montón de archivos vcf (v4.1) con variaciones estructurales de un montón de organismos no modelo (es decir, no hay variantes conocidas). Descubrí que hay bastantes herramientas para manipular archivos vcf como VCFtools, el paquete R vcfR o la biblioteca de Python PyVCF. Sin embargo, ninguno de ellos parece proporcionar un resumen rápido, algo como (preferiblemente categorizado también por tamaño):

  type countDEL xINS yINV z ....  

¿Hay alguna herramienta o función que haya pasado por alto que produzca resúmenes de este estilo?

Sé que el archivo vcf es solo un archivo de texto plano y si analizaré REF y Columnas ALT Debería poder escribir un script que haga el trabajo, pero esperaba poder evitar escribir mi propio analizador.

--- editar ---

Hasta ahora parece que la única herramienta que tiene como objetivo hacer resúmenes (respuesta @gringer) no está funcionando en vcf v4.1. Otras herramientas proporcionarían una solución solo parcial al filtrar cierto tipo de variante. Por lo tanto, acepto mis propias soluciones de analizador perl / R, hasta que haya una herramienta funcional para estadísticas de vcf con variantes estructurales .

¿Qué tipo de resumen estás buscando? ¿Solo recuentos sin procesar (de número de inserciones / eliminaciones)? Solo tengo curiosidad, ya que creo que personalmente me interesaría saber cómo esa variación se propaga realmente a lo largo de un genoma / contig / secuencia espacialmente, ¿sería eso de interés para usted?
Estoy comparando un par de especies relativas. Sin embargo, todavía no he identificado regiones homólogas (la anotación está en construcción hasta ahora), por lo tanto, el posicionamiento de los SV en este momento no es muy relevante. Sin embargo, lo que puedo hacer ya es ver si hay al menos una diferencia en los conteos / tamaños de SV entre especies.
Quizás incluso mejor que el resumen sería obtener el tipo de SV y su tamaño para cada llamada de SV. Entonces podría un resumen en R sería trivial.
Probablemente ya haya visto esto, pero `vcftools` tiene un hermano llamado` bcftools`, que tiene [una función de consulta, que permite a los usuarios consultar un VCF / BCF para extraer campos e información y generar su propio formato] (https : //samtools.github.io/bcftools/bcftools-man.html#query). Puede que no haga exactamente lo que desea, pero puede acercarlo mucho (¿lo suficiente como para necesitar un poco de procesamiento posterior en R?).
Genial, eso suena prometedor. Le echaré un vistazo.
Cinco respuestas:
#1
+8
gringer
2017-05-29 10:16:36 UTC
view on stackexchange narkive permalink

De acuerdo con la página de manual de bcftools , es capaz de producir estadísticas usando el comando bcftools stats . Al ejecutar esto yo mismo, las estadísticas se parecen a lo que está pidiendo:

  # Este archivo fue producido por bcftools stats (1.2-187-g1a55e45 + htslib-1.2.1-256-ga356746) y se puede trazar usando plot-vcfstats. # La línea de comando era: bcftools stats OVLNormalised_STARout_KCCG_called.vcf.gz ## Definición de conjuntos: # ID [2] id [3] nombres de archivo separados por tabulacionesID 0 OVLNormalised_STARout_KCCG_called.vcf.gz , Números de resumen: # SN [2] id [3] tecla [4] valorSN 0 número de muestras: 108SN 0 número de registros: 333SN 0 número de no ALT: 0SN 0 número de SNP: 313SN 0 número de MNP: 0SN 0 número de indeles: 20SN 0 número de otros: 0SN 0 número de sitios multialélicos: 0SN 0 número de sitios SNP multialélicos: 0 # TSTV, transiciones / transversiones: # TSTV [2] id [3] ts [4] tv [5 ] ts / tv [6] ts (1ª ALT) [7] tv (1ª ALT) [8] ts / tv (1ª ALT) TSTV 0 302 11 27,45 302 11 27,45 # SiS, Estadísticas Singleton: ... # IDD, Distribución InDel: # IDD [2] id [3] length (deleti ons negativos) [4] countIDD 0-9 1IDD 0-2 4IDD 0-1 6IDD 0 1 4IDD 0 2 1IDD 0 3 3IDD 0 4 1 # ST, Tipos de sustitución: # ST [2] id [3] type [4] countST 0 A>C 2ST 0 A>G 78ST 0 A>T 2ST 0 C>A 5ST 0 C>G 0ST 0 C>T 66ST 0 G>A 67ST 0 G>C 0ST 0 G>T 1ST 0 T>A 1ST 0 T>C 91 0 T>G 0 # DP, distribución Profundidad # DP [2] Identificación [3 ] bin [4] número de genotipos [5] fracción de genotipos (%) [6] número de sitios [7] fracción de sitios (%) DP 0 >500 0 0.000000 333 100.000000  
Hago casi exactamente lo que me gustaría ver. Pero aparentemente está hecho para SNV, no para SV. Etiquetó incorrectamente los puntos de interrupción como indels.
Está anotando lo que está explícitamente en el archivo VCF, y eso es SNV (e INDEL). Si desea un análisis de variantes estructurales (es decir, a una escala mayor que los nucleótidos individuales), deberá usar algo que haga más que un resumen del archivo VCF. Las inversiones, las eliminaciones a gran escala y los puntos de interrupción no forman parte del formato VCF "estándar".
Son parte de 4.1. Se describe en la documentación: https://samtools.github.io/hts-specs/VCFv4.1.pdf
Bien, ya veo. Supongo que eso lo convierte en un error de bcftools entonces: https: //github.com/samtools/bcftools/issues
Ok, confirmado: bcftools NO admite variantes estructurales https://github.com/samtools/bcftools/issues/623
#2
+6
Kamil S Jaron
2017-05-29 11:22:59 UTC
view on stackexchange narkive permalink

Las soluciones de "mi propio analizador". La información que estaba buscando en parte de la columna INFO , es decir, en las variables SVLEN y SVTYPE.

SV muy rápido tipos + recuentos (por @ user172818 en commnent):

  zcat var.vcf.gz | perl -ne 'imprimir "$ 1 \ n" si / [; \ t] SVTYPE = ([^; \ t] +) /' | ordenar | uniq -c  

tipos SV bastante lentos + conteos + tamaños:

  SV_colnames <- c ('CHROM', 'POS', 'ID', 'REF', 'ALT', 'CUAL', 'FILTRO', 'INFO', 'FORMATO', 'MUESTRA1') ssplit <- function (s, split = '=') {unlist (strsplit (s, split = split))} # nota, las letras mayúsculas solo respetan las convenciones de nomenclatura originales del archivo VCF getSVTYPE <- function (info) {ssplit (grep ("SVTYPE", info, value = T)) [2]} getSVLEN <- function (info ) {SVLEN <- ssplit (grep ("SVLEN", info, value = T)) ifelse (length (SVLEN) == 0, NA, as.numeric (SVLEN [2]))} load_sv <- función (archivo) {vcf_sv_table <- read.table (archivos, stringsAsFactors = F) COLNAMES (vcf_sv_table) <- SV_colnames # # vcf_sv_table posible filtrado <- vcf_sv_table [vcf_sv_table $ FILTRO == 'PASS',] vcf_sv_table_info <- strsplit (vcf_sv_table $ INFO, ' ; ') vcf_sv_table $ SVTYPE <- unlist (lapply (vcf_sv_table_info, getSVTYP E)) vcf_sv_table $ SVLEN <- unlist (lapply (vcf_sv_table_info, getSVLEN)) return (vcf_sv_table)} sv_df <- load_sv ('my_sv_calls.vc_df>) tabla (sv>
Si solo desea contar SVTYPE, esta línea de comando será más rápida: `zcat var.vcf.gz | perl -ne 'imprimir "$ 1 \ n" si / [; \ t] SVTYPE = ([^; \ t] +) /' | ordenar | uniq -c`. Los archivos SV son pequeños, por lo que su Rscript está bien, pero en general, R es muy lento para el procesamiento de texto. Se prefiere Perl / Python / etc cuando se trata de VCF enormes.
La solución R deja las puertas abiertas para un juego adicional con el tipo de longitud / posición / sv, pero definitivamente tiene razón sobre el rendimiento.
#3
+2
eastafri
2017-05-29 10:23:32 UTC
view on stackexchange narkive permalink

Es posible que desee probar Bio-VCF. De la descripción del autor

Bio-vcf es un analizador, filtro y conversor VCF de nueva generación. Bio-vcf no solo es muy rápido para datos de todo el genoma (WGS), sino que también viene con un lenguaje de filtrado, evaluación y reescritura realmente agradable y puede generar cualquier tipo de datos textuales, incluido el encabezado VCF y el contenido en RDF y JSON.

Además, es un analizador muy rápido. El DSL de reescritura puede ayudarlo a personalizar sus filtros y necesidades.

¿Tiene Bio-vcf alguna funcionalidad para la variación estructural como solicitó el OP? Un escaneo rápido de la página de GitHub no revela nada. Si bien muchas de estas herramientas son útiles para filtrar SNV e indels cortos, la mayoría no logran hacer nada útil para los SV.
#4
+2
Medhat Helmy
2019-09-17 01:19:02 UTC
view on stackexchange narkive permalink

Para las variantes estructurales, creo que puede usar SURVIVOR como SURVIVOR stats que fue diseñado específicamente para tal propósito (estadísticas en el archivo SV).

#5
+1
Panwen Wang
2019-10-22 04:01:05 UTC
view on stackexchange narkive permalink

Supongamos que tiene un campo INFO llamado SVTYPE para indicar el tipo de variante de estructura.

Así es como puede obtener las estadísticas fácilmente:

  1. Instalar vcfstats : pip install vcfstats
  2. Defina una macro para extraer la información de SVTYPE , digamos en svtype.py :
  from vcfstats.macros import cat @ catdef SVTYPE (variant): return variant.INFO [ 'SVTYPE']  
  1. Genere las estadísticas:
  vcfstats --vcf <your vcf file> \ --outdir <tel directorio de salida> \ --macro svtype.py \ --formula 'COUNT (1) ~ SVTYPE [DEL, INS, INV ...]' \ # svtypes deseados - title 'Número de tipos de SV'  

Luego, en el directorio de salida, tendrá el archivo de estadísticas llamado Number_of_SV_Types.txt y un gráfico para él también.

Consulte los detalles en https: // git hub.com/pwwang/vcfstats

La ventaja de esta solución es que utiliza un analizador estándar. La desventaja es que requiere varios pasos (a diferencia de `zcat var.vcf.gz | perl -ne 'print" $ 1 \ n "if / [; \ t] SVTYPE = ([^; \ t] +) /' | sort | uniq -c` que hará el resumen de una vez).


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