专利名称：基于rmi指数的高通量测序数据分析方法
技术领域：
本发明属于生物技术领域，涉及新一代高通量测序技术数据分析。
背景技术：
本发明是一种适用于高通量测序技术—Solexa [Bennett, S.，SolexaLtd. Pharmacogenomics, 2004. 5 (4) :p. 433-8]的新分析方法。本方法同时适用于转录组 测 序一RNA-seq[ffold, B. and R.M.Myers, Sequencecensus methods for functional genomics. Nat Methods, 2008. 5 (1) :p. 19-21])和免疫共沉淀测序一ChIP-seq[Johnson， D. S. , et al. , Genome-wide mapping of in vivo protein—DNA interactions. Science, 2007.316(5830) :p. 1497-502]两种实验方法产生的数据。高通量测序技术是近两年来生物技术领域的重要突破，新一代的测序技术将传统 的Sanger测序效率提高了数百倍，同时价格也大大下降。高通量测序技术的出现使得许多 的极有前景的生物医药应用成为可能1，癌症基因组。2，个性化医疗与诊断。3，药物靶标 筛选。高通量测序技术能否在这些领域的取得进展，其关键在于分析方法及软件的创新。本 专利提出了一种新的分析方法，可广泛用于高通量测序技术的数据分析。目前为止，广泛引用的高通量测序平台主要有三种1，Illumina公司的Genome Analyzer (又称为Solexa)测序平台。2，ABI公司的SoLiD测序平台。3，Roche公司的454 测序平台。三种测序平台分别基于不同的测序技术并各有特点，但目前最为广泛应用的是 Solexa 平台。Solexa平台将待扩增并测序的DNA固定于固体表面，使用BridgePCR amplification 扩增 DNA 片段，并使用 reverse dye terminator 技术进行测序。Solexa 平 台运行一次成本约8000美元，可产生 40，000, 000左右的35_70bp的序列数据。Solexa 平台成本远低于454平台(以每bp花费计)，并且不存在SoLiD技术所存在的G/C偏差的 问题，因此在生物学研究领域得到广泛引用。目前Solexa技术主要有两部分应用1，RNA-Seq，即转录组测序。将细胞或组织内 的mRNA反转录为eDNA后，进行扩增并输入Solexa平台测序，得到的结果进行分析后可以 得到mRNA的表达量。RNA-seq技术由于拥有精确定量和高灵敏度的特点，被认为将会很快 取代Microarray技术。2，ChIP-Seq，即免疫共沉淀测序技术。这项技术可以定位转录因子 (transcription factor)与DNA的结合位点(binding site)而在生物医学研究中被广泛 应用。目前针对Solexa技术平台的数据分析软件有如下几类1，序列对位软件， 将 Solexa 测序的 reads 在基因组上快速定位[Langmead, B.，et al.，Ultrafast and memory-efficient alignment of short DNA sequencesto the human genome. Genome Biol，2009. 10(3) :p. R25]。2，RNA-seq分析软件，根据RNA-seq的数据确定每个基因的表 达量[Mortazavi, A. , et al. , Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods, 2008. 5 (7) :p. 621-8]。3，ChlP-seq 分析软件，将 ChlP-seq 的结果角军析为转录因子结合位点(transcription factor binding site) [Rozowsky, J. , et al., PeakSeq enables systematic scoring of ChlP—seqexperiments relative to controls. Nat Biotechnol,2009. 27(1) :p. 66-75]。本专利致力于后两类应用，并提出了全新的分析 思路以提高分析结果的质量。
发明内容
本发明基于目前的Solexa测序技术，找到了一种新的可以定义差异表达和转录 因子结合位点的分析方法，相对于其他分析方法大幅提高了分析精度。本方法的步骤如 下(1)获取Solexa测序序列，将所有序列对参照基因组使用ELAND软件进行对位 (Alignment)。将无法对位的序列(如测序质量太低的序列)丢弃。对于有多重对位的序 列(multiple hits)保留分数最高或并列最高的对位结果。(2)将获得的对位序列文件转化成转化成RMI (Read Mass Index) Score。RMI的计 算方法如下RMI = (Read Coverage/Mappability)^Adjustment其中Read Coverage为该位点被测序的次数，我们可以利用对位文件直接计算出 精确到每bp的Read Coverage。Mappability表示该区段在零假设下被随机序列覆盖的理论值。该理论值及其 分布取决于参考基因组，无法用理论公式计算，但我们可以利用参考基因组(Reference Genome)进行计算机模拟计算得出。其计算的方法是将参考基因组拆分为35bp (或者 70bp，取决于Solexa测序的长度)的小段，以Ibp为步长，将每一个理论上可能的区段都取 出，然后将所有的小段都对原基因组对位。如此得到的对位结果即为Mappability的理论 分布。显然，基因组中的重复序列的Mappbility将比唯一的序列Mappability高，这也是 我们在计算RMI时要对Mappability进行校正的原因。Adjustment为针对该次测序的校正。Solexa测序的过程中存才测序误差，因此 并非所有的序列都可以完美对位(perfect match)到参考基因组上。有一些序列将有Ibp 的误差(lbp mismatch)，另有一些序列有2bp的误差(多于2bp误差的序列将不予考虑)。 本方法对有mismatch的序列有一定罚分，即认为这些序列的可信度比完美对位(perfect match)的序列要低。经过试验，本方法将lbp mismatch的序列可信度设为50 %，2bp mismatch的序列可信度设为25%。(3)经过以上步骤之后，我们得到了全基因组范围内的RMIindex。接下来的步骤 是计算RMI的理论分布。接下来将分为两种情况讨论A，RNA-seq分析。B，Chip-seq分析。(A)RNA-Seq分析。RNA-seq分析相对较为简单。一般来说，我们的实验设计为对 比两个样本，或者一系列时间序列的样本互相比较。我们通过步骤(2)已经得到实际RMI 分布，现在需要计算的是RMI的理论分布。在给定区段内，该分布将是一个二项分布 在理论情况下，每个基因的表达量已知(假定为Ε)，该分布中η为所测序列的总数(Total number of reads)，ρ为该区段的序列数的理论值 (MappabiIity/ Σ Mappabi 1 ity)*η*Ε。然而基因有各自不同的表达量，实验之前无法得知。 因此我们选取实验中的一个样本，计算出E和相应的ρ，定义它为标准值。其他的样本计算 出的数值则与该标准值进行统计检验。若出现显著差异，则表明存在显著差异。为了进行 以上统计检验，我们需要知道该分布的标准差。我们知道二项分布的标准差为P (I-P)。公 式中最后一个参数k为此处实际的序列数(number of observed reads)。为简化后面的计 算，当序列数较大时(大于等于30)，二项分布可以用正态分布模拟 该分布的平均值和标准差和前一个分布相同，μ = (MappabiIity/ Σ MappabiIity)σ "2 = μ *(1-μ )。如前所述，得到RMI的理论分布之后，假设我们要对两个样本的RNA-seq结果进行 差异基因表达分析，则我们设定其中一个样本为标准，方差已知(理论值)，所以用另一样 本进行标准ζ检验，得到p-value。如果有多个生物重复(biological r印licates)，则可 以进行t检验。(4)为控制假阳性结果(false positive results)，在进行基因组范围实验时，多 重校正必不可少。这里我们将对步骤(3)得到的p-value进行Bonferroni多重检验校正。 假定我们以步长X，区间1对整个基因组进行t检验，假定基因组长度为L。则我们进行的 检验总数为N= (L-y)/x。根据此检验次数，我们认为只有p-value小于0.05/N的结果方 为显著结果。(B)ChlP-seq分析。ChlP-seq的分析原理与RNA-seq基本相同。唯一不同的步骤 是ChIP-seq实验以Input样本与ChIP_ed样本进行比较。因此我们以Input样本为标准， 设定其在某区段的表达量为E，相应的二项分布中的参数为P，其标准差为P* (I-P)。其后的 步骤与RNA-seq相同。本方法的特征本方法以RMI指数为基础，设计了 一种新的分析方法，适用于RNA-seq和ChlP-seq 两大类高通量测序实验的分析。该方法比现存的其他方法更具弹性，并有更高的精确性。 具体体现在1)分析RNA-seq结果时，不是基因为单位，而是以区段为单位。区段的大小可 以任意定义，因此本方法可以同时检测差异表达与可变剪切，这是本方法与其他方法的显 著区别。2)在分析ChlP-seq序列时，只要做少许优化，(比如，将结合位点定义为可以是 p-value最小化的区间)，则本方法将给出精确到碱基的峰值，将比其他方法更加精确。这 是RMI方法比其他ChlP-seq分析方法的区别所在。以上是对本发明的描述而非限定，基于本发明思想的其它实施方式，均在本发明 的保护范围之中。
权利要求
基于RMI指数的高通量测序分析新方法，其特点是基于一个新的指数(RMIRead Mass Index)，对高通量测序结果进行快速准确的分析。该方法的特征在于有如下步骤步骤1获取转录本高通量测序(RNA-seq)或染色体免疫共沉淀高通量测序(ChIP-seq)数据。步骤2根据高通量测序信息，并基于物种基因组序列信息进行校正，估算RMI的经验分布。步骤3利用RMI经验分布和测序数据，鉴定差异表达区段或差异结合峰。
2.RMI的计算方法为RMI = (Read Coverage/Mappability)^Adjustment其中Read Coverage为该位点被测序的次数Mappability表示该区段在零假设下被随机序列覆盖的理论值。RMI的计算方法如下RMI = (Read Coverage/Mappability)^Adjustment其中Read Coverage为该位点被测序的次数，我们可以利用对位文件直接计算出精确 到每 bp 的 Read Coverage。Mappability表示该区段在零假设下被随机序列覆盖的理论值。计算的方法是将参考 基因组拆分为η个碱基(η值为不同的测序仪器读取序列的长度)的小段，以1碱基为步长， 将每一个理论上可能的区段都取出，然后将所有的小段都对原基因组对位。如此得到的对 位结果即为Mappability的理论分布。Adjustment为针对该次测序的校正。即对有不匹配的序列有一定罚分，即认为这 些序列的可信度比完美对位的序列要低。经过试验，本方法将1个碱基不匹配的序列 adjustment值设为50%，2个碱基不匹配的序列adjustment值设为25%。
3.在高通量测序仪器Solexa，454(GS-FLX) ,SOLiD和Polonator中均可以使用RMI指 数进行数据分析。
全文摘要
本发明基于一个新的指数RMI(Read Mass Index)，对高通量测序结果进行快速准确的分析。该方法同时适用于转录本高通量测序(RNA-seq)和染色体免疫共沉淀高通量测序(ChIP-seq)两种方法产生的数据进行分析。
文档编号C12Q1/68GK101886114SQ200910051239
公开日2010年11月17日 申请日期2009年5月14日 优先权日2009年5月14日
发明者周杰, 曾华宗 申请人:上海聚类生物科技有限公司