用奇异值分解的方法对基因芯片数据进行主成分分析

生物医学信号处理特别培养计划用奇异值分解的方法对基因芯片数据进行主成分分析（21142010袁辉）一．原理简介对基因阵列进行主成分分析，即提取主模式是一件十分有意义的工作。这里我们讨论用奇异值分解的方法对基因阵列进行模式提取。通过奇异值分解，可以根据奇异值的大小和采用的阈值提取出主要成分，画出图形，进而分析生理意义。在进行奇异值分解之前，对原始数据进行标准化预处理是十分必要的。通过下面的分析也可以看出做不做预处理对结果的影响是十分大的。之后为了证明该方法的可行性，我们对人工模拟的数据（即用正弦加噪（200行）和指数加噪（200行）数据及高斯噪声（1600行）组成基因阵列），在采样点为14点的情况下，用同样方法提取主模式，可以很明显的提取出正弦模式和指数模式，这证明该方法在有噪干扰的情况下同样可以对主成分进行提取，符合基因阵列的实际情况（有噪）。但是在对实际的基因阵列进行处理时，发现结果并不理想（主成分模式并不很有规律），原因还要进一步探讨。最后，对模拟数据进行归类分析，通过对每个transcriptionalrespone和提取模式求相关运算得到相关系数，然后画出散点图，可以看出散点分成三部分，和模拟数据实际情况相符。二．流程1生物医学信号处理特别培养计划三．结果及分析1．对人工数据进行处理：图6图7分析：图6中（b）和（c）为具有规律的模式，（b）有些像正弦波，（c）有些像余弦波，但没有指数模式的波形。（d）图则比较杂乱，是噪声模式。2生物医学信号处理特别培养计划图7的横坐标transcriptionalrespone和图6中model1（b图中的模式）的相关系数，纵轴为其与model2（c图中的模式）的相关系数。散点分为三部分，大多数点都比较零散的分布于原点附近，说明它们和两种模式的相关性都不强，为噪声。分布在右下角的一群散点，与model1相关性高，而右上角的一群仍与model1相关性高，不同的只是它们的相关性和model2有正负之分，相关程度也差不多。根据图6和图7的结果，我个人认为并不是很理想，因为根据结果只能说明在原始数据中有两中有明显规律的模式，但只能提取出正弦模式，并不能提取出衰减变化的指数模式，虽然图7表明有两群独立的点分别和两种模式相关，即它们具有相似的变化规律，但是结果则更像是一组正弦模式和一组余弦模式，看不出有指数模式。2．对生物数据进行处理图8取对数并进行标准化处理的结果图9未经预处理的结果图10未取对数但进行标准化处理的结果图11取对数减均值未除标准差的结果分析：看出是否进行预处理对结果影响很大。取不取对数影响不大。图8中（b）近似看作指数衰落的过程，（c）近似看作正弦模式，（d）没什么规律。总的来说结果并不是很理想。四．讨论1．对于奇异值分解的原理的理解是是否能够很好的用其解决问题的关键，但我对于奇异值分解的理解是一种矩阵分解的方法，可以提取出矩阵A的奇异值，对于奇异值的理解是它是矩阵AAH的特征值的平方根，根据矩阵谱分解的知识可知，矩阵A可以分解为一系列（个数为秩）幂等矩阵的和而特征值为这些矩阵的系数，可以认为特征值的大小表征了其3生物医学信号处理特别培养计划对应的分量在A中的比重。故用奇异值分解的方法可以提取出主成分的个数（选取合适的阈值的前题下），但是论文中是如何将经过奇异值分解得出的S和V阵赋予生物意义的，和V阵的行为什么就能表示对应的模式就不得而知了。同时，我们知道奇异值分解并不是唯一的，即V阵并不是唯一的（不知matlab中的算法是如何选取s和v阵的）。那么在v阵不唯一的情况下，还能否用其行表示基因的模式？下面有一些结果，表明用同样的数据进行奇异值分解，结果是不同的，但之间只是存在相位的变化。为此，还有待对奇异值分解原理的探究。图12图13图14图15图16图172．对图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图233．对数据预处理的讨论观察图8至图11可以发现是否进行预处理，如何预处理之间影响着数据处理的结果。而数据预处理过程中将transcriptionalrespone中的行变为均值为零和方差为1的数据。按道理来说，这只是将数据平移后放大或缩小倍数的过程，从本质上并没有改变数据的分布关系。但是这却很影响结果。到底预处理有没有一个统一的方法？在结果不同的情况下，以何标准来对结果进行选择何评价？我个人认为应该有一个对比的标准（例如用其他方法进行主成分分析得到的结果）进行比较，这样才能够对奇异值分解的方法作出客观合理的评价。不然对结果进行枉自推断，将其赋予生物意义，是没有太大意义的，也并不能说明该方法的好坏。五．附录一1．处理人工数据的程序A．生成数据矩阵子程序sdata.mpoints=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);fori=1:200ys(i,:)=f1(a1(i),t1);ye(i,:)=f2(b1(i),t1);endys=ys+noise1;ye=ye+noise1;6生物医学信号处理特别培养计划data=[ys;ye;noise];clearif1f2t1stepastepbnoise1noiseysyea1b1B．数据预处理子程序spreprocess.mn=size(data);fori=1:n(1)mid=mean(data(i,:));normdata(i,:)=(data(i,:)-mid);endcleariclearmidC．奇异值分解子程序sdosvd.msum=0;[eigenassay,singularvalue,eigengene]=svd(normdata,0);%svd的紧凑格式eigengene=eigengene';%A=usv'fori=1:n(2)singularv(i)=singularvalue(i,i);sum=sum+singularv(i);endsingularv=singularv/sum;clearsumi=1:n(2);subplot(2,2,1)bar(i,singularv)xlabel('component')ylabel('relativevariance')title('(a)')clearitime=0:10:times;subplot(2,2,2)plot(time,eigengene(1,:),'o-')xlabel('time(min)')title('(b)firsteigengene')subplot(2,2,3)plot(time,eigengene(2,:),'o-')xlabel('time(min)')7生物医学信号处理特别培养计划title('(c)secondeigengene')subplot(2,2,4)plot(time,eigengene(3,:),'o-')xlabel('time(min)')title('(d)thirdeigengene')D．散点图子程序scatter.maa=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('correlationwithexpmodel')ylabel('correlationwithsinmodel')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.mn=size(data);j=1;fori=1:n(1)test=all(data(i,:));iftest==1%当数组中元素全为非零元素时返回1，否则返回0data(i,:)=log(data(i,:));mid=mean(data(i,:));dev=STD(data(i,:));normdata(j,:)=(data(i,:)-mid)./dev;j=j+1;elsei=i+1;endendcleariclearmidcleardevclearjcleartest8生物医学信号处理特别培养计划C．奇异值分解子程序dosvd.msum=0;[eigenassay,singularvalue,eigengene]=svd(normdata,0);eigengene=eigengene';%A=usv'fori=1:14singularv(i)=singularvalue(i,i);sum=sum+singularv(i);endsingularv=singularv/sum;clearsumi=1:14;subplot(2,2,1)bar(i,singularv)xlabel('component')ylabel('relativevariance')title('(a)')clearitime=[0306090120150180210240270300330360390];subplot(2,2,2)plot(time,eigengene(1,:),'o-')xlabel('time(min)')title('(b)firsteigengene')subplot(2,2,3)plot(time,eigengene(2,:),'o-')xlabel('time(min)')title('(c)secondeigengene')subplot(2,2,4)plot(time,eigengene(3,:),'o-')xlabel('time(min)')title('(d)thirdeigengene')9