Abstract

Sekwencjonowanie o dużej przepustowości, zwłaszcza eksomów, jest popularnym narzędziem diagnostycznym, ale trudno jest określić, które narzędzia są najlepsze w analizie tych danych. W tym badaniu wykorzystujemy wyniki NIST Genome in a Bottle jako nowy zasób do walidacji naszego potoku analizy eksomów. Używamy sześciu różnych alignerów i pięciu różnych wywoływaczy wariantów, aby określić, który z 30 potoków osiąga najlepsze wyniki na ludzkim egzomie, który został użyty do wygenerowania listy wariantów wykrytych przez konsorcjum Genome in a Bottle. Spośród tych 30 potoków stwierdziliśmy, że Novoalign w połączeniu z GATK UnifiedGenotyper wykazuje najwyższą czułość, utrzymując jednocześnie niską liczbę wyników fałszywie dodatnich dla SNV. Widać jednak, że indele nadal są trudne do wykrycia przez każdy potok, ponieważ żadne z narzędzi nie osiągnęło średniej czułości wyższej niż 33% lub pozytywnej wartości predykcyjnej (PPV) wyższej niż 53%. Wreszcie, zgodnie z oczekiwaniami, stwierdzono, że alignery mogą odgrywać równie istotną rolę w wykrywaniu wariantów, jak same callery wariantów.

1. Wstęp

W ciągu ostatnich kilku lat dokonano wielu postępów w technologiach sekwencjonowania o wysokiej przepustowości. Dzięki tym postępom możliwe jest obecnie wykrycie dużej liczby potencjalnych wariantów wywołujących choroby, a w kilku przypadkach dane z sekwencjonowania następnej generacji (NGS) zostały nawet wykorzystane do celów diagnostycznych. Wynika to częściowo z rozwoju technologii sekwencjonowania w ciągu ostatnich kilku lat, ale także z liczby ulepszeń wprowadzonych do różnych narzędzi bioinformatycznych używanych do analizy gór danych produkowanych przez instrumenty NGS .

W przypadku poszukiwania mutacji u pacjenta, typowym przepływem pracy jest sekwencjonowanie ich egzomu za pomocą sekwenatora Illumina, wyrównanie surowych danych do ludzkiego genomu referencyjnego, a następnie zidentyfikowanie pojedynczych wariantów nukleotydowych (SNV) lub krótkich wstawek i delecji (indeli), które mogłyby prawdopodobnie spowodować lub wpłynąć na fenotyp zainteresowania. Podczas gdy jest to dość proste, decyzja o wyborze najlepszych narzędzi do użycia na każdym etapie analizy nie jest prosta. Istnieje duża liczba narzędzi, które są wykorzystywane w różnych etapach pośrednich, ale dwa najważniejsze etapy w całym procesie to wyrównanie surowych odczytów do genomu, a następnie wyszukiwanie wariantów (tj. SNV i indeli). W tym badaniu chcemy pomóc współczesnym bioinformatykom, wyjaśniając prawidłową kombinację narzędzia do wyrównywania krótkich odczytów i narzędzia do wywoływania wariantów w celu przetwarzania danych z sekwencjonowania eksomów wytwarzanych przez instrumenty NGS.

W przeszłości przeprowadzono wiele takich badań, ale wszystkie one miały wady w takiej czy innej formie. Idealnie byłoby mieć listę każdego znanego wariantu zawartego w próbce, tak że kiedy potok narzędzi analitycznych jest uruchomiony, można go przetestować, aby wiedzieć z pewnością, że działa poprawnie. Jednak w przeszłości taka lista nie istniała, więc walidacja musiała być przeprowadzana za pomocą mniej kompletnych metod. W niektórych przypadkach walidacja była przeprowadzana poprzez generowanie symulowanych danych, tak aby stworzyć zestaw znanych prawdziwych wyników pozytywnych (TP) i prawdziwych wyników negatywnych (TN). Chociaż w ten sposób wygodnie jest uzyskać listę każdego TP i TN w zbiorze danych, to jednak nie jest to dokładna reprezentacja biologii. Inne metody walidacji potoków wywoływania wariantów obejmują użycie tablic genotypowych lub sekwencjonowania Sangera w celu uzyskania listy TP i fałszywie pozytywnych (FP). Mają one tę zaletę, że dostarczają biologicznie potwierdzonych wyników, ale mają też tę wadę, że nie są wyczerpujące ze względu na ograniczoną liczbę miejsc na tablicach genotypowych i zaporowy koszt walidacji metodą Sangera, gdy wykonuje się ją tysiące razy. Wreszcie, żadne z tych badań nie miało na celu przyjrzenia się wpływowi, jaki na wywoływanie wariantów miał aligner krótkich odczytów. W związku z tym nie można było niezależnie ocenić wpływu wydajności alignera.

