Abstract

Die Hochdurchsatz-Sequenzierung, insbesondere von Exomen, ist ein beliebtes Diagnoseinstrument, aber es ist schwierig zu bestimmen, welche Instrumente am besten für die Analyse dieser Daten geeignet sind. In dieser Studie nutzen wir die Ergebnisse des NIST Genome in a Bottle als neue Ressource für die Validierung unserer Exom-Analyse-Pipeline. Wir verwenden sechs verschiedene Aligner und fünf verschiedene Varianten-Caller, um festzustellen, welche der insgesamt 30 Pipelines am besten auf einem menschlichen Exom funktioniert, das zur Erstellung der Liste der vom Genome in a Bottle Consortium entdeckten Varianten verwendet wurde. Von diesen 30 Pipelines haben wir festgestellt, dass Novoalign in Verbindung mit GATK UnifiedGenotyper die höchste Sensitivität aufweist und gleichzeitig eine niedrige Anzahl von falsch-positiven SNVs aufweist. Es ist jedoch offensichtlich, dass Indels nach wie vor für jede Pipeline schwierig zu handhaben sind, da keines der Tools eine durchschnittliche Sensitivität von mehr als 33 % oder einen positiven prädiktiven Wert (PPV) von mehr als 53 % erreichte. Schließlich wurde wie erwartet festgestellt, dass Aligner eine ebenso wichtige Rolle bei der Variantenerkennung spielen können wie Varianten-Caller selbst.

1. Hintergrund

In den letzten Jahren wurden zahlreiche Fortschritte bei den Hochdurchsatz-Sequenzierungstechnologien erzielt. Dank dieser Fortschritte ist es nun möglich, eine große Anzahl potenziell krankheitsverursachender Varianten zu erkennen, und in einigen wenigen Fällen wurden Next Generation Sequencing (NGS)-Daten sogar für diagnostische Zwecke verwendet. Dies ist zum Teil auf die Entwicklungen bei den Sequenzierungstechnologien in den letzten Jahren zurückzuführen, aber auch auf die zahlreichen Verbesserungen bei den verschiedenen bioinformatischen Werkzeugen, die zur Analyse der von den NGS-Geräten erzeugten Datenmengen verwendet werden.

Bei der Suche nach Mutationen bei einem Patienten besteht ein typischer Arbeitsablauf darin, das Exom mit einem Illumina-Sequenzer zu sequenzieren, die Rohdaten am menschlichen Referenzgenom auszurichten und dann einzelne Nukleotidvarianten (SNVs) oder kurze Insertionen und Deletionen (Indels) zu identifizieren, die möglicherweise den gewünschten Phänotyp verursachen oder beeinflussen könnten. Während dies relativ einfach ist, ist die Entscheidung über die besten Werkzeuge für die einzelnen Phasen der Analyse-Pipeline nicht einfach. Es gibt eine Vielzahl von Tools, die in verschiedenen Zwischenschritten eingesetzt werden, aber die beiden wichtigsten Schritte im gesamten Prozess sind das Alignment der Rohdaten zum Genom und die anschließende Suche nach Varianten (d. h. SNVs und Indels). In dieser Studie wollen wir den Bioinformatikern von heute helfen, indem wir die richtige Kombination von Short-Read-Alignment-Tool und Variantenaufruf-Tool für die Verarbeitung von Exom-Sequenzierungsdaten, die von NGS-Instrumenten erzeugt wurden, aufzeigen.

In der Vergangenheit wurde eine Reihe dieser Studien durchgeführt, die jedoch alle in der einen oder anderen Form Nachteile hatten. Idealerweise sollte man über eine Liste aller bekannten Varianten in einer Probe verfügen, so dass man bei der Ausführung einer Pipeline von Analysewerkzeugen diese testen kann, um mit Sicherheit zu wissen, dass sie korrekt arbeitet. In der Vergangenheit gab es jedoch keine solche Liste, so dass die Validierung mit weniger vollständigen Methoden durchgeführt werden musste. In einigen Fällen wurde die Validierung durch die Erzeugung simulierter Daten durchgeführt, um einen Satz bekannter wahrer Positiver (TP) und wahrer Negativer (TN) zu erstellen. Auf diese Weise erhält man zwar eine Liste aller TP und TN im Datensatz, aber die Biologie wird dadurch nur unzureichend wiedergegeben. Andere Methoden zur Validierung von Variantenaufruf-Pipelines umfassen die Verwendung von Genotypisierungs-Arrays oder Sanger-Sequenzierung, um eine Liste von TPs und Falsch-Positiven (FP) zu erhalten. Diese Methoden haben den Vorteil, dass sie biologisch validierte Ergebnisse liefern, haben aber auch den Nachteil, dass sie aufgrund der begrenzten Anzahl von Spots auf Genotypisierungs-Arrays und der unerschwinglichen Kosten der Sanger-Validierung, wenn diese Tausende von Malen durchgeführt wird, nicht umfassend sind. Schließlich wurde in keiner dieser Studien untersucht, wie sich der Short-Read-Aligner auf die Variantenbestimmung auswirkt. Folglich konnte der vorgelagerte Effekt der Alignerleistung nicht unabhängig bewertet werden.

