众所周知，samtools加bcftools的找变异流程的运行速度是非常慢的，如果是全基因组，可能得耗费三五天。可以说是 比已经慢的发指的gatk流程还磨人！
但是，我们经常对某些细胞系进行有针对性的设计变异，比如BAF155的R1064K呀，H3F3A的K27呀，那我我们拿到高通量测序数据的时候，就肯定希望可以快速的看看这个基因是否被突变成功了。假设，我们已经得到了所有样本的sort好的bam文件，想看看自己设计的基因突变是否成功了，可以有针对性的只call 某个基因的突变！
samtools加bcftools流程
代码如下：
代码语言：javascript复制grep H3F3A ~/reference/gtf/gencode/protein_coding.hg19.position
# samtools, bcftools的参数需要自己看帮助文档哦。
samtools mpileup -r chr1:226249552-226259702 -ugf ~/reference/genome/hg19/hg19.fa *sorted.bam | bcftools call -vmO z -o H3F3A.vcf.gz
gunzip H3F3A.vcf.gz
~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old H3F3A.vcf >H3F3A.annovar
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg19 --geneanno --outfile H3F3A.anno H3F3A.annovar ~/biosoft/ANNOVAR/annovar/humandb/
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg19 --dbtype knownGene --geneanno --outfile H3F3A.anno H3F3A.annovar ~/biosoft/ANNOVAR/annovar/humandb/需要自己制作好基因的起始终止坐标文件，这样就可以找到自己的基因的位置，比如我的H3F3A是chr1:226249552-226259702，用bcftoolls简单的call variation即可，得到的vcf文件用annovar注释一下，看看是否在自己设计的蛋白质的某个位点的氨基酸！
