結(jié)合WGCNA鑒定與豬肌纖維和肌內(nèi)脂肪相關(guān)的中樞基因

劉娟，王舒，左周，路暢，楊陽，蔡春波，趙燕，郭曉紅，曹果清，李步高，高鵬飛

（山西農(nóng)業(yè)大學(xué) 動(dòng)物科學(xué)學(xué)院，山西 太谷030801）

豬肉是人類蛋白質(zhì)和脂肪的主要來源之一，約占世界肉類消費(fèi)的30%。影響豬肉品質(zhì)的因素有遺傳、營養(yǎng)、應(yīng)激等，評定肉質(zhì)的感官指標(biāo)有肉色、pH值、嫩度、系水力、風(fēng)味等［1］。肌纖維組織學(xué)特性包括肌纖維直徑、肌纖維密度、肌纖維類型、肌節(jié)長度、結(jié)締組織和肌內(nèi)脂肪含量，是評價(jià)肉品質(zhì)的重要經(jīng)濟(jì)性狀。肌纖維直徑（Diameter of muscle fiber，DiaMF）與肌肉嫩度相關(guān)，直徑越小，肉質(zhì)越嫩［2，3］。肌纖維密度（Density of muscle fiber，DenMF）和肌內(nèi)脂肪（Intramuscular fat，IMF）含量影響豬肉風(fēng)味、嫩度和多汁性［4］。研究表明，肌內(nèi)脂肪含量低于2.5%的肉類會影響適口性，較高的肌內(nèi)脂肪含量通常被認(rèn)為對肉質(zhì)特征具有積極的影響［5，6］。馬身豬是一個(gè)中國北方豬品種，先前研究分析了馬身豬和大白豬之間與生長速度和肉質(zhì)相關(guān)的表型差異。馬身豬比大白豬具有較低的生長速度和更好的肉質(zhì)，馬身豬的肌內(nèi)脂肪含量高于大白豬［7］。進(jìn)行馬身豬和大白豬的轉(zhuǎn)錄組分析，闡明2個(gè)豬種不同生長發(fā)育階段的分子機(jī)制，發(fā)現(xiàn)馬身豬和大白豬的主要區(qū)別在于氨基酸代謝，脂肪酸代謝和骨骼肌發(fā)育［8］。

高通量測序技術(shù)能夠準(zhǔn)確獲得生物體或組織中特定時(shí)期的基因表達(dá)特征，近年來，國內(nèi)外研究者利用RNA-Seq技術(shù)從轉(zhuǎn)錄組水平探析了不同品種肌肉發(fā)育［9］、肌內(nèi)脂肪沉積［10］的基因表達(dá)模式，構(gòu)建了相關(guān)基因調(diào)控網(wǎng)絡(luò)，發(fā)現(xiàn)并鑒定了PPARGC1B、TRIM 63、CSRP3、ACTN2、MYL 1、MYH 6等一系列基因是影響肌肉生長和脂肪沉積的關(guān)鍵基因［11］，胰島素、絲裂原活化蛋白激酶（MAPK）、Wnt等信號通路能夠顯著影響肌纖維類型和肌肉發(fā)育［12］。WGCNA是一種應(yīng)用于基因組學(xué)、蛋白質(zhì)組學(xué)和代謝組學(xué)中描述相應(yīng)基因、蛋白或代謝物表達(dá)相關(guān)性的系統(tǒng)生物學(xué)分析方法［13］。Matthew等［14］首次將WGCNA應(yīng)用于代謝組學(xué)數(shù)據(jù)，基于WGCNA分析番茄代謝組學(xué)數(shù)據(jù)，確定了與番茄成熟和遺傳背景相關(guān)性狀的3個(gè)主要代謝模塊。Priscila等［15］利用WGCNA揭示了miRNA在內(nèi)洛爾牛脂肪酸組成中的潛在調(diào)控作用；宋萬勤等［16］使用從歐洲生物信息研究所獲得的RNA-Seq數(shù)據(jù)鑒定與間充質(zhì)干細(xì)胞成骨分化有關(guān)的長非編碼RNA，發(fā)現(xiàn)32種新穎的LncRNA是 包 括miR-689，miR-640，miR-601和miR-544等在內(nèi)的miRNA前體。

