Samtools+bcftools Call SNP

  1. bwa mem算法 比对human_g1k_v37.fasta参考序列

  2. 用samtools sort或者picard.jar SortSam对比对后的sam进行默认排序,然后转化为bam文件

  3. 用picard MarkDuplicates对bam文件进行mark duplicates,以免PCR重复reads对后续call snp造成影响,没必要去除,可参照https://www.biostars.org/p/3917/

  4. 用samtools merge将多个lane合并(如果是多个lane的WGS的话),生成merged_marked.bam文件
    (以上步骤和原始数据可参看GATK best practices(Pre-processing+Variant discovery)

  5. 使用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 
      

    注:

    &&与& 区别:两者都表示“与”运算,但是&&运算符第一个表达式不成立的话,后面的表达式不运算,直接返回。而&对所有表达式都得判断。
    || 与| 区别:两者都表示“或”运算,但是||运算符第一个表达式成立的话,后面的表达式不运算,直接返回。而|对所有表达式都得判断。