CN101974623A - 差异表达基因的检测方法 - Google Patents

Abstract

本发明提供一种差异表达基因的检测方法，包括如下步骤：确定样本的数量；通过数据读取单元，从基因芯片中获取基因表达值；通过数据转换单元进行数据转换；通过数据处理单元，根据样本的不同条件，按双样本检验单元、单样本检验单元、多样本检验单元进行数据处理；通过分布识别单元进行分析，判断样本检验结果的分布，从而确定p值；通过比较单元判断基因表达值是否存在差异，得到检测结果。本发明检测灵敏度高，假阳性低，有效性好。

Description

差异表达基因的检测方法

技术领域

本发明涉及生物领域中基因芯片分析检测技术，具体涉及一种差异表达基因的检测方法。

背景技术

DNA微阵列(DNA microarray)，也叫基因芯片，是最近数年发展起来的一种能快速、高效检测DNA片段序列、基因表达水平的新技术。它将数目从几百个到上百万个不等的称之为探针的核苷酸序列固定在小的(约1cm²)玻璃或硅片等固体基片或膜上，该固定有探针的基片就称之为DNA微阵列。根据核苷酸分子在形成双链时遵循碱基互补原则，就可以检测出样本中与探针阵列中互补的核苷酸片段，从而得到样本中关于基因表达的信息，这就是基因表达谱，因此基因表达谱可以用一个矩阵或一个向量来表示，矩阵或向量元素的数值大小即该基因的表达水平。

随着大规模基因表达谱(Gene expression profile，或称为基因表达分布图)技术的发展，各类病人的特异组织基因表达分布图都可以参照各种组织的正常的基因表达。从DNA芯片所测量的成千上万个基因中，找出决定样本类别的一组基因“标签”，即“信息基因”(informative genes)是正确识别肿瘤类型、给出可靠诊断和简化实验分析的关键所在，同时也为抗病药物的研制提供了一种新的途径。通常由于基因数目很大，在判断疾病基因标签的过程中，需要剔除掉大量“无关基因”，从而大大缩小需要搜索的致病基因范围。在删除大量“无关基因”时，需要利用各种检验方法。国内外对于两总体样本实验，当前检验方法需要样本独立的假设，两样本T检验以及基于T检验改进检验方法如：SAM，贝叶斯检验，都仅仅考虑了数据来自正态分布，它仅仅检验均值是否有差异。而F检验以及其改进检验都仅仅检验方差是否有差异。如果需要检验分布是否相同，可以用非参数检验如Wilcoxon检验，但是要检验分布是否相同在当前还没有参数检验方法。对于多总体样本实验，当前方法是方差分析(ANOVA)或者非参数Kruskal-wallis检验。ANOVA只考虑到多个总体样本的均值，忽视偏度、峰度而且无法考虑到方差、偏度、峰度是否相同。非参数Kruskal-wallis检验虽没有分布假设，但是它只利用了秩(排序后序列号大小)的信息，非参数方法是无法用参数方法或者没有参数方法的时候才考虑。

在实际中，系统是非线性和相互联系的系统。在非线性系统，各种参数(如：均值、方差、偏度和峰度)是相互影响，特别是偏度和峰度被作为非线性重要指标[5]并且在一个非常弱的系统中偏度和峰度也将被保持甚至是扩大。在非线性系统中，即使输入信号为正态分布，它的输出信号也将不再是正态分布。而基因之间是相互作用的，这种作用可以是正或负反馈循环，同时基因相互作用是非线性的。因此，偏度和峰度是不能忽视的。

当前的检验方法不仅忽视偏度和峰度从而忽略了基因之间的非线性相互作用，而且不能同时检测均值、方差、偏度和峰度等是否有差异。但在实际中要同时检验均值、方差、偏度和峰度等更高中心矩的差异性几乎是不可能的，(1)多个中心矩无法进行合并，因为偏度和峰度之间不独立，(2)即使合并了，其分布也相当复杂，很难写出分布函数从而计算p值，(3)即使其分布可以找到，在实际中因为样本量比较小得到具有偏的中心矩估计，从而使最终检验结果无法接受。

本发明克服现有技术的以上缺点，提供一种差异表达基因的检测方法，检测结果表明本发明检测的假阳性低，有效性高，可以检测到其他方法不能检测得出的差异表达基因。

发明内容

本发明提供一种差异表达基因的检测方法，其特征在于，包括如下步骤：

(1)确定样本的数量；所述样本包括病例组样本、对照组样本、或多总体中每个总体的样本；

(2)根据样本的数量，通过数据读取单元，从基因芯片中获取基因表达值；

(3)先通过2进制-10进制数据转换单元，再通过以2为底对数的数据转换单元，将所述基因表达值转换成矩阵X_i，其中，矩阵X_i是第1到第k次方矩阵，

