Articles

A Comparison of Variant Calling Pipelines Using Genome in a Bottle as a Reference

Abstract

ハイスループットのシーケンサー、特にエクソームは人気の診断ツールだが、どのツールがこのデータの解析に最も適しているかを判断することは困難である。 本研究では、エクソーム解析パイプラインの検証のための新しいリソースとして、NIST Genome in a Bottleの結果を使用する。 6種類のアライナーと5種類のバリアントコーラーを使用し、Genome in a Bottle Consortiumが検出したバリアントリストの作成に使用したヒトエクソームにおいて、全30種類のうちどのパイプラインが最も優れた性能を示すかを判断しています。 これら30のパイプラインのうち、GATK UnifiedGenotyperと組み合わせたNovoalignは、SNVの偽陽性を低く抑えながら最高の感度を示すことが分かりました。 しかし、どのツールも平均33%以上の感度、53%以上のPPV(Positive Predictive Value)を達成できず、どのパイプラインも依然としてインデルを扱うのは困難であることが明らかになりました。 最後に、予想通り、アライナーはバリアントコーラーそのものと同様にバリアント検出において重要な役割を果たすことが判明した。 背景

過去数年間で、ハイスループットシーケンス技術に多くの進歩があった。 これらの進歩により、現在では多くの潜在的な疾患原因バリアントを検出することが可能であり、いくつかのケースでは次世代シーケンサー(NGS)データが診断目的に使用されているほどである。 これは、過去数年間のシーケンス技術の発展によるものですが、NGS装置によって生成された膨大なデータを解析するために使用される様々なバイオインフォマティクスツールに多くの改良が加えられたことにも起因しています。

患者の変異を検索する場合、典型的なワークフローは、イルミナシーケンサーで患者のエクソームを配列し、生データをヒト参照ゲノムに合わせ、目的の表現型の原因または影響を及ぼす可能性のある一塩基変異(SNV)または短い挿入および欠失(indel)を特定することです。 これは非常に簡単なことですが、解析パイプラインの各段階で使用する最適なツールを決定することは容易ではありません。 様々な中間ステップで使用されるツールは多数ありますが、プロセス全体の中で最も重要な2つのステップは、生リードをゲノムにアラインメントすることと、バリアント(SNVsやindel)を検索することです。 本研究では、NGS装置で作成されたエクソームシーケンスデータを処理するためのショートリードアライメントツールとバリアントコーリングツールの正しい組み合わせを明らかにすることで、今日のバイオインフォマティシャンに役立つことを目的としています

過去にも多くの研究が行われてきましたが、いずれも何らかの形で欠点がありました。 理想的には、解析ツールのパイプラインを実行するときに、それが正しく機能していることを確実に知るためにテストできるように、サンプルに含まれるすべての既知のバリアントのリストがあるべきです。 しかし、過去にはそのようなリストは存在しなかったため、完全ではない方法で検証を行う必要がありました。 例えば、シミュレーションデータを作成し、既知の真陽性(TP)と真陰性(TN)のセットを作成することで検証を行うことができた。 この方法では、データセット内のすべてのTPとTNをリストアップすることができるが、生物学を正確に表現することはできない。 バリアントコールパイプラインを検証する他の方法としては、ジェノタイピングアレイやサンガーシーケンスを用いて、TPとFPのリストを得ることができます。 これらは生物学的に検証された結果を得られるという利点があるが、ジェノタイピングアレイのスポット数が限られていることや、サンガーバリデーションを何千回も行うとコストが高くなることから、包括的でないという欠点もある。 最後に、これらの研究はいずれも、ショートリードのアライナーがバリアントコールに与える影響を調べることを目的としていない。 本研究では、NIST主導のGenome in a Bottle (GiaB) Consortiumによって実験的に検証されたユタ州の匿名女性(被験者ID:NA12878、もともとは1000ゲノムプロジェクトで配列決定)のバリアントリストを利用しています。 このバリアントリストは、5つのシーケンサーから得られた14の異なるデータセットを統合して作成され、私たちのエクソーム解析パイプラインによって生成されたバリアントリストを検証することを可能にします 。 本研究の新規性は、アライナーとバリアントコーラーの適切な組み合わせを、包括的かつ実験的に決定されたバリアントデータセットに対して検証したことにあります。 NIST-GiaB.

