基因网络的重构

《基因网络的重构》由会员分享，可在线阅读，更多相关《基因网络的重构（14页珍藏版）》请在装配图网上搜索。

1、基因调控网络的重构及病毒感染的致病机制摘要一个基因的表达受其他基因的影响，而这个基因又影响其他基因的表达，这 种相互影响相互制约的关系构成了复杂的基因调控网络。几乎所有的细胞活动和 功能都受基因网络调控。本文通过对基因调控网络的重构来推测出病毒感染的致 病机制。由于所给的基因表达谱数据海量，所以首先对基于时间序列给出的272组基 因表达谱数据进行降维处理。为了减小噪声以及缺失值度对实验精度的干扰，在 实验前剔除四组噪声较大或有缺失的数据。对数据进行降维具体采用聚类分析 法：首先将8个健康志愿者的基因表达谱数据导入SPSS中，对数据进行K-means 聚类，将11966个基因聚为20类，分别找出

2、与各聚类中心距离最小的基因将其 作为特征基因。找出 20 个特征基因为(6002_at,5042_at,23378_at,54440_at,667_at,6192_at,84663_at,10098_at, 2671_a t,3117_a t,157680_a t,1372_a t,83440_a t,54915_a t,25880_a t,8846_a t, 10321_a t,7503_a t,8678_a t,9582_a t)。对于问题一运用线性组合模型。对于经聚类分析处理所得的20个特征基因, 分别写出它与其他特征基因以及包括其本身的线性回归模型。基于特征基因的基 因表达谱数据利用M

3、ATLAB求出其相关系数矩阵，进而判断各基因之间的调控关 系。通过此调控关系，运用MATLAB程序构建出基因调控网络。线性组合模型基 于时间序列基因表达数据，能够更准确的描述基因调控网络的动态性。对于问题二，将9个感染病毒志愿者与另8个健康志愿者在相应时刻的总体 平均基因表达数据进行比较，并根据9个志愿者出现症状的时刻不同，将其分为 5组，同样将健康志愿者也分为5组，通过定义模型来取出影响较大的基因，然后 对照基因表对基因作用的描述去寻求该重要蛋白。关键词：基因调控网络重构 K-means聚类 特征基因 线性组合模型一、问题重述通过基因之间的相互调控，生物体可以实现细胞的生长，器官的发育、以及

4、 免疫等各种生物机能。随着测序技术的发展，产生了越来越多的高通量实验数据。基于这些实验数据重建基因调控网络（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个样

5、本。基因表达谱数据见附件1,其中前8 个为未出现严重感染症状的数据，后9个为出现严重感染症状的数据。（其中行 代表探针号,对应着不同的基因；列为各个个体血液样本在各个时间节点的数据） 个体出现感染症状的时间节点示意图见附件2。问题：1）根据实验数据重构基因调控网络；2）通过比较出现感染症状的志愿者和健康志愿者的样本数据，试确定病毒感染 人体后导致志愿者是否会出现严重临床症状的重要蛋白。二、问题分析一个基因的表达受其他基因的影响，而这个基因又影响其他基因的表达，这 种相互影响相互制约的关系构成了复杂的基因调控网络。几乎所有的细胞活动和 功能都受基因网络调控。本文需要根据基因表达谱数据重构基因调控

6、网络以及对 导致出现严重临床症状的蛋白质进行研究。对于问题一，首先对数据进行聚类分析得到20个特征基因，然后建立线性 组合模型，基于基因表达谱数据利用MATLAB求出其相关系数矩阵，进而判断各 基因之间的调控关系，通过此调控关系，运用MATLAB程序构建出基因调控网络。对于问题二，由题目中所给的信息知9个健康志愿者在接种病毒后出现了相 应的临床症状，所以其体内部分蛋白肯定发生了变化。可以将出现症状后的数据 与前8个没有出现症状的相比较，如果其中某些基因表达谱的数据与正常基因表 达数据有较大的波动，则认为其感染病毒后体内出现了重要蛋白。三、模型假设与约定1假设不考虑个体之间的差异2假设所给的数据

