加權基因共表達網絡分析鑒定阻塞性睡眠呼吸暫停的關鍵基因

曹媛媛 蔡昕添 汪思敏 洪 靜 努爾古麗·買買提 李南方

阻塞性睡眠呼吸暫停(obstructive sleep apnea, OSA)是一種以夜間睡眠打鼾、呼吸暫停和日間嗜睡為特征的睡眠呼吸紊亂疾病，可導致間歇性低氧血癥、高碳酸血癥和睡眠結構紊亂[1]。與OSA相關的主要臨床風險是多器官系統損害，如心腦血管疾病、代謝綜合征和認知功能障礙等[1～3]。OSA在普通人群中的發生率呈逐年上升趨勢，據估計，中度至重度OSA在男性中高達49.7%，在女性中高達23.4%[4]。然而，絕大多數OSA患者(70%～90%)沒有得到及時診治，造成了沉重的健康和社會經濟負擔[5]。因此，研究OSA的病因和發病機制，尋找早期診斷指標和治療靶點至關重要。

據報道，OSA有較強的遺傳影響，患者一級親屬的風險增加了1.5倍以上[6]。衡量OSA嚴重程度的呼吸暫停低通氣指數(AHI)中約35%～40%的變異可以用遺傳因素來解釋[7]。目前，基于高通量測序的全基因組DNA微陣列已經成為研究復雜疾病遺傳學的一種有效且相對經濟的工具[8]。雖然已經發現了一些OSA分子標志物，但由于OSA的異質性及其復雜的病理生理狀況，單個基因并不能準確地代表OSA的特征[7，9]。與側重于單個基因的差異表達分析不同，共表達網絡分析通過以無監督的方式識別協同表達的基因模塊，為理解疾病的發病機制和治療干預機會提供了新的見解[10,11]。它已被成功地應用于嚴重哮喘、慢性阻塞性肺疾病(COPD)及癌癥等多種生物學過程的研究，在識別候選生物學標志物和治療靶點方面被證明是相當有效的[11,12]。

本研究應用加權基因共表達網絡分析(WGCNA)，識別與OSA相關的共表達模塊，對其生物學功能和通路進行注釋；鑒定關鍵模塊中的樞紐基因并深入分析，為尋找與OSA發病相關的潛在靶基因提供理論依據。

材料與方法

1.數據的獲?。簭腘CBI的GEO數據庫(http:∥www.ncbi.nlm.nih.gov/geo/)中下載OSA的全基因組表達數據集GSE135917進行初步分析，該數據集基于GPL6244平臺，來源于10例OSA患者和8例正常對照者的皮下脂肪組織樣本，其中包含了每例樣本對應的年齡、性別、體重指數(BMI)及疾病狀況的臨床信息[13]。另一基于GPL6244平臺的獨立數據集GSE38792，包含10例OSA患者和8例正常對照者的內臟脂肪組織樣本，用于后續驗證[14]。

2.數據的處理與篩選：應用R語言及其程序包對含有原始數據(.CEL文件)進行預處理、歸一化和質量控制，采用RMA法進行背景校正、分位數歸一化和中值拋光。根據注釋平臺，將探針ID轉換為基因符號。利用Limma包，以經典貝葉斯t檢驗進行OSA與正常組之間的差異表達分析，以Benjamini-Hochberg法校正P值為錯誤發現率(FDR)。本研究以|log2Fold Change|≥0.585且FDR<0.05為標準篩選出差異表達基因(DEGs)用于共表達網絡的構建。

3.加權基因共表達網絡的構建：通過R語言中的WGCNA包構建DEGs的共表達網絡[11]。利用hclust函數進行聚類分析，剔除數據集中的離群樣本，利用pickSoftThreshold函數確定合適的軟閾值(β)，得到擬合指數R2>0.8的近似無尺度網絡分布。利用blockwiseModules函數進行一步法網絡構建和模塊檢測，生成最小模塊大小為30，合并切割高度為0.25的共表達基因模塊與拓撲重疊矩陣(TOM)。對每個模塊進行主成分分析，以第一主成分計算基因模塊的特征值(module eigengenes, MEs)。引入上述臨床表型，計算MEs與各臨床表型之間的相關系數，識別與OSA顯著相關的基因模塊。模塊內分析基因表達與表型的相關性(gene significance, GS)、與模塊的相關性(module membership, MM)，以篩選關鍵樞紐基因。

4.關鍵模塊的功能富集分析：選取與OSA顯著相關的基因模塊，利用Cytoscape 3.7.2軟件中ClueGo插件進行基因本體論(GO)注釋、京都基因與基因組百科全書(KEGG)通路分析(kappa score=0.4,P≤0.01)。利用STRING 11.0數據庫進行模塊內基因的蛋白-蛋白相互作用(PPI)分析(combined score≥0.7)，利用Cytoscape軟件中cytoHubba插件的Degree法可視化排名前30位的基因(Top30)。

