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

相关文章

来自专栏生信技能树

【直播】我的基因组59:CNV初步探索

好久不见,基因组直播又来了。这篇推送是对SNV进行一个初步探索。 单纯的一个样本来找CNV,总是不太准确的,但还是那句话,毕竟是自己的基因组,硬着头皮也要上。当...

39613
来自专栏深度学习入门与实践

【原】Spark之机器学习(Python版)(二)——分类

  写这个系列是因为最近公司在搞技术分享,学习Spark,我的任务是讲PySpark的应用,因为我主要用Python,结合Spark,就讲PySpark了。然而...

2306
来自专栏吉浦迅科技

让NVIDIA Jetson AGX Xavier火力全开的秘密

之前我们写过让Jetson TX2火力全开的秘密,让大家知道命令行工具nvpmodel能够定义一组参数,从而有效地定义给定功率的性能。

2.1K3
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 1

基因组浏览器是高通量测序分析的一个重要的可视化工具。相比于最终提供的表格,基因组浏览器可以提供更多的信息,如直观展示突变位点、查看有无新转录本或新的可变剪接形式...

1222
来自专栏Albert陈凯

分布式系统(Distributed System)资料

分布式系统(Distributed System)资料 希望转载的朋友,你可以不用联系我.但是一定要保留原文链接,因为这个项目还在继续也在不定期更新.希望看到文...

9188
来自专栏我是攻城师

如何合理的控制solr查询的命中的数量和质量?

3855
来自专栏生信技能树

【直播】我的基因组47:测序深度和GC含量的关系

在前面我们用 ChIP-seq 的分析方法可视化了一下我的 WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关。 比如在TSS...

6049
来自专栏数据小魔方

R语言可视化——图表美化与套用主题(下)

昨天的分享跟大家简单介绍了关于柱形图图表元素美化的思路,今天接着分享关于套用主题。 因为单独使用代码来调整单个图表元素,实在是太费劲了,更何况图表的细节元素有那...

2766
来自专栏生信技能树

【直播】我的基因组47:测序深度和GC含量的关系

在前面我们用 ChIP-seq 的分析方法可视化了一下我的 WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关。 比如在TSS...

37911
来自专栏听雨堂

振幅和成交量的关系

用广晟有色的历史数据,用sklearn进行回归,数据如下: ? 假设每日振幅和成交量以及价格是有关系的,于是构造: # coding=utf-8 from pa...

1988

扫码关注云+社区