macでインフォマティクス

macでインフォマティクス

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

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

 

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

 

 

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

 

9/14追記

 

A、gatkを使う方法。

SelectVariants

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

 

以下の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つ以上のサンプルに対応可。

brewでgatkを入れている人はjava -jar GenomeAnalysisTK.jarの部分をgatkだけに変えて上記コマンドを打つ。

 

 

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以上を残す。

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

  

 

 

 

 

 

高度なフィルタリングでないなら、grepを使い絞り込むこともできる。

例えば、以下のようなMT(下側)にだけ検出される塩基置換を検出したいとする。

f:id:kazumaxneo:20170914140452j:plain

この部位のVCFファイルを確認する(青線部分)。

f:id:kazumaxneo:20170914140633j:plain

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はベン図を書くために共通/非共通の数を知りたい場合などにも活用可能です。


  

引用

 

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