macでインフォマティクス

macでインフォマティクス

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

バリアントコール結果のVCFフォーマット

2018 10/25 追記

2019 8/9 コマンド追記

2020 10/14 vcflibのインストールコマンドを最後に追記

 

次世代データからリファンレンスゲノムの変異検出を行うと、Variant Call Format(VCF)という形式で出力されることが多い。VCFの詳細はsamtoolsのVCFフォーマットオフィシャルページに書いてあるが、そのフォーマットについてもう少し噛み砕いて説明を残しておく。

 

The Variant Call Format (VCF) Version 4.2 Specification

f:id:kazumaxneo:20180323191132j:plain

 

What is a VCF and how should I interpret it?

https://gatkforums.broadinstitute.org/gatk/discussion/1268/what-is-a-vcf-and-how-should-i-interpret-it

 

以下はGATK haplotypecallerで変異検出して出力されたVCFファイルのコメント1行と変異コールの1行を表示したものである。 

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample_1 
chr 7438 . T C 1086.77 PASS AC=2;AF=1.00;AN=2;DP=25;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=43.50;QD=34.24;SOR=1.136 GT:AD:DP:GQ:PL 1/1:0,25:25:75:1115,75,0

 

カラムを左から見ていく。まず1−5列目

#CHROM POS ID REF ALT 
chr 7438 . T C

#CHROM リファレンスfastaのヘッダー名

#POS ポジション

#ID ID

#REF リファレンス塩基

#ALT リードから読み取ったリファレンスと異なる塩基(上だとC -> Tに置換)

 

6列目

QUAL
1086.77

quality scoreの値。この値は検出ソフトによって定義が変わるので、ツール間で比較することはできない。gatkのhaplotypercallerの場合、マッピングのクオリティを表すmapping quality scoreをリードの数だけ足した合計値に近い値になるが(すなわち大きいほどエラーの可能性は低い)、もう少し複雑な計算をされているしい。いずれにせよカバレッジによって値は大きく変わる。例えばカバレッジが1ならどんなにクオリティが高くても2桁にしかならない。繰り返し配列を含む領域では5、6桁以上になったりするので、この値は直接フィルタリングに使えるのもではないらしい。7列目にはこの値をリード数で正規化したQD値が記載されている。

mapping qualityについてはこちらを参照。

 

7列目

FILTER
PASS

多くのidnel解析ソフトではフィルリングが終わっているとPASSと付く。フィルタリングを通らないと、例えばGATKでは自分が決めた名前 (e.g., basic_snp_filter) と付く。またフィルタリングが実行されていないと "."(ピリオド)になっている。

VCFフォーマットの定義では1-7列目まで必須となっている。

 

8列目

INFO
AC=2;AF=1.00;AN=2;DP=25;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=43.50;QD=34.24;SOR=1.136

追加情報。 "英数字混じりの文字=数字"で状態を表現し、それをセミコロンで区切るように規定されている。標準で以下のような内容がある。snpEffなどで検出した変異のアノテーション情報などもここに記載される。

 

  • AA ancestral allele 先祖型の対立遺伝子
  • AC allele count in genotypes, for each ALT allele, in the same order as listed
  • AF allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes
  • AN total number of alleles in called genotypes
  • BQ RMS base quality at this position
  • CIGAR cigar string describing how to align an alternate allele to the reference allele
  • DB dbSNP membership
  • DP combined depth across samples, e.g. DP=154
  • END end position of the variant described in this record (esp. for CNVs)
  • H2 membership in hapmap2
  • MQ RMS mapping quality, e.g. MQ=52
  • MQ0 Number of MAPQ == 0 reads covering this record
  • NS Number of samples with data
  • SB strand bias at this position
  • SOMATIC indicates that the record is a somatic mutation, for cancer genomics
  • VALIDATED validated by follow-up experiment

 GATKの場合、上のようにQD(Qual By Depth.)、FS(Fisher Strand)、SORなどのパラメータも追加されている。

 

