Questão:
Existe uma maneira fácil de criar um resumo de um arquivo VCF (v4.1) com variações estruturais?
Kamil S Jaron
2017-05-29 02:40:14 UTC
view on stackexchange narkive permalink

Eu tenho um monte de arquivos vcf (v4.1) com variações estruturais de um monte de organismos não-modelo (ou seja, não há variantes conhecidas). Descobri que existem várias ferramentas para manipular arquivos vcf como VCFtools, pacote R vcfR ou biblioteca python PyVCF. No entanto, nenhum deles parece fornecer um resumo rápido, algo como (de preferência categorizado por tamanho também):

  type countDEL xINS yINV z ....  

Existe alguma ferramenta ou função que esqueci que produz resumos deste estilo?

Eu sei que o arquivo vcf é apenas um arquivo de texto simples e se irei dissecar REF e Colunas ALT Devo ser capaz de escrever um script que fará o trabalho, mas esperava poder evitar escrever meu próprio analisador.

--- editar ---

Até agora, parece que a única ferramenta que visa fazer resumos (resposta @gringer) não está funcionando no vcf v4.1. Outras ferramentas forneceriam solução apenas parcial, filtrando certo tipo de variante. Portanto, eu aceito minhas próprias soluções perl / R de analisador, até que haja uma ferramenta funcional para estatísticas de vcf com variantes estruturais .

