【直播】我的基因组81：看看我的vcf文件的vaf分布情况

这一讲中，我们对VCF中的"VAF"简单的来看一起，如果你对VCF文件还不了解的话，那你就要自我批评一下了。在基因组直播刚开始的时候，我还专门对VCF文件进行了简述。【直播】我的基因组28-必须要理解vcf格式记录的变异位点信息. 今天不说别的，我们专门对看一下VAF的分布情况。

VAF",就是variant allele frequency 或者 variant allele fraction

对于NGS测序数据来说，就是跟参考基因不同的reads与总的测序reads的比值。

一般在VCF文件里面，会有DP4这个信息，可以很容易算出vaf值，如下;

得到上面数据的代码是：

首先是shell

cat  autochr.highQuali.varType | perl -alne '{next if /^#/;/DP4=(.*?);.*VARTYPE=(.*?)s/;print "$F[0],$1,$2"}'>DP4.stat

然后是R

a=read.csv('DP4.stat',stringsAsFactors = F,header = F)
colnames(a)=c('chr','ref_f','ref_r','alt_f','alt_r','type')
a=a[grepl('chr',a[,1]),]
## Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases
head(a)
#lapply(2:5, function(i){ a[,i]=as.numeric(a[,i]) })
a[,2]=as.numeric(a[,2])
a[,3]=as.numeric(a[,3])
a[,4]=as.numeric(a[,4])
a[,5]=as.numeric(a[,5])
a$vaf=(a$alt_f+a$alt_r)/(a$alt_f+a$alt_r+a$ref_f+a$ref_r )
table(a[,c(1,6)])
snp=subset(a,type=='SNP')
head(snp)
hist(snp$vaf)
indel=subset(a,type!='SNP')
head(indel)
hist(indel$vaf)

正常人的二倍体基因组位点只有杂合或者纯合两种情况，对于纯合那么vaf必然是1，对于杂合，必然是0.5。但是现实测序得到的结果远比这要复杂，尤其是测序深度不够的时候。因为测序本身具有随机性，而且还有很多系统误差。理想情况也只能像是扔硬币。

我的vcf文件里面所有的snp突变位点的vaf值分布如下：

可以看到纯合位点和杂合位点有一个很明显的分界线，就是我们通常说的二八原则咯。

对杂合位点来说，理论上跟扔硬币一样，是概率事件。

还有，我的vcf文件里面所有的indel突变位点的vaf值分布如下：

一般来说，DP4只要比对之后很容易从bam文件里面算出来(samtools mpileup命令即可)，其实最好的情况下不需要各种call variation的软件了，简单的判断语句就知道各个位置是不是变异了，是纯合呢还是杂合。但是我们说过，实际的测序结果往往是很复杂的，在很多位点，普通的判断语句并不适用，即使是主流的variation caller的表现也往往不能统一。

而文献里面对TCGA里面的癌症样本的somatic mutation的vaf统计如下：

Figure 7 Distribution of the Variant Allele Fraction (VAF) of somatic mutations in one sample of lung adenocarcinoma from the TCGA study .

文章来源：Computational methods and resources for the

interpretation of genomic variants in cancer

可以看出tumor里面的vaf分布其实已经不再是扔硬币那样的概率了，对于杂合位点来说。

原因很多，首先tumor不一定是单纯的二倍体了，其次tumor样品一般来说本身异质性高，而我们测序是混合多个细胞的，有一些突变有一些并不突变。而且纯合的somatic mutation几乎没有，因为somatic mutation是tumor过滤了normal后留下来的变异位点，不是遗传多样性，突变这个过程既然是后天产生的，就很难保证取样部分的几百万个细胞全部突变了。

With cancer data it is important to look at the allele frequency in the sample. Most cancer samples are a mixture of non-cancerous cells mixed with cancerous cells that are clonal expansions of beneficial mutations (to the cancer). So, as you say, a 0.5 frequency indicates the site is heterozygous in the individual. A lower frequency might suggest it is a tumor mutation that has swept through the tumor cells, and an even lower frequency suggests it is a clonal subpopulation.

http://www.nature.com/nature/journal/v481/n7381/full/nature10762.html