数学建模华中赛B题优秀论文

1、精选文库第八届华中地区大学生数学建模邀请赛承诺书我们仔细阅读了第八届华中地区大学生数学建模邀请赛的竞赛细则。我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的 , 如果引用别人的成果或其他公开的资料（包括网上查到的资料） ，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们的参赛报名号为：参赛队员(签名) ：队员 1：队员 2：队员 3：武汉工业

2、与应用数学学会第八届华中地区大学生数学建模邀请赛组委会-精选文库第八届华中地区大学生数学建模邀请赛编号专用页选择的题号：B参赛的编号：（以下内容参赛队伍不需要填写）竞赛评阅编号：-精选文库第八届华中地区大学生数学建模邀请赛题目：基因调控网络的重构及病毒感染的致病机制【摘要】一个基因的表达受其他基因的影响 ,而这个基因又影响其他基因的表达 ,这种相互影响相互制约的关系构成了复杂的基因调控网络。 基因调控网络的研究是从基因之间相互作用的角度揭示复杂的生命现象，是当前生物信息学研究的前沿。疾病的发病因素和原理，对于医疗领域有着十分重要的作用。这不仅仅能够让更多的患者免受病痛的困扰，还能促进人类医学史

3、的进步。 所以根据基因数据谱来重构基因调控网络，以及某个疾病症状产生的原因的研究具有很大的意义。本文对基因调控网络的重构以及导致严重临床症状的蛋白质进行了研究和推测。由于所给的基因数据谱（附录一）十分庞大，所以首先要对数据进行降维处理。本题基于时间序列给出了 272 组基因数据，为了减小噪声以及缺失值对实验精度的干扰， 在实验前对四组噪声较大或有缺失的数据进行剔除。 具体的降维方式采用了多元统计法中的主成分分析和聚类分析： 先对这一万多个数据做主成分分析，从这一万多个数据中，通过线性变化选出了 1000 个左右的重要变量来组成新的样本。既降低了数据的处理难度，又尽量保持了新数据和原数据相比，尽

4、可能保持原数据的信息。然后用 spss 两阶聚类法粗略地对要聚类的数目进行一个估计，根据此估计用 K-means 算法对数据进行处理， 得到相应的 30 组数据。对这 30 组数据建立模型，来重构基因调控网络。本文中采用的模型是线性回归模型，并对它的合理性， 以及相较贝叶斯网络作了对比。 最后依据所得到的系数矩阵进行基因网络图的绘制与呈现。问题二在第一问的基础上，寻求导致产生严重临床症状的蛋白质。根据附录二给出的个体出现感染症状时间节点示意图， 1 代表此志愿者在该时间节点表现出了临床症状， 0 则表示没有，这是一个二分类。本题采用逻辑回归模型，利用 LR 分类器模型去寻找该重要蛋白质。用 2

5、68 组数据，其中每一个基因视为该组数据的一个属性，对这些基因进行 LR 分类，并得到相应的系数矩阵。然后对系数矩阵进行分析， 取出影响比较大的几个基因， 然后对照基因表对基因作用的描述去寻求该重要蛋白。 本题最终找出四个导致志愿者产生严重的临床症状的蛋白质。所有代码实现，以及每次得到的系数矩阵均在附录中给出。关键词：线性回归模型，基因调控网络重构，多元统计法，主成分分析，聚类分析，逻辑回归（ LR 分类器模型）-精选文库1. 问题重述通过基因之间的相互调控， 生物体可以实现细胞的生长， 器官的发育、 以及免疫等各种生物机能。随着测序技术的发展，产生了越来越多的高通量实验数据。基于这些实验数据

6、重建基因调控网络（ Gene regulatory networks, GRNs），对于深入了解生物机能的实现过程具有重要作用。生物实验中， 在 17 个健康志愿者鼻内接种流感病毒H3N2/Wisconsin，其中9 个人出现了严重的感染症状， 另外的 8 个人没有出现症状。 接种后，每隔大约8 h 从血液中采集样本测量基因表达谱数据， 实验数据一共有 16 个时间点 ( 单位 : h) ，包括 baseline ( -24) ，0， 5， 12，21，29， 36，45，53，60，69， 77，84，93，101，108，共 268 个样本。基因表达谱数据见附件1，其中前 8 个为未出现严

