对感兴趣的基因call variation

这是生信技能树创始人jimmy师兄关于针对特定基因做call variation的笔记，坦白讲，我看了好几遍才看懂个大概，菜鸟真煎熬，分享出来，希望帮到更多的人：

有这个需求，是因为我们经常对某些细胞系进行有针对性的设计变异，比如BAF155的R1064K呀，H3F3A的K27呀，那我我们拿到高通量测序数据的时候，就肯定希望可以快速的看看这个基因是否被突变成功了。

现在比对几乎不耗费什么时间了，但是得到的sam要sort的时候还是蛮耗费时间的，更耗费时间的就是全基因组时间的GATK流程了。假设，我们已经得到了所有样本的sort好的bam文件，想看看自己设计的基因突变是否成功了，可以有针对性的只call 某个基因的突变！

比如韩国人的基因组数据比对文件如下：

接近60G的数据，可以很简单的代码直接指定对某个基因进行找变异，代码如下：

grep H3F3A ~/reference/gtf/gencode/protein_coding.hg19.position

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注释一下，看看是否在自己设计的蛋白质的某个位点的氨基酸！

PS:需要自己安装annovar，可以看jimmy师兄以前的博客！

是不是看明白了呀？

生信技能树论坛