加權(quán)基因共表達(dá)網(wǎng)絡(luò)在肌少癥基因篩選中的應(yīng)用研究

黃玲莉

（上海理工大學(xué)健康科學(xué)與工程學(xué)院，上海 200093）

0 引言

肌少癥是一種隨著年齡增長(zhǎng)而導(dǎo)致進(jìn)行性肌肉量減少、肌肉力量下降或軀體功能減退的老年綜合征?；加屑∩侔Y的個(gè)體會(huì)因骨骼肌衰退、骨—肌單位萎縮，容易發(fā)生跌倒造成骨折，并且該癥狀與認(rèn)知障礙、呼吸系統(tǒng)疾病和心血管疾病相關(guān)，會(huì)導(dǎo)致患者生活質(zhì)量下降。研究表明，從40歲至80歲，人體的骨骼肌總量將下降30%～50%；60歲后，肌肉功能每年將下降3%；80歲后，患肌少癥的概率高達(dá)67.1%。

隨著中國(guó)進(jìn)入老齡化社會(huì)，肌少癥的預(yù)防和治療已受到了人們的廣泛關(guān)注，但目前人們對(duì)肌少癥的認(rèn)識(shí)仍然較少，發(fā)病機(jī)理尚不明晰。因此，迫切需要挖掘肌少癥的相關(guān)基因用于研究其發(fā)病機(jī)理，并為預(yù)防、診斷和治療該疾病提供參考。

1 研究現(xiàn)狀

隨著步入大數(shù)據(jù)時(shí)代，高通量測(cè)序技術(shù)、網(wǎng)絡(luò)方法和計(jì)算機(jī)技術(shù)快速發(fā)展，從計(jì)算角度挖掘疾病的生物標(biāo)志物或治療靶點(diǎn)已成為目前的研究熱點(diǎn)。王莉華等通過(guò)差異基因分析法對(duì)肌少癥進(jìn)行分析，但該方法忽略了基因間的內(nèi)部相關(guān)性，即具有相同或相似表達(dá)的基因可能具有相似功能。Shin等利用TargetScan算法從Dlk1-Dio3 miRNA中挖掘有關(guān)肌少癥的靶點(diǎn)基因，但該算法的處理對(duì)象為單個(gè)miRNA，無(wú)法同時(shí)考慮全基因組信息。為此，Zillikens等提出了全基因組關(guān)聯(lián)分析方法（Genome-wide Association Study，GWAS）以解決TargetScan算法中存在的問(wèn)題，該方法從基因HSD17B11、VCAN、ADAMTSL3、IRS1、FTO中發(fā)現(xiàn)了有關(guān)肌少癥的5個(gè)新位點(diǎn)，但GWAS會(huì)隨著測(cè)序數(shù)據(jù)量的增加，影響計(jì)算速度。同時(shí)，GWAS存在多重檢驗(yàn)問(wèn)題，會(huì)導(dǎo)致篩選的基因不完全是致病基因。

然而，加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（Weighted Gene Coexpression Network Analysis，WGCNA）的提出為解決以往研究存在的問(wèn)題提供了新的思路。該方法是一種基于拓?fù)渚仃囁鶚?gòu)建的共表達(dá)網(wǎng)絡(luò)，使用無(wú)監(jiān)督聚類法確定簇的數(shù)據(jù)挖掘方法。相較于上述方法，WGCNA可從整體基因組進(jìn)行研究，通過(guò)確定基因表達(dá)水平間的關(guān)聯(lián)及生物功能的相似性，將基因集分為一個(gè)個(gè)模塊以挖掘關(guān)鍵基因，這樣既能縮小篩選范圍，減少計(jì)算量，又能較好地控制多重檢驗(yàn)的假陽(yáng)性率，提高識(shí)別準(zhǔn)確率。

目前，WGCNA已被廣泛應(yīng)用于腫瘤、帕金森、糖尿病等疾病的研究當(dāng)中，其穩(wěn)定性和有效性已在多項(xiàng)研究中得以證實(shí)。例如，Su等使用WGCNA篩選卵巢癌的關(guān)鍵基因。Haase等將WGCNA應(yīng)用于研究Rett綜合癥。而本文首次運(yùn)用WGCNA分析肌少癥，以期為肌少癥的致病機(jī)制及治療研究提供新的生物標(biāo)志物或候選靶點(diǎn)分子。