ここまでが必須のフィールドであるが、gatkでは9、10列目にサンプルのgenotype、カバレッジなどの情報を記している(詳細はGATKのhttp://gatkforums.broadinstitute.org/gatk/discussion/1268/what-is-a-vcf-and-how-should-i-interpret-itを参照)。以下のように9-10列目は連動している。

FORMAT       sample1
GT:AD:DP:GQ:PL 1/1:0,25:25:75:1115,75,0

セミコロンで5つの情報を表している。順番に説明する。

GT: genotype。0/0ならリファレンスのホモ、1/1ならコールされた塩基のホモ、0/1ならヘテロと解釈されたことになる。上だと1/1なのでALTのホモになる。ポリプロイドゲノムだと他の数値も出るらしい。

AD: allele depth。フィルター無し状態でのカバレッジ。バリアントコーラーをパスしていないリードも含む。ただし冗長にマッピングされたリードなどは含まれない。REF、ALTの順番に並んでいる。上の状態なら0,25なのでREF=0、ALT=25(リンク)。

DP: filtered depth。フィルターあり状態でのカバレッジ。バリアントコーラーをPASSしたリードのカバレッジ。ただし、ADと違って冗長にマッピングされた意味の薄いリードのカバレッジも含んでいる(リンク)。

GQ: Quality of the assigned genotype。Phread-scaleでのジェのタイピングの信頼値で次のPLと関係している。上限は99に規定されており、高いほどホモ、ヘテロコールの信頼度が高いらしい。

PL: Normalized" Phred-scaled likelihoods of the possible genotypes.正規化したもっともらしさの尺度らしい。GTは大抵0/0から1/1の間を取るので、PLは大抵0になるよう補正されている。

 

注: vcfファイルのheaderにも説明は記載されている。

追記

 C++のVCF管理ライブラリvcflib(Github)のコマンドvcffilterを使い、例えば

 DP>5、MQ>30、QD>20でフィルタリングすると、freebayesやsamtools、varscanなどでコールしたレアバリアントの大半は除去されます(*1)。

vcffilter -f "DP > 5 & MQ > 30 & QD > 20" file.vcf >filtered.vcf

#さらにhomozygoteのコールだけ抽出。
cat filtered.vcf | vcffilter -g "GT = 1|1" | vcffixup - | vcffilter -f "AC > 0" >results.vcf

#DP combined depth across samples, e.g. DP=154
#MQ RMS mapping quality, e.g. MQ=52
#QD Variant call confidence normalized by depth of sample reads supporting a variant
  •  -g, --genotype-filter    specifies a filter to apply to the genotype fields of records

Filters are specified in the form "<ID> <operator> <value>:

  • -f "DP > 10" # for info fields
  • -g "GT = 1|1" # for genotype fields
  • -f "CpG" # for 'flag' fields

Operators can be any of: =, !, <, >, |, &

 

format変換

bcftools view

先にvcf <=> bcf変換などができるbcftools viewをまとめておく(他にも機能あり)。

#bcftoolsのinstall
#bioconda (link)
conda install -c bioconda -y bcftools
#vcf => bcf
bcftools view -Ou input.vcf -o output.bcf

#bcf => vcf
bcftools view -Ov input.bcf -o output.vcf

#vcf => vcf.gz
bcftools view -Oz input.vcf -o output.vcf.gz

#bcf => vcf.gz
bcftools view -Oz input.bcf -o output.vcf.gz

#vcf => compressed bcf
bcftools view -Ob input.vcf -o output.bcf.gz

#bcf => compressed bcf
bcftools view -Ob input.bcf -o output.bcf.gz

#indexing(tabix)installするにはhtslibをビルドするかcondaを使う(*1)
tabix -p vcf output.vcf.gz
=> output.vcf.gz.tbiができる

#BGZF compressin and create index (bgzip)installするにはhtslibをビルドするかcondaを使う(*2)
bgzip -i output.bcf
=> output.bcf.gz.gbiができる(-@ 4だと4スレッド使用)
  • -O <b|u|z|v>       b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]

 

 

