CN103164631B - 一种智能协同表达基因分析仪 - Google Patents

Abstract

本发明公开了一种智能协同表达基因分析仪，包括基因芯片和单片机，利用基因芯片采集技术获取生物样本的基因表达谱，应用单片机嵌入式分析技术获取协同表达的基因集。单片机有四个模块组成：基因芯片表达谱读取模块；协同表达基因提取模块；存储模块；输出模块；协同表达基因提取模块的提取协同表达基因过程包含核函数选择、协同免疫克隆Memetic核双聚类算法、获取有重叠的双聚类和协同表达的基因集四部分。这种智能协同表达基因分析仪可从生物样本的基因表达谱中提取\表达趋势一致的基因集和反向表达相关的基因集。通过对这些共同表达基因的寻找，不仅可以对基因的功能研究给予提示，还可以对基因调控途径和调控网络的研究给予启发。

Description

一种智能协同表达基因分析仪

技术领域

本发明涉及一种智能协同表达基因分析仪，主要用到基因芯片非线性特征提取技术和单片机嵌入式核双聚类算法应用。

背景技术

基因芯片技术将大量探针分子固定于支持物上后与标记的样品分子进行杂交，通过检测每个探针分子的杂交信号强度，进而获取样品分子的数量和序列信息，为寻找新的功能基因发挥着重要作用。基因芯片现已泛化到生物芯片、微阵列、DNA芯片，甚至蛋白芯片。

表达谱基因芯片的出现为检测整个基因组的表达情况提供了一个极为有利的工具。一次微阵列实验能获得细胞在某一条件下的全基因组表达数据，包含成千上万个基因在细胞中的相对或绝对丰度，构成了一个数据矩阵集，行向量代表基因；列向量代表某一条件各基因的表达水平。对于基因表达谱数据的分析任务就是从数据矩阵中找出具有相似表达模式的基因(高表达基因集、低表达基因集或者相反表达基因集)。这些具有相似表达模式的基因可能具有共同的特征，如共同的调节元件，共有的生物功能，或者共有的细胞起源等。通常某一特定途径中的基因或者受相同环境变化影响的基因应当是共调控的，并且具有相似的表达模式，通过对这些共同表达基因的寻找，不仅可以对基因的功能研究给予提示，还可以对基因调控途径和调控网络的研究给予启发。基因表达谱协同表达基因的提取，对生物医学临床诊断、药物疗效判断、揭示疾病发生机制都有重要的指导意义。

应用于基因芯片表达数据的传统聚类算法可将基因表达数据看成聚类的对象,将表达模式类似的基因,也就是共表达的基因，归入同一聚类中。传统聚类是一种全局意义上的聚类，它是以全部样本作为特征，或以全部基因作为特征进行的聚类，只能在数据矩阵的行或者列某一方向上进行信息相关性探索，它不能发现某些条件子集下的局部结构。所以，传统聚类方法难以找出在某些样本下参与调控的基因聚类和与某些基因相关联的样本。

基因表达谱分析的双聚类方法可在数据矩阵的行和列两个方向上同时聚类，能够发现基因芯片数据中隐藏的大量有用的局部模式，获得基因的共调控证据。为寻找这些信息，Cheng和Church于2000年给出了双聚类的定义，双聚类的目的就是在基因表达数据矩阵中寻找满足条件的子矩阵，使得子矩阵中基因集在对应的条件集上表达波动一致，反之亦然。不同的双聚类算法采用不同的度量方式，所以能找到的双聚类有很大差别。近年来，双聚类局部模式的搜索算法主要是在以下四个方面进行扩展：1)通过传统聚类分别对矩阵的行和列进行聚类，然后合并这些中间结果得到最终双聚类，这类算法无法完全脱离聚类的全局性，不能很好地寻找局部模式。2)基于贪心迭代搜索方法寻找双聚类:将双聚类问题视为一个优化问题,运用贪心策略通过逐步删除可以使子矩阵的平均平方残基降低的行和列，得到一个最初的双聚类，此算法每次只能找到一个双聚类,且很难找到相互之间有重叠区的双聚类，结果的随机性很大。3)采用穷举策略寻找双聚类算法，多数为穷举小的子矩阵然后合并这些子矩阵的过程。这个过程常常借助图或者树的方式进行剪枝或者添加。双聚类分析实质上是一个NP-hard问题。所以，使用穷举策略的双聚类算法虽然能够找到较优的双聚类，算法的时间复杂度会随矩阵规模的增大而呈指数增长，其缺陷显而易见，很多算法对初始矩阵的大小有限制。4)将双聚类问题转化成数学或其他模型，应用各种方法在矩阵中寻找有规律的子矩阵。然而，一种数学模型只对应一种或少数的双聚类类型，算法实现较复杂，双聚类结果类型单一。双聚类问题的复杂度依赖于实际的问题构建和计算方法，特别是构建优化函数来评价给定双聚类的质量，几乎所有的计算方法进行双聚类其复杂度是NP完全的。协同免疫算法是一种进化计算技术，是模拟自然界生物免疫进化过程和机制求解实际问题的一类自组织、自适应和自学习的一种人工智能技术，双聚类算法结合协同免疫算法，可通过协同自适应搜索优化的方式，找到最优双聚类。

