前文回顾

1. GATK官方教程 / 概述及工作前的布置

2. GATK教程 / 体细胞短变异检测 (SNV+InDel)流程概览

3. GATK教程 / 变异检测前的数据预处理

4. GATK / 体细胞短变异检测工具Mutect2的使用

5. Mutect2案例 / 有或没有正常样本配对的肿瘤体细胞变异检测结果对比

问题概述

体细胞突变结果一般记录在VCF文件中，如果不可视化和绘图，人是没办法阅读的。

突变注释格式 (Mutation Annotation Format, MAF)，通常用于存储检测到的体细胞突变，并可用于后续maftools R包。从而在基因水平上观察不同肿瘤样本的突变图景，结果类似这样：

https://qiita.com/kojix2/items/738ade193e0718509c46

https://bioconductor.org/packages/release/bioc/html/maftools.html

这就需要将上游分析获得的VCF文件，注释、整理成标准的Maf格式。需要注意：

① VCF注释 与 ② 注释后整理成maf文件格式，二者可能涉及两类不同的工具(那么能否完美衔接，就是一个潜在的问题)。另外，VCF注释有时需要区分其记录的变异是胚系变异，还是体细胞变异。例如GATK Funcotator注释时，对这两类VCF文件是严格区分的 (调用了完全不同的数据库，当然基因组版本是一样的，只是涉及肿瘤驱动突变的功能注释时，调用的数据库可能有所区别)。

① VCF注释。由于VCF文件注释工具很多，如：SnpEff、ANNOVAR、GATK Funcotator、VEP、CADD等，不同软件依赖不同的运行环境和注释数据库。

注释后产生的、新的VCF文件，虽然格式上大部分一致，但在FORMAT列中 (如：AD、DP、GQ、AF)，这些子类的个数及排序各有各的不同。这种不同其实不是VCF注释带来的 (不同软件的注释结果一般都会放在INFO列，且注释内容的格式编排不同)，而是由于VCF在产生时使用的不同程序所带来的，如：GATK HaplotypeCaller、GATK Mutect2、vcftools、bcftools)。这导致已注释的VCF文件，发现其FORMAT列可能非常不同，如：

GATK HaplotypeCaller (Germline) + bcftools/vcftools

GT:AD:DP:GQ:PL  1/1:0,64:64:99:2117,192,0

GATK Mutect2 (Somatic)

GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:0,19:0.953:19:0,8:0,6:0,18:0,0,18,1

snpeffToMaf.pl脚本的示例VCF：

GT:PL:DP:AD   0/0:0,255,255:598:596,1

因此，需要使用不同的方式整理成标准的Maf文件。

② 注释后整理成maf文件格式。在整理成maftools所要求的Maf格式时，需要使用不同的工具，如snpeffToMaf.pl、GATK Funcotator (没错，此工具既能注释VCF，也能将注释结果转为Maf文件；ANNOVAR应该也可以)。

VCF to Maf的技术路径

ANNOVAR的VCF注释+Maf文件整理，网上有很多资料，我个人以前常用，后来舍弃 (非黑，原因很复杂，主要是个人习惯问题)，所以这里不再介绍。

SnpEff对VCF文件注释后，需要配合使用snpeffToMaf.pl转为MAF格式，再使用maftools绘图。但上面介绍过，不同软件获得的VCF文件的FORMAT列的内容、排序很不同，比如AF值一般只有Mutect2程序提供。因此，snpeffToMaf.pl脚本需要不断被修改，以适应不同形式的FORMAT，生成标准MAF文件。

用户自己改就比较麻烦，且容易出错，因此不建议自行修改snpeffToMaf.pl。最好是使用VCF注释工具自动输出的Maf文件。

下面是GATK Funcotator (注释数据库下载方式见文末)注释VCF、并自动输出Maf格式的参考代码 (使用了rush并行计算)：

