数学建模华中赛b题优秀.docx
《数学建模华中赛b题优秀.docx》由会员分享，可在线阅读，更多相关《数学建模华中赛b题优秀.docx（20页珍藏版）》请在冰豆网上搜索。
数学建模华中赛b题优秀
第八届华中地区大学生数学建模邀请赛
承诺书
我们仔细阅读了第八届华中地区大学生数学建模邀请赛的竞赛细则。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。
如有违反竞赛规则的行为，我们将受到严肃处理。
我们的参赛报名号为：
参赛队员（签名）：
队员1：
队员2：
队员3：
武汉工业与应用数学学会
第八届华中地区大学生数学建模邀请赛组委会
第八届华中地区大学生数学建模邀请赛
编号专用页
选择的题号：
B
参赛的编号：
（以下内容参赛队伍不需要填写）
竞赛评阅编号：
第八届华中地区大学生数学建模邀请赛
题目：
基因调控网络的重构及病毒感染的致病机制
【摘要】
一个基因的表达受其他基因的影响,而这个基因又影响其他基因的表达,这种相互影响相互制约的关系构成了复杂的基因调控网络。
基因调控网络的研究是从基因之间相互作用的角度揭示复杂的生命现象，是当前生物信息学研究的前沿。
疾病的发病因素和原理，对于医疗领域有着十分重要的作用。
这不仅仅能够让更多的患者免受病痛的困扰，还能促进人类医学史的进步。
所以根据基因数据谱来重构基因调控网络，以及某个疾病症状产生的原因的研究具有很大的意义。
本文对基因调控网络的重构以及导致严重临床症状的蛋白质进行了研究和推测。
由于所给的基因数据谱（附录一）十分庞大，所以首先要对数据进行降维处理。
本题基于时间序列给出了272组基因数据，为了减小噪声以及缺失值对实验精度的干扰，在实验前对四组噪声较大或有缺失的数据进行剔除。
具体的降维方式采用了多元统计法中的主成分分析和聚类分析：
先对这一万多个数据做主成分分析，从这一万多个数据中，通过线性变化选出了1000个左右的重要变量来组成新的样本。
既降低了数据的处理难度，又尽量保持了新数据和原数据相比，尽可能保持原数据的信息。
然后用spss两阶聚类法粗略地对要聚类的数目进行一个估计，根据此估计用K-means算法对数据进行处理，得到相应的30组数据。
对这30组数据建立模型，来重构基因调控网络。
本文中采用的模型是线性回归模型，并对它的合理性，以及相较贝叶斯网络作了对比。
最后依据所得到的系数矩阵进行基因网络图的绘制与呈现。
问题二在第一问的基础上，寻求导致产生严重临床症状的蛋白质。
根据附录二给出的个体出现感染症状时间节点示意图，1代表此志愿者在该时间节点表现出了临床症状，0则表示没有，这是一个二分类。
本题采用逻辑回归模型，利用LR分类器模型去寻找该重要蛋白质。
用268组数据，其中每一个基因视为该组数据的一个属性，对这些基因进行LR分类，并得到相应的系数矩阵。
然后对系数矩阵进行分析，取出影响比较大的几个基因，然后对照基因表对基因作用的描述去寻求该重要蛋白。
本题最终找出四个导致志愿者产生严重的临床症状的蛋白质。
所有代码实现，以及每次得到的系数矩阵均在附录中给出。
关键词：
线性回归模型，基因调控网络重构，多元统计法，主成分分析，聚类分析，逻辑回归（LR分类器模型）
1.问题重述
通过基因之间的相互调控，生物体可以实现细胞的生长，器官的发育、以及免疫等各种生物机能。
随着测序技术的发展，产生了越来越多的高通量实验数据。
基于这些实验数据重建基因调控网络（Generegulatorynetworks,GRNs），对于深入了解生物机能的实现过程具有重要作用。
生物实验中，在17个健康志愿者鼻内接种流感病毒H3N2/Wisconsin，其中9个人出现了严重的感染症状，另外的8个人没有出现症状。
接种后，每隔大约8h从血液中采集样本测量基因表达谱数据，实验数据一共有16个时间点（单位:
h），包括baseline（-24），0，5，12，21，29，36，45，53，60，69，77，84，93，101，108，共268个样本。
基因表达谱数据见附件1，其中前8个为未出现严重感染症状的数据，后9个为出现严重感染症状的数据。
（其中行代表探针号，对应着不同的基因；列为各个个体血液样本在各个时间节点的数据）个体出现感染症状的时间节点示意图见附件2。
问题：
1）根据实验数据重构基因调控网络；
2）通过比较出现感染症状的志愿者和健康志愿者的样本数据，试确定病毒感染人体后导致志愿者是否会出现严重临床症状的重要蛋白。
2.问题分析
一个基因的表达受其他基因的影响,而这个基因又影响其他基因的表达,这种相互影响相互制约的关系构成了复杂的基因调控网络。
更一般些,几乎所有的细胞活动都被基因网络所控制。
生命是存储并加工信息的复杂系统,孤立地研究单个基因及其表达往往不能确切地反映生命现象本身的内在规律。
因此,需要从复杂系统的角度研究基因网络。
对于问题一，考察我们如何根据已有的基因表达谱（附录一）去重构基因调控网络，从而推断调控网络各节点之间潜在的调控关系。
考虑“反向分析法”来重构基因调控网络，常见的基因调控网络模型有布尔网络模型、线性组合模型和贝叶斯网络模型等等。
然而题目所给的数据集十分庞大，如果直接将这一万个基因全部带入模型，那么计算量是惊人的。
所以需要用到多元统计方法中的主成分分析和聚类分析去实现降维的操作。
对于问题二，在已经重构好的基因网络的基础上寻找导致病毒感染人体以后导致志援者是否产生严重临床症状的蛋白质。
首先我们要对数据进行分析，寻找与染病相关系数大的基因，然后依据附录一的sheet2中对于基因的描述去进一步确定关键蛋白质。
3.模型假设
针对本问题，建立如下合理假设：
（1）题目所给数据准确可靠；
（2）假设不考虑个体差异性；
（3）基因表达呈高斯分布；
4.符号说明
表示第n个基因基于时间序列的第m组数据；
表示一个基因；
为回归系数；
代表基因X在时间点t具有的表达值；
为常数；
为误差项。
5.问题一的建模与算法实现求解
5.1数据的分析
问题一需要根据所给的基因表达谱数据来重构基因调控网络，附录一中的sheet1中给出了17个志愿者体内的10000种基因，随着注入病毒后的时间变化而出现的数值变化。
由于数据集过大，所以第一步要做的就是对这一万种基因进行筛选降维操作。
只选取部分具有代表性的数据代入模型，从而减少计算量。
对于数据的处理部分，采用多元统计中的常用方法，主成分分析和聚类分析。
5.2数据预处理
5.2.1数据处理方法选择
由于这道题目的数据量庞大，所以，如何筛选数据就成了很重要的一步。
我们这里采取先对10000组数据做主成分分析，形成1000组新变量，再对这些新变量进行聚类分析，进一步降维。
5.2.2主成分分析
（1）主成分分析的基本思想：
主成分分析的基本思想是通过构造10000个基因初始数据的适当的线性组合，以产生一系列互不相关的新变量，从中选出少数几个新变量并使它们尽可能多地包含原先所有基因的信息（降维），从而使得用这几个新变量替代原变量分析问题成为可能。
即在尽可能少丢失信息的前提下从所研究的
个变量中求出几个新变量，它们能综合原有变量的信息，相互之间又尽可能不含重复信息。
（2）主成分分析的实现：
设有
个样品，
个变量（指标）的数据矩阵。
本题中n=10000，表示10000种基因；m=268，表示基于时间序列的基因数据变化指标。
寻找
个新变量
，使得
1、
2、
彼此不相关
主成分的系数向量
的分量
刻划出第
个变量关于第
个主成分的重要性。
可以证明，若
为
维随机向量，它的协方差矩阵
的
个特征值为
，相应的标准正交化的特征向量为
，则
的第
主成分为
。
称
为主成分
的贡献率，
为主成分
的累计贡献率，它表达了前
个主成分中包含原变量
的信息量大小，通常取
使累计贡献率在85%以上即可。
当然这不是一个绝对不变的标准，可以根据实际效果作取舍，例如当后面几个主成分的贡献率较接近时，只选取其中一个就不公平了，若都选入又达不到简化变量的目的，那时常常将它们一同割舍。
计算步骤如下：
1、由已知的原始数据矩阵
计算样本均值向量
；
其中
2、计算样本协方差矩阵
其中
3、把原始数据标准化，即
，记
。
形成样本相关矩阵
；
4、求
的特征根
及相应的标准正交化的特征向量
，可得主成分为
。
（3）主成分分析降维结果
用Matlab实现以上算法（代码见附录），实现结果如下：
图5.1主成分分析结果
如图可见是一个1000组新的变量，由于数据集比较大，在这里只截出一部分。
下面再对这1000组新变量做聚类分析处理。
5.2.3聚类分析
（1）聚类分析的基本思想：
聚类（clustering），简单的讲就是将一个给定的数据集分成若干个不同簇的过程。
聚类算法中的簇指的是数据对象的集合，且这种数据对象集合必须满足条件：
同一簇中的数据对象间具有较大的相似性，而不同簇中的数据对象间具有较小的相似性。
聚类的主要指导思想就是尽可能使同一簇内对象相似度达到最大，且不同簇间对象相异度达到最大。
（2）K-means算法：
首先从含有n个数据对象的数据集中随机选择K个数据对象作为初始中心然后计算每个数据对象到各中心的距离，根据最近邻原则，所有数据对象将会被划分到离它最近的那个中心所代表的簇中，接着分别计算新生成的各簇中数据对象的均值作为各簇新的中心，比较新的中心和上一次得到的中心，如果新的中心没有发生变化，则算法收敛，输出结果，如果新的中心和上一次的中心相比发生变化，则要根据新的中心对所有数据对象重新进行划分。
直到满足算法的收敛条件为止。
（3）K-means算法的实现：
①从含有1000个数据对象的基因表达谱（附录一）中随机选择100个数据对象作为初始的聚类中心；
②两基因数据的相似度可通过计算两个基因数据的欧式距离来得到,再根据最近邻原则将数据对象逐个划分到离其最近的聚类中心所代表的簇中，计算误差平方和准则函数E的值；
③更新聚类中心，即分别计算各个簇中所有数据对象的均值作为各个簇的新的中心，以新的聚类中心来计算误差平方和准则函数E的值；
将步骤③计算得到的E值和前一次计算得到的E值进行比较，若两者差值的绝对⑤值小于等于预先设定的阈值，即聚类准则函数收敛，则转步骤（5），否则转步骤②；
⑤输出K个聚类。
5.2.4数据降维后的结果显示
将10000个基因数据通过主成分分析和聚类分析缩小到30组特征数据。
图5.2各类中基因数目分布
5.3模型一的分析
目前重构基因调控网络的模型主要有贝叶斯网络，线性组合模型等，本题我们采用的是线性回归模型，下面先对基因调控网络与线性回归模型的对应关系作分析，以证明模型的可行性：
假设基因
受基因
.....
调控，则认为相应的基因表达值，有以下线性关系：
式中
代表基因
在t时刻的表达水平，
....
为常数。
若基因A，B，C之间的真实调控关系如图所示，其中A，B，C代表基因，而边代表调控关系，比如从A到B有一条有向边，代表了基因A对基因B有调控关系。
A
BC
其对应的线性模型为：
5.4模型一的建立
5.4.1模型说明：
采用多元线性回归模型
设因变量y与自变量
的线性回归模型为
（2-1）
其中，
，
，....,
是p个未知参数，称为回归系数。
y称为因变量（被解释变量），而
，
，...，
是p歌可以精确测量并可以控制的一般变量，称为自变量（解释变量）。
上面的公式就称为多元线性回归模型。
为随机误差项，对随机误差项我们常假定
（2-2）
称
为理论回归方程。
如果有n个样本数据（
），i=1，2，...，n，则样本数据代入（2-3）式可表示为：
写成矩阵形式为：
5.4.2本题中的基因调控网络与线性回归模型的对应关系：
17位志愿者，每位志愿者体内的P个基因数据，这里P=11961。
根据时间序列排布，一共16个连续采样时间点，获得272组观测数据。
其中有四组错误数据，排除后剩余268组数据
，这里t=1，2，....,268，其中，xti，表示t时刻时基因i（i代表一个基因）的表达值。
对于基因i，i=1，2，...,p,若假设它与全部基因（包括它本身在内）存在线性关系，则它的线性回归模型为
（3-5）
代入268组样本数据，得：
写成矩阵的形式表示为：
也即
式中，
矩阵X是一个（T-1）*P矩阵，称X为回归设计矩阵或者设计矩阵。
对于
所有P个基因的全模型，可以表示成统一的矩阵形式如下
Y=XB+E
其中，
。
在处理基因调控网络重构问题时，实际上是逐一寻找对基因i,i=1,2,...,p起调控作用的调控基因的集合，而不是从基因i,i=1,2,...,p调控哪些其他基因的角度出发。
这从网络的拓扑结构上来讲，我们是在逐个寻找网络中每个顶点（基因）的入边的集合（即其他基因对该顶点对应基因的调控关系），找到了每个顶点的入边集合后，由顶点和边组成的整个网络的拓扑结构就确定了。
以上线性回归算法用Matlab实现，重构出的基因调控网络如下所示：
图5.3基因网络重构效果图
6.问题二的建模与算法实现求解
6.1模型选择
问题二：
在已经重构好的基因网络的基础上寻找导致病毒感染人体以后导致支援者是否产生严重临床症状的蛋白质。
假设是否展现严重的临床症状为最后的因变量，那么由于这个因变量是二分类，所以考虑用逻辑回归模型去对基因对是否产生症状的关系做一个分析。
通过逻辑回归我们可以得到一个基因对于严重临床症状的相关性的大小，以及对产生临床症状是起着正向的作用还是抑制的作用。
我们要对所得的数据进行分析，寻找与染病相关系数大的基因，然后去依据基因表达谱中的数据中对于基因的描述去进一步确定关键蛋白质。
6.2模型简述
该模型的对输入没有严格的要求，而且模型简单直观，又容易解释.而且该模型并不容易产生过拟合。
而当前题目中只有产生和没有产生临床症状两种不同的结果，所以采用此模型去寻求重要的蛋白质是一个可以考虑的方法。
定义一种概率函数
（6-1）
对
进行线性回归。
对于上述模型可以采用最大似然估计方法对其回归参数进行估计。
既用总体的分布密度或概率分布的表达式及样本所提供的信息求位置参数估计两的一种方法。
6.3模型建立
注射后不同人随时间是否产生症状不同。
回归的因变量是二分类。
因此采用的是LR分类器的方法去求解。
将每一个志愿者的每一个时间节点作为一组数据，而每一个基因的值就是该组数据的一个属性。
如果一个志愿者在某时间点产生了严重的临床症状，则因变量结果为1，相反没有产生严重的临床症状则为0。
本次实验有17个志愿者，每个志愿者有16个时间节点的采样数据。
又由于有四个数据有较大误差，故对268组数据进行LR分类。
得到每一个基因与最终是否产生临床症状这个因变量的相关性。
通过相关性的分析便可以得到能够推动产生临床症状的关键基因，从而推出关键蛋白。
将通过逻辑回归所得到的系数矩阵绘制如下：
图6.1各基因与临床症状的相关性
6.4模型估计结果
设定一个阙值，将大于阙值的基因数据拿出来，对所拿出来的基因数据进行分析。
满足要求的基因在基因表达谱中的编号如下：
图6.2基因在基因谱中的位置
对应基因表达谱，将上述位置的基因编号拿出，方便后续数据分析，如下：
图6.3基因以及基因的指针编号
对照基因表达谱中的数据得到了四个对于最终的因变量影响较大的基因，在谱中寻找所对应的蛋白质的信息如下表：
Gene_ID
Symbol
Description
Location
Type（s）
1742
DLG4
discs,largehomolog4（Drosophila）
PlasmaMembrane
kinase
3669
ISG20
interferonstimulatedexonucleasegene20kDa
Nucleus
enzyme
27074
LAMP3
lysosomal-associatedmembraneprotein3
PlasmaMembrane
64135
IFIH1
interferoninducedwithhelicaseCdomain1
Nucleus
enzyme
表6.1对临床症状产生重要的蛋白质
综上，结果认为导致志愿者产生严重的临床症状的蛋白质主要有四种，其标志依次为DLG4,ISG20,LAMP3,IFIH1
7.模型的分析及改进
7.1对线性回归模型拟合优度的分析
假设Ki为第I个类用线性回归重构基因网络的拟合优度，那么该网络的平均拟合优度
（7-1）
带入数据计算得到，平均拟合优度大约在0.8左右，所以该网络具有一定的误差，但是大致对于一个基因网络还是有一定真是的反应。
7.2改进
在使用线性回归模型区重构基因网络的时候，虽然不容易产生拟合过度的情况，但是线性回归和逻辑回归都只能对变量只见一个线性的关系具有良好的刻画。
对于完全不相关的两组变量的时候数据的误差就会增大。
所以可结合非线性动态贝叶斯网络来重构基因网络，并且在逻辑回归之前尽早的取出不相关变量，对变量进行分析。
8.附录
图8.1基因网络系数矩阵1
图8.2基因网络系数矩阵2
图8.3基因网络系数矩阵3
图8.4基因网络系数矩阵4
代码：
LinerRegress.m
for（i=1:
30）
tep=[];
forj=1:
i-1
tep=[tep,new_gene（new_represent（j）,:
）'];
end
forj=i+1:
30
tep=[tep,new_gene（new_represent（j）,:
）'];
end
[b,bint,r,rint,stats]=regress（new_gene（i,:
）',tep）;
new_ratio（i,1:
29）=b;
assess（i,1:
4）=stats;
End
Netplot.m
%函数名netplot
%使用方法输入请helpnetplot
%无返回值
%函数只能处理无向图
%作者：
tiandsp
%最后修改：
2012.12.26
functionnetplot（A,flag）
%调用方法输入netplot（A,flag），无返回值
%A为邻接矩阵或关联矩阵
%flag=1时处理邻接矩阵
%flag=2时处理关联矩阵
%函数只能处理无向图
ifflag==1%邻接矩阵表示无向图
ND_netplot（A）;
return;
end
ifflag==2%关联矩阵表示无向图
[mn]=size（A）;%关联矩阵变邻接矩阵
W=zeros（m,m）;
fori=1:
n
a=find（A（:
i）~=0）;
W（a
（1）,a
（2））=1;
W（a
（2）,a
（1））=1;
end
ND_netplot（W）;
return;
end
functionND_netplot（A）
[nn]=size（A）;
w=floor（sqrt（n））;
h=floor（n/w）;
x=[];
y=[];
fori=1:
h%使产生的随机点有其范围，使显示分布的更广
forj=1:
w
x=[x10*rand
（1）+（j-1）*10];
y=[y10*rand
（1）+（i-1）*10];
end
end
ed=n-h*w;
fori=1:
ed
x=[x10*rand
（1）+（i-1）*10];
y=[y10*rand
（1）+h*10];
end
plot（x,y,'r*'）;
title（'网络拓扑图'）;
fori=1:
n
forj=i:
n
ifA（i,j）~=0
c=num2str（A（i,j））;%将A中的权值转化为字符型
text（（x（i）+x（j））/2,（y（i）+y（j））/2,c,'Fontsize',10）;%显示边的权值
line（[x（i）x（j）],[y（i）y（j）]）;%连线
end
text（x（i）,y（i）,num2str（i）,'Fontsize',14,'color','r'）;%显示点的序号
holdon;
end
end
end
end
9.参考文献
[1]虞慧婷,吴骋,柳伟伟,付旭平,贺佳.基因调控网络模型构建法.第二军
医大学学报.2006,27（7）:
737~740
[2]StyczynskiMP,StephanopoulosG.Overviewofcomputationalmethodsfortheinferenceofgeneregulatorynetworks,computersandchemicalengineering[J].2005,29:
519-534.
[3]BansalAK,KoradiaV.Theroleofreverseengineeringinthedevelopmentofgenericformulations,pharmaceuticaltechnology[J].PharmaceatTechnol,2005,29：
50254.
[4]Wyrick,J.J，R.A.Young.Decipheringgeneexpressionregulatorynetworks.CurrOpinGenetDev2002,12
（2）:
130-6[5]徐肖江,王连水,丁达夫.从酵母表达时间序列估计基因调控网络[J].生物化学与生物物理报,2003,35:
7072716.
升级会员 