对RNA测序平台Illumina基因分析器的通道影响的统计分析 毕业论文.doc

对rna测序平台illumina基因分析器的通道影响的统计分析摘要下一代测序技术以它低廉的成本及准确的产出结果，在生物信息界及医学制药界获得突飞猛进的发展。rna测序技术作为下一代测序技术在转录组分析的应用技术，广泛应用于转录组学分析研究。illumina genome analyzer是提供rna测序的平台，但测序平台并不是很完美，产生的数据不是绝对精确。本文将分析由于测序机器本身而产生的技术影响，特别是在illumina测序过程中的通道分组带来的数据产生的影响。进而，可以对测序平台进行改进与完善。关键词：下一代测序技术，rna测序技术，illumina测序平台，通道影响 abstractthe next generation of sequencing technology is emerged and highly developed as its low cost and the accurate output results, in the biological information industry and medical pharmaceutical industry. rna sequencing technology as the application of next generation of sequencing technology in transcriptomic analysis, has widely been used in transcriptome research. illumina genome analyzer is a machine which serves as a platform to produce rna sequencing data, but each sequencing platform is not perfect at all, thus the resulting data is not absolutely accurate. this paper will analyze the data variation, that exists because of the effect from the machine itself during the sequencing progress, especially between lanes when sequencing proceed by illumina. for what has been done above, then we can develop and perfect the sequencing platform.keywords: next-generation sequence, rna sequence, illunima genome analyzer, lane effect目 录摘要2abstract3引言51简介61.1 生物背景61.2下一代测序（next-generation sequence）61.3 rna测序（rna sequence）61.4 illumina genome analyzer72统计模型92.1列联表独立性检验92.2超几何分布模型102.3 poisson分布模型112.4 poisson广义线性模型122.4.1线性模型定义122.4.2最小二乘法132.4.3对数线性模型的应用133总结与展望15参考文献16后记与感谢17 引言近年来，下一代测序技术（next-generation sequencing，ngs），又叫深度测序技术(deep sequencing)或高通量测序技术(high-throughput sequencing)，以它低廉的成本及准确的产出结果，在生物信息界及医学制药界获得突飞猛进的发展。而它产生的海量数据使该技术不得不面临信息学难题，有效地处理和分析测序数据变得举足轻重。作为最新发展起来的有巨大市场空间的新兴技术，高通量测序技术无疑面临着机遇和挑战。当然，作为一个新兴技术，高通量测序技术的缺点尚未被人们所了解。产生的数据可供统计学家分析关于技术性的和生物性的可再生能力，错误率和测量偏度等。在这里我们将讨论某些方面的分析和数据处理是必须作为研究新一代测序技术的表达数据的考虑因素。rna测序（rna-sequencing ，rna-seq）作为下一代测序技术在转录组分析的应用技术，广泛应用于转录组学分析研究。它曾被喻为“转录组学中革命性的工具”。相比其它测序技术，如微阵列技术（microarray），rna-seq技术无需预先针对已知序列设计探针，就能在单个实验中以较短的时间同步测量成千上万的基因表达，同时避免产生背景噪声及交互影响。现今很多公司在测序平台的开发和完善上都有了很大进展，目前，以illumina/solexa genome analyzer, helicos heliscope及solid等平台为主。虽然现在的测序平台已经较为完善，但仍不是很完美。因此，产生的数据有较高的精确度及可靠性外，还是由各因素导致的一定的偏差。本文中，我们感兴趣的是在测序过程中对数据的随机分组是否会对测序结果产生影响。本文将以illumina测序平台的原理为例。探讨测序过程中对通道进行分类时产生的数据变量间的关系，即基因的表达数据是否是独立不受分组影响。当然，通道的分组不是唯一的。测序过程中还有流动槽之间的分组，和分析通道分组影响一样，我们进一步也可以对流动槽分组带来的影响进行假设检验。但在本文中不就此展开。其中，通道分组影响和流动槽分组影响均是由于机器本身原因产生的变量，它们有一个专有名词，技术变量（technical variation）。我们知道，测序过程中在技术因素带来的对数据结果的影响之外，还有分段、转录等过程产生的生物变量（biological variation）。科学家在以往的工作研究中已经发现生物变量带来的影响远远大于技术变量。技术变量虽然十分微小甚至在海量数据下可以忽视，但是对技术变量的研究能使我们在下一步的工作中对测序平台进行进一步完善与改进，从而获得更精确的数据。1简介1.1 生物背景遗传信息是由脱氧核糖核酸（dna）分子编码的，它决定了在一个有机体从祖先那继承活细胞的功能。一个dna分子由四种核苷酸（nucleotides , nt）组成，记为a，c，g和t。其中a-t及g-c配对，被称为碱基对（base pairs , bp），是双螺旋结构的dna的基本元素。我们常说的生物术语“基因表达”指的正是分子生物学中的中心法则（the central dogma），在该过程中，一个基因的dna序列作为模板产生核糖核酸（rna）分子，然后在翻译过程中，rna也为细胞合成蛋白质的模板。当一个特定的基因在这一过程出现，我们称基因被表达。随着测量由蛋白质编码的rna即信使rna（mrna）的高通量机器设备的发展，基因表达的研究有了突破性进展。通过测量从单个细胞或一组细胞内的不同的基因的mrna数量，生物学家可以观测到细胞正在生成的过程中的基因信号。通过比较不同的细胞样本中的基因信号（如肿瘤样本和正常组织），我们又可以确认基因表达的差异。事实上，一个基因不是单独的整体，为了保证细胞的复杂性功能，许多基因可以编码多个蛋白质结构。这是通过选择性剪接机制而实现的。在剪接时，基因序列连续的部分可以包含或不包含在最终生成的mrna中。其中包含的部分一般称之为外显子（exons），而不被包含的部分称之为内含子（introns）。不同的外显子子集形成特定的基因，可以连结起来成为一个基因片段（fragment），或称为异构体（isoform）。1.2下一代测序（next-generation sequence）由于dna和rna分子的核苷酸排序不同，决定了遗传信息的不同。因此，近几年，用于测量和确定dna或rna核苷酸顺序的测序技术成为了在分子生物学中最具影响力和建设性意义的工具手段之一。下一代测序技术（next-generation sequencing，ngs），又叫深度测序技术(deep sequencing)或高通量测序技术(high-throughput sequencing)，具有高效的产出和低廉成本等优势，在生物信息界及医学制药界扮演着越来越重要的角色，越来越多的生物学家选择它们作为研究的主要工具，包括转录和结构变化等。最近，几个测序平台也迅速已经发展起来，相比传统的sanger测序技术，它们产生的数据具有高的精确度与可靠性，但其产生的海量数据使该技术不得不面临着数据处理的难题。作为最新发展起来的有巨大市场空间的新兴技术，高通量测序技术无疑面临着机遇和挑战。1.3 rna测序（rna sequence）目前，下一代测序系统除了像传统测序方法应用于基因组的测序或重测序外，还用于研究转录组的rna测序（rna-seq）。rna测序（rna-sequencing ，rna-seq），又叫做转录组测序。利用下一代高通量测序技术对从mrna逆转录生成的cdna进行测序，通过统计相关读段数(reads)计算出不同基因的mrna片段的表达量，我们称之为mrna测序（mrna-seq）。同样原理，各种类型的转录本都可以由下一代测序技术进行深度测序，统称作rna-seq。目前，rna测序广泛应用于酵母，老鼠甚至人类的转录组学分析研究。它曾被喻为“转录组学中革命性的工具”。转录组（transcriptome）即某一特定组织或细胞在复制阶段下转录出来的所有转录基因的集合，因此它是基因的功能和结构研究的基础，能作为分析特定生物学过程以及疾病发生过程中的分子机理的重要工具。rna测序技术能在单核苷酸水平下对任意物种的整体转录活动进行检测，同时，发现新的转录本，精确识别剪切位点。首先，rna分子从特定的细胞中分离，接着被随机打断分成不同长度的片段，然后rna进行转录过程，逆转录成cdna，下一个阶段cdna以特定的长度进行扩张，得到的样本被放入测序仪器进行测序。长期以来被人们熟知的基因表达分析技术是微阵列技术（microarray），而该技术存在的测量转录水平的缺陷，包括由于探针区别上的不足而缺乏绝对的量化，由于饱和度和背景因素而缺乏表达水平的高低分辨率，以及由于需要预先定义的探针序列而缺乏测量转录事件的灵活性。而这些缺点，随着高通量测序技术的发展被轻而易举地避免。rna-seq技术无需预先针对已知序列设计探针，就能在单个实验中以较短的时间同步测量成千上万的基因表达，提供更精确的数字化信号及更高的测序通量，同时避免产生背景噪声及交互影响。rna测序技术已广泛应用于生物学研究、医学研究、临床研究和药物研发等。其流程简单如图1所示。1.4 illumina genome analyzer测序dna和rna分子序列的方法已经发展了30多年。目前，涌现了许多新一代高通量测序平台，如illumina公司的solexa技术和abi公司的solid技术等为标志的新一代测序技术。在此，特以illumina公司的solexa技术为例。illumina公司目前大量使用的测序仪是genome analyzer（ga），其采用边合成边测序（sequencing by synthesis，sbs）的原理。该合成在某一基底上进行，一般地称之为流动槽（flow cell）。在流动槽上，大量的合成反应能够同步地被观察。因此，这些方法的应用领域不仅局限在基因组的重排序，还拓展到了检测基因表达（rna-seq）等领域。下面我们简单对illumina genome analyzer做一下介绍。首先，dna双链分解，适配器与片段末端连接，片段材料通过适配器进行pcr（polymersase chain reaction）扩增。接着执行选择长度的步骤，选择大约200nt（nucleotide）的片段。将需要测序片段放置到流动槽（flow cell）上，作为填装的补充物也同时放置到光滑的表面上。一旦单个分子附加到流动槽，桥式扩增将在原位置上代替pcr扩增，在样本链基础上生成高度相似的多个副本，称为聚类（cluster）。随着补充物和dna聚合酶的添加，测序就开始了。补充物序列是在添加一种溶液和四个碱基上以一次一个单位生成的，每个碱基被标记上不同萤光材料并且包含一个可逆的终端，因此在每次循环中单个模板链只能附加有一个分子。同一聚类的模板副本由聚合酶扩增应该按阶段进行，并且标记上荧光信号的模板序列的聚类可以进行基本的结合，以此可被扫描机器检测到。在第一轮碱基加法和扫描后，可逆的终端从每一个结合的核苷酸中移除，并且这个过程不断重复。在特定的循环次数运行后，基本上保证了所有读段（reads）在最后具有相同的长度。每一次illumina基因组分析器的运行产生了一个流动槽的数据，并且每个流动槽包括八个区域，称之为通道（lanes），其中一组作为对照组。如下图2所示。例如，人的肝脏细胞和肾脏细胞的信使rna从一个个体中分离，并进行分段排序，然后由illumina基因分析器进行测序。目前，每个通道上大约可产生1000万的有效读段。这个数字可能会有较大变化，这取决于输入资料的质量、运行过程中使用的特定的测序化学药剂、或其他未知因素的影响。当然，对dna片段进行头尾两端测序也是可以实现的，即所谓的“双端（paired-end）”测序。2统计模型就如何对rna测序产生的数据进行分析，我们关心的是特定异构体的基因表达读段数，因此我们要对大规模的rna序列建立统计模型的估计。针对illumina平台测序仪器产生的数据，我们不仅要考虑到测序过程中，由于转录过程而产生的一些不可避免的生物性影响因素（biological variation），而且我们需要知道每个测序仪器都不是完美的，在测序仪器运行时会产生因为机器本身造成的误差变量(technical variation)。我们必须将这两个变化分开考虑并进行研究，以对观察结果深入探讨。在以往的研究中，人们已经发现生物性影响因素所带来的差异是显著的，而技术性变量的影响因素往往很小而被人们忽视。下文我们将着重讨论技术性变量在测序过程中带来的影响。以上，本文已经介绍了illumina平台对rna测序的流程。例如，mrna分别从一个人的肝脏细胞和肾脏细胞中分离，随机分段和排序后，用illumina基因分析器进行测序，rna样本分成7个通道，加上一组为对照组。在进行大量实验后，我们发现illumina平台产生的数据相比microarray，具有十分高的可信度和可再生性，但是，我们还发现即使同样的样本，在同样的集合下被放置于不同的通道上进行测试的过程中，总有基于样本误差的系统性的差别。为了方便进一步研究，我们称这样的偏差或差别为通道影响（lane effect）。由于不同的通道分类而产生的通道影响也正是因为机器本身造成的误差变量，即技术变量（technical variation）。为了研究通道影响，我们可以采取多种方法进行检测，当然每个可行性方法都有弊端，但是我们可以通过多方研究对通道影响进行判断。2.1列联表独立性检验首先，我们只做两个通道的考虑，由通道一和通道二测序后产生的数据之间是否存在关系或是互相独立。我们将用fisher精确概率检验来做通道影响的假设检验问题。对测序产生的大量的基因读段数进行独立性检验，我们可以采用列联表来做假设估计。列联表，又称交互分类表，是观测数据按两个或更多属性分类时所列出的频数表，其交互分类的目的是将两个变量分组，然后比较各组的分布情况，以分析判断变量间的关系。因为我们只考虑两个通道的情况，可以简单的做22列联表。如下表格所示。通道1通道2横向加和读段数基因a读段数其他基因读段数纵向加和读段数其中，表示第个通道上基因产生的读段数，. 表示第个通道上所有基因的总读段数，即, . 表示基因在两个通道上的读段数，即 , . 表示所有基因在两个通道上的总读段数，即.这里，为了方便起见，我们记基因a为基因1，其他基因为基因2，统称为基因。我们假定两个通道间不存在通道影响。即关心的假设检验问题为h0：=1 vs h1：1其中 ，表示测序过程中基因读段数的真实比例。接下去，我们要对原假设进行检验。因为我们的目的是利用列联表进行独立性检验，即变量间是相互独立的，因此我们可以利用独立性的一些性质。在此，我们采用广义似然比检验。我们用表示基因在通道上表达的概率。其中，。因为独立性我们可知，而由观测数据，.则，最大似然比易求得，其中，自由度.卡方的值可以查表得到，将与其比较，当以上不等式成立时，即大于时，则表示接受原假设，即变量间互相独立不受影响，也就是说测序过程中的通道的分类对测序结果没有影响。2.2超几何分布模型另一种方法，我们进行的是将完全相同的样本在同一集合下放置于两个通道上进行反复测序。我们假设通道是相互独立的，即产生的读段数来自两个独立样本，由此我们可建立超几何分布的模型。我们假定通道的分组对数据没有影响，即我们建立假设检验h0为不存在通道影响，基于这样的假设检验，我们可以描述为通道一上的基因读段数与两个通道上随机取样的样本读段数相同。首先，我们做如下定义，表示样本在通道上的读段数；表示样本在两个通道和上的总读段数，即；表示在总读段数为的实验下，样本的读段数来自通道的随机变量。基于上述定义，我们可以将假设检验用简单的数学公式表示成 vs ，其中，.原假设成立的条件是等于的概率，即服从超几何分布。由此，我们可以计算每个基因表达的p值。原假设成立的情况下，p值的分布是均匀分布。我们以在观测值下单边的p值为例。.我们发现在同样的样本在两个通道上测序过程中，只有低于0.5%的基因表达的p值足够小而拒绝原假设，即只有极少量的基因因为通道分组的不同而有不同的表达。2.3 poisson分布模型然后，我们又进一步同步地对多个通道进行测序，记通道总数为l。我们要计算的落入每一个基因的外显子区域的读段的数量和使用一个联合poisson模型来研究读段数的联合分布。可以证明，在这个猜想下的对数似然函数参数的联合分布总是凹的，所以最大似然估计（mle），可以很容易地用简单的迭代法求得。在建立统计模型之前，我们先定义下列符号：样本中基因在第条通道上的读段数；该读段数服从poisson分布。：样本在通道上产生读段的总比率；即.：样本在通道上的读段来自基因的比率，其中.则，我们可以得出服从的poisson分布的均值.我们建立不存在通道影响的假设检验，由上述定义，我们可以将原假设表示为在考虑多项分布总体时，即，我们采取拟合优度检验。原假设成立的条件是服从自由度为的分布。其中，表示样本中基因在通道上的读段数；表示读段数来自通道的概率；表示所有读段数的总和。2.4 poisson广义线性模型在实际实验中，不是所有的基因读段数都严格服从poisson分布，我们发现有少部分基因在超出了poisson分布的假设，我们称之为超poisson分布变量（extra-poisson variation）。在上述的poisson分布的基础上，又可以拓展到poisson广义线性模型（generalized linear model, glm），以此计算每个基因表达的p值。我们依然要进行同一个样本在不同通道上基因不同表达的实验。在介绍线性生成模型之前，我们简单介绍一下统计学中线性模型的概念及应用。2.4.1线性模型定义定义1（线性模型）：为一个集合，对所有，线性函数。令的分布函数为，其均值等于且方差为。令为来自定义域的数列。定义为下来自的随机样本，. 则为个观察量，其中且，.注：线性模型的统计意义是的具体形式是未知的，即参数 和的数值是未知的。虽然名称为线性模型，但函数的形式未必是线性的，它还可以是别的形式，如时，或者. 但是有时我们不确定哪种形式的线性模型更优化或建模时产生的误差最小。我们可以计算其相关系数或剩余标准差来判定。相关系数，其越大表示越拟合，即线性模型越好；剩余标准差，其和相关系数是相对的，越小则模型越好。2.4.2最小二乘法在对线性模型进行了定义，我们自然想要知道怎么确定线性模型，即如何求得参数 和的值。目前，我们采取的大多数方法为最小二乘法估计（least-squares estimation, lse）。首先，我们要估计线性模型的误差，我们可以将y与x的关系表示为，其中，.定义2（最小二乘法）令，. 为满足定义1所示线性模型的对观测值. 则 和的值的最小二乘法估计即满足下式的 和的最小值。基于上式，我们只要对 和求偏导，令其等于零，便可求出 和的最小值。我们得到, . 其中，2.4.3对数线性模型的应用目前，我们对线性模型和估计线性模型函数参数的最小二乘法进行了介绍。而对于本文关心的如何对通道影响进行估计。我们在上节讨论的poisson分布的基础上，对表示样本在通道上的读段来自基因的比率建立线性模型。我们知道，则，线性模型. 它是2.4.1介绍的线性模型的变型。在大量的数据为依据上，我们发现上式poisson广义线性模型是拟合最优的。其中表示基因在所有通道上的平均表达水平，表示通道对基因表达水平的标准均值的影响。则定义原假设为不存在通道差异而成为影响因素的假设检验可表示为 vs 在进行假设检验时，我们采用广义似然比检验，其中， 为在全空间下的最大似然估计（mle），为在假设空间下的最大似然估计（mle）。由于poisson分布的性质是变量的方差与均值相等，即，则假设检验成立的条件是lrt服从自由度为1的分布。3总结与展望下一代测序技术是近几年高速发展的新兴技术，在医学、制药及生物领域有十分大的发展空间，但大部分人对这年轻的科技尚不熟悉，本文首先简单对下一代测序技术，rna测序以及广泛运用的测序平台illumina genome analyzer进行介绍。然后我们关心如何分析illumina测序平台产生的数据，实验发现，经过该平台测序得出的数据具有相对高的可靠性，但是还是不可避免地会有一些偏差，我们又对该偏差进行了分析。以illumina测序平台的流程为例，mrna从细胞分离出来后随机分段并反转录成cdna，然后由测序平台测序，测序过程中片段被放置于流动槽，每个流动槽包括8组通道，其中一组为对照组。本文关心的是测序过程中的通道的分组是否会对测序数据的结果有影响，即通道影响（lane effect）。本文做了三种假设检验方法对测序过程中的通道影响进行估计。首先，我们想到采用最基本也是最基础的fisher精确概率检验，直接对两个通道上基因的读段数据进行分析，以列联表的格式直观呈现数据并进行独立性检验。其次，我们同样先从简单的两个通道考虑，将同样的样本以同样的集合放入通道，我们假设通道上基因的读段数服从超几何分布而进行概率分析。然后，我们将两个通道延拓为多个通道，对于多个通道的测序数据，我们又进行了通道上基因的读段数服从泊松分布的假设，同时进行拟合优度检验。最后，我们又考虑到可以用线性生成泊松分布更好地解决问题。通过大量实验得出的数据分析显示，通道的不同对测序基因表达的影响是很小的，即技术变量（technical variation）虽然存在，但在海量的数据测序过程中，几乎是可以忽略的。以此，我们也证明了下一代测序技术的日渐完善，各测序平台进行测序过程中产生的数据在一定程度上来说有相当高的精确度和可靠性。当然，不同的通道分类只是众多技术变量的一种，此外不同的流动槽测序可能产生的差异也属于由于机器本身而产生的测序结果的误差，即技术变量。下一代测序技术的飞速发展及广泛的应用，势必市场对其的期望值会更高，在未来的研究中，我们大可在技术变量上多做文章，以更进一步地提高和完善rna测序平台的准确性。参考文献1auer p l, doerge r