W tym badaniu mamy przewagę listy wariantów dla anonimowej kobiety z Utah (ID obiektu: NA12878, pierwotnie sekwencjonowanej w ramach projektu 1000 Genomes ), która została eksperymentalnie zatwierdzona przez prowadzone przez NIST konsorcjum Genome in a Bottle (GiaB). Ta lista wariantów została stworzona poprzez integrację 14 różnych zbiorów danych z pięciu różnych sekwenatorów i pozwala nam na walidację dowolnej listy wariantów wygenerowanej przez nasze potoki analizy egzomu. Nowością tej pracy jest walidacja właściwej kombinacji alignerów i wywoływaczy wariantów na podstawie wszechstronnego i eksperymentalnie określonego zbioru danych wariantów: NIST-GiaB.

Do przeprowadzenia naszej analizy wykorzystamy jeden z zestawów danych egzomowych pierwotnie użytych do stworzenia listy NIST-GiaB. Wybraliśmy tylko jeden z oryginalnych egzomów wygenerowanych przez Illumina TruSeq, ponieważ chcieliśmy przedstawić standardowy scenariusz przypadku użycia dla kogoś, kto chce przeprowadzić analizę NGS, i podczas gdy sekwencjonowanie całego genomu nadal spada w cenie, sekwencjonowanie egzomów jest nadal popularną i realną alternatywą. Ważne jest również, aby zauważyć, że według Bamshad i wsp. obecnie oczekiwana liczba SNV na egzom europejsko-amerykański wynosi 20 283 ± 523 . Pomimo tego, całkowita liczba SNV znalezionych na liście NIST-GiaB, które potencjalnie mogą występować w zbiorze danych egzomu TruSeq, wyniosła 34 886, co jest znacznie wyższe niż oczekiwano. Jest to prawdopodobnie spowodowane faktem, że podczas gdy zestaw exome został użyty do wygenerowania danych NIST-GiaB, został on również uzupełniony o sekwencjonowanie całego genomu.

Na koniec rozważyliśmy dużą liczbę alignerów i wywoływaczy wariantów, ale ostatecznie wybraliśmy 11 narzędzi opartych na powszechności, popularności i przydatności do naszego zbioru danych (np. SNVMix, VarScan2 i MuTect nie zostały użyte, ponieważ są przeznaczone do stosowania na próbkach pochodzących z nowotworów). Nasza analiza obejmuje porównanie sześciu alignerów (Bowtie2 , BWA sampe , BWA mem , CUSHAW3 , MOSAIK i Novoalign) oraz pięciu wywoływaczy wariantów (FreeBayes , GATK HaplotypeCaller, GATK UnifiedGenotyper , SAMtools mpileup i SNPSVM ). W tym badaniu próbujemy również określić, jak duży wpływ, jeśli w ogóle, ma aligner na wywołanie wariantu i które alignery działają najlepiej, gdy używamy normalnej próbki egzomu Illumina. Według naszej wiedzy, jest to pierwszy raport, który waliduje wszystkie możliwe kombinacje (w sumie 30 potoków) szerokiej gamy alignerów i wywoływaczy wariantów.

2. Metody

2.1. Datasets

