CN105117617B - 一种用于筛选环境敏感性生物分子的方法 - Google Patents

Abstract

本发明公开了一种用于环境敏感性生物分子的筛选方法，包括以下步骤：组学数据加载和标准化处理，计算特征的四分位距和设定阈值，利用此阈值进行四分位距算法过滤，得到差异表达特征，然后，用方差分析算法对标准化后的数据集进行过滤，得到差异表达特征。结合以上两种算法，匹配共同差异表达变化的特征，进行随机森林算法排序，得到重要性的差异表达特征。在此基础上，通过聚类分析，确定环境敏感的分子标志物。本发明提供的用于筛选环境敏感的分子标志物组合算法，降低了计算时间、提高了准确性，给出了特征的重要性排序，为生物学家快速定位到环境敏感的靶分子，揭示生物的响应机理，进行早期预警和防治，提供了高效便捷的数据处理工具。

Description

一种用于筛选环境敏感性生物分子的方法

技术领域

本发明属于生物信息技术领域，涉及一种用于挖掘高通量生物信息数据中对环境具有敏感性的生物分子的方法。

背景技术

在环境毒理学领域中，生物体在有害因素(物理、化学、或者生物等因素)中暴露会引起机体的各种变化，包括生理、生化、免疫、细胞和遗传等方面的改变。随着分子生物学实验技术的发展，许多研究者已开始在分子水平上研究有害环境下生物分子的变化，识别早期生物学效应，达到对危害因素的早期预警，从而更全面地进行有害因素的风险预警，同时也可用于毒理学机制和防治效果的研究。

在分子生物学实验技术中，伴随着高通量检测技术的发展，环境毒理学研究越来越重视获取不同环境处理下的基因组、转录组、蛋白质组和代谢组等组学高通量数据，并由此从系统生物学的角度全面评估有害环境下敏感的生物分子，解析受环境胁迫的基因调控网络，为分析环境毒理学机制和防治提供重要的研究手段。目前常用的生物组学技术受人为操作以及多次实验数据整合的影响，数据的重复性不高、存在噪声以及系统偏差，预测结果往往存在由于方法而引起的无法消除的误差。因此，需要发展多种有效的分析方法，尤其是能够处理多个数据集、对噪声不敏感的统计方法，以提取海量数据中蕴含的有用信息。本领域一个前沿的研究方法是，通过整合不同实验来源的生物组学数据，减少由于样本量小造成的误差，寻找更加实用的分子标记物。

通过比较大样本的生物组学数据，寻找显著差异的分子集合，可预测潜在的分子标记物。寻找差异表达分子的计算方法很多，最直接的方法就是测量变化倍数，即计算两个样本间同一分子的表达量之比。尽管变化倍数的方法直观，但是该方法没有考虑到噪声以及生物学的可变性，带有明显的任意性，造成显著的假阳性和假阴性。比如，转录因子，其差异倍数小于2时，实际上也具有很高的生物学效应(Hu et al.,2009；Sleno and Emili,2008)。此外，倍数变化法筛选到的特征的数量一般较大，对于进一步研究生物分子的功能，带来了很大的不便。

因此，特征选择成为了差异分子筛选关键的一步，是筛选出有价值的敏感分子，进行有针对性和高效研究的重要手段。根据特征选择过程与分类器设计的关系，特征选择算法一般分为过滤法(filter)、绕封法(wrapper)和嵌入法(embedded)(Saey et al.,2007)。其中，filter法独立于分类器，wrapper法和Embedded法与分类器结合使用。目前人们已经研究了多种特征选择算法，如IQR和ANOVA法属于filter法，已经成功地被运用到基于DNA芯片数据的敏感基因的筛选(Barrell etal.,2008；Barrell etal.,2009)。随机森林(RandomForests)属于embedded法，该方法是基于决策树(decision tree)的分类器集成算法，最早由Leo Breiman提出(Breiman 2001)。通过随机置换分离点的值判断其对分类效果的影响从而决定被置换基因的重要性，并对所有基因的重要性进行排序。

