基于位置權重矩陣的核小體識別及功能分析

歲品品，邢旭東，王　宏，崔　穎

( 哈爾濱醫(yī)科大學生物信息科學與技術學院， 哈爾濱 150081)

?

基于位置權重矩陣的核小體識別及功能分析

歲品品，邢旭東，王宏，崔穎*

( 哈爾濱醫(yī)科大學生物信息科學與技術學院， 哈爾濱 150081)

摘要：為研究高通量的人類CD4+T細胞的核小體定位模式，使用迭代算法對核小體定位模式進行分類，并利用位置權重矩陣方法分別構建穩(wěn)定核小體定位序列、動態(tài)核小體定位序列和連接區(qū)序列模型，通過十倍交叉驗證評估模型性能，并與Segal方法與彎曲度方法進行比較，發(fā)現位置權重矩陣方法在敏感性、精度和準確性方面都具有一定優(yōu)越性。同時采用滑窗法在全基因組選取候選序列進行核小體識別，挖掘核小體定位相關基因，并進行基因生物學進程功能富集分析，發(fā)現穩(wěn)定與動態(tài)核小體、真實與潛在核小體對應的基因所參與調控的生物學過程各有不同，但也有一些生物學過程為不同類別核小體所共有，例如對細胞內大分子的調控功能。

關鍵詞：核小體定位；位置權重矩陣；基因功能富集分析

核小體是真核生物染色質的基本組成單位。真核細胞內大約75%～90%的DNA與組蛋白相互纏繞。核小體在DNA序列上的確切位置對DNA序列參與的生物學功能有重要的影響。核小體定位在基因轉錄調控、DNA復制與修復、可變剪接等基本生命過程中都扮演著重要的角色[1-2]，然而，在全基因組上核小體的精確位置卻不是一成不變的，即在DNA序列上核小體定位呈現動態(tài)性，其定位過程、在轉錄過程中的調控功能非常復雜[3-4]。隨著測序技術的快速發(fā)展，ChIP-chip與ChIP-seq等高通量技術已經繪制出核小體定位圖譜，為研究核小體定位及其功能奠定了基礎。同時，這些高分辨率的核小體定位圖譜給采用生物信息方法預測活體內核小體定位提供了豐富的數據樣本，已經開發(fā)出來多種預測算法識別核小體定位。

本文利用位置權重矩陣構建核小體定位模型，在全基因組上識別核小體定位，挖掘核小體定位基因，通過基因富集分析挖掘到核小體參與調控的生物功能，這有利于加強人們對核小體在全基因組的定位模式的全面認識，能夠增加對核小體生物功能的了解，對核小體定位機制的研究以及核小體與基因調控的關系有一定的指導作用。

1材料和方法

1.1數據來源與處理

人類全基因組數據來源于UCSC[5]，并計算全基因組中四種堿基的背景頻率。人類核小體定位數據來自于Dustin E. Schones和 Kairong Cui所做的工作[6]，人類全基因組基因定位數據從ensemble數據庫中[7]。將每條染色體的每一個基因信息分類統(tǒng)計，找出基因的起止位置。根據Segal等人的工作，下載到Segal模型所用的酵母核小體數據[8]。

1.2方法

為了減少DNA序列本身的堿基偏好性對模型的影響，本文把位置頻率矩陣轉換為位置權重矩陣(Position Weight Matrix)。通過引入全基因組背景頻率bi(i∈{A，G，C，T})來消除DNA序列本身堿基組成的偏好性根據公式(1)構建位置權重矩陣模型元素：

(1)

(2)

(3)

其中，qi,j是堿基i在核小體序列第j(j=(1,2,3...147))個位置出現的頻率。元素值Si,j表示堿基i(A,C,G,T)在核小體序列第j位置上的權重值，根據核小體序列集合，獲得位置權重矩陣模型S，為4×147的矩陣，根據公式(3)計算候選序列與模型的相似性，計算其相似性得分，得分越高說明相似性越強。本文利用上述方法分別對穩(wěn)定核小體、動態(tài)核小體、連接區(qū)序列構建位置權重矩陣模型，并分別計算每條候選序列與三個模型的相似性得分，相似性得分最高者判斷為相應模型對應的集合。

2結果

2.1核小體定位模式

