基于數(shù)據(jù)庫(kù)分析牙髓干細(xì)胞轉(zhuǎn)錄組差異性和功能傾向性研究*

譚思逸 許翀 周珊羽 劉軍

先前的研究已經(jīng)揭示了在成年牙齒中仍能再生的干細(xì)胞群，即出生后人類(lèi)牙髓中的血管周?chē)浪韪杉?xì)胞（dental pulp stem cells,DPSCs），但從老年人牙齒中分離的DPSCs的多能性明顯降低［1］。由此可見(jiàn)，DPSCs 對(duì)形成牙本質(zhì)、維持牙齒營(yíng)養(yǎng)與美觀，提高蓋髓術(shù)成功率尤為重要［2］。

由上可知，人類(lèi)的DPSCs 參與成牙本質(zhì)細(xì)胞產(chǎn)生的牙本質(zhì)結(jié)構(gòu)的再生［3-8］。對(duì)DPSCs再生能力的基本分子機(jī)制理解并掌握，能讓我們有可能利用這些功能用于再生醫(yī)學(xué)和組織工程等未來(lái)細(xì)胞治療方法。從牙髓中分離出的干細(xì)胞已經(jīng)成功在體外實(shí)驗(yàn)分化為成脂肪、軟骨、成骨牙源性譜系［9-11］。隨著高通量測(cè)序技術(shù)的發(fā)展，利用高通量轉(zhuǎn)錄組學(xué)技術(shù)對(duì)DPSCs 的整體基因表達(dá)水平分析得以實(shí)現(xiàn)。機(jī)器學(xué)習(xí)方法如隨機(jī)森林等算法的應(yīng)用，有助于對(duì)高維數(shù)據(jù)的分析和挖掘［12］。本研究對(duì)基因表達(dá)數(shù)據(jù)庫(kù)GEO的GSE123973數(shù)據(jù)集里的DPSCs和骨髓間充質(zhì)干細(xì)胞（bone marrow mesenchymal stem cells,BMMSCs）的RNA-Seq 數(shù)據(jù)進(jìn)行分析，找到DPSCs 的差異高表達(dá)基因，通過(guò)蛋白互作網(wǎng)絡(luò)在差異性高表達(dá)基因集里找到起關(guān)鍵作用的核心基因，并通過(guò)GO和Reactome對(duì)核心基因進(jìn)行功能注釋?zhuān)愿玫卦忈孌PSCs 的特性及其涉及的主要分子作用通路。

1.材料與方法

1.1 RNA-Seq 數(shù)據(jù)檢索 本研究在研究初始階段對(duì)Gene Expression Omnibus（GEO）數(shù)據(jù)庫(kù)里DPSC相關(guān)的高通量基因測(cè)序數(shù)據(jù)集進(jìn)行檢索。

本研究數(shù)據(jù)集的納入標(biāo)準(zhǔn)是：研究?jī)?nèi)容是關(guān)于DPSC與其他來(lái)源的干細(xì)胞比較；考慮到DPSC的異質(zhì)性，我們排除了DPSC樣本小于等于3個(gè)的小樣本研究；使用的是基于測(cè)序技術(shù)獲得的RNA-Seq數(shù)據(jù)集。

截至2021年1月，我們?cè)贕EO 數(shù)據(jù)庫(kù)檢索得到涉及DPSC 的數(shù)據(jù)集有9個(gè)。基于以上標(biāo)準(zhǔn)，僅數(shù)據(jù)集GSE123973 滿足以上三個(gè)要求。所以本研究只對(duì)GSE123973數(shù)據(jù)集進(jìn)行分析。

1.2 RNA-Seq 數(shù)據(jù) 獲取本研究通過(guò)R 包GEOquery［13］從數(shù)據(jù)庫(kù)中獲得GSE123973 數(shù)據(jù)集里的原始RNA-Seq 數(shù)據(jù)，該數(shù)據(jù)集共10個(gè)DPSCs樣本和3個(gè)BMSCs 樣本。10個(gè)DPSCs 樣本中，包含7個(gè)原代培養(yǎng)的DPSCs 樣本和3個(gè)商品DPSCs樣本。原代培養(yǎng)的DPSCs 樣本從人第三磨牙中分離牙髓，提取干細(xì)胞并原代培養(yǎng)后進(jìn)行測(cè)序。3個(gè)人源DPSCs 商業(yè)細(xì)胞樣本購(gòu)于LONZA。人源BMSCs 商業(yè)細(xì)胞樣本購(gòu)于LONZA。測(cè)序樣本的基因比對(duì)參考為人類(lèi)基因組GRCH37。

