楊 敏，張 靜

(云南大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院統(tǒng)計系，昆明650091)

基因表達是指基因在生物體內(nèi)的轉(zhuǎn)錄、剪接、翻譯以及轉(zhuǎn)變成具有生物活性的蛋白質(zhì)分子之前的所有加工過程，基因的轉(zhuǎn)錄調(diào)控是基因表達調(diào)控中最重要的過程，正確的轉(zhuǎn)錄調(diào)控使得整個生物體內(nèi)的能量和資源得以正確的分配［1］。由于調(diào)控元件(或稱模體)是一些具有保守性功能的片段，在基因的長期進化中具有不變的趨勢，因此可以從調(diào)控元件的角度來分析基因表達的差異性。不同組織所使用的調(diào)控元件既有相同的，也有不同的?？赡軒讉€組織共同使用一些調(diào)控元件，也可能某個組織單獨使用一些調(diào)控元件。Yu［2］等人基于轉(zhuǎn)錄因子之間的相互作用，識別了人類30個組織特異性基因的調(diào)控模塊，開發(fā)并制作了TiGER數(shù)據(jù)庫［3］。用這個數(shù)據(jù)庫可以方便地查詢各個組織的調(diào)控模塊，但是不同組織可能共同使用了一些調(diào)控元件，因此，還缺少各個組織所特有的調(diào)控元件信息，即對組織基因調(diào)控的差異性缺乏研究。本文利用傳統(tǒng)識別調(diào)控元件的方法［4］，并應(yīng)用泊松分布和主成分分析提取出各組織的調(diào)控元件，再基于這些調(diào)控元件，用統(tǒng)計方法統(tǒng)計出現(xiàn)次數(shù)具有顯著性特征的特有模體，基于這些特有模體，尋找控制各組織基因調(diào)控的元件。

1 方法

1.1 提取過表達模體的方法

提取過表達模體的方法分為兩步:第一步是采用泊松分布計算出各個模體出現(xiàn)的概率;第二步是把模體作為變量，每個組織中的每條基因啟動子序列作為實驗條件，通過主成分分析，確定一組“主要模體”作為過表達模體。

1.1.1 泊松分布

van Helden提出泊松分布來比較兩條序列的相似性［5］。他認為模體在某條基因序列上出現(xiàn)是服從泊松分布的，可用泊松分布來度量每個模體出現(xiàn)的概率。由于模體是一些具有保守性的功能片段，在基因的長期進化中具有不變的趨勢，因此通過泊松分布度量的每個模體出現(xiàn)的概率大小可以作為模體功能大小的估計。

為不同模體的出現(xiàn)打分，需要定義被考慮的模體在基因組中出現(xiàn)的概率。把每個組織的全部基因作為背景序列來計算背景頻率。頻率表示給定模體i在組織j中出現(xiàn)的頻率，表示模體i在組織j中出現(xiàn)次數(shù)的期望，L表示啟動子序列的長度。則有

模體i在第k條啟動子序列上的出現(xiàn)次數(shù)服從泊松分布，泊松分布的分布函數(shù)F( x，)表示期望值為時，觀察到出現(xiàn)次數(shù)≤x的概率。定義第k條基因啟動子序列上模體i至少出現(xiàn)次的概率為

采用MATLAB軟件編程即可得到每個模體在每條序列中出現(xiàn)的概率。

1.1.2 主成分分析PCA

主成分分析的主要思想是以某些線性組合(主成分)來表示原始數(shù)據(jù)，再從這些線性組合中盡快提取原始數(shù)據(jù)的信息［6］。給定n維空間中的m個點，尋求一個 n×n維的矩陣 W，使得 Y=［y1，y2，…，ym］=WTX，同時滿足新坐標(biāo)系下各維之間數(shù)據(jù)的相關(guān)性最?。?］。

主成分分析的一般步驟為［8］:

假設(shè)數(shù)據(jù)為 X=［x1，x2，…，xi，…，xm］，維數(shù)為n，在下列所有運算中均有 i∈[1，n]，j∈[1，m]。

(3)計算協(xié)方差矩陣 Sn×m，當(dāng) Pocc*Pd時，S［a，

(4)對協(xié)方差矩陣進行特征分析，并將特征值由大到小排列:λ1＞λ2＞…＞λd，對應(yīng)的特征向量也作相應(yīng)排列。

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

(6)進一步分析每個主成分對信息的貢獻，確定d。令 λi表示第 i(i=1，2，…，d)個特征值，定義第i個主成分的貢獻率為則有前 r個主成分的累計貢獻率為一般要求累計貢獻率達到70%以上。