7、重感染症状的数据， 后 9 个为出现严重感染症状的数据。 （其中行代表探针号，对应着不同的基因； 列为各个个体血液样本在各个时间节点的数据）个体出现感染症状的时间节点示意图见附件2。问题：1）根据实验数据重构基因调控网络；2）通过比较出现感染症状的志愿者和健康志愿者的样本数据，试确定病毒感染人体后导致志愿者是否会出现严重临床症状的重要蛋白。2. 问题分析一个基因的表达受其他基因的影响, 而这个基因又影响其他基因的表达, 这种相互影响相互制约的关系构成了复杂的基因调控网络。更一般些, 几乎所有的细胞活动都被基因网络所控制。生命是存储并加工信息的复杂系统, 孤立地研究单个基因及其表达往往不能确切地

8、反映生命现象本身的内在规律。因此 , 需要从复杂系统的角度研究基因网络。对于问题一，考察我们如何根据已有的基因表达谱（附录一）去重构基因调控网络，从而推断调控网络各节点之间潜在的调控关系。考虑“反向分析法”来重构基因调控网络， 常见的基因调控网络模型有布尔网络模型、线性组合模型和贝叶斯网络模型等等。 然而题目所给的数据集十分庞大，如果直接将这一万个基因全部带入模型， 那么计算量是惊人的。 所以需要用到多元统计方法中的主成分分析和聚类分析去实现降维的操作。对于问题二， 在已经重构好的基因网络的基础上寻找导致病毒感染人体以后导致志援者是否产生严重临床症状的蛋白质。首先我们要对数据进行分析， 寻找与

9、染病相关系数大的基因， 然后依据附录一的sheet2中对于基因的描述去进一步确定关键蛋白质。-精选文库3. 模型假设针对本问题，建立如下合理假设：（ 1）题目所给数据准确可靠；（ 2）假设不考虑个体差异性；（ 3）基因表达呈高斯分布 ；4. 符号说明X n,m 表示第 n 个基因基于时间序列的第m 组数据；A K 表示一个基因；i 为回归系数；at, X 代表基因X 在时间点t 具有的表达值；b1, b2 ,b3 为常数；1, 2 为误差项。5. 问题一的建模与算法实现求解5.1 数据的分析问题一需要根据所给的基因表达谱数据来重构基因调控网络，附录一中的 sheet1 中给出了 17 个志愿者

10、体内的 10000 种基因，随着注入病毒后的时间变化而出现的数值变化。 由于数据集过大， 所以第一步要做的就是对这一万种基因进行筛选降维操作。 只选取部分具有代表性的数据代入模型， 从而减少计算量。 对于数据的处理部分，采用多元统计中的常用方法，主成分分析和聚类分析。5.2 数据预处理数据处理方法选择由于这道题目的数据量庞大，所以，如何筛选数据就成了很重要的一步。我们这里采取先对10000 组数据做主成分分析，形成1000 组新变量，再对这些新变量进行聚类分析，进一步降维。主成分分析（ 1）主成分分析的基本思想：主成分分析的基本思想是通过构造 10000 个基因初始数据的适当的线性组合，以产生

11、一系列互不相关的新变量， 从中选出少数几个新变量并使它们尽可能多地包含原先所有基因的信息 （降维），从而使得用这几个新变量替代原变量分-精选文库析问题成为可能。即在尽可能少丢失信息的前提下从所研究的 m 个变量中求出几个新变量，它们能综合原有变量的信息，相互之间又尽可能不含重复信息。（ 2）主成分分析的实现：设有 n 个样品， m 个变量（指标）的数据矩阵。 本题中 n=10000，表示 10000 种基因； m=268，表示基于时间序列的基因数据变化指标。x11x12Lx1mx(1)X n mx21x22Lx2 mx(2)LLLLMxn1xn 2Lxnmx( n )寻找 k 个新变量 y1

12、, y2 ,L, yk (km) ，使得1、 yl al 1x1 al 2 x2 Lalm xm , (l1,2,L, k)2、 y1 , y2 ,L yk 彼此不相关主成分的系数向量 al(al 1, al 2 ,L , alm ) 的分量 alj 刻划出第 j 个变量关于第 l 个主成分的重要性。可以证明，若 x( x1 , x2 ,L , xm )T 为 m 维随机向量，它的协方差矩阵 V 的 m 个特征值为 12Lm0 ，相应的标准正交化的特征向量为u1,u2 ,L , um ，则x ( x1 , x2 ,L, xm )T 的第 i 主成分为 yiuiT x (i1,2,L , m)