Que tipo de resumo você está procurando? Apenas contagens brutas (do número de inserções / exclusões)? Estou apenas curioso, pois acho que pessoalmente estaria interessado em saber como essa variação é realmente espalhada ao longo de um genoma / contig / sequência espacialmente, isso seria do seu interesse?
Estou comparando algumas espécies relativas. No entanto, ainda não identifiquei regiões homólogas (a anotação ainda está em construção), portanto, o posicionamento dos SVs não é muito relevante no momento. No entanto, o que já posso fazer é verificar se existe pelo menos uma diferença nas contagens / tamanhos dos SVs entre as espécies.
Talvez até melhor do que o resumo seria obter o tipo de SV e seu tamanho para cada chamada de SV. Então eu poderia um resumo em R seria trivial.
Você provavelmente já viu isso, mas `vcftools` tem um irmão chamado` bcftools`, que tem [uma função de consulta, que permite aos usuários consultar um VCF / BCF para extrair campos e informações e gerar seu próprio formato] (https : //samtools.github.io/bcftools/bcftools-man.html#query). Pode não fazer exatamente o que você deseja, mas pode chegar muito perto (o suficiente para talvez apenas precisar de um pouco de pós-processamento em R?).
Ótimo, isso parece promissor. Eu vou dar uma olhada.
Cinco respostas:
#1
+8
gringer
2017-05-29 10:16:36 UTC
view on stackexchange narkive permalink

De acordo com a página do manual bcftools , ele é capaz de produzir estatísticas usando o comando bcftools stats . Executando isso sozinho, as estatísticas parecem com o que você está pedindo:

  # Este arquivo foi produzido por bcftools stats (1.2-187-g1a55e45 + htslib-1.2.1-256-ga356746) e pode ser plotado usando plot-vcfstats. # A linha de comando era: bcftools stats OVLNormalised_STARout_KCCG_called.vcf.gz ## Definição de conjuntos: # ID [2] id [3] nomes de arquivo separados por tabulaçãoID 0 OVLNormalised_STARout_KCCG_called.vcf.gz , Números de resumo: # SN [2] id [3] chave [4] valorSN 0 número de amostras: 108SN 0 número de registros: 333SN 0 número de não-ALTs: 0SN 0 número de SNPs: 313SN 0 número de MNPs: 0SN 0 número de indels: 20SN 0 número de outros: 0SN 0 número de sites multialélicos: 0SN 0 número de sites SNP multialélicos: 0 # TSTV, transições / transversões: # 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, estatísticas de singleton: ... # IDD, Distribuição InDel: # IDD [2] id [3] comprimento (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 substituição: # ST [2] id [3] tipo [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 91ST 0 T>G 0 # DP, distribuição Profundidade # DP [2] id [3 ] bin [4] número de genótipos [5] fração de genótipos (%) [6] número de locais [7] fração de locais (%) DP 0 >500 0 0,000000 333 100,000000  
Eu faço quase exatamente o que gostaria de ver. Mas aparentemente é feito para SNVs não, SVs. Ele rotulou erroneamente os pontos de interrupção como indels.
Ele está anotando o que está explicitamente no arquivo VCF, que são SNVs (e INDELs). Se você quiser uma análise de variante estrutural (ou seja, em uma escala maior do que um único nucleotídeo), você precisará usar algo que faça mais do que um resumo do arquivo VCF. Inversões, exclusões em grande escala e pontos de interrupção não fazem parte do formato VCF "padrão".
Eles fazem parte do 4.1. Ele está descrito na documentação: https://samtools.github.io/hts-specs/VCFv4.1.pdf
Certo, entendo. Acho que isso o torna um bug do bcftools: https: //github.com/samtools/bcftools/issues
Ok, confirmado: bcftools NÃO suporta variantes estruturais https://github.com/samtools/bcftools/issues/623
#2
+6
Kamil S Jaron
2017-05-29 11:22:59 UTC
view on stackexchange narkive permalink

As soluções "meu próprio analisador". A informação que procuro na coluna INFO , nomeadamente nas variáveis ​​ SVLEN e SVTYPE.

SV muito rápido tipos + contagens (por @ user172818 no comnent):

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

tipos SV bastante lentos + contagens + tamanhos:

  SV_colnames <- c ('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1') ssplit <- function (s, split = '=') {unlist (strsplit (s, split = split))} # note, as letras maiúsculas apenas respeitam as convenções de nomenclatura originais do filegetSVTYPE VCF função <- (info) {ssplit (grep ("SVTYPE", info, valor = T)) [2]} getSVLEN <- função (info ) {SVLEN <- ssplit (grep ("SVLEN", info, valor = T)) ifelse (comprimento (SVLEN) == 0, NA, as.numeric (SVLEN [2]))} load_sv função <- (arquivo) {vcf_sv_table <- read.table (arquivos, stringsAsFactors = F) COLNAMES (vcf_sv_table) <- SV_colnames # possível filtragem # vcf_sv_table <- vcf_sv_table [vcf_sv_table $ FILTER == '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.vcf') tabela (sv_code> SVPE $)
Se você deseja apenas contar SVTYPE, esta linha de comando será mais rápida: `zcat var.vcf.gz | perl -ne 'imprimir "$ 1 \ n" if / [; \ t] SVTYPE = ([^; \ t] +) /' | classificar | uniq -c`. Os arquivos SV são minúsculos, então seu Rscript é bom, mas geralmente, R é muito lento para processamento de texto. Perl / Python / etc é o preferido quando você está lidando com VCFs enormes.
A solução R deixa as portas abertas para brincadeiras adicionais com o tipo length / position / sv, mas você está definitivamente certo sobre a performance.
#3
+2
eastafri
2017-05-29 10:23:32 UTC
view on stackexchange narkive permalink

Você pode querer experimentar Bio-VCF. Da descrição do autor

Bio-vcf é um analisador, filtro e conversor de VCF de nova geração. Bio-vcf não é apenas muito rápido para dados de todo o genoma (WGS), ele também vem com uma linguagem muito boa de filtragem, avaliação e reescrita e pode produzir qualquer tipo de dados textuais, incluindo cabeçalho VCF e conteúdo em RDF e JSON.

Além disso, é um analisador muito rápido. A reescrita DSL pode ajudá-lo a personalizar sua filtragem e necessidades.

O Bio-vcf possui alguma funcionalidade para variação estrutural conforme solicitado pelo OP? Uma rápida verificação da página do GitHub não revela nada. Embora muitas dessas ferramentas sejam úteis para filtrar SNVs e indels curtos, a maioria não faz nada de importante para os SVs.
#4
+2
Medhat Helmy
2019-09-17 01:19:02 UTC
view on stackexchange narkive permalink

Para variantes estruturais, acho que você pode usar SURVIVOR como SURVIVOR stats , que foi projetado especificamente para essa finalidade (estatísticas no arquivo SV).

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

Digamos que você tenha um campo INFO chamado SVTYPE para indicar o tipo de variante da estrutura.

Veja como obter as estatísticas facilmente:

  1. Instalar vcfstats : pip install vcfstats
  2. Defina uma macro para extrair as informações de SVTYPE , digamos em svtype.py :
  from vcfstats.macros import cat @ catdef SVTYPE (variante): return variant.INFO [ 'SVTYPE']  
  1. Gere as estatísticas:
  vcfstats --vcf <your vcf file> \ --outdir <the diretório de saída> \ --macro svtype.py \ --formula 'COUNT (1) ~ SVTYPE [DEL, INS, INV ... svtypes title 'Número de tipos de SV'  

Então, no diretório de saída você terá o arquivo de estatísticas chamado Number_of_SV_Types.txt e um gráfico para ele também.

Verifique os detalhes em https: // git hub.com/pwwang/vcfstats

A vantagem desta solução é que ela usa um analisador padrão. A desvantagem é que requer várias etapas (ao contrário de `zcat var.vcf.gz | perl -ne 'print" $ 1 \ n "if / [; \ t] SVTYPE = ([^; \ t] +) /' | sort | uniq -c` que fará o resumo de uma vez).


Estas perguntas e respostas foram traduzidas automaticamente do idioma inglês.O conteúdo original está disponível em stackexchange, que agradecemos pela licença cc by-sa 3.0 sob a qual é distribuído.
Loading...