第五章  新一代测序NGS：重测序的回帖和变异鉴定

5.1  新一代测序

• 从二十世纪前，人类认识到DNA的重要性后，一直以来将测序----确定一个特定DNA分子的序列----作为理解生命的重要方法。
• 真正可以大规模运用的核酸测序方法，是1977年由英国生物化学家Frederick Sanger提出并实现.
• Sanger测序法的广泛应用使得大规模测定基因组序列成为可能，并为人类最终在20世纪完成人类基因组草图的绘制奠定了坚实的基础。
• 2005年454技术发布为标志，新一代测序技术开始登上舞台。
• 与经典的Sanger Sequencing相比，新一代测序技术可以产生更多的reads，从而得到更高的测序深度，因此也常常被称为深度测序。（Next Generation Sequencing/ Deep Sequencing）
• 另一方面，相比于Sanger Sequencing，新一代测序产生的平均reads，长度更短，却错误率更高，从而给后续的生物信息学分析提出了更大的挑战。
• 除ABI公司的SOLiD之外，深度测序仪产生的数据通常被保存为FASTO格式；其中对于每条read，在具体的核酸序列信息之外，还包括了每个碱基对应的质量信息。
• Quality: Given p = the probability of a base calling is wrong, its Quality Score can be written as 
• Q值进一步按照转换标，被编码为ASCII码字符，保存在Fastq文件中。
• 经验规则：在实际分析工作中，常常会讲质量分数小于20（即错误概率大于0.01的碱基）认为是不可靠的；如果这样的碱基超过read长度的20%，即会考虑将此条read丢弃。
• 为了克服reads过短带来的问题，新一代测序广泛运用了pared-end reads，也即同时对一个较长的片段两端进行测序。这时，相应的reads名字会在最后分别加上/1和/2以示区别。
• 除了基因组DNA测序外，新一代测序技术还被用来研究表观遗传学修饰、RNA转录组以及蛋白-DNA相互作用等重要的生物学问题。
• 通过将个体基因组（individual genome）测序产生得到的DNA reads map到参考基因组，可以有效地发现不同个体（individuals）之间存在的遗传差异。在此基础上，通过与特定的表型差异相联系，可以开展关联分析（Association Study），研究表型差异的遗传学基础。从而，为研究遗传疾病的机制，进而探索后续的诊断和治疗方案提供重要的线索。
• RNA-seq是利用深度测序来研究转录组的技术，通过系统测序细胞中表达的全部转录本，RNA-seq技术使得研究人员可以快速确定转录组，进而鉴定存在的可变剪切体(Alternative splicing isoform). 与此同时，通过统计每个基因(locus)对应的reads数目（number of mapped reads），可以估计基因的表达水平，进而进行差异表达分析、聚类分析等统计计算以确定与给定生物过程相关的基因。
• 与RNA-seq类似，ChIP-Seq也是利用深度测序相关技术来研究转录调控的技术，ChIP-Seq通过深度测序技术，测定与特定抗体结合的DNA序列，进而腿短Protein-DNA相互作用。
• 通过选择不同的抗体，ChIP-Seq技术既可以用来检测转录因子的结合点，也可以用于检测特定染色质修饰区域。因此该技术在转录调控乃至表观遗传学研究上均获得了广泛应用。

5.2  序列回帖和变异鉴定

