基于轉錄調控模體的人不同組織基因差異性的統計分析

楊 敏，張 靜

(云南大學數學與統計學院統計系，昆明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.