杨 敏，张 静

(云南大学数学与统计学院统计系，昆明650091)

基因表达是指基因在生物体内的转录、剪接、翻译以及转变成具有生物活性的蛋白质分子之前的所有加工过程，基因的转录调控是基因表达调控中最重要的过程，正确的转录调控使得整个生物体内的能量和资源得以正确的分配［1］。由于调控元件(或称模体)是一些具有保守性功能的片段，在基因的长期进化中具有不变的趋势，因此可以从调控元件的角度来分析基因表达的差异性。不同组织所使用的调控元件既有相同的，也有不同的。可能几个组织共同使用一些调控元件，也可能某个组织单独使用一些调控元件。Yu［2］等人基于转录因子之间的相互作用，识别了人类30个组织特异性基因的调控模块，开发并制作了TiGER数据库［3］。用这个数据库可以方便地查询各个组织的调控模块，但是不同组织可能共同使用了一些调控元件，因此，还缺少各个组织所特有的调控元件信息，即对组织基因调控的差异性缺乏研究。本文利用传统识别调控元件的方法［4］，并应用泊松分布和主成分分析提取出各组织的调控元件，再基于这些调控元件，用统计方法统计出现次数具有显著性特征的特有模体，基于这些特有模体，寻找控制各组织基因调控的元件。

1 方法

1.1 提取过表达模体的方法

提取过表达模体的方法分为两步:第一步是采用泊松分布计算出各个模体出现的概率;第二步是把模体作为变量，每个组织中的每条基因启动子序列作为实验条件，通过主成分分析，确定一组“主要模体”作为过表达模体。

1.1.1 泊松分布

van Helden提出泊松分布来比较两条序列的相似性［5］。他认为模体在某条基因序列上出现是服从泊松分布的，可用泊松分布来度量每个模体出现的概率。由于模体是一些具有保守性的功能片段，在基因的长期进化中具有不变的趋势，因此通过泊松分布度量的每个模体出现的概率大小可以作为模体功能大小的估计。

为不同模体的出现打分，需要定义被考虑的模体在基因组中出现的概率。把每个组织的全部基因作为背景序列来计算背景频率。频率表示给定模体i在组织j中出现的频率，表示模体i在组织j中出现次数的期望，L表示启动子序列的长度。则有

模体i在第k条启动子序列上的出现次数服从泊松分布，泊松分布的分布函数F( x，)表示期望值为时，观察到出现次数≤x的概率。定义第k条基因启动子序列上模体i至少出现次的概率为

采用MATLAB软件编程即可得到每个模体在每条序列中出现的概率。

1.1.2 主成分分析PCA

主成分分析的主要思想是以某些线性组合(主成分)来表示原始数据，再从这些线性组合中尽快提取原始数据的信息［6］。给定n维空间中的m个点，寻求一个 n×n维的矩阵 W，使得 Y=［y1，y2，…，ym］=WTX，同时满足新坐标系下各维之间数据的相关性最小［7］。

主成分分析的一般步骤为［8］:

假设数据为 X=［x1，x2，…，xi，…，xm］，维数为n，在下列所有运算中均有 i∈[1，n]，j∈[1，m]。

(3)计算协方差矩阵 Sn×m，当 Pocc*Pd时，S［a，

(4)对协方差矩阵进行特征分析，并将特征值由大到小排列:λ1＞λ2＞…＞λd，对应的特征向量也作相应排列。

(5)取前 d 个特征值∧d=diag［λ1，λ2，…λd］和特征向量 Wd=［w1，w2，…，wd］，主成分分析可以由X在Wd上投影得到，即:Y=，原始数据重建为:X=WY+。

(6)进一步分析每个主成分对信息的贡献，确定d。令 λi表示第 i(i=1，2，…，d)个特征值，定义第i个主成分的贡献率为则有前 r个主成分的累计贡献率为一般要求累计贡献率达到70%以上。

n个变量的m个观察值，形成一个n×m的数据矩阵，n通常比较大。由于每个新变量是原有变量的线性组合，体现原有变量的综合效果，具有一定的实际意义。采用PCA的方法提取过表达模体，把模体i作为变量，每个组织中的每条基因启动子序列作为实验条件，通过主成分分析，确定一组“主要模体”作为过表达模体。

1.2 不同组织基因调控模体使用的差异性分析

1.2.1 Wilcoxon 秩和检验

Wilcoxon秩和检验［9］并不要求数据满足某种分布假设且这种方法也非常适合小样本数据集。根据基因表达数据的大小排序，然后得到数据的秩，再利用数据的秩而不是数据本身计算基因的秩和统计量。Wilcoxon秩和检验是一种建立在二项分布理论基础上的总体分布位置差异非参数检验方法。