近年来核方法被用在传统的聚类分析中，核聚类的主要思想是首先通过一个非线性映射，将输入空间中的数据点映射到高维特征空间中，通过选取合适的Mercer核函数代替非线性映射的内积，在特征空间中进行聚类。核的聚类方法比经典的传统聚类方法有较大的改进，通过非线性映射增加了数据点线性可分的概率，即扩大数据类之间的差异，能较好地分辨、提取并放大有用的特征。当各类样本的边界是线性不可分以及类分布为非椭圆分布时，经典聚类算法失效的情况下，核聚类算法常常能得到更为准确的聚类。

发明内容

本发明通过把核化思想应用于双聚类，将一元核聚类思想拓展应用于二元核双聚类，通过提出核双聚类，应用核函数非线性映射，在高维空间的扩大数据类之间的差异，来更好地分辨、同时提取行列双向特征。提出并行协同免疫克隆Memetic核双聚类(CICMKB)实现算法，通过采用并行进化计算优化策略，在群体中反复应用选择、变异和交叉等操作，可同时搜索到既具有多样性又能达到全局最优的多个双聚类，用于解决共表达基因选取的难题。目前，基于核聚类算法思想的相关专利有：孙周宝和韩立新公开的一种基于流行学习的基因表达数据的聚类方法(CN:201110112132.X)。缑水平，焦李成等公开了一种基于NJW谱聚类标记的图像分割方法，然而，上述所有的研究都是针对传统的一元聚类算法--核方法的拓展。

本发明的目的是为解决生物样本经过基因表达芯片技术处理后，智能获取共表达基因的筛选难题。通过基因芯片组合单片机，应用单片机中的基因分析模块，来获取生物样本基因表达谱中共表达基因。

本发明的一种智能协同表达基因分析仪，包括基因芯片和单片机，所述单片机有四个模块组成：一种智能协同表达基因分析仪，包括基因芯片和单片机，其特征是所述单片机有四个模块组成：

(a)基因芯片表达谱读取模块，用于获取基因表达谱数据，其中包括生物细胞数据、生物组织数据、生物血样数据；

(b)协同表达基因提取模块，用于基因表达谱分析，提取协同表达基因；

(c)存储模块，对协同表达基因集的分析结果进行保存；

(d)输出模块，用于分析的结果输出；

所述协同表达基因提取模块的提取协同表达基因过程包含核函数选择、协同免疫克隆Memetic核双聚类算法、获取有重叠的双聚类和协同表达的基因集四部分，具体依次包括以下步骤：

(1)将基因芯片表达谱读取模块读取的基因表达谱数据，组成基因表达谱数据集A，所述基因表达谱数据集A是一个二维表格数据集合，对数据集A应用9近邻法进行缺失值填补；

(2)首先从核函数库列表中选择核函数，默认为高斯核函数，其核参数为1；下面(2)(3)两部分是一个循环执行，直到输出最优核双聚类集合；

(3)应用协同免疫克隆Memetic核双聚类算法对基因表达谱数据集A获取K个双聚类；所述的σ-核双聚类定义如下：

设一个n×m二维关系表达实数矩阵A=X×Y={x_ij}(i∈[1,n],j∈[1,m])，其中X代表数据纪录行{x₁,x₂,...,x_n},Y为对应的属性列{y₁,y₂,...,y_m}，x_ij为表达数据矩阵A中的元素。若x_ij通过某映射函数φ(x)，被投射到高维特征空间F。在高维特征空间F中，设I、J分别为X、Y的子集，则对指定的子矩阵B=I×J具有以下总体核行方差RVAR，即为所有记录X的方差平均，和平均核平方残差MSR如下：

$\mathrm{RVAR}(I,J)=\frac{1}{\left|I\right|\left|J\right|}\underset{i\∈I,j\∈J}{\mathrm{\Σ}}K(x_{\mathrm{ij}},x_{\mathrm{ij}})+\frac{1}{{\left|I\right|\left|J\right|}^2}(\frac{1}{\left|J\right|}-2)\underset{i\∈I}{\mathrm{\Σ}}\left(\underset{j\∈J,v\∈J}{\mathrm{\Σ}}K(x_{\mathrm{ij}},x_{\mathrm{iv}})\right)$

$\mathrm{MSR}(I,J)=\mathrm{RVAR}(I,J)+\frac{2}{{\left|I\right|}^2\left|J\right|}\underset{j\∈J}{\mathrm{\Σ}}\left(\underset{j\∈J,u\∈I}{\mathrm{\Σ}}K\left(x_{\mathrm{ij}}{,x}_{\mathrm{uj}}\right)\right)$

$-\frac{2}{{\left|I\right|\left|J\right|}^2}\underset{i\∈I,j\∈J}{\mathrm{\Σ}}\underset{u\∈I,v\∈J}{\mathrm{\Σ}}K(x_{\mathrm{ij}},x_{\mathrm{uv}})$