13、。muiT x (ikm称 i/j 为主成分 yi1,2,L, m) 的贡献率，j /j 为主成分j1j 1j 1y1 , y2 ,Lyk的累计贡献率， 它表达了前 k 个主成分中包含原变量x1, x2 ,L , xm 的信息量大小，通常取 k 使累计贡献率在 85%以上即可。当然这不是一个绝对不变的标准，可以根据实际效果作取舍， 例如当后面几个主成分的贡献率较接近时， 只选取其中一个就不公平了， 若都选入又达不到简化变量的目的， 那时常常将它们一同割舍。计算步骤如下：1、由已知的原始数据矩阵Xn m 计算样本均值向量 ? x( x1 , x2 ,L , xm )T ；其中 xi1n1,2,L

14、 , m)xij (in j 1?11 ( sij ) ? ( ij )2、计算样本协方差矩阵Vnn其中 sij( xli xi )( xljx j )(i, j 1,2,L , m)l1-精选文库3、把原始数据标准化，即%xijxj ，记% 。形成样本相关矩阵xijX n m( xij )jj? %T %；R X X?的特征根 1 2Lm0 及相应的标准正交化的特征向量4、求 Ru1, u2 ,L , um ，可得主成分为 yiuiT x (i 1,2,L , m) 。（ 3）主成分分析降维结果用 Matlab 实现以上算法（代码见附录），实现结果如下：图 5.1 主成分分析结果如图可见是一

15、个1000 组新的变量，由于数据集比较大，在这里只截出一部分。下面再对这 1000 组新变量做聚类分析处理。聚类分析（ 1）聚类分析的基本思想 ：聚类（ clustering ），简单的讲就是将一个给定的数据集分成若干个不同簇的过程。聚类算法中的簇指的是数据对象的集合， 且这种数据对象集合必须满足条件：同一簇中的数据对象间具有较大的相似性， 而不同簇中的数据对象间具有较小的相似性。聚类的主要指导思想就是尽可能使同一簇内对象相似度达到最大，且不同簇间对象相异度达到最大。（ 2） K-means算法：-精选文库首先从含有 n 个数据对象的数据集中随机选择 K个数据对象作为初始中心然后计算每个数据对

16、象到各中心的距离， 根据最近邻原则， 所有数据对象将会被划分到离它最近的那个中心所代表的簇中， 接着分别计算新生成的各簇中数据对象的均值作为各簇新的中心， 比较新的中心和上一次得到的中心， 如果新的中心没有发生变化， 则算法收敛， 输出结果，如果新的中心和上一次的中心相比发生变化，则要根据新的中心对所有数据对象重新进行划分。 直到满足算法的收敛条件为止。（ 3） K-means算法的实现：从含有 1000 个数据对象的基因表达谱（附录一）中随机选择 100 个数据对象作为初始的聚类中心；两基因数据的相似度可通过计算两个基因数据的欧式距离来得到, 再根据最近邻原则将数据对象逐个划分到离其最近的聚

17、类中心所代表的簇中，计算误差平方和准则函数 E 的值；更新聚类中心，即分别计算各个簇中所有数据对象的均值作为各个簇的新的中心，以新的聚类中心来计算误差平方和准则函数E 的值；4 将步骤计算得到的 E 值和前一次计算得到的E 值进行比较，若两者差值的绝对值小于等于预先设定的阈值，即聚类准则函数收敛，则转步骤（5) ，否则转步骤；输出 K 个聚类。数据降维后的结果显示将 10000 个基因数据通过主成分分析和聚类分析缩小到30 组特征数据。图 5.2 各类中基因数目分布-精选文库5.3 模型一的分析目前重构基因调控网络的模型主要有贝叶斯网络， 线性组合模型等， 本题我们采用的是线性回归模型， 下面

