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

相关文章

来自专栏斑斓

大数据 | Spark中实现基础的PageRank

吴军博士在《数学之美》中深入浅出地介绍了由Google的佩奇与布林提出的PageRank算法,这是一种民主表决式网页排名技术。书中提到PageRank的核心思想...

3128
来自专栏生信技能树

点点鼠标就能完成的NMR代谢组学实战

代谢组学几乎完全不涉及生物信息学最核心的序列比对,包括武汉中科院数物所波谱国重实验室等主流科研机构都是利用化学计量学和多元统计分析方法,对通过核磁共振(NMR)...

763
来自专栏生信技能树

使用ESTIMATE来对转录组表达数据根据stromal和immune细胞比例估算肿瘤纯度

ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues usin...

1072
来自专栏专知

【专知-Deeplearning4j深度学习教程02】用ND4J自己动手实现RBM: 图文+代码

【导读】主题链路知识是我们专知的核心功能之一,为用户提供AI领域系统性的知识学习服务,一站式学习人工智能的知识,包含人工智能( 机器学习、自然语言处理、计算机视...

40510
来自专栏生信技能树

【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapmap计划数据比较

hapmap计划的人群分布结果和千人基因组计划的分布结果来分析是一样的!【直播】我的基因组55:简单的PCA分析千人基因组的人群分布 这两个计划里面收集的样本的...

3926
来自专栏生信技能树

【直播】我的基因组48:我可能测了一个假的全基因组

背景知识 男性只有一条X染色体和一条Y染色体,所以,理论上它们上面的SNV都应该是纯合的! X,Y除了同源区域外,其它地方差异很大。所以在女性样本里面即使是混入...

27912
来自专栏生信宝典

2018 升级版Jaspar数据库

R包ggseqlogo 绘制seq logo图和Seq logo 在线绘制工具—Weblogo介绍了如何用R脚本和在线工具绘制seq logo图,用于展现转录因...

692
来自专栏生信宝典

NGS基础 - FASTQ格式解释和质量评估

FASTQ文件格式和命名 高通量测序之后用于下游分析的数据一般存储在FASTQ文件中。为了节省空间,又不影响下游使用,也一般用gzip压缩的格式。 单端测序每个...

2075
来自专栏冰霜之地

Threes-AI 玩小三传奇 (上)

1 个月前和另外二位小伙伴一起参加了一个 AI 的比赛。虽然比赛结果不理想,至少我享受到了编程过程中的乐趣。从这次比赛中让我认识到 Go 除了写服务端,写游戏模...

632
来自专栏生信技能树

使用CGP数据库的表达矩阵进行药物反应预测

主页: CGP website 是 Genomics of Drug Sensitivity of Cancer (GDSC)计划的数据

811

扫描关注云+社区