2019 11/19 コマンドエラー修正（A、Bの比較なのにCの表記がある）
2020 10/20 インストール追記
以前ショートリードからindelとSNVを検出するワークフローを紹介した。
複数サンプルがある場合、上記のような方法でVCFファイルを出力した後、サンプル間で共通のSNPs、サンプルごとに固有のSNPsなどを絞り込む必要が出てくるシチュエーションは多いと思われる。やり方は複数存在するが、今回は２つ紹介する。
9/14追記
A、gatkを使って行う
SelectVariants
以下の２ステップで固有/共通変異を識別する。
１、比較するvcfファイルを結合（２サンプル）
java -jar GenomeAnalysisTK.jar -T CombineVariants -R ref.fa -genotypeMergeOptions UNIQUIFY -V sample1.vcf -V sample2.vcf -o union.vcf
-V: vcfファイルの指定。
genotypeMergeOptions: サンプルが同じ場合に区別するために必要。INFOカラムにsuffixがつく。詳細は公式HP参考。
(Determines how we should merge genotype records for samples shared across the ROD files (UNIQUIFY| PRIORITIZE|UNSORTED|REQUIRE_UNIQUE)
-Vを増やすことで３つ以上のサンプルに対応可。
brewでgatkを入れている人はjava -jar GenomeAnalysisTK.jarの部分をgatkだけに変えて上記コマンドを打つ。
２、共通する変異の抽出
java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -o intersect.vcf
３、固有の変異（= 全変異 ー 共通する変異）
java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -invertSelect -o invertselect.vcf
-invertSelect: 共通しない変異
INFOの右端にサンプル名が付いており、それでどちらのサンプルか判別可能になっている。この情報を指標に簡単なスクリプトで分離することもできる。
１のステップで結合するVCFを追加することで３つ以上のデータの共通、固有変異の抽出も可能であるし、VCFフォーマットで記述されてるデータなら、他のindel toolで検出した変異をマージすることも簡単にできると思われる。
追記; gatk、freebayes、svabaはmegeできることを確認したが、platypusやFermikitなど一部のツールではVCFのフォーマットに独自の部分があり、mergeできなかった。
B、bcftoolsを使って行う
bcftoolsのisecを使う。やり方は公式ページのisecの説明部分に書かれている。
bcftoolsを使うには、vcfファイルが圧縮され、indexがつけられている必要がある。
準備１ bgzipで圧縮する。
bgzip -c sampleA.vcf > sampleA.vcf.gz
bgzip -c sampleB.vcf > sampleB.vcf.gz
bgzip -c sampleC.vcf > sampleC.vcf.gz
-cをつけて実行することで、元のファイルを保持したままsample.vcf.gzができる。
準備２ インデックスファイル作成
bcftools index sampleA.vcf.gz
bcftools index sampleB.vcf.gz
bcftools index sampleC.vcf.gz
これで準備ができた。いかに４つの例を紹介する。
１、AとBの比較
bcftools isec sampleA.vcf.gz sampleB.vcf.gz -p output
-p: 出力ディレクトリ（必須）
outputディレクトリにRWADMEと４つのファイルができる。仕分け条件はREADMEに書いてあるが、
・１つ目: sampleAに固有の変異（アレル）
・２つ目: sampleBに固有の変異（アレル）
・３つ目: sapleAから探したAとB共通の変異（アレル）
・４つ目: sapleBから探したAとB共通の変異（アレル）
３と４はどちらのVCFをベースに作るかの違いで、検出される部位の数は完全に同じである。VCFの由来かは４つのVCFファイルそれぞれのINFOフィールドの右端に記載されている。
２、AとBとCを比較してA特異的な変異
bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -C
-C: output positions present only in the first file but missing in the others
-Cをつけると、最初のsampleA.vcf特異的な変異だけが出力される。
３、AとBとC共通の変異
bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -n=3
-n: output positions present in this many (=), this many or more (+), this many or fewer (-), or the exact same (~) files
-nをつけることで比較するVCFファイルの数を調節する。
-n=3 → ３つのVCFファイル
-n+3 → ３つ以上のVCFファイル (３を含む）
-n-3 → ３つ以下のVCFファイル(３を含む）
上記の場合、-n=3でも-n+3でも同じ結果になる。
４、A、B、C固有の変異
bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -n-1
-n-1をつけることで、sampleA.vcf、sampleB.vcf、sampleC.vcf固有の変異が３つのファイルそれぞれに出力される。
それ以外にもたくさんオプションはある。
・出力する変異の種類を調節するには、-cオプションを使う。
-c all: 全変異
-c snps: 塩基置換のみ
-c indels: indelのみ
・出力ファイルを制限するには-wを使う。例えばsampleAのみの結果だけで良いなら
-w1
・変異の種類を細かく調節するなら-e（trueを排除）と-i（trueを残す）を使う。
-e'MAF[0]<0.05' select rare variants at 5% cutoff（変異率5%以下の部位は排除。-iなら55以下だけ残す）。
-i'QUAL>10' QUALが１０以上を残す。
など無数にある。他は公式ページを参照。
高度なフィルタリングでないなら、grepを使い絞り込むこともできる。
例えば、以下のようなMT（下側）にだけ検出される塩基置換を検出したいとする。
この部位のVCFファイルを確認する（青線部分）。
VCFは前もってマージしてある。
（↓実行したコマンド）
gatk -T CombineVariants -R ref.fasta -genotypeMergeOptions UNIQUIFY -V WT_SNPs.vcf -V MT_SNPs.vcf -o SNPsunion.vcf
GT:AD:DP:GQ:PL フィールドを拡大。
GT:AD:DP:GQ:PL ./. 1/1:0,12:12:36:450,36,0
WTでは変異が全くコールされていないため、./.になっている。MTでは塩基置換が100%おきているため0,12（REF=0、ALT=12）になっている。MTにだけおきている変異をGT:AD:DP:GQ:PL ./.の行を検索するだけで絞り込める。
ホモ(GT=1/1)
grep -n "GT:AD:DP:GQ:PL\s\./\.\s1/1" SNPs_invertselect.vcf > MT_homo.vcf
ヘテロ（GT=0/1）
grep -n "GT:AD:DP:GQ:PL\s\./\.\s0/1" SNPs_invertselect.vcf > MT_hetero.vcf
VCFtoolsのvcf-isecでも同等のことができます。VCFtoolsはベン図を書くために共通/非共通の数を知りたい場合などにも活用可能です。
追記
10Mb以下のスモールゲノムならbreseqを使うのも手です。
引用
How to compare two vcf files of two cultivars? — GATK-Forum
GATK4のインストール
#bioconda gatk4(link) gatk4-sparkもある
conda install -c bioconda gatk4
#初回はgatk-registerを走らせて実行ファイルのパスを叩く
gatk-register <path>/<to>/<your>/GenomeAnalysisTK.jar