Ludzki genom referencyjny hg19 został pobrany z przeglądarki UCSC (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/) i został użyty do wykonania alignerów. Ludzki egzom, SRR098401, został pobrany z Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra). Do celów anotacji i kalibracji użyto dbSNP137 bez miejsc po wersji 129, HapMap 3.3, Human Omni 2.5 BeadChip oraz listy zestawów indeli Mills i 1000 G gold standard (wszystkie z ftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz).

2.2. The Pipeline

Figura 1 pokazuje przepływ pracy zastosowany w tym badaniu, który jest podobny do tego przedstawionego w przewodniku Best Practices opracowanym przez The Broad Institute . Obejmuje on szereg kroków mających na celu zapewnienie najwyższej jakości plików z dopasowaniami, a także kilka kolejnych, gwarantujących prawidłowe nazwanie wariantów. Po pierwsze, surowe odczyty zostały wyrównane do hg19, a następnie z wyrównania usunięto duplikaty PCR. Następnie, aby ułatwić identyfikację indeli w dalszej części potoku, wykonano wyrównanie odczytów wokół indeli. Ostatnim krokiem przetwarzania wyrównania było przeprowadzenie rekalibracji wyników jakości bazy, co pomaga zniwelować nieodłączną stronniczość i niedokładność wyników uzyskiwanych przez sekwenatory. Niestety, pomimo tych kroków, stopień wyrównania każdego z alignerów był znacznie niższy niż oczekiwano, więc aby to zrównoważyć, użyto zestawu narzędzi fastx do odfiltrowania niskiej jakości odczytów (Tabela 1). Odczyty niskiej jakości zdefiniowano jako te odczyty, które miały przynajmniej połowę swoich wyników jakości poniżej 30. Po przetworzeniu wyrównania przeprowadzono wywoływanie wariantów i filtrowanie wariantów.

Rysunek 1

Schemat zastosowanego potoku analizy danych. Aby zapewnić najwyższą jakość dopasowań, odczyty są najpierw filtrowane, a następnie dopasowywane do ludzkiego genomu referencyjnego, hg19. Następnie usuwane są duplikaty PCR, odczyty są wyrównywane wokół domniemanych indeli, a wyniki jakości baz są rekalibrowane. Wreszcie, warianty są wywoływane i walidowane względem listy wariantów NIST-GiaB.

Sześć narzędzi użytych do generowania dopasowań to Bowtie2, BWA mem, BWA sampe, CUSHAW3, MOSAIK i Novoalign, a pięć narzędzi użytych do generowania wariantów to FreeBayes, GATK HaplotypeCaller, GATK UnifiedGenotyper, SAMtools mpileup i SNPSVM, jak widać w Tabeli 2.

2.3. Filtrowanie

Raw data was acquired from the SRA (SRR098401), split with fastq-dump, and filtered using the fastx toolkit. W szczególności, fastq-dump używał flag -split_files i -split_spot, a fastq_quality_filter był uruchamiany z następującymi argumentami: -Q 33 -q 30 -p 50. Następnie odczyty zostały odpowiednio sparowane za pomocą niestandardowego skryptu.

2.4. Wyrównywanie

Wyrównywacze używały domyślnych argumentów z wyjątkiem sytuacji, gdy użyto argumentu threads, gdy był on dostępny. Użyte polecenia są następujące.

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

Aby zapewnić prawidłowe obliczenie głębokości pokrycia, użyto modułu Picard Tools CalculateHsMetrics z następującymi argumentami:(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.bedWarto zauważyć, że aby ten moduł działał, do pliku z łóżkiem egzomu TruSeq musi być dołączony nagłówek z pliku SAM alignment. Ponadto kolumna 6 musi zostać przeniesiona do kolumny 4, a kolumna 5 musi zostać usunięta z pliku złoża TruSeq.

2.6. Przetwarzanie plików wyrównania

Przetwarzanie plików wyrównania (pliki SAM/BAM) wymagało następujących kroków dla wszystkich alignerów:(1)Konwersja SAM do BAM za pomocą SAMtools view:(a)samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Sortowanie plików BAM za pomocą modułu 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=coordinate(3)Usuwanie duplikatów PCR za pomocą modułu 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 dodane do plików wyrównania przy użyciu modułu Picard Tools, AddOrReplaceReadGroups:(a)java -jar bin/AddOrReplaceReadGroups.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.dups_removed.bam O=alignments/NA12878.ALN.RG.bam SO=coordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)Wyrównanie wokół indeli przy użyciu modułów GATK RealignerTargetCreator i 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)Rekalibracja bazy przy użyciu modułów GATK BaseRecalibrator i 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. Wywoływanie wariantów