i表示第i个样本，x表示数据转换后的值，k表示最高原点矩的次数，w表示获取样本的个数；

(4)通过数据处理单元进行样本检验；

A)当所述样本包括病例组样本和对照组样本时，通过双样本检验单元进行SWang检验，所述SWang检验是指

其中，n是病例组样本的数量，m是对照组样本的数量，k是最高原点矩的次数，

是病例组样本的基因表达值转换后的原点矩向量，

是对照组样本的基因表达值转换后的原点矩向量，D₁是病例组样本的基因表达值转换后的从第1次方到第k次方之间离差阵，D₂是对照组样本的基因表达值转换后的从第1次方到第k次方之间离差阵；

B)当所述样本仅包括病例组样本时，通过单样本检验单元进行SWang检验，所述SWang检验是指

其中，S₀是已知分布的原点矩向量；

C)当所述样本为多总体样本时，通过多样本检验单元进行SWang检验，所述SWang检验是指其中，d为多总体的数量，d_i为各个总体中的样本数量，sd为d个总体的样本数量之和，sd＝d₁+d₂+…+d_p；

(5)根据分布识别单元进行分析，判断上述样本检验结果的分布，从而确定p值；

(6)通过比较单元得到检测结果，当所述p值小于0.05时，则判断所述样本之间的基因表达值存在差异。

本发明中所述步骤(2)中基因芯片包括实验制得的基因芯片。

本发明中所述步骤(4)中的k值最大为4。

本发明中所述步骤(5)中SWang检验符合F分布或皮尔逊分布。

本发明中所述步骤(6)的样本之间的基因表达值的差异，可以通过公知基因数据库或实验来验证。

本发明的多总体，是指来自多于两个总体，且每个总体的样本量可以相同也可以不相同。

本发明提供一种差异表达基因的检测方法，是通过检验样本之间的基因分布是否相同来判断其是否属于差异表达基因。而要检验样本之间的基因分布是否相同，则需要同时检验其原点矩是否相等。其工作原理在于，当F_X(x)和F_Y(y)分别是累积分布函数，同时这两个累计分布函数的所有原点矩都存在且是有限的，要F_X(u)与F_Y(u)对所有的u相等当且仅当E(X^r)＝E(Y^r)对任意r＝0，1，2，…都相等，其中，E(X^r)是原点矩。

本发明通过利用中心矩与原点矩之间相互关系进行变换，从而只需要同时检验多个原点矩是否有差异就可以来判断中心矩是否有差异，其原因在于，(1)原点矩在样本比较小时也可以无偏的，(2)他们联合分布可以得到，同时也可以得到证明。本发明还利用矩阵把多个原点矩结合到一起，这样可以写出分布也可以把相关统计量的分布转化到熟悉的F分布，从而可以得出其p值。

基因相互之间的作用可以有正或负反馈循环的调节，同时基因相互作用是非线性的。如存在基因A，B，C三个基因，A对B、C有调节作用，而B对C有调节作用，C对A有调节作用。实际上，A基因发生变化，而B、C也将发生相应的变化，同时B的变化将使C也相应变化，而后C也发生变化。这样的变化最后达到一个稳定状态。非线性就是说如果A变化了2，而B和C相应变化0.125和8，这是三次方非线性，但是基因达到20000多个，这种非线性之间关系非常复杂，要找到他们非线性函数相当困难。而在过去研究证明了偏度和峰度可以检测非线性是否存在和测量非线性大小。本发明SWang检验检测非线性系统信号的多个原点矩差异性，是基于非线性相互关联的考虑，其中，峰度和偏度是测量非线性大小的重要指标。本发明可以判断在非线性系统信号中的数据是否来自某一个特定分布。首先根据分布计算出原点矩，然后计算出数据原点矩，然后检验二者是否相等来判断数据是否来自同一特定分布。反之，也可以检验相关数据是否来自于某一个特定分布，从而判断相关数据的特定原点矩是否相等。

本发明通过数据处理单元进行样本检验。本发明利用矩阵广义逆而不是矩阵逆。本发明SWang检验方法中，(D₁+D₂)^-表示广义逆，因为在(D₁+D₂)的行列式不为0时，广义逆与其逆相同，但是在行列式为0时，以广义逆因为其逆不存在。

本发明是基于基因表达数据的分布存在各个原点矩。所述SWang检验与F分布存在对应关系，即SWang～F(k，n+m-k-1)，其自由度分别为k和n+m-k-1。因此，本发明通过分布识别单元，根据F分布，判断SWang检验结果的分布，从而确定p值。

在双样本情况下，假设第一组有b样本，记为X＝{x₁₁，x₁₂，…，x_1b}，第二组有c样本Y＝{y₁₁，y₁₂，…，y_1c}，所用数据的分布是皮尔逊家族分布，皮尔逊家族分布中包括有均值、方差、偏度和峰度。