隨機(jī)森林（Random forest，RF）［17］是一種機(jī)器學(xué)習(xí)領(lǐng)域的集成方法，它利用隨機(jī)重采樣技術(shù)和節(jié)點(diǎn)隨機(jī)分裂技術(shù)構(gòu)建多棵決策樹，通過投票得到最終分類結(jié)果。程淑萍等［19］基于ncRNA和蛋白質(zhì)的序列信息提取特征，訓(xùn)練包括隨機(jī)森林算法在內(nèi)的3個(gè)機(jī)器學(xué)習(xí)模型，預(yù)測ncRNA與蛋白質(zhì)的相互作用，預(yù)測的準(zhǔn)確率較高，證明機(jī)器學(xué)習(xí)方法能夠預(yù)測ncRNA與蛋白質(zhì)是否存在相互作用。

課題組前期以馬身豬和大白豬為研究對象，對不同生長階段肌肉組織進(jìn)行轉(zhuǎn)錄組測序，構(gòu)建了肌肉組織中19條長非編碼RNA（Long noncoding RNA，LncRNA）539條mRNA組成的分子調(diào)控網(wǎng)絡(luò)［20］。本研究采用WGCNA分析技術(shù)并結(jié)合隨機(jī)森林算法，對2個(gè)品種肌肉轉(zhuǎn)錄組測序數(shù)據(jù)深入分析，以期獲得影響豬肌肉發(fā)育及肌內(nèi)脂肪沉積的關(guān)鍵基因。

1 材料和方法

1.1 實(shí)驗(yàn)動(dòng)物與樣品采集

所有動(dòng)物手術(shù)都嚴(yán)格按照《世界醫(yī)學(xué)協(xié)會道德守則》（http：//ec.europa.eu/environment/chemicals/lab_animals/legislation_en.htm）進(jìn)行，實(shí)驗(yàn)方案經(jīng)山西農(nóng)業(yè)大學(xué)動(dòng)物倫理委員會批準(zhǔn)。本研究選用馬身豬和大白豬為研究對象，飼養(yǎng)于山西省大同市種豬場，按照豬飼養(yǎng)標(biāo)準(zhǔn)（NY/T 65-2004）進(jìn)行飼喂。分別在1、30、60、90、120、150、180 d隨機(jī)選擇3頭體重相近的去勢公豬屠宰，采集位于最后肋骨處的背最長肌樣品，液氮保存?zhèn)溆?。按照背最長肌纖維方向垂直切割，采集1 cm3肌肉組織放入4%多聚甲醛固定。

1.2 肌纖維的檢測直徑、密度和肌內(nèi)脂肪含量的測定

肌肉組織在4%甲醛中固定，經(jīng)乙醇脫水、二甲苯透明后用石蠟包埋，6μm切片后采用蘇木精—伊紅（H&E）染色。每個(gè)采樣點(diǎn)平行樣品不少于9張切片，每張切片取3～5張圖像，采用馬尾松三色染色測定肌纖維的直徑和密度。在10×10倍顯微鏡下，對每個(gè)樣品作20個(gè)視野大方格內(nèi)肌纖維根數(shù)的測定（大方格面積約為0.25 mm2），換算后得到每平方毫米纖維根數(shù)，得出肌纖維密度。在10×40倍顯微鏡下，采用不規(guī)則多邊形對每個(gè)樣本進(jìn)行100根肌纖維橫斷面積的描畫，統(tǒng)計(jì)出所有的橫截面積，取其平均值，得到肌纖維面積。根據(jù)肌纖維面積計(jì)算，假設(shè)肌纖維呈圓柱形，由圓面積公式求得肌纖維直徑［21］。采用索氏脂肪抽提法分析肌內(nèi)脂肪含量［22］。

1.3 測序數(shù)據(jù)來源與分析方法

數(shù)據(jù)來源于本課題組前期豬肌肉轉(zhuǎn)錄組數(shù)據(jù)（NCBI登錄號為SRP068558（https：//trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi？login=pda）），共18個(gè)樣本，包含馬身豬和大白豬2個(gè)品種3個(gè)日齡（1、90、180 d）的測序數(shù)據(jù)。使用R包limma［23］進(jìn)行差異基因分析，DEG的閾值設(shè)置為|log2（fold-change）|＞0.8且P＜0.05。

1.4 WGCNA分析