Domyślne argumenty były używane dla każdego wywoływacza wariantów, chyba że zawierał on flagę „threads” lub „parallel”, w którym to przypadku również była ona używana. Dodatkowo, indele były wywoływane oddzielnie od SNV, jeśli było to możliwe. Konkretnie, użyte polecenia są następujące.

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.vcfZ powodu braku wymaganych flag CIGAR w pliku wyrównania, SNPSVM nie zdołał wywołać wariantów dla CUSHAW3, a SAMtools mpileup nie mógł wywołać wariantów na wyrównaniach MOSAIK z tego samego powodu. Ponadto, ze względu na fakt, że SNPSVM wykrywa tylko SNV, nie odnotowano żadnych indeli dla tego programu.

2.8. Filtracja wariantów

Filtracja różniła się w zależności od używanego wywoływacza wariantów. W przypadku GATK HaplotypeCaller i GATK UnifiedGenotyper, moduły GATK, VariantRecalibrator i ApplyRecalibration, zostały użyte do filtrowania SNV przy użyciu HapMap 3.3, Omni 2.5 SNP BeadChip i dbSNP 137 bez danych 1000 Genome jako zestawów treningowych. Dla SNPSVM użyto wyników QUAL ≥ 4 i wartości DP ≥ 6. Dla FreeBayes i SAMtools, QUAL scores ≥ 20 i DP values ≥ 6 were used.

2.9. Variant Comparison

Do porównania wariantów użyto USeq 8.8.1 do porównania SNV wspólnych dla wszystkich zestawów danych. Do porównania indeli użyto narzędzia vcflib – vcfintersect. Plik TruSeq hg19 exome bed truseq_exome_targeted_regions.hg19.bed.chr, uzyskany 11 grudnia 2013 r., został użyty w celu ograniczenia porównań do lokalizacji, które mogły być wychwycone przez zestaw exome pull down użyty w sekwencjonowaniu SRR098401. Plik ten można uzyskać od firmy Illumina tutaj: http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. Aby zapewnić, że warianty były reprezentowane identycznie między różnymi zestawami wywołań, użyto narzędzia vcflib vcfallelicprimitives do wstępnego przetworzenia plików vcf.

2.10. Obliczenia statystyczne

True Positive (TP). Jest to mutacja, która została wykryta przez testowany potok i znajduje się na liście NIST-GiaB.

False Positive (FP). Jest to mutacja, która została wykryta przez testowany potok, ale jest to mutacja, która nie istnieje na liście NIST-GiaB.

True Negative (TN). Jest to mutacja, która nie została wykryta przez testowany potok i jest jedną z tych, które nie istnieją na liście NIST-GiaB.

Fałszywy Negatywny (FN). Jest to mutacja, która nie została wykryta przez testowany potok, ale która istnieje na liście NIST-GiaB:

3. Wyniki i Dyskusja

3.1. Prefiltrowanie wariantów

Podczas wykonywania analizy wariantów jedną z wielu pułapek, które należy wziąć pod uwagę, jest przestrzeń sekwencji egzomu (określona przez zestaw do przechwytywania egzomu) i to, jak może ona wpłynąć na wyniki analizy. W tym przypadku dysponowaliśmy pojedynczym egzomem (SRR098401), który został pobrany przy użyciu zestawu do pozyskiwania egzomów Illumina TruSeq i zsekwencjonowany na aparacie HiSeq 2000. Mając to na uwadze, chcieliśmy się upewnić, że mierzymy zdolność narzędzi bioinformatycznych do wykonywania swojej pracy, a nie to, jak dobrze działał zestaw do pozyskiwania egzomu Illumina TruSeq. Oznacza to, że chcemy tylko próbować wywoływać warianty, które powinny być obecne w eksonach, jak zdefiniowano w zestawie do ściągania. Z tego powodu używamy pliku bed dostarczonego przez Illumina, a nie ogólnego pliku bed z anotacją, na przykład RefSeq dla hg19. Stwierdziliśmy, że dla tego konkretnego osobnika, zgodnie z listą NIST-GiaB, powinno być w sumie 34 886 SNV i 1 473 indeli w regionach zdefiniowanych przez plik łóżka TruSeq.

