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
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
引數通過一個表示式來指定輸出格式,其中變數的寫法如下
-
%CHROM
代表VCF檔案中染色體那一列,其他的列,比如POS
,ID
,REF
,ALT
,QUAL
,FILTER
也是類似的寫法 -
[]
對於FORMAT欄位的資訊,必須要中括號括起來 -
%SAMPLE
代表樣本名稱 -
%GT
代表FORMAT欄位中genotype的資訊 -
\t
代表製表符分隔 -
\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>
....
----------------------------------------------------------------------------------------------
長流程
#一開始寫好引用,方便以後
ACC=AF086833
ebola=/vol2/agis/xiaoyutao_group/liuyunze/project/ebola
REF=$ebola/ref/$ACC.fa
SRR=SRR1553500
BAM=$ebola/align/$SRR.bam
R1=$ebola/raw/${SRR}_1.fastq
R2=$ebola/raw/${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
VARI=$ebola/variant
##bwa比對,samtools排序並構建索引,為了之後更快呼叫比對檔案
mkdir -p $ebola/align && cd $ebola/align
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM
mkdir -p $VARI
samtools faidx $REF
##第一種方法:bcftools召喚變異
samtools mpileup -uvf $REF $BAM | bcftools call -vm -Oz > bcftools.vcf.gz
##第二種方法:freebayes
freebayes -f $REF $BAM > $ebola/align/freebayes.vcf
##第三種方法:GATK(版本:4.0.7.0)
#注意:在使用GATK之前,需要先建立兩個參考基因組的索引檔案.dict和.fai【具體參見https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference】
#.dict中包含了基因組中contigs的名字,也就是一個字典;
#構建.dict檔案(原來要使用picard的CreateSequenceDictionary模組,但是現在gatk整合了此模組,可以直接使用)
gatk CreateSequenceDictionary -R $REF -O $ebola/ref/$ACC.dict
#.fai也就是fasta index file,索引檔案,可以快速找出參考基因組的鹼基
#構建
samtools faidx $REF
#gatk開始:
#必選 -I -O -R,代表輸入、輸出、參考
#接下來可以按照字母順序依次寫出來,這樣比較清晰
#-bamout:將一整套經過gatk程式重新組裝的單倍體基因型(haplotypes)輸出到檔案
#-stand-call-conf :低於這個數字的變異位點被忽略,可以設成標準30(預設是10)
gatk HaplotypeCaller -R $REF -I $BAM -O $ebola/align/HaplotypeCaller.vcf \
-bamout $ebola/align/$SRR.gatk.bam \
-stand-call-conf 30
# gatk用時3.95 minutes.
#<gatk補充>GATK進行變異檢測的時候,是按照染色體排序順序進行的,並非多條染色體並行檢測的。因此,如果樣本資料量比較大的話,一般多個染色體並行。
<