对于均方残差阈值σ≥0，如果子矩阵B_IJ满足MSR(I,J)≤σ，则称B_IJ子矩阵为一个σ-核双聚类，这里K(x_ij,x_uv)为核函数，通过选择不同的核函数，对 $\left\{\begin{array}{c}\min f\left(x\right)=\frac{\mathrm{MSR}(I,J)}{\mathrm{\&sigma;gRVAR}(I,J)}+\frac{\left|X\right|\left|Y\right|}{\left|I\right|\left|J\right|}\\ \mathrm{sudject\; to\; MSR}(I,J)<\mathrm{\&sigma;},\mathrm{\&sigma;}\&GreaterEqual;0\end{array}\right.$ 来求解优化问题，找到各种可重叠双聚类；

其中，|X|、|Y|分别为表达实数矩阵A的行数、列数；|I|、|J|分别表示所求双聚类的实数矩阵X、Y子集的行数和列数；

(3.1)初始抗体群：在一个L维空间中，初始化生成K个差异种群G_i，与K个双聚类对应；每个子种群G_i有100个初始抗体{z_i1,z_i2,...,z_ij}，i=1,2...,K，j=1,2,...100；定义每个抗体z_ij的编码为二进制编码，其长度为L=m+n，这里n和m分别为二维数据整体行的个数和列的个数；如果抗体z_ij某一位置为l，则意味着相应的行或者列包括在K个双聚类中；初始化设置迭代终止均方残差MSR阈值σ。

(3.2)设计抗体适应度函数affinity(z_ij)，计算每个子种群G_i中每个抗体z_ij的亲和度值，这里又称适应度值(i=1,2...,K,j=1,2,...100)：

$\mathrm{affinity}\left(z_{\mathrm{ij}}\right)=\frac{\mathrm{MSR}(I,J)}{\mathrm{\&sigma;gRVAR}(I,J)}+\frac{\left|X\right|\left|Y\right|}{\left|I\right|\left|J\right|}$

其中，I,J分别为动态所求的子矩阵行个数和列个数，RVAR为动态核子矩阵总体行方差和MSR为动态核子矩阵平均平方残差。

(3.3)抗体促进和抑制：各个子种群G_i并行计算，计算子种群中每个抗体z_ij的抗体浓度：

$C\left(z_{\mathrm{ij}}\right)=\frac{1}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}|\mathrm{affinity}\left(z_{\mathrm{ij}}\right)-\mathrm{affinity}\left(z_{\mathrm{ik}}\right)|}$

其中子种群G_i中抗体浓度C(z_ij)与适应度affinity(z_ij)相关，N为当前子种群G_i中的抗体个数；

在各子种群内，根据抗体的浓度调节机制，由抗体的期望繁殖率p_j对抗体有选择地进行克隆复制；

设在当前子种群G_i中，抗体z_ij的期望繁殖率p_j是由适应度概率p_fj和浓度抑制概率p_dj两部分组成：

p_j=μ·p_fj+(1-μ)P_dj

其中，μ是常数调节因子，取值为0.6；

适应度概率p_fj为： $p_f\left(z_{\mathrm{ij}}\right)=\frac{\mathrm{affinity}\left(z_{\mathrm{ij}}\right)}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{affinity}\left(z_{\mathrm{ik}}\right)};$

抗体浓度抑制概率p_dj为： $p_d\left(z_{\mathrm{ij}}\right)=\frac{1}{C\left(z_{\mathrm{ij}}\right)\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{C\left(z_{\mathrm{ik}}\right)}},$ 抗体的浓度高则减小该个体的选择概率，反之则增加该个体的选择概率；

(3.4)免疫选择克隆就是将免疫选择算子选中的抗体个体进行复制，于是被选中的抗体z_ij的克隆复制个数为：

$N_c=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{round}(v\&CenterDot;p_j)$

其中，round表示取整；v>1是参数因子，这里取值为当前种群个数的2倍；p_j是期望繁殖率；

(3.5)子种群内抗体间交叉，新抗体产生；每个抗体被选择交叉的概率为P_c；

P_c=P₁·r^t/T

这里P₁是一个预先设定的概率，r∈[0,1]是一个变异因子常数，在这里可设为一个[0,1]的随机数，它起着调整交叉概率；t为当前演化代数，T为最大可迭代代数；依据概率P_c从当代种群中选取若干个体，按照交叉算子