5.關鍵基因的識別與驗證：從關鍵模塊Top30基因構建的PPI網絡中，依據節點度值選取排名前3位的基因為關鍵基因。整合GSE135917和GSE38792數據集，利用ggpubr、ggplot2包對關鍵基因在OSA和正常組織的表達再次驗證，利用pROC包繪制ROC曲線評估關鍵基因的診斷價值和預測OSA的最佳截斷值。

結 果

1.芯片數據的基本信息：經過數據預處理，從18個組織樣本中共獲得23281個基因表達值。以|log2Fold Change|≥0.585且FDR<0.05為標準，OSA與正常組相比，共篩選到3425個DEGs，其中上調基因497個，下調基因2928個?；鹕綀D顯示了DEGs的分布，熱圖顯示了DEGs與樣本的雙向分層聚類結果。為降低噪音及計算機的運行負荷，將3425個DEGs用于下一步共表達網絡的構建(圖1)。

圖1 GSE135917差異基因的篩選 A. 火山圖(紅點表示上調的差異基因；藍點表示下調的差異基因；灰點表示無差異的基因)B. 熱圖(紅色表示差異基因的高表 達；藍色表示差異基因的低表達)

2.加權基因共表達網絡的構建：通過WGCNA算法，依據無尺度網絡分布擬合，選取 β=18作為本數據集的軟閾值，并計算基因間的鄰接矩陣與拓撲重疊TOM，基于TOM構建基因間的分層聚類樹，動態剪切樹法合并MEs相似度較高的模塊，最終把基因聚類成3個模塊,即Turquoise模塊(2345個基因)、Blue模塊(220個基因)、Grey模塊(3個基因)(圖2A)，其中將不能聚類到任何模塊的基因歸于Grey模塊，在后續分析中將其移除。進一步的共表達模塊與臨床表型的相關性熱圖分析(圖2B)顯示，Turquoise模塊與OSA相關性最強(r=-0.98，P=0.000)，以其作為關鍵模塊進行GS與MM分析，Turquoise模塊內基因與臨床表型相關性良好，呈明顯線性相關(圖2C)，分布在右上角的基因既與其他基因關聯性高又與OSA的發病有密切聯系，有助于疾病關鍵基因的識別。因此，筆者用STRING數據庫對Turquoise模塊構建蛋白相互作用網絡，并用Cytoscape軟件可視化節點度最高的前30個基因，結果發現，該模塊存在多個樞紐基因，如SLC2A2、PRL、SST等(圖2D)。

3.關鍵模塊的功能富集分析：利用Cytoscape中的ClueGo插件對Turquoise模塊內基因進行GO和KEGG富集分析，結果以P≤0.01為入選標準。該模塊的基因功能主要注釋于GO:0004984嗅覺受體活性，GO:0043227膜結合細胞器，GO:0005654核質，GO:0045184蛋白質定位。另外，KEGG通路分析顯示該模塊基因主要富集于嗅覺轉導通路(hsa04740)和神經活性配體-受體相互作用通路(hsa04080)(圖3)。

4.關鍵基因的識別與驗證：從Top30基因構建的PPI網絡中，依據節點度值選取排名前3位的基因SLC2A2、PRL、SST作為后續分析與驗證的關鍵基因。筆者整合GSE135917和GSE38792數據集的樣本(包括20例OSA組織和16例正常組織)，通過分析各關鍵基因的表達水平發現，與正常組織比較，SLC2A2、PRL、SST在OSA組織中的表達明顯降低(P均<0.01)，與芯片分析結果一致(圖4A)。隨后，進行ROC分析以確定3個關鍵基因的診斷價值及在基因表達水平預測OSA的最佳截斷值(圖4B，表1)，結果提示這3個關鍵基因很可能與OSA的發生、發展有重要聯系。

圖2 加權基因共表達網絡的構建 A.分層聚類樹與共表達模塊;B.模塊與臨床表型的相關性熱圖;C.Turquoise模塊內基因與臨床表型數據關聯性; D.Turquoise模塊內節點度最高的30個基因構建的蛋白相互作用網絡圖

圖3 Turquoise模塊的功能富集分析 A.GO富集分析;B.KEGG信號通路分析

圖4 關鍵基因的驗證 A.GSE135917與GSE38792中關鍵基因的表達水平;B.關鍵基因的ROC曲線

表1 關鍵基因的ROC曲線分析

討 論

