突变需要学习的知识点

学习Jimmy直播我的基因组系列

SNV（经常会被混着SNP来称呼，我们先不纠结这个细节）

1. 通常一个人的全基因组测序数据可以挖掘到四百万个SNVs(跟参考基因组不一样的单碱基位点)，还有五十万的indels(insertions or deletions),但是得到的数据通常是以vcf文件格式给出的(自行搜索什么是vcf格式)，比如下面：

2.首先记住一个很重要的知识点，变异是相对的！变异不等于突变。 简单说一下什么是找变异，变异跟突变有什么区别呢？举个栗子：有国际组织规定了人类的参考基因组（如UCSC,ENSEMBL,NCBI等，前面帖子都有讲)，就是 AAAAA(这里简化一下，就5个碱基，其实人类基因组多达30亿个) 。现在通过给自己测序得知，我与之对应的是AGCAA，那么我相比国际基因组来说，就是2个变异位点，位于基因组的坐标2和3，但是它们还不能说就是突变。

如第二位碱基，虽然我的是G，参考基因组是A，但是全球已经测序了几百万人，而我查看了他们的测序结果，其中99万人都是G，这说明是参考基因组出现了问题，可能是国际组织当年恰好选择了一个人是A，所以就规定第二个碱基是A。所以虽然我用软件找到了我的这个位点相对于参考基因组是来说，是一个变异，但是这恰好是好事，完全不用担心，我们也不需要用突变这个单词来描述它！

那么接下来看第3位碱基，同样，国际组织规定了是A，而我却测了个C，但是全球已经公布的一百万人里面99.999万人都跟参考一样，就是A。有一个人和参考基因组对应的碱基不一样，不一样的那个人是个有病的患者，这个时候，你就惨了，这个变异，就是突变了！

很多变异其实只是造成人种多样性的原因，是构成人独特性的基础，而那些跟疾病相关的变异，我们通常就会叫做是突变！

因我只举了2个极端的例子，所以大家可能会误以为，跟大多数人一样，就没事了！其实也并不是这样，一般来说，在正常人的数据库里面出现了5%的变异就可以认为没什么大的危害，而且变异还可以分成germline、somatic、de novo等情况，如果是特定性的针对某种疾病还可以找driver的mutation.

vcf就是一个尤为重要的基本数据格式，里面描述的是变异位点的信息，是一个表格形式，不同的表头下面有着不同的信息。VCF格式本来由千人基因组计划提出来，方便描述他们找到的海量(当时是海量)变异位点。本质上也是个文本文件而已，普通编辑器打开即可。但是它对每一行每一列有具体的定义，包括文件最前面一些#开头的注释信息(这个非常重要，后面每一个位点的描述的tag都在这个注释信息里面可以找到)

4.vcf文件的正文部分

vcf的正文部分，必须要有的是前面8列，一般来说可以有10列，分别是：

-1. #CHROM

2. POS

3. ID

4. REF（参考序列的碱基）

5. ALT（Variant的碱基）

6. QUAL