n個變量的m個觀察值，形成一個n×m的數(shù)據(jù)矩陣，n通常比較大。由于每個新變量是原有變量的線性組合，體現(xiàn)原有變量的綜合效果，具有一定的實際意義。采用PCA的方法提取過表達模體，把模體i作為變量，每個組織中的每條基因啟動子序列作為實驗條件，通過主成分分析，確定一組“主要模體”作為過表達模體。

1.2 不同組織基因調(diào)控模體使用的差異性分析

1.2.1 Wilcoxon 秩和檢驗

Wilcoxon秩和檢驗［9］并不要求數(shù)據(jù)滿足某種分布假設(shè)且這種方法也非常適合小樣本數(shù)據(jù)集。根據(jù)基因表達數(shù)據(jù)的大小排序，然后得到數(shù)據(jù)的秩，再利用數(shù)據(jù)的秩而不是數(shù)據(jù)本身計算基因的秩和統(tǒng)計量。Wilcoxon秩和檢驗是一種建立在二項分布理論基礎(chǔ)上的總體分布位置差異非參數(shù)檢驗方法。

設(shè) X1，…，Xm和 Y1，…，Yn是分別來自具有F( x－ μ1)和F( x－ μ2)連續(xù)分布形式的獨立總體的兩個隨機樣本，假設(shè)

將樣本X1，…，Xm和Y1，…，Ym混合在一起，并按從小到大的順序排列起來。每一個觀測值在混合排列中都有自已的秩。令Ri為Yi在這N=m+n個數(shù)中的秩(即Yi是第Ri小的)。令I(lǐng)m和In分別表示兩樣本的指標(biāo)集，則

同樣，對于X樣本也可以得到WX:

稱WY或WX為Wilcoxon秩和統(tǒng)計量。令

WXY表示在混合樣本中所有X的觀測值小于Y的觀測值的個數(shù)。則有

因此WXY(或WYX)也可作為上述檢驗問題的檢驗統(tǒng)計量，它們稱為Manm-Whitney統(tǒng)計量。當(dāng)WXY很小時可拒絕零假設(shè)。

當(dāng)m，n均較大時(大于10)，在H0假設(shè)下，可用正態(tài)分布近似。此時WXY漸近服從均值為，方差為的正態(tài)分布。因此可通過標(biāo)準正態(tài)分布作檢驗，其檢驗統(tǒng)計量為:

在許多情況下，數(shù)據(jù)中有相同的數(shù)字，稱為結(jié)(tie)。結(jié)中數(shù)字的秩為它們按升冪排列后位置的平均值。對于打結(jié)的情況，此時大樣本近似用的Z值應(yīng)修正為:

這里τi為結(jié)統(tǒng)計量，而t為結(jié)的個數(shù)。

對于顯著性水平α，當(dāng)p值＜α?xí)r，拒絕H0;否則不能拒絕。

1.2.2 基因表達的差異性分析

本文主要考慮不同組織所使用的調(diào)控元件的差異情況。具體做法是:把某個組織的過表達模體與其他與之有顯著性差異的所有組織的過表達模體相比，提取出那些在這個組織中出現(xiàn)，而在其他與之有顯著性差異的所有組織中不出現(xiàn)的模體作為測試集S。并統(tǒng)計出各個模體出現(xiàn)的次數(shù)。用超幾何分布(Hypergeometric distribution)計算模體在S中出現(xiàn)的概率:

其中N表示與這個組織有顯著性差異的組織個數(shù)，T表示某個模體mi在測試集S中出現(xiàn)的次數(shù)，把與這個組織有顯著性差異的組織分為兩類，一類是具有模體mi的組織，另一類是不具有模體mi的組織，前者組織個數(shù)記為N－n，后者記為n。pt值越小，說明模體mi是某個組織特有的模體的可能性越大。設(shè)定閾值b，當(dāng)pt＜b時，認為模體mi在測試集S中出現(xiàn)的頻率顯著高于其他模體的出現(xiàn)頻率，稱模體mi是這個組織的特有模體。然后分析這些特有模體的特征，并分析它們在啟動子序列中出現(xiàn)的位置，以此為基礎(chǔ)分析各組織使用的模體參與情況。

2 樣本和啟動子序列的選取

2.1 樣本

本文樣本包括人的HK(管家)基因以及30組TSP(組織特異性)基因。人HK基因序列數(shù)據(jù)從HuGEIndex數(shù)據(jù)庫中獲得。該數(shù)據(jù)庫給出了451條HK基因的id，剔除無內(nèi)含子和線粒體中的基因，共得到412條HK基因啟動子序列。30組TSP基因，來自Yu等［2］文獻中采用聚類方法所得。根據(jù)文獻中分析得到的30個組織特異性基因的id號，從UCSC數(shù)據(jù)庫提取人30個組織特異性基因的基因序列(見表1)。