然而，这些算法均有一定的优缺点，如filter法的优势在于能很快地排除大量的非关键的噪声特征，缩小特征子集的搜索范围，适合作为特征的预选器；但该法筛选到的特征规模一般比较大，其中仍包含一些明显的噪声特征，给寻找关键性的特征和分类带来很大的阻碍。wrapper和embedded法在速度上比filter法慢，但其考虑了特征之间的关系，造成所选的优化特征子集的规模相对要小得多，非常有利于关键特征的识别和筛选。因此，如何组合这些算法使其更适合筛选环境特异性的敏感分子标志物，以及对环境敏感性分子按照重要性进行排序，成为一个亟需要解决的问题。

发明内容

鉴于现有技术存在的上述问题，本发明提供一种既要考虑特征之间的关系、特征的重要性，又要减少程序运行时间的组合的特征选择算法，以用于筛选环境敏感的生物分子。

本发明通过以下技术方案实现：

一种用于环境敏感性生物分子的筛选方法，包括以下步骤：

①.读入相应生物样本表达谱数据矩阵，矩阵中的行属性为组学数据中的特征或者变量，列属性表示拟分析环境处理的条件；所述表达谱数据矩阵中的数据为对应环境处理条件下特征的绝对表达量；

②.步骤①得到的表达谱数据进行标准化处理；所述标准化处理方法为均值标准化处理；

③.计算步骤②得到的标准化的表达谱矩阵中每个特征的IQR的大小；所述IQR为统计学参数四分位距；

④.将步骤③得到的所有特征的IQR大小进行正态分布拟合分析：将正态分布90％累积概率对应的IQR设定为相应的阈值；

⑤.将步骤②得到的标准化的表达谱矩阵进行IQR算法过滤：将IQR与步骤④得到的相比较，并将IQR大于阈值的特征记录，得到IQR算法过滤下的差异表达特征；

⑥.将步骤②得到的标准化的表达谱矩阵中的列属性按照不同的处理条件进行重新分组；

⑦.将⑥得到的表达谱矩阵进行ANOVA算法过滤：将p设定为0.05，得到ANOVA过滤下的差异表达特征；所述的ANOVA为统计学中方差分析；所述p为组间差异性水平；

⑧.将步骤⑤和⑦中得到的差异表达特征相互匹配，得到共同差异表达的特征；

⑨.将步骤⑧得到共同差异表达特征矩阵进行Random Forest算法排序，得到特征的重要性：将重要性大于1.5的特征记录，得到重要的响应环境处理的特征矩阵；所述RandomForest表示随机森林算法，是生物信息学中一种用于筛选变量的数据挖掘算法；

⑩.将步骤⑨得到的响应环境处理的特征矩阵，进行聚类分析，得到响应环境处理的敏感分子标志物；所述聚类为无监督学习中的系统聚类。

进一步，在上述技术方案中，具体还包括以下步骤：步骤①中，所述的组学数据包括基因组，蛋白质组，代谢组等数据，其对应的特征分别为基因，蛋白质分子，代谢物等；所述的环境处理条件包括对照组和胁迫环境处理组；所述的组学数据样本量大于等于3；所述的环境处理条件大于等于2，其中包括1组对照处理。所述的组学数据的筛选方法为一种生物信息学组合特征选择算法。

进一步，在上述技术方案中，所述步骤②还包括，在中值标准化，离差标准化，标准差标准化和对数Logistic的标准化方法中任取一种方法；

进一步，在上述技术方案中，所述步骤④还包括，正态分布的检验；所述正态分布的检验方法为Kolmogorov-Smirnov检验；

进一步，在上述技术方案中，所述步骤⑤还包括，按照IQR过滤下的差异表达特征的顺序导出基于IQR的差异表达谱矩阵；所述的IQR算法的阈值，可以根据需要进行调整。

进一步，在上述技术方案中，所述步骤⑦还包括，按照ANOVA过滤下的差异表达特征的顺序导出基于ANOVA的差异表达谱矩阵；所述的ANOVA算法的p值，可以根据需要，在满足要求p≤0.05的要求内进行调整。

进一步，在上述技术方案中，所述步骤⑧还包括，按照共同差异表达的特征的顺序导出共同的差异表达谱矩阵；

进一步，在上述技术方案中，所述步骤⑨还包括以下步骤，

按照步骤⑥进行列属性分组；

对步骤⑧和上文所述的共同差异表达特征数据，调整Random Forest中的参数，得到最低的分类错误率，确定最佳的算法执行参数；所述的分类错误率为OOB错误率；