解析には、もともとNIST-GiaBリストの作成に使用されたエクソームデータセットのうちの1つを使用します。 全ゲノムシーケンスの価格は下がり続けていますが、エクソームシーケンスはまだ人気があり、実行可能な代替手段です。 また、Bamshadらによると、現在、ヨーロッパ系アメリカ人のエクソームあたりのSNVの予想数は20,283 ± 523であることに注目することが重要である。 にもかかわらず、NIST-GiaBリストで見つかった、TruSeqエクソームデータセットに存在する可能性のあるSNVの総数は34,886であり、予想よりも大幅に多くなっています。 これは、NIST-GiaBデータの生成にエクソームキットが使用された一方で、全ゲノムシーケンスによって補完されたという事実によるものと思われます。

最後に、多数のアライナーとバリアントコーラーを検討しましたが、普及率、人気、および当社のデータセットとの関連性に基づいて、最終的に11のツールを選択しました(たとえば、SNVMix、BarScan2、および MuTectは腫瘍由来試料への使用が目的なので使用されていない)。 解析自体は、6つのアライナー(Bowtie2、BWA sampe、BWA mem、CUSHAW3、MOSAIK、Novoalign)と5つのバリアントコーラー(FreeBayes、GATK HaplotypeCaller、GATK UnifiedGenotyper、SAMtools mpileup、SNPSVM )と比較するものである。 また、本研究では、アライナーがバリアントコールにどの程度影響するのか、また、通常のイルミナエクソームサンプルを使用した場合にどのアライナーが最も良いパフォーマンスを示すのかを明らかにしようと試みています。 私たちの知る限り、これは幅広いアライナーとバリアントコーラーのすべての可能な組み合わせ(合計30パイプライン)を検証した最初の報告である。 データセット