表1 各組織的基因條數(shù)Table 1 Number of genes in each tissues

2.2 啟動子序列的選取

研究表明，基因上游的轉(zhuǎn)錄調(diào)控位點一般位于翻譯起始位點(ATG)上游1 000 bp區(qū)域內(nèi)，目前大部分工作主要集中在基因上游區(qū)域。因此采取基因上游1 000 bp區(qū)域作為啟動子序列。在下載的7 261條基因序列中剔除啟動子序列相同的基因序列，共得到4 672條啟動子序列。本文主要考察6核苷酸。

3 潛在的轉(zhuǎn)錄調(diào)控模體

3.1 提取過表達模體

其中 mi(i=1，2，…，4 096)表示模體，gj(j=1，2，…，t)表示組織中的基因啟動子序列。為了消除概率很小甚至為零的數(shù)據(jù)的不利影響，對4 096個模體的出現(xiàn)概率進行加權(quán)平均，即對于模體mi，在組織中出現(xiàn)的概率為:

其中nij(j=1，2，…，t)表示模體i在第j條基因序列上出現(xiàn)的數(shù)目，Ni表示模體i在這個組織中出現(xiàn)的數(shù)目，且有Ni=ni1+ni2+…+nit。

本文提取出現(xiàn)率 pocc=1－p′i＜0.05的模體進行主成分分析，取累計貢獻率在70%以上的模體，作為過表達模體，各個組織中提取的過表達模體數(shù)目見表2。

表2 各組織中應(yīng)用PCA提取的過表達模體數(shù)目Table 2 Number of over-represented motifs extracted by principal components analysis

提取出的過表達模體的準確率，可以通過與TRANSFAC數(shù)據(jù)庫［10］進行比對。模體識別準確率的計算為識別正確的模體的個數(shù)除以提取的過表達模體的個數(shù)。通過搜集整理得到TRANSFAC數(shù)據(jù)庫中人的213個轉(zhuǎn)錄因子。與TRANSFAC數(shù)據(jù)庫比對，準確率最低是87.50%，最高可達到95.98%(表3)，說明提取出的過表達模體具有生物學(xué)上的意義。

表3 各組織的過表達模體與TRANSFAC數(shù)據(jù)庫比對的情況Table 3 Comparison of over-represented motifs with TRANSFAC database

3.2 過表達模體的特征分析

針對過表達模體中堿基的出現(xiàn)情況，將模體進行適當(dāng)分類［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

從上圖知，各個組織的調(diào)控元件中CG/AT_lack模體比率大致相同，而AT_rich模體、CG_rich模體比率變化大，說明其主要差異在于AT_rich模體、CG_rich模體。雖然CG/AT_lack模體在各個組織中的比重都遠遠高于AT_rich模體、CG_rich模體，但在各個組織中基本都出現(xiàn)，而在各個組織中AT_rich模體、CG_rich模體出現(xiàn)不盡相同。由此可以推測控制基因表達差異性的元件是AT_rich模體或CG_rich模體。從以上分析可以看出，HK基因的調(diào)控模體中AT_rich模體的比率高于CG_rich模體比率，而大部分TSP基因(除骨、大腦、宮頸、卵巢、胃、胸腺、舌頭外)的調(diào)控模體中CG_rich模體的比率高于AT_rich模體的比率。由此可以推測HK基因和TSP基因的調(diào)控元件特征不同。即HK基因的調(diào)控元件偏向AT_rich模體，而 TSP基因的調(diào)控元件偏向CG_rich模體。

4 各組織調(diào)控元件使用的差異性分析

一個組織的表現(xiàn)形式不同于另一個組織，雖然是由眾多的因素造成，但是否與每個組織自身特有的調(diào)控元件有關(guān)呢?為此我們考查了不同組織組織調(diào)控元件使用的差異性情況。

4.1 組織間調(diào)控元件的顯著性差異結(jié)果

計算每兩個組織調(diào)控元件的Z統(tǒng)計量，并根據(jù)此統(tǒng)計量檢驗每兩個組織的調(diào)控元件是否具有顯著性差異。結(jié)果顯示，除組織次級神經(jīng)系統(tǒng)、睪丸的調(diào)控元件與HK基因的調(diào)控元件不能判斷具有顯著性差異外，其他TSP基因的調(diào)控元件都與HK基因的調(diào)控元件具有顯著性差異。

4.2 調(diào)控元件的差異性分析