追記

Variant annotations in VCF format

http://snpeff.sourceforge.net/VCFannotationformat_v1.0.pdf

 

VCFファイルはbedtoolsでマージしたり編集することができます。bedtoolsについては下記にまとめています。

 サンプル間の共通する変異を抽出するなら以下を参照してください。

変異が起きた遺伝子のアノテーションをつけるなら以下を参照してください。

変異の種類をカウントしたりBEDに変換するなら以下を参照してください。

VCFのフィルタリング

高速なパーサ

vcfにgenotype情報を追記する(そのあとfilteringする)。

要約統計

 視覚化


引用

 VCFフォーマット4.2の仕様

https://samtools.github.io/hts-specs/VCFv4.2.pdf

 

 VCFフォーマット wikiまとめ

 

Biostars VCFのindex

 

追記 2018/04

1000genomesによると、VCFなどのフォーマットの策定と今後の拡大は、Global Alliance for Genomics and Health Dataが担っているらしい(100genomeリンク)。

Global AllianceのHPに、SAM/BAMやCRAM、VCFのDOcumentがある。

https://www.ga4gh.org/ga4ghtoolkit/genomicdatatoolkit/

VCF4.3も決まっているらしい。

http://samtools.github.io/hts-specs/VCFv4.3.pdf

=> 2020年現在すでに使用されている。

 

こちらもどうぞ

What is a VCF?

https://wabi-wiki.scilifelab.se/pages/viewpage.action?pageId=19988483

 

Variant Finding

http://bioinformatics.ucdavis.edu/docs/2015-september-workshop/_downloads/Tu-Variant-Finding-MB.pdf

http://bioinformatics.ucdavis.edu/docs/2015-september-workshop/_downloads/Tu-Variant-Finding-MB.pdf

 

vcflibのインストール

git clone --recursive https://github.com/vcflib/vcflib.git
cd vcflib
make -j
cd bin/

./vcffilter

# ./vcffilter 

usage: ./vcffilter [options] <vcf file>

 

options:

    -f, --info-filter     specifies a filter to apply to the info fields of records,

                          removes alleles which do not pass the filter

    -g, --genotype-filter specifies a filter to apply to the genotype fields of records

    -k, --keep-info       used in conjunction with '-g', keeps variant info, but removes genotype

    -s, --filter-sites    filter entire records, not just alleles

    -t, --tag-pass        tag vcf records as positively filtered with this tag, print all records

    -F, --tag-fail        tag vcf records as negatively filtered with this tag, print all records

    -A, --append-filter   append the existing filter tag, don't just replace it

    -a, --allele-tag      apply -t on a per-allele basis.  adds or sets the corresponding INFO field tag

    -v, --invert          inverts the filter, e.g. grep -v

    -o, --or              use logical OR instead of AND to combine filters

    -r, --region          specify a region on which to target the filtering, requires a BGZF

                          compressed file which has been indexed with tabix.  any number of

                          regions may be specified.

 

Filter the specified vcf file using the set of filters.

Filters are specified in the form "<ID> <operator> <value>:

 -f "DP > 10"  # for info fields

 -g "GT = 1|1" # for genotype fields

 -f "CpG"  # for 'flag' fields

 

Operators can be any of: =, !, <, >, |, &

 

Any number of filters may be specified.  They are combined via logical AND

unless --or is specified on the command line.  Obtain logical negation through

the use of parentheses, e.g. "! ( DP = 10 )"

 

For convenience, you can specify "QUAL" to refer to the quality of the site, even

though it does not appear in the INFO fields.

 

 

参考

*1

https://www.biostars.org/p/51439/