• reads mapping是指将测序得到的DNA片段----即reads----定位在基因组上。
• 通过reads mapping，在克服了深度测序产生的Reads过短导致的技术困难的同时，也方便将测序得到的数据与前期研究产生的注释结果进行整合。
• reads mapping往往被作为深度测序数据分析的第一步，其质量的好坏，以及速度的快慢，都会直接对后续的分析工作产生影响。
• reads mapping可以看成双序列对比，但是与之前提到的不同，首先是长度上，两者的长度有着跨数量级的差异；reads的长度通常不超过100bp，而参考基因组却通常在上百Mb；其次是数据量，深度测序产生的reads动辄达到几百Gb，相当于几十个人类基因组；在双序列对比中，我通常假定输入序列本身不会有错，但目前深度测序产生的reads质量却参差不齐，错误率较高（1e-2~1e-5）
• reads mapping需要一个混合了全局对比和局部对比的混合型的alignment（hybrid alignment）。
• “Global” for short sequence(i.e. NGS read)
• But "Local" for long sequence(i.e. Reference Genome)
• In particular, the surrounding "overhang" gaps should be not penalized.
• 借鉴BLAST的思路，采用Seeding-and-extending的策略。具体来说，我们首先通过对基因组建立索引，从而将Read快速定位，而后再通过标准的动态规划算法，构建最终的alignment。
• 所谓索引本质上就是对数据的分组，基于对原始数据中的key应用索引函数处理。
• 好的索引函数应当尽力确保划分后的数据组大小接近，这样可以更好地降低平均的检索时间。
• Hash是一种常用的索引方法。
• 通过对原始数据中的key应用哈希函数计算出的值作为索引地址，在理想情况下，只需常数时间，即O(1)，即可完成数据的查找。
• reads mapping过程中，为了提高灵敏度，通常会允许若干碱基的错配。
• 抽屉原理：桌上有十个苹果，要把这十个苹果放到九个抽屉里，无论怎样放，我们会发现至少会有一个抽屉里面放不少于两个苹果。这一现象就是我们所说的“抽屉原理”。
• 例如，依据抽屉原理，将read划分成3个不重叠的块，则可确保在错配不超过两个的情况下至少有一个片段可以作为与基因组完全匹配的种子来查询索引，从而大大提高reads mapping的速度。
• 如果搜索的时候，允许较多的错配位点，则需要把序列分解成很多小的区块，那就回导致性能的急剧下降。
• 2009年开始，数据压缩算法中常用的前缀树（prefix tree）和后缀树（suffix tree）开始被应用于reads mapping。这些数据结构通过合并共享子串（shared substring）降低内存消耗，并且可以方便的实现对最长公共子串的查找。
• 目前流行的BOWTIE、BWA、SOAP3等工具均采用了基于后缀的Burrows-Wheeler transform(BWT).
• 与BLAST不同的是，对于一个错配位点，我们需要考虑是测序错误引起假象（artifact）的可能性。
• 2008年，Li Heng与Richard Durbin提出了mapping quality分数的概念，来处理假象的问题。
• 在实际数据分析中，更多的用Mapping Quality而非序列对比分数来筛选真正reads mapping的数目。
• 根据遗传变异的尺寸，可以将之分为在单个碱基水平的单核苷酸变异SNV（Single Nucleotide Variant）和涉及到多个碱基的结构变异SV（Structural Variation）.
• SNV可以划分为碱基替换SNP（最常涉及到的遗传变异）和indel
• SV可以划分为大尺度的插入/删除（Large-scale insertion/deletion）;倒转，Inversion; 异位，Translocation;拷贝数变异， Copy Number Variation(CNV)等多种类型。
• SNP calling只是确定哪个基因组位点存在变异，并不涉及对应位点的基因型。
• Genotyping则是在SNP calling的基础上，进一步确定变异位点的基因型，包括是纯合还是杂合等等。
• 实践中，基于概率的SNP calling和Genotyping模型更为常用。

5.3  序列回帖和变异鉴定的分析演示

• 以人类的线粒体基因组为例。
• BWA分析的第一步是首先将线粒体基因组进行索引。
• 在索引之后，我们将分别对高通量测序得到的双端的reads进行mapping。
• 在mapping之后得到的是.sai文件，这种文件暂时不能打开，需要进行最后一步BWA的处理。
• 比对生成的是BAM格式的文件。
• .bam文件是二进制的文件，需要使用samtools view来进行观看。
• 在bam格式中，每一行都表示了一个reads的比对结果，包括reads的名字，它每一个比对到线粒体染色体的名字，以及坐标，以及其他的一些信息，但是这样的信息时很难直观地进行寻找突变位点的分析的，所以接下来还会进行别的处理。
• 处理之前，首先对这些BAM文件进行排序。
• 排序之后，将对它进行索引。
• 然后使用GATK提供的indel realignment工具对于原始的比对结果进行微调（原始的比对可能在某些indel附近存在错误，所以需要微调）。
• GATK的工作分两步完成。
• 第一步：使用RealignerTargetCreator工具，寻找一些需要realignment的位点。
• 第二步：使用IndelRealigner来对BAM文件进行处理。
• 上述两步工作结束之后，还将使用GATK提供的Base[Quality]Recalibration来对测序得到的base的quality来进行调整，这一步同样是分为两部分。
• 第一步：提供一个用于训练的已知的变异位点的数据集合；
• 第二步：执行TableRecalibration工具的命令，最终来生成我们所需要的final.bam格式的文件。
• 现在可以进行variance calling的工作。有多种工具。此处演示两种。
• 第一种是samtools提供的mpileup；第二种是GATK提供的UnifiedGenotyper.
• samtools的操作，我们首先需要对fasta文件建议另外一个samtools所需要的索引。
• 然后运行samtools mpileup这个命令，在这个命令中，我们提供了原始的基因组的序列以及我们经过处理后得到比对后的结果BAM格式的文件。然后产生了.VCF的文件。
• GATK提供的UnifiedGenotyper操作类似。而得到的VCF文件与samtools的基本相同，只是在header上有所区别。找到的variant的位置，是有所区别的。具体差别则需要具体的实践操作，在实践中寻找究竟哪一个才是更适合你的。