Odkąd odfiltrowaliśmy warianty, które nie znajdowały się w regionach zdefiniowanych przez plik łóżka egzomu Illumina TruSeq, przeszliśmy od setek tysięcy prawdopodobnych wariantów (dane nie pokazane) do, średnio, około 23 000 wariantów (SNV i indeli) na rurociąg (Tabela 3). Jest to ważny krok dla badaczy, od którego powinni zacząć, ponieważ znacznie zmniejsza przestrzeń wyszukiwania potencjalnie interesujących wariantów.

3.2. Raw Variant Results Compared to GiaB

Jednym z aspektów, który chcieliśmy zrozumieć podczas tego porównania, było znaczenie filtrowania wariantów wykrytych przez te narzędzia. Powodem tego jest to, że w idealnym przypadku chcielibyśmy mieć tak wysoki poziom czułości, jak to tylko możliwe, aby interesujące nas mutacje nie zostały utracone w procesie filtrowania. It therefore behoves us to determine whether or not this step is necessary and to what degree it is necessary, since it is clear from the NIST-GiaB results and the Bamshad et al. review that sensitivity could be an issue.

As we can see in Table 3, filtering is needed more for some variant callers than for others when it comes to SNVs, and it is absolutely necessary for indels. W większości przypadków liczba wariantów TP jest bliska lub wyższa niż oczekiwana przez nas liczba około 20 000, ale z drugiej strony, w niektórych przypadkach liczba FP jest bardzo wysoka.

Jasno widać, że istnieje wiele różnic w liczbach generowanych przez każdy potok. Można jednak znaleźć w nich pewne cechy wspólne, które prawdopodobnie wynikają z algorytmicznego pochodzenia każdego narzędzia. FreeBayes generuje zarówno największą liczbę niefiltrowanych wariantów, jak i największą liczbę FP. Jest prawdopodobne, że widzimy tego rodzaju wydajność tylko dzięki temu, że chociaż nie jest to jedyny program do wywoływania wariantów oparty na wnioskowaniu bayesowskim, jest on wyjątkowy w swojej interpretacji dopasowań. Oznacza to, że jest to wywoływacz oparty na haplotypach, który identyfikuje warianty w oparciu o sekwencję samych odczytów, a nie ich wyrównanie, tak jak działa UnifiedGenotyper firmy GATK.

Dodatkowo widzimy, że wyrównywacze oparte na Burrows-Wheeler działają bardzo podobnie do siebie: Bowtie2, BWA mem i BWA sampe osiągają podobne wyniki w całej rozciągłości. Można przypuszczać, że wynika to prawdopodobnie z faktu, że wszystkie te narzędzia wykorzystują podobne algorytmy podczas wykonywania swoich zadań. Obserwacja ta jest poparta faktem, że MOSAIK (gapped alignments using the Smith-Waterman algorithm) i CUSHAW3 (a hybrid seeding approach) mają bardzo różne algorytmy bazowe i w konsekwencji dają bardzo różne wyniki.

Ta różnica w wynikach korelująca z różnymi algorytmami jest najlepiej widoczna w wynikach SNPSVM. Spośród wywoływaczy wariantów, jest to jedyny, który wykorzystuje maszyny wektorów wspierających i budowanie modeli do generowania wywołań SNV. Wydaje się, że chociaż ma on wadę polegającą na tym, że nie jest tak czuły jak inne metody, to korzysta z tego, że jest niezwykle dokładny niezależnie od używanego alignera. This suggests that one is able to skip the filtering step altogether when using this variant caller.