ヒト参照ゲノムhg19はUCSCブラウザ(http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/)からダウンロードし、アライメントに使用した。 ヒトエキソーム SRR098401 は Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) からダウンロードした。 アノテーションとキャリブレーションのために、バージョン129以降のサイトを含まないdbSNP137、HapMap 3.3、Human Omni 2.5 BeadChip、Mills and 1000 G gold standard indel set listを用いた(すべてftp://ftp.broadinstitute.org/distribution/gsa/gatk_resources.tgz)。

2.2. パイプライン

図1は、本研究で使用したワークフローを示しており、これはThe Broad Instituteが作成したBest Practicesガイドに概説されているものと同様である。 これは、作成されたアラインメントファイルが最高品質であることを保証するためのいくつかのステップと、バリアントが正しくコールされることを保証するためのいくつかのステップを含んでいる。 まず、生のリードを hg19 にアライメントし、PCR の重複をアライメントから削除しました。 次に、パイプラインの後半でインデルの同定を助けるために、インデルの周囲でリードの再アラインメントが行われました。 アライメント処理の最後のステップは、塩基品質スコアの再キャリブレーションステップで、シーケンサーが発行するスコアの固有の偏りや不正確さを改善するのに役立ちます。 残念ながら、これらのステップにもかかわらず、各アライナーのアライメント率は予想よりもかなり低かったため、これを相殺するために、fastx toolkitを使用して低品質リードをフィルタリングしました(表1)。 低品質リードは、品質スコアの少なくとも半分が30以下のリードと定義されました。 アライメント処理後、バリアントコールとバリアントフィルタリングを行いました。

Aligner % read aligned (unfiltered)

% read aligned (filtered)

BWAメモリー89

47.69

Average depth of coverage
Bowtie2 89.73 98.73 47.97
BWA mem 92.91 99.85 46.97 BWA mem 98.73 47.97
BWA sampe 85.95 97.49 46.67
CUSHAW3 85.85.85.00 99.81
mosaik 85.68 96.22 45.67
Mosaik 99.81 47.69 47.69 47.6914
Novoalign 82.21 94.20 45.62
Table 1
フィルタリング・リードと非フィルター・リードに対するアライメントのパーセンテージです。 平均カバー深度は、フィルター付きリードで作成したアライメントファイルのものです。
図1
用いたデータ解析パイプラインの概念図です。 最高品質のアラインメントを作成するために、まずリードをフィルタリングし、次にヒトの参照ゲノムであるhg19にアラインメントする。 次に、PCR重複を除去し、推定されるインデルの周囲にリードをアライメントし、塩基品質スコアを再キャリブレーションします。 最後に、Variantをコールし、NIST-GiaB list of variantsと照合して検証しています。

アライメントに使用したツールは、Bowtie2、BWA mem、BWA sampe、CUSHAW3、MOSAIK、Novoalignの6つで、バリアント生成にはFreeBayes、GATK HaplotypeCaller、GATK UnifiedGenotyper、SAMtools mpileup、SNPSVMの5つが使用されており、表2に見られるとおり、このツールで生成したバリアントも、BWA mem、MOSAIK、SNPSVM、SAMtools mpileupの5つである。

0.7.5a

2.2.3

2.7-2

0.1.19

Tool Type Version Reference
Bowtie2 Aligner 2.X.S.A. Bowtie2 Aligner
Version
BWA sampe Aligner
BWA mem Aligner 0.7.5a
CUSHAW3 Aligner 3.0.3
MOSAIK Aligner
Novoalign Aligner 3.3.0.02.07 N/A
FreeBayes Genotyper v9.0 FreeBayes Genotyper50579.2-19-g011561f
GATK HaplotypeCaller Genotyper N/A
GATK UnifiedGenotyper Genotyper 2.2.0 2.7-2 2.7-2 2.7-2 2.7-2 2.7-2 2.7-2
SAMtools mpileup Genotyper
SNPSVM Genotyper 0.2.1 0.01
Table 2
これらは30種類のパイプライン(アライナー6種類バリアントコーラー5種類)で構成され11種類のツールが使用されています。 再現性を確保するためにソフトウェアのバージョンも記載しています。
2.3. Filtering

SRA (SRR098401) からRawデータを取得し、fastq-dumpで分割し、fastx toolkitでフィルタリングを行った。 具体的には、fastq-dumpは-split_filesと-split_spotフラグを使用し、fastq_quality_filterは以下の引数で実行しました。 -Q 33 -q 30 -p 50.

2.4. Aligning

Aligner は、sreads 引数を使用する場合を除き、デフォルトの引数を使用しました(利用可能な場合)。 使用したコマンドは以下の通りです。 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_Filtered.fastq alignments/NA12870…sam

(1)bwaaln-t10ゲノム.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 > alignment/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

<850>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.MOSAIKBuild -q raw_data/read_2_filtered.firstq -q2 raw_data/read_2_filtered.fastqmkb -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

Depth of Coverageを正しく計算するために、Picard ToolsのモジュールCalculateHsMetricsを以下の引数で使用した: (1)java -jar CalculateHsMetrics.jar I=NA12878.ALN.BQSR.bam O=ALN.O.log R=genome/hg19.O.log R=genome.BQSR.jarを使用した。fa TI=genome/truseq_exome.bed BI=genome/truseq_exome.bed VALIDATION_STRINGENCY=SILENT PER_TARGET_COVERAGE=ALN.ptc.bed このモジュールを機能させるには、TruSeq exome bedファイルにSAM alignment fileのヘッダーを前置しなければならないということが重要であり、このファイルは、そのためのものです。 さらに、6列目は4列目に移動させ、5列目はTruSeqベッドファイルから削除する必要があります。 アライメントファイルの処理

アライメントファイル(SAM/BAMファイル)の処理は、全てのアライナーにおいて、以下の手順が必要でした:(1)SAMtools viewによるSAM→BAM変換:(a) samtools view -bS alignments/NA12878.ALN.sam -o alignments/NA12878.ALN.bam(2)Picard ToolsモジュールSortSamによるBAMファイルのソート:(a)java -jar bin/SortSam.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.bam OUTPUT=alignments/NA12878.ALN.sorted.bam (2)BAM file sorting using the Picard Tools module of SortSam:(a)java -jar bin/SortSam.jarbam SORT_ORDER=coordinate(3)PCR Picard ToolsのモジュールMarkDuplicatesを用いた重複除去:(a)java -jar bin/MarkDuplicates.jar VALIDATION_STRINGENCY=SILENT I=alignments/NA12878.ALN.sorted.bam O=alignments/NA12878.ALN.dups_removed.Bamba (a)Picardツールズモジュールを用いた重複の除去(b)(b)MarkDuplicatesを用いた重複の除去。bam REMOVE_DUPLICATES=true M=alignments/metrics(4)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 O=ALIGNUMN.ALN.RG.Bam VALIDATION_STRINGENCY=SILENTbam SO=coordinate RGID=NA12878 RGLB=NA12878 RGPL=illumina RGPU=NA12878 RGSM=NA12878 CREATE_INDEX=true(5)GATK モジュール RealignerTargetCreator と 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.java -X:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.java -X:-DoEscapeAnalysis-jarjar -T IndelRealigner -R genome/hg19.fa -I alignments/NA12878.ALN.RG.bam -known genome/mills.vcf -o alignments/NA12878.ALN.indels.RG.bam -known genome/mills.vcf -o alignmentments/NA12878.ALN.indels.bam -maxReadsForRealignment 100000 -maxReadsInMemory 1000000 -targetIntervals tmp/ALN.intervals(6)GATK モジュール BaseRecalibrator と PrintReads によるベースリキャリブレーション:(a)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -T BaseRecalibrator -R genome/hg19.fa -I alignments/NA12878.ALN.indels.Bam -known genome/hg19.fa -I alignments/NA12878.ALN.indels.bam -targetIntervals tmp/ALN.intervals -maxReadsForRealignment 100000 -maxReadsInMemory1000000bam -knownSitesgenomedbsnp_137.hg19.excluding_sites_after_129.only_standard_chroms.vcf -o tmp/NA12878.ALN.grp(b)java -XX:-DoEscapeAnalysis -jar bin/GenomeAnalysisTK.jar -DoEscapeAnalysis.jar -XX:-DoEscapeAnalysis.jar -XX:-DoEscapeAnalysis.jar -XX:-DoEscapeAnalysis.jarjar -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. Variant calling

各Variant Callerの引数は、”threads “または “parallel “フラグが含まれていない場合はデフォルトのものを使用し、その場合も同様に使用しました。 また、indelは可能な限りSNVとは別に呼び出された。 具体的には、以下のようなコマンドを使用した。 FreeBayes

(1)freebayes -f genome/hg19.fa -i -X -u -v vcf_files/NA12878.ALIGNER.freebayes.raw.snv.vcf alignments/NA12878.ALIGNER.BQSR.Freebayes -f genome/hg19.fa -f freebayes.raw.snv.vacf 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.Bam -Na1878.ALIGNER.UB.raw.bam -Na1878.BQSR.bam -Na1878.ALIGNER.BQSR.bam -o vcf_files/NA12878.BQSR.bamvcf -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.BQSR.bam

(1.1)samtools Mpileup -uf genome/hg19.faalignments/NA12878.ALIGNER.BQSR.bamALIGNER.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.java (2)java (1)java (1)jar tools/NA12878.ALIGNER.BQSR.bamSNPSVM.raw.vcfアラインメントファイルに必要なCIGARフラグが存在しないため、SNPSVMはCUSHAW3のバリアントをコールできず、SAMtools mpileupは同じ理由でMOSAIKアラインメントのバリアントをコールできなかった。 また、SNPSVMはSNVのみを検出するため、このプログラムではindelは報告されていない。

2.8. Variant Filtration

使用するVariant Callerにより、Filtrationが異なる。 GATK HaplotypeCallerとGATK UnifiedGenotyperの場合、GATKモジュールのVariantRecalibratorとApplyRecalibrationを使用して、HapMap 3.3, Omni 2.5 SNP BeadChip, 1000 GenomeデータなしのdbSNP 137を学習セットとしてSNVをフィルタリングしていました。 SNPSVMでは、QUALスコア≥4、DP値≥6を使用しました。 FreeBayesとSAMtoolsでは、QUALスコア≧20、DP値≧6を使用した。 Variant比較

Variant比較では、全データセットに共通するSNVの比較にUSEQ 8.8.1を使用した。 indelの比較には、vcflibツールのvcfintersectを使用した。 SRR098401のシーケンスで使用したエクソームプルダウンキットで捕捉できる位置に比較を限定するために、2013年12月11日に取得したTruSeq hg19 exome bed file truseq_exome_targeted_regions.hg19.bed.chr を使用しました。 このファイルは、イルミナからこちらで入手できます。 http://support.illumina.com/sequencing/sequencing_kits/truseq_exome_enrichment_kit/downloads.ilmn. 異なるコールセット間でバリアントが同一に表現されるように、vcflibツールvcfallelicprimitivesを使用してvcfファイルの前処理を行った

2.10. 統計計算

True Positive (TP)。 テストするパイプラインで検出された変異で、NIST-GiaB listに存在する変異です。 テストされたパイプラインによって検出された変異ですが、NIST-GiaBリストには存在しない変異です。

True Negative (TN)(トゥルーネガティブ)。 テスト中のパイプラインで検出されなかった変異で、NIST-GiaBリストには存在しないものです。

False Negative (FN)(偽陰性)。 テストされたパイプラインでは検出されなかったが、NIST-GiaBリストには存在する変異。 Prefiltering Variants

バリアント解析を行う際、考慮しなければならない多くの落とし穴の1つが、エクソーム配列空間(エクソームキャプチャーキットによって定義される)とそれが解析結果にどのように影響するのかです。 今回のケースでは、Illumina TruSeq exomeキットを使用して抽出され、HiSeq 2000でシーケンスされた単一のexome(SRR098401)を持っていました。 このことを念頭に置き、バイオインフォマティクスツールの能力を測定するのであって、Illumina TruSeqエクソームキャプチャーキットがどれだけ機能したかを測定するのではないことを確認したかったのです。 つまり、プルダウンキットで定義されたエクソンに存在すると思われるバリアントのみをコールすることを試みたいのです。 このため、汎用のアノテーションベッドファイル、例えばhg19用のRefSeqではなく、Illuminaが提供するベッドファイルを使用します。 NIST-GiaBリストによると、この特定の個人については、TruSeqベッドファイルによって定義された領域内に合計34,886個のSNVと1,473個のインデルが存在するはずであることがわかりました

一度Illumina TruSeqエクソームベッドファイルによって定義された領域に位置しない変異をフィルタリングしたところ、パイプラインごとに数十万の推定変異(データなし)から平均約23000変異(SNVおよびインデル)へと変化しました(表3)。 これは、潜在的に興味深いバリアントの検索空間を大幅に削減するため、研究者が最初に行う重要なステップです。

440

Bowtie2

360

の3つ。

1.1,370

Aligner Genotyper Raw TP SNVs Raw FP SNVs Raw TP indels Raw FP indels
Bowtie2 FreeBayes 23,985 73,473 806 2,482
Bowtie2 GATK HC 21,631 273 771 1,103
Bowtie2 GATK UG 25.1 Bowtie2 UG 1,103 1,103136 2,276 420
Bowtie2 mpileup 21.1 420 2,276 420
440 440 1,030 734 1,414
Bowtie2 SNPSVM 17.0 1,030 1,040 1,414
47
BWA mem FreeBayes 23.1 FreeBayes 18,256 785 2,088
BWA mem GATK HC 21,707 367 779 1,348
BWA mem GATK UG 21.0 BWA mem GATK UG 1,348 1,348 213 402 408
BWA mem mpileup 25.1,081 2,129 761 1,772
BWA mem SNPSVM 17.0 BWA Mem 1,129 1,672 1,772> 65
BWA sampe FreeBayes 23,789 27,143 737 1,872
BWA sampe GATK HC 21,878 263 758 1,161
BWA sampe GATK UG 22,153 321 394 385
BWA sampe mpileup 25,206 2,205 684 1,401
BWA sampe SNPSVM 18,017 78
CUSHAW3 FreeBayes 23,191 53,525 624 3,310
cushaw3 gatk hc 19,673 14,814 751 4,727
cushaw3 gatk ug 19,113 13,184 1,005
CUSHAW3 mpileup 22,171 9,694 681 1,983
CUSHAW3 SNPSVM
MOSAIK FreeBayes 23,373 39,203 783 3,359
mosaik gatk hc 13,528 111 500 458
mosaik gatk ug 17,147 76 392 284
MOSAIK mpileup
mosaik snpsvm 14,586 8
Novoalign FreeBayes 22,794 2.970 678 1,554
Novoalign GATK HC 21,407 473 779
Novoalign GATK UG 21,113 144 387 365
Novoalign mpileup 24,512 1,861 773 1,781
Novoalign SNPSVM 17,109 164
表3
30パイプラインの生の変異統計(SNVとインデル含む)。
3.2. Raw Variant Results Compared to GiaB

この比較を行う際に理解したかったことの1つは、これらのツールで検出されたバリアントをフィルタリングすることの重要性です。 その理由は、理想的には、関心のある変異がフィルタリングの過程で失われないように、できるだけ高いレベルの感度を持ちたいと思うからです。 NIST-GiaB の結果と Bamshad らのレビューから、感度が問題になりうることは明らかであるため、このステップが必要かどうか、どの程度必要かを判断することは有益であると思われます。 ほとんどの場合、TP variantの数は私たちが予想した約2万個に近いかそれ以上ですが、一方で、FPの数が非常に多いケースもあります。 しかし、各ツールのアルゴリズム的な起源に由来すると思われる、数値の共通点を見つけることができます。 FreeBayesは、フィルタリングされていないバリアントの数とFPの数の両方を最も多く生成しています。 これは、ベイズ推定に基づく唯一のvariant callerではないものの、アラインメントの解釈においてユニークであるため、このツールからこのようなパフォーマンスを得ることができたと思われます。 つまり、これはハプロタイプベースのコーラーで、アライメントではなく、リードの配列そのものに基づいてバリアントを識別します。 Bowtie2、BWA mem、BWA sampeは、全体的に同じような結果になっています。 これは、これらのツールはすべて、指定されたタスクを実行する際に類似のアルゴリズムを利用するという事実によるものだろうと推測されます。 この観察は、MOSAIK (Smith-Waterman アルゴリズムを使用したギャップ アラインメント) および CUSHAW3 (ハイブリッド シーディング アプローチ) が両方とも非常に異なる基本アルゴリズムを持ち、その後非常に異なる結果を生成するという事実によってサポートされています。 バリアントコールャーのうち、サポートベクターマシンとモデル構築を利用して SNV コールを生成するのは、このツールだけです。 他の手法に比べて感度が低いという欠点がありますが、使用するアライナーに関係なく、非常に正確であるという利点があるように思われます。 4092>

インデルに関しては、このタイプの変異をうまく処理するものとして、他のどのアライナーも際立っているようには見えません。 実際、FP の数を見ると、インデルの識別の精度に最も大きな役割を果たすのはバリアントコーラーであることが明らかです。 また、CUSHAW3 + SNPSVM、MOSAIK + SAMtools mpileupのいずれのパイプラインのデータも、アライメントファイルにはバリアントコーラーが下流で機能するために必要なCIGAR文字列が含まれていないため、存在しないのです。 最後に、SNPSVMのindelデータがないのは、このツールがSNVの同定にのみ使用されているためである。 GiaBと比較したフィルタリング結果

表4のように、標準的なフィルタリングの実施により、各パイプラインで多くのFP SNVを除去することができた。しかし、これらのフィルタはほとんどの場合、SNV検出には少し強すぎるが、インデルには十分厳しくないようである。 このことは、各データセットで報告されたFPの数の違いを見れば明らかです。 例えば、Freebayesを用いたBowtie2では、72,570個のFP SNVが除去されましたが(98%の削減)、1,736個のFP indelしか除去されませんでした(70%の削減)。 また、使用したフィルターはパイプラインに依存し、ほとんどの場合、各パイプライン内でSNVとインデルのFPを同様に減少させたことに留意する必要があります。 例外は、CUSHAW3アライメントから同定されたバリアントの数が、他のアライメントと比較された場合です。 は、全体的にTP SNVの数が少なく、FP SNVの数が多く、フィルタリング後に1,000以上のFPインデルを生成する唯一のアライナーである。

1856 1856

BWA sampe

の3つ。

Aligner Genotyper Filtered TP SNVs Filtered FP SNVs Filtered TP indels Filtered FP indels
Bowtie2 FreeBayes 17,504 903 481 746
Bowtie2 GATK HC 17.1.2,330 29 648 687
Bowtie2 GATK UG 19,937 49 395 338
Bowtie2 mpileup 17.1.1,049 153 402 541
Bowtie2 SNPSVM 13.X.Y. Bowtie2 541
402 541 8
BWA mem FreeBayes 17,376 347 461 739
BWA mem GATK HC 19,388 302 689 860
BWA mem GATK UG 20,000 48 397 355
BWA mem mpileup 17,070 57 403 606
BWAのmem SNPSVM 15,060 10
BWA sampe FreeBayes 17,435 450 443 647
BWA sampe GATK HC 19,438 214 630 725
BWA sampe GATK UG 19.1,557 27 384 336
BWA sampe mpileup 17.1 BWA sampe 336 111 387 518
BWA sampe SNPSVM 15.0 BWA sampe SNPSVM BWA sampe 10
CUSHAW3 FreeBayes 16.1 FREEBAYS FREEBAYS 7,627 362 1,294
cushaw3 gatk hc 16,590 2,195 665 1,551
cushaw3 gatkug 17,939 2,202 357 545
CUSHAW3 mpileup 15,942 4,029 368 796
cushaw3 snpsvm
MOSAIK FreeBayes 17.177 679 458 645
mosaik gatk hc 11.X.Y. 1857 1858 679 645 1857616 33 426 255
mosaik gatk ug 16.0 1857> 1857> 1857.0423 42 381 224
MOSAIK mpileup
mosaik snpsvm 4,727 3
Novoalign FreeBayes 16.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X.X,658 219 384 559
Novoalign GATK HC 19,406 385 702
Novoalign GATK UG 20,521 46 386 315
Novoalign mpileup 16,493 62 396 579
Novoalign SNPSVM 14.1.2,451 18
表4
30パイプラインのフィルタリングバリアント統計(SNVとインデル含む)。

フィルタリングによってTP variantの数が大幅に減少することを考えると、CUSHAW3やFreeBayesを用いたパイプラインを除き、希少でインパクトの強いvariantを検索する場合はこのステップをスキップすることが賢明かもしれません。 その代わりに、統計学ではなく生物学に基づいたフィルタリング処理に時間をかけるとよいでしょう。 例えば、スプライスサイト変異、フレームシフトを引き起こすインデル、トランケーション変異、ストップロス変異、または関心のある表現型に生物学的に関連することが知られている遺伝子の変異など、影響が大きいと思われる小さな変異のリストを特定することに時間を費やす方がより理にかなっているかもしれません

3.4. 平均TPRと感度

表5に見られるように、CUSHAW3を除く各ツールのSNVのPositive Predictive Value(PPV)は91%~99.9%であるが、平均感度は非常に低い(50%程度)。 この相違はいくつかの理由が考えられますが、最も可能性が高いのは、エクソン間のカバレッジの深さが変化していることです。 SNV感度が低いことに加えて、インデルの感度も低い(約30%)ことがわかりますが、インデルのPPVはかなり低くなっています(35.86%から52.95%)。 これは、非常に短いインデルは従来のNGSでは検出が困難であること、異なるバリアントコーラーによるインデルの表現が、2つのインデルが異なると主張するツールを引き起こすこと、またはアライメントツールが同じインデルを異なる表現で作成することなどの理由によるものと考えられます。

32.79% 32.80% 32.80% 32.80%

50.77%

32.2%

50.76%

Tool Average SNV PPV Average SNV sensitivity Average indel PPV Average indel sensitivity
Bowtie2 98.69% 49.19% 45.45% 32.69%
BWA mem 99.9% 32.69% 32.49% 32.69% 32.49% 50.96% 43.24% 33.10%
BWA sampe 99.1% 99.09% 50.85% 45.31% 31.30%
cushaw3 80.9% 35.1% 35.2% 35.2%
48.08% 29.50% 29.74%
mosaik 98.51% 35.79% 52.95% 28.63%
Novoalign 99.17% 50.1% 50.18% 44.55% 31.70%
FreeBayes 90.95% 51.00% 35.0%
FreeBayes 31.70% 31.70%
32.79%
GATK HaplotypeCaller 97.05% 51.03% 43.17% 31.86% 32.79%
GATK UnifiedGenotyper 97.93% 52.12% 31.9% 31.1% 31.20% 32.1% 32.2%
SAMtools mpileup 94.99% 39.15% 31.30%
SNPSVM 99.92% 50.20%31.30% N/A

表5

ツールごとの平均PPV(陽性予測値)と感度を示した。

おそらく、両方のタイプの変異について最も可能性が高い説明は、深さの問題です。 どのバリアント解析研究でもそうですが、カバレッジの深さが増すと感度が上がりますが、エクソームキャプチャーキットがエクソンを均一にプルダウンできないため、良好なカバレッジの深さを保証することは不可能です 。 さらに、単一のエクソームキャプチャーキットがすべてのエクソンをカバーすることはない。 実際、全ゲノムシーケンスでは、エクソンよりも低い平均深度でバリアント解析を行った方が、その深度が均一であるため、良好な結果が得られることが示されている。 したがって、NIST-GiaBリストがエクソームおよびゲノムシーケンスデータのコンパイルから作成されたという事実によって、多数のバリアントが欠落している可能性があります。 最終的に、適切な感度を得るためには、最終的に全ゲノムシーケンスを行う必要があるが、それは現在ほとんどのラボにとってコスト的に不可能である。 幸い、このコストは下がり続けており、近い将来、エクソーム解析からより完全な全ゲノム解析へと徐々に移行することが予想される

3.5. 深さの関数としての感度

感度はツールの最も重要な性能指標の1つを反映し、ほとんどのツールは50%以上の感度を達成するのに苦労するので、我々はさらに深さがバリアントコール感度にどのように影響するかを調べたいと思います。 パイプライン、バリアントコーラー、アライナーについて、様々なツールの組み合わせを検討し、最適なものを決定しました。 Figure 2 では、Variant Caller と aligner の組み合わせのうち、感度と偽陽性率(FPR)により決定されるベスト 5 を採用しました。 つまり、FP SNVの数が最も少なく、TP SNVの数が最も多いものを選びました。 すぐに目につくのは、感度が予想より低いことです。 すべてのパイプラインがほぼ同じレベルで機能しており、約150倍の深さに到達するまでにほとんどのバリアントを同定しています。これは、この深さで十分である可能性が高く、バリアントの数が不足しているのは、特定のエクソンが平均よりも低いカバレッジであることが原因である可能性を示しています。 図3(b)に示すように、使用するアライナーに関係なく、その優れた性能を示しています。

図2
上位5パイプラインの深さの関数としての感度。 上位5つのパイプラインは、すべてのSNVの深さを感度に対してプロットしてここに示されている。
(a)
(a)
(b)
(b)
(a)
(a)(b)
(b)
図3
トップアライナーとトップバリアントコーラーにおける深さの関数としての感度。 (a) トップアライナーであるBWA memとすべてのvariant callerの組み合わせで、すべてのSNVの深さを感度に対してプロットした結果。 (b)全てのアライナーとペアになったトップバリアントコーラー、GATK UnifiedGenotyperの感度に対してプロットされた全てのSNVの深さの結果

トップ5のパイプラインを見ることに加えて、我々は、全てのアライナーとペアになったトップバリアントコーダ(図3(a))と同様に、全てのアライナーとペアになったベストバリアントコーダ(図3(b)で同じ解析を行うことは有用だと判断しました。) パイプラインと同様に、最良のアライナーは、最も多くのTP SNVを生成し、最も少ないFP SNVを生成するもの(この場合はBWA mem)として同定されました。 最適なアライメントを使用したにもかかわらず、バリアントコーラー間でかなり大きな差が見られますが、これは採用するアルゴリズムの違いに起因すると思われます(図3(a))。 しかし、最も性能の良いVariant CallerであるGATK UnifiedGenotyperの場合、上位4つのアライナーの間でばらつきが少なく、CUSAHW3とMOSAIKを除いて、ほとんどの状況でかなり良い性能であることを示しています。 上位パイプライン間の共有バリアント

最後に、異なるパイプライン間でバリアントコールセットがどのようにユニークであるかを知りたいと思いました。 これを行うために、再び上位5つのバリアントコールパイプラインに焦点を当てました。 Bowtie2 + UnifiedGenotyper、BWA mem + UnifiedGenotyper、BWA sampe + HaplotypeCaller、BWA sampe + UnifiedGenotyper、Novoalign + UnifiedGenotyperの5つのパイプラインです。 図4からわかるように、合計22,324個の異なるバリアントのうち、15,489個のSNV(70%)が共有されており、当該5種類のパイプラインの間には大きな重複があることがわかります。 しかし、これは5つのパイプラインのうち4つのパイプラインがバリアントコーラーとしてUnifiedGenotyperを使用していることが大きな原因であると言うこともできます。 この考え方は、パイプラインに固有のバリアント数が最も多い367個は、BWA sampeとHaplotypeCallerの組み合わせに属するという事実からも裏付けられます。 また、ユニークなSNVの2番目に多い数もBWA sampeアライナーに属していることは注目に値するので、ユニークなSNVの高い数は、バリアントコーラーよりもアライナーに起因するのが良い可能性がある。

図4
上位5パイプラインによって識別されるSNVの交差点。

4.結論

テストした30種類のパイプラインの中で、NovoalignとGATK UnifiedGenotyperはSNVのFP数を低く抑えながら最高の感度を示すことが分かりました。 アライナーの中では、BWA memが一貫して最高のパフォーマンスを示しましたが、使用するバリアントコーラーによって結果は大きく変わりました。 当然ながら、最も優れたバリアントコーラーであるGATK UnifiedGenotyperは、使用するアライナーに関わらず、ほぼ同様の結果を得ることができた。 しかし、どのパイプラインも平均感度33%以上、PPV53%以上を達成することができず、依然としてインデルの扱いが難しいことが読み取れます。 インデルの検出性能が低いことに加え、変異の種類に関わらず、感度は本論文で概説したすべてのパイプラインで問題になっています。 NA12878のエクソームで予想されるSNVの数は34,886ですが、上位5つのパイプラインで特定されたすべてのバリアントのユニオンを使用しても、特定された最大の数は非常に少なく(22,324)、その結果、NA12878のエクソームでは、SNVの数は減少しました。 4092>

情報公開

Adam CornishはChittibabu Gudaの研究室の大学院生で、コンピュータ科学とゲノミクスのトレーニングを受けています。 Chittibabu Guda(准教授)は、分子生物学と計算生物学の学際的なバックグラウンドを持っています。

利益相反

著者らは、利益相反を認識していない。

著者貢献

アダム・コーニッシュは、研究の設計、すべての分析の実施、図の作成、および論文の執筆を担当した。 Chittibabu Gudaは、論文の改善に関する重要なフィードバックと分析自体に関する意見を提供し、論文を徹底的に編集した。 4092>

謝辞

この研究は、ネブラスカ大学医療センター(UNMC)からChittibabu Guda氏への開発資金により支援されました。 Novoalignの試用版を提供してくださった制作者の方々、本研究を円滑に進めるためのインフラを支援してくださったUNMCのBioinformatics and Systems Biology Core facilityに感謝します