牙周病的表觀遺傳學分析：生物信息學研究

王田田 劉香瓊 李思敏

牙周炎是一種主要由菌斑生物膜導致的多因素疾病，該病可導致牙齒支持組織的進行性破壞[1]。遺傳和表觀遺傳學因素可通過調節牙周致病菌引起的宿主免疫反應，影響牙周炎的進展[2]。為闡明牙周炎的病因機制，確定牙周炎發病過程中遺傳和表觀遺傳因素的改變顯得非常必要。

非編碼RNA，作為表觀遺傳學的重要調控因素，有兩種主要類型，微小RNA（miroRNAs，miRNAs）和長鏈非編碼RNA（long noncoding RNAs，lncRNAs）。長鏈非編碼RNA，微小RNA，可與信使RNA，三者之間相互作用，形成調節性競爭性內源性RNA網絡[3]。微小RNA含有22個堿基對，能夠負向調控信使RNA的表達[4]。最近使用芯片分析的研究[5，6]表明，牙周炎患者的炎癥牙齦組織和健康牙齦組織之間在微小RNA的表達譜上存在很大差異。長鏈非編碼RNA，作為非編碼RNA的一種，包含200個以上的堿基對，可以通過正向或者負向的方式來調控編碼蛋白質的信使RNA的表達[7]。研究表明長鏈非編碼RNA在慢性牙周炎患者的炎癥牙齦組織中差異性表達[8]。

應用RNA芯片分析方法研究各種RNA在牙周炎中的表達，所得到的結果存在很大的異質性[9]。因此，有必要通過生物信息學的技術，將這些既往研究中的所有的芯片數據集進行整合，從而預測出重要的與牙周炎發病風險有關的生物學標記物。近些年來，已有生物信息學研究[10]確定了與牙周炎發病過程相關的重要的基因、蛋白質，以及生物學通路。然而，基于競爭性內源性RNA理論，對牙周炎發病過程中競爭性內源性RNA網絡的研究尚無報道。

本研究采用生物信息學方法，基于牙周炎的高通量RNA測序及芯片數據，通過整合長鏈非編碼RNA，微小RNA以及信使RNA的表達譜數據，旨在構建參與牙周炎發病機制的競爭性內源性RNA（ceRNA）網絡。

資料和方法

數據采集與差異表達分析：從基因表達數據庫（GEO）中下載與牙周炎樣品匹配的基因表達數據和miRNA表達譜數據。將炎癥性牙齦樣品作為實驗組，健康牙齦樣品作為對照組。使用R語言的limma包中的Linear模型，對牙周炎樣本和牙周健康樣本之間差異表達的mRNA和miRNA進行差異表達分析。P值＜0.05，FDR＜0.05，|log2FC（倍數變化）|≥0.58的mRNA和miRNA分別被認為是顯著差異表達的基因（DEG），以及顯著差異表達的miRNA（DE-miRNA）。

1.RNA互作對數據的提取

從TarBase v7.0，miRTarBase Release 6.0，和miRecords三個數據集，提取已獲得實驗驗證的324507對miRNA-mRNA互作對。從miRcode和starBase database v2.0數據集，提取出1756對lncRNA-miRNA互作對，整合后得到DE-miRNA的mRNA靶標。

2.差異表達基因的共表達分析

為構建共表達網絡和確定潛在的風險基因，利用R語言進行加權基因共表達網絡分析。首先，獲得P值小于0.05的顯著差異表達基因的表達譜，并對基因之間的相互關系的鄰接矩陣進行計算。基于公式（1），對基因互作對間的連接系數進行計算。

χi，χj：基因i和基因j的表達載體；cor：兩個基因載體的皮爾遜相關系數；β：皮爾遜相關系數的指數變換，將皮爾遜相關系數進一步覆蓋到連接系數αij。

基于網絡的拓撲性質，通過加權基因共表達網絡分析，獲得共表達網絡的模塊。該網絡分析了兩種互作的基因以及與其互作的其他基因之間的關系。根據以下公式（2），這兩種互作的基因將連接系數αij覆蓋到加權系數wij上。基于加權系數，獲得了顯著的基因模塊。

