|
% read aligned (filtered)
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メモリー89
BWA sampe |
85.95 |
97.49 |
46.67 |
CUSHAW3 |
85.85.85.00 |
99.81 |
47.69
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つである。
|
Tool |
Type |
Version |
Reference |
|
Bowtie2 |
Aligner |
2.X.S.A. |
Bowtie2 |
Aligner |
Version |
|
BWA sampe |
Aligner |
0.7.5a
|
BWA mem |
Aligner |
0.7.5a |
|
CUSHAW3 |
Aligner |
3.0.3 |
|
MOSAIK |
Aligner |
2.2.3
|
Novoalign |
Aligner |
3.3.0.02.07 |
N/A |
FreeBayes |
Genotyper |
v9.0 |
FreeBayes |
Genotyper50579.2-19-g011561f |
|
GATK HaplotypeCaller |
Genotyper |
2.7-2
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 |
0.1.19
|
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)。 これは、潜在的に興味深いバリアントの検索空間を大幅に削減するため、研究者が最初に行う重要なステップです。
|
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 |
440 |
1,030 |
734 |
1,414 |
Bowtie2 |
SNPSVM |
17.0 |
1,030 |
1,040 |
1,414 |
Bowtie2
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 |
360
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 |
の3つ。
– |
– |
mosaik |
snpsvm |
14,586 |
8 |
– |
Novoalign |
FreeBayes |
22,794 |
2.970 |
678 |
1,554 |
Novoalign |
GATK HC |
21,407 |
473 |
779 |
1.1,370
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インデルを生成する唯一のアライナーである。
|
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 |
1856 1856
111 |
387 |
518 |
BWA sampe |
SNPSVM |
15.0 |
BWA sampe |
SNPSVM |
BWA sampe
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 |
の3つ。
– |
– |
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つのインデルが異なると主張するツールを引き起こすこと、またはアライメントツールが同じインデルを異なる表現で作成することなどの理由によるものと考えられます。
|
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% |
32.79% 32.80% 32.80% 32.80%
GATK UnifiedGenotyper |
97.93% |
50.77%
52.12% |
31.9% |
31.1% |
31.20% |
32.1% |
32.2% |
32.2%
SAMtools mpileup |
94.99% |
50.76%
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)
(b)
(a) (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に感謝します
。 |