將18個(gè)樣本隨機(jī)分為訓(xùn)練集和測試集，其中訓(xùn)練集含10個(gè)樣本，測試集含有8個(gè)樣本。訓(xùn)練集數(shù)據(jù)采用R-Studio 4.0.0軟件的WGCNA包［13］進(jìn)行分析。首先計(jì)算所有基因的表達(dá)相關(guān)系數(shù)，根據(jù)基因間相關(guān)性確定軟閾值（power），構(gòu)建基因聚類樹。根據(jù)基因表達(dá)相關(guān)系數(shù)確定表達(dá)模塊，每個(gè)表達(dá)模塊最少基因數(shù)設(shè)定為30個(gè)。將每種基因表達(dá)模塊作為整體對象，計(jì)算模塊間的相關(guān)系數(shù)，當(dāng)模塊間相關(guān)系數(shù)高于0.85，則合并相關(guān)模塊。

1.5 GO和KEGG富集分析

基因功能注釋分析采用DAVID［24］（http：//david.abcc.ncifcrf.gov/）在線數(shù)據(jù)庫進(jìn)行。提取GO分析和KEGG富集分析的結(jié)果，分別基于閾值P＜0.05和P＜0.01進(jìn)行分析。

1.6 模塊核心基因的篩選及隨機(jī)森林法驗(yàn)證

核心基因由模塊內(nèi)基因連通性決定。根據(jù)基因的表達(dá)趨勢計(jì)算皮爾遜相關(guān)系數(shù)，取相關(guān)系數(shù)絕對值大于0.2的基因進(jìn)行共表達(dá)分析。利用cytohubba插件［25］篩選每個(gè)模塊的核心基因，使用Cytoscape（version 3.5.1）［26］對 模 塊 構(gòu) 建 互 作 網(wǎng)絡(luò)。利用R語言Random Forest模塊［17］進(jìn)行隨機(jī)森林分析，計(jì)算模塊核心基因的平均基尼指數(shù)減少值（Mean Decrease in the Gini index，MDG），篩選hub基因。采用Graphpad 8.0軟件可視化作圖。

1.7 qPCR檢測hub基因發(fā)育性表達(dá)變化

采用qPCR技術(shù)檢測hub基因在不同月齡2個(gè)品種背最長肌的發(fā)育性表達(dá)變化。在NCBI數(shù)據(jù)庫中獲得篩選得到的hub基因序列，通過Primer3plus在線軟件進(jìn)行引物設(shè)計(jì)，選取18S作為檢測內(nèi)參基因，引物信息詳見表1。采用Prime ScriptTM RT reagent Kit with gDNA Eraser試劑盒（TaKaRa，大連，中國）合成cDNA，反轉(zhuǎn)錄時(shí)總RNA濃度為500 ng·μL-1。使用SYBR Green RTPCR（Taκara，Dalian，China）進(jìn)行qPCR檢測。每個(gè)樣本進(jìn)行3次技術(shù)重復(fù)，分別以無模板和無逆轉(zhuǎn)錄酶的體系作為陰性對照。反應(yīng)程序如下：預(yù)變性階段，95℃30 s；循環(huán)階段，95℃30 s，60℃30 s，設(shè)定循環(huán)數(shù)為45；溶解階段，95℃30 s，60℃1 min，95℃30 s。

表1 關(guān)鍵基因的引物列表Table 1 The list of quantitative primers of hub genes

1.8 統(tǒng)計(jì)分析

qPCR結(jié)果采用2-△△Ct法計(jì)算hub基因的相對表達(dá)量。運(yùn)用SPSS 22.0軟件進(jìn)行基因表達(dá)差異分析，顯著水平選擇0.05和0.01。采用Graphpad 8.0軟件進(jìn)行基因表達(dá)量作圖。

2 結(jié)果與分析

2.1 WGCNA分析結(jié)果

利用R語言WGCNA包計(jì)算訓(xùn)練集樣本中基因表達(dá)相關(guān)系數(shù)并聚類，表達(dá)趨勢性相似的基因合并為一組，構(gòu)建基因共表達(dá)網(wǎng)絡(luò)（圖1）。訓(xùn)練集樣本聚為兩大類，馬身豬和大白豬1日齡樣本聚為第一大類，90日齡和180日齡樣本聚為第二大類，每一大類下各品種分別聚一簇。樣本聚類樹下方為豬背最長肌肌纖維直徑、密度和肌內(nèi)脂肪含量的檢測熱圖，可以看出2個(gè)品種豬背最長肌肌纖維直徑隨著日齡均顯著增加，肌纖維密度逐漸降低，肌內(nèi)脂肪含量增加。同一日齡2個(gè)品種間相比，馬身豬的肌纖維密度顯著大于大白豬，肌內(nèi)脂肪含量也顯著高于大白豬。

