VCFやBCF を扱う bcftools

2020 11/17 追記

2021 4/244 追記、5/24 docker imageのリンク追加、6/2 callコマンド追記、9/17 論文引用、10/1 追記

2023/07/24 mpileup修正

bcftoolsは変異をコールしてバリアントコールフォーマットのVCFを出力したり、VCFやBCF（VCFのバイナリーフォーマット）を操作するツール。多様なコマンドから成る。samtoolsの論文で発表された（論文より　"The SAMtools package consists of two key components samtools and bcftools."）。

HP

manual

 インストール

Ubuntu18.04のdockerイメージでテストした（ホストOS；macos10.14）。

 本体　Github

#bioconda (link)
mamba install -c bioconda -y bcftools

#docker image(biocontainers)
#v1.9
docker pull biocontainers/bcftools:v1.9-1-deb_cv1
#v1.5
docker pull biocontainers/bcftools:v1.5_cv3

#docker image(mgibio/bcftools-cwl)
#v1.12
docker pull mgibio/bcftools-cwl:1.12

> bcftools

# bcftools 

Program: bcftools (Tools for variant calling and manipulating VCFs and BCFs)

Version: 1.8 (using htslib 1.8)

Usage:   bcftools [--version|--version-only] [--help] <command> <argument>

Commands:

 -- Indexing

    index        index VCF/BCF files

 -- VCF/BCF manipulation

    annotate     annotate and edit VCF/BCF files

    concat       concatenate VCF/BCF files from the same set of samples

    convert      convert VCF/BCF files to different formats and back

    isec         intersections of VCF/BCF files

    merge        merge VCF/BCF files files from non-overlapping sample sets

    norm         left-align and normalize indels

    plugin       user-defined plugins

    query        transform VCF/BCF into user-defined formats

    reheader     modify VCF/BCF header, change sample names

    sort         sort VCF/BCF file

    view         VCF/BCF conversion, view, subset and filter VCF/BCF files

 -- VCF/BCF analysis

    call         SNP/indel calling

    consensus    create consensus sequence by applying VCF variants

    cnv          HMM CNV calling

    csq          call variation consequences

    filter       filter VCF/BCF files using fixed thresholds

    gtcheck      check sample concordance, detect sample swaps and contamination

    mpileup      multi-way pileup producing genotype likelihoods

    roh          identify runs of autozygosity (HMM)

    stats        produce VCF/BCF stats

 Most commands accept VCF, bgzipped VCF, and BCF with the file type detected

 automatically even when streaming from a pipe. Indexed VCF and BCF will work

 in all situations. Un-indexed VCF and BCF and streams will work in most but

 not all situations.

> bcftools stats

$ bcftools stats

About:   Parses VCF or BCF and produces stats which can be plotted using plot-vcfstats.

         When two files are given, the program generates separate stats for intersection

         and the complements. By default only sites are compared, -s/-S must given to include

         also sample columns.

Usage:   bcftools stats [options] <A.vcf.gz> [<B.vcf.gz>]

Options:

        --af-bins <list>               allele frequency bins, a list (0.1,0.5,1) or a file (0.1\n0.5\n1)

        --af-tag <string>              allele frequency tag to use, by default estimated from AN,AC or GT

    -1, --1st-allele-only              include only 1st allele at multiallelic sites

    -c, --collapse <string>            treat as identical records with <snps|indels|both|all|some|none>, see man page for details [none]

    -d, --depth <int,int,int>          depth distribution: min,max,bin size [0,500,1]

    -e, --exclude <expr>               exclude sites for which the expression is true (see man page for details)

    -E, --exons <file.gz>              tab-delimited file with exons for indel frameshifts (chr,from,to; 1-based, inclusive, bgzip compressed)

    -f, --apply-filters <list>         require at least one of the listed FILTER strings (e.g. "PASS,.")

    -F, --fasta-ref <file>             faidx indexed reference sequence file to determine INDEL context

    -i, --include <expr>               select sites for which the expression is true (see man page for details)

    -I, --split-by-ID                  collect stats for sites with ID separately (known vs novel)

    -r, --regions <region>             restrict to comma-separated list of regions

    -R, --regions-file <file>          restrict to regions listed in a file

    -s, --samples <list>               list of samples for sample stats, "-" to include all samples

    -S, --samples-file <file>          file of samples to include

    -t, --targets <region>             similar to -r but streams rather than index-jumps

    -T, --targets-file <file>          similar to -R but streams rather than index-jumps

    -u, --user-tstv <TAG[:min:max:n]>  collect Ts/Tv stats for any tag using the given binning [0:1:100]

        --threads <int>                number of extra decompression threads [0]

    -v, --verbose                      produce verbose per-site and per-sample output