With regard to indels, no aligner seems to stand out among the rest as one that handles this type of mutation well. W rzeczywistości, patrząc na liczbę FP, jasne jest, że to wywoływacz wariantu odgrywa największą rolę w dokładności identyfikacji indeli. Dodatkowo, nie ma danych ani dla CUSHAW3 plus SNPSVM ani MOSAIK plus SAMtools mpileup pipelines z powodu plików wyrównania nie zawierających niezbędnych ciągów CIGAR dla wywoływaczy wariantów, aby funkcjonować downstream. Wreszcie, powodem braku danych indel dla SNPSVM jest fakt, że narzędzie to jest używane wyłącznie do identyfikacji SNV.

3.3. Przefiltrowane wyniki w porównaniu z GiaB

Jak widać w Tabeli 4, standardowe praktyki filtrowania zdołały usunąć dużą liczbę FP SNV dla każdego potoku; wydaje się jednak, że filtry te są nieco zbyt agresywne w większości przypadków dla wykrywania SNV, ale nie dość rygorystyczne dla indeli. Staje się to oczywiste, gdy spojrzymy na różnice w liczbie FP odnotowanych w każdym zestawie danych. Na przykład Bowtie2 z Freebayes pozwala na usunięcie 72 570 FP SNV (redukcja o 98%), ale tylko 1 736 FP indeli (redukcja o 70%). Należy również zauważyć, że zastosowane filtry były zależne od potoku i w większości przypadków w ramach każdego potoku dały podobną redukcję FP SNV i indeli. Jedynym wyjątkiem jest liczba wariantów zidentyfikowanych z dopasowania CUSHAW3 w porównaniu z innymi dopasowaniami: Ogólnie liczba TP SNV jest niższa, liczba FP SNV jest wyższa i jest to jedyny aligner, który produkuje ponad 1000 FP indeli po odfiltrowaniu.

Zważywszy na fakt, że filtrowanie znacząco zmniejsza liczbę wariantów TP, rozsądnie byłoby, z wyjątkiem potoków wykorzystujących CUSHAW3 i FreeBayes, pominąć ten krok przy poszukiwaniu rzadkich, wysoce wpływowych wariantów. Zamiast tego można poświęcić więcej czasu na proces filtrowania, który opiera się raczej na biologii niż na statystyce. Na przykład, bardziej sensowne może być zainwestowanie czasu w zidentyfikowanie małej listy wariantów, które prawdopodobnie będą miały duży wpływ: mutacje w miejscu splotu, indele, które powodują zmiany ramek, mutacje truncacji, mutacje stop-loss lub mutacje w genach, które są biologicznie istotne dla interesującego nas fenotypu.

3.4. Średnia TPR i czułość

Jak widać w tabeli 5, pozytywna wartość predykcyjna (Positive Predictive Value, PPV) dla każdego narzędzia, z wyjątkiem CUSHAW3, waha się od 91% do 99,9% dla SNV, ale średnia czułość jest bardzo niska (około 50%). Ta rozbieżność może wynikać z wielu powodów, ale najbardziej prawdopodobnym jest zmienna głębokość pokrycia eksonów. Widzimy, że oprócz niskiej czułości SNV, czułość indeli jest niska (około 30%); jednakże PPV dla indeli jest znacznie niższa (35.86% do 52.95%). Może to wynikać z jednego z następujących powodów: bardzo krótkie indele są trudne do wykrycia przez konwencjonalne NGS , reprezentacja indeli przez różne variant callers może powodować, że narzędzia nieprawidłowo twierdzą, że dwa indele są różne, lub narzędzia wyrównujące produkują różne reprezentacje tego samego indelu .