wij：代表與基因i相互作用的基因，與基因j相互作用的基因兩者之間的重復。

3.共表達網絡中風險基因模塊的確定

為了篩選出特定的基因模塊，構建了評分分類器。對每個模塊進行差異表達分析，分析了在牙周炎標本和牙周健康標本之間差異表達的每個基因的表達水平，P值被設定為整合表達評分。此外，我們還計算了模塊中每個基因的度（degree）與其-log10（P值）之間的關系。之后，可以篩選出明顯富集的基因模塊。

4.牙周炎ceRNA網絡的構建和分析

從TarBase，miRTarBase和miRecords三個數據庫中獲取到差異表達的9498對mRNA-miRNA相互作用關系對。從miRcode和starBase v2.0兩個數據庫中下載了1726對DEmiRNA-lncRNA的互作對。通 過 將DEG-mRNA-miRNA互 作 對 和DEmiRNA-lncRNA互作對整合到一起，構建了競爭性內源性RNA（ceRNA）網絡。

5.分析結果預測精確度的評估

為了評估本研究中分析方法的預測精確度，我們應用了含有核函數的支持向量機（SVM）。通過計算受試者工作特征（ROC）曲線（AUC）下的面積，預測模型的診斷能力得以評估。通過使用R包中的ROCR.61，整合假陽性率和敏感度的數據，繪制出ROC曲線。此外，為了可視化地顯示分類器的性能，我們進行了無監督的分層聚類分析。通過將完整的聯結和皮爾遜相關性作為距離測量，實現表達譜的分層聚類。

結 果

1.差異表達的mRNA，miRNAs和lncRNAs的確定

通過進行差異表達分析，我們確定了704個差異表達基因和47個差異表達的miRNA。通過使用TarBase，miRTarBase和miRecords這三個數據庫，我們獲得了3282對差異表達的miRNA與其靶基因的互作對，這些互作對是由37個差異表達的miRNA和2650個靶基因構成的。

2.共表達分析

通過差異表達分析，獲得P值＜0.05的基因并對其進行共表達分析。為了建立無標度網絡，我們選擇了相鄰矩陣（β）的權重，此外我們還對不同β值（范圍：1-20）的基因的相鄰矩陣的相關系數和平均值進行了計算。當軟權重為12時，無標度R∧2趨于穩定。

3.顯著模塊的確定

在核查已構建網絡的參數后，通過對加權共表達網絡進行建模，將基因分為不同模塊。從加權關聯網路中，檢索出11個具有不同剪切靜態顏色的基因模塊。隨機抽取出100個基因之后，對所有基因模塊的基因相關性進行計算。每個基因的度k和-log10（p）之間的相關性被加權，加權結果顯示這四個模塊（綠色，紫色，藍綠色，黃色）中的基因與牙周炎的發病機制顯著相關（P＜0.05，圖1b-e）。

4.ceRNA網絡的構建

在所有11個模塊（圖1a）中，綠色和藍綠色模塊的顯著性最高，因此被用于構建ceRNA網絡（圖2a，b）。綠色模塊包含2種DE-miRNA（miR-125a-3p和miR-142-3p）和三種DEGs（HSPA4L，EVI2B和PANK3）（圖2a）。藍綠色模塊含有3種DE-miRNAs（miR-142-3p，miR-125a-3p和miR-200a），和三種DEGs（YOD1，ITGAL，和CTNNBIP1）（圖2b）。在ceRNAs網絡中，3種lncRNA（TUG1，FGD5-AS1，MALAT1）能夠同時直接靶向兩種微小RNA，包括miR-125a-3p和miR-142-3p。沒有一種lncRNA能夠同時直接靶向miR-125a-3p和miR-200a。

圖1 共表達網絡中的模塊

圖2 綠色模塊（a）和藍綠色模塊（b）的競爭性內源性網絡

5.使用兩種額外的獨立隊列進行進一步的基因驗證

為了獨立驗證我們預測出來的這6種基因的診斷能力，基于GPL6244這一平臺，我們在GSE33774這一數據集（含有7個牙周炎樣品和8個牙周健康對照樣本）中，對這些基因的鑒別性能進行了檢測。基于這六種基因的表達譜，我們進行了SVM分析。分析結果（圖3）顯示，HSPA4L基因（被hsa-miR-125a-3p所靶向）的ROC為0.9643，是這6個基因ROC中最高的。這表明，HSPA4L是在牙周炎樣本與牙周健康對照樣本中最為顯著性差異表達的基因。