> bcftools convert

# bcftools convert

About:   Converts VCF/BCF to other formats and back. See man page for file

         formats details. When specifying output files explicitly instead

         of with <prefix>, one can use '-' for stdout and '.' to suppress.

Usage:   bcftools convert [OPTIONS] <input_file>

VCF input options:

   -e, --exclude <expr>        exclude sites for which the expression is true

   -i, --include <expr>        select sites for which the expression is true

   -r, --regions <region>      restrict to comma-separated list of regions

   -R, --regions-file <file>   restrict to regions listed in a file

   -s, --samples <list>        list of samples to include

   -S, --samples-file <file>   file of samples to include

   -t, --targets <region>      similar to -r but streams rather than index-jumps

   -T, --targets-file <file>   similar to -R but streams rather than index-jumps

VCF output options:

       --no-version               do not append version and command line to the header

   -o, --output <file>            output file name [stdout]

   -O, --output-type <b|u|z|v>    b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]

       --threads <int>            number of extra output compression threads [0]

GEN/SAMPLE conversion (input/output from IMPUTE2):

   -G, --gensample2vcf <...>   <prefix>|<gen-file>,<sample-file>

   -g, --gensample <...>       <prefix>|<gen-file>,<sample-file>

       --tag <string>          tag to take values for .gen file: GT,PL,GL,GP [GT]

       --chrom                 output chromosome in first column instead of CHROM:POS_REF_ALT

       --sex <file>            output sex column in the sample-file, input format is: Sample\t[MF]

       --vcf-ids               output VCF IDs in second column instead of CHROM:POS_REF_ALT

gVCF conversion:

       --gvcf2vcf              expand gVCF reference blocks

   -f, --fasta-ref <file>      reference sequence in fasta format

HAP/SAMPLE conversion (output from SHAPEIT):

       --hapsample2vcf <...>   <prefix>|<hap-file>,<sample-file>

       --hapsample <...>       <prefix>|<hap-file>,<sample-file>

       --haploid2diploid       convert haploid genotypes to diploid homozygotes

       --sex <file>            output sex column in the sample-file, input format is: Sample\t[MF]

       --vcf-ids               output VCF IDs instead of CHROM:POS_REF_ALT

HAP/LEGEND/SAMPLE conversion:

   -H, --haplegendsample2vcf <...>  <prefix>|<hap-file>,<legend-file>,<sample-file>

   -h, --haplegendsample <...>      <prefix>|<hap-file>,<legend-file>,<sample-file>

       --haploid2diploid            convert haploid genotypes to diploid homozygotes

       --sex <file>                 output sex column in the sample-file, input format is: Sample\t[MF]

       --vcf-ids                    output VCF IDs instead of CHROM:POS_REF_ALT

TSV conversion:

       --tsv2vcf <file>        

   -c, --columns <string>      columns of the input tsv file [ID,CHROM,POS,AA]

   -f, --fasta-ref <file>      reference sequence in fasta format

   -s, --samples <list>        list of sample names

   -S, --samples-file <file>   file of sample names

実行方法

stats

VCF/BCFの要約統計 

1つのVCFファイルだけ提供すると、非参照アレル頻度、デプス分布、品質およびサンプルごとのカウント、シングルトン統計などが出力される。2つのVCFファイルを提供した場合は、コンコーダンス（非参照対立遺伝子頻度別遺伝子型コンコーダンス、サンプル別遺伝子型コンコーダンス、非参照不一致）や相関関係などの統計量も出力される。サイトごとの不一致（PSD）も--verboseモードで表示される。

#one raw vcf file
bcftools stats input.vcf > stats

#two gzipped vcf file (only 2)
bcftools stats A.vcf.gz B.vcf.gz > stats

#視覚化（matplotlibが必要）
plot-vcfstats stats -p plot-vcfstats-test

#multiqcを使うことで要約統計(3sample 以上もOK）
bcftools stats sampleA.vcf > sampleA_bcftools_stats
bcftools stats sampleB.vcf > sampleB_bcftools_stats
bcftools stats sampleC.vcf > sampleC_bcftools_stats
multiqc .

 出力 

 multiqc (1サンプル)

annotate

