macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

サンプル間で共通する変異と固有の変異を抽出する

 

以前ショートリードからindelとSNVを検出するワークフローを紹介した。

 

 

複数サンプルがある場合、上記のような方法でVCFファイルを出力した後、サンプル間で共通のSNPs、サンプルごとに固有のSNPsなどを絞り込む必要が出てくるシチュエーションは多いと思われる。やり方は複数存在するが、今回は2つ紹介する。

 

A、gatkを使う方法。

以下の質問

How to compare two vcf files of two cultivars? — GATK-Forum

でGATKの開発チームがgatkを使ったやり方を指南している。以下の2ステップで固有/共通変異を識別する。

 

1、vcfファイルの結合(2サンプル)

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を増やすことで3つ以上のサンプルに対応可。

 

2、共通する変異の抽出

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -o intersect.vcf

 

3、固有の変異(= 全変異 ー 共通する変異)

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -invertSelect -o invertselect.vcf

-invertSelect: 共通しない変異

INFOの右端にサンプル名が付いており、それでどちらのサンプルか判別可能になっている。この情報を指標に簡単なスクリプトで分離することもできる。

 

1のステップで結合するVCFを追加することで3つ以上のデータの共通、固有変異の抽出も可能であるし、VCFフォーマットで記述されてるデータなら、他のindel toolで検出した変異をマージすることも簡単にできると思われる。

追記; gatk、freebayes、svabaはmegeできることを確認したが、platypusやFermikitなど一部のツールではVCFのフォーマットに独自の部分があり、mergeできなかった。

 

 

 

B、bcftoolsを使う方法

bcftoolsのisecを使う。やり方は公式ページのisecの説明部分に書かれている。

BCFtools manual page

 

bcftoolsを使うには、vcfファイルが圧縮され、indexがつけられている必要がある。

準備1 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ができる。

 準備2 インデックスファイル作成 

bcftools index sampleA.vcf.gz
bcftools index sampleB.vcf.gz
bcftools index sampleC.vcf.gz

 これで準備ができた。いかに4つの例を紹介する。

 

1、AとBの比較

bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output

 -p: 出力ディレクトリ(必須)

 

outputディレクトリにRWADMEと4つのファイルができる。仕分け条件はREADMEに書いてあるが、

・1つ目: sampleAに固有の変異(アレル)

・2つ目: sampleBに固有の変異(アレル)

・3つ目: sapleAから探したAとB共通の変異(アレル)

・4つ目: sapleBから探したAとB共通の変異(アレル)

3と4はどちらのVCFをベースに作るかの違いで、検出される部位の数は完全に同じである。VCFの由来かは4つのVCFファイルそれぞれのINFOフィールドの右端に記載されている。

 

2、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特異的な変異だけが出力される。

 

3、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 → 3つのVCFファイル

-n+3 → 3つ以上のVCFファイル (3を含む)

-n-3  → 3つ以下のVCFファイル(3を含む)

上記の場合、-n=3でも-n+3でも同じ結果になる。

 

4、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固有の変異が3つのファイルそれぞれに出力される。

 

それ以外にもたくさんオプションはある。

・出力する変異の種類を調節するには、-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が10以上を残す。

など無数にある。他は公式ページを参照。

 

 

 2つの方法を紹介した。他にもvcf-isecでも同等のことができる。

 

ベン図を書くために共通/非共通の数を知りたければ、下記URLのアメリエフさんのブログが非常に参考になります。VCFtoolsのコマンドで共通/非共通の変異を調べる例が説明されています。

アメリエフのブログ | 20121106