Parte 1: Visão geral do método¶
Tradução assistida por IA - saiba mais e sugira melhorias
Chamada de variantes é um método de análise genômica que visa identificar variações em uma sequência genômica em relação a um genoma de referência. Aqui vamos usar ferramentas e métodos projetados para chamar variantes germinativas curtas, ou seja, SNPs e indels, em dados de sequenciamento de genoma completo.

Um pipeline completo de chamada de variantes normalmente envolve muitas etapas, incluindo mapeamento para a referência (às vezes chamado de alinhamento genômico) e filtragem e priorização de variantes. Para simplificar, neste curso vamos focar apenas na parte de chamada de variantes.
Métodos¶
Vamos mostrar duas maneiras de aplicar chamada de variantes a amostras de sequenciamento de genoma completo para identificar SNPs e indels germinativos. Primeiro começaremos com uma abordagem por amostra simples que chama variantes independentemente de cada amostra. Depois mostraremos uma abordagem de chamada conjunta mais sofisticada que analisa múltiplas amostras juntas, produzindo resultados mais precisos e informativos.
Antes de mergulharmos na escrita de qualquer código de fluxo de trabalho para qualquer uma das abordagens, vamos testar os comandos manualmente em alguns dados de teste.
Conjunto de dados¶
Fornecemos os seguintes dados e recursos relacionados:
- Um genoma de referência consistindo de uma pequena região do cromossomo humano 20 (de hg19/b37) e seus arquivos acessórios (índice e dicionário de sequência).
- Três amostras de sequenciamento de genoma completo correspondentes a um trio familiar (mãe, pai e filho), que foram reduzidas a uma pequena fatia de dados no cromossomo 20 para manter os tamanhos de arquivo pequenos. Estes são dados de sequenciamento Illumina de leituras curtas que já foram mapeados para o genoma de referência, fornecidos no formato BAM (Binary Alignment Map, uma versão comprimida de SAM, Sequence Alignment Map).
- Uma lista de intervalos genômicos, ou seja, coordenadas no genoma onde nossas amostras têm dados adequados para chamar variantes, fornecida no formato BED.
Software¶
As duas principais ferramentas envolvidas são Samtools, um kit de ferramentas amplamente usado para manipular arquivos de alinhamento de sequências, e GATK (Genome Analysis Toolkit), um conjunto de ferramentas para descoberta de variantes desenvolvido no Broad Institute.
Essas ferramentas não estão instaladas no ambiente GitHub Codespaces, então vamos usá-las via contêineres obtidos através do serviço Seqera Containers (veja Hello Containers).
Dica
Certifique-se de estar no diretório nf4-science/genomics para que a última parte do caminho mostrada quando você digita pwd seja genomics.
1. Chamada de variantes por amostra¶
A chamada de variantes por amostra processa cada amostra independentemente: o chamador de variantes examina os dados de sequenciamento de uma amostra por vez e identifica posições onde a amostra difere da referência.
Nesta seção testamos os dois comandos que compõem a abordagem de chamada de variantes por amostra: indexar um arquivo BAM com Samtools e chamar variantes com GATK HaplotypeCaller. Estes são os comandos que vamos encapsular em um fluxo de trabalho Nextflow na Parte 2 deste curso.
- Gerar um arquivo de índice para um arquivo de entrada BAM usando Samtools
- Executar o GATK HaplotypeCaller no arquivo BAM indexado para gerar chamadas de variantes por amostra em VCF (Variant Call Format)
Começamos testando os dois comandos em apenas uma amostra.
1.1. Indexar um arquivo de entrada BAM com Samtools¶
Arquivos de índice são uma característica comum dos formatos de arquivo de bioinformática; eles contêm informações sobre a estrutura do arquivo principal que permite que ferramentas como GATK acessem um subconjunto dos dados sem ter que ler o arquivo inteiro. Isso é importante por causa de quão grandes esses arquivos podem ficar.
Arquivos BAM são frequentemente fornecidos sem um índice, então o primeiro passo em muitos fluxos de trabalho de análise é gerar um usando samtools index.
Vamos baixar um contêiner Samtools, iniciá-lo interativamente e executar o comando samtools index em um dos arquivos BAM.
1.1.1. Baixar o contêiner Samtools¶
Execute o comando docker pull para baixar a imagem do contêiner Samtools:
Saída do comando
1.20--b5dfbd93de237464: Pulling from library/samtools
6360b3717211: Pull complete
2ec3f7ad9b3c: Pull complete
7716ca300600: Pull complete
4f4fb700ef54: Pull complete
8c61d418774c: Pull complete
03dae77ff45c: Pull complete
aab7f787139d: Pull complete
4f4fb700ef54: Pull complete
837d55536720: Pull complete
897362c12ca7: Pull complete
3893cbe24e91: Pull complete
d1b61e94977b: Pull complete
c72ff66fb90f: Pull complete
0e0388f29b6d: Pull complete
Digest: sha256:bbfc45b4f228975bde86cba95e303dd94ecf2fdacea5bfb2e2f34b0d7b141e41
Status: Downloaded newer image for community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
community.wave.seqera.io/library/samtools:1.20--b5dfbd93de237464
Se você não baixou esta imagem antes, pode levar um minuto para completar. Uma vez concluído, você tem uma cópia local da imagem do contêiner.
1.1.2. Iniciar o contêiner Samtools interativamente¶
Para executar o contêiner interativamente, use docker run com as flags -it.
A opção -v ./data:/data monta o diretório local data dentro do contêiner para que as ferramentas possam acessar os arquivos de entrada.
Você notará que seu prompt muda para algo como (base) root@a1b2c3d4e5f6:/tmp#, indicando que você está agora dentro do contêiner.
Verifique que você pode ver os arquivos de dados de sequência em /data/bam:
Com isso, você está pronto para testar seu primeiro comando.
1.1.3. Executar o comando de indexação¶
A documentação do Samtools nos fornece a linha de comando para executar para indexar um arquivo BAM.
Precisamos apenas fornecer o arquivo de entrada; a ferramenta gerará automaticamente um nome para a saída anexando .bai ao nome do arquivo de entrada.
Execute o comando samtools index em um arquivo de dados:
O comando não produz nenhuma saída no terminal, mas você deve agora ver um arquivo chamado reads_mother.bam.bai no mesmo diretório que o arquivo de entrada BAM original.
Conteúdo do diretório
Isso completa o teste da primeira etapa.
1.1.4. Sair do contêiner Samtools¶
Para sair do contêiner, digite exit.
Seu prompt deve agora estar de volta ao que era antes de você iniciar o contêiner.
1.2. Chamar variantes com GATK HaplotypeCaller¶
Queremos executar o comando gatk HaplotypeCaller no arquivo BAM que acabamos de indexar.
1.2.1. Baixar o contêiner GATK¶
Primeiro, vamos executar o comando docker pull para baixar a imagem do contêiner GATK:
Saída do comando
Algumas camadas mostram Already exists porque são compartilhadas com a imagem do contêiner Samtools que baixamos anteriormente.
4.5.0.0--730ee8817e436867: Pulling from library/gatk4
6360b3717211: Already exists
2ec3f7ad9b3c: Already exists
7716ca300600: Already exists
4f4fb700ef54: Already exists
8c61d418774c: Already exists
03dae77ff45c: Already exists
aab7f787139d: Already exists
4f4fb700ef54: Already exists
837d55536720: Already exists
897362c12ca7: Already exists
3893cbe24e91: Already exists
d1b61e94977b: Already exists
e5c558f54708: Pull complete
087cce32d294: Pull complete
Digest: sha256:e33413b9100f834fcc62fd5bc9edc1e881e820aafa606e09301eac2303d8724b
Status: Downloaded newer image for community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
community.wave.seqera.io/library/gatk4:4.5.0.0--730ee8817e436867
Isso deve ser mais rápido que o primeiro download porque as duas imagens de contêiner compartilham a maioria de suas camadas.
1.2.2. Iniciar o contêiner GATK interativamente¶
Inicie o contêiner GATK interativamente com o diretório de dados montado, assim como fizemos para Samtools.
Seu prompt muda para indicar que você está agora dentro do contêiner GATK.
1.2.3. Executar o comando de chamada de variantes¶
A documentação do GATK nos fornece a linha de comando para executar para realizar chamada de variantes em um arquivo BAM.
Precisamos fornecer o arquivo de entrada BAM (-I) assim como o genoma de referência (-R), um nome para o arquivo de saída (-O) e uma lista de intervalos genômicos para analisar (-L).
No entanto, não precisamos especificar o caminho para o arquivo de índice; a ferramenta procurará automaticamente por ele no mesmo diretório, baseado na convenção estabelecida de nomenclatura e colocalização.
O mesmo se aplica aos arquivos acessórios do genoma de referência (arquivos de índice e dicionário de sequência, *.fai e *.dict).
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_mother.bam \
-O /data/vcf/reads_mother.vcf \
-L /data/ref/intervals.bed
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.vcf -L /data/ref/intervals.bed
00:27:50.687 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
00:27:50.854 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.858 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
00:27:50.858 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
00:27:50.858 INFO HaplotypeCaller - Executing as root@a1fe8ff42d07 on Linux v6.10.14-linuxkit amd64
00:27:50.858 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
00:27:50.859 INFO HaplotypeCaller - Start Date/Time: February 8, 2026 at 12:27:50 AM GMT
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.859 INFO HaplotypeCaller - ------------------------------------------------------------
00:27:50.861 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
00:27:50.861 INFO HaplotypeCaller - Picard Version: 3.1.1
00:27:50.861 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
00:27:50.862 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
00:27:50.863 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
00:27:50.864 INFO HaplotypeCaller - Deflater: IntelDeflater
00:27:50.864 INFO HaplotypeCaller - Inflater: IntelInflater
00:27:50.864 INFO HaplotypeCaller - GCS max retries/reopens: 20
00:27:50.864 INFO HaplotypeCaller - Requester pays: disabled
00:27:50.865 INFO HaplotypeCaller - Initializing engine
00:27:50.991 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
00:27:51.016 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
00:27:51.029 INFO HaplotypeCaller - Done initializing engine
00:27:51.040 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
00:27:51.042 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
00:27:51.042 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
00:27:51.046 INFO HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
00:27:51.063 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
00:27:51.085 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
00:27:51.086 INFO IntelPairHmm - Available threads: 10
00:27:51.086 INFO IntelPairHmm - Requested threads: 4
00:27:51.086 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
00:27:51.128 INFO ProgressMeter - Starting traversal
00:27:51.136 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
00:27:51.882 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
00:27:52.969 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
00:27:52.971 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1145.7
00:27:52.971 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
00:27:52.976 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.003346916
00:27:52.976 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.045731709
00:27:52.977 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
00:27:52.981 INFO HaplotypeCaller - Shutting down engine
[February 8, 2026 at 12:27:52 AM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=203423744
A saída de log é muito verbosa, então destacamos as linhas mais relevantes no exemplo acima.
Os arquivos de saída, reads_mother.vcf e seu arquivo de índice, reads_mother.vcf.idx, são criados dentro do seu diretório de trabalho no contêiner.
O arquivo VCF contém as chamadas de variantes, como veremos em um minuto, e o arquivo de índice tem a mesma função que o arquivo de índice BAM, permitir que ferramentas busquem e recuperem subconjuntos de dados sem carregar o arquivo inteiro.
Como VCF é um formato de texto e este é um arquivo de teste pequeno, você pode executar cat reads_mother.vcf para abri-lo e visualizar seu conteúdo.
Se você rolar de volta até o início do arquivo, encontrará um cabeçalho composto de muitas linhas de metadados, seguido por uma lista de chamadas de variantes, uma por linha.
Conteúdo do arquivo (resumido)
No exemplo de saída acima, destacamos a última linha do cabeçalho, que fornece os nomes das colunas para os dados tabulares que seguem. Cada linha de dados descreve uma possível variante identificada nos dados de sequenciamento da amostra. Para orientação sobre como interpretar o formato VCF, veja este artigo útil.
1.2.4. Mover os arquivos de saída¶
Qualquer coisa que permaneça dentro do contêiner ficará inacessível para trabalhos futuros.
O arquivo de índice BAM foi criado diretamente no diretório /data/bam no sistema de arquivos montado, mas não o arquivo VCF e seu índice, então precisamos mover esses dois manualmente.
Conteúdo do diretório
Uma vez feito isso, todos os arquivos estão agora acessíveis no seu sistema de arquivos normal.
1.2.5. Sair do contêiner GATK¶
Para sair do contêiner, digite exit.
Seu prompt deve voltar ao normal. Isso conclui o teste de chamada de variantes por amostra.
Escreva como um fluxo de trabalho!
Sinta-se à vontade para seguir direto para a Parte 2 se quiser começar a implementar esta análise como um fluxo de trabalho Nextflow. Você só precisará voltar para completar a segunda rodada de testes antes de seguir para a Parte 3.
2. Chamada conjunta em uma coorte¶
A abordagem de chamada de variantes que acabamos de usar gera chamadas de variantes por amostra. Isso é bom para olhar variantes de cada amostra isoladamente, mas produz informações limitadas. Muitas vezes é mais interessante olhar como as chamadas de variantes diferem entre múltiplas amostras. O GATK oferece um método alternativo chamado chamada conjunta de variantes para este propósito.
A chamada conjunta de variantes envolve gerar um tipo especial de saída de variante chamado GVCF (para Genomic VCF) para cada amostra, depois combinar os dados GVCF de todas as amostras e executar uma análise estatística de 'genotipagem conjunta'.

O que é especial sobre o GVCF de uma amostra é que ele contém registros resumindo estatísticas de dados de sequência sobre todas as posições na área alvo do genoma, não apenas as posições onde o programa encontrou evidência de variação. Isso é crítico para o cálculo de genotipagem conjunta (leitura adicional).
O GVCF é produzido pelo GATK HaplotypeCaller, a mesma ferramenta que acabamos de testar, com um parâmetro adicional (-ERC GVCF).
Combinar os GVCFs é feito com GATK GenomicsDBImport, que combina as chamadas por amostra em um armazenamento de dados (análogo a um banco de dados).
A análise de 'genotipagem conjunta' propriamente dita é então feita com GATK GenotypeGVCFs.
Aqui testamos os comandos necessários para gerar GVCFs e executar genotipagem conjunta. Estes são os comandos que vamos encapsular em um fluxo de trabalho Nextflow na Parte 3 deste curso.
- Gerar um arquivo de índice para cada arquivo de entrada BAM usando Samtools
- Executar o GATK HaplotypeCaller em cada arquivo de entrada BAM para gerar um GVCF de chamadas de variantes genômicas por amostra
- Coletar todos os GVCFs e combiná-los em um armazenamento de dados GenomicsDB
- Executar genotipagem conjunta no armazenamento de dados GVCF combinado para produzir um VCF em nível de coorte
Agora precisamos testar todos esses comandos, começando com a indexação de todos os três arquivos BAM.
2.1. Indexar arquivos BAM para todas as três amostras¶
Na primeira seção acima, indexamos apenas um arquivo BAM. Agora precisamos indexar todas as três amostras para que o GATK HaplotypeCaller possa processá-las.
2.1.1. Iniciar o contêiner Samtools interativamente¶
Já baixamos a imagem do contêiner Samtools, então podemos iniciá-lo diretamente:
Seu prompt muda para indicar que você está dentro do contêiner, com o diretório de dados montado como antes.
2.1.2. Executar o comando de indexação em todas as três amostras¶
Execute o comando de indexação em cada um dos três arquivos BAM:
samtools index /data/bam/reads_mother.bam
samtools index /data/bam/reads_father.bam
samtools index /data/bam/reads_son.bam
Conteúdo do diretório
Isso deve produzir os arquivos de índice no mesmo diretório que os arquivos BAM correspondentes.
2.1.3. Sair do contêiner Samtools¶
Para sair do contêiner, digite exit.
Seu prompt deve estar de volta ao normal.
2.2. Gerar GVCFs para todas as três amostras¶
Para executar a etapa de genotipagem conjunta, precisamos de GVCFs para todas as três amostras.
2.2.1. Iniciar o contêiner GATK interativamente¶
Já baixamos a imagem do contêiner GATK anteriormente, então podemos iniciá-lo diretamente:
Seu prompt muda para indicar que você está dentro do contêiner GATK.
2.2.2. Executar o comando de chamada de variantes com a opção GVCF¶
Para produzir um VCF genômico (GVCF), adicionamos a opção -ERC GVCF ao comando base, que ativa o modo GVCF do HaplotypeCaller.
Também mudamos a extensão do arquivo de saída de .vcf para .g.vcf.
Tecnicamente isso não é um requisito, mas é uma convenção fortemente recomendada.
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_mother.bam \
-O reads_mother.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_mother.bam -O reads_mother.g.vcf -L /data/ref/intervals.bed -ERC GVCF
16:51:00.620 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
16:51:00.749 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.751 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
16:51:00.751 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
16:51:00.751 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
16:51:00.751 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
16:51:00.752 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 4:51:00 PM GMT
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - ------------------------------------------------------------
16:51:00.752 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
16:51:00.753 INFO HaplotypeCaller - Picard Version: 3.1.1
16:51:00.753 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
16:51:00.753 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
16:51:00.754 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
16:51:00.754 INFO HaplotypeCaller - Deflater: IntelDeflater
16:51:00.754 INFO HaplotypeCaller - Inflater: IntelInflater
16:51:00.754 INFO HaplotypeCaller - GCS max retries/reopens: 20
16:51:00.754 INFO HaplotypeCaller - Requester pays: disabled
16:51:00.755 INFO HaplotypeCaller - Initializing engine
16:51:00.893 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
16:51:00.905 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
16:51:00.910 INFO HaplotypeCaller - Done initializing engine
16:51:00.912 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
16:51:00.917 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
16:51:00.919 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
16:51:00.919 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
16:51:00.923 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
16:51:00.923 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
16:51:00.933 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
16:51:00.945 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
16:51:00.945 INFO IntelPairHmm - Available threads: 4
16:51:00.945 INFO IntelPairHmm - Requested threads: 4
16:51:00.945 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
16:51:00.984 INFO ProgressMeter - Starting traversal
16:51:00.985 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
16:51:01.452 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
16:51:02.358 INFO HaplotypeCaller - 7 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
7 total reads filtered out of 1867 reads processed
16:51:02.359 INFO ProgressMeter - 20_10037292_10066351:13499 0.0 35 1529.5
16:51:02.359 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
16:51:02.361 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0022800000000000003
16:51:02.361 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.061637120000000004
16:51:02.361 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
16:51:02.362 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 4:51:02 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=257949696
Isso cria o arquivo de saída GVCF reads_mother.g.vcf no diretório de trabalho atual no contêiner, assim como seu arquivo de índice, reads_mother.g.vcf.idx.
Se você executar head -200 reads_mother.g.vcf para visualizar as primeiras 200 linhas do conteúdo do arquivo, verá que é muito mais longo que o VCF equivalente que geramos na primeira seção, e a maioria das linhas parece bem diferente do que vimos no VCF.
Conteúdo do arquivo (resumido)
| reads_mother.g.vcf | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 | |
Destacamos mais uma vez a última linha do cabeçalho, assim como as três primeiras chamadas de variantes 'propriamente ditas' no arquivo.
Você notará que as linhas de chamada de variantes estão intercaladas com muitas linhas não-variantes, que representam regiões não-variantes onde o chamador de variantes não encontrou evidência de variação. Como mencionado brevemente acima, isso é o que é especial sobre o modo GVCF de chamada de variantes: o chamador de variantes captura algumas estatísticas descrevendo seu nível de confiança na ausência de variação. Isso torna possível distinguir entre dois casos muito diferentes: (1) há dados de boa qualidade mostrando que a amostra é homozigota-referência, e (2) não há dados bons suficientes disponíveis para fazer uma determinação de qualquer forma.
Em um GVCF como este, normalmente há muitas dessas linhas não-variantes, com um número menor de registros de variantes espalhados entre elas.
2.2.3. Repetir o processo nas outras duas amostras¶
Agora vamos gerar GVCFs para as duas amostras restantes executando os comandos abaixo, um após o outro.
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_father.bam \
-O reads_father.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_father.bam -O reads_father.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:28:30.677 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:28:30.801 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.803 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:28:30.804 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:28:30.804 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:28:30.804 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:28:30.804 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:28:30 PM GMT
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.804 INFO HaplotypeCaller - ------------------------------------------------------------
17:28:30.805 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:28:30.805 INFO HaplotypeCaller - Picard Version: 3.1.1
17:28:30.805 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:28:30.806 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:28:30.806 INFO HaplotypeCaller - Deflater: IntelDeflater
17:28:30.807 INFO HaplotypeCaller - Inflater: IntelInflater
17:28:30.807 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:28:30.807 INFO HaplotypeCaller - Requester pays: disabled
17:28:30.807 INFO HaplotypeCaller - Initializing engine
17:28:30.933 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:28:30.946 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:28:30.951 INFO HaplotypeCaller - Done initializing engine
17:28:30.953 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:28:30.957 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:28:30.959 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:28:30.960 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:28:30.963 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:28:30.963 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:28:30.972 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:28:30.987 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:28:30.987 INFO IntelPairHmm - Available threads: 4
17:28:30.987 INFO IntelPairHmm - Requested threads: 4
17:28:30.987 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:28:31.034 INFO ProgressMeter - Starting traversal
17:28:31.034 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:28:31.570 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:28:32.865 INFO HaplotypeCaller - 9 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
9 total reads filtered out of 2064 reads processed
17:28:32.866 INFO ProgressMeter - 20_10037292_10066351:13338 0.0 38 1245.2
17:28:32.866 INFO ProgressMeter - Traversal complete. Processed 38 total regions in 0.0 minutes.
17:28:32.868 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0035923200000000004
17:28:32.868 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.10765202500000001
17:28:32.868 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.03 sec
17:28:32.869 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:28:32 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=299892736
gatk HaplotypeCaller \
-R /data/ref/ref.fasta \
-I /data/bam/reads_son.bam \
-O reads_son.g.vcf \
-L /data/ref/intervals.bed \
-ERC GVCF
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar HaplotypeCaller -R /data/ref/ref.fasta -I /data/bam/reads_son.bam -O reads_son.g.vcf -L /data/ref/intervals.bed -ERC GVCF
17:30:10.017 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:30:10.156 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.159 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:30:10.159 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
17:30:10.159 INFO HaplotypeCaller - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:30:10.159 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:30:10.159 INFO HaplotypeCaller - Start Date/Time: February 11, 2026 at 5:30:09 PM GMT
17:30:10.159 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - ------------------------------------------------------------
17:30:10.160 INFO HaplotypeCaller - HTSJDK Version: 4.1.0
17:30:10.160 INFO HaplotypeCaller - Picard Version: 3.1.1
17:30:10.161 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:30:10.161 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:30:10.161 INFO HaplotypeCaller - Deflater: IntelDeflater
17:30:10.162 INFO HaplotypeCaller - Inflater: IntelInflater
17:30:10.162 INFO HaplotypeCaller - GCS max retries/reopens: 20
17:30:10.162 INFO HaplotypeCaller - Requester pays: disabled
17:30:10.162 INFO HaplotypeCaller - Initializing engine
17:30:10.277 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:30:10.290 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:30:10.296 INFO HaplotypeCaller - Done initializing engine
17:30:10.298 INFO HaplotypeCallerEngine - Tool is in reference confidence mode and the annotation, the following changes will be made to any specified annotations: 'StrandBiasBySample' will be enabled. 'ChromosomeCounts', 'FisherStrand', 'StrandOddsRatio' and 'QualByDepth' annotations have been disabled
17:30:10.302 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
17:30:10.303 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
17:30:10.304 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
17:30:10.307 INFO HaplotypeCallerEngine - Standard Emitting and Calling confidence set to -0.0 for reference-model confidence output
17:30:10.307 INFO HaplotypeCallerEngine - All sites annotated with PLs forced to true for reference-model confidence output
17:30:10.315 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
17:30:10.328 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
17:30:10.329 INFO IntelPairHmm - Available threads: 4
17:30:10.329 INFO IntelPairHmm - Requested threads: 4
17:30:10.329 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
17:30:10.368 INFO ProgressMeter - Starting traversal
17:30:10.369 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
17:30:10.875 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
17:30:11.980 INFO HaplotypeCaller - 14 read(s) filtered by: MappingQualityReadFilter
0 read(s) filtered by: MappingQualityAvailableReadFilter
0 read(s) filtered by: MappedReadFilter
0 read(s) filtered by: NotSecondaryAlignmentReadFilter
0 read(s) filtered by: NotDuplicateReadFilter
0 read(s) filtered by: PassesVendorQualityCheckReadFilter
0 read(s) filtered by: NonZeroReferenceLengthAlignmentReadFilter
0 read(s) filtered by: GoodCigarReadFilter
0 read(s) filtered by: WellformedReadFilter
14 total reads filtered out of 1981 reads processed
17:30:11.981 INFO ProgressMeter - 20_10037292_10066351:13223 0.0 35 1302.7
17:30:11.981 INFO ProgressMeter - Traversal complete. Processed 35 total regions in 0.0 minutes.
17:30:11.983 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0034843710000000004
17:30:11.983 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.048108363
17:30:11.983 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.02 sec
17:30:11.984 INFO HaplotypeCaller - Shutting down engine
[February 11, 2026 at 5:30:11 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=226492416
Uma vez concluído, você deve ter três arquivos terminando em .g.vcf no seu diretório atual (um por amostra) e seus respectivos arquivos de índice terminando em .g.vcf.idx.
Conteúdo do diretório
Neste ponto, chamamos variantes em modo GVCF para cada uma de nossas amostras de entrada. É hora de seguir para a chamada conjunta.
Mas não saia do contêiner! Vamos usar o mesmo na próxima etapa.
2.3. Executar genotipagem conjunta¶
Agora que temos todos os GVCFs, podemos experimentar a abordagem de genotipagem conjunta para gerar chamadas de variantes para uma coorte de amostras. É um método de duas etapas que consiste em combinar os dados de todos os GVCFs em um armazenamento de dados, depois executar a análise de genotipagem conjunta propriamente dita para gerar o VCF final de variantes chamadas conjuntamente.
2.3.1. Combinar todos os GVCFs por amostra¶
Esta primeira etapa usa outra ferramenta GATK, chamada GenomicsDBImport, para combinar os dados de todos os GVCFs em um armazenamento de dados GenomicsDB. O armazenamento de dados GenomicsDB é uma espécie de formato de banco de dados que serve como armazenamento intermediário para as informações de variantes.
gatk GenomicsDBImport \
-V reads_mother.g.vcf \
-V reads_father.g.vcf \
-V reads_son.g.vcf \
-L /data/ref/intervals.bed \
--genomicsdb-workspace-path family_trio_gdb
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar GenomicsDBImport -V reads_mother.g.vcf -V reads_father.g.vcf -V reads_son.g.vcf -L /data/ref/intervals.bed --genomicsdb-workspace-path family_trio_gdb
17:37:07.569 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:37:07.699 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.702 INFO GenomicsDBImport - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:37:07.702 INFO GenomicsDBImport - For support and documentation go to https://software.broadinstitute.org/gatk/
17:37:07.703 INFO GenomicsDBImport - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:37:07.703 INFO GenomicsDBImport - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:37:07.704 INFO GenomicsDBImport - Start Date/Time: February 11, 2026 at 5:37:07 PM GMT
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.704 INFO GenomicsDBImport - ------------------------------------------------------------
17:37:07.706 INFO GenomicsDBImport - HTSJDK Version: 4.1.0
17:37:07.706 INFO GenomicsDBImport - Picard Version: 3.1.1
17:37:07.707 INFO GenomicsDBImport - Built for Spark Version: 3.5.0
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:37:07.709 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:37:07.710 INFO GenomicsDBImport - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:37:07.710 INFO GenomicsDBImport - Deflater: IntelDeflater
17:37:07.711 INFO GenomicsDBImport - Inflater: IntelInflater
17:37:07.711 INFO GenomicsDBImport - GCS max retries/reopens: 20
17:37:07.711 INFO GenomicsDBImport - Requester pays: disabled
17:37:07.712 INFO GenomicsDBImport - Initializing engine
17:37:07.883 INFO FeatureManager - Using codec BEDCodec to read file file:///data/ref/intervals.bed
17:37:07.886 INFO IntervalArgumentCollection - Processing 6369 bp from intervals
17:37:07.889 INFO GenomicsDBImport - Done initializing engine
17:37:08.560 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:37:08.561 INFO GenomicsDBImport - Vid Map JSON file will be written to /tmp/family_trio_gdb/vidmap.json
17:37:08.561 INFO GenomicsDBImport - Callset Map JSON file will be written to /tmp/family_trio_gdb/callset.json
17:37:08.561 INFO GenomicsDBImport - Complete VCF Header will be written to /tmp/family_trio_gdb/vcfheader.vcf
17:37:08.561 INFO GenomicsDBImport - Importing to workspace - /tmp/family_trio_gdb
17:37:08.878 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.359 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.487 INFO GenomicsDBImport - Importing batch 1 with 3 samples
17:37:09.591 INFO GenomicsDBImport - Done importing batch 1/1
17:37:09.592 INFO GenomicsDBImport - Import completed!
17:37:09.592 INFO GenomicsDBImport - Shutting down engine
[February 11, 2026 at 5:37:09 PM GMT] org.broadinstitute.hellbender.tools.genomicsdb.GenomicsDBImport done. Elapsed time: 0.03 minutes.
Runtime.totalMemory()=113246208
Tool returned:
true
A saída desta etapa é efetivamente um diretório contendo um conjunto de diretórios ainda mais aninhados contendo os dados de variantes combinados na forma de múltiplos arquivos diferentes. Você pode explorar, mas rapidamente verá que este formato de armazenamento de dados não é destinado a ser lido diretamente por humanos.
Dica
O GATK inclui ferramentas que tornam possível inspecionar e extrair dados de chamada de variantes do armazenamento de dados conforme necessário.
2.3.2. Executar a análise de genotipagem conjunta propriamente dita¶
Esta segunda etapa usa ainda outra ferramenta GATK, chamada GenotypeGVCFs, para recalcular estatísticas de variantes e genótipos individuais à luz dos dados disponíveis em todas as amostras da coorte.
Saída do comando
Using GATK jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar GenotypeGVCFs -R /data/ref/ref.fasta -V gendb://family_trio_gdb -O family_trio.vcf
17:38:45.084 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/conda/share/gatk4-4.5.0.0-0/gatk-package-4.5.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:38:45.217 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.220 INFO GenotypeGVCFs - The Genome Analysis Toolkit (GATK) v4.5.0.0
17:38:45.220 INFO GenotypeGVCFs - For support and documentation go to https://software.broadinstitute.org/gatk/
17:38:45.220 INFO GenotypeGVCFs - Executing as root@be1a0302f6c7 on Linux v6.8.0-1030-azure amd64
17:38:45.220 INFO GenotypeGVCFs - Java runtime: OpenJDK 64-Bit Server VM v17.0.11-internal+0-adhoc..src
17:38:45.221 INFO GenotypeGVCFs - Start Date/Time: February 11, 2026 at 5:38:45 PM GMT
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - ------------------------------------------------------------
17:38:45.221 INFO GenotypeGVCFs - HTSJDK Version: 4.1.0
17:38:45.222 INFO GenotypeGVCFs - Picard Version: 3.1.1
17:38:45.222 INFO GenotypeGVCFs - Built for Spark Version: 3.5.0
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:38:45.222 INFO GenotypeGVCFs - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:38:45.223 INFO GenotypeGVCFs - Deflater: IntelDeflater
17:38:45.223 INFO GenotypeGVCFs - Inflater: IntelInflater
17:38:45.223 INFO GenotypeGVCFs - GCS max retries/reopens: 20
17:38:45.223 INFO GenotypeGVCFs - Requester pays: disabled
17:38:45.223 INFO GenotypeGVCFs - Initializing engine
17:38:45.544 INFO GenomicsDBLibLoader - GenomicsDB native library version : 1.5.1-84e800e
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field InbreedingCoeff - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAC - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.561 INFO NativeGenomicsDB - pid=221 tid=222 No valid combination operation found for INFO field MLEAF - the field will NOT be part of INFO fields in the generated VCF records
17:38:45.577 INFO GenotypeGVCFs - Done initializing engine
17:38:45.615 INFO ProgressMeter - Starting traversal
17:38:45.615 INFO ProgressMeter - Current Locus Elapsed Minutes Variants Processed Variants/Minute
17:38:45.903 WARN InbreedingCoeff - InbreedingCoeff will not be calculated at position 20_10037292_10066351:3480 and possibly subsequent; at least 10 samples must have called genotypes
GENOMICSDB_TIMER,GenomicsDB iterator next() timer,Wall-clock time(s),0.07757032800000006,Cpu time(s),0.07253379200000037
17:38:46.421 INFO ProgressMeter - 20_10037292_10066351:13953 0.0 3390 252357.3
17:38:46.422 INFO ProgressMeter - Traversal complete. Processed 3390 total variants in 0.0 minutes.
17:38:46.423 INFO GenotypeGVCFs - Shutting down engine
[February 11, 2026 at 5:38:46 PM GMT] org.broadinstitute.hellbender.tools.walkers.GenotypeGVCFs done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=203423744
Isso cria o arquivo de saída VCF family_trio.vcf no diretório de trabalho atual no contêiner, assim como seu índice, family_trio.vcf.idx.
É outro arquivo razoavelmente pequeno, então você pode executar cat family_trio.vcf para visualizar o conteúdo do arquivo, e rolar para baixo para encontrar as primeiras linhas de variantes.
Conteúdo do arquivo (resumido)
Destacamos mais uma vez a última linha do cabeçalho, que marca o início dos dados de chamada de variantes.
Isso parece similar ao VCF que geramos anteriormente, exceto que desta vez temos informações em nível de genótipo para todas as três amostras. As últimas três colunas no arquivo são os blocos de genótipo para as amostras, listadas em ordem alfabética de seu campo ID, como mostrado na linha de cabeçalho destacada.
Se olharmos os genótipos chamados para nosso trio familiar de teste para a primeira variante, vemos que o pai é heterozigoto-variante (0/1), e a mãe e o filho são ambos homozigotos-variante (1/1).
Essa é, em última análise, a informação que estamos procurando extrair do conjunto de dados!
2.3.3. Mover os arquivos de saída¶
Como observado anteriormente, qualquer coisa que permaneça dentro do contêiner ficará inacessível para trabalhos futuros. Antes de sairmos do contêiner, vamos mover os arquivos GVCF, o VCF final multi-amostra e todos os seus arquivos de índice manualmente para o sistema de arquivos fora do contêiner. Dessa forma, teremos algo para comparar quando construirmos nosso fluxo de trabalho para automatizar todo esse trabalho.
Conteúdo do diretório" hl_lines="14-19 22-23
data
├── bam
│ ├── reads_father.bam
│ ├── reads_father.bam.bai
│ ├── reads_mother.bam
│ ├── reads_mother.bam.bai
│ ├── reads_son.bam
│ └── reads_son.bam.bai
├── ref
│ ├── intervals.bed
│ ├── ref.dict
│ ├── ref.fasta
│ └── ref.fasta.fai
├── samplesheet.csv
└── vcf
├── family_trio.vcf
├── family_trio.vcf.idx
├── reads_father.g.vcf
├── reads_father.g.vcf.idx
├── reads_mother.g.vcf
├── reads_mother.g.vcf.idx
├── reads_mother.vcf
├── reads_mother.vcf.idx
├── reads_son.g.vcf
└── reads_son.g.vcf.idx
Uma vez feito isso, todos os arquivos estão agora acessíveis no seu sistema de arquivos normal.
2.3.4. Sair do contêiner GATK¶
Para sair do contêiner, digite exit.
Seu prompt deve voltar ao normal. Isso conclui o teste manual dos comandos de chamada conjunta de variantes.
Conclusão¶
Você sabe como testar os comandos de indexação Samtools e chamada de variantes GATK em seus respectivos contêineres, incluindo como gerar GVCFs e executar genotipagem conjunta em múltiplas amostras.
O que vem a seguir?¶
Faça uma pausa, depois siga para a Parte 2 para aprender como encapsular esses mesmos comandos em fluxos de trabalho que usam contêineres para executar o trabalho.