7、都准确无误3假设各基因之间的调控关系是为线性的4假设得出的系数矩阵为正数是存在促进关系，为负数时存在抑制关系，零则表 示基因之间没有调控作用5假设误差项 是服从期望为0方差为b 2的正态分布函数四、符号说明及名词定义1.2.3.4.5.6.7.x , i = 1 ,2，/2=0 ;1 , 2 ,itA , i = 1,2,3,4,5ia ,i = 1,2,3,4,5iB , i = 1,2,3,4,5ib , i = 1,2,3,4,5i0 , i = 1,2,20is , i = 1,2,20基因i在t时刻的表示值出现症状的志愿者的组号各组对应基因总体平均表达谱数据 没有出现症状的志愿者的组

8、别各组对应基因总体平均表达谱数据线性组合模型对应的系数矩阵线性回归模型的误差五、模型的建立与求解5.1问题一的模型基因表达谱K-means 聚类基因组聚类中心特征基因线性组合模型基因调控网络图一流程图数据的处理基于题目提供的实验数据及说明重建基因调控网络，首先我们要分析从接种 病毒H3N2/Wisconsin各时间段内均没有出现症状的8个志愿者体内采集的血液 样本所测量的基因表达谱数据，即我们开始需要做的就是对这些数据进行筛选处 理，找出有代表性的特征基因以及其在志愿者体内各个时间段的变化数据。对于这些数据的筛选处理，我们使用的是Kmeans聚类方法。首先，将数 据从EXCEL导入SPSS中，

9、后对其进行聚类。根据上述所提到的聚类方法，将其 中具有相同或者类似作用的基因聚在一类，这里我们将这11966个基因分成20 类。由此程序运行之后，可以得到一系列的图表，如各个聚类所包含的基因数目 表，如下：类别代码12345678910案例个数1272156834317165715910751类别代码11121314151617181920案例个数9831890498700958805317951表一 各聚类所包含基因数目表每个聚类中的案例数图二聚类分布图还有每个基因所在的的聚类及与其聚类中心的距离合并表（附件）。通过后 面的这一张表格，我们找出与各聚类中心距离最小的基因，将其作为各聚类的特

10、征基因。因此，我们得到了 20个特征基因，基因的探针号及其所对应的序号如 下：基因探针6002_at5042_at23378_at54440_at667_at6192_at序号123456基因探针84663_at10098_at2671_at3117_at157680_at1373_at序号789101112基因探针83440_at54915_at25880_at8846_at10321_at7503_at序号131415161718基因探针8678_at9582_at序号1920表二基因探针号对应序号表模型的建立与求解通过上面的聚类方法，我们得到了 20个特征基因。对于这些基因，我们需 要通

11、过题目所给的数据求出各个基因之间的联系，即基因与基因之间所具有的调 控关系。由于实验中每隔一定时间都会对这些志愿者的血液做采样测量工作，从 而得到每个时间点上基因表达谱的数据。所以对于这样一个基于时间序列的基因 表达数据，我们可以建立一个线性组合模型。由这20个基因及实验所得到的这 些基因的16组（16个连续采样时间点）观测数据（基因的表达值） 讣xt2,xt,2051二1,2,16,其中，卞1 = 12，20表示基因i在t时刻的表示值。 对于基因i,二匕2,，20，若假设它与其他全部基因（包括其本身）存在线性 关系，则它的线性回归模型为：X =卩 X + 卩 x H卩 x +i=1,2；,s

12、tii1 t-1,1i 2 t-1,2is t-1,si此模型称为关于基因i的全模型。因题目中有8个志愿者没有出现症状，且于每个志愿者，都对其各个时间段 的血液进行了采样测量。若对这些志愿者进行逐个分析的话，则所需要构建的方 程组的数目很多。由于每个志愿者在同一个时刻其体内的基因表达谱数据都大径 相同，所以我们可以用其均值来代替。这样我们就得到了这8个志愿者在这16 个时间段里面的平均基因表达谱数据，并将其汇总成表，下表为部分数据。Hour 基因BaselineHour 00Hour 005Hour 0126002 at5.7351465165.6745275315.7173740335.72

13、34783785042 at13.1120851913.0997956913.0939249613.167355323378 at8.0904983798.1320371818.1058126588.23630825154440 at11.8324369911.8240585211.6782359311.79153007667 at4.8448019134.8395511814.895735834.8586254356192 at11.8374999711.980239112.0161526112.1392600784663 at7.9564622967.7873894917.70919020

14、67.67417541310098 at9.7811938559.5923621469.4855521069.662105892671 at6.457163146.4561789296.4713552756.4814495973117 at10.1267395410.2246935410.3288319610.55447331157680 at7.0150843857.1058961287.0345967756.9665245811373 at4.0528464234.0331791874.0617862894.01160286683440 at10.4049415510.478871210.

