samtools 之前博文已经介绍过一些常用的方法。本篇主要说下如何利用samtools 和 bcftools来call snp。

和其他工具一样，bam文件都要经过处理（另见博文）。假如对C17样本进行call snp, 数据为：

LC17-1_L002.sorted.rmp.rg.recal.bam

LC17-2_L006.sorted.rmp.rg.recal.bam

LC17-3_L002.sorted.rmp.rg.recal.bam

RC17-1_L003.sorted.rmp.rg.recal.bam

RC17-2_L004.sorted.rmp.rg.recal.bam

RC17-3_L004.sorted.rmp.rg.recal.bam

数据准备好后，执行命令：

samtools mpileup -P ILLUMINA -f RAP_cDAN.fasta -EgD \

LC17-1_L002.sorted.rmp.rg.recal.bam LC17-2_L006.sorted.rmp.rg.recal.bam \

LC17-3_L002.sorted.rmp.rg.recal.bam RC17-1_L003.sorted.rmp.rg.recal.bam \

RC17-2_L004.sorted.rmp.rg.recal.bam RC17-3_L004.sorted.rmp.rg.recal.bam \

> samtools_result.bcf

命令解释：

mpileup 是samtools中call snp的工具。

-P 指platform, 现在短reads测序一般是ILLUMINA。

-f 后跟参考序列，序列文件必须提前建好index。

-E， Extended BAQ（base alignment quality） computation, 如果有的话，会提高检测出MNPs的灵敏度，当然会轻微的减下特异度。

-g Compute genotype likelihoods and output them in the binary call format(BCF).

-D Output per-sample read depth

> 是将结果保存到samtools_result.bcf文件中

最终得到的samtools_result.bcf 是二进制文件，到此完成了call snp的第一步。

得到bcf文件以后，第二步执行命令：

bcftools view -cNegv samtools_result.bcf > samtools_result.vcf

命令解释：

veiw 是bcftools中主要的方法，‘Convert between BCF and VCF ， call variant candidates and estimate allele frequencies.’

-c  Call variants using Bayesian inference .

-N Skip sites where the REF field is not A/T/G/C. 一般的参考基因组序列都是由四种碱基组成，不知道还能有什么？ 难道是没出来部分的N 么 ？ 也可以不加这个参数，

   我测试过，这种情况非常非常少。

-e 其实也可以不加， 因为如果前面有-c 那么-e就默认被使用了。‘Perform max-likelihood inference only,including estimating the site allele frequency, testing Hardy-Weinberg equlibrium and testing associations with LRT’.

-g Call per-sample genotypes at variant sites.(用-c的方法)

-v output variant site only .我们当然只关心编译位点

> 将结果保存到samtools_result.vcf文件中。 是vcf格式的，具体关于vcf文件格式解释见其他博文。