是x₁₁，x₁₂，…，x_1b的偏度，

是x₁₁，x₁₂，…，x_1b的峰度。

是y₁₁，y₁₂，…，y_1c的偏度，是y_i1，y_i2，…，y_ic的峰度。

对x₁₁，x₁₂，…，x_1m来自正态分布，正态分布的均值为mu和标准差为sigma。首先，利用公式(x_ij-mu)/sigma对数据进行标准化，但是在标准化后偏度和峰度将保持不变.因此，偏度和峰度与正态分布的参数mu和sigma是独立的。因为

是正态分布mu的估计同时是正态分布sigma的平方的诡计。最后，偏度和峰度是与

和

独立的。

由于直接利用中心阶矩会在检验中遇到一些问题：比如，中心阶矩在样本很小时，通过样本估计会产生偏差；中心阶矩的联合分布，相当复杂，无法推导；利用所有中心阶矩没有比较好的方法进行同时检验。本发明的检测方法采用中心阶矩与原点矩变换，

$r_o=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^o---\left(1\right)$

r₁＝0                              (2)

$r_2=-\frac{S_1^2}{n^2}+\frac{S_2}{n}---\left(3\right)$

$r_3=\frac{{6S}_1^3}{n^3}-\frac{3S_1{S}_3}{n^2}\frac{S_3}{n}---\left(4\right)$

$r_4=-\frac{{3S}_1^4}{n^4}+\frac{{6S}_1^2{S}_2}{n^3}-\frac{{4S}_1{S}_3}{n^2}+\frac{S_4}{n}---\left(5\right)$

$r_5=\frac{{4S}_1^5}{n^5}-\frac{{10S}_1^3{S}_2}{n^4}+\frac{{10S}_1^2{S}_3}{n^3}-\frac{10S_1{S}_4}{n^2}+\frac{S_5}{n}---\left(6\right)$

……

是原点矩，o＝1，2，3，...，k，r₁，r₂，r₃，r₄，和r₅是中心阶矩。同时，原点矩可以在样本很小的情况下得到无偏估计，这时样本只要大于4。

根据Cramer theory，如果X来自正态分布，均值向量为μ和协方差矩阵为∑，则

其中g是一个影射

如偏微分函数在μ∈R^d领域是连续的。则

根据Cramer theory，假设前8阶矩都存在，根据中心极限定理得到：

$\left(\begin{array}{c}{S}_1\\ S_2\\ S_3\\ S_4\end{array}\right)\stackrel{L}{~}N(\left(\begin{array}{c}{\mathrm{\μ}}^1\\ {\mathrm{\μ}}^2\\ {\mathrm{\μ}}^3\\ {\mathrm{\μ}}^4\end{array}\right),\left(\begin{array}{cccc}\mathrm{Var}\left(X\right)& \mathrm{Cov}(X^2,X)& \mathrm{Cov}(X^3,X)& \mathrm{Cov}(X^4,X)\\ \mathrm{Cov}(X^2,X)& \mathrm{Var}\left(X^2\right)& \mathrm{Cov}(X^3,X^2)& \mathrm{Cov}(X^4,X^2)\\ \mathrm{Cov}(X^3,X)& \mathrm{Cov}(X^3,X^2)& \mathrm{Var}\left(X^3\right)& \mathrm{Cov}(X^4,X^3)\\ \mathrm{Cov}(X^4,X)& \mathrm{Cov}(X^4,X^2)& \mathrm{Cov}(X^4,X^3)& \mathrm{Var}\left(X^4\right)\end{array}\right))---\left(7\right)$

Hoteling T²检验是用来检验多元正态分布均值向量是否相等，同时检验多个相关变量均值是否相等。该检验不仅考虑到每个变量均值的差异，还考虑到变量之间的相关系数。

X来自多元正态分布，其均值向量为

其p×p协方差矩阵为∑且∑未知，Y来自多元正态分布，其均值向量为

其p×p协方差矩阵为∑。则Hoteling T²检验为

其中A₁A₂是离差矩阵。

本发明中，原点矩样本均值向量

和

可以通过以下进行估计：

${\stackrel{\‾}{X}}_i=\left[\begin{array}{c}\stackrel{\‾}{x_i^1}\\ \stackrel{\‾}{x_i^2}\\ .\\ .\\ .\\ \stackrel{\‾}{x_i^k}\end{array}\right],$ ${\stackrel{\‾}{Y}}_i=\left[\begin{array}{c}\stackrel{\‾}{y_i^1}\\ \stackrel{\‾}{y_i^2}\\ .\\ .\\ .\\ \stackrel{\‾}{y_i^k}\end{array}\right]---\left(11\right)$

其中

得到，原点矩样本离差矩阵D_xi和D_Yi