利用迭代匹配算法，即將休眠狀態(tài)下的核小體起止位置與激活狀態(tài)下的核小體起止位置進行匹配，獲得4種核小體定位模式。(1)如果核小體定位在激活狀態(tài)下相對于休眠狀態(tài)未發(fā)生任何位置的改變，則該核小體定位定義為完全穩(wěn)定模式(Completely Stable Mode，CSM);(2)如果核小體定位在激活狀態(tài)下相對于休眠狀態(tài)向左或向右移動小于147 bp，則定義為滑動模式(Shift Mode，SM);(3)如果核小體定位在激活狀態(tài)下相對于休眠狀態(tài)向左或向右移動超過147 bp，則定義為完全動態(tài)核小體定位(Completely Dynamic Mode, CDM);(4)如果核小體定位在激活狀態(tài)下相對于休眠狀態(tài)下無核小體定位，則定義為核小體缺失模式(Delete Mode，DEM)(見圖1)。

圖1　核小體定位穩(wěn)定模式和動態(tài)模式Fig.1　The nucleosome position of stable and dynamic pattern

本文分析核小體定位4種模式的DNA序列，發(fā)現CSM和SM存在很大相似性，其模式在DNA序列上的位置變化相對較小，因此將CSM和SM歸為穩(wěn)定模式(Stable Model, SM)，而CDM和DM模式中核小定位呈現非常大的動態(tài)性，兩者可能是同時相互協調發(fā)揮調控作用，因此將CDM和DEM歸為動態(tài)模式(Dynamic Model, DM)。獲得穩(wěn)定模式核小體約53.21%，動態(tài)模式核小體定位約46.79%，核小體定位的多種模型可能和具體的生物過程有關。

2.2模型建立與模型比較

分別構建穩(wěn)定核小體位置權重矩陣(Stable Nucleosome Position Weight Matrix, SNPM)、動態(tài)核小體位置權重矩陣(Dynamic Nucleosome Position Weight Matrix, DNPM)和連接序列位置權重矩陣(Linker Sequence Position Weight Matrix, LSM)，并使用Wilcoxon-test檢驗三個模型間的差異是否具有顯著性，對三個模型中4種堿基在1到147位置的權重差異性如表1和圖2所示，結果表明兩兩模型間4種堿基的差異具有顯著性，即三個模型之間存在顯著差異，此差異為利用三模型識別核小體提供依據。

那天碰巧星期五，臨下班時歐陽鋒接到錢多多的一個電話。錢多多在電話中邀歐陽鋒去鴻運酒樓吃海鮮。歐陽鋒顯得有些猶豫。電話那頭的錢多沒等歐陽鋒回絕，說，歐陽鋒，你不會連這點面子也不給吧？——你在單位門口等著，我馬上讓司機過去接你！

表1　模型間的差異性檢驗結果

圖2　模型間4種堿基對應位置元素值對比Fig.2　The compare of 4 types of bases according to related position in different models