# 注释数据库：
# funcotator_dataSources.v1.7.20200521s
cat ${result}/metadata.txt | cut -f 1 | sort -u | \rush -j 6 "gatk Funcotator -R ${ref} \-V ${result}/{}.mt2.clean.vcf.gz \--output-file-format MAF \-O ${result}/{}.mt2.clean.funcotator.maf \--data-sources-path  ~/gatk_ex/funcotator/funcotator_dataSources.v1.7.20200521s/ \--ref-version hg38"

Maf格式长什么样子？

① snpeffToMaf.pl输出的Maf文件 (表头/列的名称)：

Hugo_Symbol     Entrez_Gene_Id  Center  CDS_Change

Chromosome      Start_Position  End_Position    Strand

Variant_Classification  Variant_Type    Reference_Allele

Tumor_Seq_Allele1  Tumor_Seq_Allele2  Tumor_Sample_Barcode

Protein_Change  i_TumorVAF_WU   i_transcript_name

下面这3列，不在GATK Funcotator注释结果中；审查其内容后发现是有的，只是列的名称不同：

CDS_Change  |  i_TumorVAF_WU  |  i_transcript_name

i_TumorVAF_WU    i_transcript_name

1.00000             ENST00000335137.3

0.98824             ENST00000335137.3

0.38350             MA0139.1

0.35156             ENST00000635509.1

0.64286             ENST00000635509.1

0.34286             ENST00000429505.5

VAF：变异等位基因频率 (Variant allele frequency, VAF)

② gatk Funcotator输出的Maf文件 (表头/列的名称)：

Genome_Change/Annotation_Transcript/cDNA_Change/Codon_Change/Protein_Change

g.chr1:3728151C>T    ENST00...4.5    c.1008C>T    c.(1006-1008)gcC>gcT   p.A336A

g.chr1:3728190T>C    ENST00...4.5    c.1047T>C    c.(1045-1047)caT>caC    p.H349H

注意：下划线数字的倍数关系 (约3倍，即3联体密码子)

gatk Funcotator注释的Maf文件一共217列：

1  Hugo_Symbol

2  Entrez_Gene_Id

3  Center

4  NCBI_Build

5  Chromosome

6  Start_Position

7  End_Position

8  Strand

9  Variant_Classification

10  Variant_Type

11  Reference_Allele

12  Tumor_Seq_Allele1

13  Tumor_Seq_Allele2

14  dbSNP_RS

15  dbSNP_Val_Status

16  Tumor_Sample_Barcode

17  Matched_Norm_Sample_Barcode

18  Match_Norm_Seq_Allele1

19  Match_Norm_Seq_Allele2

20  Tumor_Validation_Allele1

21  Tumor_Validation_Allele2

22  Match_Norm_Validation_Allele1

23  Match_Norm_Validation_Allele2

24  Verification_Status

25  Validation_Status

26  Mutation_Status

27  Sequencing_Phase

28  Sequence_Source

29  Validation_Method

30  Score

31  BAM_File

32  Sequencer

33  Tumor_Sample_UUID

34  Matched_Norm_Sample_UUID

35  Genome_Change

36  Annotation_Transcript

37  Transcript_Strand

38  Transcript_Exon

39  Transcript_Position

40  cDNA_Change

41  Codon_Change

42  Protein_Change

43  Other_Transcripts

44  Refseq_mRNA_Id

45  Refseq_prot_Id

46  SwissProt_acc_Id

47  SwissProt_entry_Id

48  Description

49  UniProt_AApos

50  UniProt_Region

51  UniProt_Site

52  UniProt_Natural_Variations

53  UniProt_Experimental_Info

54  GO_Biological_Process

55  GO_Cellular_Component

56  GO_Molecular_Function

57  COSMIC_overlapping_mutations

58  COSMIC_fusion_genes

59  COSMIC_tissue_types_affected

60  COSMIC_total_alterations_in_gene

61  Tumorscape_Amplification_Peaks

62  Tumorscape_Deletion_Peaks

63  TCGAscape_Amplification_Peaks

64  TCGAscape_Deletion_Peaks

65  DrugBank

66  ref_context

67  gc_content