15、4123292510.4601505354915 at9.7052267129.7740230979.7184067479.86327548325880_at7.6387606537.6930863877.674135087.7041838738846 at8.4630967218.5314779938.5418774058.50012686210321 at5.2876413025.30073825.2075033954.8078558757503 at3.9801596843.9568949513.9318787073.9267675898678 at9.10426669.15475856

16、29.0754370039.1052674689582 at7.1622454557.0923982147.047716487.022314327表三部分平均基因表达数据综上，经过运算和整理，最终我们得到了 16个样本数据，对于基因i,我 们将这16个样本数据带入上述的全模型中，得到的方程组可以表示为：X = P X + P XX +2ii1 11 i 2 12is 1si1x = P x + P x +P X +s 3ii1 21 i2 22is 2si2X = P X + P X+P X +SJ mii1 m1,1i2 m1,2is ml,si,ml若令：y = (x ,x ,x )t x

17、 = (x ,x ,x )ti2i 3imii1i 2 im _1,i，P= (P , P，P )t，i = 1,Si i1 i2is则上述的方程组可表示为:y =0 x + 0 x + + 0 x+ iii1 1i 2 2is s也可看作：其中：y =i二 X 0 +iiX = x, x1 2,x =sx11x21x12x22x1sx2sx1- m-1,1xm-1,2xm1-1将上面EXCEL中的平均值表导入MATLAB中，利用相应的函数命令，求出基 因i所对应的系数矩阵0 i二12,20的值，如下：i下表为基因i的部分系数矩阵：0.00000.00000.00000.00000.00000

18、.00000.00000.00000.29990.00000.00000.00000.00000.0000-1.99650.00000.00000.00000.51973.17490.03040.98151.89497.91090.0000-0.24420.00000.00000.0000-1.01340.00000.00000.00000.76140.00000.00000.0000-3.4107-1.0881-1.04390.03290.0124-0.1295-0.6203-0.06340.00000.00000.0000-0.0494-0.0471-0.3509-0.10330.0296

19、-0.87930.3642-3.0061-0.0770-0.1476-0.30520.24620.18700.8681-0.10203.01290.58530.46561.9047-0.25510.26153.46220.7409-10.7046-0.15550.2314-0.05450.1586-0.1705-1.0649-0.38404.3068-0.04100.0000-0.65302.01050.52160.00000.00004.85090.00000.83930.00000.00000.00000.00003.38860.0000-0.26750.4764-0.1972-0.380

20、5-0.29561.24641.37913.13800.00000.00000.00000.00000.00000.00000.00000.0000-0.23070.0000-0.1546-0.9423-0.14711.64192.4032-8.75020.3333-0.69690.33750.35940.7069-1.6878-5.19372.8845-0.0087-0.03250.1065-0.0905-0.0396-0.4041-0.10590.4277-0.5995-0.1416-0.0051-1.5385-0.80060.3475-1.0721-5.47940.57670.2111-

21、0.2934-2.84750.35742.64410.0000-3.7557-0.07140.69470.55370.3532-0.29790.00700.7536-2.1191表三基因系数矩阵表所以，通过求出的系数矩阵的值，我们可以初步地了解到每个基因之间的调 控关系。用MATLAB中的程序（见附录），画出相应的基因调控网络。图二基因调控网络无向图由于网络的复杂性用MATLAB做出的无向图看起来杂乱无章，我们尝试加以 改进。图是用netplot函数做的无向图，我们尝试做一下改进：1：将无向图改为有向图；2：将边的权值改为-1与1 （权值大于0的改为1,小于0改为-1）然后用MATLAB程序做

22、成图（代码见附录）。图四基因调控网络有向图5.2问题二的模型问题的分析对于问题二，由题目中所给的信息知9个健康志愿者在接种病毒后出现了相 应的临床症状，其体内部分蛋白肯定出现了变化，所以我们可以将出现症状后的 数据与前8个没有出现症状的相比较，如果其中某些基因表达谱的数据与正常基 因表达数据有较大的波动，则认为其感染病毒后体内出现了重要蛋白。模型的建立与求解由于接种病毒后出现症状的志愿者人数与没有出现症状的人数不一致，比较 起来较为麻烦，而各志愿者在同一时刻体内的基因表达数据大致一样，所以我们 用这9个志愿者在出现症状之后的总体平均基因表达数据与另8个人的在相应时 刻的总体平均数据进行比较。由