$\left\{\begin{array}{c}{G}_1^{\mathrm{new}}={\mathrm{\&omega;}}_1\&CenterDot;G_1+(1-{\mathrm{\&omega;}}_1)\&CenterDot;G_2\\ G_2^{\mathrm{new}}={\mathrm{\&omega;}}_2\&CenterDot;G_2+(1-{\mathrm{\&omega;}}_2)\&CenterDot;G_1\end{array}\right.$

进行变异，其中G₁,G₂为从种群中随机选择的两个父个体，为通过交叉运算子运算后产生的子代对应新个体；ω₁，ω₂为[0,1]上随机选取的参数；

(3.6)在各子种群内，以突变概率P_m，执行抗体间的变异，其中

$P_m=P_2{e}^{-r_{2}t/T}$

这里，P₂代表了P_m的初始值，P₂设置为0.5；r₂为抗体突变参数，t为当前演化代数，T为最大可迭代代数；依据概率P_m从中选取若干个体，按照变异算子

$V^{\&prime;}=\left\{\begin{array}{c}\mathrm{int}(V+(b_{\sup }-V){[e^{-r_2}\&CenterDot;(1-\frac{t}{T})]}^2),\mathrm{sign}=0\\ \mathrm{int}(V-(V-b_{\inf }){[r^{-r_2}\&CenterDot;(1-\frac{t}{T})]}^2),\mathrm{sign}=1\end{array}\right.$

进行变异；

其中，V'为变异后的参数；V是选中的变异参数，这里V=0.4；sign随机取0或者1；b_sup=1和b_inf=0分别为参数取值的上界和下界；int()为取整；

(3.7)亲和度计算：重新计算子种群内每个抗体适应度affinity(z_ij)；

(3.8)最优抗体：将不同子种群中的适应度函数affinity(z_ij)最小的抗体作为局部最优抗体保留；

(3.9)抗体记忆池：子种群间筛选的最优抗体，存储入抗体记忆池中，形成精英抗体群；

(3.10)对精英抗体群中每个抗体，采用单纯形法进行Memetic局部搜索，获取每个抗体的局部最优值，这些最优值更新精英抗体记忆池；所述局部最优值是指抗体的适应度函数affinity(z_ij)稳定于某个值，且适应度函数affinity(z_ij)小于种群中的局部最优个体；

(3.11)当算法已达到最大进化代数或者前一次迭代与当前迭代精英抗体亲和度整体均值的差小于常数eps=0.0001且MSR(I,J)≤σ，算法终止，输出前K个精英抗体，为最优K个双聚类；否则，算法转向下一步；

(3.12)各个子种群群体更新，每个子种群重新恢复POPSIZE个种群大小，算法转向步骤(3.2)，重新循环执行，直到输出当前核函数下最优K个双聚类，并存入存储器当前核函数下的最优K个双聚类的每个平均核平方残差MSR(I,J)，它的总体核行方差RVAR(I,J)，和它的行数I与列数J。

算法转向步骤(2)，选择另外一个核函数，从(2)-(3)重新计算另外一个核函数下的最优K个双聚类，并存入存储器另外一个核函数下的最优K个双聚类的每个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；直到所有的核函数全被选择运行完；

(4)根据最优K个双聚类获得K个协同表达基因集；依据存储器中每个核函数下的最优K个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；在多个核函数间，选择具有最低平均均方残差且双聚类行列大小大的核函数，所对应的K个双聚类为输出的K个双聚类，其选择计算标准为：

$\mathrm{sum}\_\mathrm{fitness}=\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\frac{\mathrm{MS}{R}_i(I,J)}{\mathrm{\&sigma;gRVA}{R}_i(I,J)}+\frac{K}{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}(I_i+J_i)}$

其中，i=1,2,...,K表示在某个核函数下的K个双聚类；选择最小的sum_fitness所对应的核函数下的K个双聚类为输出的K个双聚类，其对应K个协同表达基因集，包括：1)相似表达的基因集，共同高表达或共同低表达，变化几乎相同基因集；2)表达趋势一致的基因集，共同高表达或共同低表达，变化一致基因集；3)反向表达相关的基因集，基因表达总相反；输出协同表达的基因集。

作为优选的技术方案：

如上所述的一种智能协同表达基因分析仪，所述的另外一个核函数是对称正定的距离测度函数，为线性核函数、Tricube函数和Epanechnikov函数。

如上所述的一种智能协同表达基因分析仪，所述的高斯核函数为其核参数δ设为1，x_ij,x_uv为基因表达谱数据集A中的两个元素，K(x_ij,x_uv)为核函数，exp表示e指数函数。

本发明不仅有方法创新而且具有实用价值，通过引进了单片机嵌入式系统分析手段，让方法架构与单片机上，使得应用更加直观化，用户可直接通过后台的显示器便能查看基因表达谱芯片中基因的共表达信息，结果一目了然。本发明利用基因芯片采集技术获取生物样本的基因表达谱，应用单片机嵌入式分析技术获取协同表达的基因集。

有益效果

本发明的一种智能协同表达基因分析仪，具有以下优点：

(1)在高维特征空间中获取双聚类，增加二维行列数据点可分的概率，扩大数据类之间的差异，能更好地分辨、提取并放大有用的特征，同时解决传统聚类随着维数增长复杂度迅速上升而导致聚类算法的性能下降的问题，这在一般双聚类模型中无法获得。

(2)本发明提出非线性核双聚类的定义和协同免疫克隆Memetic核双聚类的获取，可在数据矩阵的行和列两个方向上同时聚类，采用并行进化计算优化策略，可同时搜索到既具有多样性又能达到全局最优的多个双聚类，获取各种协同表达基因集，一般算法不具有这个特性。

(3)基因分析模块中核双聚类可根据不同的核函数的选择，演变成各种双聚类模型，包括线性和非线性双聚类，这在其他双聚类方法中无法获得。

(4)本发明分析仪通过软硬件结合的手段，将共表达基因直接以可视化图形化方法展现，使得应用更加直观化，让使用者可直接通过后台的显示器观察基因表达谱芯片中基因的共表达信息，结果一目了然，当前还未有此类硬件。

附图说明

图1描述了智能协同表达基因分析仪基本技术流程

图2描述了基因分析模块中协同免疫克隆Memetic核双聚类算法的流程

具体实施方式

下面结合具体实施方式，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

下面结合啤酒酵母细胞周期表达值基因芯片，获取数据中的局部模式共表达基因调控信息，来进一步阐述本发明图1中智能协同表达基因分析仪的基因分析模块，其基本技术流程如下：

(1)生物样本：收集啤酒酵母细胞不同时期的细胞生物样本。

(2)基因芯片：提取啤酒酵母细胞不同时期的细胞生物样本的mRNA，拟转录为cDNA，并且标记，将其点阵到基因芯片上。

(3)获取基因表达谱：利用基因芯片技术，通过杂交、激光扫描，获取基因芯片中探针所代表的基因的荧光强度，将其转换成啤酒酵母细胞周期表达的基本数据集(rawdata)，通过移除重复出现的探针数据；移除标记数据；并以中位数对荧光强度的数据进行标准化(Normalized)校正，转换基因芯片中基本数据集，为啤酒酵母细胞周期的基因表达谱，这里包括2884个基因在17个条件下的基因表达水平值。(2)-(3)过程是基因芯片技术过程，不是本发明的重点。

