【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏数据库

速来围观!——三种NCBI常见数据库

在微生物测序分析中,常常需要对未知的核酸或蛋白序列进行物种,功能或类别注释。注释方法种类较多,其中最常用的是与一些标准数据库进行相似性搜索,也就是序列比对。因此...

19310
来自专栏生信技能树

Variant 分析阶段小结1-基础碎碎念

所谓遗传变异是生物体内遗传物质发生变化而造成的可以遗传给后代的变异,这些变异导致了生物在不同水品上体现出遗传的多样性。生物信息学中各种基因组研究的基础就是遗传变...

1163
来自专栏生信技能树

【直播】我的基因组78:简单解析一下蛋白编码基因的测序深度及覆盖度

上一讲中,我们对蛋白的编码基因的测序深度和覆盖度进行了统计,其中有的覆盖度很高,有的覆盖度却又很低,针对这个统计出的测序深度及覆盖度,我们就可以做一些简单的统计...

3587
来自专栏生信技能树

【直播】我的基因组79:为什么这些基因的覆盖度如此之低?

在之前,我们计算了每个基因的GC含量以及基因长度,也要samtools计算了每个基因的覆盖度以及平均测序深度,还有基因的内部测序深度差异值(S值)。也在上一讲中...

26810
来自专栏生信技能树

单细胞转录组3大R包之monocle2

主要是针对单细胞转录组测序数据开发的,用来找不同细胞类型或者不同细胞状态的差异表达基因。分析起始是表达矩阵,作者推荐用比较老旧的Tophat+Cufflinks...

1.6K9
来自专栏Y大宽

ToppGene Suite中文使用指南

2007.12:Improved human disease candidate gene prioritization using mouse phenoty...

923
来自专栏大数据挖掘DT机器学习

Python NLTK自然语言处理:词干、词形与MaxMatch算法

CSDN:白马负金羁 自然语言处理是计算机科学领域与人工智能领域中的一个重要方向。自然语言工具箱(NLTK,Natural Language Toolkit)...

4175
来自专栏生信技能树

第1篇:ATAC-seq的背景介绍以及与ChIP-Seq的异同

ATAC-seq(Assay for Transposase-Accessible Chromatin with high throughput sequenc...

1182
来自专栏数说戏聊

11.散点图&折线图&饼图1.散点图2.折线图饼图

以一个变量为横坐标,另一个变量为纵坐标,利用散点(坐标点)的分布形态反映变量关系的图形。

661
来自专栏新智元

【TensorFlow1.2.0版发布】14大新功能,增加Intel MKL集成

【新智元导读】TensorFlow 今天发布最新版 1.2.0,公布了14大最新功能。新智元带来最新介绍,包括 API 的重要变化、contrib API的变化...

3249

扫码关注云+社区