1.3 差異基因篩選 對(duì)基因數(shù)據(jù)分布進(jìn)行初步分析后，使用R 包DESeq2［14］進(jìn)行基因表達(dá)數(shù)據(jù)的標(biāo)準(zhǔn)化和差異性分析。使用Wald test得出基因表達(dá)的差異是否滿足原假設(shè)。通過(guò)多重檢驗(yàn)BH 法對(duì)P值進(jìn)行了校正。

1.4 蛋白互作網(wǎng)絡(luò)的構(gòu)建 通過(guò)STRING 數(shù)據(jù)庫(kù)構(gòu)建DPSCs 差異表達(dá)基因的蛋白質(zhì)相互作用（protein-protein interaction networks，PPI）網(wǎng)絡(luò)［15］。并對(duì)DPSCs 差異表達(dá)基因相關(guān)的蛋白之間的相互作用網(wǎng)絡(luò)進(jìn)行研究。

1.5 核心基因的功能注釋 通過(guò)蛋白互作網(wǎng)絡(luò)獲得核心基因（圖3）后，對(duì)核心基因進(jìn)行分子通路相關(guān)性功能注釋。應(yīng)用GO［16］和Reactome［17］DPSCs核心基因集的功能傾向性。

1.6 關(guān)鍵調(diào)控基因分析 對(duì)方法1.3 所獲得的差異基因進(jìn)行調(diào)控基因分析。通過(guò)對(duì)差異基因中的轉(zhuǎn)錄因子作為隨機(jī)森林網(wǎng)絡(luò)的調(diào)控基因，構(gòu)建決策樹(shù)。通過(guò)對(duì)決策樹(shù)中的性能和所含差異基因的P值進(jìn)行打分，分值越高的決策樹(shù)被認(rèn)為是越關(guān)鍵的調(diào)控基因。基于該排名，應(yīng)用Cytoscape 3.7 對(duì)排名第一的關(guān)鍵調(diào)控基因網(wǎng)絡(luò)進(jìn)行可視化。

圖1 研究分析流程圖

2.結(jié)果

2.1 RNA-Seq數(shù)據(jù)預(yù)處理 通過(guò)GEOquery共獲得13個(gè)樣本，DPSCs組10個(gè)樣本和BMMSCs組3個(gè)樣本。本研究分析使用的參照基因組序列為Ensembl GRCh37，每個(gè)樣本包含63677個(gè)基因數(shù)據(jù)。對(duì)中位數(shù)計(jì)數(shù)低于3的基因進(jìn)行過(guò)濾，獲得20394個(gè)基因。DPSCs的差異表達(dá)將在這20394個(gè)基因中進(jìn)行分析得出。

2.2 DPSCs大部分基因表達(dá)模式異于BMMSCs通過(guò)過(guò)濾低表達(dá)基因后，對(duì)基因數(shù)據(jù)進(jìn)行差異表達(dá)分析，經(jīng)過(guò)多重檢驗(yàn)后，獲得8200個(gè)差異表達(dá)基因（FDR＜0.01）。取其中組間差異性最大的100個(gè)基因進(jìn)行熱圖聚類(lèi)和可視化，發(fā)現(xiàn)DPSCs和BMMSCs各自得到明確的聚類(lèi)（圖2）。

圖2 DPSCs與BMMSCs差異表達(dá)基因熱圖

2.3 差異基因?qū)?yīng)的蛋白互作網(wǎng)絡(luò)投射出核心基因 由于差異表達(dá)基因較多，在進(jìn)行基因的功能注釋前，需要找到關(guān)鍵的調(diào)控基因，以減少過(guò)多基因造成的背景噪音影響。把差異基因輸入STRING數(shù)據(jù)庫(kù)進(jìn)行蛋白互作網(wǎng)絡(luò)的構(gòu)建，得到差異基因相關(guān)的蛋白互作網(wǎng)絡(luò)。選取高互作得分的節(jié)點(diǎn)（confidence＞0.7）和應(yīng)用MCL聚類(lèi)（inflation parameter為5）結(jié)合，過(guò)濾無(wú)任何連接的節(jié)點(diǎn)，得到核心基因網(wǎng)絡(luò)（圖3）。

圖3 DPSCs核心基因的蛋白互作網(wǎng)絡(luò)

2.4 功能注釋顯示DPSCs 干細(xì)胞功能差異性基于前期蛋白互作網(wǎng)絡(luò)鑒定到的核心基因集，通過(guò)Fish精確試驗(yàn)從GO數(shù)據(jù)庫(kù)中得到其顯著相關(guān)的注釋信息，從而對(duì)核心基因集進(jìn)行功能分類(lèi)注釋。把核心基因改變的倍數(shù)作為參數(shù)，得到主要GO條目及其相關(guān)基因的關(guān)系圖，體現(xiàn)了其干細(xì)胞功能差異性（圖4）。本研究也應(yīng)用Reactome數(shù)據(jù)庫(kù)進(jìn)行通路研究，包括經(jīng)典的代謝通路、信號(hào)轉(zhuǎn)導(dǎo)、基因轉(zhuǎn)錄調(diào)控、細(xì)胞凋亡與疾病。把核心基因改變的倍數(shù)作為參數(shù)，得到主要Reactome條目及其相關(guān)基因的關(guān)系圖（圖5）。

