1 min read

bcftools使用笔记

什么是vcf?全称是variant call format,是用来记录SNP,INDEL和SV变异信息的文本。在GATK软件中得到最好的支持,当然SAMtools得到的结果也是VCF格式,和GATK的VCF格式有点差别。可参考以下链接文件,深入学习。

1.查看vcf文件的头部信息

bcftools view -h shrimp.vcf.gz

2.提取指定个体产生新vcf文件

如果只想分析vcf文件中部分样本,可以利用view命令中的—S选项过滤个体。

bcftools view -S selectedinds.txt shrimp.vcf.gz -Ov > shrimp104.vcf

-S选项后边跟一个文本文件,每一行为保留个体的ID编号。如果样品少,也可以 在-S 后边直接跟样品的ID号。>为重定向符号,表示把保留的个体信息存到shrimp.vcf文件中;-Ov表示输出未经压缩的vcf文件。

3.压缩vcf为gz压缩文件

bcftools一般要求读入文件为gz格式。而且gz格式的文件,压缩幅度非常大,便于传输。因此可以利用view命令把vcf文件压缩为gz格式。

bcftools view shrimp.vcf -Oz -o shrimp.vcf.gz --threads 8

-Oz表示输出格式为压缩文件gz格式;-o后边跟压缩文件名字;–threads 8 表示通过8个线程并行压缩shrimp.vcf文件。

4.提取等位基因和基因型信息

参考这个链接Extracting allele, Genotype from VCF file.

用到的命令为“query”,说明书中该命令的功能为“transform VCF/BCF into user-defined formats”,即把VCF/BCF文件转换为用户定义格式。

bcftools query -f '%CHROM %ID %POS %REF %ALT [ %TGT]\n' query.vcf.gz -o query.extract.txt

上述命令提取vcf文件中染色体、基因型等信息,输出为空格分隔的文本文件。

其中:

  • %CHROM 染色体列
  • %ID 变异位点名称
  • %POS 变异位点位置
  • %REF 参考等位基因
  • %ALT 变异等位基因
  • %TGT 字符格式如A/G的基因型;%GT为0/1格式的基因型

也可以输出为英文逗号分隔的csv文件。若想用tab分隔,符号是bcftools query -f ‘%CHROM,%ID,%POS,%REF,%ALT[,%TGT]’ query.vcf.gz -o query.extract.csv`

5.变异位点的基本统计分析

stats命令用于统计VCF文件的基本信息,比如突变位点的总数,不同类型突变位点的个数等。用法如下:

bcftools stats view.vcf > view.stats

然后会生成一个名为view.stats的文本文件。

然后调用一下命令,进行可视化输出:

plot-vcfstats view.stats -p output

多个结果文件保存在output文件夹下。其中summary.pdf文件,包括了汇总分析的结果。主要包括:

  • 转换(Transition:嘌呤(purine)转换A(腺嘌呤)-G(鸟嘌呤);嘧啶(pyrimidine)转换C(胞嘧啶)-T(胸腺嘧啶))和颠换(Transversion:A-C;A-T;G-C;G-T)类型分布统计。尽管颠换类型要多于转换,但在基因组上,转换的实际发生频率要高于颠换类型。 一般转换与颠换的发生比例为2:1。

需要注意,plot-vcfstats脚本在bcftools安装目录misc文件夹下,这是一个perl脚本,会调用python3绘图模块Matplotlib。如果没有安装该模块,可以通过pip3命令pip3 install -U matplotlib进行安装。如果没有pip3,通过sudo apt-get install python3-pip命令安装。

plot-vcfstats生成pdf文件还需要pdf-latex,如果系统没有安装latex,通过sudo apt-get install texlive-full进行安装。

上述命令在ubuntu 18.04系统中执行。

6. 替换染色体名称

需要用到annotate命令中的–rename-chrs参数。

命令形式:

bcftools annotate --rename-chrs NewChrName.txt old.xxx.vcf.gz -Oz -o new.xxx.gz --threads n

其中NewChrName.txt文件中包括了旧和新染色体名称的对应关系。--threads可以设置多线程加快新vcf文件的生成速度。

可以通过bioawk命令来提取旧染色体名字。

bioawk -c vcf ' { print $chrom } ' old.vcf.gz | sort -n | uniq > ChrName.txt