2019 7/29 condaインストール追記
2021 6/23 phyloコマンド追記
2023/09/02 追記
集団および量的遺伝学は、集団内の個体がどのように異なるかを調べる。これらの違いを特定することにより、様々な分析を行うことができる。例えば、遺伝子変異は、表現型の基礎を特定し、進化論的な質問に答えるために、または法医学を促進するために使用され得る。個々の群またはサンプル群にわたる遺伝的変異を表すための標準としての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
以下のような機能がある。
公式より転載。
インストール
VCF-kit has been tested with Python 3.6.
依存
- bwa (v 0.7.12)
- samtools (v 1.3)
- bcftools (v 1.3)
- blast (v 2.2.31+)
- muscle (v 3.8.31)
GitHub - AndersenLab/VCF-kit: VCF-kit: Assorted utilities for the variant call format
本体はpipで導入でき、依存も全てbrewで導入できる。ここではcondaでpython3.10の環境を作って導入する。
mmaba create -n vk python=3.10 -y
conda activate vk
pip install matplotlib
pip install yahmm
pip install numpy
pip install VCF-kit
#python以外の依存はbrewで導入する。
brew install bwa samtools bcftools blast muscle
#Bioconda(link) テストしてません
mamba install -c bioconda -y vcfkit
> vk --help
$ vk --help
usage:
vk <command> [<args>...]
vk setup
vk -h | --help
vk --version
commands:
calc
call
filter
genome
hmm
phylo
primer
rename
tajima
vcf2tsv
ラン
一般にmulti-VCFの入力を前提としている。1000genomesのphase3のVCFを使ってみる（リンク）。chr１つだと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
- n - number of variants/records with the genotype distribution shown.
- ref - number of homozygous reference genotypes
- het - number of heterozygous genotypes
- alt - number of homozygous alternative genotypes (biallelic variants only)
- 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の３つ以下のコールを出力。
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>
複数データベースが見つかった。上から１番目以外が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
- --prefix 先頭
- --suffix 末尾
- --subst 置換
３つは同時に使える。柔軟なリネームが可能になっている。
vcf2tsv（リンク）:表形式（行列）に変換。wideかlongは必須。
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（リンク）: VCFから系統解析を行う。
#fasta-alignment出力(vk phylo fastaコマンド)
vk phylo fasta input.vcf.gz
#multi sample VCF => tree (nj|upgma, ここではupgma)
vk phylo tree upgma --plot multi-sample.vcf.gz
*muscleを使う。PATHが通っていること。
検証できなかったが、他にもいくつかコマンドがある。詳細は各リンクから確認してください。
- 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)
