【直播】我的基因组25：用bcftools来call variation

在这里，我们可以考虑比较对3种不同的bam文件来分别call variation的差异，探索对bam文件的不同过滤模式对snp calling的影响，分别是，原始比对的bam文件，去除低质量和多比对还有PCR重复的bam，以及用GATK进行realign的bam文件。

首先采取最经典的软件来做这个工作，后面我们会使用其它软件并且比较：

代码如下：

samtools mpileup -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta \

P_jmzeng.filter.rmdup.bam  |bcftools call -vmO z -o P_jmzeng.rmdup.bcftools.vcf.gz

在这里  \ 是换行符 在脚本里面很多时候一条命令很长  包含数据的绝对路径以及命令的选项  比如GATK 的很多命令就很长 加入换行符还是一行命令 在脚本中看起来更加清晰

所以只需要理解samtools mpileup 和bcftools call的参数的意义就可以了

首先我们来看samtools的mpileup这个命令我选择的参数

很明显，我就选择了3个参数，-f 来输入有索引文件的fasta参考序列； -g 输出到bcf格式给bcftools软件来做input，而-u就是说不要压缩，因为我们要通过管道直接把数据给bcftools的。这个mpileup以前为pileup，如果大家看到的教程比较陈旧，会有这样的疑问。

接下来我们看bcftools：

这次，我选择的参数有点多。

我-O选择的是z格式输出，没有添加多线程，所以对43G的bam文件处理了30个小时。

-v 就是说，我只输出有变异的位点，而不是全基因组坐标的输出。

-m比较难以理解，应该是一个位点，可以允许多种突变情况吧。

同时也可以选择多线程计算  因为全基因组的计算量很大,可以指定 --threads 24 便是使用24线程 。在先跑通最简单命令的基础之上，不要一成不变，可以通过读软件的文档来理解参数的意义，通过适当的调整参数，来更好地解决你的计算任务，但是参数不要随便改，因为不够robust的算法在面对不同复杂度的数据的时候，很有可能会带给你错误的结果

参考：

http://www.htslib.org/doc/samtools.html

https://samtools.github.io/bcftools/bcftools.html

http://cbsu.tc.cornell.edu/lab/doc/Variant_workshop_Part1.pdf

文：Jimmy、阿尔的太阳

图文编辑：吃瓜群众