2 基于WGCNA的基因識(shí)別

本文利用WGCNA篩選肌少癥基因，具體實(shí)驗(yàn)過(guò)程包含以下6個(gè)部分：①對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理；②依據(jù)基因間的關(guān)聯(lián)性構(gòu)建共表達(dá)網(wǎng)絡(luò)；③利用共表達(dá)網(wǎng)絡(luò)對(duì)基因集進(jìn)行模塊識(shí)別；④將劃分后的模塊與外部性狀進(jìn)行關(guān)聯(lián)分析，以確定關(guān)鍵模塊；⑤在關(guān)鍵模塊中篩選肌少癥的樞紐基因；⑥對(duì)樞紐基因進(jìn)行功能驗(yàn)證。

2.1 數(shù)據(jù)預(yù)處理

本文使用的肌少癥相關(guān)基因表達(dá)數(shù)據(jù)和臨床信息（GSE111016）來(lái)源于美國(guó)國(guó)立生物信息中心的GEO（Gene Expression Omnibus）數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/gds/），該數(shù)據(jù)庫(kù)是目前國(guó)際上最全面的高通量基因表達(dá)公共數(shù)據(jù)庫(kù)，其中包含了20例肌少癥患者和20例健康對(duì)照者，然而發(fā)現(xiàn)數(shù)據(jù)庫(kù)中的樣本3存在數(shù)據(jù)缺失問(wèn)題，需要將其剔除。

圖1為通過(guò)類平均法計(jì)算樣本間的距離所得到的樣本聚類樹(shù)。由圖1可見(jiàn)，將聚類高度設(shè)為77時(shí)，樣本11、樣本13和樣本18為離群樣本，需要先將其剔除。然后，對(duì)剩余的基因表達(dá)數(shù)據(jù)進(jìn)行l(wèi)og2轉(zhuǎn)換，并進(jìn)行TMM歸一化處理，最終得到36例臨床信息的16 879個(gè)基因表達(dá)數(shù)據(jù)用于后續(xù)研究分析。

Fig.1 Sample hierarchical clustering tree圖1 樣本分層聚類樹(shù)

2.2 共表達(dá)網(wǎng)絡(luò)構(gòu)建

基因共表達(dá)網(wǎng)絡(luò)關(guān)系的構(gòu)建通常由鄰接矩陣定義，該矩陣的元素大小表示兩個(gè)節(jié)點(diǎn)間的連接強(qiáng)度。本研究在構(gòu)建鄰接矩陣前，先按式（1）、式（2）定義一個(gè)中間變量，即相似度矩陣S。

其中，

x

，

x

表示任意兩個(gè)基因。

然后，利用閾值化方法將相似度矩陣S轉(zhuǎn)化為鄰接矩陣A，傳統(tǒng)鄰接矩陣計(jì)算公式為：

其中，τ為所取閾值。

但由于該閾值法為硬閾值，無(wú)法反映共表達(dá)信息的連續(xù)性，易丟失部分信息，并且該方法只能表示基因間是否存在聯(lián)系，無(wú)法表示關(guān)系的強(qiáng)弱。為此，本文通過(guò)引入軟閾值

β

對(duì)相似度進(jìn)行冪運(yùn)算，計(jì)算公式如式（5）：

其中，

β

≥1。

接下來(lái)，使用R軟件的pick Soft Threshold（）函數(shù)對(duì)軟閾值進(jìn)行篩選，并對(duì)篩選出的軟閾值進(jìn)行無(wú)尺度拓?fù)錂z查。該操作決定了所構(gòu)建網(wǎng)絡(luò)是否符合無(wú)尺度網(wǎng)絡(luò)，冪運(yùn)算使得之前關(guān)系緊密的基因不受影響或者影響較小，而相關(guān)性較弱的基因相關(guān)性則會(huì)明顯下降，進(jìn)而表示關(guān)系的強(qiáng)弱。

考慮到兩個(gè)基因間除直接影響外，還會(huì)受到其它基因的間接調(diào)控，本文根據(jù)拓?fù)渲丿B公式將鄰接矩陣A轉(zhuǎn)化為拓?fù)渲丿B矩陣Ω。該矩陣既能在網(wǎng)絡(luò)拓?fù)渌缴象w現(xiàn)基因間的強(qiáng)度大小，還能減少噪音和假陽(yáng)性，使結(jié)果更穩(wěn)健。計(jì)算公式如式（6）-式（8）所示：

