CN103177197A - 基于高通量测序检测差异表达与可变剪切分析的方法 - Google Patents

Abstract

本发明基于目前的Solexa测序技术，找到了一种新的可以定义差异表达和转录因子结合位点的分析方法，相对于其他分析方法大幅提高了分析精度。该方法同时适用于转录本高通量测序(RNA-seq)和染色体免疫共沉淀高通量测序(ChlP-seq)两种方法产生的数据进行分析。

Description

基于高通量测序检测差异表达与可变剪切分析的方法

技术领域

本发明属于生物技术领域，涉及新一代高通量测序技术数据分析。

背景技术

本发明是一种适用于高通量测序技术的新分析方法。高通量测序技术是近两年来生物技术领域的重要突破，新一代的测序技术将传统的Sanger测序效率提高了数百倍，同时价格也大大下降。高通量测序技术的出现使得许多的极有前景的生物医药应用成为可能：1，癌症基因组。2，个性化医疗与诊断。3，药物靶标筛选。高通量测序技术能否在这些领域的取得进展，其关键在于分析方法及软件的创新。本专利提出了一种新的分析方法，可广泛用于高通量测序技术的数据分析。

Solexa平台将待扩增并测序的DNA固定于固体表面，使用Bridge PCR amplification扩增DNA片段，并使用reverse dyeterminator技术进行测序。Solexa平台运行一次成本约8000美元，可产生～40,000,000左右的35-70bp的序列数据。Solexa平台成本远低于454平台(以每bp花费计)，并且不存在SoLiD技术所存在的G/C偏差的问题，因此在生物学研究领域得到广泛引用。

目前Solexa技术主要有两部分应用：1，RNA-Seq，即转录组测序。将细胞或组织内的mRNA反转录为cDNA后，进行扩增并输入Solexa平台测序，得到的结果进行分析后可以得到mRNA的表达量。RNA-seq技术由于拥有精确定量和高灵敏度的特点，被认为将会很快取代Microarray技术。2，ChIP-Seq，即免疫共沉淀测序技术。这项技术可以定位转录因子(transcription factor)与DNA的结合位点(binding site)而在生物医学研究中被广泛应用。

目前针对Solexa技术平台的数据分析软件有如下几类：1，序列对位软件，将Solexa测序的reads在基因组上快速定位。2，RNA-seq分析软件，根据RNA-seq的数据确定每个基因的表达量。3，ChIP-seq分析软件，将ChIP-seq的结果解析为转录因子结合位点(transcription factor binding site)[Rozowsky，J.，et al.，PeakSeqenables systematic scoring of ChIP-seq experiments relative to controls.Nat Biotechnol，2009.27(1)：p.66-75]。本专利致力于后两类应用，并提出了全新的分析思路以提高分析结果的质量。

发明内容

本发明基于目前的Solexa测序技术，找到了一种新的可以定义差异表达和转录因子结合位点的分析方法，相对于其他分析方法大幅提高了分析精度。本方法的步骤如下：

(1)获取Solexa测序序列，将所有序列对参照基因组使用ELAND软件进行对位(Alignment)。将无法对位的序列(如测序质量太低的序列)丢弃。对于有多重对位的序列(multiple hits)保留分数最高或并列最高的对位结果。

(2)将获得的对位序列文件转化成转化成RMI(Read MassIndex)Score。RMI的计算方法如下：

RMI＝(Read Coverage/Mappability)＊Adjustment

其中Read Coverage为该位点被测序的次数，我们可以利用对位文件直接计算出精确到每bp的Read Coverage。

Mappability表示该区段在零假设下被随机序列覆盖的理论值。该理论值及其分布取决于参考基因组，无法用理论公式计算，但我们可以利用参考基因组(Reference Genome)进行计算机模拟计算得出。其计算的方法是：将参考基因组拆分为35bp(或者70bp，取决于Solexa测序的长度)的小段，以1bp为步长，将每一个理论上可能的区段都取出，然后将所有的小段都对原基因组对位。如此得到的对位结果即为Mappability的理论分布。显然，基因组中的重复序列的Mappbility将比唯一的序列Mappability高，这也是我们在计算RMI时要对Mappability进行校正的原因。

Adjustment为针对该次测序的校正。Solexa测序的过程中存才测序误差，因此并非所有的序列都可以完美对位(perfect match)到参考基因组上。有一些序列将有1bp的误差(1bp mismatch)，另有一些序列有2bp的误差(多于2bp误差的序列将不予考虑)。本方法对有mismatch的序列有一定罚分，即认为这些序列的可信度比完美对位(perfect match)的序列要低。经过试验，本方法将1bp mismatch的序列可信度设为50％，2bp mismatch的序列可信度设为25％。