Prawdopodobnie najbardziej prawdopodobnym wyjaśnieniem dla obu typów mutacji jest kwestia głębokości. Jak w przypadku każdego badania analizy wariantów, wzrost głębokości pokrycia prowadzi do wzrostu czułości, ale niemożliwe jest zagwarantowanie dobrej głębokości pokrycia ze względu na niezdolność zestawów do przechwytywania eksomu do równomiernego ściągania eksonów. Dodatkowo, żaden pojedynczy zestaw do wychwytywania eksomów nie obejmuje każdego eksonu. Rzeczywiście, wykazano, że analiza wariantów sekwencjonowania całego genomu przy średniej głębokości mniejszej niż egzom działa lepiej ze względu na jednorodność wspomnianej głębokości. Jest więc prawdopodobne, że brakuje dużej liczby wariantów ze względu na fakt, że lista NIST-GiaB została utworzona z kompilacji danych z sekwencjonowania egzomicznego i genomicznego. Ostatecznie, aby osiągnąć właściwą czułość, trzeba będzie wykonać sekwencjonowanie całego genomu, ale jest to obecnie nieopłacalne dla większości laboratoriów. Na szczęście koszt ten nadal spada i wkrótce zobaczymy stopniowe przejście od analizy egzomu do bardziej kompletnej analizy całego genomu.

3.5. Sensitivity as a Function of Depth

Ponieważ czułość odzwierciedla jeden z najważniejszych wskaźników wydajności narzędzia, a większość narzędzi walczy o osiągnięcie czułości wyższej niż 50%, chcielibyśmy dalej zbadać jak głębokość wpływa na czułość wywoływania wariantów. Przyjrzeliśmy się wielu różnym kombinacjom narzędzi, aby ustalić, jakie są najlepsze potoki, narzędzia do wywoływania wariantów i alignery. Na Rysunku 2 wybraliśmy pięć najlepszych kombinacji narzędzi do wywoływania wariantów i wyrównywania, określonych na podstawie ich czułości i współczynnika fałszywych wyników pozytywnych (FPR). Oznacza to, że wybraliśmy te, które miały największą liczbę wywołanych SNV TP i najmniejszą liczbę SNV FP. Podczas inspekcji, rzeczą, która wyróżnia się natychmiast jest to, że czułość jest niższa niż oczekiwano. Wszystkie potoki osiągają mniej więcej ten sam poziom: identyfikują większość wariantów do czasu osiągnięcia głębokości około 150x, co wskazuje, że głębokość ta jest prawdopodobnie wystarczająca, a liczba brakujących wariantów jest prawdopodobnie spowodowana tym, że niektóre eksony mają niższe niż przeciętne pokrycie. Należy zauważyć, że cztery z pięciu najlepiej działających potoków mają GATK UnifiedGenotyper jako swój wywoływacz wariantów, wykazując jego wyższą wydajność niezależnie od użytego alignera, jak pokazano na Rysunku 3(b).

Rysunek 2

Czułość jako funkcja głębokości dla pięciu najlepszych potoków. Pięć najlepszych rurociągów pokazano tutaj z głębokością każdego SNV wykreśloną względem czułości.

(a)

(b)

(a)
(b)

Oprócz przyjrzenia się pięciu najlepszym potokom, stwierdziliśmy, że warto przeprowadzić tę samą analizę na najlepszym alignerze połączonym z każdym wywoływaczem wariantów (Rysunek 3(a)), jak również na najlepszym wywoływaczu wariantów połączonym z każdym alignerem (Rysunek 3(b)). Podobnie jak w przypadku rurociągów, najlepszy aligner został zidentyfikowany jako ten, który wytworzył największą liczbę TP SNV i najmniejszą liczbę FP SNV – w tym przypadku BWA mem. Pomimo posiadania najlepszego wyrównania, nadal obserwujemy dość dużą różnicę między wywoływaczami wariantów, co prawdopodobnie wynika z różnych algorytmów, jakie stosują (Rysunek 3(a)). Jednakże, w przypadku najlepiej działającego wywoływacza wariantów, GATK UnifiedGenotyper, występuje mniejsza zmienność wśród czterech najlepszych wyrównywaczy, wskazując, że działa on dość dobrze w większości sytuacji, z wyjątkami w postaci CUSAHW3 i MOSAIK.

3.6. Wspólne warianty wśród najlepszych potoków

