bcftools使用總結

使用bcftools進行SNP calling

bcftools也能夠進行SNP calling。在以前的版本中，一般都是和samtools的mpileup命令結合使用, 命令以下

samtools mpileup -uf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf

因爲samtools和bcftools更新得都很快，只要有一個版本不對，採用上面的pipeline就會報錯。爲了減小版本不合適帶來的問題，bcftools的開發團隊將mpileup這個功能添加到了bcftools中。

在最新版的bcftools 中，只須要使用bcftools這一個工具就能夠實現SNP calling， 用法以下

bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa | bcftools call -mv -o raw.vcf

--fasta-ref參數指定參考序列的fasta文件，mpileup.bam是輸入文件，一般都是GATK 標準預處理流程獲得的bam文件。

須要注意的是mpileup命令雖然也會輸出VCF格式的文件，可是並不直接進行snp calling。下面的命令能夠生成VCF格式的文件

bcftools mpileup mpileup.1.bam --fasta-ref mpileup.ref.fa >mpileup.vcf

直接生成的VCF文件內容以下

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT HG00100
17 1 . A <*> 0 . DP=5; PL
17 2 . A <*> 0 . DP=5; PL
17 3 . G <*> 0 . DP=5; PL
17 4 . C <*> 0 . DP=5; PL
17 5 . T <*> 0 . DP=5; PL

裏面的每一條記錄並非一個SNP位點，而是染色體上每一個鹼基的比對狀況的彙總。這種信息官方稱之爲genotype likelihoods。

call命令纔是真正的執行SNP calling的程序，基本用法以下

bcftools call mpileup.vcf -c  -v -o variants.vcf

在進行SNP calling 時，必須選擇一種算法，有兩種calling算法可供選擇，分別對應-c和-m參數。-c參數對應consensus-caller算法， -m參數對應multiallelic-caller算法，後者更適合多種allel和罕見變異的calling。

-v參數也是經常使用參數，做用是隻輸出變異位點的信息，若是一個位點不是snp/indel, 不會輸出。

注：新版本bcftools中 call命令可替代view命令

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

其餘命令

1. annotate

annotate命令有兩個用途，第一個是用於註釋VCF文件，用法以下

bcftools annotate -a db.vcf -c ID,QUAL,+TAG view.vcf -o annotate.vcf

-a參數指定註釋用的數據庫文件，格式能夠是VCF, BED, 或者是\t分隔的自定義文件。在\t分隔的自定義文件中，必須包含CHROM, POS字段；-c參數指定將數據庫的哪些信息添加到輸出文件中。

第二個用途是編輯VCF文件，好比去除其中的某些註釋信息，或者去除某些樣本，用法以下

bcftools annotate -x ID,INFO/DP,FORMAT/DP  view.vcf -o remove.id.vcf

-x參數表示去除VCF文件中的註釋信息，能夠是其中的某一列，好比ID, 也能夠是某些字段，好比INFO/DP，多個字段的信息用逗號分隔；去除以後，這些信息所在的列並不會去除，而是用.填充。

2. concat

concat命令能夠將多個VCF文件合併爲一個VCF文件，要求輸入的VCF文件必須是排序以後的，若是包含多個樣本的信息，樣本的順序也必須一致。經典的應用場景包括合併不一樣染色體上的VCF文件，合併SNP和INDEL 兩種類型的VCF文件，用法以下

bcftools concat merge.2.a.vcf.gz merge.3.a.vcf.gz -o -o merge.vcf

還須要注意一點，輸入的VCF文件必須是通過bgzip壓縮的文件。

3. merge

merge命令也是用於合併VCF文件，主要用於將單個樣本的VCF文件合併成一個多個樣本的VCF文件。用法以下

bcftools merge merge.a.vcf.gz merge.b.vcf.gz -o merge.vcf

該命令要求輸入文件必須是通過bgzip壓縮的文件， 並且還須要有.tbi的索引。

4. isec

isec用於在多個VCF文件之間取交集，差集，並集等操做，經典的應用場景是對多種軟件的SNP calling 結果進行venn 分析。用法以下

bcftools isec A.vcf.gz B.vcf.gz -p dir

默認參數就是取交集，更多高級用法請參考幫助文檔。

5. stats

stats命令用於統計VCF文件的基本信息，好比突變位點的總數，不一樣類型突變位點的個數等。用法以下

bcftools stats view.vcf >  view.stats

輸出文件中記錄了不少類型的統計數據，重點介紹如下幾種

基本信息：

SN 0 number of samples: 3
SN 0 number of records: 15
SN 0 number of no-ALTs: 1
SN 0 number of SNPs: 11
SN 0 number of MNPs: 0
SN 0 number of indels: 3
SN 0 number of others: 0
SN 0 number of multiallelic sites: 1
SN 0 number of multiallelic SNP sites: 0

顛換和轉換的比例：

# TSTV, transitions/transversions:
# TSTV  [2]id  [3]ts  [4]tv  [5]ts/tv  [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV  0  8  3  2.67  8  3  2.67

Indel 長度分佈:

# IDD, InDel distribution:
# IDD [2]id [3]length (deletions negative) [4]count
IDD 0 -2 1
IDD 0 1 2
IDD 0 3 1