18、先对基因调控网络与线性回归模型的对应关系作分析，以证明模型的可行性：假设基因 A k 受基因 A K1 ,A k2，. A km 调控，则认为相应的基因表达值，有以下线性关系：at ,kb1at 1,k1 b2at 1,k 2 . bmat 1,kmk式中 at ,k 代表基因 A K 在 t 时刻的表达水平，b1 . bm 为常数。若基因 A ，B，C 之间的真实调控关系如图所示，其中A ，B，C 代表基因，而边代表调控关系，比如从A 到 B 有一条有向边，代表了基因A 对基因 B 有调控关系。ABC其对应的线性模型为：at ,cb1at 1,a b2 at1,b1at ,bb3at 1,a

19、25.4 模型一的建立模型说明：采用多元线性回归模型设因变量y 与自变量 x1 , x2 ,., xp ,的线性回归模型为y1x 12 x 2 . p x p（2-1）其中，1 ，2 ，., p 是 p个未知参数，称为回归系数。y 称为因变量（被解释变量），而 x1 ， x 2 ，.， x p 是 p 歌可以精确测量并可以控制的一般变量，称为自变量（解释变量）。上面的公式就称为多元线性回归模型。为随机误差项，对随机误差项我们常假定-精选文库E()0var( )2（2-2）称 E( y)1x1 2 x2 . p xp为理论回归方程。如果有 n 个样本数据（ xi1 , xi 2 ,., xip

20、; yi ），i=1，2，.，n，则样本数据代入（ 2-3）式可表示为：写成矩阵形式为：y X本题中的基因调控网络与线性回归模型的对应关系：17 位志愿者，每位志愿者体内的 P 个基因数据，这里 P=11961。根据时间序列排布，一共 16 个连续采样时间点，获得 272 组观测数据。其中有四组错误数据，排除后剩余268 组数据 x t1 , xt 2 ,.， xtp ，这里 t=1 ，2， .,268，其中， xti，表示 t 时刻时基因 i （i 代表一个基因 ) 的表达值。对于基因 i ，i=1 ，2，.,p,若假设它与全部基因 （包括它本身在内） 存在线性关系， 则它的线性回归模型为（

21、 3-5）代入 268 组样本数据，得：写成矩阵的形式表示为：-精选文库也即yi xii式中，矩阵 X 是一个（ T-1 ）*P 矩阵，称X 为回归设计矩阵或者设计矩阵。对于所有 P 个基因的全模型，可以表示成统一的矩阵形式如下Y=XB+E其中， Yy1, y2 ,., yp T , B1, 2 ,.,p T , E1 , 2 ,., P T 。在处理基因调控网络重构问题时，实际上是逐一寻找对基因i,i=1,2,.,p 起调控作用的调控基因的集合， 而不是从基因 i,i=1,2,.,p 调控哪些其他基因的角度出发。 这从网络的拓扑结构上来讲， 我们是在逐个寻找网络中每个顶点 ( 基因 ) 的入

22、边的集合 ( 即其他基因对该顶点对应基因的调控关系 ) ，找到了每个顶点的入边集合后， 由顶点和边组成的整个网络的拓扑结构就确定了。以上线性回归算法用 Matlab 实现，重构出的基因调控网络如下所示：-精选文库图 5.3 基因网络重构效果图6. 问题二的建模与算法实现求解6.1 模型选择问题二：在已经重构好的基因网络的基础上寻找导致病毒感染人体以后导致支援者是否产生严重临床症状的蛋白质。 假设是否展现严重的临床症状为最后的因变量，那么由于这个因变量是二分类， 所以考虑用逻辑回归模型去对基因对是否产生症状的关系做一个分析。 通过逻辑回归我们可以得到一个基因对于严重临床症状的相关性的大小，以及对

23、产生临床症状是起着正向的作用还是抑制的作用。我们要对所得的数据进行分析， 寻找与染病相关系数大的基因， 然后去依据基因表达谱中的数据中对于基因的描述去进一步确定关键蛋白质。6.2 模型简述该模型的对输入没有严格的要求，而且模型简单直观，又容易解释 . 而且该模型并不容易产生过拟合。 而当前题目中只有产生和没有产生临床症状两种不同的结果，所以采用此模型去寻求重要的蛋白质是一个可以考虑的方法。定义一种概率函数Pr(Y 1| X1 x1, X2 x2Xn1（ 6-1）xn)11 a*e b1X1bnX n对 log() 进行线性回归。对于上述模型可以采用最大似然估计方法对其回归参数进行估计。 既用总