按照重要差异表达的特征的顺序导出重要的差异表达谱矩阵；

所述的Random Forest算法中的重要性阈值，可以根据需要进行调整。

进一步，在上述技术方案中，所述步骤⑩还包括，按照聚类分析的结果导出响应环境处理的敏感分子标志物。

与现有技术相比，本发明具有以下有益效果：

1.本发明提供了一种筛选环境敏感性生物分子的方法，组合了特征选择算法中过滤法(IQR和ANOVA算法)和封装法(Random Forest算法)的优势，可比单一的特征选择算法更好地筛选环境敏感性生物分子。

2.本发明通过生物信息学方法编写程序，可快速准确地从高通量的数据中筛选环境敏感性生物分子，操作简便，与手工分析相比，可大大减少工作量。

3.本发明能对这些特异性变化的生物分子进行关联水平的排序，为进一步研究其生物学功能和获取生物标志物提供敏感的靶分子。

附图说明

图1为本发明方法的实现流程图；

图2为IQR算法计算每个特征的IQR分布图；

图3为IQR算法计算每个特征的IQR分布的正态分布拟合图；

图4为IQR与ANOVA算法筛选到的特征比较图；

图5为Random Forest算法下参数ntree与分类错误率图；

图6为Random Forest算法筛选到的重要性的特征图；

图7为重要性特征的聚类分析图。

具体实施方式

下面主要结合附图和具体实施例对本发明的构建和结果进行详细描述。但是本领域的普通技术人员将会理解，下列实施例仅用于说明本发明，但不应视为限制本发明的范围。

下面以线虫全基因表达组为例(类似的信息在其它相关的实验中也容易获得)，对本发明进行详细说明。

线虫全基因表达组芯片：上海康成生物公司

线虫全基因组表达组芯片数据，包括九组数据。其中，该数据集包含三种线虫类型，分别为野生型，dys-1突变体，和ced-1突变体。所述三种线虫暴露在空间协同、空间辐射与地面对照等三种不同的环境处理条件下。

空间协同组(SF)：在神舟-8号飞船空间飞行环境下，线虫固定在SIMBOX盒中，该盒置于飞船内固定。空间飞行16.5天，待飞船返回后7h内收集线虫，提取得到总RNA。利用Invitrogen Superscript试剂盒合成DS-cDNA，并对其进行单色荧光标记(NimbleGen)。使用NimbleGen芯片杂交系统和Axon GenePix 4000B芯片扫描进行荧光杂交和对荧光结果扫描。最后利用NimbleScan software(version 2.5)以及Agilent GeneSpring GX software(version 11.5.1)对芯片结果进行读取和分析，获取在空间协同条件下的线虫全基因组的绝对表达量数据；

空间辐射组(SC)：在神舟-8号飞船空间飞行环境下，线虫固定在SIMBOX盒中，该盒置于飞船内带有1g离心机的设备中。在空间飞行16.5天，待飞船返回后7h内收集线虫，提取得到总RNA，按上述空间协同组相同的方法，获得在空间辐射条件下的线虫全基因组的绝对表达量数据。所述方法中带有1g离心机的设备是为了对抗微重力的影响，获得只有空间辐射的环境。

地面对照组(GC)：在地面条件下同步培养16.5天后，收集线虫，提取得到未经空间飞行环境处理的线虫总RNA，按上述空间协同组相同的方法，获得在地面条件下的线虫全基因组的绝对表达量数据，作为对照。

按照下述方法在线虫全基因组表达谱中，筛选对空间环境敏感的生物分子(基因)：

(1)图1是本发明筛选环境敏感分子标志物的方法流程图。已知线虫全基因表达谱数据矩阵。该表达谱矩阵中，行表示基因，共有18186个基因；列表示处理条件，共有9组，分别为SF(野生型)，SF(dys-1突变体)，SF(ced-1突变体)，SC(野生型)，SC(dys-1突变体)，SC(ced-1突变体)和GC(野生型)，GC(dys-1突变体)，GC(ced-1突变体)。所述的不同处理条件下括号的内容表示线虫的类型。

(2)使用read.table命令加载步骤(1)中全基因组表达谱矩阵到R语言平台。