圖3 ROC曲線（a），預測出的六種風險基因的熱圖（b）

討 論

本研究圍繞三種類型的RNA（信使RNA，微小RNA，長鏈非編碼RNA），對牙周炎的遺傳和表觀遺傳學機制進行剖析。六種信使RNA（HSPA4L，CTNNBIP1，ITGAL，YOD1，PANK3，和EVI2B）被 預測為牙周炎的生物標記物。HSPA4L，熱沖擊蛋白家族的一員，在牙周炎發病機制中發揮雙重作用：一方面可以維持牙齦健康，防止受到口腔微生物的破壞[11]；另一方面，HSPA在某種條件下會引發牙齦炎，加重后向牙周炎轉化，進而增加牙槽骨的損失[12]。CTNNBIP1，又名β-鏈互作蛋白1，可抑制核內β-鏈蛋白的轉錄能力[13]。核內β-鏈蛋白的激活能夠通過抑制非經典的Wnt通路，促進牙周韌帶干細胞的增殖，并導致其成骨分化的減弱[14]。ITGAL，整合素亞基αL，可編碼整合素淋巴細胞功能相關抗 原-1（LFA-1/CD11a）的 整 合 素αL鏈[15]。LFA-1/CD11a是一種重要的免疫調節細胞因子，可通過介導T細胞的召集和保留，調節牙周炎的免疫過程[16]。YOD1，泛蛋白酵素，在牙周炎的樣本中表達下調，可被miR-30b所靶向。MiR-30家族能夠負向地調控骨髓間充質干細胞的骨形態蛋白-2引起的成骨分化[17]。PANK3，泛酰酸激酶3，已被證明可被miR-142-3p所靶向；MiR-142-3p能夠調節脂多糖引起的牙周炎癥[9]。EVI2B，同向性病毒結合位點2B，已被證明在多種微生物感染的大鼠骨模型中表達量增加[18]。

三種微小RNA（miR-125a-3p，miR-200a，和miR-142-3p）被預測可作為牙周炎的風險生物標記物。miR-125a-3p來源于pre-miR-125a的3'末端，在牙周炎樣本中表達上調。MiR-125a的上調可通過調節TRAF6/NFATc1/miR-125a調節反饋回路，抑制破骨細胞的生成，進而抑制牙周炎的進展[19]。miRNA-200a通過減少牙槽骨骨髓基質細胞（BMSCs）的成骨分化，從而抑制牙周骨再生[20]。MiR-142-3p是介導破骨細胞死亡的RANKL的獨立誘導物，可通過調節Wnt通路促進成骨細胞的分化[21]。

三種長鏈非編碼RNA（MALAT1，TUG1，和FGD5-AS1）被預測可作為牙周炎的風險生物標記物。MALAT1，肺腺癌轉移相關轉錄子1，能夠與轉錄因子NF-κB相互作用，調節脂多糖引發的炎癥反應[22]。NF-κB信號通路可以調節RANK配體引起的破骨細胞生成，進而引起牙周炎的骨吸收[23]。TUG1，牛磺酸上調基因1，能夠通過促進轉化生長因子TGF-β，從而抑制細胞的增殖[24]。TGF-β信號通路在牙周炎的發病過程中，參與了牙齦成纖維細胞與牙齦上皮細胞之間的相互作用[25]。FGD5已被證明能夠在牙周炎癥進展過程中，通過調節血管內皮生長因子VEGF的促血管生成作用，來促進血管網絡的延伸[26]。

結論：六種mRNA（HSPA4L，PANK3，YOD1，CTNNBIP1，EVI2B，和ITGAL），三 種miRNA（miR-125a-3p，miR-200a，和miR-142-3p），三種lncRNA（MALAT1，TUG1，和FGD5-AS1）參與了牙周炎的競爭性內源性RNA網絡。這些結果有助于解釋牙周炎的發病機制，為未來研究提供了潛在的治療靶標。