专利名称：确定差异表达基因的混合线性模型两步筛选法的制作方法
技术领域：
本发明属于大规模基因表达数据处理技术领域，特别涉及一种高效准确鉴别差异表达基因的方法。
背景技术：
在患病人群中，某些基因发生突变和功能失调等原因会导致其本身基因表达水平的变异，这样，通过基因芯片技术分析和比较患病人群和正常人群中基因表达谱，就有可能发现在两种人群中一些表达不一致(即有差异)的基因。这些基因可能与某特定疾病有关，因此，可用来做进一步功能分析、药物设计和临床诊断等应用。同样，也可通过基因芯片技术分析和比较服药人群和正常人群以鉴别某特定药物的药效等等。
随着基因芯片制作技术的提高和基因芯片价格的下降，该技术在基因功能分析、致病机制研究、药物毒性分析、中药药理机制和药效成分研究及医疗诊断等方面得到了广泛应用。准确分析基因芯片测定结果的需求大增。然而，现有的分析基因芯片测定结果的数据处理方法所获得的结果不很可靠，因为基于基因芯片技术的大规模基因表达谱测定过程中存在多种误差来源，如芯片制作点样差异，探针标记杂交温度和湿度变化，不能准确地鉴别真正由于考察因素如药物处理或疾病等引起差异表达的基因，难以辨别真实的基因表达变化和随机的实验误差，容易导致假阴性和假阳性。根据包含大量假阴性和假阳性的基因得到的聚类分析难以得到具有实际生物学意义的科学结论。
因此，本领域中迫切需要开发出一种能克服上述缺陷，降低分析结果的假阳性概率的方法。
发明内容
为解决上述问题，本发明提供了一种鉴别差异表达基因的方法，其特征在于，该方法通过两个步骤鉴定差异表达基因先将芯片数据通过噪音过滤消除大的试验系统误差，然后通过单基因模型初步判断符合显著性标准A的差异表达基因；其次，用多基因模型分析这些初定的差异表达基因，选出符合显著性标准B的基因；其中显著性标准A和B之间的关系是A/B≥5，以便控制假阳性。
在一个较佳的实施方案中，在所述用多基因模型分析差异表达基因的步骤中，剖分出基因主效应及基因与各项考察因素的互作效应，用MINQUE法估计各项效应的方差和协方差分量，用AUP法预测随机效应，在较严格的标准下鉴定差异表达基因控制假阳性。
在另一方案中，在单基因模型分析步骤中，用MINQUE法估计各项变异来源的方差分量σu(l)2，用Jackknife抽样法估算方差分量估计值的抽样方差；针对每个基因检验该基因的处理效应的方差是否为0；如果拒绝第l基因的H0，则该基因被初步判定为可能的差异表达基因，对应的该基因的观察值将进入多基因模型进一步分析。
在另一方案中，在多基因模型中用MINQUE法估计模型的各项方差分量，用AUP法预测随机效应大小；Jackknife抽样技术结合t检验用来检验多基因模型包括的可能的差异表达基因是否真正为处理引起的差异表达，也就是基因与处理的互作效应是否与0有显著差异；对于第l个基因至少有一个基因与处理互作效应不为0，则将基因l作为差异表达基因。
具体而言，本发明的方法包括以下步骤1)通过噪音过滤消除大的试验系统误差；2)通过单基因模型检测考察因素方差变异的大小在一个比较宽松的标准下初步筛选差异表达基因；3)用多基因模型分析这些数目相对较少的初定的差异表达基因，可以剖分出基因主效应及基因与各项考察因素的互作效应，用MINQUE法估计各项效应的方差和协方差分量，用AUP法预测随机效应，在较严格的标准下鉴定差异表达基因控制假阳性；4)采用Jacknife重复抽样技术估算检验统计量的显著性；5)将基因和考察因素的互作效作为鉴定差异表达基因的具体指标；6)对基因主效应的无偏估计值和基因与处理因素互作效应的无偏预测值进行进一步的聚类分析，得到相对真实的具有生物学意义的科学结论。
本发明通过对大规模基因芯片表达谱数据，首先应用本发明所提出的单基因模型完成差异表达基因鉴别的初筛选，然后再对初选到的基因用多基因模型进行分析，根据基因和考察因素的互作效应这个指标判定差异表达基因。实验证明，该方法的分析结果非常可靠，能大大降低分析结果的假阳性概率，大幅度提高大规模基因表达数据分析的效率。
本发明的经济效益和社会效益分析，基因芯片分析发现的差异表达基因可进一步用于基因功能分析和药物设计等领域，经济效益和社会效益非常显著。在实际应用过程中，必须尽量确保这些发现的基因是具有真实的差异表达(即不是假阳性的基因)。对假阳性的基因进行功能分析会造成资源的浪费，同样也会导致药物设计的失败和临床诊断的错误等。在一般基因芯片试验中，通常能发现与研究目标(如寻找某种癌症致病机理或某种药物药效等)有关的差异表达基因约200个，本发明能在大多情形下有效降低假阳性的基因至少5％至10％(即10至20个基因)，这样一次应用试验可节省费用＝10～20个×10万元/基因＝100～200万元。同时，本发明能更有效地发现新的传统方法未能发现的基因2～5％(即4至10个基因)，这样一次应用试验可增加经济效益＝4～10×20万元/基因＝80～200万元。累计节省费用和经济效益随着基因芯片技术的推广和应用而不断增加。
附图简述
图1.三种方法鉴定差异表达基因的判错率(FDR)及功效(power)比较，图中虚线表示判错率，实线表示功效，三种方法为混合线性模型(圆圈)，无过滤的t-检验(三角)及变异倍数过滤加t-检验(方块)。VGT/VP为基因与处理的互作方差占总表型变异的比例，其值从0.1变化到0.9。EQUAL为效应A，D，T，GA，GD和ε的方差相同，即VA∶VD∶VT∶VGA∶VGD∶Vε＝1∶1∶1∶1∶1∶1；ARRAYDOM表示效应A和GA在剩余的表型变异中占绝对优势，即(VA+VGA)/(VP-VGT)＝0.9，VD∶VT∶VGD∶Vε＝1∶1∶1∶1；DYEDOM表示效应D和GD在剩余的表型变异中占绝对优势，即(VD+VGD)/(VP-VGT)＝0.9，VA∶VT∶VGA∶Vε＝1∶1∶1∶1；TREATDOM表示效应T在剩余的表型变异中占绝对优势，即VT/(VP-VGT)＝0.9，VA∶VD∶VGA∶VGD∶Vε＝1∶1∶1∶1∶1。
图2显示了三种方法鉴定差异表达基因的判错的基因个数(FN)比较。图中三种方法为混合线性模型(圆圈)，无过滤的t-检验(三角)及变异倍数过滤加t-检验(方块)。
图3示意性地显示了本发明整个方法的分析流程。
具体实施例方式
本发明克服了现有分析数据处理方法的缺点，设计了一种基于混合线性模型的两步法分析基因芯片数据的方法。该方法可以有效地将基因表达水平根据变异来源的不同剖分为几个组成部分。其中根据不同的研究试验设计，模型可以做适当的调整。通过两个步骤鉴定差异表达基因首先，将经过数据标准化的芯片数据通过噪音过滤消除大的试验系统误差，然后在一个比较宽松的标准下通过单基因模型初步判断差异表达基因；其次，用多基因模型分析这些初定的差异表达基因以便在较严的标准下控制假阳性。用最小范数二次无偏估计方法(MINQUE法)估计各项效应的方差和协方差分量，用调整预测方法(AUP法)预测随机效应。例如对于一个不同药物处理对基因表达影响的研究，噪音过滤模型为yijkl＝μ+Ai+Tj+Dk+εijkl，其中yijkl是实验获得的第i张芯片上的第j种药物处理的标记为第k种荧光的第l个基因的表达量；μ是所有基因的平均表达量；Ai是芯片效应，Ai～(0，σA2)；Tj是处理效应，Tj～(0，σT2)；Dk是荧光标记效应，Dk～(0，σD2)；εijkl是随机误差，εijkl～(0，σε2)。通过最小二乘估计法可获得以上几个参数的估计值 和 由公式rijkl=yijkl-(μ^+A^i+T^j+D^k)]]>计算得到该模型的残差可作为单基因模型的观察值。对於第l个基因的单基因模型可以表示为rijkl＝μl+Ail+Tjl+Dkl+γijkl，其中μi、Ail、Tjl、Dkl和γijkl分别是第l个基因的平均表达量、芯片效应、处理效应、荧光标记效应和随机误差。用MINQUE(1)法可估计各项方差分量σu(l)2。用Jackknife抽样法估算方差分量估计值的抽样方差。针对每个基因检验该基因的药物处理效应的方差是否为0即H0:σT(l)2=0,]]>H1:σT(l)2≠0.]]>如果拒绝第l基因的H0，则该基因被初步判定为可能的差异表达基因，对应的该基因的观察值将进入多基因模型进一步分析。多基因模型可以表述为yijkl＝μ+Gl+Ai+Tj+Dk+GAli+GTlj+GDlk+εijkl，其中Gl是第l个基因的平均表达量，GAli是第l个基因与第i张芯片的互作效应，GTlj是第l个基因与第j种药物处理的互作效应，GDlk是第l个基因与第k种荧光标记的互作效应，其他参数定义如同噪音过滤模型。用MINQUE(1)法可估计模型的各项方差分量。用AUP(Adjusted UnbiasedPrediction)法预测随机效应大小。Jackknife抽样技术结合t检验用来检验多基因模型包括的可能的差异表达基因是否真正为药物处理引起的差异表达，也就是基因与药物处理的互作效应是否与0有显著差异。对于第l个基因至少有一个基因与处理互作效应不为0，则基因l为差异表达基因。两个步骤中的筛选标准可以根据用户的具体追求目标来确定如果用户追求低假阳性，则采用严格的筛选标准(即较小的显著性测验水平)；如果用户追求低假阴性，则采用宽松的筛选标准(即较大的显著性测验水平)。在我们为用户提供的与本发明配套的软件中可供用户选择的显著性水平有0.05，0.01，0.005，0.001，0.0005，0.0001和0.00005。在下面的实例数据分析中两个步骤采用的显著性水平分别为0.0005和0.0001。
蒙特卡罗模拟结果表明混合线性模型两步法在绝大多数情况下优于传统的t检验方法(图1和2)。当A，D，T，GA，GD和ε的方差相等(EQUAL)时，混合线性模型两步法检验差异表达基因的功效比t检验法高，同时判错率也比t检验法低。当A和GA效应的方差在除GT效应外的剩余表型变异中占绝对优势(ARRAYDOM)时，混合线性模型两步法获得比t检验法显著高的功效，同时判错率与t检验法相当或比其低。当D和GD效应方差占绝对优势(DYEDOM)时，混合线性模型两步法仍可获得比t检验法显著高的功效；当VGT/VP大于0.3，判错率比t检验法略高。当T效应的方差在除GT效应外的剩余表型变异中占绝对优势(TREATDOM)时，t检验法的功效比混合线性模型两步法高，但其付出的代价是异常高的判错率。在上述四种情况下，混合线性模型两步法鉴定差异表达基因时判错的基因个数总是比t检验法少，尤其是当T效应的方差在除GT效应外的剩余表型变异中占绝对优势的情况，t检验法鉴定出的差异表达基因中大约2500到3000个是假阳性或假阴性，而混合线性模型两步法鉴定出的差异表达基因中只有4到40个是假阳性或假阴性。
采用基因和处理因素的互作效应作为鉴定差异表达基因的具体指标，大大降低了传统的基於比值法或t测验法的假阴性和假阳性。对基因主效应的无偏估计值和基因与处理因素互作效应的无偏预测值进行聚类可以获得具有统计学和生物学意义的结果。发明人同时用实例数据检验该方法的有效性和实用性。采用的实例数据来源于斯坦福大学的DNA芯片开放数据库(Stanford Microarray Database；http://genome-www5.stanford.edu/；芯片编号为11374、11333、11339、11323、11375和11342)。该实例是用来研究模式生物拟南芥的信使RNA的降解模式。采用本发明提出的方法找到了一些用传统的数据处理方法不能发现的新基因。本发明可以供广大生物医学科研机构和基因芯片公司及制药企业分析基因芯片实验数据，准确可靠地鉴定差异表达基因。
下面将根据具体实施例来进一步描述本发明。然而，应当理解，本发明的方法并不仅限于该实施例。
数据来源于斯坦福大学的DNA芯片开放数据库(Stanford Microarray Database；http//genome-www5.stanford.edu/；芯片编号为11374、11333、11339、11323、11375和11342)。该实例目的是研究模式生物拟南芥的信使RNA的降解模式，共有11521个克隆。试验者做了三次独立的3’-脱氧腺苷处理试验，每个试验样本分别在两张芯片杂交，其中一张芯片上红色荧光标记未处理的样品，绿色荧光标记3’-脱氧腺苷处理了120分钟的样品，另外一张芯片正好相反，绿色荧光标记未处理的样品，红色荧光标记3’-脱氧腺苷处理了120分钟的样品。原文采用t检验法鉴定差异表达基因，即信使RNA不稳定的基因。在显著性水平为0.0001，共检测到100个基因的信使RNA不稳定。
采用本发明提出的方法我们找到了一些用传统的数据处理方法不能发现的新基因。具体过程如下
针对这批数据，我们所用的噪音过滤模型为yijkl＝μ+Ai+Tj+Dk+εijkl，其中yijkl是实验获得的第i张芯片上的第j种处理的标记为第k种荧光的第l个基因的表达量；μ是所有基因的平均表达量；Ai是芯片效应，Ai～(0，σA2)；Tj是处理效应，Tj～(0，σT2)；Dk是荧光标记效应，Dk～(0，σD2)；εijkl是随机误差，εijkl～(0，σε2)；i为1，2，3，4，5或6，即芯片编号；j为1或2，分别代表未处理或3’-脱氧腺苷处理了120分钟的样品；k为1或2，分别代表红色荧光标记或绿色荧光标记；l为基因编号，从1到11521。通过最小二乘估计法获得以上几个参数的估计值 和 由公式rijkl=yijkl-(μ^+A^i+T^j+D^k)]]>计算得到该模型的残差作为单基因模型的观察值。对於第l个基因的单基因模型为rijkl＝μl+Ail+Tjl+Dkl+γijkl，其中μl、Ail、Tjl、Dkl和γijkl分别是第l个基因的平均表达量、芯片效应、处理效应、荧光标记效应和随机误差。用MINQUE(1)法估计各项方差分量σu(l)2。用Jackknife抽样法估算方差分量估计值的抽样方差。针对每个基因检验该基因的处理效应的方差是否为0即H0:σT(l)2=0,]]>H1:σT(l)2≠0.]]>如果拒绝第l基因的H0，则该基因被初步判定为可能的差异表达基因，对应的该基因的观察值进入多基因模型进一步分析。我们采用的多基因模型为yijkl＝μ+Gl+Ai+Tj+Dk+GAli+GTlj+GDlk+εijkl，其中Gl是第l个基因的平均表达量，GAli是第l个基因与第i张芯片的互作效应，GTlj是第l个基因与第j种药物处理的互作效应，GDlk是第l个基因与第k种荧光标记的互作效应，其他参数定义如同噪音过滤模型。用MINQUE(1)法估计模型的各项方差分量。用AUP(Adjusted UnbiasedPrediction)法预测随机效应大小。用Jackknife抽样技术结合t检验检验多基因模型包括的可能的差异表达基因是否真正为3’-脱氧腺苷处理引起的差异表达，也就是基因与处理的互作效应是否与0有显著差异。对于第l个基因至少有一个基因与处理互作效应不为0，则基因l为差异表达基因。为了与原作者的分析结果比较，我们在单基因模型和多基因模型中采用的显著性水平分别为0.0005和0.0001。
运用混合线性模型两步法，共检测出90个基因在3’-脱氧腺苷处理了120分钟后有明显的信使RNA降解现象，其中的51个基因用t检验法也检测到了，另外的39个基因仅用混合线性模型两步法检测出了(见表1)。AA395830和N37328是位点At1g13260上的DNA结合蛋白基因RAV1的两个表达序列标签(ESTs)，H77088和T04337是位点At2g22430上的同源异形域转录因子基因(ATHB-6)的两个EST。混合线性模型两步法检测出这4个EST在3’-脱氧腺苷处理后都发生了降解，而t检验法只鉴别出N37328和T04337。AA720100，AA720105和T76004都来自于位点At4g30440的核苷糖差向异构酶基因；T20600，N65459和T75944都来自于位点At4g31500的细胞色素P450单氧酶基因。t检验法只鉴别出AA720100和T20600是不稳定的，混合线性模型两步法检测出AA720105，T76004，N65459和T75944是容易降解的基因。T20543，AA720239和AA720240这3个ESTs都来自于位点At5g64260上的基因，它们仅被混合线性模型两步法检测出是不稳定基因AA067525和AA067498来自于位点At5g37770的钙调蛋白2基因；AA597715和H36178来自于位点At5g61590上的似乙烯响应元素组合因子基因；AA597849和T46143来自于位点At1g72450上的基因。两种方法分别在每个基因上找到一个转录物。但是t检验法没有象混合线性模型两步法那样找到来自于同一个基因的多个转录物。因为同一个基因编码的探针应该具有相似的表达模式，所以显而易见在同一个基因中同时发现多个不稳定的转录物是非常可能的。从这一点来说，混合线性模型两步法鉴定出了较多的合理的不稳定转录物。
另外，在信使RNA翻译过程中，polyA尾巴对于增强信使RNA的稳定性和保证信使RNA的正常功能具有重要作用。研究表明缺少polyA尾巴的组蛋白信使RNA的半衰期少于30分钟。混合线性模型两步法发现了两个与组蛋白有关的ESTs(H76940和AA720291)，但t检验法没有鉴定出来。
表1.用混合模型方法鉴别具有不稳定转录本的拟南芥(Arabidopsis thaliana)基因(AtGUTs)
*混合模型方法和t-检验方法将ESTs鉴别成差异表达的基因
权利要求
1.一种鉴别差异表达基因的方法，其特征在于，该方法包括先将芯片数据通过噪音过滤消除大的试验系统误差，然后通过单基因模型初步判断符合显著性标准A的差异表达基因；其次，用多基因模型分析这些初定的差异表达基因，选出符合显著性标准B的基因；其中显著性标准A和B之间的关系是A/B≥5，以便控制假阳性；采用Jackknife重复抽样技术估算检验统计量的显著性，在单基因模型中将处理考察因素的方差作为初定差异表达基因的具体指标，在多基因模型中将基因和考察因素的互作效应作为具体指标来鉴定差异表达基因。
2.根据权利要求1所述的方法，其特征在于，在用多基因模型分析差异表达基因的步骤中，剖分出基因主效应及基因与各项考察因素的互作效应，用MINQUE法估计各项效应的方差和协方差分量，用AUP法预测随机效应，在较严格的标准下鉴定差异表达基因控制假阳性。
3.根据权利要求1所述的方法，其特征在于，在单基因模型分析步骤中，用MINQUE法估计各项变异来源的方差分量σu(l)2，用Jackknife抽样法估算方差分量估计值的抽样方差；针对每个基因检验该基因的处理效应的方差是否为0；如果拒绝第l基因的H0，则该基因被初步判定为可能的差异表达基因，对应的该基因的观察值将进入多基因模型进一步分析。
4.根据权利要求1所述的方法，其特征在于，在多基因模型中用MINQUE法估计模型的各项方差分量，用AUP法预测随机效应大小；Jackknife抽样技术结合t检验用来检验多基因模型包括的可能的差异表达基因是否真正为处理引起的差异表达，也就是基因与处理的互作效应是否与0有显著差异；对于第l个基因至少有一个基因与处理互作效应不为0，则将基因l作为差异表达基因。
5.根据权利要求1所述的方法，其特征在于，该方法包括以下步骤1)通过噪音过滤消除大的试验系统误差；2)通过单基因模型检测考察因素方差变异的大小在一个比较宽松的标准下初步筛选差异表达基因；3)用多基因模型分析这些数目相对较少的初定的差异表达基因，可以剖分出基因主效应及基因与各项考察因素的互作效应，用MINQUE法估计各项效应的方差和协方差分量，用AUP法预测随机效应，在较严格的标准下鉴定差异表达基因控制假阳性；4)采用Jacknife重复抽样技术估算检验统计量的显著性；5)将基因和考察因素的互作效应作为鉴定差异表达基因的具体指标；6)对基因主效应的无偏估计值和基因与处理因素互作效应的无偏预测值进行进一步的聚类分析，得到相对真实的具有生物学意义的科学结论。
全文摘要
一种鉴别差异表达基因的方法，其特征在于，该方法通过两个步骤鉴定差异表达基因先将芯片数据通过噪音过滤消除大的试验系统误差，然后通过单基因模型初步判断符合显著性标准A的差异表达基因；其次，用多基因模型分析这些初定的差异表达基因，选出符合显著性标准B的基因；其中显著性标准A和B之间的关系是A/B≥5，以便控制假阳性。采用Jackknife重复抽样技术估算检验统计量的显著性，在单基因模型中将处理考察因素的方差作为初定差异表达基因的具体指标，在多基因模型中将基因和考察因素的互作效应作为鉴定差异表达基因的具体指标。本发明的方法能够提供可靠的分析结果，大大降低分析结果的假阳性概率，并大幅度提高大规模基因表达数据分析的效率。
文档编号C12Q1/68GK1786192SQ20041008920
公开日2006年6月14日 申请日期2004年12月8日 优先权日2004年12月8日
发明者朱军, 陆燕, 李亦学 申请人:李亦学, 朱军, 陆燕