设 X1，…，Xm和 Y1，…，Yn是分别来自具有F( x－ μ1)和F( x－ μ2)连续分布形式的独立总体的两个随机样本，假设

将样本X1，…，Xm和Y1，…，Ym混合在一起，并按从小到大的顺序排列起来。每一个观测值在混合排列中都有自已的秩。令Ri为Yi在这N=m+n个数中的秩(即Yi是第Ri小的)。令Im和In分别表示两样本的指标集，则

同样，对于X样本也可以得到WX:

称WY或WX为Wilcoxon秩和统计量。令

WXY表示在混合样本中所有X的观测值小于Y的观测值的个数。则有

因此WXY(或WYX)也可作为上述检验问题的检验统计量，它们称为Manm-Whitney统计量。当WXY很小时可拒绝零假设。

当m，n均较大时(大于10)，在H0假设下，可用正态分布近似。此时WXY渐近服从均值为，方差为的正态分布。因此可通过标准正态分布作检验，其检验统计量为:

在许多情况下，数据中有相同的数字，称为结(tie)。结中数字的秩为它们按升幂排列后位置的平均值。对于打结的情况，此时大样本近似用的Z值应修正为:

这里τi为结统计量，而t为结的个数。

对于显著性水平α，当p值＜α时，拒绝H0;否则不能拒绝。

1.2.2 基因表达的差异性分析

本文主要考虑不同组织所使用的调控元件的差异情况。具体做法是:把某个组织的过表达模体与其他与之有显著性差异的所有组织的过表达模体相比，提取出那些在这个组织中出现，而在其他与之有显著性差异的所有组织中不出现的模体作为测试集S。并统计出各个模体出现的次数。用超几何分布(Hypergeometric distribution)计算模体在S中出现的概率:

其中N表示与这个组织有显著性差异的组织个数，T表示某个模体mi在测试集S中出现的次数，把与这个组织有显著性差异的组织分为两类，一类是具有模体mi的组织，另一类是不具有模体mi的组织，前者组织个数记为N－n，后者记为n。pt值越小，说明模体mi是某个组织特有的模体的可能性越大。设定阈值b，当pt＜b时，认为模体mi在测试集S中出现的频率显著高于其他模体的出现频率，称模体mi是这个组织的特有模体。然后分析这些特有模体的特征，并分析它们在启动子序列中出现的位置，以此为基础分析各组织使用的模体参与情况。

2 样本和启动子序列的选取

2.1 样本

本文样本包括人的HK(管家)基因以及30组TSP(组织特异性)基因。人HK基因序列数据从HuGEIndex数据库中获得。该数据库给出了451条HK基因的id，剔除无内含子和线粒体中的基因，共得到412条HK基因启动子序列。30组TSP基因，来自Yu等［2］文献中采用聚类方法所得。根据文献中分析得到的30个组织特异性基因的id号，从UCSC数据库提取人30个组织特异性基因的基因序列(见表1)。

表1 各组织的基因条数Table 1 Number of genes in each tissues

2.2 启动子序列的选取

研究表明，基因上游的转录调控位点一般位于翻译起始位点(ATG)上游1 000 bp区域内，目前大部分工作主要集中在基因上游区域。因此采取基因上游1 000 bp区域作为启动子序列。在下载的7 261条基因序列中剔除启动子序列相同的基因序列，共得到4 672条启动子序列。本文主要考察6核苷酸。

3 潜在的转录调控模体

3.1 提取过表达模体

其中 mi(i=1，2，…，4 096)表示模体，gj(j=1，2，…，t)表示组织中的基因启动子序列。为了消除概率很小甚至为零的数据的不利影响，对4 096个模体的出现概率进行加权平均，即对于模体mi，在组织中出现的概率为:

其中nij(j=1，2，…，t)表示模体i在第j条基因序列上出现的数目，Ni表示模体i在这个组织中出现的数目，且有Ni=ni1+ni2+…+nit。

本文提取出现率 pocc=1－p′i＜0.05的模体进行主成分分析，取累计贡献率在70%以上的模体，作为过表达模体，各个组织中提取的过表达模体数目见表2。

表2 各组织中应用PCA提取的过表达模体数目Table 2 Number of over-represented motifs extracted by principal components analysis

提取出的过表达模体的准确率，可以通过与TRANSFAC数据库［10］进行比对。模体识别准确率的计算为识别正确的模体的个数除以提取的过表达模体的个数。通过搜集整理得到TRANSFAC数据库中人的213个转录因子。与TRANSFAC数据库比对，准确率最低是87.50%，最高可达到95.98%(表3)，说明提取出的过表达模体具有生物学上的意义。