7. FILTER [来自于##FILTER]

8. INFO

9. FORMAT

10. 可能会有样本的名称本

CHROM 和 POS：参考序列名和variant的位置；如果是INDEL的话，位置是INDEL的第一个碱基位置。

ID：variant的ID。比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’."表示其为一个novel variant。

REF 和 ALT：参考序列的碱基 和 Variant的碱基。

QUAL：Phred格式(Phred_scaled)的质量值，表 示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。

FILTER：使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为”PASS”;若variant不可靠，则该项不为”PASS”或”.”。

INFO：这一行是variant的详细信息，内容很多，以下再具体详述。

FORMAT 和 TTG11B：这两行合起来提供了’TTG11B′这个sample的基因型的信息。’TTG11B′代表这该名称的样品，是由BAM文件中的@RG下的 SM 标签决定的。

image

前面7列都很简单，顾名思义，分别就是该变异位点位于参考基因组的哪条染色体，哪个位置，是否被一下数据库给标记了ID(通常说的是dbSNP)，该位置的参考基因组是什么碱基，这个变异位点变异成了什么碱基。找到这个变异的软件给它的质量值是多少，是否合格。下面这个表格里面我们可以看到第十列就是'realign'，可以看到比对时候@RG留下来的sam的样本名称，就可以知道这个vcf是经过realign的那个bam里面call出来的突变。

vcf只学七列是远远不够的，我们有必要下功夫把较为复杂的第8列和第9列的内容好好学习一下！

image

第8列 INFO 就非常复杂了，该列信息最多了，看起来是一列，但是里面可以无限包容，可以根据字段拆分成多列，都是以 “TAG=Value”,并使用”;”分隔的形式。其中很多的TAG含义在VCF文件的头部注释信息##INFO中已给出。

通常我们熟悉的tag有：

AC，AF 和，AN[A开头的多和等位基因有关]：

AC(Allele Count) 表示该Allele的数目；

AF(Allele Frequency) 表示Allele的频率； 

AN(Allele Number) 表示Allele的总数目。

对于1个diploid sample[二倍体样本]而言

则基因型 0/1 表示sample为杂合子，Allele数为1(双倍体的sample在该位点只有1个等位基因发生了突变)，Allele的频率为0.5(双倍体的 sample在该位点只有50%的等位基因发生了突变)，总的Allele为2； 基因型 1/1 则表示sample为纯合的，Allele数为2，Allele的频率为1，总的Allele为2。

DP：reads覆盖度。是一些reads被过滤掉后的覆盖度。[注意，第八列和第九列都有DP，都表示该位点覆盖深度的信息，但是详细意义可能是不同的大家可以探究一下，在head里面就可以找到相应信息]

Dels：Fraction of Reads Containing Spanning Deletions。进行SNP和INDEL calling的结果中，有该TAG并且值为0表示该位点为SNV，没有则为INDEL。[这个值很重要，可以根据这个tag分离indel和snv]

image

如果你觉得call变异的软件默认给出的tag不符合你的要求，你可以继续用其它软件在该列里面不停的增加tag，我见过给该列直接添加到180个tag的，我们后面主要讲如何来添加tag。

有了这8列，已经是标准的vcf文件了，但是大家肯定会奇怪，还没有关于这个位点的基因型，测序深度的描述的信息。

这就是属于后面的第9列FORMAT规定的了，如果有多个样本，就会按照第九列的格式不停的增加下去。

第九列可以是GT,DP,FT,GL,PL,GP等等，都可以在该vcf文件的表头里面找到关于它们的解释。前面所讲的 ##FORMAT 表头部分 便是对第九列的解释

第九列相对于第八列来说没有那么复杂的信息，数据格式是比较固定的，其中包含的信息也很重要，主要是某一个特定位点基因型，测序深度的描述，因此有必要弄清楚。

image image

第9列数据，包含两列内容，两列内容是对应的，前者为格式，后者为格式对应的数据。

GT：样品的基因型（genotype）。两个数字中间用’/"分 开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。因此： 0/0 表示sample中该位点为纯合的，和ref一致； 0/1 表示sample中该位点为杂合的，有ref和variant两个基因型； 1/1 表示sample中该位点为纯合的，和variant一致。

AD 和 DP：AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid中则是用逗号分割的两个值，前者对应ref基因型，后者对应variant基因型； DP（Depth）为sample中该位点的覆盖度。

GQ：基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值，表示在该位点该基因型存在的可能性；该值越高，则Genotype的可能性越 大；计算方法：Phred值 = -10 * log (1-p) p为基因型存在的概率。

PL：指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。和之前不一致，该值越大，表明为该种基因型的可能 性越小。 Phred值 = -10 * log (p) p为基因型存在的概率。

最需要理解的就是DP4和GT了：

image

第十列的话就是样本的信息 可以在比对的时候使用@RG来做一个标记