$D_{\mathrm{Xi}}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i){(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i)}^{\′},$ $D_{\mathrm{Yi}}=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}(Y_{\mathrm{ij}}-{\stackrel{\‾}{Y}}_i){(Y_{\mathrm{ij}}-{\stackrel{\‾}{Y}}_i)}^{\′}---\left(12\right)$

其中 $X_{\mathrm{ij}}={[x_{\mathrm{ij}}^1,x_{\mathrm{ij}}^2,...,x_{\mathrm{ij}}^k]}^{\′},$ $Y_{\mathrm{ij}}={[y_{\mathrm{ij}}^1,y_{\mathrm{ij}}^2,...,y_{\mathrm{ij}}^k]}^{\′}.$

同时，可以得到原点矩样本协方差矩阵

可以得到，原点矩样本相关系数矩阵

其中v_ij是其中o＝1，2，3，…，k，q＝1，2，3，…，k。

本发明中的

是正态分布，其中均值为0，方差为

$(\frac{1}{n}+\frac{1}{m})\mathrm{\Σ}---\left(15\right)$

因此，

$\sqrt{\frac{\mathrm{nm}}{m+n}}(\stackrel{\‾}{X}-\stackrel{\‾}{Y})---\left(16\right)$

是多元正态分布，其中均值为0向量，协方差为∑。

$A_1=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(X_{\left(j\right)}-\stackrel{\‾}{X}){(X_{\left(j\right)}-\stackrel{\‾}{X})}^{\′},$ $A_2=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(Y_{\left(j\right)}-\stackrel{\‾}{Y}){(Y_{\left(j\right)}-\stackrel{\‾}{Y})}^{\′}---\left(17\right)$

A₁和A₂是Wishart分布，其自由度都为n-1，协方差矩阵是∑。

因为Wishart分布可加性。

∑W_i～W(h，∑)，W_i～W(H_i，∑)，h＝∑H_i(18)

因此，

A₁+A₂～W(2n-2，∑)(19)

另一方面，T²是Hotelling’sT²分布，其自由度分别为p和n+m-2。根据T²分布和F分布，本发明中可以转换T²检验为F检验。

$F=\frac{(n+m-2)-p+1}{(n+m-2)*p}{T}^2~F(p,n+m-p-1)---\left(20\right)$

本发明还可以采用SWang(h，k)检验，其中h是最小阶矩，k是最大阶矩。SWang(h，k)可以利用来自h阶矩到k阶矩的信息，例如SWang(1，3)利用来自第一阶矩到第三阶矩的信息进行检验。因此，阶矩需要多少依赖分布和样本个数。如果k等于1、SWang(1，1)是T检验的平方，就是F检验，但是k大于2。

SWang(h，k)

$\mathrm{SWang}=\frac{(n+m-2)-(k-1+1)+1}{n+m-2}{T}^2$

