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文件格式解释见其他博文。
如何用 samtools 和 bcftools call snp,布布扣,bubuko.com
如何用 samtools 和 bcftools call snp
原文地址：http://www.cnblogs.com/freemao/p/3763949.html