鹼基替換類型：

# ST, Substitution types:
# ST [2]id [3]type [4]count
ST 0 A>C 0
ST 0 A>G 0
ST 0 A>T 0
ST 0 C>A 1
ST 0 C>G 0
ST 0 C>T 4
ST 0 G>A 1
ST 0 G>C 1
ST 0 G>T 1
ST 0 T>A 0
ST 0 T>C 3
ST 0 T>G 0

輸出文件能夠用於plot-vcfstats命令，進行可視化, 這個腳本位於bcftools安裝目錄的misc目錄下。用法以下

plot-vcfstats view.stats -p output

-p參數指定輸出結果的目錄，這個腳本依賴latex 生成pdf 文件，因此係統上的latext 必定要安裝好。

輸出目錄下文件不少，詳細列表以下

├── counts_by_af.indels.dat
├── counts_by_af.snps.dat
├── depth.0.dat
├── depth.0.pdf
├── depth.0.png
├── indels.0.dat
├── indels.0.pdf
├── indels.0.png
├── plot.py
├── plot-vcfstats.log
├── substitutions.0.pdf
├── substitutions.0.png
├── summary.aux
├── summary.log
├── summary.pdf
├── summary.tex
├── tstv_by_af.0.dat
└── tstv_by_qual.0.dat

主要看summary.pdf文件就能夠了，該文件包含了不少信息

1.不一樣類型的突變位點彙總

2.插入缺失長度分佈圖

3.測序深度分佈

4.鹼基轉換類型分佈

6. index

index命令用於對VCF文件創建索引，要求輸入的VCF文件必須是使用bgzip壓縮以後的文件，支持.csi和.tbi兩種索引，默認狀況下創建的索引是.csi格式， 用法以下

bgzip view.vcf
bcftools index view.vcf.gz

運行成功後，會生成索引文件view.vcf.gz.csi。若是須要創建.tbi格式的索引，用法以下

bcftools index -t view.vcf.gz

tbi索引文件爲view.vcf.gz.tbi。

7.  view

view命令能夠用於VCF和BCF格式的轉換，用法以下

bcftools view view.vcf.gz -O u -o view.bcf

-O參數指定輸出文件的類型，b表明壓縮後的BCF文件，u表明未經壓縮的BCF文件，z表明壓縮後的VCF文件，v表明未經壓縮的VCF文件；-o參數指定輸出文件的名字。

還能夠根據樣本篩選VCF文件，用法以下

bcftools view view.vcf.gz -s NA00001,NA00002  -o subset.vcf

-s參數指定想要保留的樣本信息，多個樣本用逗號分隔。

若是樣本名稱添加了^前綴，表明去除這些樣本，好比-s ^NA00001,NA00002，這個寫法表示從VCF文件中去除NA00001,NA00002這兩個樣本的信息。

還能夠過濾突變位點，過濾的條件很是多，能夠根據突變位點的類型，基因型類型等等條件進行過濾，詳細的參數能夠參考軟件的幫助文檔，這裏只作一個基本示例

bcftools view view.vcf.gz -k -o known.vcf

-k參數表示篩選已知的突變位點，即ID那一列值不是.的突變位點。

8. query

query命令也是用於格式轉換，和view命令不一樣，query經過表達式來指定輸出格式，可定製化程度更高。用法以下

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' view.vcf.gz

-f參數經過一個表達式來指定輸出格式，其中變量的寫法以下

1. %CHROM
   表明VCF文件中染色體那一列，其餘的列，好比POS, ID, REF, ALT, QUAL, FILTER也是相似的寫法

2. []
   對於FORMAT字段的信息，必需要中括號括起來

3. %SAMPLE
   表明樣本名稱

4. %GT
   表明FORMAT字段中genotype的信息

5. \t
   表明製表符分隔

6. \n
   表明新的一行

    

輸出文件以下

11 2343543 A . NA00001=0/0 NA00002=0/0 NA00003=0/0
11 5464562 C T NA00001=./. NA00002=./. NA00003=./.
20 76962 T C NA00001=0/1 NA00002=1/1 NA00003=1/1

更多變量的寫法請參考官方文檔。

9. sort

sort 命令用於對VCF文件排序， 按照染色體位置進行排序，用法以下

bcftools sort view.vcf.gz -o sort.view.vcf

10. reheader

reheader命令有兩個用途，第一用途用於編輯VCF文件的頭部，第二個用途用於替換VCF文件中的樣本名。

替換樣本的用法以下

bcftools reheader -s sample.file view.vcf -o new.sample.vcf

-s參數指定須要替換的樣本名，內容以下

NA00001 NA1
NA00002 NA2
NA00003 NA3

第一列表明VCF文件中原始的樣本名稱，第二列表明替換後的樣本名稱，兩類之間用空格分隔，須要注意的是，樣本名不容許有空格。

編輯VCF文件的頭部的用法以下

bcftools reheader -h header.file  view.vcf -o new.header.vcf

-h參數指定新的header文件，內容以下

##fileformat=VCFv4.3
##reference=file:///seq/references/1000Genomes-NCBI37.fasta
##contig=<ID=11,length=135006516>
##contig=<ID=20,length=63025520>
....