表3 各组织的过表达模体与TRANSFAC数据库比对的情况Table 3 Comparison of over-represented motifs with TRANSFAC database

3.2 过表达模体的特征分析

针对过表达模体中碱基的出现情况，将模体进行适当分类［11］。如果6核苷酸中有4个或4个以上的碱基是A或T，认为该6核苷酸富含A、T，称为AT_rich模体;同样如果6核苷酸中有4个或4个以上的碱基是C或G，认为该6核苷酸富含 C、G，称为CG_rich模体;其它既非AT_rich模体又非CG_rich模体，称为CG/AT_lack模体。

HK基因启动子中碱基A+T含量和C+G含量分别为52%和48%，即在HK基因启动子序列中，A+T含量较高。提取出的过表达模体中富含A、T的模体比率高于富含C、G模体的比率。各组织的过表达模体特征见表4。

表4 各组织中模体的碱基使用情况Table 4 Base usage of motifs in each tissues

把上表做成图形以便于更直观的观察(见图1)。

图1 各组织中各富含模体的比率Fig.1 The rate of rich motifs in each tissues

从上图知，各个组织的调控元件中CG/AT_lack模体比率大致相同，而AT_rich模体、CG_rich模体比率变化大，说明其主要差异在于AT_rich模体、CG_rich模体。虽然CG/AT_lack模体在各个组织中的比重都远远高于AT_rich模体、CG_rich模体，但在各个组织中基本都出现，而在各个组织中AT_rich模体、CG_rich模体出现不尽相同。由此可以推测控制基因表达差异性的元件是AT_rich模体或CG_rich模体。从以上分析可以看出，HK基因的调控模体中AT_rich模体的比率高于CG_rich模体比率，而大部分TSP基因(除骨、大脑、宫颈、卵巢、胃、胸腺、舌头外)的调控模体中CG_rich模体的比率高于AT_rich模体的比率。由此可以推测HK基因和TSP基因的调控元件特征不同。即HK基因的调控元件偏向AT_rich模体，而 TSP基因的调控元件偏向CG_rich模体。

4 各组织调控元件使用的差异性分析

一个组织的表现形式不同于另一个组织，虽然是由众多的因素造成，但是否与每个组织自身特有的调控元件有关呢?为此我们考查了不同组织组织调控元件使用的差异性情况。

4.1 组织间调控元件的显著性差异结果

计算每两个组织调控元件的Z统计量，并根据此统计量检验每两个组织的调控元件是否具有显著性差异。结果显示，除组织次级神经系统、睾丸的调控元件与HK基因的调控元件不能判断具有显著性差异外，其他TSP基因的调控元件都与HK基因的调控元件具有显著性差异。

4.2 调控元件的差异性分析

以HK基因为例，HK基因与除次级神经系统、睾丸这两个组织的其他28个组织都具有显著性差异，而通过统计发现 ATTTCG、GGCATA、GGTATA、GTATAG、TATAGT和TATGAC这6个模体在测试集S中都出现了28次，说明这6个模体与在其他28个组织中都不出现的，而只出现在HK基因中，因此，我们可以推测这6个模体很可能是HK基因的特有调控元件。这6个模体所对应的转录因子是HTF。再分析这6个模体，发现它们中有5个是AT_rich模体。若把启动子序列分为五个区域:0～200 bp、200～400 bp、400～600 bp、600～800 bp和 800～1 000 bp。统计发现，这六个模体在启动子序列上的分布情况如下:

表5 HK基因中特有模体在启动子序列上的位置分布Table 5 Location of specific motifs in the promoter sequences of HK genes

由上表知，特有模体的位置分布也具有一定的偏向性，在0～200 bp上分布最密集，其次在400～600 bp上也有分布，而其他序列区域里则分布得少，甚至没有。

对每个组织仿效HK基因，对得到的测试集进行分析并确定这个组织的特有模体。结果显示各个组织的特有模体大部分都不相同，其中不乏共同使用特有模体的。如肾脏与小肠共同使用特有模体CGCATG，这也能从一个侧面得到证实:肾脏和小肠［12］都具有消化功能。膀胱与子宫共同使用特有模体CCGACG，有临床试验表明:子宫切除术后，由于支配膀胱的神经损伤和膀胱解剖位置的改变，常常引起膀胱功能的障碍［13］。再如肌肉、骨髓与肾脏共同使用特有模体GCGCTA，这是由于体内肌肉组织代谢产物肌酐，经血循环到达肾脏，从肾小球滤过后从尿中排泄，因此可以通过测定血清肌酐浓度判断肾小球滤过功能。同时，肾脏分泌促进红细胞生成素，促进骨髓造血，生成红细胞。各组织特有模体的模体特征都偏向于AT_rich模体或CG_rich模体。