(4)在单片机中，基因芯片表达谱读取模块：通过计算机程序，获取啤酒酵母细胞周期的基因表达谱数据的二维表达实数矩阵A，数据矩阵A中有34个丢失值，这里采用9近邻法进行缺失值填补。

(5)在单片机中，协同表达基因提取模块的提取协同表达基因过程包含核函数选择、协同免疫克隆Memetic核双聚类算法、获取有重叠的双聚类和协同表达的基因集四部分。

(5.1)核函数选择:首先从核函数库列表中选择核函数，首选高斯核函数，其核参数为1；核函数库中是常用的核函数可以选择对称正定的核函数为距离测度,如线性核函数、

首选选择高斯核函数 $K\left(x_{\mathrm{ij},}{x}_{\mathrm{uv}}\right)=\exp (-\frac{{\left|\right|x_{\mathrm{ij}}-x_{\mathrm{uv}}\left|\right|}^2}{{\mathrm{\&delta;}}^2}),$ 其核参数δ设为1。

(5.2)应用协同免疫克隆Memetic核双聚类算法对基因表达谱数据集A获取当前核函数下的K个双聚类；

设置协同免疫克隆Memetic算法的初始参数，如表2所示。输入迭代终止均方残差阈值σ=300和最大双聚类个数K=15，依据图2描述的协同免疫克隆Memetic核双聚类算法，获取在均方残差阈值σ下的15个最大双聚类，对应的15个协同表达基因集。

表2协同免疫克隆Memetic算法的参数设置

(5.3)获取有重叠的双聚类:算法转向步骤(5.1)，选择另外一个核函数，例如：选择Tricube核函数返回步骤(5.2)，重新计算Tricube核函数下最优K个双聚类，并存入存储器这最优K个双聚类的每个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；

(5.1)-(5.3)实际是一个循环，直到分别获取所有不同核函数下的核双聚类集合；

(5.4)协同表达基因提取:根据不同核函数下核双聚类集合，获得最优K个协同表达基因集；依据存储器中每个核函数下的最优K个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；在多个核函数间，选择具有最低平均均方残差且双聚类行列大小大的核函数，所对应的K个双聚类为输出的K个双聚类，其选择计算标准为：

$\mathrm{sum}\_\mathrm{fitness}=\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\frac{\mathrm{MS}{R}_i(I,J)}{\mathrm{\&sigma;gRVA}{R}_i(I,J)}+\frac{K}{\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}(I_i+J_i)}$

其中，i=1,2,...,K表示在某个核函数下的K个双聚类；选择最小的sum_fitness所对应的核函数下的K个双聚类为输出的K个双聚类，其对应K个协同表达基因集，包括：1)相似表达的基因集，共同高表达或共同低表达，变化几乎相同基因集；2)表达趋势一致的基因集，共同高表达或共同低表达，变化一致基因集；3)反向表达相关的基因集，基因表达总相反；输出协同表达的基因集。

(6)存储模块，对最优K=15个协同表达基因集的结果进行保存；

(7)输出这K=15个最大双聚类，即为酵母数据集所挖掘出的在不同之间段中共表达的基因集。这些随着不同时间协同表达调控基因含义可在GeneOntology水平上详细解释。

图2详细描述了基因分析模块中协同免疫克隆Memetic核双聚类算法流程。

协同免疫克隆Memetic核双聚类算法是根据变异和突变原理，并行产生多个数据集。以每个数据点积开始搜索，但由于拟牛顿梯度下降法（爬山法）对初始点非常敏感，往往陷入局部极值,采用进化计算非梯度下降手段使解“跳”到爬山法可求解范围。在爬山法可求解范围应用单纯形法进行局部搜索，找到局部极值。汇聚局部最优极值点，构成精英机制，寻找全局最优极值（全局极值可能存在多个）。

协同免疫克隆Memetic算法采用协同双层进化机制，上层协作优化空间实现进化协作过程；底层种群空间实现免疫克隆Memetic算法，并根据进化代数自适应调节其参数。在协同免疫克隆Memetic核双聚类算法中，首先定义抗原，抗原识别，抗原和抗体之间的亲和性。抗原对应于算法所要解决的问题，这里是指待优化的目标函数式。抗原识别对应于算法所要解决问题的可行解，亲和度函数是评价函数，反应抗原和抗体之间的亲和性，是解与目标函数的匹配程度。

图2中描述的协同免疫克隆Memetic核双聚类算法具体步骤如下：

(1)初始抗体群：在一个L维空间中，初始化生成K个差异种群G_i，与K个双聚类对应；每个子种群G_i有100个初始抗体{z_i1,z_i2,...,z_ij}，i=1,2...,K，j=1,2,...100；定义每个抗体z_ij的编码为二进制编码，其长度为L=m+n，这里n和m分别为二维数据整体行的个数和列的个数；如果抗体z_ij某一位置为l，则意味着相应的行或者列包括在K个双聚类中；初始化设置迭代终止均方残差MSR阈值σ。

(2)设计抗体适应度函数affinity(z_ij)，计算每个子种群G_i中每个抗体z_ij的亲和度值，这里又称适应度值(i=1,2...,K,j=1,2,...100)：

$\mathrm{affinity}\left(z_{\mathrm{ij}}\right)\frac{\mathrm{MSR}(I,J)}{\mathrm{\&sigma;gRVAR}(I,J)}\frac{\left|X\right|\left|Y\right|}{\left|I\right|\left|J\right|}$

