【直播】我的基因组82:如何对maf格式的突变文件统计vaf
前面我们已经详细讲解了vaf,而且也对我的vcf文件进行了统计。【直播】我的基因组81:看看我的vcf文件的vaf分布情况
文末我提到了一篇文章里面画出了TCGA的一个tumor样本的somatic mutation的vaf分布!
后台很多小伙伴留言那个该如何做。
其实非常简单的。
前面我们介绍了TCGA的maf文件的下载。
我们随便看一个癌种吧!
先在shell里面
cut -f 40-42 TCGA.BLCA.mutect.somatic.maf >TCGA.BLCA.counts
然后在R里面:
a=read.table("TCGA.BLCA.counts",header = T,stringsAsFactors = F)
a$vaf=a$t_alt_count/a$t_depth
hist(a$vaf,breaks = 0:100/100)
就出图啦:
虽然很粗糙,但是意思已经明白了。
但是可以看出tumor里面的vaf分布其实已经不再是扔硬币那样的概率了,对于杂合位点来说。
原因很多,首先tumor不一定是单纯的二倍体了,其次tumor样品一般来说本身异质性高,而我们测序是混合多个细胞的,有一些突变有一些并不突变。
而且纯合的somatic mutation几乎没有,因为somatic mutation是tumor过滤了normal后留下来的变异位点,不是遗传多样性,突变这个过程既然是后天产生的,就很难保证取样部分的几百万个细胞全部突变了。
为什么取对TCGA下载的maf文件取第40~42列呢,也是经过搜索学习的。
MAF格式定义如下:
https://software.broadinstitute.org/software/igv/MutationAnnotationFormat
https://wiki.nci.nih.gov/display/TCGA/Mutation+Annotation+Format+%28MAF%29+Specification
但是它只规定了前34列,后面所有的列,都是可以自由扩展的。
而我们下载的maf文件共有114 columns,如下:
Hugo_Symbol
Entrez_Gene_Id
Center
NCBI_Build
Chromosome
Start_Position
End_Position
Strand
Variant_Classification
Variant_Type
Reference_Allele
Tumor_Seq_Allele1
Tumor_Seq_Allele2
dbSNP_RS
dbSNP_Val_Status
Tumor_Sample_Barcode
Matched_Norm_Sample_Barcode
Match_Norm_Seq_Allele1
Match_Norm_Seq_Allele2
Tumor_Validation_Allele1
Tumor_Validation_Allele2
Match_Norm_Validation_Allele1
Match_Norm_Validation_Allele2
Verification_Status
Validation_Status
Mutation_Status
Sequencing_Phase
Sequence_Source
Validation_Method
Score
BAM_File
Sequencer
Tumor_Sample_UUID
Matched_Norm_Sample_UUID
HGVSc
HGVSp
HGVSp_Short
Transcript_ID
Exon_Number
t_depth
t_ref_count
t_alt_count
n_depth
n_ref_count
n_alt_count
all_effects
Allele
Gene
Feature
Feature_type
Consequence
cDNA_position
CDS_position
Protein_position
Amino_acids
Codons
Existing_variation
ALLELE_NUM
DISTANCE
TRANSCRIPT_STRAND
SYMBOL
SYMBOL_SOURCE
HGNC_ID
BIOTYPE
CANONICAL
CCDS
ENSP
SWISSPROT
TREMBL
UNIPARC
RefSeq
SIFT
PolyPhen
EXON
INTRON
DOMAINS
GMAF
AFR_MAF
AMR_MAF
ASN_MAF
EAS_MAF
EUR_MAF
SAS_MAF
AA_MAF
EA_MAF
CLIN_SIG
SOMATIC
PUBMED
MOTIF_NAME
MOTIF_POS
HIGH_INF_POS
MOTIF_SCORE_CHANGE
IMPACT
PICK
VARIANT_CLASS
TSL
HGVS_OFFSET
PHENO
MINIMISED
ExAC_AF
ExAC_AF_AFR
ExAC_AF_AMR
ExAC_AF_EAS
ExAC_AF_FIN
ExAC_AF_NFE
ExAC_AF_OTH
ExAC_AF_SAS
GENE_PHENO
FILTER
src_vcf_id
tumor_bam_uuid
normal_bam_uuid
GDC_Validation_Status
GDC_Valid_Somatic
而能够统计vaf的列,我已经标记了红色,就是40~42列咯。