1. 程式人生 > >bcftools使用總結

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>
....

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

長流程

#一開始寫好引用,方便以後
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進行變異檢測的時候,是按照染色體排序順序進行的,並非多條染色體並行檢測的。因此,如果樣本資料量比較大的話,一般多個染色體並行。<