其中，I,J分别为动态所求的子矩阵行个数和列个数，RVAR为动态核子矩阵总体行方差和MSR为动态核子矩阵平均平方残差。

(3)抗体促进和抑制：各个子种群G_i并行计算，计算子种群中每个抗体z_ij的抗体浓度：

$C\left(z_{\mathrm{ij}}\right)=\frac{1}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}|\mathrm{affinity}\left(z_{\mathrm{ij}}\right)-\mathrm{affinity}\left(z_{\mathrm{ik}}\right)|}$

其中子种群G_i中抗体浓度C(z_ij)与适应度affinity(z_ij)相关，N为当前子种群G_i中的抗体个数；

在各子种群内，根据抗体的浓度调节机制，由抗体的期望繁殖率p_j对抗体有选择地进行克隆复制；

设在当前子种群G_i中，抗体z_ij的期望繁殖率p_j是由适应度概率p_fj和浓度抑制概率p_dj两部分组成：

p_j=μ·p_fj+(1-μ)P_dj

其中，μ是常数调节因子，取值为0.6；

适应度概率p_fj为： $p_f\left(z_{\mathrm{ij}}\right)=\frac{\mathrm{affinity}\left(z_{\mathrm{ij}}\right)}{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{affinity}\left(z_{\mathrm{ik}}\right)};$

抗体浓度抑制概率p_dj为：抗体的浓度高则减小该个体的选择概率，反之则增加该个体的选择概率；

(4)免疫选择克隆就是将免疫选择算子选中的抗体个体进行复制，于是被选中的抗体z_ij的克隆复制个数为：

$N_c=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{round}(v\&CenterDot;p_j)$

其中，round表示取整；v>1是参数因子，这里取值为当前种群个数的2倍；p_j是期望繁殖率；

(5)子种群内抗体间交叉，新抗体产生；每个抗体被选择交叉的概率为P_c；

P_c=P₁·r^t/T

这里P₁是一个预先设定的概率，r∈[0,1]是一个变异因子常数，在这里可设为一个[0,1]的随机数，它起着调整交叉概率；t为当前演化代数，T为最大可迭代代数；依据概率P_c从当代种群中选取若干个体，按照交叉算子