24、体的分布密度或概率分布的表达式及样本所提供的信息求位置参数估计两的一种方法。6.3 模型建立注射后不同人随时间是否产生症状不同。回归的因变量是二分类。 因此采用的是-精选文库LR 分类器的方法去求解。将每一个志愿者的每一个时间节点作为一组数据，而每一个基因的值就是该组数据的一个属性。 如果一个志愿者在某时间点产生了严重的临床症状，则因变量结果为 1，相反没有产生严重的临床症状则为 0。本次实验有 17 个志愿者，每个志愿者有 16 个时间节点的采样数据。 又由于有四个数据有较大误差，故对 268 组数据进行 LR分类。得到每一个基因与最终是否产生临床症状这个因变量的相关性。 通过相关性的分析便

25、可以得到能够推动产生临床症状的关键基因， 从而推出关键蛋白。 将通过逻辑回归所得到的系数矩阵绘制如下：图 6.1 各基因与临床症状的相关性6.4 模型估计结果设定一个阙值，将大于阙值的基因数据拿出来，对所拿出来的基因数据进行分析。满足要求的基因在基因表达谱中的编号如下：图 6.2 基因在基因谱中的位置对应基因表达谱，将上述位置的基因编号拿出，方便后续数据分析，如下：图 6.3 基因以及基因的指针编号对照基因表达谱中的数据得到了四个对于最终的因变量影响较大的基因，在谱中寻找所对应的蛋白质的信息如下表：-精选文库Gene_IDSymbolDescriptionLocationType(s)1742

26、DLG4discs, large homolog 4Plasmakinase(Drosophila)Membrane3669ISG20interferon stimulatedNucleusenzymeexonuclease gene 20kDa27074lysosomal-associated membranePlasmaLAMP3protein 3Membrane64135interferon induced withNucleusenzymeIFIH1helicase C domain 1表 6.1对临床症状产生重要的蛋白质综上，结果认为导致志愿者产生严重的临床症状的蛋白质主要有四种，其

27、标志依次为DLG4,ISG20,LAMP3,IFIH17. 模型的分析及改进7.1 对线性回归模型拟合优度的分析假设 Ki 为第 I 个类用线性回归重构基因网络的拟合优度，那么该网络的平均拟合优度i mK(Ki / m)i 1（7-1）带入数据计算得到，平均拟合优度大约在 0.8 左右，所以该网络具有一定的误差，但是大致对于一个基因网络还是有一定真是的反应。7.2 改进在使用线性回归模型区重构基因网络的时候， 虽然不容易产生拟合过度的情况，但是线性回归和逻辑回归都只能对变量只见一个线性的关系具有良好的刻画。对于完全不相关的两组变量的时候数据的误差就会增大。 所以可结合非线性动态贝叶斯网络来重构

28、基因网络，并且在逻辑回归之前尽早的取出不相关变量，对变量进行分析。-精选文库8. 附录图 8.1 基因网络系数矩阵1图 8.2基因网络系数矩阵2图 8.3基因网络系数矩阵3-精选文库图 8.4基因网络系数矩阵4代码： LinerRegress.mfor(i=1:30)tep=;for j=1:i-1tep=tep,new_gene(new_represent(j),:)'endfor j=i+1:30tep=tep,new_gene(new_represent(j),:)'endb,bint,r,rint,stats=regress(new_gene(i,:)',tep

29、);new_ratio(i,1:29)=b;assess(i,1:4)=stats;EndNetplot.m%函数名 netplot%使用方法输入请help netplot%无返回值%函数只能处理无向图%作者： tiandsp%最后修改：functionnetplot(A,flag)%调用方法输入netplot(A,flag)，无返回值%A为邻接矩阵或关联矩阵%flag=1时处理邻接矩阵%flag=2时处理关联矩阵%函数只能处理无向图ifflag=1%邻接矩阵表示无向图ND_netplot(A);return;-精选文库endifflag=2%关联矩阵表示无向图m n=size(A);%关联矩阵变邻接矩阵W=zeros(m,m);fori=1:na=find(A(:,i)=0);W(a(1),a(2)=1;W(a(2),a(1)=1;endND_netplot(W);return;endfunctionND_netplot(A)n n=size(A);w=floor(sqrt(n);h=floor(n/w);x=;y=;fori=1:h%使产生的