肿瘤外显子数据处理系列教程（六）vcf文件的注释及ANNOVAR的使用
大家好，非常抱歉，上一周因为事务繁忙来不及更新，这周我们继续。
上一节我们讲到了，GATK的最佳实践，并且拿到了mutect2的vcf文件。这一周，我们先对vcf文件结合bam文件进行一定的探索，然后对其进行注释，顺便学习一下王凯老师的ANNOVAR软件的使用。
vcf载入IGV
同样的，vcf文件要载入igv进行可视化之前也要构建索引，构建的方法很简单，以case1_biorep_A_techrep_filter.vcf为例：
bgzip -c case1_biorep_A_techrep_filter.vcf >case1_biorep_A_techrep_filter.vcf.gz
bcftools index -t case1_biorep_A_techrep_filter.vcf.gz
然后将vcf.gz文件及其索引下载到本地，存放在同一个目录下，直接把vcf.gz文件载入igv。还有我们之前番外篇的小bam文件small.bam也一起载入，定位到17号，因为我们之前提取的小bam文件是17号染色体的。
我们从vcf文件中随机挑选一个突变位点进行查看，如：
CHROM POS ID REF ALT QUAL FILTER INFO FORMAT case1_biorep_A_techrep case1_germline
chr17 78526659 . G C . PASS DP=357;ECNT=1;NLOD=34.92;N_ART_LOD=-2.076e+00;POP_AF=1.000e-05;P_CONTAM=0.00;P_GERMLINE=-5.877e+01;TLOD=85.25 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:SA_MAP_AF:SA_POST_PROB 0/1:187,34:0.176:84,19:103,15:40,30:150,136:60:18:0.152,0.121,0.154:3.206e-03,0.076,0.921 0/0:116,0:0.021:59,0:57,0:37,0:151,0:0:0
在igv中搜索框输入对应的坐标即可定位到该位点，并将光标定在图示位置，即可查看该位点的详细信息，与vcf文件中记录的信息相对应，同时我们可以看到，在bam文件中的比对结果：
为了更加直观的看突变位点，我们在bam文件的reads位置右键选择按照碱基排序：
排序后可以看到，在参考基因组hg38上，改位点的碱基是G，而在该样本中是部分为C，有一条reads为A，一条reads为N，后两者测到的reads数太少，可信度都比较低，因此在vcf文件中没有被记录。感兴趣的朋友还可以进一步探索vcf文件和bam文件之间的关系，我们继续下一步。
ANNOVAR的使用
关于ANNOVAR的使用，其实曾老师写过很多教程，如：ANNOVAR软件用法还可以更复杂：http://www.bio-info-trainee.com/4007.html，这里我们就不再展开了，直接使用就好。
这里我们直接采用联合注释的方法：假如只有一个vcf文件，那可以用下面的方法做注释：
~/biosoft/ANNOVAR/annovar/table_annovar.pl case1_biorep_A_techrep_filter.vcf ~/biosoft/ANNOVAR/annovar/humandb/ \
-buildver hg38 \
-out case1_biorep_A_techrep_filter \
-remove \
-protocol knownGene,clinvar_20170905 \
-operation g,f \
-nastring . \
-vcfinput
但是我们这里有多个样本，所以要做批处理，同样的我们需要一个config配置文件，如下：
cat config
case1_biorep_A_techrep_filter
case1_biorep_B_filter
case1_biorep_C_filter
case1_techrep_2_filter
......
case6_biorep_A_techrep_filter
case6_biorep_B_filter
case6_biorep_C_filter
case6_techrep_2_filter
批处理的脚本是：
cat annovar.sh
cat config|while read id
do
~/biosoft/ANNOVAR/annovar/table_annovar.pl ${id}.vcf ~/biosoft/ANNOVAR/annovar/humandb/ \
-buildver hg38 \
-out ${id} \
-remove \
-protocol knownGene,clinvar_20170905 \
-operation g,f \
-nastring . \
-vcfinput
done
注释结果
每个样本的输出结果有3个文件，我们关心的是下面两个文件，同样以这个样本为例case1_biorep_A_techrep：
case1_biorep_A_techrep_filter.hg38_multianno.txt
case1_biorep_A_techrep_filter.hg38_multianno.vcf
第一个文件就记录了突变注释信息，如突变位点所在的基因、是nonsynonymous或者synonymous等
Chr Start End Ref Alt Func.knownGene Gene.knownGene GeneDetail.knownGene ExonicFunc.knownGene AAChange.knownGene CLINSIG CLNDBN CLNACC CLNDSDB CLNDSDBID Otherinfo
chr1 6146376 6146376 G T exonic CHD5 . nonsynonymous SNV CHD5:uc001amb.3:exon11:c.C1638A:p.N546K,CHD5:uc057btb.1:exon11:c.C1638A:p.N546K . . . . . 0.25 . 75 chr1 6146376 . G T . PASS DP=184;ECNT=1;NLOD=22.56;N_ART_LOD=-1.888e+00;POP_AF=1.000e-05;P_CONTAM=0.00;P_GERMLINE=-4.514e+01;TLOD=5.41 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:SA_MAP_AF:SA_POST_PROB 0/1:97,3:0.058:48,1:49,2:37,38:154,210:60:20:0.020,0.030,0.030:0.018,3.223e-03,0.979 0/0:75,0:0.020:45,0:30,0:36,0:154,0:0:0
chr1 6461445 6461445 G T exonic TNFRSF25 . nonsynonymous SNV TNFRSF25:uc001ana.4:exon5:c.C694A:p.R232S,TNFRSF25:uc001anf.4:exon9:c.C1132A:p.R378S,TNFRSF25:uc001ang.4:exon9:c.C1108A:p.R370S,TNFRSF25:uc001ane.4:exon10:c.C1243A:p.R415S,TNFRSF25:uc001anh.4:exon10:c.C1270A:p.R424S . . . . . 0.25 . 9
chr1 6461445 . G T . PASS DP=21;ECNT=1;NLOD=2.70;N_ART_LOD=-1.049e+00;POP_AF=1.000e-05;P_CONTAM=0.00;P_GERMLINE=-1.853e+00;TLOD=6.33 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:SA_MAP_AF:SA_POST_PROB 0/1:6,2:0.278:3,2:3,0:29,39:204,184:60:14:0.00,0.253,0.250:0.023,0.025,0.952 0/0:9,0:0.126:3,0:6,0:24,0:169,0:0:0
chr1 31756671 31756671 C A exonic BAI2 . nonsynonymous SNV BAI2:uc057ego.1:exon1:c.G130T:p.D44Y,BAI2:uc057egq.1:exon1:c.G145T:p.D49Y,BAI2:uc001btq.2:exon2:c.G130T:p.D44Y,BAI2:uc010ogn.3:exon2:c.G175T:p.D59Y,BAI2:uc010ogq.3:exon2:c.G166T:p.D56Y,BAI2:uc031twa.2:exon2:c.G268T:p.D90Y,BAI2:uc057egm.1:exon2:c.G130T:p.D44Y,BAI2:uc001btn.4:exon4:c.G166T:p.D56Y,BAI2:uc001bto.4:exon4:c.G166T:p.D56Y . . . . . 0.25 . 33
......
简单统计一下可以发现，总共有236个突变位点，235位点都在外显子上，然而，却有一个位点是在intergenic，有点诡异，我们后面有时间再讨论：
cat case1_biorep_A_techrep_filter.hg38_multianno.txt|grep exonic|less -S|wc
# 235 7036 122942
cat case1_biorep_A_techrep_filter.hg38_multianno.txt|grep -v exonic
# Chr Start End Ref Alt Func.knownGene Gene.knownGene GeneDetail.knownGene ExonicFunc.knownGene AAChange.knownGene CLINSIG CLNDBN CLNACC CLNDSDB CLNDSDBID Otherinfo
# chr19 22857426 22857426 G - intergenic ZNF99,CTD-2579N5.3 dist=73319;dist=49660 . . . . . . . 0.25 . 60 chr19 22857425 . TG T . PASS DP=271;ECNT=1;NLOD=18.05;N_ART_LOD=-1.794e+00;POP_AF=1.000e-05;P_CONTAM=0.00;P_GERMLINE=-2.893e+01;RPA=2,1;RU=G;STR;TLOD=124.43 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:SA_MAP_AF:SA_POST_PROB 0/1:148,46:0.255:67,21:81,25:33,33:153,153:60:18:0.222,0.222,0.237:0.013,9.937e-03,0.977 0/0:60,0:0.024:31,0:29,0:32,0:148,0:0:0
下一周我们将用maftools来对ANNOVAR的注释结果进一步分析。
但是，这个教程并没有回答我们生信技能树今天推出的那个问题，不同转录组流程结果到底该如何比较，里面有TP53这个基因的某个蛋白的杂合突变 R280T ，需要找到坐标，首先需要知道转录本，有一个网页工具可以继续转换，不知道这个ANNOVAR能否做到？
如果你完全看不懂这个系列教程，那么下面的课程你可能会需要！
1
全国巡讲第17站
2
全国巡讲第18站
１ | 生信－R语言入门 |
２ | GEO数据库挖掘 |
３ | 生信－LINUX基础 |
如果你精选10篇我们生信技能树2019对你帮助最大的推文教程，发到我邮箱 jmzeng1314@163.com 并且写出你的故事，就有惊喜哦！