macでインフォマティクス

macでインフォマティクス

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

VCFのユーティリティツール VCF-kit

 

 集団および量的遺伝学は、集団内の個体がどのように異なるかを調べる。これらの違いを特定することにより、様々な分析を行うことができる。例えば、遺伝子変異は、表現型の基礎を特定し、進化論的な質問に答えるために、または法医学を促進するために使用され得る。個々の群またはサンプル群にわたる遺伝的変異を表すための標準としてのVCFフォーマット(Danecek et al。、2011)の開発により、変異データを用いた遺伝子解析のための多数のツールの開発が促進された。その例には、遺伝子型の帰属(Browning and Browning、2007)、変異のアノテーション(Pedersen et al、2016)、変異型予測(Cingolani et al、2012)および集団遺伝解析(Pfeifer et al、2014)が含まれる。利用可能なツールはあるが、著者らの調べではVCFファイルを直接扱えることができる独自のツールセットが必要であることがわかった。その結果、著者らは、一連のツールをPythonで記述されたコマンドラインベースのプログラム、VCF-kitに組み込んだ。

 VCF-kitはVCFファイルを直接扱い、独自のさまざまな分析を実行する。VCF-kitには、バリアント検証のためのプライマー生成、樹形図作成、連鎖解析におけるシーケンスデータからの遺伝子型帰属、その他のツールを含むVCFファイルの処理と分析に不可欠なユーティリティが含まれている。

 

マニュアル

http://vcf-kit.readthedocs.io/en/latest/?badge=latest

以下のような機能がある。

f:id:kazumaxneo:20180326113244j:plain

公式より転載。 

 

インストール

依存

  • bwa (v 0.7.12)
  • samtools (v 1.3)
  • bcftools (v 1.3)
  • blast (v 2.2.31+)
  • muscle (v 3.8.31)

本体 Github 

GitHub - AndersenLab/VCF-kit: VCF-kit: Assorted utilities for the variant call format

本体はpipで導入でき、依存も全てbrewで導入できる。

pip install matplotlib
pip install yahmm
pip install numpy
pip install VCF-kit

#python以外の依存はbrewで導入する。
brew install bwa samtools bcftools blast muscle

> vk --help

$ vk --help

 

usage:

  vk <command> [<args>...] 

  vk setup

  vk -h | --help

  vk --version

 

commands:

  calc

  call

  filter

  geno

  genome

  hmm

  phylo

  primer

  rename

  tajima

  vcf2tsv

 

 

 

ラン

一般にmulti-VCFの入力を前提としている。1000genomesのphase3のVCFを使ってみる(リンク)。chr1つだと1GBくらいあるので、tabixで必要な領域だけダウンロードする。

tabix -p vcf  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz chr1:1-1000000

 

下の説明ではあまり統一できていないが、基本出力は非圧縮のvcfになる。データが多いなら、パイプでbcftoolsにつなぎ、gz出力する。

| bcftools view -O z > output.vcf.gz

 

 

calcリンク: VCFの遺伝子型/対立遺伝子の頻度または数を調べる。 

サンプル間で共有/固有のhomozygous genotypes。

vk calc sample_hom_gt input.vcf.gz

 genotypeの分布。 

vk calc genotypes input.vcf.gz > output
  • --frequency   calculate the frequencies of each genotype.

出力は行列(表)形式になる。

$ vk calc genotypes out.vcf.gz 

n ref het alt mis

1192 4 2 0 0

790 2 4 0 0

728 0 0 6 0

452 0 2 4 0

443 0 6 0 0

379 0 4 2 0

359 0 0 4 2

357 2 2 2 0

356 0 0 2 4

  1.  n - number of variants/records with the genotype distribution shown.
  2. ref - number of homozygous reference genotypes
  3. het - number of heterozygous genotypes
  4. alt - number of homozygous alternative genotypes (biallelic variants only)
  5. mis - number of missing genotypes

アレルの頻度

vk calc genotypes input.vcf.gz > output

$ vk calc spectrum out.vcf.gz 

n alt_allele_freq

1443 1.0

452 0.833333333333

119 0.75

548 0.666666666667

1080 0.5

973 0.333333333333

141 0.25

1192 0.166666666667

 

 

callリンク: サンガーシーケンス(AB1 / ABI、FASTQ、またはFASTAファイルを許容)から特定されたバリアントを、NGS解析から特定されたVCF内に存在するバリアントと比較する。 

vk call DL238_sanger.AB1 --ref=WBcel235 --all-sites illumina_seq.vcf.gz
  • --all-sites    Shows all calls from Sanger and VCF.
  • --vcf-targets   Only show variants present in the VCF.

出力は公式マニュアルで確認してください(リンク)。

 

 

filterリンク:homozygous reference, heterozygous, homozygous alternative、またはVCF由来の変異体の数または頻度に基づいて、バリアントをハードフィルタリングまたはソフトフィルタリングする。ここでの出力はgz圧縮にしている。