Na koniec chcieliśmy się dowiedzieć, jak unikalne są zestawy wywołań wariantów między różnymi potokami. Aby to zrobić, ponownie skupiliśmy się na pięciu najlepszych potokach wywoływania wariantów: Bowtie2 plus UnifiedGenotyper, BWA mem plus UnifiedGenotyper, BWA sampe plus HaplotypeCaller, BWA sampe plus UnifiedGenotyper oraz Novoalign plus UnifiedGenotyper. Jak widać na Rysunku 4, istnieje duża ilość nakładania się pomiędzy pięcioma różnymi potokami, z 15,489 SNV (70%) współdzielonymi z 22,324 odrębnych wariantów. Jednakże, można również stwierdzić, że jest to w dużej mierze spowodowane tym, że cztery z pięciu potoków używają UnifiedGenotyper jako swojego wywoływacza wariantów. Tę tezę potwierdza fakt, że największa liczba wariantów unikalnych dla potoku, 367, należy do kombinacji BWA sampe plus HaplotypeCaller. Warto również zauważyć, że druga najwyższa liczba unikalnych SNV również należy do alignera BWA sampe, więc możliwe, że wysoka liczba unikalnych SNV jest lepiej przypisana alignerowi niż variant caller.

Rysunek 4

Przecięcie SNV zidentyfikowanych przez pięć najlepszych rurociągów.

4. Wnioski

Zauważyliśmy, że spośród trzydziestu różnych testowanych pipelin Novoalign plus GATK UnifiedGenotyper wykazywały najwyższą czułość przy zachowaniu niskiej liczby FP dla SNV. Spośród alignerów, BWA mem konsekwentnie osiągał najlepsze wyniki, ale wyniki nadal różniły się znacznie w zależności od zastosowanego wywoływacza wariantów. Naturalnie, wynika z tego, że najlepszy wywoływacz wariantów, GATK UnifiedGenotyper, w większości przypadków dawał podobne wyniki niezależnie od użytego alignera. Jest jednak oczywiste, że indele są nadal trudne do opanowania przez każdy potok, ponieważ żaden z nich nie osiągnął średniej czułości wyższej niż 33% lub PPV wyższej niż 53%. Poza niską ogólną wydajnością, jaką widzimy w wykrywaniu indeli, czułość, niezależnie od typu mutacji, stanowi problem dla każdego potoku opisanego w niniejszym dokumencie. Spodziewana liczba SNV dla egzomu NA12878 wynosi 34 886, ale nawet przy użyciu unii wszystkich wariantów zidentyfikowanych przez pięć najlepszych potoków, największa zidentyfikowana liczba była bardzo niska (22 324). Wydaje się, że choć wciąż bardzo użyteczna analiza egzomu ma swoje ograniczenia, nawet jeśli chodzi o coś tak pozornie prostego jak wykrywanie SNV.

Ujawnienie

Adam Cornish jest studentem studiów magisterskich w laboratorium Chittibabu Guda z wykształceniem w dziedzinie informatyki i genomiki. Chittibabu Guda (profesor nadzwyczajny) posiada interdyscyplinarne wykształcenie w dziedzinie biologii molekularnej i obliczeniowej. Od 2001 roku opublikował wiele metod obliczeniowych o różnorodnych zastosowaniach w badaniach biomedycznych.

Konflikt interesów

Autorzy nie są świadomi żadnych konkurencyjnych interesów.

Wkład autorów

Adam Cornish zaprojektował badanie, przeprowadził wszystkie analizy, wykonał ryciny i napisał pracę. Chittibabu Guda dostarczył istotnych informacji zwrotnych na temat poprawek w pracy i wkładu w same analizy oraz dokładnie zredagował pracę. Wszyscy autorzy przeczytali i zatwierdzili ostateczną wersję pracy.

Podziękowania

Ta praca była wspierana przez fundusze rozwojowe dla Chittibabu Guda z University of Nebraska Medical Center (UNMC). Autorzy pragną podziękować twórcom Novoalign za udostępnienie oprogramowania w wersji testowej oraz Bioinformatics and Systems Biology Core facility w UNMC za wsparcie infrastrukturalne ułatwiające przeprowadzenie tych badań.