【直播】我的基因组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

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-24

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏一个爱瞎折腾的程序猿

sqlserver使用存储过程跟踪SQL

USE [master] GO /****** Object: StoredProcedure [dbo].[sp_perfworkload_trace_s...

2060
来自专栏芋道源码1024

熔断器 Hystrix 源码解析 —— 断路器 HystrixCircuitBreaker

本文主要基于 Hystrix 1.5.X 版本 1. 概述 2. HystrixCircuitBreaker 3. HystrixCircuitBreaker....

5327
来自专栏大内老A

The .NET of Tomorrow

Ed Charbeneau(http://developer.telerik.com/featured/the-net-of-tomorrow/) Exciti...

31710
来自专栏Ceph对象存储方案

Luminous版本PG 分布调优

Luminous版本开始新增的balancer模块在PG分布优化方面效果非常明显,操作也非常简便,强烈推荐各位在集群上线之前进行这一操作,能够极大的提升整个集群...

3145
来自专栏跟着阿笨一起玩NET

c#实现打印功能

2762
来自专栏飞扬的花生

jsencrypt参数前端加密c#解密

      写程序时一般是通过form表单或者ajax方式将参数提交到服务器进行验证,如何防止提交的请求不被抓包后串改,虽然无法说绝对安全却给非法提交提高了难度...

3869
来自专栏陈仁松博客

ASP.NET Core 'Microsoft.Win32.Registry' 错误修复

今天在发布Asp.net Core应用到Azure的时候出现错误InvalidOperationException: Cannot find compilati...

4848
来自专栏pangguoming

Spring Boot集成JasperReports生成PDF文档

由于工作需要,要实现后端根据模板动态填充数据生成PDF文档,通过技术选型,使用Ireport5.6来设计模板,结合JasperReports5.6工具库来调用渲...

1.2K7
来自专栏魂祭心

原 canvas绘制clock

4074
来自专栏落花落雨不落叶

canvas画简单电路图

61811

扫码关注云+社区