VCF/BCFのカラムを編集するコマンド。カラムを追加したり消したりもできる。

bcftools annotate

merge　

重複していないサンプルセットの複数のVCF/BCFファイルをマージして、1つのマルチサンプルファイルを作成する。サンプル名がすべてのファイルでユニークでないとエラーになる（--force-samplesをつけると強制処理）。 vcftoolsのvcf-mergeコマンドに相当。

bcftools merge

concat　

VCF/BCFファイルを縦に結合する。。入力ファイルは、chとpositionでソートされてている必要がある。vcftoolsのvcf-concatコマンドに相当。

bcftools concat

isec   

VCFファイルの交叉、和（union）、余集合を作成する。

bcftools isec

plugin   user-defined plugins

様々なユーティリティのための共通フレームワーク。プラグインは、通常のコマンドと同じように使用できるが、名前の前に「+」が付いている。

bcftools plugin

 
reheader    

VCF/BCFファイルのヘッダを変更し、サンプル名を変更する。

bcftools reheader

sort

VCF/BCFをソートする。

bcftools sort

convert　

VCF/BCFを他のツールがサポートするフォーマットに変換したり、他のフォーマットからVCF/BCFに変換

#gVCF => VCF
bcftools convert input.gvcf -f hg38.fa --gvcf2vcf -o out.vcf 

# 23andme results => VCF
bcftools convert -c ID,CHROM,POS,AA -s SampleName -f 23andme-ref.fa --tsv2vcf 23andme.txt -Oz -o out.vcf.gz

• --gvcf2vcf　convert gVCF to VCF. Only sites with FILTER set to "PASS" or "." will be expanded.
• -G　convert IMPUTE2 output to VCF.
• -g　convert from VCF to gen/sample format used by IMPUTE2 and SHAPEIT. --hapsample2vcf　convert from hap/sample format to VCF.
• -H　convert from hap/legend/sample format used by IMPUTE2 to VCF
• -h　convert from VCF to hap/legend/sample format used by IMPUTE2 and SHAPEIT.
• --tsv2vcf　convert from TSV (tab-separated values) format (such as generated by 23andMe) to VCF.

call

SNP/indel コール

#samtoolsでpileupしてbcf出力（-uでBCF, -gでgVCF）, bcftools callのmultiallelic-caller でコール（consensus-callerよりも推奨）。（性染色体やオルガネラゲノムも含めて）すべてのサイトが2倍体であると仮定
samtools mpileup -uf ref.fa map.bam | bcftools call -mv > var.raw.vcf
#追記　bcftools pileupになった(2023/07/23)
bcftools mpileup -Ou -f ref.fa map.bam | bcftools call -mv -Ov > var.raw.vcf


#bcftools filterコマンドで低クオリティなコール（QUAL<20）とディープカバレッジな領域（DP>100、平均リードの深さの約2倍を推奨）からのコールをフィルタリング
bcftools filter -s LowQual -e '%QUAL<20 || DP>100' var.raw.vcf > var.flt.vcf

• -c, --consensus-caller   the original calling method (conflicts with -m)
• -m, --multiallelic-caller   alternative model for multiallelic and rare-variant calling (conflicts with -c)
• -s, --soft-filter <string>   annotate FILTER column with <string> or unique filter name ("Filter%d") made up by the program ("+")
• -e, --exclude <expr>   exclude sites for which the expression is true (see man page for details)

bcftoolsを使ったコンセンサスコールも参照（link）

view

サブセット、フィルタリング、変換する。

#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だと４スレッド使用)

#特定の個人だけ取り出す
bcftools view -O v -o output.vcf -s NA19062,NA19436 ALL.chrMT.phase3.genotypes.vcf.gz

以前の記事も参照

bcftools query

以前の記事を参照

bcftools consensus

以前の記事を参照

bcftools mpileup

以前の記事を参照

他にも様々なコマンドがあります。マニュアルを確認して下さい。

引用

A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data

Heng Li

Bioinformatics. 2011 Nov 1; 27(21): 2987–2993. Published online 2011 Sep 8. doi: 10.1093/bioinformatics/btr509 PMCID: PMC3198575

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3198575/

2021

Twelve years of SAMtools and BCFtools

Petr Danecek, James K Bonfield, Jennifer Liddle, John Marshall, Valeriu Ohan, Martin O Pollard, Andrew Whitwham, Thomas Keane, Shane A McCarthy, Robert M Davies, Heng Li

Gigascience. 2021 Feb 16;10(2)

関連