圖4 DPSC核心基因集參與的GO條目及網(wǎng)絡(luò)

圖5 DPSC核心基因集參與的Reactome條目及網(wǎng)絡(luò)

2.5 DPSCs 干細(xì)胞核心基因調(diào)控網(wǎng)絡(luò) 進(jìn)一步對(duì)以上重要差異基因進(jìn)行調(diào)控網(wǎng)絡(luò)的分析，獲得排名第一的調(diào)控基因DMBX1，由圖6 左所示，與BMMSCs 樣本相比，這個(gè)關(guān)鍵調(diào)控基因在DPSCs 樣本中有明顯上調(diào)的差異表達(dá)。進(jìn)一步應(yīng)用網(wǎng)絡(luò)可視化軟件cytoscape 對(duì)這個(gè)關(guān)鍵調(diào)控基因網(wǎng)絡(luò)的連接關(guān)系進(jìn)行可視化，得到其調(diào)控網(wǎng)絡(luò)（圖6 右）：

圖6 DPSC核心調(diào)控基因DMBX1的差異表達(dá)及其調(diào)控網(wǎng)絡(luò)

3.討論

不少前期研究關(guān)注人類(lèi)脫落乳牙干細(xì)胞（SHEDs）與DPSCs 的功能比較。在血管生成功能方面，Zhu SY等的研究發(fā)現(xiàn)SHEDs和DPSCs顯示出與周細(xì)胞相似分化模式，為DPSCs替代周細(xì)胞的臨床應(yīng)用奠定了理論基礎(chǔ)[18]。Wang HH等的研究發(fā)現(xiàn)，在早期和晚期傳代時(shí)，SHEDs比DPSCs表現(xiàn)出更高的增殖率和成骨分化能力，其結(jié)果可能為選擇安全有效的用于再生醫(yī)學(xué)和治療的SHEDs和DPSCs提供指導(dǎo)[19]。盡管這些研究讓口腔科學(xué)研究受益匪淺，但系統(tǒng)地比較DPSCs和其他MSCs的研究還相對(duì)較少。前期研究表明，與DPSCs相似，BMMSCs的再生能力也隨著年齡的增長(zhǎng)會(huì)下降。老年人的骨髓間充質(zhì)干細(xì)胞的增殖和分化能力與年輕人相比顯著下降[20]。另外，BMMSCs可作為不同環(huán)境背景的干細(xì)胞對(duì)照，因此，BMMSCs是研究DPSCs差異性特別是口腔組織分化功能的良好對(duì)照。

前期研究多聚焦在DPSCs和其他來(lái)源的干細(xì)胞的差異表達(dá)基因的數(shù)量和特征上。本研究主要以基因互作網(wǎng)絡(luò)為切入點(diǎn)，通過(guò)對(duì)數(shù)據(jù)庫(kù)中的DPSCs和BMMSCs轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行比較分析，獲得DPSCs的差異表達(dá)基因。這些差異基因在DPSCs中大部分已有表達(dá)，而在BMMSCs卻處于靜止表達(dá)階段，說(shuō)明相對(duì)于BMMSCs，DPSCs更像是一群處于分化階段的細(xì)胞。應(yīng)用差異表達(dá)基因構(gòu)建蛋白互作網(wǎng)絡(luò)，篩選核心基因。通過(guò)選取高得分的節(jié)點(diǎn)和應(yīng)用MCL聚類(lèi)分析結(jié)合，有助于挖掘核心的調(diào)控基因，這些基因及相關(guān)表達(dá)蛋白可能是DPSCs生長(zhǎng)、分化過(guò)程中受影響或者發(fā)揮激活、抑制作用的關(guān)鍵位點(diǎn)。其中獲得關(guān)鍵節(jié)點(diǎn)如細(xì)胞分裂周期蛋白7(cell division cyck 7，CDC7)，苯甲亞胺醇1（budding uninhibited by benzimidazoles 1,BUB1）等。對(duì)核心基因進(jìn)行功能注釋?zhuān)@得與BMMSCs相比DPSCs的基因功能傾向性。GO功能注釋顯示DPSCs差異性GO分類(lèi)主要集中在染色體區(qū)（GO:0098687，chromosomal region；GO:0000775,chromosome，centromeric region；GO:0000793，condensed chromosome），生長(zhǎng)特異性（GO:0006260，DNA replication；GO:0006261，DNA-dependent，DNA replication），以及分化特異性（GO:0001503，ossification）。Reactome數(shù)據(jù)庫(kù)比對(duì)結(jié)果顯示，DNA復(fù)制相關(guān)的基因（R-HSA-68962，Activation of the pre-replicative complex；RHSA-69306，DNA Replication）仍是DPSCs核心基因集參與的通路。另外，細(xì)胞周期的改變也是Reactome數(shù)據(jù)庫(kù)提示的DPSCs核心基因集參與的一系列通路。通過(guò)對(duì)分析挖掘，我們發(fā)現(xiàn)，DPSCs的高表達(dá)基因主要富集在GO和Reactome功能注釋中的細(xì)胞周期（包括DNA復(fù)制）以及分化特性相關(guān)（骨細(xì)胞分化）。骨分化的影響表現(xiàn)了DPSCs的特性，證明了其具有干細(xì)胞多向分化能力的特征,在再生牙髓治療中具有重要潛力[21,22]。最后，我們應(yīng)用隨機(jī)森林算法，構(gòu)建基因調(diào)控網(wǎng)絡(luò)，獲得DPSCs的核心調(diào)控基因。本研究從差異基因整體表達(dá)情況和相關(guān)通路入手，通過(guò)對(duì)數(shù)據(jù)庫(kù)層層深入的分析，引入無(wú)偏估計(jì)的隨機(jī)森林算法處理高維度數(shù)據(jù)獲得最有代表性的轉(zhuǎn)錄因子網(wǎng)絡(luò)，模型泛化能力強(qiáng)。