5 结语

本文用简单易行的方法提取出各个组织的调控元件及特有模体。结果显示，在管家基因和30个组织特异性基因中，它们的调控元件都具有一定的偏向性，偏向于AT_rich模体或CG_rich模体，而具有CG/AT_lack模体特征的调控元件是极少的。由调控元件的这一偏向规律，在今后的调控元件识别中可以首先忽略CG/AT_lack模体，从而为识别过程缩小了范围。分析发现各组织之间有共同使用的调控元件，例如管家基因和眼睛，其调控元件个数分别为126和132，他们共同使用的调控元件个数为55。然而各个组织也有不少自己单独使用的调控元件即特有模体，这些特有模体专职调控这些组织的基因表达，模体特征都偏向于AT_rich模体或CG_rich模体，与过表达模体的偏向是一致的。由此可推测，特有模体控制着各组织基因的表达。本文的工作对于掌握人基因的转录调控机制和调控模体的作用具有一定的指导意义。

References)

［1］ 孙啸，陆祖宏，谢建明.生物信息学基础［M］.北京:清华大学出版社，2005.SUN Xiao，LU Zuhong，XIE Jianming.Foundation of Bioinformatics［M］.Beijing:Tsinghua University press，2005.

［2］ YU X，LIN J，ZACK D J，et al.Computational analysis of tissue-specific combinatorial gene regulation:predicting interaction between transcription factors in human tissues［J］.Nucleic Acids Research，2006，34(17):4925 －4936.

［3］ LIU Xiong，YU Xueping，DONALD Z，et al.TiGER:A database for tissue-specific gene expression and regulation［J］.BMCBioinformatics，2008，9:271.

［4］ 李婷婷，蒋博，汪小我，等.转录因子结合位点的计算分析方法［J］.生物物理学报，2008，24(5).LI Tingting，JIANG Bo，WANG Xiaowo，et al.The method of calculation and analysis for transcription factor binding sites ［J］.Acta Biophysica Sinica，2008，24(5):334－347.

［5］ VAN HELDEN J.Metrics for comparing regulatory sequences on the basis of pattern counts［J］.Bioinformatics，2004，20(3):399 －406.

［6］ 朱建平.应用多元统计分析［M］.北京:科学出版社，2006.ZHU Jianping.Application ofmultivariate statistical analysis［M］.Beijing:Science Press，2006.

［7］ 李刚，高政.人脸自动识别方法综述［J］.计算机应用研究，2003，(8):4 －9，40.LI Gang，GAO Zheng.A survey of automatic human face recognition ［J］.Application Research of Computers，2003，(8):4 －9，40.

［8］ 侯咏佳，方东博，袁生光，等.主成分分析算法的FPGA实现［J］.机电工程，2008，25(9):37－40.HOU Yongjia，FANG Dongbo，YUAN Shengguang，et al.The principal component analysis algorithm realized by FPGA［J］.Mechanical＆ Electrical Engineering，2008，25(9):37－40.

［9］ 吴喜之.非参数统计［M］.北京:中国统计出版社，1999.WU Xizhi.Non-parametric statistics［M］.Beijing:China Statistical Press，1999.

［10］ WINGENDER E，CHEN X，HEHL R，et al.TRANSFA:an integrated system for gene expression regulation ［J］.Nucleic Acids Res.，2000，28，316 －319.

［11］ LI Huimin，CHEN Dan，ZHANG Jing.Statistical analysis of combinatorial transcriptional regulatory motifs in human intron-containing promoter sequences［J］.Computational Biology and Chemistry，2013，43(2):35－45.

［12］庞智.小肠消化吸收营养基因的区段性表达和调节［J］.国外医学(卫生学分册)，1998，25(4):230 －234.PANG Zhi.The digestion and absorption section expression and regulation of nutrient gene in small intestinal［J］.Foreign Medical Sciences(Health Sciences)，1998，25(4):230－234.

［13］张文淼，施铮铮，吴雪清，等.子宫切除术后膀胱功能障碍患者的尿动力学分析［J］.中华妇产科杂志，2005，40(11):778－779.ZHANG Wenmiao，SHI Zhengzheng，WU Xueqing，et al.Analysis of hysterectomy urinary dynamics in patients with bladder dysfunction ［J］.Chinese Journal of Obstetrics and Gynecology，2005，40(11):778－779.