圖1 樣本聚類圖和性狀熱圖Fig.1 Cluster diagram of samplesand heatdiagram of traits

采用平均鏈接層次聚類將表達(dá)模式相似的差異基因進(jìn)行分組聚類，建立第一層模塊聚類樹狀圖（圖2A Dynamic Tree Cut），計(jì)算模塊間的相關(guān)系數(shù)，建立模塊特征值聚類樹狀圖（圖2B），當(dāng)模塊相關(guān)系數(shù)大于0.85，軟閾值大于5時(shí)，對模塊進(jìn)行二次合并，建立第二層模塊聚類樹狀圖（圖2A Merged dynamic）。本研究中7262條差異基因首先獲得19個(gè)表達(dá)模塊，根據(jù)模塊相關(guān)性優(yōu)化后最終獲得16個(gè)模塊基因集，之后的分析按照合并后的模塊基因集進(jìn)行。16個(gè)模塊平均含有450條基因，藍(lán)色（blue）模塊中含有的基因數(shù)最多（2960條），灰色60（grey60）和淺綠色（lightgreen）模塊基因最少，僅含50條。

圖2 差異表達(dá)基因基于相關(guān)系數(shù)聚類樹狀圖Fig.2 Dendrogram of all differentially expressed genes clustered based on correlation coefficient

2.2 表型相關(guān)模塊篩選

通過計(jì)算模塊特征與表型間的相關(guān)系數(shù)，獲得模塊性狀相關(guān)性熱圖（圖3）。肌內(nèi)脂肪與MEblack模塊成極顯著正相關(guān)（R=0.78，P=0.008），與MEsalmon模塊（R=-0.73，P=0.02）呈顯著負(fù)相關(guān)。肌纖維密度與MEblack模塊（R=-0.78，P=0.007）和MElightgreen模塊（R=-0.66，P=0.04）存在顯著負(fù)相關(guān)。只有MEmagenta模塊與肌纖維直徑存在極顯著正相關(guān)（R=0.9，P=0.0004）。

圖3 模塊特征基因與肌肉性狀的相關(guān)性熱圖Fig.3 Heat map of correlation between module characteristic genes and muscle traits

2.3 GO和KEGG富集結(jié)果

對MEblack、MEsalmon、MElightgreen和MEmagenta模塊基因集分別進(jìn)行功能富集（圖4）。GO分析結(jié)果表明，MEblack模塊基因功能富集于生物合成、分子修飾及細(xì)胞內(nèi)各種代謝進(jìn)程，MElightgreen模塊基因功能富集于核酸代謝、蛋白分子復(fù)合物組裝及細(xì)胞內(nèi)代謝進(jìn)程，MEsalmon模塊基因功能富集于細(xì)胞質(zhì)、細(xì)胞器、細(xì)胞骨架、轉(zhuǎn)錄因子結(jié)合及一些細(xì)胞代謝進(jìn)程，MEmagenta模塊基因顯著富集于細(xì)胞運(yùn)動(dòng)、細(xì)胞形態(tài)發(fā)生等進(jìn)程。

圖4 GO和KEGG富集圖Fig.4 GO and KEGG enrichment plots

MEblack和MElightgreen模塊與肌纖維密度存在極顯著的負(fù)相關(guān)，這2個(gè)模塊共同富集于代謝和大分子修飾等進(jìn)程，主要包含總代謝進(jìn)程（GO：0008152），細(xì)胞代謝（GO：0044237）、初級 代謝（GO：0044238）等進(jìn)程。MEblack模塊與肌內(nèi)脂肪含量存在極顯著的正相關(guān)，其中碳水化合物代謝（GO：0005975）、小分子代謝（GO：0044281）、氧化還原活性（GO：0016491）等生物進(jìn)程可能是顯著影響肌內(nèi)脂肪沉積的關(guān)鍵分子進(jìn)程。MEmagenta模塊與肌纖維直徑存在顯著的正相關(guān)，主要集中在 細(xì) 胞 骨 架（GO：0007010）、細(xì) 胞 質(zhì)（GO：0005829）和細(xì)胞器（GO：0006996）等生物進(jìn)程。有趣的是MEsalmon模塊沒有與其它3個(gè)模塊共有的生物進(jìn)程，該模塊與肌內(nèi)脂肪含量存在極顯著負(fù)相關(guān)，主要集中于激酶活性（GO：0016301）、離子結(jié)合（GO：0043167）、細(xì)胞及胞內(nèi)組成結(jié)構(gòu)形態(tài)發(fā)生（GO：0000902）等生物進(jìn)程，說明這些進(jìn)程對肌內(nèi)脂肪沉積存在負(fù)調(diào)控作用。

