Samtools+bcftools Call SNP

bwa mem算法比对human_g1k_v37.fasta参考序列
用samtools sort或者picard.jar SortSam对比对后的sam进行默认排序，然后转化为bam文件
用picard MarkDuplicates对bam文件进行mark duplicates，以免PCR重复reads对后续call snp造成影响，没必要去除，可参照https://www.biostars.org/p/3917/
用samtools merge将多个lane合并（如果是多个lane的WGS的话），生成merged_marked.bam文件（以上步骤和原始数据可参看GATK best practices（Pre-processing+Variant discovery））
使用samtools mpileup + bcftools call SNP
- 先使用samtools mpileup将bam文件生成bcf文件（二进制文件）
  
  mpileup是samtools中call snp的工具，可以不使用-g参数，则会生成一个文本格式的文件，我们可以看到参考序列上每个碱基的比对结果：
  总共6列，分别是参考序列名（染色体），位置，参考序列的碱基，比对上的reads数，比对情况，比对上的碱基的质量
  
  生成bcf文件的主要参数：
  
  -g Compute genotype likelihoods and output them in the binary call format (BCF)(计算genotype likelihoods并以bcf格式输出)
  
  -f The faidx-indexed reference file in the FASTA format(用samtools faidx对参考序列建索引.fai文件，用其他软件建索引也可以的)
  
  -u Generate uncompressed VCF/BCF output(输出不压缩的bcf文件)
  
  -t Comma-separated list of FORMAT and INFO tags to output (case-insensitive)（按照自己的需求，输出vcf格式中常见的各种tags，比如-t DP输出vcf中的DP值；-t SP输出类似GATK call snp 出的VCF的GQ值？）
```
  samtools mpileup -go samtools.bcf -f ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -t DP -t SP merged_marked.bam
```
  后来发现，这里-t DP -t SP加不加都行，因为最后的bcftools call snp时，都会有DP和SP的tag
- 再用bcftools将bcf文件生成常用的vcf格式文件
  
  bcftools call snp的主要参数如下：
  
  -O Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v)(一般可以选择z,压缩的vcf格式)
  
  --threads Number of output compression threads to use in addition to main thread（主要用于压缩时的多线程，在-O b/z时使用）
  
  -m alternative modelfor multiallelic and rare-variant calling designed to overcome known limitations in -c calling model(一种更优的新算法，克服了旧算法的缺陷)
  
  -v output variant sites only(只输出突变位点信息)
  
  -f --format-fields 可以在输出的VCF文件中增加新的tag，例如增加GQ
  
  -V, --skip-variants 可以跳过输出snps|indels（比如我想只输出snps，则跳过indels）
  
  下面我就用一个常见的输出：
```
  bcftools call -vmO z -o bcftools_raw.vcf.gz samtools.bcf
```
  当然上述两步可以合起来运行（但是不能提交后台运行，反正会报错。。。）
```
  samtools mpileup -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -t DP -t SP merged_marked.bam |bcftools call -vmO z -o bcftools_raw.vcf.gz
```
- 用bcftools过滤不可靠的位点
  
  bcftools filter的主要参数如下：
  
  -e --exclude 主要用表达式方式去除过匹配上的位点，这个参数很关键，很多过滤都依靠这个表达式
  
  -g, --SnpGap 过滤INDEL附近的snp位点，比如-SnpGap 5则过滤INDEL附近5个碱基距离内的SNP
  
  -G, --IndelGap 过滤INDEL附近的INDEL位点。。。应该是这个意思
  
  -o, --output 输出文件的名称
  
  -O, --output-type 输出的格式，一般z和v都行
  
  -s, --soft-filter 将过滤掉的位点用字符串注释
  
  -S, --set-GTs set genotypes of failed samples to missing value (.) or reference allele (0)（将不符合要求的个体基因型改为./.）
  
  下面就简单的过滤QUAL小于10，DP值小于5，INDEL附近的位点
```
  bcftools filter -O v -o bcftools_filter.vcf -s LOWQUAL -e 'QUAL<10 || FMT/DP <5' --SnpGap 5 --set-GTs . bcftools_raw.vcf.gz
```
- 提取过滤后的SNP位点
```
  bcftools view -v snps bcftools_filter.vcf >bcftools_snp_filter.vcf
```
  或者在vcf文件中INFO列里，如果是INDEL的话，会标注出INDEL，因此提取SNP也可以：
```
  grep -v 'INDEL' bcftools_fiter.vcf >bcftools_snp_filter.vcf
```
注：
>&&与& 区别：两者都表示“与”运算，但是&&运算符第一个表达式不成立的话，后面的表达式不运算，直接返回。而&对所有表达式都得判断。
|| 与| 区别：两者都表示“或”运算，但是||运算符第一个表达式成立的话，后面的表达式不运算，直接返回。而|对所有表达式都得判断。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处