In dieser Studie haben wir den Vorteil einer Variantenliste für eine anonyme Frau aus Utah (Subjekt-ID: NA12878, ursprünglich sequenziert für das 1000 Genomes Projekt), die experimentell durch das vom NIST geleitete Genome in a Bottle (GiaB) Konsortium validiert wurde. Diese Variantenliste wurde durch die Integration von 14 verschiedenen Datensätzen von fünf verschiedenen Sequenzierern erstellt und ermöglicht es uns, jede von unseren Exomanalyse-Pipelines erstellte Variantenliste zu validieren. Die Neuheit dieser Arbeit besteht darin, die richtige Kombination von Alignern und Varianten-Callern anhand eines umfassenden und experimentell ermittelten Variantendatensatzes zu validieren: NIST-GiaB.

Zur Durchführung unserer Analyse verwenden wir einen der Exom-Datensätze, die ursprünglich zur Erstellung der NIST-GiaB-Liste verwendet wurden. Wir haben uns für einen der ursprünglich von Illumina TruSeq generierten Exome entschieden, weil wir ein Standardanwendungsszenario für jemanden bieten wollten, der eine NGS-Analyse durchführen möchte, und obwohl die Preise für die Ganzgenomsequenzierung weiter sinken, ist die Exomsequenzierung immer noch eine beliebte und praktikable Alternative. Es ist auch wichtig zu wissen, dass laut Bamshad et al. die erwartete Anzahl von SNVs pro europäisch-amerikanischem Exom derzeit 20.283 ± 523 beträgt. Trotzdem betrug die Gesamtzahl der in der NIST-GiaB-Liste gefundenen SNVs, die potenziell im TruSeq-Exom-Datensatz vorhanden sind, 34 886, was deutlich höher ist als erwartet. Dies ist wahrscheinlich auf die Tatsache zurückzuführen, dass das Exom-Kit zwar für die Generierung der NIST-GiaB-Daten verwendet wurde, aber auch durch die Sequenzierung des gesamten Genoms ergänzt wurde.

Schließlich haben wir eine große Anzahl von Alignern und Varianten-Callern in Betracht gezogen, aber letztendlich die 11 Tools auf der Grundlage ihrer Verbreitung, Popularität und Relevanz für unseren Datensatz ausgewählt (z. B. wurden SNVMix, VarScan2 und MuTect nicht verwendet, da sie für die Verwendung bei von Tumoren stammenden Proben bestimmt sind). Unsere Analyse selbst umfasst den Vergleich von sechs Alignern (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK und Novoalign) und fünf Varianten-Callern (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup und SNPSVM ). In dieser Studie versuchen wir auch herauszufinden, ob und wie stark sich der Aligner auf das Varianten-Calling auswirkt und welche Aligner bei Verwendung einer normalen Illumina-Exom-Probe am besten abschneiden. Unseres Wissens ist dies der erste Bericht, der alle möglichen Kombinationen (insgesamt 30 Pipelines) einer breiten Palette von Alignern und Varianten-Callern validiert.

2. Methoden

2.1. Datensätze

Das menschliche Referenzgenom hg19 wurde vom UCSC-Browser (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) heruntergeladen und für die Durchführung der Alignments verwendet. Das menschliche Exom, SRR098401, wurde aus dem Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) heruntergeladen. Für die Annotation und Kalibrierung wurden dbSNP137 ohne Sites nach Version 129, HapMap 3.3, Human Omni 2.5 BeadChip und Mills und 1000 G Goldstandard Indel Set Listen verwendet (alle von ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. Die Pipeline

Abbildung 1 zeigt den in dieser Studie verwendeten Arbeitsablauf, der dem im Leitfaden „Best Practices“ des Broad Institute skizzierten Arbeitsablauf ähnelt. Er umfasst eine Reihe von Schritten, um sicherzustellen, dass die erzeugten Alignment-Dateien von höchster Qualität sind, sowie mehrere weitere Schritte, um zu gewährleisten, dass die Varianten korrekt benannt werden. Zunächst wurden die Rohdaten an hg19 ausgerichtet, und dann wurden PCR-Duplikate aus dem Alignment entfernt. Um die Indel-Identifizierung später in der Pipeline zu erleichtern, wurde anschließend ein Read-Realignment um Indels herum durchgeführt. Der letzte Schritt der Alignment-Verarbeitung bestand in der Rekalibrierung des Basenqualitäts-Scores, der dazu beiträgt, die inhärenten Verzerrungen und Ungenauigkeiten der von den Sequenzierern vergebenen Scores auszugleichen. Leider war die Alignment-Rate der einzelnen Aligner trotz dieser Schritte deutlich niedriger als erwartet. Um dies auszugleichen, wurde das fastx-Toolkit zum Herausfiltern von Reads geringer Qualität verwendet (Tabelle 1). Als minderwertige Reads wurden solche definiert, bei denen mindestens die Hälfte der Qualitätswerte unter 30 lag. Nach der Alignment-Verarbeitung wurden Variantenaufrufe und Variantenfilterung durchgeführt.

Abbildung 1

Schematische Darstellung der verwendeten Datenanalyse-Pipeline. Um sicherzustellen, dass Alignments von höchster Qualität erstellt werden, werden die Reads zunächst gefiltert und dann an das menschliche Referenzgenom hg19 angeglichen. Anschließend werden PCR-Duplikate entfernt, Reads um mutmaßliche Indels herum aligniert und die Basenqualitätswerte neu kalibriert. Schließlich werden die Varianten aufgerufen und anhand der NIST-GiaB-Liste der Varianten validiert.

Die sechs Tools, die zur Erzeugung von Alignments verwendet wurden, waren Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK und Novoalign, und die fünf Tools, die zur Erzeugung von Varianten verwendet wurden, waren FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup und SNPSVM, wie in Tabelle 2 zu sehen ist.

2.3. Filterung

Rohdaten wurden von der SRA (SRR098401) erfasst, mit fastq-dump aufgeteilt und mit dem fastx-Toolkit gefiltert. Bei fastq-dump wurden die Flags -split_files und -split_spot verwendet, und fastq_quality_filter wurde mit den folgenden Argumenten ausgeführt: -Q 33 -q 30 -p 50. Anschließend wurden die Reads mit einem benutzerdefinierten Skript korrekt gepaart.

2.4. Aligning

Die Aligner verwendeten Standardargumente, außer wenn ein Thread-Argument verwendet wurde, sofern verfügbar. Die verwendeten Befehle sind wie folgt.

2.4.1. Bowtie2

(1)bowtie2 -p 10 -x $INDEX -1 raw_data/read_1_filtered.fastq -2 raw_data/read_2_filtered.fastq -S alignments/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 -in alignments/NA12878.MOSAIK.mkb -out alignments/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 Depth of Coverage Calculation

Um eine korrekte Berechnung der Abdeckungstiefe zu gewährleisten, wurde das Picard Tools-Modul CalculateHsMetrics mit den folgenden Argumenten verwendet:(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.bedEs ist wichtig zu beachten, dass der TruSeq-Exome-Bettdatei der Header der SAM-Alignment-Datei vorangestellt werden muss, damit dieses Modul funktioniert. Außerdem muss Spalte 6 in Spalte 4 verschoben werden, und Spalte 5 muss aus der TruSeq-Bettdatei entfernt werden.

2.6. Verarbeitung der Ausrichtungsdateien

Die Verarbeitung der Ausrichtungsdateien (SAM/BAM-Dateien) erforderte für alle Aligner die folgenden Schritte:(1)SAM zu BAM Konvertierung mit SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Sortierung der BAM-Dateien mit dem Picard Tools-Modul SortSam:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam SORT_ORDER=coordinate(3)PCR-Duplikatentfernung mit dem Picard Tools-Modul MarkDuplicates:(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=Zuordnungen/NA12878.ALN.sorted.bam O=Zuordnungen/NA12878.ALN.dups_removed.bam REMOVE_DUPLICATES=true M=alignments/metrics(4)Mit dem Picard Tools-Modul AddOrReplaceReadGroups zu Alignment-Dateien hinzugefügte Lesegruppen:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=koordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Realignment um Indels unter Verwendung der GATK-Module RealignerTargetCreator und IndelRealigner:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T RealignerTargetCreator -R genome/hg19.fa -I alignments/NA12878.ALN.RG.bam -known genome/mills.vcf -o tmp/ALN.intervals(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T IndelRealigner -R genome/hg19.fa -I alignments/NA12878.ALN.RG.bam -known genome/mills.vcf -o alignments/NA12878.ALN.indels.bam -maxReadsForRealignment 100000 -maxReadsInMemory 1000000 -targetIntervals tmp/ALN.intervals(6)Basenrekalibrierung mit den GATK-Modulen BaseRecalibrator und 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. Variantenaufruf

Die Standardargumente wurden für jeden Variantenaufrufer verwendet, es sei denn, er enthielt ein „threads“- oder „parallel“-Flag; in diesem Fall wurde auch dieses verwendet. Außerdem wurden Indels nach Möglichkeit getrennt von SNVs aufgerufen. Im Einzelnen wurden die folgenden Befehle verwendet:

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 && 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.vcfDa die erforderlichen CIGAR-Flags in der Alignment-Datei nicht vorhanden waren, konnte SNPSVM keine Varianten für CUSHAW3 aufrufen, und SAMtools mpileup konnte aus demselben Grund keine Varianten für MOSAIK-Alignments aufrufen. Da SNPSVM nur SNVs erkennt, wurden auch keine Indels für dieses Programm gemeldet.

2.8. Variantenfilterung

Die Filterung variierte je nach verwendetem Varianten-Caller. Im Falle von GATK HaplotypeCaller und GATK UnifiedGenotyper wurden die GATK-Module VariantRecalibrator und ApplyRecalibration verwendet, um SNVs zu filtern, wobei HapMap 3.3, der Omni 2.5 SNP BeadChip und dbSNP 137 ohne 1000 Genome-Daten als Trainingssätze verwendet wurden. Für SNPSVM wurden QUAL-Scores ≥ 4 und DP-Werte ≥ 6 verwendet. Für FreeBayes und SAMtools wurden QUAL-Scores ≥ 20 und DP-Werte ≥ 6 verwendet.

2.9. Variantenvergleich

Für den Variantenvergleich wurde USeq 8.8.1 verwendet, um SNVs zu vergleichen, die zwischen allen Datensätzen gemeinsam sind. Zum Vergleich von Indels wurde das vcflib-Tool vcfintersect verwendet. Die TruSeq hg19-Exom-Bed-Datei truseq_exome_targeted_regions.hg19.bed.chr, die am 11. Dezember 2013 erhalten wurde, wurde verwendet, um Vergleiche auf Stellen zu beschränken, die mit dem Exom-Pull-Down-Kit erfasst werden konnten, das bei der Sequenzierung von SRR098401 verwendet wurde. Diese Datei kann hier von Illumina bezogen werden: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Um sicherzustellen, dass die Varianten in den verschiedenen Call-Sets identisch dargestellt werden, wurde das vcflib-Tool vcfallelicprimitives zur Vorverarbeitung der vcf-Dateien verwendet.

2.10. Statistische Berechnungen

True Positive (TP). Es handelt sich um eine Mutation, die von der getesteten Pipeline erkannt wurde und die in der NIST-GiaB-Liste enthalten ist.

Falsch positiv (FP). Es handelt sich um eine Mutation, die von der getesteten Pipeline erkannt wurde, aber nicht in der NIST-GiaB-Liste vorhanden ist.

True Negative (TN). Es handelt sich um eine Mutation, die von der getesteten Pipeline nicht erkannt wurde und die in der NIST-GiaB-Liste nicht vorhanden ist.

Falsch Negativ (FN). Es handelt sich um eine Mutation, die von der getesteten Pipeline nicht erkannt wurde, die aber in der NIST-GiaB-Liste vorhanden ist:

3. Ergebnisse und Diskussion

3.1. Vorfilterung von Varianten

Bei der Durchführung von Variantenanalysen ist eine der vielen Fallstricke, die berücksichtigt werden müssen, der Exom-Sequenzraum (wie vom Exom-Capture-Kit definiert) und wie er die Analyseergebnisse beeinflussen kann. In diesem Fall hatten wir ein einzelnes Exom (SRR098401), das mit dem Illumina TruSeq Exom-Kit extrahiert und auf einem HiSeq 2000 sequenziert wurde. Vor diesem Hintergrund wollten wir sicherstellen, dass wir die Fähigkeit der Bioinformatik-Tools messen, ihre Aufgaben zu erfüllen, und nicht, wie gut das Illumina TruSeq Exome Capture Kit funktioniert. Das heißt, wir wollen nur versuchen, Varianten zu nennen, die in den vom Pull-Down-Kit definierten Exons vorhanden sein sollen. Aus diesem Grund verwenden wir die von Illumina bereitgestellte Bettdatei und nicht eine allgemeine Annotationsbettdatei, z. B. RefSeq für hg19. Wir fanden heraus, dass es für dieses bestimmte Individuum gemäß der NIST-GiaB-Liste insgesamt 34.886 SNVs und 1.473 Indels innerhalb der von der TruSeq-Bettdatei definierten Regionen geben sollte.

Nachdem wir Varianten herausgefiltert hatten, die sich nicht in den von der Illumina TruSeq Exome-Bettdatei definierten Regionen befanden, gingen wir von Hunderttausenden von mutmaßlichen Varianten (Daten nicht gezeigt) auf durchschnittlich etwa 23.000 Varianten (SNVs und Indels) pro Pipeline über (Tabelle 3). Dies ist ein wichtiger Schritt für Forscher, da es den Suchraum für potenziell interessante Varianten erheblich reduziert.

3.2. Rohvariantenergebnisse im Vergleich zu GiaB

Ein Aspekt, den wir bei diesem Vergleich verstehen wollten, war die Bedeutung der Filterung von Varianten, die von diesen Tools erkannt wurden. Der Grund dafür ist, dass man idealerweise eine möglichst hohe Empfindlichkeit haben möchte, damit die interessierenden Mutationen nicht im Filterungsprozess verloren gehen. Es ist daher ratsam festzustellen, ob und in welchem Maße dieser Schritt notwendig ist, da aus den NIST-GiaB-Ergebnissen und dem Bericht von Bamshad et al. hervorgeht, dass die Sensitivität ein Problem darstellen könnte.

Wie aus Tabelle 3 hervorgeht, ist die Filterung bei einigen Varianten-Callern für SNVs notwendiger als für andere, und bei Indels ist sie absolut notwendig. In den meisten Fällen liegt die Zahl der TP-Varianten nahe bei oder höher als die von uns erwartete Zahl von etwa 20.000, aber andererseits ist die Zahl der FPs in einigen Fällen sehr hoch.

Die von den einzelnen Pipelines generierten Zahlen variieren sehr stark. Man kann jedoch einige Gemeinsamkeiten in den Zahlen finden, die wahrscheinlich auf die algorithmischen Ursprünge der einzelnen Tools zurückzuführen sind. FreeBayes erzeugt sowohl die größte Anzahl ungefilterter Varianten als auch die höchste Anzahl von FPs. Es ist wahrscheinlich, dass wir diese Art von Leistung nur bei diesem Tool sehen, da es zwar nicht der einzige Varianten-Caller ist, der auf Bayes’scher Inferenz basiert, aber einzigartig in seiner Interpretation von Alignments ist. Das heißt, es handelt sich um einen Haplotyp-basierten Caller, der Varianten auf der Grundlage der Sequenz der Reads selbst und nicht des Alignments identifiziert, wie es der UnifiedGenotyper von GATK tut.

Außerdem sehen wir, dass die Burrows-Wheeler-basierten Aligner sehr ähnlich abschneiden: Bowtie2, BWA mem und BWA sampe erzielen durchweg ähnliche Ergebnisse. Man könnte vermuten, dass dies wahrscheinlich darauf zurückzuführen ist, dass alle diese Tools ähnliche Algorithmen verwenden, wenn sie ihre Aufgabe erfüllen. Diese Beobachtung wird durch die Tatsache gestützt, dass MOSAIK (Gapped Alignments unter Verwendung des Smith-Waterman-Algorithmus) und CUSHAW3 (ein hybrider Seeding-Ansatz) beide sehr unterschiedliche Algorithmen zugrunde liegen und folglich sehr unterschiedliche Ergebnisse liefern.

Dieser Unterschied in den Ergebnissen, die mit verschiedenen Algorithmen korrelieren, zeigt sich am besten bei den SNPSVM-Ergebnissen. Von den Variantenaufrufern ist dies der einzige, der Support-Vektor-Maschinen und Modellbildung einsetzt, um SNV-Aufrufe zu generieren. Es hat den Anschein, dass diese Methode zwar den Nachteil hat, dass sie nicht so empfindlich ist wie andere Methoden, jedoch den Vorteil hat, dass sie unabhängig von dem verwendeten Aligner extrem genau ist. Dies deutet darauf hin, dass man den Filterungsschritt ganz überspringen kann, wenn man diesen Varianten-Caller verwendet.

In Bezug auf Indels scheint kein Aligner unter den anderen hervorzustechen, der diese Art von Mutation gut handhabt. Betrachtet man die Anzahl der FPs, so wird deutlich, dass der Variantenrufer die größte Rolle bei der Genauigkeit der Indel-Identifizierung spielt. Darüber hinaus gibt es weder für CUSHAW3 plus SNPSVM noch für MOSAIK plus SAMtools Mpileup-Pipelines Daten, da die Alignment-Dateien nicht die erforderlichen CIGAR-Zeichenfolgen enthalten, damit die Varianten-Caller nachgeschaltet funktionieren können. Schließlich gibt es keine Indel-Daten für SNPSVM, weil dieses Tool ausschließlich zur Identifizierung von SNVs verwendet wird.

3.3. Gefilterte Ergebnisse im Vergleich zu GiaB

Wie in Tabelle 4 dargestellt, gelingt es den Standard-Filterverfahren, eine große Anzahl von FP SNVs für jede Pipeline zu entfernen; es scheint jedoch, dass diese Filter in den meisten Fällen für die SNV-Erkennung etwas zu aggressiv, für Indels jedoch nicht streng genug sind. Dies wird deutlich, wenn man sich die Unterschiede in der Anzahl der FPs ansieht, die in jedem Datensatz gemeldet werden. Bei Bowtie2 mit Freebayes werden beispielsweise 72.570 FP SNVs entfernt (eine Reduktion von 98%), aber nur 1.736 FP Indels (eine Reduktion von 70%). Es ist auch anzumerken, dass die verwendeten Filter pipeline-abhängig waren und größtenteils innerhalb jeder Pipeline zu einer ähnlichen Reduzierung der SNV- und Indel-FPs führten. Die einzige Ausnahme ist die Anzahl der Varianten, die aus den CUSHAW3-Alignments im Vergleich zu anderen Alignments identifiziert wurden: Insgesamt ist die Zahl der TP-SNVs geringer, die Zahl der FP-SNVs höher, und es ist der einzige Aligner, der nach dem Filtern mehr als 1.000 FP-Indels produziert.

In Anbetracht der Tatsache, dass das Filtern die Anzahl der TP-Varianten erheblich reduziert, könnte es ratsam sein, mit Ausnahme der Pipelines, die CUSHAW3 und FreeBayes verwenden, diesen Schritt bei der Suche nach seltenen, hochwirksamen Varianten zu überspringen. Stattdessen könnte man mehr Zeit auf einen Filterungsprozess verwenden, der eher auf Biologie als auf Statistik basiert. So kann es beispielsweise sinnvoller sein, Zeit in die Identifizierung einer kleinen Liste von Varianten zu investieren, die wahrscheinlich hochwirksam sind: Spleißstellenmutationen, Indels, die Frameshifts verursachen, Trunkierungsmutationen, Stop-Loss-Mutationen oder Mutationen in Genen, von denen bekannt ist, dass sie biologisch relevant für den interessierenden Phänotyp sind.

3.4. Durchschnittliche TPR und Sensitivität

Wie aus Tabelle 5 hervorgeht, liegt der positive Vorhersagewert (PPV) für jedes Tool, mit Ausnahme von CUSHAW3, zwischen 91 % und 99,9 % für SNVs, aber die durchschnittliche Sensitivität ist sehr niedrig (etwa 50 %). Diese Diskrepanz könnte auf eine Reihe von Gründen zurückzuführen sein, der wahrscheinlichste ist jedoch die unterschiedliche Tiefe der Abdeckung der Exons. Wir sehen, dass neben der geringen SNV-Sensitivität auch die Indel-Sensitivität gering ist (ca. 30 %); der PPV für Indels ist jedoch deutlich niedriger (35,86 % bis 52,95 %). Dies könnte auf einen der folgenden Gründe zurückzuführen sein: Sehr kurze Indels sind mit herkömmlichen NGS schwer zu erkennen, die Darstellung von Indels durch verschiedene Varianten-Caller kann dazu führen, dass Tools fälschlicherweise behaupten, dass zwei Indels unterschiedlich sind, oder Alignment-Tools erzeugen unterschiedliche Darstellungen desselben Indels.

Die wahrscheinlichste Erklärung für beide Arten von Mutationen ist die Frage der Tiefe. Wie bei jeder Studie zur Variantenanalyse führt eine Erhöhung der Abdeckungstiefe zu einer Erhöhung der Sensitivität, aber es ist unmöglich, eine gute Abdeckungstiefe zu garantieren, da die Exome Capture Kits nicht in der Lage sind, Exons gleichmäßig zu erfassen. Außerdem deckt kein einziges Exom-Capture-Kit jedes Exon ab. In der Tat hat sich gezeigt, dass die Variantenanalyse bei der Sequenzierung des gesamten Genoms mit einer durchschnittlichen Tiefe, die geringer ist als die eines Exoms, aufgrund der Einheitlichkeit dieser Tiefe besser funktioniert. Daher ist es wahrscheinlich, dass eine große Anzahl von Varianten fehlt, da die NIST-GiaB-Liste aus einer Zusammenstellung von Exom- und Genomsequenzierungsdaten erstellt wurde. Um eine angemessene Sensitivität zu erreichen, wird man letztendlich eine Ganzgenomsequenzierung durchführen müssen, aber das ist für die meisten Labors derzeit zu teuer. Glücklicherweise sinken diese Kosten weiter, und wir werden bald eine allmähliche Verlagerung von der Exomanalyse zur vollständigeren Ganzgenomanalyse erleben.

3.5. Sensitivität in Abhängigkeit von der Tiefe

Da die Sensitivität eine der wichtigsten Leistungskennzahlen eines Tools darstellt und die meisten Tools Mühe haben, eine Sensitivität von mehr als 50 % zu erreichen, wollten wir weiter untersuchen, wie sich die Tiefe auf die Sensitivität der Variantenerkennung auswirkt. Wir untersuchten eine Reihe verschiedener Kombinationen von Tools, um die besten Pipelines, Varianten-Caller und Aligner zu ermitteln. Für Abbildung 2 wurden die fünf besten Kombinationen von Variantenaufrufern und Alignern anhand ihrer Empfindlichkeit und Falsch-Positiv-Rate (FPR) ermittelt. Das heißt, wir haben diejenigen ausgewählt, die die höchste Anzahl von TP SNVs und die niedrigste Anzahl von FP SNVs aufwiesen. Bei der Inspektion fällt sofort auf, dass die Empfindlichkeit geringer ist als erwartet. Alle Pipelines schneiden in etwa gleich gut ab: Sie identifizieren die meisten ihrer Varianten, wenn eine Tiefe von etwa 150x erreicht ist, was darauf hindeutet, dass diese Tiefe wahrscheinlich ausreicht und dass die Anzahl der fehlenden Varianten wahrscheinlich darauf zurückzuführen ist, dass bestimmte Exons eine unterdurchschnittliche Abdeckung aufweisen. Man beachte, dass vier der fünf leistungsstärksten Pipelines GATK UnifiedGenotyper als Varianten-Caller verwenden, was seine überlegene Leistung unabhängig vom verwendeten Aligner zeigt, wie in Abbildung 3(b) dargestellt.

Abbildung 2

Sensitivität als Funktion der Tiefe für die fünf besten Pipelines. Die fünf besten Pipelines sind hier dargestellt, wobei die Tiefe jedes SNVs gegen die Empfindlichkeit aufgetragen ist.

(a)

(b)

(a)
(b)

Abbildung 3

Sensitivität als Funktion der Tiefe für den Top-Aligner und den Top-Varianten-Caller. (a) Ergebnisse für die Tiefe jeder SNV aufgetragen gegen die Empfindlichkeit für den Top-Aligner, BWA mem, gepaart mit jedem Varianten-Caller. (b) Ergebnisse für die Tiefe jeder SNV, aufgetragen gegen die Sensitivität für den besten Variantenaufrufer, GATK UnifiedGenotyper, gepaart mit jedem Aligner.

Neben der Betrachtung der fünf besten Pipelines haben wir festgestellt, dass es sinnvoll wäre, die gleiche Analyse für den besten Aligner, der mit jedem Variantenaufrufer gekoppelt ist (Abbildung 3(a)), sowie für den besten Variantenaufrufer, der mit jedem Aligner gekoppelt ist (Abbildung 3(b)), durchzuführen. Wie bei den Pipelines wurde der beste Aligner als derjenige identifiziert, der die höchste Anzahl von TP-SNVs und die niedrigste Anzahl von FP-SNVs produziert – in diesem Fall BWA mem. Trotz des besten Alignments, mit dem wir arbeiten können, sehen wir immer noch einen ziemlich großen Unterschied zwischen den Variantenaufrufern, was wahrscheinlich auf die unterschiedlichen Algorithmen zurückzuführen ist, die sie verwenden (Abbildung 3(a)). Im Fall des leistungsstärksten Variantenaufrufers, GATK UnifiedGenotyper, scheint es jedoch weniger Unterschiede zwischen den vier besten Alignern zu geben, was darauf hindeutet, dass er in den meisten Situationen recht gut abschneidet, mit Ausnahme von CUSAHW3 und MOSAIK.

3.6. Gemeinsame Varianten in den Top-Pipelines

Zuletzt wollten wir wissen, wie einzigartig die Variantenaufrufe in den verschiedenen Pipelines sind. Zu diesem Zweck konzentrierten wir uns erneut auf die fünf führenden Pipelines für Variantenaufrufe: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper, und Novoalign plus UnifiedGenotyper. Wie in Abbildung 4 zu sehen ist, gibt es eine große Menge an Überschneidungen zwischen den fünf verschiedenen Pipelines, wobei 15.489 SNVs (70 %) von insgesamt 22.324 verschiedenen Varianten gemeinsam genutzt werden. Man könnte jedoch auch argumentieren, dass dies größtenteils darauf zurückzuführen ist, dass vier der fünf Pipelines den UnifiedGenotyper als ihren Varianten-Caller verwenden. Diese Annahme wird durch die Tatsache bestätigt, dass die größte Anzahl von Varianten, die nur in einer Pipeline vorkommen, nämlich 367, zu der Kombination BWA sampe plus HaplotypeCaller gehört. Es ist auch erwähnenswert, dass die zweithöchste Anzahl einzigartiger SNVs ebenfalls zum BWA sampe Aligner gehört, so dass es möglich ist, dass die hohe Anzahl einzigartiger SNVs eher dem Aligner als dem Varianten-Caller zuzuschreiben ist.

Abbildung 4

Die Schnittmenge der von den fünf führenden Pipelines identifizierten SNVs.

4. Schlussfolgerungen

Wir fanden heraus, dass unter den dreißig verschiedenen getesteten Pipelines Novoalign plus GATK UnifiedGenotyper die höchste Sensitivität aufwies und gleichzeitig eine niedrige Anzahl von FP für SNVs beibehielt. Von den Alignern schnitt BWA mem durchweg am besten ab, aber die Ergebnisse variierten je nach verwendetem Varianten-Caller dennoch stark. Daraus folgt natürlich, dass der beste Variantencaller, GATK UnifiedGenotyper, unabhängig vom verwendeten Aligner meist ähnliche Ergebnisse lieferte. Es ist jedoch leicht ersichtlich, dass Indels immer noch für jede Pipeline schwierig zu handhaben sind, da keine der Pipelines eine durchschnittliche Sensitivität von mehr als 33 % oder einen PPV von mehr als 53 % erreicht. Zusätzlich zu der geringen Gesamtleistung, die wir bei der Erkennung von Indels feststellen, ist die Sensitivität, unabhängig vom Mutationstyp, ein Problem für alle in dieser Arbeit vorgestellten Pipelines. Die erwartete Anzahl von SNVs für das Exom von NA12878 beträgt 34.886, aber selbst bei Verwendung der Vereinigung aller von den fünf besten Pipelines identifizierten Varianten war die größte identifizierte Anzahl sehr niedrig (22.324). Es scheint, dass die Exomanalyse, obwohl sie immer noch sehr nützlich ist, ihre Grenzen hat, selbst wenn es um etwas so scheinbar Einfaches wie die SNV-Erkennung geht.

Bekanntgabe

Adam Cornish ist Doktorand im Labor von Chittibabu Guda und hat eine Ausbildung in Informatik und Genomik. Chittibabu Guda (außerordentlicher Professor) hat einen interdisziplinären Hintergrund in Molekular- und Computerbiologie. Er hat seit 2001 eine Reihe von Berechnungsmethoden mit einer Vielzahl von Anwendungen in der biomedizinischen Forschung veröffentlicht.

Interessenkonflikt

Die Autoren sind sich keiner konkurrierenden Interessen bewusst.

Beitrag der Autoren

Adam Cornish hat die Studie entworfen, alle Analysen durchgeführt, die Abbildungen erstellt und die Arbeit geschrieben. Chittibabu Guda lieferte wesentliche Rückmeldungen zu Verbesserungen an der Arbeit und Beiträge zu den Analysen selbst und redigierte die Arbeit gründlich. Alle Autoren haben die endgültige Fassung gelesen und gebilligt.

Danksagungen

Diese Arbeit wurde durch die Entwicklungsgelder für Chittibabu Guda vom University of Nebraska Medical Center (UNMC) unterstützt. Die Autoren möchten den Entwicklern von Novoalign dafür danken, dass sie die Software als Testversion zur Verfügung gestellt haben, sowie der Bioinformatik- und Systembiologie-Kerneinrichtung am UNMC für die Infrastrukturunterstützung, die diese Forschung ermöglicht hat.