在DPSCs和BMMSCs差異表達(dá)基因的蛋白互作網(wǎng)絡(luò)的研究中，我們發(fā)現(xiàn)了網(wǎng)絡(luò)關(guān)鍵節(jié)點(diǎn)如CDC7，BUB1等（圖4，圖5）。CDC7參與有絲分裂細(xì)胞周期G1/S轉(zhuǎn)變和DNA復(fù)制啟動(dòng)。BUB1是紡錘體組裝檢查點(diǎn)最主要的組件蛋白之一，它作為監(jiān)控機(jī)制的“感受器蛋白”，影響細(xì)胞的增殖、分化等過(guò)程。這些結(jié)果提示雖然DPSCs和BMMSCs都是干細(xì)胞，但其二者的細(xì)胞周期有所不同，主要表現(xiàn)在G1期。細(xì)胞周期功能是細(xì)胞進(jìn)行細(xì)胞分裂、分化的基礎(chǔ)，其功能的核心包括細(xì)胞周期蛋白、細(xì)胞周期蛋白依賴(lài)性激酶等[23]。越來(lái)越多的證據(jù)表明[24]，不同的細(xì)胞有不同的細(xì)胞周期運(yùn)行機(jī)制，尤其體現(xiàn)在不同的干細(xì)胞上，不同的干細(xì)胞具備其獨(dú)特細(xì)胞周期模式。并且，干細(xì)胞的自我更新和多向分化潛能等特性和功能都與其獨(dú)特的細(xì)胞周期運(yùn)行方式有關(guān)[25]。

最后，有別于以往研究從基因表達(dá)差異為依據(jù)的基因重要性排序方式，我們應(yīng)用隨機(jī)森林算法計(jì)算核心轉(zhuǎn)錄因子及其調(diào)控網(wǎng)絡(luò)在DPSCs和BMMSCs之間的整體差異性。我們發(fā)現(xiàn)轉(zhuǎn)錄因子DMBX1是DPSCs最關(guān)鍵的轉(zhuǎn)錄調(diào)節(jié)因子。DMBX1是中樞神經(jīng)系統(tǒng)發(fā)育相關(guān)的調(diào)控基因，還能調(diào)節(jié)細(xì)胞周期[26]。其在DPSCs的高表達(dá)考慮與牙髓神經(jīng)的發(fā)育相關(guān)。另外，我們還發(fā)現(xiàn)其目標(biāo)基因不僅包括細(xì)胞周期基因如周期蛋白依賴(lài)性激酶樣2 (CDKL2)，還包括HLA II類(lèi)組織相容性抗原HLA-DOB，說(shuō)明該調(diào)控基因在DPSCs中很有可能具有免疫調(diào)控功能。

綜上，本研究通過(guò)DPSCs與BMMSCs的差異基因調(diào)控網(wǎng)絡(luò)分析，以轉(zhuǎn)錄因子和受調(diào)控基因作為節(jié)點(diǎn)，以調(diào)控關(guān)系作為邊的有向網(wǎng)絡(luò)，對(duì)轉(zhuǎn)錄基因與目標(biāo)基因之間的關(guān)系網(wǎng)絡(luò)進(jìn)行精確的量化，從而更深入地闡述DPSCs與BMMSCs間的差異生物學(xué)整體特征。