對于KEGG通路富集分析發(fā)現(xiàn)，MEblack模塊 顯 著 富 集 于cGMP-PKG（ko04022）、MAPK（ko04010）、基 質(zhì) 粘 附（ko04510）和 代 謝 途 徑（ko01100）等信號通路。lightgreen模塊基因顯著富集于配體識別（ko04371）、C型凝集素受體（ko04625）、cAMP（ko04024）以 及 生 物 合 成（ko04722）等信號通路。MEmagenta模塊基因極顯 著 富 集 于cGMP-PKG（ko04022）、Wnt（ko04310）、Rap1（ko04015）、cAMP（ko04024）、配體識別（ko04371）、甲狀腺激素（ko04919）、C型凝集素受體（ko04625）、液體剪切應(yīng)力和動(dòng)脈粥樣硬化（ko05418）等信號通路。salmon模塊基因顯著富集于液體剪切應(yīng)力和動(dòng)脈粥樣硬化（ko05418）、cAMP（ko04024）、Wnt（ko04310）等信號通路。

2.4 共表達(dá)子網(wǎng)絡(luò)的構(gòu)建及hub基因的篩選

根據(jù)基因聚類結(jié)果，MEblack、MElightgreen、MEsalmon和MEmagenta4個(gè)模塊分別包含156條、50條、83條和113條差異基因，利用Cytoscape分別構(gòu)建4個(gè)模塊基因互作網(wǎng)絡(luò)。圖5A為MEblack模塊共表達(dá)網(wǎng)絡(luò)，156個(gè)基因共形成11800條連接，根據(jù)cytohubba程序獲得各基因間的連通度，最 終 確 定CENPI、MAP3K 3、XRCC6、FAM 131B、RIC8B、PKP2和PLEKHA 5等基因?yàn)樵撃K的hub基因。圖5B、圖5C和圖5D分別是MElightgreen、MEsalmon和MEmagenta模塊基因形成的表達(dá)網(wǎng)絡(luò)。與MEblack模塊類似，基于基因連通度結(jié)果，最終確定SNRNP48和HK3為MElightgreen模塊的hub基因，F(xiàn)AM 76B、TNXB、AGTRAP和CSDC2為MEsalmon模塊hub基因，AKT 3、CDK 16、RAB21、CBX7、ZNF792和PARP6為MEmagenta模塊hub基因。

圖5 模塊基因的共表達(dá)網(wǎng)絡(luò)Fig.5 The coexpression network of module genes

2.5 隨機(jī)森林對核心基因的檢驗(yàn)結(jié)果

根據(jù)WGCNA分析結(jié)果，共獲得19個(gè)核心基因。將這些核心基因作為樣本特征屬性，建立隨機(jī)森林分類模型對訓(xùn)練集樣本進(jìn)行深度學(xué)習(xí)，通過測試集對模型進(jìn)行檢驗(yàn)，從而確定顯著影響肌肉肌纖維直徑、肌纖維密度和肌內(nèi)脂肪沉積3個(gè)性狀的關(guān)鍵基因。檢驗(yàn)結(jié)果見表2，訓(xùn)練集中共有10個(gè)樣本，只有1 d樣本出現(xiàn)預(yù)測錯(cuò)誤，90 d和180 d的樣本全部預(yù)測準(zhǔn)確，預(yù)測準(zhǔn)確率為90%；測試集預(yù)測結(jié)果與訓(xùn)練集相似，1 d樣品出現(xiàn)預(yù)測錯(cuò)誤，其余樣本均準(zhǔn)確預(yù)測，預(yù)測準(zhǔn)確率為87.5%。

表2 訓(xùn)練集和測試集的分類結(jié)果Table 2 Classification results of training set and test set