OSA發病機制復雜，涉及多個器官系統的病理生理變化[2]。既往研究多數從單個基因出發，僅能對生物學過程做出局部解釋，WGCNA通過構建基因間的鄰接矩陣、拓撲重疊TOM，識別高度協同變化的基因模塊，并能結合臨床信息，分析模塊與臨床表型的相關性，充分利用基因組大數據信息，對其進行整體全面探索[11]。本研究利用WGCNA法篩選與OSA顯著相關的基因模塊，并對模塊內基因進行功能富集分析，這些基因主要與嗅覺轉導、神經活性配體-受體相互作用等密切相關。隨后對關鍵模塊構建蛋白相互作用網絡，將樞紐基因可視化后識別出3個關鍵基因，通過在另一張芯片上再次驗證，發現SLC2A2、PRL、SST在OSA組織中表達均下調，進一步的ROC曲線分析顯示這3個關鍵基因可能是OSA潛在的生物學標志物。目前3個基因在OSA中尚未有相關報道，但仍有證據提示它們與OSA存在潛在關聯。

SLC2A2是溶質載體家族2成員，編碼葡萄糖轉運子樣蛋白GLUT2，分布于肝、腎、腸、胰島β細胞和中樞神經系統，促進葡萄糖在質膜上的被動轉運，在控制機體葡萄糖穩態中起重要作用[15]。在芬蘭糖尿病預防研究中，SLC2A2的單核苷酸多態性(SNPs)與糖耐量受損向2型糖尿病的轉化相關，而且這種關聯與體重變化無關[16]。Borglykke等[17]在評估46個2型糖尿病相關基因變異的單獨效應和累積效應時發現，只有SCL2A2的小等位基因與心血管事件發生風險增加顯著相關，而且這種關聯與基線糖尿病狀態無關。本研究結果表明，SLC2A2在OSA患者中的表達水平降低，且具有良好的識別能力(AUC=0.9594)。目前OSA對心血管疾病、代謝綜合征和神經精神障礙的不利影響日益受到關注，一些基因變異參與OSA發病機制，且可能與OSA相關疾病的發生也存在因果關系[7]。因此，研究SLC2A2變異是否同時與OSA和共病疾病(如糖尿病、心血管疾病)相關，即SLC2A2變異可能通過不同的機制同時影響這兩種表型，可能具有重要意義。

催乳素(PRL)是一種由腦垂體分泌的蛋白質激素，參與泌乳、生殖、血管生成、免疫反應和滲透調節等多種生物學過程[18]。PRL還可影響睡眠結構，PRL基因缺陷小鼠的快速眼動睡眠比野生型小鼠減少[19]。此外，在高催乳素血癥患者中觀察到代謝改變、體重易增加，這些患者在其催乳素正?；篌w重可減輕[20]。PRL可通過調節LPL活性和脂肪生成，減少脂聯素在人體脂肪組織中的釋放來調節能量代謝[21]。在一項針對早發和病態肥胖的全基因組關聯研究(GWAS)中發現，PRL基因附近的變異與常見肥胖和BMI變異相關[22]。Nilsson等[21]在芬蘭西部的一項大規模人群研究中成功復制了這種關聯，發現PRL基因附近的變異與男性肥胖的增加有關。本研究中，與正常組織比較，OSA患者PRL基因的表達顯著降低。眾所周知，肥胖使OSA的發病風險增加10～14倍，有望解釋高達40%的AHI變異，鑒定決定“中間表型”的基因可能有助于識別OSA的易感基因[7]。因此，有必要進一步研究PRL變異是否通過睡眠結構改變或中間表型(如肥胖)導致OSA。

生長抑素(SST)是一種環肽激素，影響生長激素的釋放和胃腸功能。SST可通過調節胃腸運動、腸道養分吸收和能量平衡來影響機體生長和體重[23]。此外，作為中樞神經系統中廣泛存在的神經遞質或調質，SST在突觸可塑性和神經元活性的微調中發揮作用[24]。研究表明，SST+/nNOS+神經元功能障礙可能導致與慢波活動產生障礙和認知受損相關的病理生理變化，包括阿爾茨海默病(AD)、癲癇、精神分裂癥和創傷性腦損傷等[25]。衰老早期大腦中SST表達下調導致腦啡肽酶活性逐漸下降，導致AD患者Aβ-淀粉樣蛋白的沉積，因此，SST基因變異可能改變生長抑素的表達或功能，從而參與AD的發病過程[26]。本研究中SST在OSA樣本中的表達降低，如上所述，OSA與心腦血管病、神經精神障礙等密切相關，一些基因變異與OSA和相關疾病的發生可能均存在因果關系，故SST變異是否會通過不同的機制效應同時導致OSA和認知功能障礙，有必要進一步驗證。

綜上所述，本研究通過WGCNA方法在OSA患者皮下脂肪組織中發現了2個共表達模塊和3個關鍵基因，有助于為OSA和相關疾病的研究及其診斷治療、靶點選擇提供新的線索。但關鍵基因在OSA發生中的具體作用，仍需要開展進一步的體內、體外實驗予以深入探討。