68  CCLE_ONCOMAP_overlapping_mutations

69  CCLE_ONCOMAP_total_mutations_in_gene

70  CGC_Mutation_Type

71  CGC_Translocation_Partner

72  CGC_Tumor_Types_Somatic

73  CGC_Tumor_Types_Germline

74  CGC_Other_Diseases

75  DNARepairGenes_Activity_linked_to_OMIM

76  FamilialCancerDatabase_Syndromes

77  MUTSIG_Published_Results

78  OREGANNO_ID

79  OREGANNO_Values

80  tumor_f

81  t_alt_count

82  t_ref_count

83  n_alt_count

84  n_ref_count

85  Gencode_34_secondaryVariantClassification

86  Achilles_Top_Genes

87  ClinVar_VCF_AF_ESP

88  ClinVar_VCF_AF_EXAC

89  ClinVar_VCF_AF_TGP

90  ClinVar_VCF_ALLELEID

91  ClinVar_VCF_CLNDISDB

92  ClinVar_VCF_CLNDISDBINCL

93  ClinVar_VCF_CLNDN

94  ClinVar_VCF_CLNDNINCL

95  ClinVar_VCF_CLNHGVS

96  ClinVar_VCF_CLNREVSTAT

97  ClinVar_VCF_CLNSIG

98  ClinVar_VCF_CLNSIGCONF

99  ClinVar_VCF_CLNSIGINCL

100  ClinVar_VCF_CLNVC

101  ClinVar_VCF_CLNVCSO

102  ClinVar_VCF_CLNVI

103  ClinVar_VCF_DBVARID

104  ClinVar_VCF_GENEINFO

105  ClinVar_VCF_MC

106  ClinVar_VCF_ORIGIN

107  ClinVar_VCF_RS

108  ClinVar_VCF_SSR

109  ClinVar_VCF_ID

110  ClinVar_VCF_FILTER

111  CosmicFusion_fusion_id

112  Familial_Cancer_Genes_Synonym

113  Familial_Cancer_Genes_Reference

114  Gencode_XHGNC_hgnc_id

115  HGNC_HGNC_ID

116  HGNC_Status

117  HGNC_Locus_Type

118  HGNC_Locus_Group

119  HGNC_Previous_Symbols

120  HGNC_Previous_Name

121  HGNC_Synonyms

122  HGNC_Name_Synonyms

123  HGNC_Chromosome

124  HGNC_Date_Modified

125  HGNC_Date_Symbol_Changed

126  HGNC_Date_Name_Changed

127  HGNC_Accession_Numbers

128  HGNC_Enzyme_IDs

129  HGNC_Ensembl_Gene_ID

130  HGNC_Pubmed_IDs

131  HGNC_RefSeq_IDs

132  HGNC_Gene_Family_ID

133  HGNC_Gene_Family_Name

134  HGNC_CCDS_IDs

135  HGNC_Vega_ID

136  HGNC_OMIM_ID(supplied_by_OMIM)

137  HGNC_RefSeq(supplied_by_NCBI)

138  HGNC_UniProt_ID(supplied_by_UniProt)

139  HGNC_Ensembl_ID(supplied_by_Ensembl)

140  HGNC_UCSC_ID(supplied_by_UCSC)

141  Oreganno_Build

142  Simple_Uniprot_alt_uniprot_accessions

143  dbSNP_ASP

144  dbSNP_ASS

145  dbSNP_CAF

146  dbSNP_CDA

147  dbSNP_CFL

148  dbSNP_COMMON

149  dbSNP_DSS

150  dbSNP_G5

151  dbSNP_G5A

152  dbSNP_GENEINFO

153  dbSNP_GNO

154  dbSNP_HD

155  dbSNP_INT

156  dbSNP_KGPhase1

157  dbSNP_KGPhase3

158  dbSNP_LSD

159  dbSNP_MTP

160  dbSNP_MUT

161  dbSNP_NOC

162  dbSNP_NOV

163  dbSNP_NSF

164  dbSNP_NSM