vk filter REF --min=1 input.vcf.gz | vk filter ALT --min=1 - | bcftools view -O z > one_homozygous.vcf.gz
  • --min=    Set a minimum threshold for the specified filter. An integer is interpretted as a minimum count. A float is interpretted as a minimum frequency.
  • --max=    Set a maximum threshold for the specified filter. An integer is interpretted as a maximum count. A float is interpretted as a maximum frequency.
  • --soft-filter=   When sepcified, the filter is added to the FILTER column and variants are NOT removed.
  • --mode=(+|x)    Specifies the mode for soft-filtering. When set to +, the filter is added to existing filters in the FILTER column. When set to x the soft-filter replaces existing filters in the FILTER column.

heterozygousの3つ以下のコールを出力。

vk filter HET --max=3 input.vcf.gz > output.vcf

 

 

genomeリンク: ゲノムをダウンロードし、さらにblastとbwa、samtoolsのindexをつける。ゲノムはgzで保存される。

NCBI genome database(リンク)からゲノムデータを探す。例えばいもち病菌を検索。

vk genome --search="Magnaporthe oryzae"

$ vk genome --search="Magnaporthe oryzae"

  /Users/user/.genome

  

  Searching...

  

    assembly_accession    bioproject    organism_name                     asm_name                 ftp_path

    --------------------  ------------  --------------------------------  -----------------------  ----------------------------------------------------------------------------------------------

    GCF_000002495.2       PRJNA1433     Magnaporthe oryzae 70-15          MG8                      ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/495/GCF_000002495.2_MG8

    GCF_000856605.1       PRJNA15041    Magnaporthe oryzae virus 1        ViralProj15041           ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/856/605/GCF_000856605.1_ViralProj15041

    GCF_000874825.1       PRJNA28297    Magnaporthe oryzae virus 2        ViralProj28297           ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/874/825/GCF_000874825.1_ViralProj28297

    GCF_000887575.1       PRJNA51685    Magnaporthe oryzae chrysovirus 1  ViralMultiSegProj51685   ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/887/575/GCF_000887575.1_ViralMultiSegProj51685

    GCF_000915895.1       PRJNA231682   Magnaporthe oryzae chrysovirus 3  ViralMultiSegProj231682  ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/915/895/GCF_000915895.1_ViralMultiSegProj231682

    GCF_000930815.1       PRJNA272442   Magnaporthe oryzae RNA virus      ViralProj272442          ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/930/815/GCF_000930815.1_ViralProj272442

    GCF_001028965.1       PRJNA286258   Magnaporthe oryzae virus 3        ViralProj286258          ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/028/965/GCF_001028965.1_ViralProj286258

  

  To download a genome and setup for use:

    

    vk genome ncbi --ref=<asm_name>

複数データベースが見つかった。上から1番目以外がchromosomeらしいので、これをダウンロードする。 ダウンロード時はasm_nameのカラムの名前を指定する。

vk genome ncbi --ref=MG8

データベースはncniのほかWormBase(wiki)にも対応している。ダウンロード時、chromosome名が複雑なら"Chromosome1"のようなシンプルな名前に変更される。 

 

ダウンロード済みのゲノムの一覧。

vk genome list

$ vk genome list

  /Users/user/.genome

    

    MG8

    ViralProj15041

ダウンロードパスの表示。

vk genome location

 

 

renameリンク:リネームする。

例えばGATKでコールしたVCFなので、先頭のサンプル名全てにGATK_をつける。

vk rename --prefix="GATK_" input.vcf.gz > output.vcf

 サンプルHG00097だけNEW_SAMPLE_NAMEにリネーム。

vk rename --subst=HG00097:NEW_SAMPLE_NAME input.vcf.gz > out.vcf

f:id:kazumaxneo:20180326121807j:plain

  • --prefix   先頭
  • --suffix 末尾
  • --subst 置換

3つは同時に使えるので、かなじ自由にリネームが可能になっている。

 

 

vcf2tsvリンク:表形式(行列)に変換。widelongは必須。

vk vcf2tsv wide input.vcf.gz > output.txt
  • wide   Prints a row for every CHROM:POS x # annotations (if --ANN is specified). A row will be output for every annotation.
  • long   Prints a row for every CHROM:POS x sample x annotation (if --ANN is specified) combination.
  • --print-header   Print a header row specifying column names.
  • --ANN   Parse annotation fields. Results in 1 row x (n)annotations x (n)samples when long or 1 row x (n)annotations x (n)variants when wide.

--ANNはSnpEffなどでアノテーション追記されたVCFなどに用いる。

 

phyloリンク:

 

 

検証できなかったが、他にもいくつかコマンドがある。詳細は各リンクから確認してください。

  • tajimaリンク:Tajima's D (wiki) を計算する。
  • genoリンク: FORMATフィールドフィルタ列(FT)を生成し、FILTER列に適用されたフィルタをFTフィールドに移す。
  • hmmリンク:連鎖解析のコマンド。カバレッジの低いシーケンスデータで使用するように設計されている(not as sophisticated as alternatives such as such as Impute2 or Beagle which can also be used to perform this type of analysis)。
  • primerリンク:見つかったgenoetypeについて、確認のプライマーを設計する。

 

 

引用

VCF-kit: assorted utilities for the variant call format

Cook DE, Andersen EC.

Bioinformatics. 2017 May 15;33(10):1581-1582.

 

参考ページ (Biostars)

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