(3)将步骤(2)中得到的全部基因组表达谱数据集按均值标准化方法进行处理，处理结果存储到预先定义的同样维度的矩阵变量中。

(4)调用“genefilter”软件包中的“rowQ”函数，将步骤(3)得到的标准化全基因组表达谱矩阵进行基因的IQR计算，得到18186个基因的IQR分布，结果如图2所示。从图中可以看出，很大比例的基因变动性很小(IQR接近于0)。

(5)将步骤(4)得到的每个基因的IQR分布进行正态分布的拟合分析和Kolmogorov-Smirnov检验。若符合正态分布，取正态分布90％累计概率对应的IQR设置为相应的阈值。图3结果显示，该数据集经所述检验后得到D＝0.05(p<0.01)，符合正态分布。在上述结果基础上，得到该设定的阈值为0.69。

(6)调用“genefilter”软件包，将步骤(5)得到的阈值用于IQR算法的筛选以得到差异表达的基因：将基因在不同环境处理条件下的IQR大于等于0.69的基因记录，并输出到.excel文件中，以供查阅。其中，将基因的名字和符号对应输出到该文件中。图4显示经过IQR算法过滤后得到1819个差异表达基因。

(7)调用“genefilter”软件包，将步骤(3)得到的标准化全基因组表达谱矩阵进行ANOVA算法过滤，其中该算法中Anova函数的参数p设定为0.05，调用该程序包中的“filterfun”和“genefilter”函数筛选得到差异表达基因：将基因在不同环境条件下的组间差异p小于该阈值的基因记录，并输出到.excel文件中，以供查阅。其中，将基因的名字，符号对应输出到该文件中。图4显示经过IQR算法过滤后得到1189个差异表达基因。

(8)将步骤(6)得到差异表达基因与步骤(7)中得到的差异表达基因中的名字相比较。若来自步骤(6)和(7)中基因的名字一致，将此基因记录，并输出到.excel文件中，以供查阅。其中，将基因的名字，符号对应输出到该文件中。图4显示IQR和ANOVA算法有86个重叠的差异表达基因。

(9)将步骤(8)得到的重叠差异表达基因进行生物功能富集分析。功能分析软件为DAVID，该软件中显著性富集的p值为0.05，结果如表1所示。对比传统的倍数变化法，所述的重叠差异表达基因富集了19个生物学过程，且大多与线粒体呼吸过程相关；而传统的方法仅富集了6个生物学过程，与生长或者有机酸代谢等相关。此外，结果显示：所述的重叠差异表达基因富集的生物学通路为氧化磷酸化通路，而传统的方法富集的生物学通路为色氨酸代谢通路。根据以往文献报道，空间环境下能导致生物体一系列的氧化应激反应，这与所述差异表达基因富集的生物学过程和通路一致。所述方法找到的空间环境敏感的分子也明显优于传统方法。

表1 本研究提出的方法和传统方法筛选到基因生物学过程富集分析

注：利用DAVID软件进行分析，且p值设定为小于0.05。

(10)调用“randomForest”软件包，将步骤(8)得到的重叠的差异表达基因矩阵进行Random Forest算法分类，图5为Random Forest算法中参数分类树的数目(ntree)与分类错误率图；结果显示，在默认参数条件下(ntree＝500和mtry＝9)，错误率为0，达到进一步分析的要求。所述的mtry表示为每个分支节点处随机抽取的变量个数。

(11)根据步骤(9)中设定的参数，调用“importance”函数，将步骤(8)得到的重叠的差异表达基因矩阵进行Random Forest算法重要性排序：将重要性大于1.5的重要性的差异表达基因记录，并输出到.excel文件中，以供查阅。其中，将基因的名字，符号对应输出到该文件中。图6为利用Random Forest算法排序得到30个重要性基因。

(12)将步骤(10)得到的基因表达矩阵进行聚类分析，所述的聚类分析方法为系统聚类，结果如图7所示。将图7中只在SF和SC中高表达的聚类基因进行筛选，得到9个对SF环境敏感的基因和8个对SC环境敏感的基因，结果如表2所示。

表1 空间辐射与空间辐射和微重力协同环境下的敏感基因

注：SF为空间辐射和微重力协同环境；SC为空间辐射环境；

本发明中生成的文件可以通过在设置的默认路径下进行查看。