其中，

I

表示引入第3個(gè)基因

u

之后，第

i

個(gè)基因和第

j

個(gè)基因新的連通性大小。

2.3 基因模塊識(shí)別

本文使用無(wú)監(jiān)督聚類法確定基因模塊，具體通過(guò)層次聚類算法實(shí)現(xiàn)，如公式（9）所示：

同時(shí)，利用動(dòng)態(tài)剪切樹(shù)法將層次聚類樹(shù)切分為若干模塊，每個(gè)模塊表示共表達(dá)程度較高的一類基因。

2.4 關(guān)鍵模塊確定

首先對(duì)每個(gè)模塊矩陣進(jìn)行主成分分析，得到模塊特征值。然后，提取聚類后的模塊特征值與外部性狀向量，以計(jì)算Pearson的相關(guān)系數(shù)。最后，選擇與肌少癥具有顯著關(guān)系（

p

＜0.05）的關(guān)鍵模塊進(jìn)行后續(xù)樞紐基因的篩選。同時(shí)，提取關(guān)鍵模塊內(nèi)每個(gè)基因的基因顯著性值（Gene Significance，GS）和模塊身份值（Module Membership，MM）以進(jìn)一步驗(yàn)證所選模塊的可靠性。通過(guò)分別計(jì)算GS和MM的相關(guān)系數(shù)，若其結(jié)果越顯著（

p

＜0.05），則表明所選模塊越穩(wěn)健。

2.5 樞紐基因篩選

計(jì)算關(guān)鍵模塊內(nèi)各基因的|MM|、|GS|值，可有效提高篩選樞紐基因的準(zhǔn)確性。因此，本文通過(guò)分析大量數(shù)據(jù)后，選擇|MM|＞0.8、|GS|＞0.2、

P

＜0.05為篩選條件。

2.6 基因功能驗(yàn)證

GO（Gene Ontology）是研究基因功能的重要數(shù)據(jù)庫(kù)，現(xiàn)已被廣泛用于分析、驗(yàn)證高通量基因數(shù)據(jù)?；贕O對(duì)篩選出的基因進(jìn)行富集分析，當(dāng)

P

＜0.05時(shí)表示基因組與對(duì)應(yīng)的功能具有統(tǒng)計(jì)學(xué)意義，證明所選基因是合理、有效的。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 WGCNA結(jié)果與分析

通過(guò)pick Soft Threshold（）函數(shù)得出理想軟閾值為5，其無(wú)尺度拓?fù)浞植紮z驗(yàn)結(jié)果如圖2所示。其中，

k

為連接度，橫坐標(biāo)log10（

k

）表示某基因連接數(shù)的對(duì)數(shù)，縱坐標(biāo)表示該基因出現(xiàn)概率的對(duì)數(shù)，擬合系數(shù)R=0.84＞0.8，斜率為-2.29，呈線性關(guān)系，表明所選軟閾值符合無(wú)尺度拓?fù)浞植迹鶚?gòu)建的網(wǎng)絡(luò)具有較強(qiáng)的魯棒性。

如表1所示，WGCNA將過(guò)濾后的肌少癥公共數(shù)據(jù)集聚類為29個(gè)基因模塊。其中，不同顏色代表不同基因模塊，結(jié)果跟理論預(yù)期相一致，每個(gè)模塊都包含具有相似功能的多個(gè)基因。

Fig.2 Scale free topology check with soft threshold of 5圖2 軟閾值為5的無(wú)尺度拓?fù)錂z查

Table 1 Number of module genes表1 模塊基因數(shù)

圖3為識(shí)別出的29個(gè)基因模塊與樣本表型（肌少癥和健康受試者對(duì)照性狀）關(guān)系的可視化圖（彩圖掃OSID碼可見(jiàn)）。其中，紅色表示正相關(guān)，藍(lán)色表示負(fù)相關(guān)，顏色越深則相關(guān)性越強(qiáng)，可見(jiàn)Skyblue模塊（

r

=0.53，

p

=9E-04）、Turquoise模塊（

r

=0.38，

p

=0.02）和Black模塊（

r

=0.37，

p

