Variant 分析阶段小结1-基础碎碎念
6000字，约12分钟，思考问题的熊 专栏9
遗传变异碎碎念
什么是遗传变异
所谓遗传变异是生物体内遗传物质发生变化而造成的可以遗传给后代的变异，这些变异导致了生物在不同水品上体现出遗传的多样性。生物信息学中各种基因组研究的基础就是遗传变异的研究，比如进化和各种表型的研究。
遗传变异包括单核苷酸多态性(SNP)，小片段的插入缺失(Indel)，结构变异(SV)，拷贝数变异(CNV)等等。区分这些变异简单的方法就是变了几个，其中SNP是单个核苷酸的改变，indel通常是50bp以下的变异，SV和CNV则要更长。
遗传变异如何检验
SNP 检测方法主要就是基于高通量数据中的reads在某个位点上的碱基同时结合概率统计进行检验。
基因组结构变异检测主要有四种思路：
- 双端 reads 连配 (pair-end mehtod) 首先进行mapping，然后评价双端reads的方向信息和距离是否与建库结果一致
- reads 覆盖深度 (read-depth method) 多用于检测拷贝数变异，假设比对到基因组的覆盖深度符合泊松分布，然后检测与该分布不一致的区域。
- reads 分割法 (split-read) 首先保证双端的一个read 可以完整的唯一比对，然后把另一条read进行拆分再去进行比对。
- 从头组装 (assembly) 啥也不管直接拼，拼完了再去和参考基因组比对。
相关技术及应用
目前比较常用的基组重测序技术包括全基因组重测序以及简化基因组测序(GBS，RAD-seq)。这些技术可以用来进行基因组遗传多样性的研究（不同品系中是否有不同的基因存在），群体遗传计划研究，构建遗传图谱，快速检测突变位点等等。
基本概念碎碎念
SNV 和 SNP
SNP 和 SNV 都是单碱基的突变，但是SNP 是多了一个频率属性的SNV，比如在群体中1%以上。
biallelic and multiallelic
biallelic 表示在基因组的某个位点上有两个等位基因，即可以有一个突变等位基因。换句话说这个位置上可能存在一个和参考基因组相同的碱基和一个和参考基因组不同的碱基或者是一个deletion。 multiallelic 多等位基因表示在基因组的某个位点可以观测到三个或者多个等位基因，在vcf文件中可以看到两个或者三个非参考基因组的突变。多等位基因并不常见，在各种vcf文件相关工具中，都可以统计这两种信息。
Transition vs Transversion
关于转换和颠换用下面三幅图就可以非常清楚的展示。转换(transition)则是嘌呤被嘌呤，或嘧啶被嘧啶替代，颠换(transversion)是指嘌呤与嘧啶的变化。转换突变比颠换突变更常见，且与颠换相比在氨基酸序列上产生差异的可能性更低。一般Transition 是 Transversion 数量的2倍。
What is the difference between Transition and Transversion?
参考 https://www.differencebetween.com/difference-between-transition-and-vs-transversion/
SNP 种类
全基因组SNP突变可以分成6类（C>A, C>G, C>T, A>C, A>G, A>T）。以A:T>C:G为例，此种类型SNP突变包括A>C和T>G。由于测序数据既可比对到参考基因组的正链，也可比对到参考基因组的负链，当T>C类型突变出现在参考基因组正链上，A>G类型突变即在参考基因组负链的相同位置，所以通常也会将T>C和A>G划分成一类。
SNP的可能影响
如果SNP发生在编码区，根据密码子简并性 SNP 不一定会引起编码氨基酸的改变，不引起任何变化的叫做Synonymous SNP， 而引起氨基酸变化的叫做 Non-Synonymous SNP。如果编码的某种氨基酸的密码子变成另一种，会导致多肽链的氨基酸种类和序列发生改变，这就是一个错义突变。当突变使一个编码氨基酸的密码子变成终止子时，则蛋白质合成进行到该突变位点时会提前终止，这时就是无义突变。
VCF 文件碎碎念
如果想要完完全全彻彻底底十分清楚地理解vcf文件，最好地选择是看这个28页PDF的官方说明文档。
vcf 文件由header和数据结果两部分组成，如果想初步且全面的理解vcf文件各种结果意义,认真地读一下header部分就可以。header本身是带着 ##开头的注释信息， 一方面它解释了数据结果部分含义，另一方面记录了这个vcf文件的来历，例如经过了哪些软件的处理以及所用到的参数。有时候一些公司分析的结果文件有的会删掉部分header有的不会，通过header部分，就可以看看他们是如何处理数据的。
hader
下面是两个不同软件samtools和GATK跑出来的header
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##samtoolsVersion=1.3.1+htslib-1.3.1
##samtoolsCommand=samtools mpileup …… 这里是samtools命令……下面是公司跑出来的header
##fileformat=VCFv4.1
##FILTER=<ID=LowQual,Description="Low quality">
##FILTER=<ID=SnpCluster,Description="SNPs found in clusters">
##FILTER=<ID=my_snp_filter,Description="QUAL<30||QD < 2.0 || FS > 60.0 || MQ < 40.0">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=exported,Date="Mon Dec 19 18:35:07 CST 2016",Epoch=1482143707597,Comma上面的信息显示，公司给的vcf数据版本是4.1，自己跑出来的结果是4.2。FILTER显示公司的数据已经经过了若干条件的筛选和过滤。FORMAT和后面的INFO都是结果中标记信息的解释。另外，还会有参考基因组信息。
结果信息如下
GATK结果
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 20
Chr1 23314 . G GTATC 654.77 . AC=1;AF=0.5;AN=2;BaseQRankSum=-0.399;ClippingRankSum=0;DP=43;ExcessHet=3.0103;FS=
0;MLEAC=1;MLEAF=0.5;MQ=58.99;MQRankSum=2.921;QD=34.24;ReadPosRankSum=-0.897;SOR=0.818 GT:AD:DP:GQ:PL 0/1:2,16:18:14:671,0,14samtools 结果
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT /mnt/zf/ory_reseq_HXM/rawdata/02samtools_treat/Rice_D558-R02
-I_good_R1_rmdup.bam
Chr1 31071 . A G 228 . DP=47;VDB=0.888219;SGB=-0.693146;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,21,23;MQ=60 G
T:PL 1/1:255,132,0
Chr1 31478 . C T 228 . DP=42;VDB=0.689358;SGB=-0.693146;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,21,21;MQ=60 G
T:PL 1/1:255,126,0固定信息
INFO信息
以 “TAG=Value”, 不同的tag之间使用；分隔。samtools和GATK跑出来的结果顺序并不一样，各自的tag数量也不太一样。需要小心。
仔细观察上面AC值不同的突变位点，可以体现出下面的信息：
对于二倍体样本：基因型GT 0/1 表示样本为杂合子，Allele（AC）为1（二倍体样本在该位点只有1个等位基因发生突变），Allele频率（AF）为0.5（二倍体的样本在该位点只有50%等位基因发生突变)，总Allele（AN）为2；基因型 1/1 表示样本为纯合，Allele（AC）为2，Allele频率（AF）为1，总Allele（AN）为2
不同的软件跑的tag不同，但是在自身的header里面都会有详细的解释；这一列的tag可以无限添加，比如加上各种后期的注释信息。
samtools的结果会产生一个DP，其解释是Raw read depth；
DP4的含义是Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases；
MQ0F的解释是Fraction of MQ0 reads (smaller is better)，
MQSB的意思是Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)；
RPB：Mann-Whitney U test of Read Position Bias (bigger is better)；
MQB：Mann-Whitney U test of Mapping Quality Bias (bigger is better)；
VDB checks if variant bases occur at random positions in the alignedportion of the reads.
GATK可能会产生BaseQRankSum，含义是Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities；QD含义是Variant Confidence/Quality by Depth；MLEAF指的是Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed.
FORMAT信息
对于二倍体，GT值表示这个样本在该位点所携带的两个等位基因。0表示跟REF一样；1表示表示跟ALT一样；2表示第二个ALT（如果有）。当只有一个ALT 等位基因的时候，0/0表示纯和且跟REF一致；0/1表示杂合，两个allele一个是ALT一个是REF；1/1表示纯和且都为ALT。
关于AD和DP，These two values will usually, but not always sum to the DP value. Reads that are not used for calling are not counted in the DP measure, but are included in AD.
关于两个不同的DP，The INFO DP is the combined depth over all samples. The FORMAT DP is for the sample. DP at sample level is for individual samples. At site level in INFO, it should be approximately sum of all sample level DPs
关于GQ和PL, GQ is the difference between the PL of the second most likely genotype, and the PL of the most likely genotype. the values of the PLs are normalized so that the most likely PL is always 0, so the GQ ends up being equal to the second smallest PL, unless that PL is greater than 99. Basically the GQ gives you the difference between the likelihoods of the two most likely genotypes. If it is low, you can tell there is not much confidence in the genotype.
如果理解了上面这段话的意思，就可以理解 GT:AD:DP:GQ:PL0/1:2,16:18:14:671,0,14的结果并不怎么好。