(3)经过以上步骤之后，我们得到了全基因组范围内的RMIindex。接下来的步骤是计算RMI的理论分布。接下来将分为两种情况讨论：A，RNA-seq分析。B，Chip-seq分析。

(A)RNA-seq分析。RNA-seq分析相对较为简单。一般来说，我们的实验设计为对比两个样本，或者一系列时间序列的样本互相比较。我们通过步骤(2)已经得到实际RMI分布，现在需要计算的是RMI的理论分布。在给定区段内，该分布将是一个二项分布：

$f_X\left(x\right)=P[X=k]=\left(\begin{array}{c}n\\ k\end{array}\right)p^k{(1-p)}^{n-k}\mathrm{for\; k}=\mathrm{0,1,2},...,n$

在理论情况下，每个基因的表达量已知(假定为E)，该分布中n为所测序列的总数(Total number of reads)，p为该区段的序列数的理论值：(Mappability/∑Mappability)＊n＊E。然而基因有各自不同的表达量，实验之前无法得知。因此我们选取实验中的一个样本，计算出E和相应的p，定义它为标准值。其他的样本计算出的数值则与该标准值进行统计检验。若出现显著差异，则表明存在显著差异。为了进行以上统计检验，我们需要知道该分布的标准差。我们知道二项分布的标准差为p(1-p)。公式中最后一个参数k为此处实际的序列数(number of observed reads)。为简化后面的计算，当序列数较大时(大于等于30)，二项分布可以用正态分布模拟：

$f\left(x\right)=\frac{1}{\mathrm{\σ}\sqrt{2\mathrm{\π}}}{e}^{-\frac{{(1-\mathrm{\μ})}^2}{{2n}^2}}$

该分布的平均值和标准差和前一个分布相同，μ＝(Mappability/∑Mappability)＊n＊E，σ^2＝μ＊(1-μ)。

如前所述，得到RMI的理论分布之后，假设我们要对两个样本的RNA-seq结果进行差异基因表达分析，则我们设定其中一个样本为标准，方差已知(理论值)，所以用另一样本进行标准z检验，得到p-value。如果有多个生物重复(biological replicates)，则可以进行t检验。

(4)为控制假阳性结果(false positive results)，在进行基因组范围实验时，多重校正必不可少。这里我们将对步骤(3)得到的p-value进行Bonferroni多重检验校正。假定我们以步长x，区间y对整个基因组进行t检验，假定基因组长度为L。则我们进行的检验总数为：N＝(L-y)/x。根据此检验次数，我们认为只有p-value小于0.05/N的结果方为显著结果。

(B)ChIP-seq分析。ChIP-seq的分析原理与RNA-seq基本相同。唯一不同的步骤是：ChIP-seq实验以Input样本与ChIP-ed样本进行比较。因此我们以Input样本为标准，设定其在某区段的表达量为E，相应的二项分布中的参数为P，其标准差为P＊(1-P)。其后的步骤与RNA-seq相同。

本方法的特征：

本方法以RMI指数为基础，设计了一种新的分析方法，适用于RNA-seq和ChIP-seq两大类高通量测序实验的分析。该方法比现存的其他方法更具弹性，并有更高的精确性。具体体现在：1)分析RNA-seq结果时，不是基因为单位，而是以区段为单位。区段的大小可以任意定义，因此本方法可以同时检测差异表达与可变剪切，这是本方法与其他方法的显著区别。2)在分析ChIP-seq序列时，只要做少许优化，(比如，将结合位点定义为可以是p-value最小化的区间)，则本方法将给出精确到碱基的峰值，将比其他方法更加精确。这是RMI方法比其他ChIP-seq分析方法的区别所在。

以上是对本发明的描述而非限定，基于本发明思想的其它实施方式，均在本发明的保护范围之中。

Claims (1)

1.基于高通量测序分析新方法，其特点是基于一个新的指数RMI对高通量测序结果进行快速准确的分析。该方法的特征在于有如下步骤：

步骤1：获取转录本高通量测序RNA-seq)或染色体免疫共沉淀高通量测序(ChIP-seq)数据。

步骤2：根据高通量测序信息，并基于物种基因组序列信息进行校正，估算RMI的经验分布。

步骤3：利用RMI经验分布和测序数据，鉴定差异表达区段或差异结合峰。