发明名称 一种智能协同表达基因分析仪
摘要 本发明公开了一种智能协同表达基因分析仪,包括基因芯片和单片机,利用基因芯片采集技术获取生物样本的基因表达谱,应用单片机嵌入式分析技术获取协同表达的基因集。单片机有四个模块组成:基因芯片表达谱读取模块;协同表达基因提取模块;存储模块;输出模块;协同表达基因提取模块的提取协同表达基因过程包含核函数选择、协同免疫克隆Memetic核双聚类算法、获取有重叠的双聚类和协同表达的基因集四部分。这种智能协同表达基因分析仪可从生物样本的基因表达谱中提取/表达趋势一致的基因集和反向表达相关的基因集。通过对这些共同表达基因的寻找,不仅可以对基因的功能研究给予提示,还可以对基因调控途径和调控网络的研究给予启发。
申请公布号 CN103164631A 申请公布日期 2013.06.19
申请号 CN201310130664.5 申请日期 2013.04.16
申请人 东华大学 发明人 丁永生;程丽俊;程铎辉
分类号 G06F19/10(2011.01)I 主分类号 G06F19/10(2011.01)I
代理机构 上海天翔知识产权代理有限公司 31224 代理人 吕伴
主权项 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<sub>ij</sub>}(i∈[1,n],j∈[1,m]),其中X代表数据纪录行{x<sub>1</sub>,x<sub>2</sub>,...,x<sub>n</sub>},Y为对应的属性列{y<sub>1</sub>,y<sub>2</sub>,...,y<sub>m</sub>},x<sub>ij</sub>为表达数据矩阵A中的元素;若xij通过某映射函数φ(x),被投射到高维特征空间F;在高维特征空间F中,设I、J分别为X、Y的子集,则对指定的子矩阵B=I×J具有以下总体核行方差RVAR即为所有记录X的方差平均,和平均核平方残差MSR如下:<maths num="0001"><![CDATA[<math><mrow><mi>RVAR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><mo>|</mo><mi>I</mi><mo>|</mo><mo>|</mo><mi>J</mi><mo>|</mo></mrow></mfrac><munder><mi>&Sigma;</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>I</mi><mo>,</mo><mi>j</mi><mo>&Element;</mo><mi>J</mi></mrow></munder><mi>K</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>ij</mi></msub><mo>,</mo><msub><mi>x</mi><mi>ij</mi></msub><mo>)</mo></mrow><mo>+</mo><mfrac><mn>1</mn><msup><mrow><mo>|</mo><mi>I</mi><mo>|</mo><mo>|</mo><mi>J</mi><mo>|</mo></mrow><mn>2</mn></msup></mfrac><mrow><mo>(</mo><mfrac><mn>1</mn><mrow><mo>|</mo><mi>J</mi><mo>|</mo></mrow></mfrac><mo>-</mo><mn>2</mn><mo>)</mo></mrow><munder><mi>&Sigma;</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>I</mi></mrow></munder><mrow><mo>(</mo><munder><mi>&Sigma;</mi><mrow><mi>j</mi><mo>&Element;</mo><mi>J</mi><mo>,</mo><mi>v</mi><mo>&Element;</mo><mi>J</mi></mrow></munder><mi>K</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>ij</mi></msub><mo>,</mo><msub><mi>x</mi><mi>iv</mi></msub><mo>)</mo></mrow><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0002"><![CDATA[<math><mrow><mi>MSR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow><mo>=</mo><mi>RVAR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow><mo>+</mo><mfrac><mn>2</mn><mrow><msup><mrow><mo>|</mo><mi>I</mi><mo>|</mo></mrow><mn>2</mn></msup><mo>|</mo><mi>J</mi><mo>|</mo></mrow></mfrac><munder><mi>&Sigma;</mi><mrow><mi>j</mi><mo>&Element;</mo><mi>J</mi></mrow></munder><mrow><mo>(</mo><munder><mi>&Sigma;</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>I</mi><mo>,</mo><mi>u</mi><mo>&Element;</mo><mi>I</mi></mrow></munder><mi>K</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>ij</mi></msub><mo>,</mo><msub><mi>x</mi><mi>uj</mi></msub><mo>)</mo></mrow><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0003"><![CDATA[<math><mrow><mo>-</mo><mfrac><mn>2</mn><msup><mrow><mo>|</mo><mi>I</mi><mo>|</mo><mo>|</mo><mi>J</mi><mo>|</mo></mrow><mn>2</mn></msup></mfrac><munder><mi>&Sigma;</mi><mrow><mi>i</mi><mo>&Element;</mo><mi>I</mi><mo>,</mo><mi>j</mi><mo>&Element;</mo><mi>J</mi></mrow></munder><munder><mi>&Sigma;</mi><mrow><mi>u</mi><mo>&Element;</mo><mi>I</mi><mo>,</mo><mi>v</mi><mo>&Element;</mo><mi>J</mi></mrow></munder><mi>K</mi><mrow><mo>(</mo><msub><mi>x</mi><mi>ij</mi></msub><mo>,</mo><msub><mi>x</mi><mi>uv</mi></msub><mo>)</mo></mrow></mrow></math>]]></maths>对于均方残差阈值σ≥0,如果子矩阵B<sub>IJ</sub>满足MSR(I,J)≤σ,则称B<sub>IJ</sub>子矩阵为一个σ-核双聚类,这里K(x<sub>ij</sub>,x<sub>uv</sub>)为核函数,通过选择不同的核函数,对<maths num="0004"><![CDATA[<math><mfenced open='{' close=''><mtable><mtr><mtd><mi>min</mi><mi>f</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>MSR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow><mrow><mi>&sigma;gRVAR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>|</mo><mi>X</mi><mo>|</mo><mo>|</mo><mi>Y</mi><mo>|</mo></mrow><mrow><mo>|</mo><mi>I</mi><mo>|</mo><mo>|</mo><mi>J</mi><mo>|</mo></mrow></mfrac></mtd></mtr><mtr><mtd><mi>subject to MSR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow><mo>&lt;</mo><mi>&sigma;</mi><mo>,</mo><mi>&sigma;</mi><mo>&GreaterEqual;</mo><mn>0</mn></mtd></mtr></mtable></mfenced></math>]]></maths>来求解优化问题,找到各种可重叠双聚类;其中,|X|、|Y|分别为表达实数矩阵A的行数、列数;|I|、|J|分别表示所求双聚类的实数矩阵X、Y子集的行数和列数;(3.1)初始抗体群:在一个L维空间中,初始化生成K个差异种群G<sub>i</sub>,与K个双聚类对应;每个子种群G<sub>i</sub>有100个初始抗体{z<sub>i1</sub>,z<sub>i2</sub>,...,z<sub>ij</sub>},i=1,2...,K,j=1,2,...100;定义每个抗体z<sub>ij</sub>的编码为二进制编码,其长度为L=m+n,这里n和m分别为二维数据整体行的个数和列的个数;如果抗体z<sub>ij</sub>某一位置为l,则意味着相应的行或者列包括在K个双聚类中;初始化设置迭代终止均方残差MSR阈值σ;(3.2)设计抗体适应度函数affinity(z<sub>ij</sub>),计算每个子种群G<sub>i</sub>中每个抗体z<sub>ij</sub>的亲和度值,这里又称适应度值(i=1,2...,K,j=1,2,...100):<maths num="0005"><![CDATA[<math><mrow><mi>affinity</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ij</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>MSR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow><mrow><mi>&sigma;gRVAR</mi><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>|</mo><mi>X</mi><mo>|</mo><mo>|</mo><mi>Y</mi><mo>|</mo></mrow><mrow><mo>|</mo><mi>I</mi><mo>|</mo><mo>|</mo><mi>J</mi><mo>|</mo></mrow></mfrac></mrow></math>]]></maths>其中,I,J分别为动态所求的子矩阵行个数和列个数,RVAR为动态核子矩阵总体行方差和MSR为动态核子矩阵平均平方残差;(3.3)抗体促进和抑制:各个子种群G<sub>i</sub>并行计算,计算子种群中每个抗体z<sub>ij</sub>的抗体浓度:<maths num="0006"><![CDATA[<math><mrow><mi>C</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ij</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><munderover><mi>&Sigma;</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mo>|</mo><mi>affinity</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ij</mi></msub><mo>)</mo></mrow><mo>-</mo><mi>affinity</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ik</mi></msub><mo>)</mo></mrow><mo>|</mo></mrow></mfrac></mrow></math>]]></maths>其中子种群G<sub>i</sub>中抗体浓度C(z<sub>ij</sub>)与适应度affinity(z<sub>ij</sub>)相关,N为当前子种群G<sub>i</sub>中的抗体个数;在各子种群内,根据抗体的浓度调节机制,由抗体的期望繁殖率p<sub>j</sub>对抗体有选择地进行克隆复制;设在当前子种群G<sub>i</sub>中,抗体z<sub>ij</sub>的期望繁殖率p<sub>j</sub>是由适应度概率p<sub>fj</sub>和浓度抑制概率p<sub>dj</sub>两部分组成:p<sub>j</sub>=μ·p<sub>fj</sub>+(1-μ)P<sub>dj</sub>其中,μ是常数调节因子,取值为0.6;适应度概率p<sub>fj</sub>为:<maths num="0007"><![CDATA[<math><mrow><msub><mi>p</mi><mi>f</mi></msub><mrow><mo>(</mo><msub><mi>z</mi><mi>ij</mi></msub><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>affinity</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ij</mi></msub><mo>)</mo></mrow></mrow><mrow><munderover><mi>&Sigma;</mi><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>affinity</mi><mrow><mo>(</mo><msub><mi>z</mi><mi>ik</mi></msub><mo>)</mo></mrow></mrow></mfrac><mo>;</mo></mrow></math>]]></maths>抗体浓度抑制概率p<sub>dj</sub>为:<img file="FDA00003052950500032.GIF" wi="478" he="202" />抗体的浓度高则减小该个体的选择概率,反之则增加该个体的选择概率;(3.4)免疫选择克隆就是将免疫选择算子选中的抗体个体进行复制,于是被选中的抗体z<sub>ij</sub>的克隆复制个数为:<maths num="0008"><![CDATA[<math><mrow><msub><mi>N</mi><mi>c</mi></msub><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mi>round</mi><mrow><mo>(</mo><mi>v</mi><mo>&CenterDot;</mo><msub><mi>p</mi><mi>j</mi></msub><mo>)</mo></mrow></mrow></math>]]></maths>其中,round表示取整;v&gt;1是参数因子,这里取值为当前种群个数的2倍;p<sub>j</sub>是期望繁殖率;(3.5)子种群内抗体间交叉,新抗体产生;每个抗体被选择交叉的概率为P<sub>c</sub>;P<sub>c</sub>=P<sub>1</sub>·r<sup>t/T</sup>这里P<sub>1</sub>是一个预先设定的概率,r∈[0,1]是一个变异因子常数,在这里可设为一个[0,1]的随机数,它起着调整交叉概率;t为当前演化代数,T为最大可迭代代数;依据概率P<sub>c</sub>从当代种群中选取若干个体,按照交叉算子<maths num="0009"><![CDATA[<math><mfenced open='{' close=''><mtable><mtr><mtd><msubsup><mi>G</mi><mn>1</mn><mi>new</mi></msubsup><mo>=</mo><msub><mi>&omega;</mi><mn>1</mn></msub><mo>&CenterDot;</mo><msub><mi>G</mi><mn>1</mn></msub><mo>+</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>&omega;</mi><mn>1</mn></msub><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>G</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><msubsup><mi>G</mi><mn>2</mn><mi>new</mi></msubsup><mo>=</mo><msub><mi>&omega;</mi><mn>2</mn></msub><mo>&CenterDot;</mo><msub><mi>G</mi><mn>2</mn></msub><mo>+</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><msub><mi>&omega;</mi><mn>2</mn></msub><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>G</mi><mn>1</mn></msub></mtd></mtr></mtable></mfenced></math>]]></maths>进行变异,其中G<sub>1</sub>,G<sub>2</sub>为从种群中随机选择的两个父个体,<img file="FDA00003052950500035.GIF" wi="207" he="76" />为通过交叉运算子运算后产生的子代对应新个体;ω<sub>1</sub>,ω<sub>2</sub>为[0,1]上随机选取的参数;(3.6)在各子种群内,以突变概率P<sub>m</sub>,执行抗体间的变异,其中<maths num="0010"><![CDATA[<math><mrow><msub><mi>P</mi><mi>m</mi></msub><mo>=</mo><msub><mi>P</mi><mn>2</mn></msub><msup><mi>e</mi><mrow><mo>-</mo><msub><mi>r</mi><mn>2</mn></msub><mi>t</mi><mo>/</mo><mi>T</mi></mrow></msup></mrow></math>]]></maths>这里,P<sub>2</sub>代表了P<sub>m</sub>的初始值,P<sub>2</sub>设置为0.5;r<sub>2</sub>为抗体突变参数,t为当前演化代数,T为最大可迭代代数;依据概率P<sub>m</sub>从中选取若干个体,按照变异算子<maths num="0011"><![CDATA[<math><mrow><msup><mi>V</mi><mo>&prime;</mo></msup><mo>=</mo><mfenced open='{' close=''><mtable><mtr><mtd><mi>int</mi><mrow><mo>(</mo><mi>V</mi><mo>+</mo><mrow><mo>(</mo><msub><mi>b</mi><mi>sup</mi></msub><mo>-</mo><mi>V</mi><mo>)</mo></mrow><msup><mrow><mo>[</mo><msup><mi>e</mi><msub><mrow><mo>-</mo><mi>r</mi></mrow><mn>2</mn></msub></msup><mo>&CenterDot;</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mfrac><mi>t</mi><mi>T</mi></mfrac><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>)</mo></mrow><mo>,</mo><mi>sign</mi><mo>=</mo><mn>0</mn></mtd></mtr><mtr><mtd><mi>int</mi><mrow><mo>(</mo><mi>V</mi><mo>-</mo><mrow><mo>(</mo><mi>V</mi><mo>-</mo><msub><mi>b</mi><mi>inf</mi></msub><mo>)</mo></mrow><msup><mrow><mo>[</mo><msup><mi>e</mi><msub><mrow><mo>-</mo><mi>r</mi></mrow><mn>2</mn></msub></msup><mo>&CenterDot;</mo><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mfrac><mi>t</mi><mi>T</mi></mfrac><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>)</mo></mrow><mo>,</mo><mi>sign</mi><mo>=</mo><mn>1</mn></mtd></mtr></mtable></mfenced></mrow></math>]]></maths>进行变异;其中,V'为变异后的参数;V是选中的变异参数,这里V=0.4;sign随机取0或者1;b<sub>sup</sub>=1和b<sub>inf</sub>=0分别为参数取值的上界和下界;int()为取整;(3.7)亲和度计算:重新计算子种群内每个抗体适应度affinity(z<sub>ij</sub>);(3.8)最优抗体:将不同子种群中的适应度函数affinity(z<sub>ij</sub>)最小的抗体作为局部最优抗体保留;(3.9)抗体记忆池:子种群间筛选的最优抗体,存储入抗体记忆池中,形成精英抗体群;(3.10)对精英抗体群中每个抗体,采用单纯形法进行Memetic局部搜索,获取每个抗体的局部最优值,这些最优值更新精英抗体记忆池;所述局部最优值是指抗体的适应度函数affinity(z<sub>ij</sub>)稳定于某个值,且适应度函数affinity(z<sub>ij</sub>)小于种群中的局部最优个体;(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个双聚类,其选择计算标准为:<maths num="0012"><![CDATA[<math><mrow><mi>sum</mi><mo>_</mo><mi>fitness</mi><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></munderover><mfrac><mrow><mi>MS</mi><msub><mi>R</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow><mrow><mi>&sigma;gRVA</mi><msub><mi>R</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>I</mi><mo>,</mo><mi>J</mi><mo>)</mo></mrow></mrow></mfrac></mrow><mfrac><mi>K</mi><mrow><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></munderover><mrow><mo>(</mo><msub><mi>I</mi><mi>i</mi></msub><mo>+</mo><msub><mi>J</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow></mfrac></math>]]></maths>其中,i=1,2,...,K表示在某个核函数下的K个双聚类;选择最小的sum_fitness所对应的核函数下的K个双聚类为输出的K个双聚类,其对应K个协同表达基因集,包括:1)相似表达的基因集,共同高表达或共同低表达,变化几乎相同基因集;2)表达趋势一致的基因集,共同高表达或共同低表达,变化一致基因集;3)反向表达相关的基因集,基因表达总相反;输出协同表达的基因集。
地址 201620 上海市松江区松江新城区人民北路2999号