165  dbSNP_NSN

166  dbSNP_OM

167  dbSNP_OTH

168  dbSNP_PM

169  dbSNP_PMC

170  dbSNP_R3

171  dbSNP_R5

172  dbSNP_REF

173  dbSNP_RV

174  dbSNP_S3D

175  dbSNP_SAO

176  dbSNP_SLO

177  dbSNP_SSR

178  dbSNP_SYN

179  dbSNP_TOPMED

180  dbSNP_TPA

181  dbSNP_U3

182  dbSNP_U5

183  dbSNP_VC

184  dbSNP_VP

185  dbSNP_WGT

186  dbSNP_WTD

187  dbSNP_dbSNPBuildID

188  dbSNP_ID

189  dbSNP_FILTER

190  HGNC_Entrez_Gene_ID(supplied_by_NCBI)

191  dbSNP_RSPOS

192  dbSNP_VLD

193  AS_FilterStatus

194  AS_SB_TABLE

195  AS_UNIQ_ALT_READ_COUNT

196  CONTQ

197  DP

198  ECNT

199  GERMQ

200  MBQ

201  MFRL

202  MMQ

203  MPOS

204  NALOD

205  NCount

206  NLOD

207  OCM

208  PON

209  POPAF

210  ROQ

211  RPA

212  RU

213  SEQQ

214  STR

215  STRANDQ

216  STRQ

217  TLOD

结论：gatk Funcotator 算出来的 tumor_f 列不完全等于 snpeffToMaf.pl 算出来的 i_TumorVAF_WU 列；虽然都可能是由AD (Allelic depths for the ref and alt alleles in the order listed)推导而来 (但在肿瘤体细胞突变中不完全等于：alt/(alt+ref)，如下图)。

6/36=0.1666667  |  6/42=0.1428571

4/39=0.1025641  |  4/43 [1] 0.09302326

9/70=0.1285714  |  9/79=0.1139241

158/(158+254)=0.3834951  |  9/14=0.6428571

(附)GATK Funcotator注释数据库下载

GATK包含两组预打包的数据源，允许在无需 (大量)额外配置的情况下使用Functator。这些数据源包分别对应于胚系和体细胞突变情形。

广义地说，如果你有胚系VCF，胚系数据源就是你想要开始使用的。相反，如果您有一个体细胞VCF，那么体细胞数据源就是您想要开始使用的数据源。

注释数据库下载方式 (2选1)：

1.ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/

2.gs://broad-public-datasets/funcotator/

gsutil ls gs://broad-public-datasets/

gsutil ls gs://broad-public-datasets/funcotator/

注意：上图最下方g和s的区别！

① 获取注释数据库的胚系数据源：

gsutil -m cp -r gs://broad-public-datasets/funcotator/funcotator_dataSources.v1.7.20200521g/ \funcotator/
# annotation
# --output-file-format The output file format. Either VCF, MAF, or SEG. # Please note that MAF output for germline use case VCFs is unsupported.# SEG will generate two output files: a simple tsv and a gene list. Required. Possible values: {VCF, MAF, SEG}

下载结果 (3.6 GiB)

② 获取注释数据库的体细胞数据源：

  同上，但小写 g 改为小写 s

查看README，发现二者用处完全不同

Operation completed over 81 objects/34.4 GiB.

查看、比较下载后的文件大小

另外，gsutil也可以使用一些bash命令，例如：du -sh

资料来源

Tool Documentation Index

https://gatk.broadinstitute.org/hc/en-us/articles/5358824293659--Tool-Documentation-Index#Funcotator

https://gatk.broadinstitute.org/hc/en-us/articles/5358830252187-Funcotator

https://gatk.broadinstitute.org/hc/en-us/articles/360035889931-Funcotator-Information-and-Tutorial

https://mp.weixin.qq.com/s/Y-_aDGTNUZVW-QMbhZDHXA

https://qiita.com/kojix2/items/738ade193e0718509c46

https://bioconductor.org/packages/release/bioc/html/maftools.html

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集