GATK下游VCF maf可视化

VCF可以通过不同的上游分析获得，在格式上又有一些不同，盲目转化会带来很多问题。查看VCF文件内容可知是通过GATK上游处理得到的。

less raw.vcf

##fileformat=VCFv4.1
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##GATKCommandLine.HaplotypeCaller=

为了对变异结果进一步进行可视化分析，查找了一下网上的相关资料，可以通过annovar软件进行转化，然后通过maftools R包进行可视化分析。确定了整体的分析策略之后，先安装annovar软件。

1.安装annovar

官网下载：http://annovar.openbioinformatics.org/en/latest/user-guide/download/
注册 网站发送邮件 直接邮件下载 拷贝到服务器 （需使用机构邮箱）。
拷贝到服务器，在相应的文件夹下进行解压即可。

tar -zxvf annovar.latest.tar.gz

2.VCF转换为annovar格式

### 批量转换
for i in *vcf; do perl ~/biosoft/annovar/convert2annovar.pl -format vcf4 $i -outfile $i.annovar -includeinfo; done

3.VCF变异根据基因进行注释

for i in *vcf; do \
perl ~/biosoft/annovar/table_annovar.pl $i.annovar ~/path/to/humandb/ -buildver hg19 
-out variants --otherinfo -remove -protocol ensGene -operation g -nastring NA

这个时候会遇到一个报错：

hg19_ensGeneMrna.fa does not exist. Please use 'annotate_variation.pl --downdb ensGene ~/biosoft/annovar/humandb' download the database

按照提示下载相应的注释参考数据库，提示相应的ref文件已经无法下载。重新构建ref

perl ./annotate_variation.pl --buildver hg19 --downdb seq ~/biosoft/annovar/humandb/hg19_seq
perl ./retrieve_seq_from_fasta.pl /home/data/ssy066/biosoft/annovar/humandb/hg19_ensGene.txt -seqdir ~/biosoft/annovar/humandb/hg19_seq -format ensGene -outfile ~/biosoft/annovar/humandb/hg19_ensGeneMrna.fa

重新进行注释。

4.转化为maf格式文件

目标文件夹下会出现后缀为multianno.txt的文件，导入R语言中转化为maf文件，进行可视化即可。、

library(maftools)
var.annovar.maf = annovarToMaf(annovar = "~/BN4027P6.raw.vcf.anno.hg19_multianno.txt", 
                               Center = 'CSI-NUS', refBuild = 'hg19', MAFobj = T, table = 'ensGene')
plotmafSummary(maf = var.annovar.maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)

参考内容

1.maftools : Summarize, Analyze and Visualize MAF Files (bioconductor.org)