注：彩圖見電子版(http://swxxx.alljournals.cn/index.aspx)(2016年第1期)。

2.3模型性能評估

本文采用十倍交叉驗證方法對模型的性能進行了評估,其性能評估指標為敏感性、特異性、精度和準確性如下列公式所示。

敏感性：Sensitivity=TP/(TP+FN)

(4)

特異性：Specificity=TN/(TN+FP)

(5)

精度：Precision=TP/(TP+FP)

(6)

準確性Accuracy=(TP+TN)/(TP+FN+FP+TN)

(7)

其中TP為真陽性數目、FP為假陽性數目、TN為真陰性數目和FN為假陰性數目十倍交叉驗證，并與Segal[9]和彎曲度模型比較結果如表2所示。通過文獻查找，Segal模型預測的敏感性為68.04%和(陽性)準確性42.32%。對Segal模型所用到的核小體數據進行處理，共得到60 073條釀酒酵母核小體序列和10 030條釀酒酵母連接區(qū)序列，利用位置權重矩陣方法進行模型評估，結果敏感性約為63.8%，特異度約為61.2%，精度約為90.8%，準確性約為63.5%，綜合四項評估指標，位置權重矩陣模型要優(yōu)于Segal模型，與彎曲度譜方法比較。彎曲度譜方法的敏感性為69.85%和(陽性)準確率為59.51%，本文方法敏感性為71.96%和準確性為75.40%均優(yōu)于彎曲度譜方法[10]。因此可以將位置權重矩陣方法應用到人類核小體識別當中。

在全基因組上采用滑窗法，以單堿基為步長，147 bp為窗口寬度來選取候選序列，并去掉含有“N”的候選序列，24條染色體上的總候選序列條數為28億多條，個別染色體候選序列集在硬件存儲大小達到30 G以上，這對于硬件設備是一個嚴峻的考驗。將候選序列集分別投入到模型中，根據打分公式(3)分別計算候選序列與SNPM、DNPM及LSM三個模型的相似度得分，并將候選序列歸類到相似度得分最高的模型分類中。由于候選序列是采用滑動窗口法以單堿基為步長進行提取的，這種方法使候選序列中存在非常大的數據冗余，這使模型的識別結果也存在一定的冗余，為消除這種冗余對模型識別結果的影響，本方對經模型識別后的結果進行了去冗余。去冗余方法：將每個結果中的核小體候選序列與相鄰的核小體候選序列的重疊(超過73 bp)情況合并為核小體區(qū)域，否則不合并，將此結果若核小體識別區(qū)域完全覆蓋實驗核小體定位為正確識別結果即穩(wěn)定核小體定位(Stability Nucleosome Positioning，SNP)和動態(tài)核小體定位(Dynamic Nucleosome Positioning，DNP)，否則認為識別結果為可能存在的核小體定位即潛在的核小體定位，包括潛在穩(wěn)定核小體定位(Potential Stability Nucleosome Positioning，PSNP)和潛在動態(tài)核小體定位(Potential Dynamic Nucleosome Positioning，PDNP)。

表2　模型性能比較

全基因組穩(wěn)定核小體定位識別結果達到64%以上，全基因組動態(tài)核小體定位識別結果約為60%，模型預測的潛在穩(wěn)定的核小體為35%以上，潛在的動態(tài)核小體定位結果為40%以上，此結果與模型評估的準確性基本一致，反應了模型不但有較好的發(fā)現真實核小體的能力，還可以有效地識別全基因組上潛在的核小體。

2.5挖掘核小體相關基因

為了分析核小體定位的功能，分別挖掘到核小體相關基因如圖3所示。四類核小體相關的基因集合間存在很大交疊，但各集合也有相當一部分單獨相關的基因存在。在真實核小體定位的相關基因集合中，大部分基因與真實核小體的兩種定位(真實穩(wěn)定核小體定位與真實動態(tài)核小體定位)都相關。

同樣，在潛在核小體相關基因集合中，大部分基因與潛在核小體的兩種定位(潛在穩(wěn)定核小體定位于潛在動態(tài)核小體定位)都相關，說明大部分基因可能同時受到真實核小體不同定位或者是潛在核小體不同定位的調控作用。相比較而言，穩(wěn)定核小體和潛在核小體相關基因之間的交集較小(真實穩(wěn)定核小體定位與潛在穩(wěn)定核小體定位之間，真實動態(tài)核小體定位與潛在動態(tài)核小體定位之間)，說明真實核小體和潛在核小體同時調控同一個基因的機率相對較小。

圖3　核小體定位相關基因Fig 3　The genes related to nucleosome positioning

注：SNP:真實穩(wěn)定核小體；DNP：真實動態(tài)核小體；

PSNP：潛在穩(wěn)定核小體；PDNP：潛在動態(tài)核小體。

Notes:SNP:Stability Nucleosome Positioning;DNP:Dynamic Nucleosome Positioning;PSNP:Potential Stability Nucleosome Positioning;PDNP:Potential Dynamic Nucleosome Positioning.

2.6功能富集分析

將四類核小體定位相關基因ID分別投入到DAVID9中進行Gene Ontology的Biological Process富集分析。為了使功能富集分析更加詳盡減少冗余，選擇Gene Ontology中的GOTERM_BP_4 ，顯著性閾值P=0.001。并對顯著性P值最小的前10個結果進行展示分析:

(1)如圖4(a)所示，真實穩(wěn)定核小體相關基因富集到的前10個生物學過程中涉及到細胞進程的調控(Positive regulation of cellular process、negative regulation of cellular process)、細胞內大分子調控(Biopolymer modification、cellular protein metabolic process、protein metabolic process、cellular macromolecule catabolic process)、細胞信號與通訊(Regulation of cell communication、regulation of signal transduction、intracellular signaling cascade)以及系統(tǒng)發(fā)育(Nervous system development)。

(2)如圖4(b)真實動態(tài)核小體相關基因富集到的前10個生物學過程中，涉及到細胞進程的調控(Positive regulation of cellular process、negative regulation of cellular process)、細胞內大分子調控(Biopolymer modification、cellular protein metabolic process、positive regulation of macromolecule metabolic process、protein metabolic process positive regulation of macromolecule biosynthetic process、cellular macromolecule catabolic process)、細胞信號與通訊(Regulation of cell communication、regulation of signal transduction)說明真實核小體中，穩(wěn)定核小體與動態(tài)核小體在調控功能上基本相似，除了在各生物學功能的顯著性存在一定的差異外，真實穩(wěn)定核小體還參與神經系統(tǒng)的發(fā)育。

(3)如圖4(c)所示，潛在穩(wěn)定核小體相關基因富集到的生物學過程滿足顯著性閾值的功能有9個。涉及到分化(Keratinocyte differentiation、epidermal cell differentiation、epithelial cell differentiation)、發(fā)育(Epidermis development、ectoderm development)、細胞內大分子調控(Cellular macromolecule biosynthetic process、regulation of macromolecule biosynthetic process)、轉錄事件(Transcription)、RNA代謝(Regulation of RNA metabolic process)。

(4)如圖4(d)所示，潛在動態(tài)核小體相關基因富集到的生物學過程滿足顯著性閾值的功能也只有9個。分別參與的功能為分化(Epidermal cell differentiation、keratinocyte differentiation、epithelial cell differentiation)、發(fā)育(Epidermis development、ectoderm development、organ development、tissue development)、細胞內大分子調控(Cellular macromolecule biosynthetic process、protein-lipid complex assembly)。說明真實動態(tài)核小體和潛在動態(tài)核小體除在細胞大分子調控功能上類似、其他功能有很大差異。四種核小體都參與的功能為細胞內大分子的調控功能。

圖4　基因富集分析結果Fig.4　The result of functional enrichment analysis

3結果與討論

本文通過建立位置權重矩陣模型來識別核小體定位。研究結果顯示，位置權重矩陣模型具有較高的敏感性和準確性，但假陽性率仍然比較高，原因可能是候選序列中每一條真正核小體前后都有與真正核小體相近的打分，但較高的假陽性也為挖掘試驗中沒有發(fā)現的核小體奠定了基礎。另外，通過對挖掘得到的核小體相關基因進行功能富集分析，發(fā)現穩(wěn)定與動態(tài)核小體、真實與潛在核小體對應的基因所參與調控的生物學過程各有不同，但也有一些生物學過程為不同類別核小體所共有，例如對細胞內大分子的調控功能。利用位置權重矩陣模型對在全基因組內選取的候選序列進行識別，除了發(fā)現實驗中的已經發(fā)現的真實和穩(wěn)定核小體之外，還挖掘到了一些具有核小體可能性的序列。對不同類別核小體相關的基因進行功能富集分析，發(fā)現真實與潛在、穩(wěn)定與動態(tài)核小體區(qū)域相關的基因所參與調控的生物學過程這對核小體定位機制以及核小體與基因調控的關系的研究有一定的指導意義。我們推測一方面細胞通過全基因組范圍內核小體定位模式的一致性來維持細胞的正常功能，另一方面細胞通過內部各類核小體定位模式的差異來發(fā)揮核小體的調控作用。不同的生長階段、生理條件下細胞內基因的表達水平可能存在不同，其受很多因素的調控[12]。核小體通過具體的動態(tài)位置變化來隱蔽或暴露DNA上的蛋白結合位點，這些蛋白結合位點往往與轉錄因子等和基因表達緊密相關的蛋白質相結合來調控基因表達。雖然至今仍不能確定核小體定位的動態(tài)變化是引起基因表達水平變化的決定因素，但是至少兩者之間存在著緊密的聯系，值得進一步探索。

參考文獻

[1]陳偉. 核小體定位對RNA剪接的影響及組蛋白變體的識別[D]. 呼和浩特:內蒙古大學, 2010.

CHEN Wei.The effect of nucleosome positioning on RNA splicing and the regcognition of histone rariants[D].Hohhot:Inner Mongolia University,2010.

[2] 蔡祿, 趙秀娟. 核小體定位研究進展[J]. 生物物理學報, 2009,25(6): 385-395.

CAI Lu, ZHAO Xiujuan. Advances in nucleosome positioning [J].Acta Biophysica Sinica,2009, 25(6): 385-395.

[3]SCHONES D E, CUI K, CUDDAPAH S, et al. Dynamic regulation of nucleosome positioning in the human genome[J]. Cell, 2008, 132(5): 887-898.

[4]JIANG C, PUGH B F.Nucleosome positioning and gene regulation:advances through genomics[J].Nature Reviews Genetics, 2009, 10(3): 161-172.

[5]KENT W J, SUGNET C W, FUREY T S, et al. The human genome browser at UCSC[J]. Genome Research, 2002, 12(6): 996-1006.

[6]ZHANG Y, SHIN H, SONG J S, et al. Identifying positioned nucleosomes with epigenetic marks in human from ChIP-Seq[J]. BMC Genomics, 2008, 9(1): 537.

[7]HUBBARD T J P，AKEN B L，BEAL K, et al. Ensembl 2007[J]. Nucleic Acids Research, 2007, 35(Suppl 1): D610-D617.

[8]SEGAL E, FONDUFE-MITTENDORF Y, CHEN L, et al. A genomic code for nucleosome positioning[J]. Nature, 2006, 442(7104): 772-778.

[9]GLYNN D J R , BRAD T S , DOUGLAS A H, et al. DAVID: Database for annotation,visualization,and integrated discovery[J]. Genome Biology, 2003, 4:R60(9): 54-56.

[10]SEGAL E, FONDUFE-MITTENDORF Y, CHEN L, et al. A genomic code for nucleosome positioning[J]. Nature, 2006, 442(7104): 772-778.

[11]張德金, 劉宏德, 袁志棟, 等. 基于Web技術的核小體在線預測平臺實現[J]. 微計算機信息, 2010, 26(36): 185-187.

ZHANG Dejin,LIU Hongde,YUAN Zhidong,et al.Construction of an on-line platform of predicting nucleosomes based on web techniques[J].Microcomputer Information, 2010, 26(36): 185-187.

[12] TEIF V B，VAINSHTEIN Y，CAUDRON-HERGER M， et al. Genome-wide nucleosome positioning during embryonic stem cell development[J].Nature Structural & Molecular Biology, 2012, 19(11):1185-1192.

Nucleosome positioning identification and functional analysis on the position weight matrix

SUI Pinpin，XING Xudong，WANG Hong，CUI Ying*

(CollegeofBioinformaticsScienceandTechnology,HarbinMedicalUniversity,Harbin150081,China)

Abstract：This study was based on high throughout nucleosome positioning data of CD4+T cell in human genome to investigate the model of nucleosome positioning and category the nucleosomes.We constructed three the models by using position weight matrix, including stable nucleosome model, dynamic nucleosome model and linker sequences model respectively. Ten-fold cross validation was used to evaluate the performance of the three models, and the assessment results were compared with Segal model and curvature profile method.It was found that the position weight matrix method was superior to the other two methods in terms of sensitivity, precision and accuracy. At the same time the sliding window method is adopted to select candidate sequences in the genome to identify the nucleosomes. Furthermore we mined the related genes of nucleosome positioning and completed enrichment analysis of gene functions and found that different nucleosome positioning modes involved in both a certain similarity and difference in regulation function in biological processes.Whereas there are some biological processes are co-regulated by different nucleosome positioning modes,such as regulation of macromolecule.

Keywords：Nucleosome; Position weight matrix; Enrichment analysis of gene function

中圖分類號：Q523

文獻標志碼：A

文章編號：1672-5565(2016)01-001-06

doi:10.3969/j.issn.1672-5565.2016.01.01

作者簡介：歲品品，男，本科生，研究方向：生物信息學；E-mail:1447806377@qq.com.*通信作者：崔穎 ，女，講師,研究方向：生物信息學；E-mail:cuiying204@163.com.

基金項目：黑龍江省衛(wèi)生廳科研課題資助(2013129)。

收稿日期：2015-11-08;修回日期：2015-12-11.