本发明提供的一种用于环境敏感性生物分子的筛选方法，不仅实现了快速查找环境敏感性生物分子的功能，同时也能对该生物分子与处理环境关联的程度进行排序。该方法为环境毒理学研究提供了高效便捷的数据处理工具，对比传统的手工筛选环境敏感分子的方法大大节约了计算时间并提高了准确性。同时，该方法的功能不只局限于线虫的基因表达组，也适用于其它生物的高通量数据的差异表达分子筛选，可为今后不同生物体内高通量数据挖掘等提供有力的帮助。

以上实例仅说明本发明的一种实施方式，其描述较为具体和详细，但不能因此理解为对本发明专利范围的限制。基于本发明思想的其他实施方式，均在本发明的保护范围内。

Claims (9)

1.一种用于筛选环境敏感性生物分子的方法，其特征在于包括以下步骤：

①.读入相应生物样本表达谱数据矩阵，矩阵中的行属性为组学数据中的特征或者变量，列属性表示拟分析环境处理的条件；所述表达谱数据矩阵中的数据为对应环境处理条件下特征的绝对表达量；

②.步骤①得到的表达谱数据进行标准化处理；所述标准化处理方法为均值标准化处理；

③.计算步骤②得到的标准化的表达谱矩阵中每个特征的IQR的大小；所述IQR为统计学参数四分位距；

④.将步骤③得到的所有特征的IQR大小进行正态分布拟合分析：将正态分布90%累积概率对应的IQR设定为相应的阈值；

⑤.将步骤②得到的标准化的表达谱矩阵进行IQR算法过滤：将IQR与步骤④得到的相比较，并将IQR大于阈值的特征记录，得到IQR算法过滤下的差异表达特征；

⑥.将步骤②得到的标准化的表达谱矩阵中的列属性按照不同的处理条件进行重新分组；

⑦.将⑥得到的表达谱矩阵进行ANOVA算法过滤：将p设定为0.05，得到ANOVA算法过滤下的差异表达特征；所述的ANOVA为统计学中方差分析；所述p为组间差异性水平；

⑧.将步骤⑤和⑦中得到的差异表达特征相互匹配，得到共同差异表达的特征；

⑨.将步骤⑧得到共同差异表达特征矩阵进行Random Forest算法排序，得到特征的重要性：将重要性大于1.5的特征记录，得到重要的响应环境处理的特征矩阵；所述RandomForest表示随机森林，是生物信息学中一种用于筛选变量的数据挖掘算法；

⑩.将步骤⑨得到的响应环境处理的特征矩阵，进行聚类分析，得到响应环境处理的敏感分子标志物；所述聚类为无监督学习中的系统聚类方法。

2.根据权利要求1所述的方法，其特征在于：所述的组学数据包括基因组，蛋白质组，代谢组数据；所述的环境处理条件包括对照组和胁迫环境处理组。

3.根据权利要求1所述的方法，其特征在于：所述步骤②还包括，在中值标准化，离差标准化，标准差标准化和对数Logistic的标准化方法中任取一种方法。

4.根据权利要求1所述的方法，其特征在于：所述步骤④还包括，正态分布的检验；所述正态分布的检验方法为Kolmogorov - Smirnov检验。

5.根据权利要求1所述的方法，其特征在于：所述步骤⑤还包括，按照IQR过滤下的差异表达特征的顺序导出基于IQR的差异表达谱矩阵。

6.根据权利要求1所述的方法，其特征在于：所述步骤⑦还包括，按照ANOVA过滤下的差异表达特征的顺序导出基于ANOVA的差异表达谱矩阵。

7.根据权利要求1所述的方法，其特征在于：所述步骤⑧还包括，按照共同差异表达的特征的顺序导出共同的差异表达谱矩阵。

8.根据权利要求1所述的方法，其特征在于：所述步骤⑨还包括以下步骤，

按照步骤⑥进行列属性分组；

对步骤⑧中的共同差异表达特征数据，调整Random Forest中的参数，得到最低的分类错误率，确定最佳的算法执行参数；所述的分类错误率为OOB错误率；

按照重要差异表达的特征的顺序导出重要的差异表达谱矩阵。

9.根据权利要求1所述的方法，其特征在于：所述步骤⑩还包括，按照聚类分析的结果导出响应环境处理的敏感分子标志物。