（概率论与数理统计专业论文）基因芯片数据归一化处理的几点研究.pdfVIP

摘要随着近代分子生物学实验技术和计算机技术的迅猛发展，以及人类基冈组草图绘制的顺利完成，标志着现代生命科学研究已经进入了后基因组时代。基因芯片作为一种新型的高通量的检测技术方法，可以同时测量成千上万个基因的表达水平，已成为后基因组时代研究基因与基因间相互作用的一一个强有力的工具。而在分析微阵列数据的早期阶段，遇到的一个基本问题就是数据的归一化。归一化的目的就是要消除基因芯片实验过程中系统变异对基因表达水平所带来的影响，它对基因芯片数据的后续分析起着重要的作用。使用一个恰当的归一化方法能保证归一化后的数据真实反映了生物学的意义。本文是在参, 照t l o w e s s 归一化方法以及h u a n g ，j e ta 1 ( 2 0 0 5 ) 提出的t w - s l m 归一一化方法的基础上，鉴于稳健性方面的考虑分别对它们作了一些改进，如针对l o w e s s 的局部加权最小二乘法，我们考虑改用局部加权最小一乘回归。而对于t w - s l m 归一化方法的改进，则是考虑运用加权最小，一l 乘回归来估计参数并把己估计出来的误差方差作为其权重系数。在后续的随机模拟部分，我们会对上述几种方法进行比较分析。关键词：l o w e s s 、t w - s l m 、微阵列、差异表达基因、半参数回归a bs t r a c tw i t ht h er a p i dd e v e l o p m e n to fm o l e c u l a rb i o l o g i c a la n dc o m p u t a t i o n a lt e c h -n o l o g i e s ，t h ei m p l e m e n to fh g pa n da c h i e v e m e n to fh g ds h o wt h a tm o d e ml i f es c i e n c er e s e a r c hh a se n t e r e dt h ep o s t - g e o m i ce r a t h ed n am i c r o a r r a yi sc a p a b l eo fp r o f i l i n gt h ee x p r e s s i o nl e v e l so fm a n yg e n e ss i m u l t a n e o u s l y , a n di sap r o m i s i n gt e c h n o l o g yf o rt h er e s e a r c ho fg e n ei n t e r a c t i o n ab a s i cq u e s t i o ni na n a l y z i n gc d n am i c r o a r r a yd a t ai sn o r m a l i z a t i o n t h ep u r p o s eo fn o r m a l i z a t i o ni st or e m o v es y s t e m a t i cb i a sw h i c he x i s tt h r o u g ht h ee x p e r i m e n ti nt h eo b s e r v e de x p r e s s i o nv a l u e sb ye s t a b l i s h i n gan o r m a l i z a t i o nc u r v ea c r o s st h ew h o l ed y n a m i cr a n g e i t 7 si m p o r t a n tf o rt h en e x td a t aa n a l y z i n g i nt h ea r t i c l e ，w eh a v em a d es o m ec h a n g sa n di m p r o v e m e n t st ot h el o w e s sn o r -m a l i z a t i o na n dt w - s l mn o r m a l i z a t i o np r o p o s e db yh u a n g ，j e ta 1 ( 2 0 0 5 ) r e s p e c t i v e l y , i nv i e wo fr o b u s t n e s s f o re x a m p l e ，w ec o n s i d e rl o c a l l yw e i g h t e dl e a s ta b s o l u td e v i a t i o nr e g r e s s i o ni n s t e a do fl o c a l l yw e i g h t e dl e a s ts q u a r e di nt h el o w e s sn o r -m a l i z a t i o n a l s o ，w eu s et h ev a r i a n c ee s t i m a t e db yt w - s l mn o r m a l i z a t i o na st h ew e i g h t sf o rt h en e x tw e i g h t e dl e a s ts q u a r e d w ea l s oc o n d u c ts i m u l a t i o ns t u d i e st oe v a l u a t et h ea b o v e m e n t i o n e dn o r m a l i z a t i o nm e t h o d s k e yw o r d s ：l o w e s st w - s l md i f f e r e n t i a l l ye x p r e s s i o n e dg e n e sm i c r o a r r a ys e m i p a r a m e t r i cr e g r e s s i o n第1 章引言随着生命科学进入“后基因组时代”，基因组研究的重心也逐渐转向了基因功能的研究。基因：卷片技术无疑为基因功能研究提供了一种强有力的工具，高通量特点( 该技术可在同一一块芯片上同时检测出成千上万个基因的农达水平) 使其在基因表达分析、疾病诊断和治疗、新药发现等众多领域得到广泛应用。但是如何对基因芯片应用过程中所产生的关系复杂的海量数据进行准确且合理的分析，并从中挖掘出有意义的生物信息已成为限制该技术进一步发展的主要“瓶颁”，也是当前生物信息学的主要研究内容和主要发展方向。本文主要就基因表达数据前期处理过程中数据归一化这一环节进行分析和研究。( 伍业舟，2 0 0 4 )1 1归一化方法的背景介绍微阵列技术已成为高通量检测基因表达模式的强有力工具并广泛运用在生物遗传学领域。c d n a 芯片或微阵列研究的基本步骤是：从实验样本( t e s ts a m p l e ) 和对照样本( c o n t r o ls a m p l e ) 中等量提取r n a ，随后进行反转录为c d n a ，并用两种不i 一颜色的荧光探针c y 3 ( 绿色) 、c y 5 ( 宴 3 2 色) 进行标记，之后将标记了的c d n a 混合杂交到制备好的c d n a 芯片或微阵列上，待充分杂交后，进行清洗，然后用带有两种不同波长的激光扫描仪对芯片进行扫描处理，经过扫描后的芯片可产生两张无符号的1 6 位字节的图像文件，最后对其进行图像处理和数据分析。我们所要研究的样品中某个基因的表达水平是在微阵列实验中通过实验样本和对照样本的浓度之比来确定的，而浓度之比表现为结合在芯j 卜的该基凶所对应的探针上的两种来源的c d n a 的多少，b p c y 3 和c y 5 的荧光强度之比。而实际微阵列实验巾测量到的c v 3 和c y 5 f l 勺荧光强度值受许多因素的影响，这些因素的影响会造成测量值的变异，这些变异包括系统变异和随机变异。随机变异由随机因素引起，是无法控制12的因素。引起系统变异的因素包括样本的生物学变异和系统误差，而研究人员的研究目的在于了解样本的牛物学差异，因此必须把系统变异中的系统误差部分消除，由此也就产生了数据的归一化处理。归一化实质上是通过在整个表达值范围内建立个归一化曲线( 或称归一化因子) 后，用原数据减去这个归一化因子，由此去消除实验数据中的系统误差。使用一个恰当的归一化方法，能保证归一化后的浓度比值反映的是真实的有生物学意义的基因相对表达水平。在微阵列实验中引起系统误差的因素有很多，包括两种荧光标记物的物理和化学属性的差异、扫描仪的属性设置对荧光强度测量的影响以及两种荧光标记物在标记效率上的差异、不同点样头问的差异、芯片的空问位置的不同等因素，这些因素对荧光强度或强度比的影响程度的差异将造成系统误差，使测量的荧光强度或强度比不能真实反映样本的生物学差异。所以在进行后续的数据统计分析( 如鉴别差异表达基冈的显著性分析和基冈表达模式的聚类等) 之前，对数据进行归一化处理是非常有必要而且是重要的。对于常见的c d n a 微阵列芯片。让r 和g 分别表示芯片上某个点除去了背景荧光强度之后的c y 5 ( r ) 和c y 3 ( g ) 的荧光强度值，它们代表的其实就是实验样本和对照样本中某个基因的表达水平。考虑到生物学上的意义，我们经常是对r 和g 作一个对数变换：令z = l 0 9 2 佩，y = l 0 9 2 暑。微阵列中系统偏倚是否存在、系统偏倚的模式以及对归一化效果的评价，都需要采用适当的方法对数据中的系统误差进行呈现，最为简单、便捷和直观的方法为统计图，我们可以以z 为横坐标，为纵坐标画出它们的散点图。如果数据不存在明显的系统偏倚，那么散点图上的点大致应该关于直线= 0 对称。从下面左图可以看到，在任何强度下，c y 5 的荧光强度都系统地高于c v 3 的荧光强度，若由生物学角度可以判断差异表达基因的比例很小，则数据存在明显的偏倚。从图形分布来看，这种偏倚不存在点的荧光强度的依存性。而下面右图中的系统误差则存在明显的强度依存性。( 参考李瑶，2 0 0 6 )浙江大学硕士学位论文3归一一化方法实际上就是要找到归一化冈子或者说是归一化曲线，然后用减去归一化凶了。y a n ge ta 1 ( 2 0 0 1 ) 系统地研究和分析了多种归一化方法，比如全局、强度依存偏倚以及染色互换等归一化方法。其中全局归一化方法需要假设芯片上对于所有基因来说归一一化因子是个常数并且红绿荧光强度是通过常数因子相关联的，兄= 七g 。因此全局归一化后的结果就是把强度对数比的均数归一化到0 对于有强度依存偏倚的归一化方法，y a n ge ta 1 提m 了在关于z ，的散点图中采用局部加权线性平滑( c l e v e l a n d1 9 7 9 ) ，并且把平滑后的残差作为归一化后的强度比对数。k e r re ta 1 ( 2 0 0 0 ) 提出的方差分析方法( a n o v a ) 和w 6 l f i n g e re ta 1 ( 2 0 0 1 ) 提出的混合线性模型则是把芯片以及染色的效应都考虑进去并假设常数归一化因子。而f a n , j q e ta 1 ( 2 0 0 4 ，2 0 0 5 ，2 0 0 6 ) 讨论的是利用单张芯片内的重复点样，存单张芯片内建立半线性模型( s l i m ) 。s l i m 方法要求芯片上有重复点样的基因。如果重复点样的基因个数较少，那么这些重复点样的基冈的表达水平就不能涵盖这块芯片上所有基冈的表达水平或是反映出它的空间差异。h u a n g ，j e ta 1 ( 2 0 0 5 ) 提出了一种t w - s l m 的归二化方法，他们是从半参数回归模型的角度重新考察了l o w e s s 归一化方法后得山的一种两向半线性模型。m a ，s c h u m a k e r e ta 1 ( 2 0 0 6 ) 从稳健性出发提出把加权最小绝对偏差回归和加权b o o t s t r a p 方法运用在模型s l s m 中得出了s l s m 归一化方泫。星jx = l 2 ( i o g r h o g g )x = l r 2 ( k j g r h o g g )浙江大学硕士学位论文41 2 本文提出的归一化方法的简单介绍全目前为止，已被提出了很多用于基冈芯片数据处理的归一化方法。本文主要是在基于l o w e s s 归一化方法和t w - s l m 归一化方法研究的基础上考虑其稳健性易受到异常点的影响而作了一些改进。对十运用最为广泛的l o w e s s 归一化方法，它需要至少满足以一f - - - 个生物学上的假设：( 1 ) 芯片上有差异表达基因所占的比例要很小，或者( 2 ) 在表达水平上上调表达的基冈与下调表达的基因呈对称性。l o w e s s 归一化方法实际卜就是采用了局部加权线性回归去平滑每个点，当然在c l e v e l a n d ( 1 9 7 9 ) 中提及的稳健局部加权回归运用在归一化以后，可以使l o w e s s 归一化方法对异常点更具有稳健性，其中它考虑的是运用加权最小二乘，并通过改变第二步的权重系数，叠代循环米达到稳健性，而这样就增大了计算量。在统计计算方面我们都知道最小绝对偏差( 或称最小一乘) 比最小_ 乘更具有稳健性更不易受异常点的影响。因此在本文的第二章我们将会讨论把局部加权最小一一乘| 日1 归运用在数据归。+ 化中，它是对l o w e s s 归一化方法的一种改进。h u a n g ，j e ta 1 ( 2 0 0 5 ) 提出把t w - s l m 运用在c d n a 微阵列数据的归一化中，其实是从半参数回归模型的角度重新考察了l o w e s s 归一化方法后得出的一种两向半线性模型归一化方法。它像一般的运用模型进行归一化的方法一样，不需要满足l o w e s s 归一化方法所要满足的生物学假设：( 1 ) 芯片上差异表达基因所占的比例要很小，或者( 2 ) 在表达水平上上调表达的基因与下调表达的基因呈对称性。并且它的模型假设也是很自然会满足的，因此我们在第三章提出的归一化方法中仍采用t w - s l m 模型，在原模型中参数估计采用的是最小二乘估计，非参数部分用的是b 样条估计，在本文中我们会沿用t w - s l m 的各种模型假设条件，只是在参数估计部分，会作出一些改进。考虑到稳健性方面，我们把最小二乘估计改为加权最小二乘估计，用估计出的方差的倒数作为其权重系数。我们将在第四章中利用h u a n g ，j e ta 1 ( 2 0 0 5 ) 中提及的网个模型通过随机模拟对卜述几种归一化方法进行比较和评价。第2 章l o w e s s 归一化方法和它的改进2 1l o w e s s 9 3 - - 化方法的介绍假设在实验中或研究对象中有n 块芯片、j 个基因，并且每块芯片上的每个基因都只点样一次。让尼j 和g 巧分别表示来自实验样本和对照样本里的第i 块芯片上第歹个基因的荧光强度值，按照c h e ne ta 1 ( 1 9 9 7 ) 和y a n ge ta 1 ( 2 0 0 1 ) 的做法，令表示第i 块芯片上第j 个基因的荧光强度比对数，z i j 表示的是对数荧光强度的平均值。也即有下式：y i j = l 0 9 2 瓦r i j ，= 互1l 。9 2 ( 嘞g 巧) ，江l ，n ，j = l ，z( 2 1 )l o w e s s 归一化方法是分别对每块芯片进行的，那么对于第i 块芯片( i = 1 ，n ) ，l o w e s s 归一化方法所要作的就是去甲滑拟合下面的非参数回归模型：= f i ( x i j ) + 0 ，j = 1 ，( 2 2 )这里假设六是平滑函数，e 0 是随机误差。利用c l e v e l a n d 的l o w e s s 方法后，可以算出 的l o w e s s 估计值五，之后可以通过上述模型计算出它的残差誉0 = 一五( z 玎) ，i = 1 ，n ，j = 1 ，z这些残差就被称为归一化后的数据，并被用于后续的统计分析中。所以一般来说，整个分析过程包含二步：( 1 ) 归一化；以及( 2 ) 基于归一化后的数据芘进行统计分析。比如说，在一个直接比较两个c d n a 样本的实验设计中( 即两个c d n a 样本同时杂交在一块芯片上) ，一个典型的分析方法就是第+ 一步运用l o w e s s 方法把数据归。化得到归一化后的数据，之后再基于这些数据去鉴别有差异表达的基因。在直接比较实验设计中卜述的j ：步分析方法在统计框架里可以看成足使用了两个模型。第一个就是在( 2 2 ) 中5浙江大学硕士学位论文给出的用于归一化的非参数回归模型。第二个考虑残差模型：吒= | 1 3 j + ( 2 3 )这里的岛表示的是基因j 的相对表达水平。显著性分析的目的就是要找到那些岛o 的基冈。( 2 2 ) 式和( 2 3 ) 式在这个两步分析方法中是分开讨论的两个模型，并且在分析中认为岛e 易，忽略掉它的近似误差。下面将说h ) j l o w e s s ：h - 法中是如何估计 的，也即五是如何计算得出的。参照c l e v e l a n d ( 1 9 7 9 ) 中的方法，先给出文中所用的“t r i c u b e ”权函数：w c z ，= 一l z l 3 ) 3 ：三：三：1 当i z l 0 ；3 当z o 时，( z ) 是一个非增函数；4 当1 时，w ( x ) = 0 点( ，协) ，j = 1 ，j ( 下标i = 1 省略) 。按照模型( 2 2 ) 有下式成立：协= 厂( 巧) + 弓，j = 1 ，z要想得至lj f ( x j ) 的估计值坊，也就是要去平滑点( ，协) ，由于厂光滑性的假设，因此( 巧，协) 领域内附近的点可以帮助估计彩，而它们的权重应该与它们的横坐标到巧的距离成反比。让o b 1 ，r = b j 】，增人b ，会使得用于估计奶的领域内的点增多，并且会增加平滑点的光滑性。下面给出详细的过程，对于每个基因j ，令允f 表示横坐标上z j 领域范围内离它第r i f t 的点与它的距离。w k ( x j ) 表示的是第七个点( z 缸，弧) 用于估计奶的权重系数，它由7下式给出：毗) 划( 半= 1 。以下步骤用于计算，( ) 的估计值奶，它包括局部加权回归和稳健局部加权回归：第1 步对于每一们，运用加权最小二乘法( 叫七( 巧) 作为点( z 知，纨) 的权重系数) 计算d 维多项式参数的估计值反( 巧) ，i = 0 ，d 也即：j = 。r g r o i n 叫南( 巧) ( 玑一风一f l l x k 一阮z 2 ) 2 屉= l那么对于局部加权回归而言，( 巧) 的估计值奶就可以写成：第2 步先给出b i s q u a r e 权函数b ：b c z ，= 一z 2 ) 2 ：三：三：令：弓= 协一易表示的是通过第一步之后计算出的残差。s 表示为所有i 弓i ，( j =瓯= b ( e * 七6 s ) 第3 步把第一步中用于计算厦( 巧) 的权重系数换成如( ) ，重新利用局部加权最小二乘回归法计算锄第4 步重复第2 第3 步t 次，最后算出的雪j 就是稳健局部加权回归的平滑值。在l o w e s s 方法中是通过第2 4 步来减小异常点的影响，使其方法具有稳健性。而我们都知道在统计领域里最小一乘法比最d x _ c 乘法更具有稳健性。因此本文基于这种想法所提出的稳健性改进，则是把第一步中的局部加权最小j - 乘回归改为局部加权最小一乘回归，我们将在下一节讨论研究。zz陵d 铷=如82 2对l o w e s s 9 3 - - 化方法的改进在这里我们沿用l o w e s s 归一化方法中所提及的权函数和权重系数，并考虑线性多项式啊归( 此时d = 1 ) 。| 一样考虑非参数回归模型：协= ，( q ) + 歹= l ，对于每一个j ，可以通过下列步骤计算，( 巧) 的估计值坊：1 用局部加权最小一乘回归计算线性多项式参数阮( ) ，隗( ) 。 = n r g r o i n 叫彪( 巧) i 玑一岛一岛z 七2 计算鲂= 阮( ) + f l l ( x j ) x j 。其中在第一步中关于岛( q ) ，庑( ) 的计算可以采用循环叠代的方法，如我们可以令初始值硝o ) = 0 对于第i 次叠代( i = 0 ，1 ，2 ，)( i ) 通过使目标函数壹叫南( ) i 玑一肺一硝z 詹l 达到最小来计算硝七= 1( i i ) 利用前一步已经计算出米的硝，通过使目标函数妻后( 巧) l y 七一f 1 0 ( i 一f l l x 七i 达到最小来计算p 仆。重复( i ) ( i i ) 直到满足给出的收敛条件，假设它在第k 次达到收敛条件，则岛=靥刖，麂= 硝。第3 章n 怛s l m 归一化方法和它的改进30 1t w - s l m 归一化方法的介绍h u a n g ，j ( 2 0 0 5 ) 提出把两向半线性模型( n - s l m ) 用于c d n a 微阵列数据的归一化，是基于从半参数回归模犁的角度对l o w e s s 归一化方法的重新思考。下面我们先给出半参数i 司归模型的简单介绍。3 1 1半参数回归模型假没有一批数据含有礼组观测值( 玑，鼢，乞) ，i = 1 ，n ，其中玑代表的是响应变量，- 币i ( x i ，忍) 代表的是共变量。那么半参数回归模型( s r m ) - - 般可以表示成下式：轨= f ( x t ) + 旎7 p + 旬，i = l ，n ( 3 1 )上式中的厂是一个未知函数，p 是回归参数，岛代表的是误差。上述模型在很多场合都可适用。比如说，z i 常常可以代表两种不同条件f 的两分变量( 如接受治疗与服用安慰剂这两种条件) ，并且我们所关心的只是治疗的效席卢，但是又需要协调好连续性协变量的效应。对于一个p 维的协变量甄= ( x i l ，z 咖) 7 ，常常是把，当作可加的结构来处理( h a s t i ea n dt i b s h i r a n i1 9 9 0 ) 。因此另一种可加形式的半参数回归模型可表示为：y i = f l ( x 1 ) + + 。厂p ( z 咖) + 名7 p + e i ，i = 1 ，佗( 3 2 )模型( 3 1 ) 和( 3 2 ) 是两个基本的半参数回归模型。对于模型中的参数估计部分常常有两点是需要注意的。第一，和是需要同时估计的。像下面的做法都是不合适的，例如固定p 为o ，去得- 至l j f 的估计，之后把这个，当成已知量代入到模型( 3 1 ) ，去估计p 。第二，按照半参数的信息理论，在估计p 时，需要把估计，时的不确定性考虑进去。具体的例子可以参见( b i c k e le ta 1 1 9 9 3 ) 。91 053 1 2两向半线性模型( t w - s l m )我们先通过特例( 直接比较实验) 来描述一下什么是两向半线性模型( t w - s l m ) 。在这个实验设计中，把两个标有不同荧光标记的来自不同细胞母体的c d n a 样本同时杂交到一块芯片上。让和z 巧定义为( 2 1 ) 式中的形式。则对于这个实验设计，t w - s l m 可表示为下式：y i j = 五( z 巧) + 岛+ 巧，i = 1 ，n ，j = 1 ，( 3 3 )这里的五表示的是第i 块芯片上有强度依存性的归一化曲线。风r 代表的是归一化后基凶歹的相对表达水平值。而鼠，是随机误差，均值为0 ，方差为盯磊。模型t w - s l m 可以看成是f l j l o w e s s 归一化方法中( 2 2 ) 式和( 2 3 ) 式结合而成的，把( 2 3 ) 式替换到( 2 2 ) 式里就可以得到( 3 3 ) 式。把( 2 2 ) 式和( 2 3 ) 式这两个模型结合起来考虑可以使我们同时对归一化曲线和基因效应作估计，而这样也是合理的，因为在作估计之前是并不清楚哪些基因是稳定表达的( 即岛= 0 ) ，倘若已经识别了大批稳定表达的基阗，利用这些基冈当然是可以得到渐近无偏的归一化，但实际上这种情况却是很少会发生的。在h u a n g ，j e ta 1 ( 2 0 0 5 ) 的文章中提到过之所以把( 3 3 ) 称为两向模型，是冈为它可以看成是两向方著模型半参数性质的推广。也就是说当五= q i ，i = 1 ，凡( q i 是常数参数) 时，( 3 3 ) 就变成了一个两向方差模型。另方面，t w - s l m 可以看成是标准的半参数【旦j 归模型( 3 1 ) 的延伸和扩充但又不同于( 3 1 ) ，当然它也不同于( 3 2 ) 。尤其是，在( 3 1 ) 和( 3 2 ) 中有限和无限维的参数个数是确定的、不受样本量的约束，并且他们也不包含两向方差分析模型作为子模型。相反，在模型1 、-s l m 中，有限维参数( 角，胁) 的个数是j ，它正好是估计五的样本容量，而无限维参数( ( z ) ，厶( z ) ) 的个数是几，它却是用来估计岛的样本容量。一般情况下，令z i r d 是对应着第i 块芯片的共变向量。那么模型t w - s l m 的一浙江大学硕士学位论文1 1般形式可以写成：= 六( z 巧) + 么7 岛十巧，i = 1 ，n ，j = 1 ，( 3 4 )这里的岛r d 是第j 个基因的效应，五与的含义与( 3 3 ) 中一样。共变向量z i 可以用来编译各种实验设计( k e r r , m k e ta 1 2 0 0 1 ) 。对于两个样本的直接比较实验，此时z i = l ，i = l ，n ，模型( 3 4 ) 就变成模型( 3 3 ) 。而对于借助参照样本的问接比较实验设计，则可以引进。4 个两维的共变向量z i = ( z i ，z i 2 ) 7 。当第i 块芯片是用样本一与参照样本比较时令兹= ( 1 ，0 ) 7 ，当第i 块芯片是用样本二与参照样本比较时令忍= ( 0 ，1 ) 7 。此时传= ( 岛1 ，3 j 2 ) 7 是一个两维向最，岛1 一岛z 代表的是归一化后基因j 表达水平的差值。在模型( 3 4 ) 中，很直接就能看出归一化曲线五足依赖于某个特定芯片的，如果假设五也依赖于芯片上某个点样区组，则也很容易把模型( 3 4 ) 进行推广。具体地，假设在每块芯片上有k 个点样区组，在第庇个区组上有以个基因被点样。令犰幻和锄幻分别表示第i 块芯片上第k 个区组内基因j 的荧光强度比对数和平均对数荧光强度。则模型t w - s l m 可写为：耽幻= f i k ( x i k j ) 十乏岛+ e t 幻，( 3 5 )i = 1 ，n ，k = 1 ，k ，j = 1 ，以t w - s l m 也可以推广到当芯片上有基因被重复点样的实验设计模型，一般来说，这种实验设计可以帮助提高精确性，并通过变异系数( c v ) 来评判芯片的质n ( t s e n ge ta 1 2 0 0 1 ) 。假设在第i 块芯片的第k 个区组内有以个基因被点样，其中基因j 在这个区组内被重复点样尼他次。则模型t w - s l m 可写成：y i k j ，= 知( z t 幻，) + 磊岛+ 岛幻，( 3 6 )i = l ，n ，k = 1 ，k ，j = 1 ，以，r = 1 ，尼j 七浙江大学硕上学位论文1 23 1 3t w - s l m l ) 3 - - 化方法h u a n g ，j e ta 1 ( 2 0 0 5 ) 提出用模型t w - s l m 来解决c d n a 微阵列芯片数据的归一化问题，实际上就是要对模型中的归一化曲线以及基因表达水平值进行估计，而他们采用的是半参数最小：一：乘估计( s l s e ) 并结合多项式样条的方法( s c h u m a k e r1 9 8 1 ) 来实现的。下面就来看看模型t w - s l m 中是如何进行半参数最d , - 乘估计的。令q j d 是由所有了d 维矩阵p 三( p 1 ，励) 7 开、枘, 4 4 - 疋乙g 任1 岛= 0 所组成的空间。从模型( 3 4 ) n 怛s l m 的定义中可以看出，p 能够满足窄问q 苫d 内的条件。因此可以不妨假设：，p q 三 p ：岛= o )( 3 7 )j = l倘若一。开始不满足名，岛= 0 ，则用传一冬，仇l 厂去代替岛，用五( z ) 十岛j ：，臃，去代替 ( z ) 。假设b m ? b i k , 是b 样条基底函数。并设：& 三 b i o ( x ) 三1 ，b i k ( x ) ，k = 1 ，k )( 3 8 )是由这些基底函数的所有线性组合构成的空间。并用下式来近似估计 ：风q i 。+ q 他玩七( z ) 兰b i ( x ) 7 q t &七= 1这里b t ( z ) = ( 1 ，b i l ( z ) ，b i k e ( x ) ) 7 ，而啦= ( q 0 ，c e i l ，c l i k i ) 7 是需要从数据中被估计的系数。令t 厂= ( ，厶) 和p = ( 角，励) 7 。则最小二乘( l s ) 目标函数是：njd 2 ( p ，) = 一五( z 巧) 一磊岛】2 i = lj = lh u a n g ，j e ta 1 ( 2 0 0 5 ) 提出的关于 p ， 的半参数最小二乘估计是q g x d n 鍪1 & 中使目标函数达到最小的 p ， 。也就是：( p ，) = a r gm ，i n d 2 ( p ，)( 3 9 )( p ，) q ：“n 饕1 &浙江大学硕士学位论文令嘞= ( 1 ，b 订( x i j ) ，b i k ；( x i j ) ) 7 表示的是样条基函数在z 玎，1sis 佗，l j ，处的取值。则第i 块芯”的样条基矩阵是：i 毽j耻l 、lrh现1 ( 既1 ) 玩k ( x i lb t l ( x i j ) b i k i ( x i j如果令0 ：= ( n 1 i q 。) 7 ，玑= ( 玑l ，y i g ) ，则d 2 ( p ，o l ) = d 2 ( p ，厂) 。考虑( p ，q ) 的估计，使d 2 ( p ，a ) 达到最小的问题等价于求解下面的线性方程：卢z i z i7 ) + b i & t 旎7 = y i z i 7 ，最如十矽乞= 玑i = 1i = 1t = 1让( p ，a ) 代表上述方程的解，且定义五( z ) 兰b i ( z ) 龟，i = 1 ，n 注：用目标函数d 2 ( 卢，a ) 分别对p 和o i ，i = 1 ，n 求微商，再令微商等于0 ，经过化简就可以得到上述的线性方程。虽然可以通过求解上述线性方程来估计p 和o ，但在实际计算中是利用g a u s s s e i d e l 方法( h a s t i ee ta 1 2 0 0 1 ) 来最小化目标函数，它是通过交替更新q 和p 来实现的。下面是它的计算步骤，令初始值卢( o ) = 0 ，对于第k 次叠代( 尼= 0 ，1 ，2 ，) ：第一步：求解q ( 七) 使d 2 ( p ( “，a ) 达到最小。解的具体形式是：q = ( b 鼠) 一1b ：( 玑一p ( 七) 乞) ，i = 1 ，n ( 3 1 0 )第二步：利用第一步已经解出来的q ( 七) ，令一南( z ) = b i ( z ) 7 q ：舶，求解p ( 岛+ 1 ) 使d 2 ( p ，a ( 七) ) 达到最小。解的具体形式是：( 3 1 1 )重复第一第二步直到满足给定的收敛条件。由于目标函数是严格凸的，因此算法一定收敛。假设上述算法在第k 步达到收敛准则。则岛的估计值是岛= f l ( j k ) , j =1 ，且估计的归一化曲线是：五( z ) = b i ( z ) 7 q ：k = b i ( z ) 7 ( 鼠) 一1 e ( 玑一声乙) ，i = 1 ，扎( 3 1 2 )zli |抚砷矗一兹n 曲一，魂。汹=d+岛烈令x = ( 兢，x i d ) ， ( x 。) = ( 五( z i l ) ，五( z i ，) ) 7 以及q = 鼠( e b d e 。通过( 3 1 2 ) ，则五( z i ) 的估计为：五( 鼢) = q i ( 犰一f l z d 因此归一化曲线可以看成是线性平滑算子q i 作用在玑一声乞上的结果，从玑中去除了基因的效应厉旎。相反l o w e s s 归一化方法是并没有去除基因效应的。以上就是t w - s l m 归化方法的过程，它不仅估计了归，化曲线( 这是归一化的目的) 也同时估计了基因的效应。这样就能把归一化过程的误差考虑进后续的显著性分析过程中，h u a n g ，j e ta 1 ( 2 0 0 5 ) 中还讨论了如何利用1 、- s l m 进行显著性分析，这里我们就不讨论了。3 2t w - s l m l ) 3 - - 化方法稳健性的改进在t w 二s l m 归一化方法巾的参数估计部分作者考虑的是用半参数最小二乘估计法来估计的，而一般来说最小二乘估计对于异常点不太具有稳健性，容易受到异常点的干扰，而对于c d n a 微阵列实验，影响实验结果的凶素有很多，再加上若差异表达基因的数目增多的话，则产生的数据往往会有很多异常点，这时n - s l m 归一化方法的稳健性就会受到影响。本节就是出于稳健性方面的考虑对t w - s l m 归一化方法中的参数估计部分作了相应的改进。w a n g , d l e ta 1 ( 2 0 0 5 ) 也是出于此目的提出了稳健性的t w - s l m 归一化方法。而我们考虑用半参数加权最小二二乘估计法代替最小二二乘估计以提高稳健性，我们借a ，s h u a n g g ee ta 1 ( 2 0 0 6 ) 所提出的归一化方法中用方差的倒数作为权重系数的考虑，在他的文章中提出的模型是s l s m ，包含模型( 3 3 ) 作为它的特例，并提出用加权最小绝对偏差法来作估计。当然也提到过可用加权最小二乘法来作参数估计。而本文所提出的加权最小二乘估计与m a ，s h u a n g g ee ta 1 ( 2 0 0 6 ) 中的有所不同，他们是把方差也同时作为估计的部分，这样大大增加了计算量。而我们所提出的归一。化方法实际上是分了两步。第一步：按照t w 二s l m 方法作参数估计，之后利用残差对误差方差作估计。第二步：我们利用第一步已经估计出的方差再对参数作加权最d , - 乘估计。此时我们的权重系数已知，这样既能达到稳健性又能减少计算量。下面将具体说明本文所提出的对t w - s l m 归一化方法的改进。在这里我们以模型t w - s l m 中的特例( 3 3 ) 为例加以说明，此时d = 1 ，实验设计为两样本的直接比较实验。模型为f 式：y i j = ( z 巧) + 岛+ e 订，i = 1 ，n ，j = 1 ，上式的各项说明如z 前t w - s l m 归一化方法中所1 段设的，我们所矍作的就是估计归一化曲线五和基因效应岛。第一步：利用1 、- s l m 归一化方法中已经估计出的j i ：f 和岛对误差方差吒t f 估h - 。我们先计算估计后的残差毛：毒弓= ( 一五( z 巧) 一岛) 2 ，i = l ，n ，j = 1 ，z对于如何估计吒2 j 可以考虑两种模型。( i ) 假设误差方差对于每个基因是不同的，但是不会随着基因芯片的变化而改变。也就是说，对于j = l ，z喝= 弓，i = 1 ，佗这时我们可以这样水估计窍：豸= i - 1e i 但是对于这种估计方著的做法有一个问题存在，因为在微阵列数据研究中基因的数目往往是非常巨大的，因此即使基因表达差异非常小，许多方差的估计值霹也非常小，但是偶然也会导致统计量亡( 用于显著性分析) 的值增大。一个解决的办法就是给莺增加一个适当的常数。( t u s h e r , vg e ta 1 2 0 0 1 ) 。伍) 如果像归一化曲线一样，假设误差方差存在强度依存性，并且依存性随着：卷片的变化向变化。则我们可以假设误差方差满足如下模犁：盯弓= 砰( z 巧) ，i = 1 ，几，j = 1 ，z浙江大学硕士学位论文1 6这里的盯；是光滑的正函数，这个模型把芯片的效应也考虑到方差里。由于模型中一光滑性的假设，也就是说在每个芯片上具有相似荧光强度表达水甲值的基凶他们的误差方差也相同，这个假设也是合理的，因为对于许多微阵列数据，强度比对数的变异性依赖于整个的荧光强度值。尤其表现为强度低的变异性比强度高的变异性大。我们采用非参数回归模型中平滑散点的办法对散点图( ，铭) 拟合光滑曲线。相当于拟合下列的非参数回归模型：誉弓= a ；( x i j ) + t i d ，j = l ，z我们使用( 3 1 2 ) 0 7 用来估计五的样条基矩阵，因此类似于五的估计，鳄的样条估计可以表示为：a y ( x ) = b ：( z ) ( b ：鼠) _ 1 鼠箐，这里管= ( 晶，拿易) 7 ，因此吗的估计为铭= 霹( z 巧) 。第二步：利用第一步已经估计出的误差方差铭，再次对模型t w - s l m 中的五和岛进行加权最小j ：乘估计。这时我们要考虑的加权最小二乘估计的目标函数为下式：nj1d 2 ( p ，) = 翱一五( z 巧) 一z j 2一1i 一1 1 1如果写成d 2 ( 口，q ) 就变为：令：则得到线性方程：d 2 ( 卢，q )n( 嘉) 岛l = 1 。”一b ( z 巧) 7 q i 一岛】2( 一q t 一岛) 2嘞( 一a t 一伤) = 0( 一q t 一伤) = 0= 耋1 呼 =”。巧。巧j触j似n甜n：亘。巧。虿j博。汹一生警酱俦一。虿j 触=皿。虿j 芦( p ，q ) 的估计值即为上述线性方程的解( 矽，a ) 。且定义五( z ) 三b i ( z ) 7 也在实际计算时，我们同样采用g a u s s s e i d e l y y 法通过交替更新q 和矽来求解。计算步骤如f ：同样令初始值p ( o ) = 0 ，对于第七次叠代( 后= 0 ，1 ，2 ，) ：步骤一：求解q ( 七) ，使d 2 ( p ( 舢，q ) 达到最小。解的具体形式为：q = c ；j 去嘞) _ 1c 骞南嘞c 一碍q ni = 1 棚步骤二：利用步骤一中求出的q ( 詹) 。令一七( z ) = b i ( z ) n ：m ，通过使d 2 ( p ，q ( 詹) 达到最小求出p ( 七+ 1 ) 。解的具体形式为：秽=(妻i=i扩妻i皆一，z=i、3、3重复步骤一和步骤二直到满足给定的收敛条件。假设第k 次后收敛，则：岛的估计值是岛= 磅刖，j = 1 ，z 并且估计的归一化曲线是：五( z ) = b 舡舳：k ) = b 舡) ，( ；j 专) _ 1 ( ；j 专嘞( 一岛) ) i = - 以通过第二步我们就得剑了参数的加权最小二乘估计值房和五( z ) 。当然我们可以通过它再次估计误差方差然后再返回到第二步进行计算，重复几次后得到最终的估计值。第4 章随机模拟这一章我们主要利用随机模拟实验来比较评判上述几种归一化方法。我们考虑采用均方误差( m s e ) 作为评判的标准。令a 1 和a 2 分别表示卜调和下调基因所占的百分比，并令a = a 1 + a 2 。在随机模拟中我们同样考虑以下四种模型( h u a n g , j e ta 1 ( 2 0 0 5 ) 中提到的) 。模型1 ：假设没有染色偏倚。冈此真正的归一化曲线水平为0 即：五( z ) 兰0 ，1si n 。另外，上调和下调表达基凶是均衡的，也即入1 = a 2 。模型2 ：像模型1 一样，归一化曲线 ( z ) 三0 ，1si n 。但是上调基因与下调基因所占百分比不同，我们设a 1 = 3 a 2 。模型3 ：假设存在非线性的有强度依存性的染色偏倚。上调和下调表达基因是均衡的，a 1 = a 2 。模型4 ：假设存在非线性的有强度依存性的染色偏倚。上调基因与下调基因所占百分比不同，并设入1 = 3 a 2 。( i ) 产生岛，对于大多数基因来说我们假设岛n ( o ，亏) 。并且他们所占的百分比是1 一入。对于上调基凶，我们假设岛一( p ，4 j ) ，这里p 0 。对于下调基因，我们设伤一( 一p ，弓j ) 。并且在模拟中我们设勺= 0 6 ，肛= 2 ，t d j = r v j = 1 。( i i ) 产生z 巧，我们假设z 巧一1 6 木b e t a ( 1 ，2 5 ) 。( i i i ) 产生，我们假设n ( 0 ，吒2 ) 。其中= 盯( ) ，盯( z ) = 0 3 水z - 1 3 。因此强度低的误差方差比强度高的要大。( i v ) 强度比对数值就可以如下计算：= f i ( x i j ) + 岛+ 。在模型3 和模型4 中，对于第i 块芯片上的归一化曲线，我们假设：施) = 咩爰字艇 0 , 1 6 1上式中的a i l 和a 1 2 是相互独立地从均匀分布u ( o 6 ，1 4 ) 中产生的。因此归一化曲线随着芯片的变动而变化。1 8浙江大学硕士学位论文1 9在模拟实验中，假设每次数据包含1 0 块芯片，每块芯片上含有3 0 0 0 个基因，每次模拟重复1 0 次，我们依据这1 0 次的结果计算估计值的均方误差。在上述的每个模型里我们考虑差异表达基凶所占百分比数的三种水平：a = o 0 1 ，0 0 6 ，0 1 2 。下面的四个表格列出了三种归一化方法分别在上述四个随机模拟模型下基因相对表达水平岛的估计值的均方误差的统计描述。下表中的t m 代表h u a n g ，j e ta 1 ( 2 0 0 5 ) 提出的t w - s l m 归一化方法，m t m f 表本文提出的对t w - s l m 改进后的归一化方法。表一列的是第一个模型下的结果，可以看到t m 和l o w e s s 归一化方法差不多，而m t m 归一化方法下的均方误差( m s e ) 去p - - 致的比另两种要小，特别是当存在非线性的染色偏倚，而且上调基因与下调基因所占百分比不同时，l o w e s s 归一化方法就显得不太好，并且可以看到在给出的四个模型下，改进后的m t m 归+ 化方法下的均方误差比1 、- s l m 归一化下的均方误差有了明显的减小，从图2 的箱线图中也可看至- i j m t m 方法比n - s l m 方法更稳健，而这也是我们的目的。本文中的所有模拟