bcftools使用笔记

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