以HK基因為例，HK基因與除次級神經(jīng)系統(tǒng)、睪丸這兩個組織的其他28個組織都具有顯著性差異，而通過統(tǒng)計發(fā)現(xiàn) ATTTCG、GGCATA、GGTATA、GTATAG、TATAGT和TATGAC這6個模體在測試集S中都出現(xiàn)了28次，說明這6個模體與在其他28個組織中都不出現(xiàn)的，而只出現(xiàn)在HK基因中，因此，我們可以推測這6個模體很可能是HK基因的特有調(diào)控元件。這6個模體所對應(yīng)的轉(zhuǎn)錄因子是HTF。再分析這6個模體，發(fā)現(xiàn)它們中有5個是AT_rich模體。若把啟動子序列分為五個區(qū)域:0～200 bp、200～400 bp、400～600 bp、600～800 bp和 800～1 000 bp。統(tǒng)計發(fā)現(xiàn)，這六個模體在啟動子序列上的分布情況如下:

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

由上表知，特有模體的位置分布也具有一定的偏向性，在0～200 bp上分布最密集，其次在400～600 bp上也有分布，而其他序列區(qū)域里則分布得少，甚至沒有。

對每個組織仿效HK基因，對得到的測試集進行分析并確定這個組織的特有模體。結(jié)果顯示各個組織的特有模體大部分都不相同，其中不乏共同使用特有模體的。如腎臟與小腸共同使用特有模體CGCATG，這也能從一個側(cè)面得到證實:腎臟和小腸［12］都具有消化功能。膀胱與子宮共同使用特有模體CCGACG，有臨床試驗表明:子宮切除術(shù)后，由于支配膀胱的神經(jīng)損傷和膀胱解剖位置的改變，常常引起膀胱功能的障礙［13］。再如肌肉、骨髓與腎臟共同使用特有模體GCGCTA，這是由于體內(nèi)肌肉組織代謝產(chǎn)物肌酐，經(jīng)血循環(huán)到達腎臟，從腎小球濾過后從尿中排泄，因此可以通過測定血清肌酐濃度判斷腎小球濾過功能。同時，腎臟分泌促進紅細胞生成素，促進骨髓造血，生成紅細胞。各組織特有模體的模體特征都偏向于AT_rich模體或CG_rich模體。

5 結(jié)語

本文用簡單易行的方法提取出各個組織的調(diào)控元件及特有模體。結(jié)果顯示，在管家基因和30個組織特異性基因中，它們的調(diào)控元件都具有一定的偏向性，偏向于AT_rich模體或CG_rich模體，而具有CG/AT_lack模體特征的調(diào)控元件是極少的。由調(diào)控元件的這一偏向規(guī)律，在今后的調(diào)控元件識別中可以首先忽略CG/AT_lack模體，從而為識別過程縮小了范圍。分析發(fā)現(xiàn)各組織之間有共同使用的調(diào)控元件，例如管家基因和眼睛，其調(diào)控元件個數(shù)分別為126和132，他們共同使用的調(diào)控元件個數(shù)為55。然而各個組織也有不少自己單獨使用的調(diào)控元件即特有模體，這些特有模體專職調(diào)控這些組織的基因表達，模體特征都偏向于AT_rich模體或CG_rich模體，與過表達模體的偏向是一致的。由此可推測，特有模體控制著各組織基因的表達。本文的工作對于掌握人基因的轉(zhuǎn)錄調(diào)控機制和調(diào)控模體的作用具有一定的指導(dǎo)意義。

References)

［1］ 孫嘯，陸祖宏，謝建明.生物信息學(xué)基礎(chǔ)［M］.北京:清華大學(xué)出版社，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］ 李婷婷，蔣博，汪小我，等.轉(zhuǎn)錄因子結(jié)合位點的計算分析方法［J］.生物物理學(xué)報，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］ 朱建平.應(yīng)用多元統(tǒng)計分析［M］.北京:科學(xué)出版社，2006.ZHU Jianping.Application ofmultivariate statistical analysis［M］.Beijing:Science Press，2006.

［7］ 李剛，高政.人臉自動識別方法綜述［J］.計算機應(yīng)用研究，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實現(xiàn)［J］.機電工程，2008，25(9):37－40.HOU Yongjia，F(xiàn)ANG Dongbo，YUAN Shengguang，et al.The principal component analysis algorithm realized by FPGA［J］.Mechanical＆ Electrical Engineering，2008，25(9):37－40.

［9］ 吳喜之.非參數(shù)統(tǒng)計［M］.北京:中國統(tǒng)計出版社，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］龐智.小腸消化吸收營養(yǎng)基因的區(qū)段性表達和調(diào)節(jié)［J］.國外醫(yī)學(xué)(衛(wèi)生學(xué)分冊)，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］張文淼，施錚錚，吳雪清，等.子宮切除術(shù)后膀胱功能障礙患者的尿動力學(xué)分析［J］.中華婦產(chǎn)科雜志，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.