和訳：GATKによるジェノタイプコールの非モデル生物におけるベストプラクティス
2022-06-02
2022-06-10 最終更新
GATK: the best practice for genotype calling in a non-model organism
非モデル生物での遺伝子型判定には、GATKのベストプラクティスの修正が不可欠な場合があります。この投稿では、この問題に対する私のアプローチを紹介します。
GATK (Genome Analysis Toolkit) は、様々な生物のハイスループット・シーケンス・データにおけるジェノタイプコールに最もよく使用されているソフトウェアです。GATKのベストプラクティスは、SAM、BAM、CRAM、VCF形式のシーケンスデータの様々な解析のための素晴らしいガイドとなります。しかし、GATKは主にヒトの遺伝子データを解析するために設計されており、そのパイプラインはすべてこの目的に最適化されています。ヒト以外のデータに対して同じパイプラインをそのまま使用すると、何らかの不正確な結果を招く可能性があります。特に、参照ゲノムが解析対象サンプルと同種でない場合に問題となります。 ここでは、Capsella bursa-pastorisという非モデル生物で、参照ゲノムが姉妹種しかない全ゲノムシーケンスデータに対して、GATKパイプラインで遺伝子型判定を行った結果を説明します。特定の研究事例ではありますが、私の修正点を説明することで、他の研究者がこのパイプラインを非モデル生物に採用する際に役立つと信じています。
更新：このパイプラインはより自動化することができます。詳しくはこの投稿をご覧ください。 Check sequences quality
解析の前に、配列データの品質をチェックし、何も問題がないことを確認するのは良い習慣です。この目的のためには、FASTQCを使うことをお勧めします。グラフィカルなインターフェースもありますが、私はコマンドラインを使うのが好きです。 code: shell
fastqc input.fq.gz
https://scrapbox.io/files/62982528a912ae001dc24070.png
出力はHTMLファイルで、様々な品質メトリクスがグラフィカルに表示されます。出力を説明するビデオもあります。データに大きな問題がないことを確認してください。研究者の中には、Trimomaticのようなツールを使って、低品質の配列尾部をトリミングする人もいます。昔はこのステップも有効だったと思いますが、配列品質の現状では、もう必要ないと思います。また、トリミングによって偽のバリアントが発生するという証拠もあります。 Map to the reference genome
GATKベストプラクティスでは、DNAデータにはBWAアライナーを使用することを推奨しています。BWAは、リファレンスとの乖離が少ない配列をマッピングするための、精度と速度の点でおそらく最良のソフトウェアです。非モデル生物では、一般的に参照ゲノムが姉妹種しかなく、その姉妹種がかなり乖離している可能性があります。そのため、この乖離を考慮したアライナーを使用する必要があります。 大貫注：もう非モデル生物でも十分高品質な参照ゲノム（染色体スケールの一歩手前くらいとか）をそこそこ簡単に作れる場合はあると思うので，もしそうならBWAを使っておけばよさそう。
Heliconiusの多様な配列リードをマッピングするための様々なアライナーを比較した素晴らしいブログポストがあります。この記事の著者である私も、Stampyを使うことを好んでいます。Stampyはデフォルトの設定でもそれなりに良い結果が得られますが、-substitutionrateオプションで参照元との乖離を正しく指定することで、精度を向上させることができます。 まず、参照ゲノムを準備します。
code: shell
stampy.py -G REF reference_genome.fasta.gz && stampy.py -g REF -H REF
次に、乖離を考慮した参照ゲノムにリードをマップします（今回は0.025）。
code:shell
stampy.py -t8 -g REF -h REF --substitutionrate=0.025 -o output.sam -M R1.fq.gz R2.fq.gz
研究対象種が、最も近い参照ゲノムから乖離しすぎている場合、de novoゲノムアセンブリの方が良い場合があります。ただし、乖離した参照ゲノムにマッピングする方が、ラフなde novo アセンブリー（ADD reference）よりも良い結果が得られる場合があります。
Prepare BAM files
GATKでジェノタイピングするためのBAMファイルを準備します。
SAMからBAMへの変換
上記のステップで、SAM形式のアライメントが生成されます。SAMファイルはタブ区切りのテキストファイルで、通常は非常に大きなサイズになります。スペースを節約し、ソフトウェアが大きなアライメントを扱いやすくするために、通常、SAMはバイナリ形式のBAMに変換されます。Picard Toolsは、ファイル内のアライメントされたリードを同時にソートすることもできるため、変換にお勧めのソフトウェアです。 code:shell
java -Xmx8g -jar picard.jar SortSam \
INPUT=aligned_reads.sam \
OUTPUT=sorted_reads.bam \
SORT_ORDER=coordinate
Picard Toolsでは変換に失敗する場合があります。その場合は、SAMtoolsを使うのも一つの方法です。
code: shell
samtools view -bS -@ 8 SE14_stampy0.025_DNA.sam > SE14_stampy0.025_DNA.bam
samtools sort -@ 8 SE14_stampy0.025_DNA.bam SE14_stampy0.025_DNA_sorted
samtools index SE14_stampy0.025_DNA_sorted.bam
大貫注：samtoolsに直接samをインプットして，samtools sortにbamでアウトプットさせることもできる。
BAMファイル取得後、SAMファイルを削除することができる。注意! ファイルを削除する前に、すべての変換が成功したことを確認してください。成功していなかった場合、マッピングのステップを繰り返すことになり、非常に時間がかかります。
マッピングクオリティのチェック
マッピングの品質はSAMファイルでも評価できますが、BAMファイルの方が容量を消費しないので、私はBAMファイルを解析しています。この解析には、Qualimapを好んで使用しています。Qualimapは非常に高速で、必要なグラフを含む非常に詳細なサマリーを作成することができます。 code: shell
qualimap bamqc -nt 8 -bam file.bam -outdir results_folder
https://scrapbox.io/files/6298252dac3c50001dc978d6.png
参照ゲノム全体のカバレッジ、ヌクレオチド含有量、マッピング品質の分布は、問題のある領域を特定するのに役立ちます。例えば上の図では、セントロメア近傍の領域はゲノムの他の部分よりも100倍高いカバレッジを持っているので、これらの領域で見つかったことの解釈には注意が必要です。
重複のマーク
PCR重複の可能性があるものはPicard Toolsでマークする必要があります。
code: shell
java -Xmx8g -jar picard.jar MarkDuplicates \
INPUT=sorted_reads.bam \
OUTPUT=sorted_reads_duplMarked.bam \
METRICS_FILE=sorted_reads_duplMarked.metrics \
MAX_FILE_HANDLES=15000
MAX_FILE_HANDLESはulimit -nコマンドの出力より少し小さめにしてください。
重複したリードは削除されず、必要に応じてその後の解析に含めることができる（GATKオプション：-drf DuplicateRead）ので、PCRなしのライブラリ作成手順を使用した場合でも、重複をマークすることは意味があります。
グループの追加
GATKでは、BAMファイル中にリードグループの情報が必要です。これは、サンプルを区別したり、シーケンス技術に関連するアーティファクトを検出したりするために使用されます。リードグループを追加するにはPicard Toolsを使用します。 code: shell
java -Xmx4g -jar picard.jar AddOrReplaceReadGroups \
I=sorted_reads_duplMarked.bam \
O=sorted_reads_duplMarked_readgroup.bam \
RGLB=library \
RGPL=illumina \
RGPU=barcode \
RGSM=sample
出来上がったBAMファイルにインデックスを付けます。
code: shell
samtools index SE14_stampy025_DNA_sorted_DuplMarked_readgroup.bam
インデルのリアライン
Realigner Targetを作成します（-ntオプションでマルチスレッドに対応します）。
code: shell
java -Xmx4g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T RealignerTargetCreator \
-nt 16 \
-R ~/reference/Crubella_183.fa \
-I sorted_reads_duplMarked_readgroup.bam \
-o sorted_reads_duplMarked_readgroup.intervals \
-log sorted_reads_duplMarked_readgroup.intervals.log
リアラインメントを行います（マルチスレッディングには対応していません）。
code: shell
java -Xmx4g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T IndelRealigner \
-R reference_genome.fasta \
-I sorted_reads_duplMarked_readgroup.bam \
-targetIntervals sorted_reads_duplMarked_readgroup.intervals \
-o sorted_reads_duplMarked_readgroup_realigned.bam
Base Quality Scoreの再キャリブレーション
GATKのベストプラクティスでは、Base Quality Scoreの再キャリブレーションの実行を推奨しています。この手順は、参照用トレーニングデータセットと比較することで、データの系統的なエラーを検出するものです。非モデル生物の問題は、通常、トレーニングに必要な大規模で信頼性の高い遺伝子型データベースが存在しないことです。 この問題の解決策の1つは、解析したデータと同じデータから、非常に厳しいフィルタリング基準を用いてカスタムのトレーニングデータベースを作成し、データを再キャリブレーションすることです。私はこの問題をGATKのフォーラムで取り上げ、私のアプローチはGATKチームによってサポートされましたが、実際にはうまくいきませんでした。以下のように、単に分布が低いスコアにシフトしてしまうのです。 https://scrapbox.io/files/629825310187e10023519b67.png
というわけで、大規模な学習データセットがない限り、このステップはスキップする必要がありそうです。
Genotype
BAMファイルの準備ができたので、GATKで処理し、遺伝子型データを取得することができます。
GVCFファイルの生成
各ファイルに対して、GVCFモードでHaplotypeCallerを実行します。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-nct 16 \
-R reference_genome.fasta \
-I sorted_reads_duplMarked_readgroup_realigned.bam \
-ERC GVCF \
-hets 0.015 \
-indelHeterozygosity 0.01 \
-o sorted_reads_duplMarked_readgroup_realigned.g.vcf
研究対象生物がヒト以外の場合、ヘテロ接合度の調整が必要な場合があります。デフォルトは0.001です。ヘテロ接合度が異なる場合は、-hets および -indelHeterozygosity で指定してください。過去に作成したデータから推定することができます。
近似ヘテロ接合度 = ((インデルのあるサイト数)/(全サイト数))*(1-(部位ごとの平均ホモ接合型数)/(全個体数))
このコードで部位ごとの平均ホモ接合型数を計算することができます。
code: shell
awk -F\0\/0: '!/^ *#/ {total += NF-1; count++} END { print total/count }' Indels.vcf
awk -F\0\/0: '!/^ *#/ {total += NF-1; count++} END { print total/count }' SNPs.vcf
HaplotypeCallerのアセンブリの確認
これはオプションのステップです。HaplotypeCallerは、リードの局所的な再アセンブリを行い、その後サイトのジェノタイピングを行います。このローカルに再集合されたフラグメントから BAM ファイルを作成し、IGV で品質を評価することができます。 code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-R reference_genome.fasta \
-I sorted_reads_duplMarked_readgroup_realigned.bam \
-ERC GVCF \
-hets 0.015 \
-indelHeterozygosity 0.01 \
-L scaffold_1 \
-bamout sorted_reads_duplMarked_readgroup_realigned_scaffold_1.g.vcf.bam \
-o sorted_reads_duplMarked_readgroup_realigned_scaffold_1.g.vcf
https://scrapbox.io/files/629825341b641c001dff4705.png
IGV: Original mapping (above) and HaplotypeCaller re-assembly (below)
全GVCFファイルのジェノタイピング
全てのGVCFファイルに対してGATKを実行します。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T GenotypeGVCFs \
-nt 16 \
-R reference_genome.fasta \
-V 01_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 02_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 03_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 04_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 05_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 06_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 07_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 08_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 09_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-V 10_sorted_reads_duplMarked_readgroup_realigned.g.vcf \
-hets 0.015 \
-indelHeterozygosity 0.01 \
-stand_call_conf 30.0 \
-allSites \
-o GVCFall.vcf
この際、必要に応じて-stand_call_confを調整したり、他のフィルターを追加したりしてください。
また、私は通常、-allSiteオプションですべてのサイト（すなわち、非変動および多型）をジェノタイピングし、後で変動サイトを選択します。集団遺伝学の統計量の推定には、サイト数の合計でスケーリングする必要があるためです（つまり、非可変を含む）。
SNPsとIndelsの選択
SNPを選択する。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference_genome.fasta \
-V GVCFall.vcf \
-selectType SNP \
-o GVCFall_SNPs.vcf
Indelを選択する。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference_genome.fasta \
-V GVCFall.vcf \
-selectType INDEL \
-o GVCFall_INDELs.vcf
Filter Variants
まずバリアントサイトのフィルタリングを行い、次に個々のGenotypeのフィルタリングを行っています。
VariantsフィルタリングはSNPs/Indels VCFファイルのみに適用されますが、GenotypeフィルタリングはSNPs/Indels VCFと全ゲノムVCFの両方に適用することができます。
Variant Quality Scoreの抽出
SNPの場合。
code:shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall_SNPs.vcf \
-F CHROM -F POS -F QUAL -F QD -F DP -F MQ -F MQRankSum -F FS -F ReadPosRankSum -F SOR \
--allowMissingData \
-o GVCFall_SNPs.table
Indelの場合。
code: shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall_INDELs.vcf \
-F CHROM -F POS -F QUAL -F QD -F DP -F MQ -F MQRankSum -F FS -F ReadPosRankSum -F SOR \
--allowMissingData \
-O GVCFall_INDELs.table
大貫注：-oは-Oのミスと思われ，訳では修正している。また，2022-06-10現在，--allowMissingDataはなくなっているよう。
これらのスコアは、分布の構築やフィルタリングカットオフの定義に使用されます。
Variants Scoreの診断プロットの作成
以下のRコードを実行します。
code: R
library('gridExtra')
library('ggplot2')
VCFsnps <- read.csv('GVCFall_SNPs.table', header = T, na.strings=c("","NA"), sep = "\t")
VCFindel <- read.csv('GVCFall_INDELs.table', header = T, na.strings=c("","NA"), sep = "\t")
dim(VCFsnps)
dim(VCFindel)
VCF <- rbind(VCFsnps, VCFindel)
VCF$Variant <- factor(c(rep("SNPs", dim(VCFsnps)1), rep("Indels", dim(VCFindel)1))) snps <- '#A9E2E4'
indels <- '#F4CCCA'
DP <- ggplot(VCF, aes(x=DP, fill=Variant)) + geom_density(alpha=0.3) +
geom_vline(xintercept=c(10,6200))
QD <- ggplot(VCF, aes(x=QD, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=2, size=0.7)
FS <- ggplot(VCF, aes(x=FS, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=c(60, 200), size=0.7) + ylim(0,0.1)
MQ <- ggplot(VCF, aes(x=MQ, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=40, size=0.7)
MQRankSum <- ggplot(VCF, aes(x=MQRankSum, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=-20, size=0.7)
SOR <- ggplot(VCF, aes(x=SOR, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=c(4, 10), size=1, colour = c(snps,indels))
ReadPosRankSum <- ggplot(VCF, aes(x=ReadPosRankSum, fill=Variant)) + geom_density(alpha=.3) +
geom_vline(xintercept=c(-10,10,-20,20), size=1, colour = c(snps,snps,indels,indels)) + xlim(-30, 30)
svg("Co_10accessions_FromStephen.svg", height=20, width=15)
theme_set(theme_gray(base_size = 18))
grid.arrange(QD, DP, FS, MQ, MQRankSum, SOR, ReadPosRankSum, nrow=4)
dev.off()
以下のようなプロットが生成される。
https://scrapbox.io/files/62982539f3a4a9001d170911.png
Annotations:
DP - combined depth per SNP across samples (24 samples in the case above)
QD - variant confidence standardized by depth
MQ - Mapping quality of a SNP
FS - strand bias in support for REF vs ALT allele calls
SOR - sequencing bias in which one DNA strand is favored over the other
MQRankSum - Rank sum test for mapping qualities of REF vs. ALT reads
ReadPosRankSum - do all the reads support a SNP call tend to be near the end of a read
Variantフィルタリングの適用
GATKのベストプラクティスに沿ってデータをフィルタリングし、私のデータ用に若干の調整を加える。例えば、私の場合、リードが予想以上にリファレンスから乖離しているため、MQRankSumフィルターを緩和しています。また、特定のスコアの分布に基づいてカットオフを定義し、できるだけ多くのデータを保持しつつ、最も信頼性の低いサイトを削除しています。
SNPをフィルタリングする。
code: shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R reference_genome.fasta \
-V GVCFall_SNPs.vcf \
--filterExpression "QUAL < 0 || MQ < 30.00 || SOR > 4.000 || QD < 2.00 || FS > 60.000 || MQRankSum < -20.000 || ReadPosRankSum < -10.000 || ReadPosRankSum > 10.000" \ # you define what to remove here
--filterName "my_snp_filter" \
-o GVCFall_SNPs_filter.vcf
grep -E '^#|PASS' GVCFall_SNPs_filter.vcf > GVCFall_SNPs_filterPASSED.vcf
indelをフィルタリングする。
code: shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R reference_genome.fasta \
-V GVCFall_INDELs.vcf \
--filterExpression "QUAL < 0 || MQ < 30.00 || SOR > 10.000 || QD < 2.00 || FS > 200.000 || ReadPosRankSum < -20.000 || ReadPosRankSum > 20.000" \ # you define what to remove here
--filterName "my_indel_filter" \
-o GVCFall_INDELs_filter.vcf
grep -E '^#|PASS' GVCFall_INDELs_filter.vcf > GVCFall_INDELs_filterPASSED.vcf
Variantsフィルタリング結果確認
何サイトが削除されたかを見ます。
フィルタリング前のサイト数は
code: shell
grep -vc "^#" GVCFall_INDELs.vcf
grep -vc "^#" GVCFall_SNPs.vcf
フィルタリング後に保持されているサイト数は
code: shell
grep -vc "^#" GVCFall_SNPs_filterPASSED.vcf
grep -vc "^#" GVCFall_INDELs_filterPASSED.vcf
また、失敗したフィルタがあるかどうかの確認も重要です。フィルタはよく失敗します。上記の「Variant Quality Scoresの抽出」と同様に、フィルタリングされたデータのスコアテーブルを生成します。このRコードを実行すると、すべてのフィルターで0が得られるはずです。
code: R
# SNPs
VCFsnps <- read.csv('GVCFall_SNPs_filterPASSED.table', header = T, na.strings=c("","NA"), sep = "\t")
head(VCFsnps)
sum(na.omit(VCFsnps$QD) < 2)
sum(na.omit(VCFsnps$FS) > 60)
sum(na.omit(VCFsnps$MQ) < 40)
sum(na.omit(VCFsnps$MQRankSum) < -20)
sum(na.omit(VCFsnps$SOR) > 4)
sum(na.omit(VCFsnps$ReadPosRankSum) < -10)
sum(na.omit(VCFsnps$ReadPosRankSum) > 10)
# Indels
VCFindel <- read.csv('GVCFall_INDELs_filterPASSED.table', header = T, na.strings=c("","NA"), sep = "\t")
head(VCFindel30)
sum(na.omit(VCFindel$QD) < 2)
sum(na.omit(VCFindel$FS) > 200)
sum(na.omit(VCFindel$MQ) < 40)
sum(na.omit(VCFindel$SOR) > 10)
sum(na.omit(VCFindel$ReadPosRankSum) < -20)
sum(na.omit(VCFindel$ReadPosRankSum) > 20)
いくつかのフィルタが失敗した場合は、-filterExpressionで順序を変えてみてください。
Filter genotypes
信頼度の低いバリアントサイトがすべて取り除かれたら、VCFファイルをgenotype qualityのためにfilterします。GATKはgenotypeのアノテーションスコアをいくつか提供していますが、それらはvariantサイトとnon-variantサイトで同じように使うことはできません（例えばSNPsではGQ、non-variantサイトではRGQなど）。このようなスコアを同時にフィルタリングすると、variant/non-variant sitesの比率にバイアスがかかるため、合理的ではありません。そこで、通常、depthによるフィルタリングを行います（GenotypeGVCFsの呼び出し時に、Genotype Quality Filter 30も適用されていることに注意してください）。
depth情報の抽出
全ゲノムVCFに対してのみ、深さ情報を抽出します。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall.vcf \
-F CHROM -F POS -GF GT -GF DP \
-O GVCFall.DP.table
大貫注：-oは-Oのミスと思われ，訳では修正している。
抽出した情報の可視化
GVCFall.DP.tableは通常Rでは大きすぎるので、サンプルごとに分割し、ジェノタイピングされたポジションのみを残します（!="./.")。
code: shell
for ((i=3; i<=21; i +=2)); do cut -f $i,$((i+1)) GVCFall.DP.table | awk '$1 != "./." {print $2}' > $i.DP; done
ここで、3-21は10サンプルの奇数番号です。各サンプルは2列(GT, DP)で表されるため、この番号付けが必要です。
DP分布をプロットし、カットオフを定義（サンプル数に合わせて変更）。
code: R
# define the file names list (10 samples here)
nameList <- c()
for (i in 3:21) { # 21 - odd number for 10 samples
if (i %% 2 ==1) nameList <- append(nameList, paste(i, ".DP"))
}
qlist <- matrix(nrow = 10, ncol = 3) # define number of samples (10 samples here)
qlist <- data.frame(qlist, row.names=nameList)
colnames(qlist)<-c('5%', '10%', '99%')
jpeg("GVCFall.DP.jpeg", height=1600, width=1200)
par(mar=c(5, 3, 3, 2), cex=1.5, mfrow=c(8,4)) # define number of plots for your sample
for (i in 1:31) {
DP <- read.table(nameListi, header = T) qlisti, <- quantile(DP,1, c(.05, .1, .99), na.rm=T) d <- density(DP,1, from=0, to=100, bw=1, na.rm =T) plot(d, xlim = c(0,100), main=nameListi, col="blue", xlab = dim(DP)1, lwd=2) abline(v=qlisti,c(1,3), col='red', lwd=3) }
dev.off()
write.table(qlist, "GVCFall.DP.percentiles.txt")
https://scrapbox.io/files/629825413a26e9001d6fb678.png
Example of the DP distributions
GVCFall.DP.percentiles.txt の内容を確認し、許容できるカットオフを選択することができます。私は通常、5パーセンタイル以下と99パーセンタイル以上のジェノタイプを破棄しています。
depthフィルタの適用
フィルタリングされたgenotypeにマークをつけます。
code: shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R reference_genome.fasta \
-V GVCFall.vcf \
-G_filter "DP < 6 || DP > 100" \
-G_filterName "DP_6-100" \
-o GVCFall_DPfilter.vcf
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R reference_genome.fasta \
-V GVCFall_SNPs_filterPASSED \
-G_filter "DP < 6 || DP > 100" \
-G_filterName "DP_6-100" \
-o GVCFall_SNPs_filterPASSED_DPfilter.vcf
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R reference_genome.fasta \
-V GVCFall_Indels_filterPASSED.vcf \
-G_filter "DP < 6 || DP > 100" \
-G_filterName "DP_6-100" \
-o GVCFall_Indels_filterPASSED_DPfilter.vcf
フィルタリングされたサイトをノーコールに設定します。
code: shell
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference_genome.fasta \
-V GVCFall_DPfilter.vcf \
--setFilteredGtToNocall \
-o GVCFall_DPfilterNoCall.vcf
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference_genome.fasta \
-V GVCFall_SNPs_filterPASSED_DPfilter.vcf \
--setFilteredGtToNocall \
-o GVCFall_SNPs_filterPASSED_DPfilterNoCall.vcf
java -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T SelectVariants \
-R reference_genome.fasta \
-V GVCFall_Indels_filterPASSED_DPfilter.vcf \
--setFilteredGtToNocall \
-o GVCFall_Indels_filterPASSED_DPfilterNoCall.vcf
VCF to Tab
全てのフィルターが適用されたら、アノテーション情報はもう必要ないので、信頼度の高いジェノタイプだけを残しておきます。
code: shell
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall_DPfilterNoCall.vcf \
-F CHROM -F POS -GF GT \
-o whole_Genome.table
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall_SNPs_filterPASSED_DPfilterNoCall.vcf \
-F CHROM -F POS -GF GT \
-o SNPs.table
java -Xmx8g -jar ~/Programs/GATK/GenomeAnalysisTK.jar \
-T VariantsToTable \
-R reference_genome.fasta \
-V GVCFall_Indels_filterPASSED_DPfilterNoCall.vcf \
-F CHROM -F POS -GF GT \
-o Indels.table
また、私は通常、2文字コード化された遺伝子型を1文字コード化した遺伝子型に変換しています。詳しくはこちら。 code: shell
python vcfTab_to_callsTab.py -i whole_Genome.table -o whole_Genome.tab
python vcfTab_to_callsTab.py -i SNPs.table -o SNPs.tab
全ゲノムファイルとSNPファイルをマージ（オプション）
解析内容によっては、SNPファイルと全ゲノムファイルをマージすることも可能です。この2つのファイルは、SNPファイルに対してより厳格に、それぞれ異なるフィルタリングが行われています。そのため、全ゲノムファイルに含まれる一部の遺伝子型は、偽陽性である可能性が高くなります。このような遺伝子型を排除するには、SNPsファイルに存在しない多型部位を全ゲノムファイルから削除すればよいでしょう。これは、私の merge_SNP_wholeGenome_TabFiles.py スクリプトで行うことができます。なお、2文字コード化された遺伝子型を1文字コード化した遺伝子型に変換する必要があるのは、先ほど説明したとおりです。 code: shell
python merge_SNP_wholeGenome_TabFiles.py -g whole_Genome.tab -s SNPs.tab -o merged.tab
What is next?
さて、これらの遺伝子型テーブルを他の形式に変換したり、直接解析するためのカスタムスクリプトを書くのはとても簡単です。
ご質問やご提案がありましたら、お気軽にメール（dmytro.kryvokhyzha@evobio.eu）してください。
Written on September 22, 2016