$=\frac{(n+m-2)-(k-1+1)+1}{(n+m-2)*(k-1+1)}\frac{\mathrm{nm}}{(n+m)}{({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{\left(\frac{D_1+D_2}{n+m-2}\right)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

$=\frac{(n+m-2)-(k-1+1)+1}{(n+m-2)*(k-1+1)}\frac{\mathrm{nm}}{(n+m)}(n+m-2){({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{(D_1+D_2)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

$={\frac{(n+m-k-1)\mathrm{nm}}{(n+m)k}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{(D_1+D_2)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

$\mathrm{SWang}(h,k)=\frac{(n+m-2)-(k-h+1)+1}{n+m-2}{T}^2$

$=\frac{(n+m-2)-(k-h+1)+1}{(n+m-2)*(k-h+1)}\frac{\mathrm{nm}}{(n+m)}{({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{\left(\frac{D_1+D_2}{n+m-2}\right)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

$=\frac{(n+m-2)-(k-h+1)+1}{(n+m-2)*(k-h+1)}\frac{\mathrm{nm}}{(n+m)}{(n+m-2)({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{(D_1+D_2)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

$=\frac{(n+m-k+h-2)\mathrm{nm}}{(n+m)*(k-h+1)}{({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{(D_1+D_2)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

如果k＝h＝1，SWang就是T检验的平方，其中方差为未知且相等。

$\mathrm{SWang}(h,k)=\frac{(n+m-2)-(1-1+1)+1}{n+m-2}{T}^2$

$=\frac{(n+m-2)-(1-1+1)}{(n+m-2)*(1-1+1)}\frac{\mathrm{nm}}{(n+m)}{({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)}^{\′}{\left(\frac{d_1+d_2}{n+m-2}\right)}^{-1}({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)$

$=\frac{(n+m-2)-(1-1+1)+1}{(n+m-2)*(1-1+1)}\frac{\mathrm{nm}}{(n+m)}(n+m-2){({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)}^{\′}{(d_1+d_2)}^{-1}({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)$

$=\frac{(n+m-2)\mathrm{nm}}{(n+m)*(1-1+1)}{({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)}^{\′}{(d_1+d_2)}^{-1}({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)$

$=\frac{(n+m-2)\mathrm{nm}}{(n+m)}{({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)}^{\′}{(d_1+d_2)}^{-1}({\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2)$

$={\left(\frac{{\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2}{\sqrt{\frac{\mathrm{nm}}{m+n}}\left(\frac{d_1+d_2}{n+m-2}\right)}\right)}^2$

$={\left(\frac{{\stackrel{\‾}{x}}_1-{\stackrel{\‾}{x}}_2}{\sqrt{\frac{\mathrm{nm}}{m+n}}\left(\frac{(n-1)s_1+(m-1)s_2}{n+m-2}\right)}\right)}^2$

本发明中检测多个总体样本时，假设有d个总体，第i个总体的样本量为d_i，d个总体样本量之和为sd，sd＝d₁+d₂+…+d_d。每个总体的样本进行变换

$X_i={\left[\begin{array}{cccc}{x}_{i1}& x_{i2}& ...& x_{{\mathrm{in}}_i}\\ x_{i1}^2& x_{i2}^2& ...& x_{{\mathrm{in}}_i}^2\\ ...& ...& ...& ...\\ x_{i1}^k& x_{i2}^k& ...& x_{{\mathrm{in}}_i}^k\end{array}\right]}_{k\×n_i}$

其为第i个总体多个原点矩均值向量，

${\stackrel{\‾}{X}}_i=\left[\begin{array}{c}\stackrel{\‾}{x_i^1}\\ \stackrel{\‾}{x_i^2}\\ .\\ .\\ .\\ \stackrel{\‾}{x_i^k}\end{array}\right]$

是p个总体的原点均值向量，

$\stackrel{\‾}{X}=\frac{1}{\mathrm{\Σ}{n}_i}\left[\begin{array}{c}\mathrm{\Σ\Σ}{x}_i^1\\ \mathrm{\Σ\Σ}{x}_i^2\\ .\\ .\\ .\\ \mathrm{\Σ\Σ}{x}_i^k\end{array}\right]$

$T=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}(X_{\mathrm{ij}}-\stackrel{\‾}{X}){(X_{\mathrm{ij}}-\stackrel{\‾}{X})}^{\′}$

$=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i+{\stackrel{\‾}{X}}_i-\stackrel{\‾}{X}){(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i+{\stackrel{\‾}{X}}_i-\stackrel{\‾}{X})}^{\′}$

$=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i){(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i)}^{\′}+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X}){({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X})}^{\′}$

$=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{D}_i+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{n}_i({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X}){({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X})}^{\′}$

其中， $D_i=\underset{j=1}{\overset{n_i}{\mathrm{\Σ}}}(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i){(X_{\mathrm{ij}}-{\stackrel{\‾}{X}}_i)}^{\′}$

$\mathrm{SWang}=(\frac{k-1-(d-1)+1}{2}-(\mathrm{sd}-k))\log \left(\frac{\left|\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{D}_i\right|}{|\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{D}_i+\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{d}_i({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X}){({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X})}^{\′}|}\right)$

以上表明，

因此，

则此时，SWang～χ²((k-1)(d-1))。最后，通过转变后SWang检验统计量服从卡方分布，其自由度是(k-1)(d-1)。

附图说明

图1(a)是多个样本下真实假阳性与估计假阳性的示意图。

图1(b)显示图1(a)中一个样本下真实假阳性与估计假阳性的比较示意图。

图2是RT-PCR实验结果的图谱。

图3本发明差异表达基因的检测方法的流程框图。

图4本发明检测方法中数据处理单元的流程示意图。

具体实施方式

结合以下具体实施例来进一步详细地说明本发明。本发明差异表达基因的检测方法，包括如下步骤：

(1)确定样本的数量；

样本包括病例组样本、对照组样本、或多总体中每个总体的样本。在病例组-对照组实验中，确定病例组和对照组的样本量。在只有病例组的、没有对照组的情况下，可以确定病例组样本量以及标准范围，在标准范围内的认为是正常，否则认为是生病的。标准范围的选择可以根据以往研究或者是自己的经验值。

(2)根据样本的数量，通过数据读取单元，从基因芯片中获取基因表达值；

确定待研究的针对性疾病类型，从生物实验室产生生物芯片或者NCBI网站下载对应的生物芯片数据(在生物芯片中，基因个数一般为2000到50000之间)，根据样本的数量，通过数据读取单元，获得相关基因表达值。基因芯片可以通过商业渠道购得、实验制得、或公开数据库中的基因芯片。本发明的基因芯片是指所有用于检测不同条件下基因表达谱差异的芯片杂交技术。

(3)先通过2进制-10进制数据转换单元，再通过以2为底对数的数据转换单元，将所述基因表达值转换成矩阵X_i，其中，矩阵X_i是第1到第k次方矩阵，

i表示第i个样本，x表示数据转换后的值，k表示最高原点矩的次数，w表示获取样本的个数；

(4)通过数据处理单元进行样本检验；

A)当所述样本包括病例组样本和对照组样本时，通过双样本检验单元进行SWang检验，所述SWang检验是指

其中，n是病例组样本的数量，m是对照组样本的数量，k是最高原点矩的次数，

是病例组样本的基因表达值转换后的原点矩向量，

是对照组样本的基因表达值转换后的原点矩向量，D₁是病例组样本的基因表达值转换后的从第1次方到第k次方之间离差阵，D₂是对照组样本的基因表达值转换后的从第1次方到第k次方之间离差阵；

B)当所述样本仅包括病例组样本时，通过单样本检验单元进行SWang检验，所述SWang检验是指

其中，S₀是已知分布的原点矩向量；

C)当所述样本为多总体样本时，通过多样本检验单元进行SWang检验，所述SWang检验是指

其中，d为多总体的数量，d_i为各个总体中的样本数量，sd为d个总体的样本数量之和，sd＝d₁+d₂+…+d_p；

(5)根据分布识别单元进行分析，判断样本检验结果的分布，从而确定p值；

(6)通过比较单元，分析得到检测结果，当所述p值小于0.05时，则判断所述样本之间的基因表达值存在差异。如果p值大于0.05或0.01，则认为不存在差异表达。

检测得到的差异表达基因，可以通过到NCBI网站或者通过PubMed等公共数据库中检索结果，或者通过生物实验结果(如：PCR实验)，来验证选取的基因是否存在差异表达，是否与要研究的疾病有关。

在通过数据处理单元时，根据样本的不同条件，按不同的检验单元及流程进行检验：

A：两样本情况下

(1)确认样本量m，n

(2)先对数据进行一次方，平方，三次方，四次方直到到k次方构造一个数据矩阵A

(3)计算出数据矩阵A的各个次方平均向量和各个次方协方差矩阵

(4)代入 $\mathrm{SWang}={\frac{(n+m-k-1)\mathrm{nm}}{(n+m)k}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)}^{\′}{(D_1+D_2)}^{-1}({\stackrel{\‾}{X}}_1-{\stackrel{\‾}{X}}_2)$

(5)根据F分布表或者直接利用SAS软件中PROBF求出p-value

(6)计算出联合置信域 $[{\stackrel{\‾}{S}}_i-c\sqrt{D_{1_{\mathrm{ii}}}/\left[m(m-1)\right]+D_{2_{\mathrm{ii}}}/\left[n(n-1)\right]},{\stackrel{\‾}{S}}_i+c\sqrt{D_{1_{\mathrm{ii}}}/[m(m-1)+D_{2_{\mathrm{ii}}}/[\left]n(n-1)\right]}]$

其中， $c=\sqrt{\frac{(n+m)k}{(n+m-k-1)\mathrm{nm}}\mathrm{SWang}}$

(7)根据p值进行判断，如果p值大于阀值则认为没有差异，反之，则认为有差异。

B：单样本情况下

(1)确定非线性系统产生机理，根据机理判断输出数据的分布，算出各个原点矩

组成S₀，同时确认样本量m

(2)先对数据进行一次方，平方，三次方，四次方直到到k次方构造一个数据矩阵A

(3)计算出数据矩阵A的各个次方平均向量和各个次方协方差矩阵

(4)代入 $\mathrm{SWang}={\frac{n(n-1-k+1)}{k}({\stackrel{\‾}{X}}_1-S_0)}^{\′}{\left(D_1\right)}^{-1}({\stackrel{\‾}{X}}_1-S_0)$

(5)根据F分布表或者直接利用SAS软件中PROBF求出p值

(6)计算出联合置信域 $[{\stackrel{\‾}{S}}_i-c\sqrt{D_{\mathrm{ii}}/\left[m(m-1)\right]},{\stackrel{\‾}{S}}_i+\mathrm{cc}\sqrt{D_{\mathrm{ii}}/[m(m-1)}]$

其中， $c=\sqrt{\frac{(n-1)k}{(n-k)}\mathrm{SWang}}$

(7)根据p值进行判断，如果p值大于阀值则认为没有差异，否则认为有差异。

C：检验特定原点矩

(1)确认样本量m，n，确定特定原点矩的次方数{o1，o2，o3，...}

(2)先对数据进行{o1，o2，o3，...}次方转换成一个数据矩阵A

(3)双样本按A流程，单样本按B流程

(4)判断是否有差异，如果有差异，差异值是多少。

D：当多个总体样本的情况

(1)确认总体个数为d，各个总体样本量 d_i，计算出总样本量sd＝d₁+d₂+…+d_d

(2)先对数据进行一次方，平方，三次方，四次方直到到k次方构造一个数据矩阵A

(3)计算出数据矩阵A的各个次方平均向量和各个次方离差矩阵

(4)代入 $\mathrm{SWang}=(\frac{k-1-(d-1)+1}{2}-(\mathrm{sd}-k))\log \left(\frac{\left|\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{D}_i\right|}{|\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{D}_i+\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{d}_i({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X}){({\stackrel{\‾}{X}}_i-\stackrel{\‾}{X})}^{\′}|}\right)$

(5)根据卡方分布表或者直接利用SAS软件中PROBCHI求出p值

(6)根据p值，如果p值大于阀值则认为没有差异，否则认为有差异。

实施例1：

在病例组和对照组都由相同正态分布产生随机基因表达数据，重复10000次，然后利用二项分布以0.2概率产生0或1，“1”被定为真实表达基因，“0”被定为真实的没有表达基因。最后利用F检验、T检验、本实施例中样本检验计算各个基因的统计量并计算p值，如果小于真实假阳性则认为基因差异表达，否则认为没有差异表达。计算假阳性得到估计假阳性，得到图1。在图中如果线段越趋向横坐标则说明该方法越好。

图1所示为多个样本下真实假阳性与估计假阳性，SS是样本个数，蓝色是SWang检验，黑色和灰色分别是T检验和F检验。图1表明本发明的假阳性比较低，相对与T检验和F检验非常好。各图横坐标是真实假阳性，纵坐标是估计假阳性，SS表示样本量。图1显示本发明检测的曲线总体比其他曲线更趋向横坐标。

实施例2：

生物芯片可以同时测量生物几乎所有基因的变化情况，而基因之间是相互关联的，同时这种关联是非线性的，因此基因表达谱数据的分布几乎不是正态分布。而如T检验等检验方法要求分布是正态，在生物是非常难的。本实施例不需要一定是正态分布，只要来自的分布存在原点矩就可以，非常适合这种情况下的检测。

本实施例中确定数据所在位置http://research.nhgri.nih.gov/microarray/NEJMSupplement，在这里可以下载数据，然后直接打开需要的原始数据。

其中，病例组有8个样本(BRAC2)，对照组有7个样本(BRCA1)，且基因个数为3226个。

本实施例中样本检验为样本检验的统计量服从F分布

算出自由度p值；通过p值小于0.05认为基因差异表达。如表1所示，其他方法没有检测出的差异表达基因，可以被本发明的方法检测出。本实施例还通过文献可以验证，如表2所示，还通过PCR实验加以证实，如图2所示。用RT-PCR实验证实了本发明差异表达基因的检测方法，检测到新的与breast癌症有关的基因。图2是RT-PCR实验结果的图谱，其中1表示细胞器MCF-10A，2表示细胞器SK-BR-3，3表示细胞器MCF-7，4表示细胞MDA-MB-231。以BCR基因为例加以说明，当用其他方法检测认为BCR没有差异表达时，在图谱上则不显示条带，但是本实施例中BCR在图谱上显示有亮条，从而说明本发明的检测方法可以检测出BCR存在差异表达。以IDS为例也是同样情况。

表1为比较不同检验方法下基因表达值。

基因标签为mutation_Id，gene是基因名字，SK1是对照组基因表达的偏度，SK2是疾病组基因表达的偏度，KU1是对照组基因表达的峰度，KU2是疾病组基因表达的峰度，T_T是T检验统计量，pt是T检验的统计量的P值，F_T是F检验统计量，pf是F检验统计量的P值。SWang是本发明所采用的样本检验方法，P_sw是本发明样本检测得到的P值。

表1.不同检验方法下基因表达值

以上基因在乳腺癌症研究中具有重要的作用，比如，TP53。通过表1可以得到，通过T检验并且F检验表明没有差异性表达的基因，却可以通过本发明的检测出存在差异性表达基因。在T检验和F检验的p值都大于0.05，也就是说，从均值、方差看，这些基因没有差异性，而本发明的p值都小于0.05，表明从偏度和峰度等方面这些基因存在比较大的差异，比如，病例组中偏度和峰度分别为2.2864和5.5197，对照组中偏度和峰度分别为1.1775和-.3575，偏度之差为1.11，峰度之差为5.88。以上表时，基因表达值在对照组和疾病组的分布是不同。

表2本发明检测得到的差异表达基因由相关公开文献证实的情况

实施例3：

首先，产生各个分布的数据，然后按照单样本情况下和对特定原点矩情况下，分别进行检验。

本实施例中，数据来自标准正态分布，据此，计算出正态分布的第一原点矩为0，第二原点矩为1，第二原点矩为0，第四原点矩为3。表3.1，是检验不同分布在样本为30时是否可以判断这些分布为正态分布，从而来说明本发明方法的有效性。表3.2，是检验不同分布在样本为30时，其第4和第5原点矩是否分别等于3和0。

表3.1.不同分布下的检验

通过表3.1表明，本发明可以用来检验分布是否是某一个假设的分布，在假设为标准正态分布的情况下，就可以计算出相应的各个原点矩，通过比较，表明用指数分布或者伽玛分布产生的数据拒绝了原假设(分布来自正态分布)。

表3.2.不同分布下的检验

通过表3.2表时，本发明可以用来检验特定原点矩，在假设为标准正态分布的第四和五原点矩分别为3和0，通过比较表明，用指数分布或者伽玛分布产生的数据拒绝了原假设(第四原点矩等于3同时第五原点矩等于0)。

实施例4：

假设所有总体样本来自同一分布，进行样本检验。

首先，从多个分布产生随机数，每个分布是一个总体，每个样本量为30，然后按照D。这些分布是正态分布、指数分布、伽玛分布和均匀分布。

表4在多个总体分布下，方差分析、Kruskal-walis检验和本发明SWang检验比较

通过表4表明，本发明可检测出多个总体分布是否相同。如果分布相同时，本发明检测的统计量为0，因此其p值为1，方差分析统计量为0，其p值为1，Kruskal-walis检验统计量为0，其p值也为1，完全接受多个分布总体是相同的。在分布不同时，本发明检测方法的p值小于0.05，表明是可以用来检验多个总体分布是否相同。但是方差分析检验和kruskal-walis检验统计量存在大于0.05。

比如，以下4组：“1个均匀分布(-1，1)，1个均匀分布(-0.5，0.5)，1个均匀分布(-3，3)，1个均匀分布(-5，5)”；“1个均匀分布(-1，1)，1个均匀分布(-0.9，0.9)，1个均匀分布(-0.8，0.8)，1个均匀分布(-1.1，1.1)”；“1个均匀分布(-1，1)，1个均匀分布(-0.9，0.9)，1个均匀分布(-0.8，0.8)，1个混合分布(均匀分布(-1.1，1.1)+标准正态分布)”；“1个正态分布(均值为0，方差为0.25)、1个正态分布(均值为0，方差为1)、1个正态分布(均值为0，方差为2.25)、1个正态分布(均值为0，方差为4)”；在这四组中，分布是完全不同的，但是利用方差分析和Kruskal-walis检验则认为相同，而本发明检测得到的p值小于0.05，则认为这四组数据的分布是完全不同。

Claims (5)

1.一种差异表达基因的检测方法，其特征在于，包括如下步骤：

(1)确定样本的数量；所述样本包括病例组样本、对照组样本、或多总体中每个总体的样本；

(2)通过数据读取单元，从基因芯片中获取基因表达值；

(3)先通过2进制-10进制数据转换单元，再通过以2为底对数的数据转换单元，将所述基因表达值转换成矩阵X_i，其中，矩阵X_i是第1到第k次方矩阵，i表示第i个样本，x表示数据转换后的值，k表示最高原点矩的次数，w表示获取样本的个数；

(4)通过数据处理单元进行样本检验；

A)当所述样本包括病例组样本和对照组样本时，通过双样本检验单元进行SWang检验，所述SWang检验是指

其中，n是病例组样本的数量，m是对照组样本的数量，k是最高原点矩的次数，

是病例组样本的基因表达值转换后的原点矩向量，

是对照组样本的基因表达值转换后的原点矩向量，D1是病例组样本的基因表达值转换后的从第1次方到第k次方之间离差阵，D2是对照组样本的基因表达值转换后的从第1次方到第k次方之间离差阵；

B)当所述样本仅包括病例组样本时，通过单样本检验单元进行SWang检验，所述SWang检验是指

其中，S₀是已知分布的原点矩向量；C)当所述样本为多总体样本时，通过多样本检验单元进行SWang检验，所述SWang检验是指其中，d为多总体的数量，d_i为各个总体中的样本数量，sd为d个总体的样本数量之和，sd＝d₁+d₂+…+d_p；

(5)通过分布识别单元进行分析，判断样本检验结果的分布，从而确定p值；

(6)通过比较单元得到检测结果，当所述p值小于0.05时，则判断所述样本之间的基因表达值存在差异。

2.如权利要求1所述差异表达基因的检测方法，其特征在于，所述步骤(2)中基因芯片包括实验制得的基因芯片。

3.如权利要求1所述差异表达基因的检测方法，其特征在于，所述步骤(4)中的k值最大为4。

4.如权利要求1所述差异表达基因的检测方法，其特征在于，所述步骤(5)中样本检验结果是符合F分布或卡方分布。

5.如权利要求1所述差异表达基因的检测方法，其特征在于，所述步骤(6)的样本之间的基因表达值的差异，可以通过公共基因数据库中已经发表的被证实的基因或实验来验证。

Images