肿瘤外显子数据处理系列教程（六）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对你帮助最大的推文教程，发到我邮箱 [e-mail] 并且写出你的故事，就有惊喜哦！