4.2 序列回帖和变异鉴定 NGS：Reads Mapping

Reads Mapping 是啥？

是指测序得到的DNA片段也就是Reads，定位到基因组上，通过Reads Mapping 在克服了深度测序产生的Reads过短导致的技术困难的同时，也方便利用基因组位置作为桥梁来将测序得到的数据与前期研究产生的注释结果进行有机的整合
作用

往往被作为深度测序数据分析的第一步，其质量的好坏以及速度的快慢都会对后续的分析产生影响
Mapping: Input Data

本质上还是双序列比对问题
- Reference Genome
  
  — Nucleotide
  
  — Length: Hundreds of Mb per chromosome
  
  — ～3Gb in total(for human genome)
- Reads
  
  — Nucleotide, with various qualities(relatively high error rate 1e-2 ~ 1e-5)
  
  — Length: 36~80 bp per read
  
  — Hundreds of Gbs per run
BLAST Ideas: Seeding-and-extending

通过对基因组进行索引，从而将Reads快速定位，而后再通过标准动态规划算法来构建最终的alignment
1. Find matches(seed) between the query and subject 寻找查询和主题之间的匹配
2. Extend seed into High Scoring Segment Pairs(HSPs) 将种子扩展为高分段对(HSP)
  
  — Run Smith-Waterman algorithm on the specified region only 仅在指定区域内运行Smith-Waterman算法
3. Assess the reliability of the alignment 评估对齐的可靠性
- 数据索引
  
  所谓的索引本质上就是对数据的分组，基于对原始数据中的keys应用索引函数来进行处理，可以将原始数据划分为若干更小的组，并通过将与keys相关的检索限制在某个小组里，来缩小搜索空间，从而降低搜索时间
  - 哈希
    
    例：对一段基因组序列，首先参考BLAST思路，将之切分成若干的细的片段，而后针对每一个细的片段计算其哈希值作为其在索引表中的地址，并在地址表中保存这段序列及其在基因组中的坐标，在Reads中再次出现这个片段的时候就可以在地址表中迅速找到其在基因组中出现的位置
  - 抽屉原理提高Reads Mapping速度
    
    在Reads Mapping过程中，为了提高灵敏度，通常我们会允许若干碱基的错配，这时我们就可以依据抽屉原理，将Reads划分为若干不重叠的块。
    
    例如我们将Reads划分为3个不重叠的块，那么就可以确保在错配不超过两个的情况下，至少有一个片段可以作为与基因组完全匹配的种子来查询索引，从而大大提高了Reads Mapping的速度
  - 前缀树、后缀树 —— 大大提高内存利用的效率和比对的速度
    
    早期的Reads Mapping算法ELAND、MAQ、SOAP1等等均广泛的应用了这个思路，但是如果搜索的时候允许较多的错配位点，就需要把序列分成很多小的区块，那这会导致性能急剧下降。因此，09年开始用于数据压缩算法的前缀树和后缀树开始被应用于Reads Mapping，这些数据结构通过合并共享子串降低内存消耗并且可以方便实现最长公共子串的查找。
    
    例BOWTIE、BWA、SOAP3等工具均采用了基于后缀的BWT转换，可以逐位对比对片段基因延伸，大大提高内存利用的效率和比对的速度
- 找位置
  
  通过索引找到相应的地址之后，我们就可以参照之前BLAST思路，向左右两个方向延伸扩展之后用动态规划算法来确定最后的比对
  
  与BLAST搜索不同的是，新一代测序技术的错误率很高，对于一个测序位点我们就需要考虑它是由测序错误引起的假象的可能性，因此08年提出了Mapping Quality 分数的概念来处理这一问题
  
  首先假定看到的错配都是由于测序错误引起的，这个时候的概率就是各个错配碱基所对应的概率乘积，将其做过log transformation后的数值称之为SQ。
  
  在所有的Mapping位置都计算了一个相应的SQ值之后，就可以进一步的对每个位置计算出错误mapping的概率E，这里计算E的时候是同时考虑序列的相似程度和测序质量，因此在实际的测序中，会更多的使用Mapping Quality而不是序列比对分数来筛选真正的Reads Mapping的位置
Genetic Variants —— 遗传变异

当将测序得到的Reads正确Map到基因组之后，就可以开始来鉴定遗传变异，根据遗传变异的尺寸，可以将之分为在单个碱基水平的单核苷酸变异即SNV，和涉及到多个碱基的结构变异SV
- SNV：Single Nucleotide Variant 单核苷酸变体
  
  — Substitution(SNP) 碱基替换
  
  — Indel: insertion/deletion
- Structural Variation(SV) 结构变异
  
  — Large-scale insertion/deletion 大规模插入/删除
  
  — Inversion 倒置
  
  — Translocation 易位
  
  — Copy Number Variation(CNV) 拷贝数变异
- SNP Calling is NOT Genotyping
  - "SNP calling aims to determine in which positions there are polymorphisms or in which positions at least one of the bases differs from a reference sequence”
    
    "SNP调用的目的是确定在哪些位置存在多态性，或在哪些位置至少有一个碱基与参考序列不同”
  - "Genotype calling is the process of determining the genotype for each individual and is typically only done for positions in which a SNP or a ‘variant’ has already been called.”
    
    "基因型调用是确定每个个体的基因型的过程，通常只针对已经调用了SNP或’变体’的位置。”
- Counting: an intuitive (and naïve) approach
  
  在Reads被正确Map的前提下，我们可以参考特定位点上的Map的Reads，运用经验规则来基于计数的方法同时完成SNP calling和Genotyping
  
  例：
  - Counting high-confident, non - reference allele (i.e. Quality >= 20)
    
    — Freq(突变位点的频率) <20% or > 80%: homozygous genotype 纯合子变异
    
    — Otherwise: heterozygous 杂合变异位点
  - Works well for “deeply sequenced regions” (DSR), i.e. depth > 25x
    
    对 “深度测序区域”（DSR）效果良好，即深度大于25倍
    
    — But suffer from under-calling of heterozygous genotypes for low -coverage regions
    
    但对低覆盖率区域的杂合基因型则存在低估的问题
    
    — And can’t give an objective measurement for reliability
    
    不能对可靠性进行客观的测量
- A Simple Probabilistic Model for Genotyping
  1. For a diploid genome, there will be at most two different alleles (A and a)observed at a given site:
    
    对于一个二倍体基因组，在一个给定的部位最多可以观察到两种不同的等位基因(A和a)
    
    — 3 possible genotypes: <A,A>, <A,a>, <a,a> — Number of A： k; Number of a: n-k
  2. Then, the probability for each genotypes is
  3. Bayes Formula can be further employed to calculate posterior probabilities, i.e. P(<A,A>| D), P(<a, a>|D), and P(<A,a>| D) if we can estimate the prior probabilities P(<A,A>), P(<a,a>) and P(<A,a>)