=0.02）與肌少癥具有顯著正相關(guān)性，而與對(duì)照組呈負(fù)相關(guān)。因此，可確定這3個(gè)模塊為肌少癥的關(guān)鍵模塊。此外，所選關(guān)鍵模塊的GS和MM值也表明了篩選結(jié)果的可靠性。

實(shí)驗(yàn)從3個(gè)關(guān)鍵模塊中共識(shí)別出了86個(gè)樞紐基因，并通過(guò)GO驗(yàn)證了69個(gè)，該結(jié)果與當(dāng)前文獻(xiàn)記載肌少癥是多基因參與的復(fù)雜疾病的結(jié)論相一致。實(shí)驗(yàn)結(jié)果表明，WGCNA在識(shí)別樞紐基因時(shí)，準(zhǔn)確率達(dá)到了80.2%。表2為關(guān)鍵模塊中的樞紐基因。其中，Skyblue模塊存在17個(gè)樞紐基因，驗(yàn)證成功15個(gè)；Turquoise模塊存在68個(gè)樞紐基因，驗(yàn)證成功53個(gè)；Black模塊只有1個(gè)樞紐基因并已成功驗(yàn)證。

表3為驗(yàn)證成功的基因，可發(fā)現(xiàn)大多數(shù)基因是由多個(gè)路徑共同作用，并且根據(jù)GO_ID從GO庫(kù)查找對(duì)應(yīng)的功能發(fā)現(xiàn)，肌少癥的樞紐基因主要通過(guò)調(diào)節(jié)因子活性轉(zhuǎn)錄、肽基賴氨酸修飾、組蛋白修飾等生物學(xué)過(guò)程影響肌少癥。

Fig.3 Correlation analysis between key modules and sarcopenia圖3 關(guān)鍵模塊與肌少癥的關(guān)聯(lián)分析

Table2 Hub genes in key modules and their verification information表2 關(guān)鍵模塊中的樞紐基因及其驗(yàn)證信息

3.2 實(shí)驗(yàn)對(duì)比

為驗(yàn)證算法的有效性，本文選用目前基因篩選常用的差異基因分析法進(jìn)行對(duì)比。實(shí)驗(yàn)基于R軟件的limma包，通過(guò)差異基因分析法對(duì)實(shí)驗(yàn)數(shù)據(jù)篩選肌少癥基因，并利用GO對(duì)其進(jìn)行驗(yàn)證。通過(guò)實(shí)驗(yàn)分析表明，差異基因分析法可從16 879個(gè)基因中識(shí)別出5個(gè)肌少癥基因（SCTR-AS1、ANKRD18B、ITLN1、RPS16P1、OVOS2），并成功驗(yàn)證了ITLN1（GO：0070206，

p

=0.002）、OVOS2（GO：0010951，

p

=0.03），準(zhǔn)確率為40%。WGCNA識(shí)別出了86個(gè)肌少癥相關(guān)基因，準(zhǔn)確率達(dá)到了80.2%，相較于差異基因分析法，該方法更高效，準(zhǔn)確率約提高了40%。

Table3 Results of genefun ction validation analysis表3 基因功能驗(yàn)證分析結(jié)果

4 結(jié)論

由于肌少癥是多基因參與的復(fù)雜疾病，對(duì)發(fā)病機(jī)制的探究造成了很大的困難，利用現(xiàn)代計(jì)算機(jī)技術(shù)發(fā)掘其致病基因是一種行之有效的方法。WGCNA分析算法是目前基因共表達(dá)網(wǎng)絡(luò)研究中最為常用的生物信息學(xué)分析方法，本文首次將該算法應(yīng)用于分析肌少癥，成功篩選出86個(gè)與肌少癥相關(guān)的樞紐基因，并驗(yàn)證了69個(gè)具有真實(shí)生物信息基因，相較于傳統(tǒng)方法，該方法考慮信息更全面，準(zhǔn)確率約提高了40%。

驗(yàn)證成功的基因?yàn)樘剿骷∩侔Y的發(fā)病機(jī)制提供了新的思路，也為臨床診斷、治療靶點(diǎn)提供了參考。然而，實(shí)驗(yàn)所用樣本數(shù)較少，后期將會(huì)尋找更合適的數(shù)據(jù)以進(jìn)一步提高識(shí)別準(zhǔn)確率。