【3.2】samtools 和 bcftools call snp

bcftools和samtools类似,用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件,后者为其二进制文件。官网: https://github.com/samtools/

更多资料:

Variant calling using SAMtools

1.bam文件都要经过处理

$  samtools mpileup -P ILLUMINA -f /bioinfo/data/iGenomes/Homo_sapiens/NCBI/build37.2/Sequence/BWAIndex/version0.6.0/genome.fa -EgD 199.sort.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文件中
还有一个参数-b 

如果bam文件较多,可以放到一个文件中,每行放一个,将含有bam名字的文件的文件名放到-b 后面就可以了

2.得到bcf文件以后,第二步执行命令:

$ bcftools call -c samtools_result.bcf -o samtools_result.vcf
 
Note: Neither --ploidy nor --ploidy-file given, assuming all sites are diploid

目前用的最新版1.3,跟之前的命令有些变化

报错

bcftools view -cNegv samtools_result.bcf > samtools_result.vcf

Error: Could not parse --min-ac Negv

备注:

bcftools的–max-idepth默认值是250,超过250的reads就跳过indel的计算了。把这个参数调高就可以看到插入缺失突变了。

参考资料

药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn