macでインフォマティクス

macでインフォマティクス

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

VCFtoolsでVCFを管理、編集する

VCFtoolsは、バリアントコールフォーマットのVCFファイルのマージ、ソートやフィルタリング、固有変異の抽出などができるツール。 よく使いそうなコマンドに限って紹介する。

  

マニュアル

https://vcftools.github.io/perl_module.html

インストール

git clone https://github.com/vcftools/vcftools.git
cd vcftools
./autogen.sh
./configure make
sudo make install

 

ラン

マニュアル

man vcftools

入力はgz圧縮が望ましいとされる。

bgzip my_file.vcf 
tabix -p vcf my_file.vcf.gz

 

 

vcf-merge  

2つ以上のVCFをマージする。

vcf-merge A.vcf.gz B.vcf.gz C.vcf.gz | bgzip -c > out.vcf.gz

 ポジションが同じ部位は1行にまとめられる。右端の4列を見ればどうなっているかわかる。

AとBに存在。Cにはない。

GT:GQ:DP:PL:AD 1/1:75:25:1115,75,0:0,25 1/1:60:37:795,60,0:0,37 .

Aだけに存在。

GT:GQ:DP:PL:AD 1/1:99:59:2415,181,0:0,59 . .

 

vcf-concat 

2つ以上のVCFをコメント(##)を除き結合する。(ポジションが同じ部位をマージしたりソートしたりはしない)。

vcf-concat A.vcf.gz B.vcf.gz C.vcf.gz | gzip -c > out.vcf.gz 

 

vcf-sort

chrとポジションでソートする。

vcf-sort file.vcf.gz

ソートコマンドの sort -k1,1d -k2,2n(1列目ポジションソート+2列目ポジションソート)を実行したのと同じ結果になる。

 

vcf-consensus

VCFの変異部位を元にレファレンスfastaを改変する。

cat ref.fa | vcf-consensus file.vcf.gz > out.fa 
  • --sample <name> If not given, all variants are applied
  • --haplotype <int> Apply only variants for the given haplotype (1,2)
  • --iupac-codes Apply variants in the form of IUPAC ambiguity codes

 

vcf-convert

VCFのフォーマットを3.3から4.0に変える。

zcat file.vcf.gz | vcf-convert -r reference.fa > out.vcf
  • --version <string> 4.0, 4.1, 4.2

 

vcf-isec

共通の変異や固有の変異を検出する。

2つ以上のVCFに共通して存在している変異を抽出。

vcf-isec -n +2 A.vcf.gz B.vcf.gz | bgzip -c > out.vcf.gz 
  • -n [+-=]<int> Output positions present in this many (=), this many or more (+), or this many or fewer (-) files.

Aに存在して、B、Cには存在しない変異を抽出。

vcf-isec -c A.vcf.gz B.vcf.gz C.vcf.gz | bgzip -c > out.vcf.gz
  • -c Output positions present in the first file but missing from the other files.  
  • -f  Continue even if the script complains about differing columns, VCF versions, etc.
  • -p <path> If present, multiple files will be created with all possible isec combinations. (Suitable for Venn Diagram analysis.)
  • -r <list|file> Do only the given regions (comma-separated list or one region per line in a file).
  • -w <int> In repetitive sequences, the same indel can be called at different positions. Consider records this far apart as matching (be it a SNP or an indel).

 

vcf-annotate

MQが30以上の部位を"PASS"させる。(MQリンク

vcf-annotate -f MinMQ=30 input.vcf > output.vcf
  • -f <file|list> Apply filters, list is in the format flt1=value/flt2/flt3=value/etc. If argument to -f is a file, user-defined filters be applied. See User Defined Filters below.
  • -H Remove lines with FILTER anything else than PASS or "."
  • -r <list> Comma-separated list of tags to be removed (e.g. ID,INFO/DP,FORMAT/DP,FILTER).

出力を確認すると、MQが30以上のものがフィルターをPASSしている。

f:id:kazumaxneo:20170916144847j:plain

 

デプス、アレルのデプス、マッピングクオリティなどのフィルターをかける。

vcf-annotate -f +/-a/c=3,10/d=10/-D/MinMQ=40/ input.vcf > output.vcf 

以下のようなフィルターがある。

f:id:kazumaxneo:20170916151147j:plain

上ではMinMQと書いているが、qと書いてもOK。PASSしなかった部位はその名前がつく。MinDPが通らなかったらFilterのカラムにMinDPと付く。

 

 

vcf-compare

 VCFを比較して固有の変異共通の変異の数を出力する(ベン図を書くときに使える)。

 vcf-compare A.vcf.gz B.vcf.gz C.vcf.gz > output
  • -a Ignore lines where FILTER column is anything else than PASS or '.'
  • -c <list|file> Same as -r, left for backward compatibility. Please do not use as it will be dropped in the future.
  • -d Debugging information. Giving the option multiple times increases verbosity
  • --cmp-genotypes Compare genotypes, not only positions
  • -ignore-indels Exclude sites containing indels from genotype comparison
  • -R <file> Compare the actual sequence, not just positions. Use with -w to compare indels.

VCF2ファイルを比較すると以下のような出力が得られる。

o$ cat out

# This file was generated by vcf-compare.

# The command line was: vcf-compare(v0.1.14-12-gcdb80b8) HaplotypeCaller.vcf.gz UnifiedGenotyper.vcf.gz

#

#VN 'Venn-Diagram Numbers'. Use `grep ^VN | cut -f 2-` to extract this part.

#VN The columns are: 

#VN        1  .. number of sites unique to this particular combination of files

#VN        2- .. combination of files and space-separated number, a fraction of sites in the file

VN 5 UnifiedGenotyper.vcf.gz (7.9%) #VCF1固有

VN 19 HaplotypeCaller.vcf.gz (24.7%) #VCF2固有

VN 58 HaplotypeCaller.vcf.gz (75.3%) UnifiedGenotyper.vcf.gz (92.1%) #共通

#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.

SN Number of REF matches: 58

SN Number of ALT matches: 57

SN Number of REF mismatches: 0

SN Number of ALT mismatches: 1

SN Number of samples in GT comparison: 0

 

 

vcf-to-tab

VCFをシンプルなタブ形式に変換する。

zcat file.vcf.gz | vcf-to-tab > out.tab

以下のような出力が得られる。

user$ head -6 out.tab

#CHROM  POS     REF     Z

chr     831647  C       T/T

chr     943495  G       A/A

chr     1012958 G       T/T

chr     1114645 T       TTTAATTAGGGTGATAAATTCG/TTTAATTAGGGTGATAAATTCG

chr     1204616 G       A/A

 

 

vcf-contrast  A tool for finding differences between groups of samples, useful in trio analysises, cancer genomes etc.

fill-an-ac Fill or recalculate AN and AC INFO fields.

fill-fs Annotates the VCF file with flanking sequence (INFO/FS tag) masking known variants with N's. Useful for designing primers.

fill-ref-md5 Fill missing reference info and sequence MD5s into VCF header.

vcf-annotate The script adds or removes filters and custom annotations to VCF files. To add custom annotations to VCF files, create TAB delimited file with annotations such as

vcf-fix-newlines Fixes diploid vs haploid genotypes on sex chromosomes, including the pseudoautosomal regions.

vcf-fix-ploidy Fixes diploid vs haploid genotypes on sex chromosomes, including the pseudoautosomal regions.

vcf-indel-stats Calculate in-frame ratio.

vcf-phased-join Concatenates multiple overlapping VCFs preserving phasing.

vcf-query Powerful tool for converting VCF files into format defined by the user. Supports retrieval of subsets of positions, columns and fields.

vcf-shuffle-cols Reorder columns

vcf-stats Outputs some basic statistics: the number of SNPs, indels, etc.

vcf-subset Remove some columns from the VCF file.

 

追記

変異が数十万以上あって動作が緩慢なら、最初にcyvcf2で不要な情報を除外してから使ってください。cyvcf2は2017年に発表された高速なVCFのパーサです。

 

 

引用

The variant call format and VCFtools

Petr Danecek,1,† Adam Auton,2,† Goncalo Abecasis,3 Cornelis A. Albers,1 Eric Banks,4 Mark A. DePristo,4 Robert E. Handsaker,4 Gerton Lunter,2 Gabor T. Marth,5 Stephen T. Sherry,6 Gilean McVean,2,7 Richard Durbin,1,* and 1000 Genomes Project Analysis Groupcorresponding author

Bioinformatics. 2011 Aug 1; 27(15): 2156–2158.

 

BIostars 質問 VCFのマージ

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