$\left\{\begin{array}{c}{G}_1^{\mathrm{new}}={\mathrm{\&omega;}}_1\&CenterDot;G_1+(1-{\mathrm{\&omega;}}_1)\&CenterDot;G_2\\ G_2^{\mathrm{new}}={\mathrm{\&omega;}}_2\&CenterDot;G_2+(1-{\mathrm{\&omega;}}_2)\&CenterDot;G_1\end{array}\right.$

进行变异，其中G₁,G₂为从种群中随机选择的两个父个体，为通过交叉运算子运算后产生的子代对应新个体；ω₁，ω₂为[0,1]上随机选取的参数；

(6)在各子种群内，以突变概率P_m，执行抗体间的变异，其中

$P_m=P_2{e}^{-r_{2}t/T}$

这里，P₂代表了P_m的初始值，P₂设置为0.5；r₂为抗体突变参数，t为当前演化代数，T为最大可迭代代数；依据概率P_m从中选取若干个体，按照变异算子

$V^{\&prime;}=\left\{\begin{array}{c}\mathrm{int}(V+(b_{\sup }-V){[e^{-r_2}\&CenterDot;(1-\frac{t}{T})]}^2),\mathrm{sign}=0\\ \mathrm{int}(V-(V-b_{\inf }){[r^{-r_2}\&CenterDot;(1-\frac{t}{T})]}^2),\mathrm{sign}=1\end{array}\right.$

进行变异；

其中，V'为变异后的参数；V是选中的变异参数，这里V=0.4；sign随机取0或者1；b_sup=1和b_inf=0分别为参数取值的上界和下界；int()为取整；

(7)亲和度计算：重新计算子种群内每个抗体适应度affinity(z_ij)；

(8)最优抗体：将不同子种群中的适应度函数affinity(z_ij)最小的抗体作为局部最优抗体保留；

(9)抗体记忆池：子种群间筛选的最优抗体，存储入抗体记忆池中，形成精英抗体群；

(10)对精英抗体群中每个抗体，采用单纯形法进行Memetic局部搜索，获取每个抗体的局部最优值，这些最优值更新精英抗体记忆池；所述局部最优值是指抗体的适应度函数affinity(z_ij)稳定于某个值，且适应度函数affinity(z_ij)小于种群中的局部最优个体；

(11)当算法已达到最大进化代数或者前一次迭代与当前迭代精英抗体亲和度整体均值的差小于常数eps=0.0001且MSR(I,J)≤σ，算法终止，输出前K个精英抗体，为最优K个双聚类；否则，算法转向下一步；

(12)各个子种群群体更新，每个子种群重新恢复POPSIZE个种群大小，算法转向步骤(3.2)，重新循环执行，直到输出当前核函数下最优K个双聚类，并存入存储器当前核函数下的最优K个双聚类的每个平均核平方残差MSR(I,J)，它的总体核行方差RVAR(I,J)，和它的行数I与列数J。

Claims (4)

1.一种智能协同表达基因分析仪，包括基因芯片和单片机，其特征是所述单片机有四个模块组成：

(a)基因芯片表达谱读取模块，用于获取基因表达谱数据；

(b)协同表达基因提取模块，用于基因表达谱分析，提取协同表达基因；

(c)存储模块，对协同表达基因集的分析结果进行保存；

(d)输出模块，用于分析的结果输出；

所述协同表达基因提取模块的提取协同表达基因过程包含核函数选择、协同免疫克隆Memetic核双聚类算法、获取有重叠的双聚类和协同表达的基因集四部分，具体依次包括以下步骤：

(1)将基因芯片表达谱读取模块读取的基因表达谱数据，组成基因表达谱数据集A，所述基因表达谱数据集A是一个二维表格数据集合，对数据集A应用9近邻法进行缺失值填补；

(2)首先从核函数库列表中选择核函数，默认为高斯核函数，其核参数为1；下面(2)(3)两部分是一个循环执行，直到输出最优核双聚类集合；

(3)应用协同免疫克隆Memetic的σ-核双聚类算法对基因表达谱数据集A获取K个双聚类；所述的σ-核双聚类定义如下：

设一个n×m二维关系表达实数矩阵A＝X×Y＝{x_ij}(i∈[1,n],j∈[1,m])，其中X代表数据纪录行{x₁,x₂,...,x_n},Y为对应的属性列{y₁,y₂,...,y_m}，x_ij为表达数据矩阵A中的元素；若x_ij通过某映射函数φ(x)，被投射到高维特征空间F；在高维特征空间F中，设I、J分别为X、Y的子集，则对指定的子矩阵B＝I×J具有以下总体核行方差RVAR即为所有记录X的方差平均，和平均核平方残差MSR如下：

$RVAP(I,J)=\frac{1}{\left|I\right|\left|J\right|}\underset{i\&Element;I,j\&Element;J}{\mathrm{\&Sigma;}}K(x_{ij},x_{ij})+\frac{1}{\left|I\right||J{|}^2}(\frac{1}{\left|J\right|}-2)\underset{i\&Element;I}{\mathrm{\&Sigma;}}\left(\underset{j\&Element;J,v\&Element;J}{\mathrm{\&Sigma;}}K\right(x_{ij},x_{iv}\left)\right)$

$\begin{array}{c}MSR(I,J)=RVAR(I,J)+\frac{1}{\left|I{|}^2\right|J|}\underset{j\&Element;J}{\mathrm{\&Sigma;}}\left(\underset{i\&Element;I,u\&Element;I}{\mathrm{\&Sigma;}}K\right(x_{ij},x_{uj}\left)\right)\\ -\frac{1}{\left|I\right||J{|}^2}\underset{i\&Element;I,j\&Element;J}{\mathrm{\&Sigma;}}\underset{u\&Element;I,u\&Element;J}{\mathrm{\&Sigma;}}K(x_{ij},x_{uv})\end{array}$

对于均方残差阈值σ≥0，如果子矩阵B_IJ满足MSR(I,J)≤σ，则称B_IJ子矩阵为一个σ-核双聚类，这里K(x_ij,x_uv)为核函数，通过选择不同的核函数，对

来求解优化问题，找到各种可重叠双聚类；

其中，|X|、|Y|分别为表达实数矩阵A的行数、列数；|I|、|J|分别表示所求双聚类的实数矩阵X、Y子集的行数和列数；

(3.1)初始抗体群：在一个L维空间中，初始化生成K个差异种群G_i，与K个双聚类对应；每个子种群G_i有100个初始抗体{z_i1,z_i2,...,z_ij}，i＝1,2...,K，j＝1,2,...100；定义每个抗体z_ij的编码为二进制编码，其长度为L＝m+n，这里n和m分别为二维数据整体行的个数和列的个数；如果抗体z_ij某一位置为l，则意味着相应的行或者列包括在K个双聚类中；初始化设置迭代终止均方残差MSR阈值σ；

(3.2)设计抗体适应度函数affinity(z_ij)，计算每个子种群G_i中每个抗体z_ij的亲和度值，这里又称适应度值(i＝1,2...,K,j＝1,2,...100)：

其中，I,J分别为动态所求的子矩阵行个数和列个数，RVAR为动态核子矩阵总体行方差和MSR为动态核子矩阵平均平方残差；

(3.3)抗体促进和抑制：各个子种群G_i并行计算，计算子种群中每个抗体z_ij的抗体浓度：

$C\left(z_{ij}\right)=\frac{1}{\underset{k=1}{\overset{N}{\&Sigma;}}|affinity\left(z_{ij}\right)-affinity\left(z_{ik}\right)|}$

其中子种群G_i中抗体浓度C(z_ij)与适应度affinity(z_ij)相关，N为当前子种群G_i中的抗体个数；

在各子种群内，根据抗体的浓度调节机制，由抗体的期望繁殖率p_j对抗体有选择地进行克隆复制；

设在当前子种群G_i中，抗体z_ij的期望繁殖率p_j是由适应度概率p_fj和浓度抑制概率p_dj两部分组成：

p_j＝μ·p_fj+(1-μ)P_dj

其中，μ是常数调节因子，取值为0.6；

适应度概率p_fj为： $p_f\left(z_{ij}\right)=\frac{affinity\left(z_{ij}\right)}{\underset{k=1}{\overset{N}{\&Sigma;}}affinity\left(z_{ik}\right)};$

抗体浓度抑制概率p_dj为：抗体的浓度高则减小该抗体的选择概率，反之则增加该抗体的选择概率；

(3.4)免疫选择克隆就是将免疫选择算子选中的抗体个体进行复制，于是被选中的抗体z_ij的克隆复制个数为：

$N_c=\underset{i=1}{\overset{n}{\&Sigma;}}round(v\&CenterDot;p_j)$

其中，round表示取整；v＞1是参数因子，这里取值为当前种群个数的2倍；p_j是期望繁殖率；

(3.5)子种群内抗体间交叉，新抗体产生；每个抗体被选择交叉的概率为P_c；

P_c＝P₁·r^t/T

这里P₁是一个预先设定的概率，r∈[0,1]是一个变异因子常数，在这里可设为一个[0,1]的随机数，它起着调整交叉概率；t为当前演化代数，T为最大可迭代代数；依据概率P_c从当代种群中选取若干个体，按照交叉算子

$\left\{\begin{array}{c}{G}_1^{new}={\mathrm{\&omega;}}_1\&CenterDot;G_1+(1-{\mathrm{\&omega;}}_1)\&CenterDot;G_2\\ G_2^{new}={\mathrm{\&omega;}}_2\&CenterDot;G_2+(1-{\mathrm{\&omega;}}_2)\&CenterDot;G_1\end{array}\right.$

进行变异，其中G₁,G₂为从种群中随机选择的两个父个体，为通过交叉运算子运算后产生的子代对应新个体；ω₁，ω₂为[0,1]上随机选取的参数；

(3.6)在各子种群内，以突变概率P_m，执行抗体间的变异，其中

$P_m=P_2{e}^{-r_{2}t/T}$

这里，P₂代表了P_m的初始值，P₂设置为0.5；r₂为抗体突变参数，t为当前演化代数，T为最大可迭代代数；依据概率P_m从中选取若干个体，按照变异算子

$V^{\&prime;}=\left\{\begin{array}{c}\mathrm{int}(V+(b_{\sup }-V){\&lsqb;e^{-r_2}\&CenterDot;(1-\frac{t}{T})\&rsqb;}^2),sign=0\\ \mathrm{int}(V-(V-b_{\inf }){\&lsqb;e^{-r_2}\&CenterDot;(1-\frac{t}{T})\&rsqb;}^2),sign=1\end{array}\right.$

进行变异；

其中，V'为变异后的参数；V是选中的变异参数，这里V＝0.4；sign随机取0或者1；

b_sup＝1和b_inf＝0分别为参数取值的上界和下界；int()为取整；

(3.7)亲和度计算：重新计算子种群内每个抗体适应度affinity(z_ij)；

(3.8)最优抗体：将不同子种群中的适应度函数affinity(z_ij)最小的抗体作为局部最优抗体保留；

(3.9)抗体记忆池：子种群间筛选的最优抗体，存储入抗体记忆池中，形成精英抗体群；

(3.10)对精英抗体群中每个抗体，采用单纯形法进行Memetic局部搜索，获取每个抗体的局部最优值，这些最优值更新精英抗体记忆池；所述局部最优值是指抗体的适应度函数affinity(z_ij)稳定于某个值，且适应度函数affinity(z_ij)小于种群中的局部最优个体；

(3.11)当算法已达到最大进化代数或者前一次迭代与当前迭代精英抗体亲和度整体均值的差小于常数eps＝0.0001且MSR(I,J)≤σ，算法终止，输出前K个精英抗体，为最优K个双聚类；否则，算法转向下一步；

(3.12)各个子种群群体更新，每个子种群重新恢复POPSIZE个种群大小，算法转向步骤(3.2)，重新循环执行，直到输出当前核函数下最优K个双聚类，并存入存储器当前核函数下的最优K个双聚类的每个平均核平方残差MSR(I,J)，它的总体核行方差RVAR(I,J)，和它的行数I与列数J；

算法转向步骤(2)，选择另外一个核函数，从(2)-(3)重新计算另外一个核函数下的最优K个双聚类，并存入存储器另外一个核函数下的最优K个双聚类的每个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；直到所有的核函数全被选择运行完；

(4)根据最优K个双聚类获得K个协同表达基因集；依据存储器中每个核函数下的最优K个双聚类的平均行核平方残差MSR(I,J)、总体核行方差RVAR(I,J)和它的行数I与列数J；在多个核函数间，选择具有最低平均均方残差且双聚类行列大小大的核函数，所对应的K个双聚类为输出的K个双聚类，其选择计算标准为：

其中，i＝1,2,...,K表示在某个核函数下的K个双聚类；选择最小的sum_fitness所对应的核函数下的K个双聚类为输出的K个双聚类，其对应K个协同表达基因集，包括：1)相似表达的基因集，共同高表达或共同低表达，变化几乎相同基因集；2)表达趋势一致的基因集，共同高表达或共同低表达，变化一致基因集；3)反向表达相关的基因集，基因表达总相反；输出协同表达的基因集。

2.根据权利要求1所述的一种智能协同表达基因分析仪，其特征在于，所述的另外一个核函数是对称正定的距离测度函数，为线性核函数、Tricube函数和Epanechnikov函数。

3.根据权利要求1所述的一种智能协同表达基因分析仪，其特征在于，所述的高斯核函数为其核参数δ设为1，x_ij,x_uv为基因表达谱数据集A中的两个元素，K(x_ij,x_uv)为核函数，exp表示e指数函数。

4.根据权利要求1所述的一种智能协同表达基因分析仪，其特征在于，所述的基因芯片可采集的数据包括生物细胞数据、生物组织数据、生物血样数据。