Uma comparação de condutas de chamada de variantes usando o genoma numa garrafa como referência
Abstract
Sequência de alto rendimento, especialmente de exomas, é uma ferramenta de diagnóstico popular, mas é difícil determinar quais são as melhores ferramentas para analisar estes dados. Neste estudo, nós usamos o Genoma NIST em uma garrafa como um novo recurso para validação do nosso pipeline de análise de exomas. Utilizamos seis alinhadores diferentes e cinco chamadores de variantes diferentes para determinar qual pipeline, do total de 30, executa o melhor num exoma humano que foi utilizado para ajudar a gerar a lista de variantes detectadas pelo Genoma em um Consórcio de Garrafas. Desses 30 pipelines, descobrimos que o Novoalign em conjunto com o GATK UnifiedGenotyper demonstrou a maior sensibilidade, mantendo um baixo número de falsos positivos para SNVs. No entanto, é aparente que os indels ainda são difíceis de manusear para qualquer gasoduto sem que nenhuma das ferramentas atinja uma sensibilidade média superior a 33% ou um Valor Preditivo Positivo (VPP) superior a 53%. Finalmente, como esperado, verificou-se que os alinhadores podem desempenhar um papel tão vital na detecção de variantes como os próprios chamadores de variantes.
1. Antecedentes
Nos últimos anos tem havido muitos avanços nas tecnologias de sequenciamento de alto rendimento. Devido a esses avanços, agora é possível detectar um grande número de variantes potencialmente causadoras de doenças e, em alguns casos, os dados de sequenciamento da próxima geração (NGS) têm até sido utilizados para fins de diagnóstico. Isto deve-se em parte aos desenvolvimentos nas tecnologias de sequenciamento nos últimos anos, mas também ao número de melhorias feitas nas várias ferramentas bioinformáticas utilizadas para analisar as montanhas de dados produzidos pelos instrumentos NGS .
Na busca de mutações em um paciente, um fluxo de trabalho típico é sequenciar seu exoma com um sequenciador Illumina, alinhar os dados brutos com o genoma de referência humano e, em seguida, identificar variantes de nucleotídeos simples (SNVs) ou inserções e deleções curtas (indels) que poderiam possivelmente causar ou influenciar o fenótipo de interesse . Embora isto seja bastante simples, decidir sobre as melhores ferramentas a utilizar em cada etapa do pipeline de análise não é. Há um grande número de ferramentas que são utilizadas em várias etapas intermediárias, mas as duas etapas mais importantes em todo o processo são o alinhamento das leituras brutas com o genoma e a busca de variantes (ou seja, SNVs e indels). Neste estudo, o nosso objetivo é ajudar o bioinformático de hoje elucidando a combinação correta de ferramenta de alinhamento de leitura curta e ferramenta de chamada de variante para o processamento de dados de seqüenciamento de exomas produzidos por instrumentos NGS.
Um número destes estudos foi realizado no passado, mas todos eles tiveram inconvenientes de alguma forma ou de outra. O ideal seria ter uma lista de todas as variantes conhecidas contidas em uma amostra para que quando um pipeline de ferramentas de análise for executado, você possa testá-lo para saber com certeza que está executando corretamente. No entanto, no passado não existia tal lista, por isso a validação tinha de ser feita por métodos menos completos. Em alguns casos, a validação era realizada pela geração de dados simulados de forma a criar um conjunto de verdadeiros positivos (TP) e verdadeiros negativos (TN) conhecidos. Enquanto isso convenientemente fornece uma lista de cada TP e TN no conjunto de dados, ele faz um trabalho pobre de representação precisa da biologia. Outros métodos de validação de variantes chamando pipelines incluem o uso de matrizes de genotipagem ou sequenciamento Sanger para obter uma lista de TPs e falsos positivos (FP) . Estes têm o lado positivo de fornecer resultados biologicamente validados, mas também têm o lado negativo de não serem abrangentes devido ao número limitado de pontos nas matrizes de genotipagem e ao custo proibitivo da validação de Sanger quando realizada milhares de vezes. Por último, nenhum destes estudos visou analisar o efeito que o alinhador de leitura curta teve na chamada de variantes. Consequentemente, o efeito a montante do desempenho do alinhador não pôde ser avaliado independentemente.
Neste estudo, temos a vantagem de uma lista de variantes para uma mulher anônima de Utah (ID do sujeito: NA12878, originalmente sequenciada para o projeto 1000 Genomes ) que foi validada experimentalmente pelo Consórcio Genoma em Garrafa (GiaB) liderado pelo NIST. Esta lista de variantes foi criada pela integração de 14 conjuntos de dados diferentes de cinco sequenciadores diferentes, e permite-nos validar qualquer lista de variantes geradas pelos nossos gasodutos de análise de exomas . A novidade deste trabalho é validar a combinação certa de alinhadores e chamadores de variantes em relação a um conjunto de dados de variantes abrangente e determinado experimentalmente: NIST-GiaB.
Para realizar nossa análise estaremos usando um dos conjuntos de dados de exoma originalmente usados para criar a lista NIST-GiaB. Escolhemos apenas um dos exomos gerados pelo Illumina TruSeq original porque queríamos fornecer um cenário de caso de uso padrão para alguém que deseja realizar a análise NGS, e enquanto todo o sequenciamento do genoma continua a cair no preço, o sequenciamento de exomos ainda é uma alternativa popular e viável . É também importante notar que, por Bamshad et al., atualmente o número esperado de SNVs por exoma europeu-americano é de 20.283 ± 523 . Apesar disso, o número total de SNVs encontrados na lista NIST-GiaB com potencial para existir no conjunto de dados do exome TruSeq foi de 34.886, o que é significativamente maior do que o esperado. Isto é provavelmente devido ao fato de que, enquanto o kit exome foi usado para gerar dados NIST-GiaB, ele também foi suplementado por sequenciamento de genoma inteiro.
Por último, nós consideramos um grande número de alinhadores e chamadores de variantes, mas finalmente escolhemos as 11 ferramentas baseadas na prevalência, popularidade e relevância para o nosso conjunto de dados (por exemplo, SNVMix, VarScan2 e MuTect não foram usados, pois eles são destinados para uso em amostras derivadas de tumores). Nossa própria análise envolve a comparação de seis alinhadores (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK , e Novoalign) e cinco chamadores de variantes (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup , e SNPSVM ). Neste estudo também tentamos determinar quanto de um efeito, se houver, o alinhador tem sobre a chamada de variantes e quais alinhadores têm melhor desempenho quando se utiliza uma amostra normal do exome Illumina. Para nosso conhecimento, este é o primeiro relatório que valida todas as combinações possíveis (total de 30 pipelines) de uma ampla gama de alinhadores e chamadores de variantes.
2. Métodos
2.1. Datasets
Human reference genome hg19 foi baixado do navegador da UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) e foi usado para realizar os alinhamentos. O exoma humano, SRR098401, foi baixado do Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Para fins de anotação e calibração, foram usadas as listas dbSNP137 sem sites após a versão 129, HapMap 3.3, Human Omni 2.5 BeadChip, e Mills e 1000 G gold standard indel set (todas de ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).
2.2. O Pipeline
Figure 1 mostra o fluxo de trabalho utilizado neste estudo, que é semelhante ao descrito no guia de Boas Práticas produzido pelo The Broad Institute . Isto envolve uma série de passos para garantir que os arquivos de alinhamento produzidos sejam da mais alta qualidade, bem como vários outros para garantir que as variantes sejam chamadas corretamente. Primeiro, as leituras em bruto foram alinhadas com hg19, e depois as duplicatas de PCR foram removidas do alinhamento. Em seguida, para ajudar na identificação do indel mais tarde no pipeline, a leitura do realinhamento foi feita em torno dos indels. O último passo do processamento do alinhamento foi realizar uma recalibração da pontuação de qualidade de base, o que ajuda a melhorar o viés inerente e as imprecisões das pontuações emitidas pelos sequenciadores. Infelizmente, apesar desses passos, a taxa de alinhamento de cada alinhador foi significativamente menor do que o esperado, então para compensar isso, o kit de ferramentas fastx foi usado para filtrar as leituras de baixa qualidade (Tabela 1). As leituras de baixa qualidade foram definidas como aquelas que tinham pelo menos metade de suas pontuações de qualidade abaixo de 30. Após o processamento do alinhamento, a chamada de variantes e a filtragem de variantes foram realizadas.
|
||||||||||||||||||||||||||||||||||||||||
As seis ferramentas usadas para gerar alinhamentos foram Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK e Novoalign, e as cinco ferramentas usadas para gerar variantes foram FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup, e SNPSVM, como pode ser visto na Tabela 2.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2.3. Filtragem
Os dados em bruto foram adquiridos do SRA (SRR098401), divididos com fastq-dump, e filtrados usando o kit de ferramentas fastx. Especificamente, fastq-dump usou as bandeiras -split_files e -split_spot, e fastq_quality_filter foi executado com os seguintes argumentos: -Q 33 -q 30 -p 50. Depois as leituras foram devidamente emparelhadas com um script personalizado.
2.4. Alinhando
Alinhadores usaram argumentos padrão, exceto quando um argumento de threads foi usado, quando disponível. Os comandos utilizados são os seguintes.
2.4.1. Bowtie2
(1)bowtie2 -p 10 -x $INDEX -1 raw_data/read_1_filtered.fastq -2 raw_data/read_2_filtered.fastq -alinhamentos S/NA12878.bt2.sam
2.4.2. BWA Sampe
(1)bwa aln -t 10 genome/hg19.fa raw_data/read_1_filtered.fastq > alignments/NA12878.R1.sai(2)bwa aln -t 10 genome/hg19.fa raw_data/read_2_filtered.fastq > alignments/NA12878.R2.sai(3)bwa sampe genome/hg19.fa alignments/NA12878.R1.sai alignments/NA12878.R2.sai raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > alignments/NA12878.bwa-sampe.sam
2.4.3. BWA Mem
(1)bwa mem -t 10 genome/hg19.fa raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq > alignments/NA12878.bwa-mem.sam
2.4.4. CUSHAW3
(1)cushaw3 align -r $INDEX -t 10 -o alignments/NA12878.CUSHAW3.sam -q raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq
2.4.5. MOSAIK
(1)MosaikBuild -q raw_data/read_1_filtered.fastq -q2 raw_data/read_2_filtered.fastq -st illumina -outalignments/NA12878.MOSAIK.mkb(2)MosaikAligner -em alinhamentos/NA12878.MOSAIK.mkb – sem alinhamentos/NA12878.MOSAIK -p 10 -ia genome/hg19.dat -j genome/hg19_15 -annpe tools/MOSAIK/src/networkFile/2.1.78.pe.ann -annse tools/MOSAIK/src/networkFile/2.1.78.se.ann
2.4.6. Novoalign
(1)novoalign -d $INDEX -f raw_data/read_1_filtered.fastq raw_data/read_2_filtered.fastq -o SAM -c 10 > alignments/NA12878.novoalign.sam
2.5. Alignment Proth of Coverage Calculation
Para garantir o cálculo adequado da profundidade de cobertura, o módulo CalculateHsMetrics do Picard Tools foi utilizado com os seguintes argumentos:(1)java -jar CalculateHsMetrics.jar I=NA12878.ALN.BQSR.bam O=ALN.O.log R=genome/hg19.fa TI=genome/truseq_exome.bed BI=genome/truseq_exome.bed VALIDATION_STRINGENCY=SILENT PER_TARGET_COVERAGE=ALN.ptc.bed É importante notar que o arquivo TruSeq exome bed deve ter o cabeçalho do arquivo de alinhamento SAM preparado para ele para que este módulo funcione. Além disso, a coluna 6 deve ser movida para a coluna 4, e a coluna 5 precisa ser removida do arquivo de leito TruSeq.
2.6. Processamento do arquivo de alinhamento
Processamento dos arquivos de alinhamento (arquivos SAM/BAM) requer os seguintes passos para todos os alinhadores:(1)Conversão de SAM para BAM com visualização SAMtools:(a)visualização samtools -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)BAM classificação de arquivos usando o módulo Picard Tools, SortSam:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam SORT_ORDER=coordenado(3)Remoção de duplicados PCR usando o módulo Picard Tools, MarkDuplicates:(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.sorted.bam O=alignments/NA12878.ALN.dups_removed.bam REMOVE_DUPLICATES=true M=alignments/metrics(4)Read Group added to alignment files using the Picard Tools module, AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=coordenado RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Realinhamento em torno de indels usando os módulos GATK RealignerTargetCreator e IndelRealigner:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T RealignerTargetCreator -R genoma/hg19.fa -I alignments/NA12878.ALN.RG.bam – genoma/mills.vcf conhecido -o tmp/ALN.intervals(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T IndelRealigner -R genoma/hg19.fa -I alignments/NA12878.ALN.RG.bam – genoma/mills.vcf -o alignments/NA12878.ALN.indels.bam -maxReadsForRealignment 100000 -maxReadsInMemory 1000000 -targetIntervals tmp/ALN.intervals(6)Recalibration Base using the GATK modules BaseRecalibrator and PrintReads:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T BaseRecalibrator -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -knownSitesgenomedbsnp_137.hg19.excluding_sites_after_129.only_standard_chroms.vcf -o tmp/NA12878.ALN.grp(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T PrintReads -R genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -BQSR tmp/NA12878.ALN.grp -o alignments/NA12878.ALN.BQSR.bam
2.7. Chamada de variante
Default arguments were used for each variant caller unless it contained a “threads” or “parallel” flag in which case that was used as well as. Adicionalmente, os indels foram chamados separadamente das SNVs sempre que possível. Especificamente, os comandos utilizados são os seguintes.
2.7.1. FreeBayes
(1)freebayes -f genome/hg19.fa -i -X -u -v vcf_files/NA12878.ALIGNER.freebayes.raw.snv.vcf alignments/NA12878.ALIGNER.BQSR.bam(2)freebayes -f genome/hg19.fa -I -X -u -v vcf_files/NA12878.ALIGNER.freebayes.raw.indel.vcf alignments/NA12878.ALIGNER.BQSR.bam
2.7.2. GATK HaplotypeCaller
(1)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R genome/hg19.fa -I alignments/NA12878.ALIGNER.BQSR.bam -dbsnp $DBSNP -o vcf_files/NA12878.ALIGNER.HC.raw.vcf -stand_call_conf 50
2.7.3. GATK UnifiedGenotyper
(1)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R genome/hg19.fa -nt 10 -I alignments/NA12878.ALIGNER.BQSR.bam -o vcf_files/NA12878.ALIGNER.UG.raw.snv.vcf -glm SNP -D $DBSNP(2)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T UnifiedGenotyper -R genome/hg19.fa -nt 10 -I alignments/NA12878.ALIGNER.BQSR.bam -o vcf_files/NA12878.ALIGNER.UG.raw.indel.vcf -glm INDEL -D $MILLS
2.7.4. SAMtools Mpileup
(1)samtools mpileup -uf genome/hg19.faalignments/NA12878.ALIGNER.BQSR.bam | bcftools view -bvcg – > vcf_files/NA12878.ALIGNER.mpileup.bcf &753> bcftools view vcf_files/NA12878.ALIGNER.mpileup.bcf > vcf_files/NA12878.ALIGNER.mpileup.raw.vcf
2.7.5. SNPSVM
(1)java -XX:ParallelGCThreads=10 -jar tools/SNPSVM/snpsvm.jar predict -R genome/hg19.fa -B alignments/NA12878.ALIGNER.BQSR.bam -M tools/SNPSVM/models/default.model -V vcf_files/NA12878.ALIGNER.SNPSVM.raw.vcfDue à inexistência de flags CIGAR necessários no arquivo de alinhamento, SNPSVM falhou em chamar variantes para CUSHAW3, e SAMtools mpileup não pôde chamar variantes em alinhamentos MOSAIK pelo mesmo motivo. Além disso, devido ao fato da SNPSVM detectar apenas SNVs, nenhum indels foi reportado para este programa.
2.8. Filtração de Variantes
Filtração variou dependendo da variante que está sendo utilizada. Nos casos do GATK HaplotypeCaller e GATK UnifiedGenotyper, os módulos GATK, VariantRecalibrator e ApplyRecalibration, foram usados para filtrar SNVs usando HapMap 3.3, o Omni 2.5 SNP BeadChip, e o dbSNP 137 sem dados de 1000 Genoma como conjuntos de treinamento. Para SNPSVM, foram utilizadas as pontuações QUAL ≥ 4 e os valores DP ≥ 6. Para FreeBayes e SAMtools, foram utilizados os valores de QUAL ≥ 20 e DP ≥ 6.
2.9. Comparação de Variantes
Para comparação de variantes, foi usado o USeq 8.8.1 para comparar SNVs compartilhados entre todos os conjuntos de dados. Para comparar indels, foi utilizada a ferramenta vcflib vcfintersect. O arquivo TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, obtido em 11 de dezembro de 2013, foi usado para restringir comparações a locais que poderiam ser capturados pelo kit de pull down exome usado no sequenciamento da SRR098401. Este arquivo pode ser obtido aqui da Illumina: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Para garantir que as variantes fossem representadas de forma idêntica entre diferentes conjuntos de chamadas, a ferramenta vcflib vcfallelicprimitives foi utilizada para pré-processar arquivos vcf.
2.10. Cálculos estatísticos
True Positive (TP). É uma mutação que foi detectada pelo pipeline a ser testado e que existe na lista NIST-GiaB.
False Positive (FP). É uma mutação que foi detectada pelo pipeline a ser testado, mas que não existe na lista NIST-GiaB.
Negativo Verdadeiro (TN). É uma mutação que não foi detectada pelo pipeline a ser testado e é uma mutação que não existe na lista NIST-GiaB.
Falso Negativo (FN). É uma mutação que não foi detectada pelo pipeline a ser testado mas que não existe na lista NIST-GiaB:
3. Resultados e Discussão
3.1. Pré-filtragem de Variantes
Ao realizar a análise de variantes, uma das muitas armadilhas que devem ser levadas em consideração é o espaço de sequência de exomas (conforme definido pelo kit de captura de exomas) e como isso pode afetar os resultados da análise. Neste caso, tivemos um único exome (SRR098401) que foi extraído usando o kit exome Illumina TruSeq e sequenciado em um HiSeq 2000. Com isto em mente, queríamos ter certeza de que estávamos medindo a capacidade das ferramentas bioinformáticas para fazer seu trabalho e não o quão bem o kit de captura do exome do Illumina TruSeq funcionava. Ou seja, só queremos tentar chamar variantes que devem estar presentes nos exons, como definido pelo kit pull down. Por este motivo, usamos o arquivo cama fornecido pelo Illumina, não um arquivo cama de anotação genérica, por exemplo, RefSeq para hg19. Verificamos que para este indivíduo em particular, de acordo com a lista NIST-GiaB, deveria haver um total de 34.886 SNVs e 1.473 indels dentro das regiões definidas pelo arquivo de leito TruSeq.
Após filtrarmos as variantes que não estavam localizadas nas regiões definidas pelo arquivo de leito exome do Illumina TruSeq, passamos de centenas de milhares de variantes putativas (dados não mostrados) para, em média, cerca de 23.000 variantes (SNVs e indels) por duto (Tabela 3). Este é um passo importante para os investigadores, uma vez que reduz significativamente o espaço de pesquisa de variantes potencialmente interessantes.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3.2. Resultados das Variantes Brutas Comparadas com GiaB
Um aspecto que quisemos compreender ao fazer esta comparação foi a importância de filtrar as variantes detectadas por estas ferramentas. A razão para isso é que, idealmente, gostaríamos de ter um nível de sensibilidade tão alto quanto possível para que as mutações de interesse não se percam no processo de filtragem. Portanto, cabe a nós determinar se esta etapa é ou não necessária e até que ponto ela é necessária, já que está claro a partir dos resultados do NIST-GiaB e da revisão de Bamshad et al. que a sensibilidade pode ser um problema.
Como podemos ver na Tabela 3, a filtragem é necessária mais para alguns chamadores de variantes do que para outros quando se trata de SNVs, e é absolutamente necessária para os indels. Na maioria dos casos, o número de variantes TP é próximo ou superior ao nosso número esperado de cerca de 20.000 , mas, por outro lado, em alguns casos, o número de FPs é muito alto.
Claramente há muita variação nos números gerados por cada gasoduto. No entanto, é possível encontrar alguns pontos em comum nos números que provavelmente provêm da origem algorítmica de cada ferramenta. O FreeBayes produz tanto o maior número de variantes não filtradas como o maior número de FPs. É provável que só vejamos este tipo de desempenho desta ferramenta devido ao fato de que, embora não seja a única variante que chama com base na inferência Bayesiana, ela é única na sua interpretação de alinhamentos. Ou seja, é um chamador baseado em haplótipos que identifica variantes com base na sequência das próprias leituras em vez do alinhamento, o último dos quais é como funciona o UnifiedGenotyper do GATK.
Adicionalmente vemos os alinhadores baseados em Burrows-Wheeler terem um desempenho muito semelhante um ao outro: Bowtie2, BWA mem, e BWA sampe alcançam resultados semelhantes em todos os aspectos. Pode-se supor que isso se deve provavelmente ao fato de que todas essas ferramentas utilizam algoritmos similares ao executar sua tarefa designada. Esta observação é suportada pelo fato de que MOSAIK (alinhamentos espaçados usando o algoritmo Smith-Waterman) e CUSHAW3 (uma abordagem de semeadura híbrida) têm ambos algoritmos subjacentes muito diferentes e subsequentemente produzem resultados muito diferentes.
Esta diferença nos resultados correlacionados com algoritmos diferentes é melhor vista nos resultados SNPSVM. Dos chamadores de variantes, é o único que utiliza máquinas vetoriais de suporte e construção de modelos para gerar chamadas SNV. Parece que embora tenha a desvantagem de não ser tão sensível quanto outros métodos, ela se beneficia de ser extremamente precisa, independentemente do alinhador utilizado. Isto sugere que se é capaz de saltar completamente o passo de filtragem ao usar esta variante de chamada.
No que diz respeito aos indels, nenhum alinhador parece destacar-se entre os demais como aquele que lida bem com este tipo de mutação. Na verdade, ao olhar para o número de FPs, é claro que é o chamador da variante que desempenha o maior papel na precisão da identificação do indel. Adicionalmente, não existem dados para os pipelines CUSHAW3 mais SNPSVM ou MOSAIK mais SAMtools mpileup devido aos ficheiros de alinhamento que não contêm as cordas CIGAR necessárias para que os chamadores de variantes possam funcionar a jusante. Finalmente, a razão pela qual não existem dados indel para SNPSVM é porque esta ferramenta é utilizada apenas para identificação de SNVs.
3.3. Resultados Filtrados Comparado com GiaB
Como na Tabela 4, as práticas de filtragem padrão conseguem remover um grande número de SNVs FP para cada gasoduto; no entanto, parece que estes filtros são um pouco agressivos demais na maioria dos casos para a detecção de SNV, mas não o suficiente para os indels. Isso se torna óbvio quando se observam as diferenças no número de FPs reportadas em cada conjunto de dados. Por exemplo, Bowtie2 com Freebayes vê uma remoção de 72.570 FP SNVs (uma redução de 98%), mas apenas uma remoção de 1.736 FP indels (uma redução de 70%). Também deve ser notado que os filtros utilizados foram dependentes da tubulação e, na sua maioria, dentro de cada tubulação produziu reduções semelhantes em SNV e FP indel. A única exceção aqui é o número de variantes identificadas a partir dos alinhamentos CUSHAW3, quando comparado com outros alinhamentos: Em geral, o número de SNV TP é menor, o número de SNV FP é maior e é o único alinhador que produz mais de 1.000 indels de FP após a filtragem.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Posto que a filtragem reduz significativamente o número de variantes TP, pode ser sensato, com excepção dos pipelines que utilizam CUSHAW3 e FreeBayes, saltar este passo ao procurar variantes raras e de alto impacto. Ao invés disso, pode-se gastar mais tempo em um processo de filtragem baseado em biologia do que em estatísticas. Por exemplo, pode fazer mais sentido investir tempo na identificação de uma pequena lista de variantes que provavelmente serão de alto impacto: emendas mutações no local, indels que causam emendas, mutações de truncagem, mutações stop-loss, ou mutações em genes conhecidos por serem biologicamente relevantes para o fenótipo de interesse.
3.4. TPR médio e Sensibilidade
Como pode ser visto na Tabela 5, o Valor Predictivo Positivo (VPP) para cada ferramenta, com exceção de CUSHAW3, varia de 91% a 99,9% para SNVs, mas a sensibilidade média é muito baixa (cerca de 50%). Esta discrepância pode ser devida a uma série de razões, mas a mais provável é que a profundidade de cobertura seja variável entre os exons. Podemos ver que, além da baixa sensibilidade do VNS, a sensibilidade do indel é baixa (cerca de 30%); contudo o VPP para os indels é consideravelmente mais baixo (35,86% a 52,95%). Isto pode ser devido a qualquer uma das seguintes razões: os indels muito curtos são difíceis de detectar pelo NGS convencional , a representação dos indels por diferentes chamadores de variantes pode fazer com que as ferramentas afirmem incorrectamente que dois indels são diferentes, ou ferramentas de alinhamento produzem representações diferentes do mesmo indel .
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Talvez a explicação mais provável para ambos os tipos de mutações seja a questão da profundidade. Como em qualquer estudo de análise de variantes, um aumento na profundidade de cobertura leva a um aumento na sensibilidade, mas é impossível garantir uma boa profundidade de cobertura devido à incapacidade dos kits de captura de exome de puxar para baixo os exons de forma uniforme. Além disso, nenhum kit de captura de exônomos cobre cada exônomo. De facto, foi demonstrado que a análise de variantes de sequenciamento de todo o genoma a uma profundidade média inferior a um exoma tem um melhor desempenho devido à uniformidade da referida profundidade. Assim, é provável que um grande número de variantes esteja faltando devido ao fato de que a lista NIST-GiaB foi criada a partir de uma compilação de dados de seqüenciamento exonômico e genômico. Em última análise, para se conseguir uma sensibilidade adequada, será necessário realizar o sequenciamento genômico completo, mas isso é atualmente proibitivo em termos de custo para a maioria dos laboratórios. Felizmente, este custo continua a cair, e em breve veremos uma mudança gradual da análise de exomas para a mais completa análise do genoma inteiro.
3.5. Sensibilidade como função da profundidade
Porque a sensibilidade reflete uma das métricas de desempenho mais importantes de uma ferramenta e a maioria das ferramentas luta para atingir sensibilidade superior a 50%, gostaríamos de explorar mais a profundidade que a variante afetada chama sensibilidade. Nós examinamos uma série de diferentes combinações de ferramentas para determinar quais eram os melhores pipelines, chamadores de variantes e alinhadores. Para a Figura 2, tomamos as cinco melhores combinações de chamadores e alinhadores de variantes conforme determinado pela sua sensibilidade e taxa de falsos positivos (FPR). Ou seja, seleccionámos as que tinham o maior número de SNV TP chamadas, para além do menor número de SNV FP. Na inspeção, o que se destaca imediatamente é que a sensibilidade é menor do que o esperado. Todos os gasodutos têm um desempenho aproximadamente ao mesmo nível: identificam a maioria das suas variantes quando é atingida uma profundidade de cerca de 150x, o que indica que esta profundidade é provavelmente suficiente e que o número de variantes em falta é provavelmente devido ao facto de certos exons terem uma cobertura inferior à média. Note que quatro dos cinco pipelines com melhor desempenho têm o GATK UnifiedGenotyper como seu chamador de variantes, demonstrando seu desempenho superior independentemente do alinhador utilizado como mostrado na Figura 3(b).
(a)
(b)
(a)
(b)
Além de olharmos para os cinco alinhadores superiores, determinamos que seria útil realizar a mesma análise no melhor alinhador acoplado com cada chamador de variante (Figura 3(a)), bem como o melhor chamador de variante acoplado com cada alinhador (Figura 3(b)). Como com os dutos, o melhor alinhador foi identificado como o que produziu o maior número de SNV TP e o menor número de SNV FP – neste caso, BWA mem. Apesar de ter o melhor alinhamento para trabalhar, ainda vemos uma diferença bastante grande entre os chamadores da variante, que é provavelmente atribuível aos diferentes algoritmos que eles empregam (Figura 3(a)). No entanto, no caso do chamador da variante com melhor desempenho, GATK UnifiedGenotyper, parece haver menos variação entre os quatro alinhadores superiores indicando que ele tem um desempenho bastante bom na maioria das situações, sendo as exceções CUSAHW3 e MOSAIK.
3.6. Variantes Compartilhadas entre os Topo de Dutos
Por último, queríamos saber quão únicos eram os conjuntos de chamadas de variantes entre os diferentes dutos. Para fazer isto, concentrámo-nos novamente nas cinco principais variantes que chamam pipelines: Bowtie2 mais UnifiedGenotyper, BWA mem mais UnifiedGenotyper, BWA sampe mais HaplotypeCaller, BWA sampe mais UnifiedGenotyper, e Novoalign mais UnifiedGenotyper. Como pode ser visto na Figura 4, existe uma grande quantidade de sobreposição entre os cinco diferentes pipelines em questão, com 15.489 SNVs (70%) compartilhados de um total de 22.324 variantes distintas. No entanto, também se pode argumentar que isto se deve em grande parte a quatro dos cinco gasodutos que utilizam o UnifiedGenotyper como seu chamador de variantes. Esta noção é corroborada pelo fato de que o maior número de variantes exclusivas de um gasoduto, 367, pertence à combinação BWA sampe mais HaplotypeCaller. Também vale a pena notar que o segundo maior número de SNVs exclusivos também pertence ao alinhador de SNVs do BWA, portanto é possível que o alto número de SNVs exclusivos seja melhor atribuído ao alinhador do que ao chamador de variantes.
4. Conclusões
Verificámos que entre os trinta diferentes pipelines testados Novoalign mais GATK UnifiedGenotyper demonstraram a maior sensibilidade mantendo um número baixo de FP para SNV. Entre os alinhadores, o BWA mem apresentou consistentemente o melhor desempenho, mas os resultados ainda variaram muito dependendo da variante de chamada utilizada. Naturalmente, segue-se que o melhor chamador de variantes, o GATK UnifiedGenotyper, na maioria das vezes produziu resultados semelhantes, independentemente do alinhador utilizado. No entanto, é facilmente aparente que os indels ainda são difíceis de manusear para qualquer gasoduto com nenhum dos gasodutos atingindo uma sensibilidade média superior a 33% ou um PPV superior a 53%. Além do baixo desempenho geral que vemos na detecção de indels, a sensibilidade, independentemente do tipo de mutação, é um problema para cada duto descrito neste trabalho. O número esperado de SNV para o exoma de NA12878 é de 34.886, mas mesmo quando se utiliza a união de todas as variantes identificadas pelos cinco maiores gasodutos, o maior número identificado foi muito baixo (22.324). Parece que embora ainda muito útil a análise do exome tem suas limitações mesmo quando se trata de algo aparentemente tão simples como a detecção de SNV.
Disclosure
Adam Cornish é um estudante de pós-graduação no laboratório de Chittibabu Guda com treinamento em ciência da computação e genômica. Chittibabu Guda (Professor Associado) tem uma formação interdisciplinar em biologia molecular e computacional. Ele tem publicado uma série de métodos computacionais com uma variedade de aplicações em pesquisa biomédica, desde 2001.
Conflito de Interesses
Os autores desconhecem quaisquer interesses concorrentes.
Contribuição dos autores
Adam Cornish desenhou o estudo, realizou todas as análises, fez as figuras, e escreveu o artigo. Chittibabu Guda forneceu o feedback essencial sobre as melhorias do trabalho e a contribuição das próprias análises e editou completamente o trabalho. Todos os autores leram e aprovaram o trabalho final.
Agradecimentos
Este trabalho foi apoiado pelos fundos de desenvolvimento para Chittibabu Guda do Centro Médico da Universidade de Nebraska (UNMC). Os autores gostariam de agradecer aos criadores do Novoalign por disponibilizar o software como versão de teste e ao Núcleo de Bioinformática e Biologia de Sistemas da UNMC pelo suporte de infra-estrutura que facilitou esta pesquisa.