23、附件二可知这9个志愿者出现症状的时刻有些不 同，编号为01,05,06,07的人在第45个小时出现症状，编号为13,15的人在第 101个小时出现症状，其他三个人出现症状的时间均不相同，所以将其分为5组, 分别记为A,i = 123,4,5其对应的总体平均表达数据记为a ,而将没有出现症状ii的志愿者也同样分为5组，记为B，相应的数据记为b .所以我们可以建立下面ii的模型：a - b .p = i i,i =12345 ibi若其比值的绝对值大于0.5,则认为该基因出现了重要蛋白。将整理后的数值带入上面的模型中，可以得出比值大于0.4的基因分别是：1. 3117_atHLA一DQA1 maj

24、or histocompatibility complex , class II, DQ alpha 1 Plasma membrane（质膜）t ransmembrane recep tor （胯膜受体）2. 9584_atRBM39 RNA binding motif protein 39Nucleus （细胞核）transcrip tion regula tor （转录调控因子）3. 25880_atTMEM186 transmembrane protein 186Cyto plasm（细胞质）ot her六、模型评价与改进6.1 :对问题一模型的评价优点： 模型是连续的，不需要对样本数

25、据离散化，避免了信息丢失。 基于时间序列基因表达数据，能够更准确的描述基因调控网络的动态 性。缺点：K-means聚类结果分布略有不均 基因调控关系及其复杂，往往很难通过线性关系来准确描述 算法的有效性难以验证改进：通过聚类分析将基因分为20组，在进行基因调控网络重构时，不仅可以 选取特征基因在组间进行重构，还可以在组内进行重构。即在组内再分组，不断 进行重构，最终可构建出一个较为完整的基因调控网络。而对于线性组合模型只 能对变量的一个线性关系具有良好的刻画，对不相干的两组变量刻画线性关系误 差就会增大。所以可结合非线性动态贝叶斯网络来重构基因网络。6.2：对问题二模型的评价优点：利用控制变量

26、法减少了干扰因素，利用求平均的方法消除了偶然因素，使 得找到的基因更加准确。缺点：模型比较简单，除了考虑重要蛋白开始大幅增加或减少的时间点之外，没 有考虑到基因调控的因素。七、参考文献1 刘世平，数据挖掘技术及其应用，高等教育出版社2010.1,31-37。2 刘佳宝，基于线性回归模型的基因网络重构，哈尔滨工业大学硕士学位论文, 2009.6。3 虞慧婷 吴聘 柳伟伟 付旭平 贺佳，基因调控网络构建方法，第二军医 大学学报，2006.7。4 杨斌，基因调控网络构建方法，济南大学，2007。八、附录由相关系数矩阵做出无向图的程序Netplot.mfunction netplot(A,flag)i

27、f flag=1ND_netplot(A);return;endif flag=2m n=size(A);W=zeros(m,m);for i=1:na=find(A(:,i)=0);W(a(l),a(2)=l;W(a(2),a(l)=l;endND_netplot(W);return;endND_netplot.mfunction ND_netplot(A)n n=size(A);w=floor(sqrt(n);h=floor(n/w);x=;y=;for i=1:hfor j=1:wx=x 10*rand(1)+(j-1)*10;y=y 10*rand(1)+(i-1)*10;endend

28、ed=n-h*w;for i=l:edx=x 10*rand(1)+(i-1)*10;y=y 10*rand(1)+h*10;endplot(x,y,k.);title(网络拓扑图)；for i=1:nfor j=1:nif A(i,j)=0 c=num2str(A(i,j);text(x (i)+x(j)/2,(y (i)+y(j)/2,c,Fontsize,8); line(x(i) x(j),y(i) y(j);end text(x(i),y(i),num2str(i),Fontsize,14,color,k); hold on;endend转换为无向图的程序A;IDS=A,B,C,D,E,F,G,H,T,J,K,L,M,N,O,P,Q,R,S,T; bg=biograph(A,IDS);set(bg.nodes,shape,circle,color,1,1,1,lineColor,0,0,0); set(bg,layoutType,radial);bg.showWeights=on;set(bg.nodes,textColor,0,0,0,lineWidth,2,fontsize,9); set(bg,arrowSize,12,edgeFontSize,9);get(bg.nodes,position);view(bg);