本研究中，隨機(jī)森林的預(yù)測準(zhǔn)確率達(dá)到87.5%?；?qū)颖痉诸愵A(yù)測的貢獻(xiàn)率，可根據(jù)基因的MDG篩選（圖6）。CBX7（MDG=1.58）、PLEKHA 5（MDG=0.83）、RIC8B（MDG=0.53）、FAM 131B（MDG=0.52）和PARP6（MDG=0.29）5個(gè)基因?qū)颖痉诸愗暙I(xiàn)最大，最終確定這5個(gè)基因?yàn)橛绊懠±w維發(fā)育和肌內(nèi)脂肪沉積的hub基因。

圖6 基因MDG直方圖Fig.6 Histogram plot of MDG of Gene

2.6 hub基因在豬肌肉組織中的時(shí)序性表達(dá)變化

為了驗(yàn)證hub基因在豬生長發(fā)育過程中的變化趨勢，以馬身豬和大白豬為對象，選擇1、30、60、90、120、150、180 d共7個(gè)時(shí)間點(diǎn)進(jìn)行采樣，采用qPCR檢測hub基表達(dá)（圖7）。大白豬中，CBX7、PARP6、FAM 131B和RIC8B四個(gè)基因均隨著日齡表達(dá)量逐漸升高。與大白豬相比，CBX7和PARP6在馬身豬中基因表達(dá)趨勢一致，隨著日齡表達(dá)量逐漸升高；FAM 131B和RIC8B基因出現(xiàn)先升高后降低的趨勢，在90 d時(shí)表達(dá)量最高。PLEKHA 5基因表達(dá)趨勢較特殊，在2個(gè)品種中均為先降低后升高的變化，90 d時(shí)表達(dá)量降到最低。

圖7 2個(gè)豬品種中5個(gè)關(guān)鍵基因的時(shí)序性表達(dá)Fig.7 Time-sequence expression line diagram of 5 hub genes in 2 pig breeds

3 討論

肌肉中肌纖維的直徑、密度和肌內(nèi)脂肪含量是影響豬肉風(fēng)味、嫩度、多汁性等肉質(zhì)性狀的關(guān)鍵因素［27］。我國地方品種和大白豬、長白豬等國外引進(jìn)豬品種相比，肌纖維直徑小，肌內(nèi)脂肪含量高，系水力強(qiáng)，鮮嫩多汁，肉味香濃［27，28］。近年來隨著高通量測序技術(shù)的快速發(fā)展，國內(nèi)外學(xué)者在轉(zhuǎn)錄組水平對影響肌肉發(fā)育和脂肪沉積的關(guān)鍵基因和分子調(diào)控機(jī)制進(jìn)行了大量研究，獲得了一系列關(guān)鍵基因和信號通路。

Sollero等［29］對皮奧伊州豬和長大二元豬胎兒肌肉組織轉(zhuǎn)錄本進(jìn)行分析，發(fā)現(xiàn)了能夠顯著影響胎兒期肌肉發(fā)育的1300條差異表達(dá)的mRNA。Hu等［30］對沙子嶺豬和大白豬背最長肌進(jìn)行轉(zhuǎn)錄組和蛋白組分析，獲得了488條顯著差異的基因轉(zhuǎn)錄本和23個(gè)差異蛋白，這些基因?qū)ωi肌肉脂肪酸代謝、糖酵解和肌肉發(fā)育等具有顯著的影響。最近課題組利用基因表達(dá)模式分析揭示了影響大白豬發(fā)育的基因（GHSR，EZR，F(xiàn)OXS1，DRD2，SH 3PXD2B，CSF1和TSHR），并 利 用GO和KEGG尋找了與脂肪酸合成相關(guān)的差異表達(dá)基因（ACACA，ACSF3，OXSM，CBR4和HSD17B8），發(fā)現(xiàn)肌內(nèi)脂肪在豬發(fā)育的90～180 d沉積［31］。本研究將WGCNA網(wǎng)絡(luò)分析和隨機(jī)森林深度學(xué)習(xí)相結(jié)合，從大數(shù)據(jù)多維角度探尋影響肌肉性狀的關(guān)鍵基因，最終篩選到5條hub基因和一些關(guān)鍵調(diào)控通路。李鑫等［32］基于網(wǎng)絡(luò)分析和隨機(jī)森林方法對肝癌分期進(jìn)行研究，挑選出BUB1B、CCNA 2、CCNB1、CCNB2、CDC20、MAD2L1、MCM 4、PC-NA、RFC4和TOP2A等10條對肝細(xì)胞癌有重要影響的關(guān)鍵基因。有研究者利用WGCNA挖掘牛皮癬的關(guān)鍵控制因素和驅(qū)動(dòng)因素，使用候選基因構(gòu)建基于隨機(jī)森林的二元分類器，成功將驗(yàn)證組中的正常樣本與牛皮癬樣本區(qū)分開，準(zhǔn)確度為0.95～1，篩選獲得的hub基因作為潛在的牛皮癬病發(fā)進(jìn)程的的治療靶標(biāo)［33］。

