把maf格式的somatic突变数据导入annovar去除人群频率变异 / 开普饭

首先maf格式的somatic突变数据制作成为annovar软件的输入格式：

cut -f 5-7,12,13,1,16 human_brca_all_mutect2.maf |cut -f 2-7 > 1 cut -f 5-7,12,13,1,16 human_brca_all_mutect2.maf |cut -f 1 > 2 paste 1 2 > for_annovar.input ### 共 13027 位点

然后运行annovar软件的批量注释功能

bin=/home/haitaowang/Database/hg38/Annovar/ db=/home/haitaowang/Database/hg38/Annovar/ perl $bin/table_annovar.pl for_annovar.input $db -buildver hg38 -out tmp \ -protocol refGene,cytoBand,esp6500siv2_all,exac03,gnomad_genome,cosmic70,1000g2015aug_all,clinvar_20170905 \ -operation g,r,r,f,f,f,f,f -nastring NA

因为数据库较多，所以注释耗时很长，注释结果如下：

551K Aug 17 16:52 for_annovar.input 99K Aug 17 17:10 tmp.hg38_ALL.sites.2015_08_dropped 423K Aug 17 17:10 tmp.hg38_ALL.sites.2015_08_filtered 459K Aug 17 16:57 tmp.hg38_avsnp150_dropped 174K Aug 17 16:57 tmp.hg38_avsnp150_filtered 56K Aug 17 17:11 tmp.hg38_clinvar_20170905_dropped 477K Aug 17 17:11 tmp.hg38_clinvar_20170905_filtered 35K Aug 17 17:10 tmp.hg38_cosmic70_dropped 470K Aug 17 17:10 tmp.hg38_cosmic70_filtered 667K Aug 17 16:53 tmp.hg38_cytoBand 366K Aug 17 17:10 tmp.hg38_dbnsfp33a_dropped 445K Aug 17 17:10 tmp.hg38_dbnsfp33a_filtered 115K Aug 17 16:53 tmp.hg38_esp6500siv2_all_dropped 410K Aug 17 16:53 tmp.hg38_esp6500siv2_all_filtered 423K Aug 17 16:53 tmp.hg38_exac03_dropped 307K Aug 17 16:53 tmp.hg38_exac03_filtered 286K Aug 17 16:53 tmp.hg38_genomicSuperDups 731K Aug 17 16:55 tmp.hg38_gnomad_genome_dropped 178K Aug 17 16:55 tmp.hg38_gnomad_genome_filtered 153K Aug 17 17:11 tmp.hg38_intervar_20180118_dropped 436K Aug 17 17:11 tmp.hg38_intervar_20180118_filtered 40K Aug 17 17:12 tmp.hg38_mcap_dropped 457K Aug 17 17:12 tmp.hg38_mcap_filtered 6.3M Aug 17 17:12 tmp.hg38_multianno.txt 48K Aug 17 17:12 tmp.hg38_revel_dropped 447K Aug 17 17:12 tmp.hg38_revel_filtered 68K Aug 17 17:12 tmp.invalid_input 956 Aug 17 17:12 tmp.log 208K Aug 17 16:53 tmp.refGene.exonic_variant_function 68K Aug 17 16:53 tmp.refGene.invalid_input 1.2K Aug 17 16:53 tmp.refGene.log 761K Aug 17 16:53 tmp.refGene.variant_function

简单数一下：

1454 tmp.hg38_ALL.sites.2015_08_dropped 7400 tmp.hg38_avsnp150_dropped 170 tmp.hg38_clinvar_20170905_dropped 326 tmp.hg38_cosmic70_dropped 937 tmp.hg38_dbnsfp33a_dropped 1800 tmp.hg38_esp6500siv2_all_dropped 4262 tmp.hg38_exac03_dropped 7305 tmp.hg38_gnomad_genome_dropped 1167 tmp.hg38_intervar_20180118_dropped 648 tmp.hg38_mcap_dropped 890 tmp.hg38_revel_dropped

需要一个个数据库来解读。

首先，被千人基因组计划的人群频率0.05过滤掉的坐标拿出来：

perl -alne '{print if $F[1]>0.05}' tmp.hg38_ALL.sites.2015_08_dropped > filter_by_1000g.pos

有了这些坐标，就可以去过滤我们的原始maf文件啦。

cat filter_by_1000g.pos ../human_brca_all_mutect2.maf |perl -alne '{if(/^1000/){$h{"$F[2]\t$F[3]"}=1}else{print unless exists $h{"$F[4]\t$F[5]"}}}' > filter_by_1000g.maf

同理，其它数据库也是同样的操作

perl -alne '{print if (split(",",$F[1]))[0]>0.05}' tmp.hg38_gnomad_genome_dropped > filter_by_gnomad.pos perl -alne '{print if (split(",",$F[1]))[0]>0.05}' tmp.hg38_exac03_dropped > filter_by_exac03.pos cat filter_by_exac03.pos filter_by_1000g.maf |perl -alne '{if(/^exac03/){$h{"$F[2]\t$F[3]"}=1}else{print unless exists $h{"$F[4]\t$F[5]"}}}' > filter_by_1000g_exac03.maf cat filter_by_gnomad.pos filter_by_1000g_exac03.maf|perl -alne '{if(/^gnomad_genome/){$h{"$F[2]\t$F[3]"}=1}else{print unless exists $h{"$F[4]\t$F[5]"}}}' > filter_by_1000g_exac03_gnomeAD.maf

过滤效果如下：

wc -l *.maf 9955 filter_by_1000g_exac03_gnomeAD.maf 10588 filter_by_1000g_exac03.maf 12141 filter_by_1000g.maf

过滤后的maf作为肿瘤外显子分析结果应该是会更合理。

生信技能树GATK4系列教程

GATK4的gvcf流程

你以为的可能不是你以为的

新鲜出炉的GATK4培训教材全套PPT，赶快下载学习吧

曾老湿最新私已：GATK4实战教程

GATK4的CNV流程-hg38

然后是 CNV相关工具

WES的CNV探究-conifer软件使用

单个样本NGS数据如何做拷贝数变异分析呢