Bioinfo|bedtools-操作VCF文件

参考：https://mp.weixin.qq.com/s?__biz=MzIyNDI1MzgzOQ==&mid=2650393871&idx=1&sn=44314063cd7760e36a33a89efaa1c144&chksm=f01cafefc76b26f921ebe931b196f4cad47384687169b68686ec4b2c683b784a31d421dc58b9&scene=21#wechat_redirect

 初步设想在Bioinfo板块中分享一些常见的生信分析软件的使用，原则就是有现成的轮子就不去自己造了。

    当我们得到原始的测序数据（raw data）后，经过初步的QC质控得到clean data （公司或者自己的），可以通过FastQC软件对结果进行一系列的统计和图形展示R|fastqcr QC数据处理。当确定数据得到要求后，进行mapping分析得到bam文件后，就是重头戏检测SNP了。

    本文介绍的进行变异检测的工具为BCFTOOLS， http://samtools.github.io/bcftools/bcftools.html。

bcftools — utilities for variant calling and manipulating VCFs and BCFs 。

    呐，文本就通过几个常见的参数简单介绍一下bedtools是怎么“manipulating” VCF文件的。

一 mpileup 函数

使用方法：bcftools mpileup [OPTIONS] -f ref.fa in.bam [in2.bam […]]

示例：
bcftools mpileup -C 50 -m 2 -F 0.002  -d 1000 -q 20 -Q 20 -a DP,DP4,ADR -f hg19.fa A.bam B.bam

参数说明：

-f： 指定参考基因组，后跟需要变异检测的bam文件，bam较少时可以直接写；

bcftools mpileup  -f hg19.fa A.bam B.bam

-b：后跟一个bamlist的文件；

bcftools mpileup -f hg19.fa -b bamList

-C –adjust-MQ INT  矫正的MQ值，推荐是50；

-q, -min-MQ INT （MQ质量值）# 过滤MQ质量值；

-Q, –min-BQ INT （base质量值）# 过滤base质量值；

-r, –regions 只call特定的染色体，或者点，区域

形式为CHR|CHR:POS|CHR:FROM-TO|CHR:FROM-[,…]；

示例：
bcftools mpileup -r chr1 -f hg19.fa  A.bam | bcftools call -c -v > A.chr1.vcf
bcftools mpileup -r chr1:111111-122222 -f hg19.fa  A.bam B.bam | bcftools call -c -v > A.chr1_111111-122222.vcf

-R, –regions-file FILE  当有多个region时，将region存入文件，使用-R参数，文件格式：tab分割的三列，chr start end即可 ；

-a, –annotate LIST    在INFO/FORMAT列中增加一些如DP,DP4,AD等的信息，群call可能会需要；

Comma-separated list of FORMAT and INFO tags to output.

二 bcftools call函数

使用方法：bcftools call [OPTIONS] FILE

示例：bcftools mpileup -f hg19.fa A.bam B.bam | bcftools call -c -v -r chr1:111111-122222 -s A -V indels

示例说明：

只检测A样本，且只输出chr1:111111-122222区域的SNP信息，保存到test2.vcf文件中。

参数说明：

-c  -m为两种call的方式，旧版本使用-c ，现在默认使用-m，能解决大部分任务

-V, –skip-variants snps|indels 跳过SNP或者INDEL位点

-v, –variants-only  只输出变异位点

-s, –samples LIST 只检测此处给出的样本ID （通用参数）

-S, –samples-file FILE 只对此文件中列出的样本进行检测 （通用参数）

-r  -R为通用参数，与上面用法一致

三 bcftools filter 函数

使用方法：bcftools filter [OPTIONS] FILE

示例：bcftools filter  ALL.vcf TYPE="snp" -e 'DP < 20 || MQ < 20 || %QUAL<10 || DP < 20 ||MAF < 0.05 || ' 

示例说明：群体检测结果ALL.vcf文件中，只保留 质量值大于10，MQ大于20，MAF大于等于0.05，样本总深度大于20X的过滤后的SNP数据集。

参数说明：

TYPE=”snp”：只保留SNP信息；

-s：样本过滤，bcftools filter ALL.vcf -s A

-e 可加感兴趣的各种参数，常见的如下：

群体样本总深度DP：bcftools filter  ALL.vcf TYPE=”snp” -e ‘DP < 20’ 

前两个样本的DP：bcftools filter  ALL.vcf -e ‘FORMAT/DP[0-1] < 20 ‘ 

最小等位基因频率：bcftools filter  ALL.vcf -e ‘MAF < 0.05 ‘ |ee

质量值和mapping 质量值：bcftools filter  ALL.vcf -e ‘%QUAL<10 || MQ < 20 ‘ 

BAF值（alt depth / (ref depth + alt depth )）

其他过滤参数及使用，官网链接：

http://samtools.github.io/bcftools/bcftools.html#expressions

熟练使用以上及链接中的参数，就可以办到不写py或者pl脚本，强势过滤了。

四 bcftools index 函数

bgzip 压缩 vcf 文件为 gz 文件

bgzip -c A.vcf >A.vcf.gz ；bgzip -c B.vcf >B.vcf.gz

bcftools 为 gz 文件建索引

bcftools index -t A.vcf.gz ； bcftools index -t B.vcf.gz

五  bcftools merge 函数

使用说明：将多个VCF文件，合并成一个VCF文件

bcftools merge -m snps -f PASS,. --force-samples A.vcf.gz B.vcf.gz > A_B.merge.vcf

同样不用自己写脚本合并VCF文件，省事 且较少出错。

    建议多测试，熟悉以上参数的意义以及结果的不同，找到符合自己需求的参数组合，以得到理想的结果。

好了 ，拿着bedtools 去操纵VCF文件吧。