本研究發(fā)現(xiàn)CBX7和PARP6基因歸類于MEmagenta模塊中，主要作用細(xì)胞運(yùn)動(dòng)、細(xì)胞骨架形成、細(xì)胞形態(tài)發(fā)生等功能，對肌纖維直徑有促進(jìn)作用。研究發(fā)現(xiàn)CBX7通過基因重組蛋白（Dickkopf-1，DKK-1）介導(dǎo)的Wnt/β-catenin途徑能夠抑制乳腺腫瘤［34］，膠質(zhì)瘤［35］等腫瘤細(xì)胞的增殖。結(jié)合網(wǎng)絡(luò)分析，CBX7在MEmagenta模塊中，MEmagenta模塊基因和脂肪沉積性狀高度相關(guān)，隨機(jī)森林分析發(fā)現(xiàn)，CBX7對樣本分類貢獻(xiàn)最大。本研究中發(fā)現(xiàn)CBX7基因?qū)?nèi)脂肪沉積沒有產(chǎn)生顯著影響，但有研究發(fā)現(xiàn)，在CBX7基因敲除小鼠的脂肪組織含量顯著高于野生型小鼠［36］，說明了CBX7對脂肪沉積也具有重要作用，具體作用機(jī)制有待深入研究。因此可以推斷，CBX7基因是影響脂肪沉積性狀和肌肉發(fā)育的關(guān)鍵基因。研究表明PARP6是MAPK的信號靶點(diǎn)，影響肌肉發(fā)育和脂肪沉積［37］，對大白豬皮下和肌內(nèi)脂肪組織進(jìn)行測序發(fā)現(xiàn)MAPK信號通路中PARP6顯著減少，說明該基因在脂肪細(xì)胞分化過程中具有重要的調(diào)節(jié)作用［22］。PARP6能夠與胰島素反應(yīng)性氨基肽酶（insulin-responsive amino peptidase，IRAP）和GRB14結(jié)合，調(diào)節(jié)脂肪細(xì)胞中葡萄糖轉(zhuǎn)運(yùn)體4（glucose transporter 4，GLUT 4）的活性［37］，猜測PARP6可能是影響肌肉發(fā)育和肌內(nèi)脂肪沉積的關(guān)鍵基因。

在本研究中，F(xiàn)AM 131B、RIC8B、PLEKHA 5均屬于MEblack模塊，能夠促進(jìn)肌內(nèi)脂肪沉積并同時(shí)抑制肌纖維密度增加。FAM 131B基因參與MAPK信號通路，促進(jìn)細(xì)胞增殖，機(jī)體免疫的形成。Huriye Cin等［38］發(fā)現(xiàn)FAM 131B與BRAF結(jié)合，激活MAPK信號通路，導(dǎo)致腫瘤發(fā)生。RIC8B基因編碼GEF，屬于G蛋白介導(dǎo)的信號傳導(dǎo)途徑，在細(xì)胞分裂過程中發(fā)揮重要作用［39］。研究發(fā)現(xiàn)在成骨細(xì)胞增殖過程中C/EBPβ調(diào)控RIC8B基因的高表達(dá)，而在成骨細(xì)胞分化過程中RIC8B基因表達(dá)量顯著下降［40］。PLEKHA 5在在脂質(zhì)和脂蛋白代謝通路中起到重要作用，對腫瘤細(xì)胞的增殖具有重要作用也有相關(guān)報(bào)道［41］。

4 結(jié)論

本研究通過WGCNA分析獲得了與肌纖維發(fā)育和肌內(nèi)脂肪沉積的基因集，構(gòu)建了重要模塊的共表達(dá)網(wǎng)絡(luò)，利用隨機(jī)森林算法對與性狀相關(guān)的功能基因進(jìn)行深度學(xué)習(xí)與驗(yàn)證，最終獲得5條hub基因，包括CBX7、PARP6、PLEKHA 5、FAM 131B和RIC8B。這些基因?qū)∪獍l(fā)育和脂肪代謝存在顯著作用，但具體的作用機(jī)制還需進(jìn)一步研究。