基因芯片数据的主成分分析

生物医学信号处理特别培养 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 用奇异值分解的方法对基因芯片数据进行主成分分析 （21142010 袁辉） 一． 原理简介 对基因阵列进行主成分分析，即提取主模式是一件十分有意义的工作。这里我们讨论用 奇异值分解的方法对基因阵列进行模式提取。通过奇异值分解，可以根据奇异值的大小和采 用的阈值提取出主要成分，画出图形，进而分析生理意义。在进行奇异值分解之前，对原始 数据进行 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 化预处理是十分必要的。通过下面的分析也可以看出做不做预处理对结果的影 响是十分大的。之后为了证明该方法的可行性，我们对人工模拟的数据（即用正弦加噪（200 行）和指数加噪（200 行）数据及高斯噪声（1600 行）组成基因阵列），在采样点为 14 点的 情况下，用同样方法提取主模式，可以很明显的提取出正弦模式和指数模式，这证明该方法 在有噪干扰的情况下同样可以对主成分进行提取，符合基因阵列的实际情况（有噪）。但是 在对实际的基因阵列进行处理时，发现结果并不理想（主成分模式并不很有规律），原因还 要进一步探讨。最后，对模拟数据进行归类分析，通过对每个 transcriptional respone 和提 取模式求相关运算得到相关系数，然后画出散点图，可以看出散点分成三部分，和模拟数据 实际情况相符。 二． 流程 快递问题件怎么处理流程河南自建厂房流程下载关于规范招聘需求审批流程制作流程表下载邮件下载流程设计 1 生物医学信号处理特别培养计划 三． 结果及分析 1． 对人工数据进行处理： 图 6 图 7 分析：图 6 中（b）和（c）为具有规律的模式，（b）有些像正弦波，（c）有些像余弦波，但 没有指数模式的波形。（d）图则比较杂乱，是噪声模式。 2 生物医学信号处理特别培养计划 图 7 的横坐标 transcriptional respone 和图 6 中 model1（b 图中的模式）的相关系数， 纵轴为其与 model2（c 图中的模式）的相关系数。散点分为三部分，大多数点都比较零散的 分布于原点附近，说明它们和两种模式的相关性都不强，为噪声。分布在右下角的一群散点， 与 model1 相关性高，而右上角的一群仍与 model1 相关性高，不同的只是它们的相关性和 model2 有正负之分，相关程度也差不多。根据图 6 和图 7 的结果，我个人认为并不是很理 想，因为根据结果只能说明在原始数据中有两中有明显规律的模式，但只能提取出正弦模式， 并不能提取出衰减变化的指数模式，虽然图 7 表明有两群独立的点分别和两种模式相关，即 它们具有相似的变化规律，但是结果则更像是一组正弦模式和一组余弦模式，看不出有指数 模式。 2． 对生物数据进行处理 图 8 取对数并进行标准化处理的结果 图 9 未经预处理的结果 图 10 未取对数但进行标准化处理的结果 图 11 取对数减均值未除标准差的结果 分析：看出是否进行预处理对结果影响很大。取不取对数影响不大。图 8 中（b）近似看作 指数衰落的过程，（c）近似看作正弦模式，（d）没什么规律。总的来说结果并不是很理想。 四． 讨论 1．对于奇异值分解的原理的理解是是否能够很好的用其解决问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 的关键，但我对于奇异值 分解的理解是一种矩阵分解的方法，可以提取出矩阵 A 的奇异值，对于奇异值的理解是它 是矩阵 AA H 的特征值的平方根，根据矩阵谱分解的知识可知，矩阵 A 可以分解为一系 列（个数为秩）幂等矩阵的和而特征值为这些矩阵的系数，可以认为特征值的大小表征了其 3 生物医学信号处理特别培养计划 对应的分量在 A 中的比重。故用奇异值分解的方法可以提取出主成分的个数（选取合适的 阈值的前题下），但是 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 中是如何将经过奇异值分解得出的 S 和 V 阵赋予生物意义的，和 V 阵的行为什么就能表示对应的模式就不得而知了。同时，我们知道奇异值分解并不是唯一 的，即 V 阵并不是唯一的（不知 matlab 中的算法是如何选取 s 和 v 阵的）。那么在 v 阵不唯 一的情况下，还能否用其行表示基因的模式？下面有一些结果，表明用同样的数据进行奇异 值分解，结果是不同的，但之间只是存在相位的变化。为此，还有待对奇异值分解原理的探 究。 图 12 图 13 图 14 图 15 图 16 图 17 2．对图 7 进行讨论： 4 生物医学信号处理特别培养计划 观察图 7，很奇怪的发现右边两群散点和指数模式的相关性都高，而它们与正弦模式， 分别是具有相同程度的相关性，只是一个是正相关，一群是负相关。按图 7 的结果表面两种 模式之间有一定的内在联系。就此我增加采样点数，再做试验，分别取采样点为 20 点和 40 点，得到以下结果： 图 18 采样点为 20 点的主模式图 图 19 采样点为 20 点的散点图 图 20 采样点为 40 点的主模式图 图 21 采样点为 40 点的散点图 分析：由图 14 可以发现，所提取出的两个主模式之间的波形很相似，本身就具有很强的相 关性，二者之间成反相关系，这就解释了为什么散点图会如此分别。据此结果，认为该方法 无法提取出指数模式。分析原因，认为在造数据时，选取正弦模式和指数模式，并不能很理 想的提取出两种模式。而这两种模式又是根据生物数据的特性造的，所以我个人认为奇异值 分解的方法对于处理生物数据还是存在缺陷的。 图 16 和图 17 分别是选取模式为 a*sin(2*pi*t/140)，b*t.^2 和 a*sin(2*pi*t/140)，b*t.^3 时，对主成分的提取图，正弦模式仍然可以提取出来，但是对其他模式就很难看出什么规律 了。 5 生物医学信号处理特别培养计划 图 22 图 23 3． 对数据预处理的讨论 观察图 8 至图 11 可以发现是否进行预处理，如何预处理之间影响着数据处理的结果。 而数据预处理过程中将 transcriptional respone 中的行变为均值为零和方差为 1 的数据。按道 理来说，这只是将数据平移后放大或缩小倍数的过程，从本质上并没有改变数据的分布关系。 但是这却很影响结果。到底预处理有没有一个统一的方法？在结果不同的情况下，以何标准 来对结果进行选择何评价？我个人认为应该有一个对比的标准（例如用其他方法进行主成分 分析得到的结果）进行比较，这样才能够对奇异值分解的方法作出客观合理的评价。不然对 结果进行枉自推断，将其赋予生物意义，是没有太大意义的，也并不能说明该方法的好坏。 五． 附录一 1． 处理人工数据的程序 A． 生成数据矩阵子程序 sdata.m points=14; times=points*10-1; stepa=1.5/200; stepb=4/200; a1=1.5:stepa:3-stepa; b1=4:stepb:8-stepb; t1=0:10:times; f1=inline('a*sin(2*pi*t/140)','a','t'); %f2=inline('b*exp(-t/100)','b','t'); f2=inline('b*t.^2','b','t'); noise1=normrnd(0,0.5,200,points); noise=normrnd(0,0.5,1600,points); for i=1:200 ys(i,:)=f1(a1(i),t1); ye(i,:)=f2(b1(i),t1); end ys=ys+noise1; ye=ye+noise1; 6 生物医学信号处理特别培养计划 data=[ys;ye;noise]; clear i f1 f2 t1 stepa stepb noise1 noise ys ye a1 b1 B． 数据预处理子程序 spreprocess.m n=size(data); for i=1:n(1) mid=mean(data(i,:)); normdata(i,:)=(data(i,:)-mid); end clear i clear mid C． 奇异值分解子程序 sdosvd.m sum=0; [eigenassay,singularvalue,eigengene]=svd(normdata,0); %svd 的紧凑格式 eigengene=eigengene';%A=usv' for i=1:n(2) singularv(i)=singularvalue(i,i); sum=sum+singularv(i); end singularv=singularv/sum; clear sum i=1:n(2); subplot(2,2,1) bar(i,singularv) xlabel('component') ylabel('relative variance') title('(a)') clear i time=0:10:times; subplot(2,2,2) plot(time,eigengene(1,:),'o-') xlabel('time(min)') title('(b) first eigengene') subplot(2,2,3) plot(time,eigengene(2,:),'o-') xlabel('time(min)') 7 生物医学信号处理特别培养计划 title('(c) second eigengene') subplot(2,2,4) plot(time,eigengene(3,:),'o-') xlabel('time(min)') title('(d) third eigengene') D． 散点图子程序 scatter.m aa=eigengene(1:2,:)'; vv=normdata'; corrdata=[aa,vv]; corr=corrcoef(corrdata); %matlab 中是对列求相关系数 correxp=corr(1,3:n(1)+2); corrsin=corr(2,3:n(1)+2); plot(correxp,corrsin,'o') xlabel('correlation with exp model') ylabel('correlation with sin model') 2． 处理生物数据 A． 读取数据子程序 readdata.m [genes,data(:,1),data(:,2),data(:,3),data(:,4),data(:,5),data(:,6),data(:,7),data(:,8),data(:,9),data(:,10), data(:,11),data(:,12),data(:,13),data(:,14)]=textread('Select_Elutriation.txt','%s%f%f%f%f%f%f%f %f%f%f%f%f%f%f'); B． 数据预处理子程序 preprocess.m n=size(data); j=1; for i=1:n(1) test=all(data(i,:)); if test==1 %当数组中元素全为非零元素时返回 1，否则返回 0 data(i,:)=log(data(i,:)); mid=mean(data(i,:)); dev=STD(data(i,:)); normdata(j,:)=(data(i,:)-mid)./dev; j=j+1; else i=i+1; end end clear i clear mid clear dev clear j clear test 8 生物医学信号处理特别培养计划 C． 奇异值分解子程序 dosvd.m sum=0; [eigenassay,singularvalue,eigengene]=svd(normdata,0); eigengene=eigengene';%A=usv' for i=1:14 singularv(i)=singularvalue(i,i); sum=sum+singularv(i); end singularv=singularv/sum; clear sum i=1:14; subplot(2,2,1) bar(i,singularv) xlabel('component') ylabel('relative variance') title('(a)') clear i time=[0 30 60 90 120 150 180 210 240 270 300 330 360 390]; subplot(2,2,2) plot(time,eigengene(1,:),'o-') xlabel('time(min)') title('(b) first eigengene') subplot(2,2,3) plot(time,eigengene(2,:),'o-') xlabel('time(min)') title('(c) second eigengene') subplot(2,2,4) plot(time,eigengene(3,:),'o-') xlabel('time(min)') title('(d) third eigengene') 9