利用加权基因共表达网络分析识别胆道闭锁相关的枢纽基因

周红，吴梅，李鑫

（1.宜宾市第一人民医院儿科，宜宾 644000;2.乌鲁木齐儿童医院普外科，乌鲁木齐 830000）

胆道闭锁（biliary atresia，BA）是一种婴幼儿期常见的极为严重的疾病，特点为发生于生后3 个月内的部分或者全部肝外胆道完全性纤维化梗阻，如果不治疗，将发展为肝硬化、肝衰竭，不经治疗的平均生存期在1 年左右，主要的治疗方式是肝门肠吻合术（hepatoportoenterostomy，HPE）和肝移植[1-4]。即使在出生后45 d 内进行HPE，2 岁时未进行肝移植的生存率仅为65.5%[5]。BA 面临早期诊断非常困难、手术治疗效果不佳、肝移植后并发症和免疫抑制管理困难等多种挑战。目前BA 的发病机制、病因仍不清楚[6-7]，近年来专注于鉴定BA 相关的遗传和免疫因子作用的研究成为热点。且BA 是新生儿胆汁淤积最常见的原因，BA 与非BA 的胆汁淤积性疾病在临床表现和实验室检查上非常相似，除手术探查外的方法难以鉴别[8]，因此寻求一种无创的方法以鉴别这两种疾病非常有必要。

加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）是一种分析多个样本基因表达模式的分析方法，将具有类似表达模式的基因整合到同一的模块，并分析模块与特定性状或表型之间的关联关系[9-11]。本研究拟利用WGCNA筛选出与BA 相关的基因模块及枢纽基因，以期为BA的发病机制及临床诊治带来新的理解。

1 材料与方法

1.1 数据集的收集 从美国国立生物信息技术中心（National Center for Biotechnology Information，NCBI）的基因表达综合数据库（Gene Expression Omnibus，GEO）下载GSE46960数据集，芯片测序平台为GPL6244。从该数据集中选取64 例诊断为BA 的肝脏样本（BA 组），14 例非BA 的胆汁淤积性疾病肝脏样本（Non_BA 组）以及7 名已故无肝脏疾病捐赠者的肝脏样本（Normol 组）数据进行分析。采用R 软件进行数据的处理。

1.2 加权共表达基因网络的构建 以表达量为筛选条件，选择平均表达量最高的5 000 个基因在R 软件的WGCNA 包构建加权基因共表达网络。为确保无尺度网络分别计算所有基因之间的成对Pearson 相关系数，设置阈值进行筛选，以将成对相关矩阵转换为邻近相关矩阵。采用动态混合剪切树算法标准，设置每个基因模块最少基因数目为30，并且依次计算每个模块的特征向量值，然后对模块进行聚类分析，将距离较近的模块合并为新的模块。

1.3 模块特征相关性的计算 通过WGCNA 算法计算模块基因与疾病分组表型之间的相关性，相关性的强度通过热图反应出来。当P＜0.05 时，认为单个模块与表型显著相关。选择与BA 相关系数最高的模块作为关键模块。

1.4 枢纽基因的筛选 计算每个共表达模块与基因特征值的皮尔森相关系数。选取模块身份（module membership，MM）＞0.8 且基因显著性（gene significance，GS）＞0.65 的基因为枢纽基因。

1.5 关键模块的富集分析 提取关键模块中的基因，利用R 软件ClusterProfiler 包进行基因本体论（Gene Ontology，GO）功能注释和《京都基因和基因组百科全书》（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路富集分析。P＜0.05 为差异具有统计学意义。

2 结果

2.1 样本数据筛选 经过归一化和基因ID 转换生成18 762 个基因和85 个样本的表达矩阵。选取平均表达量最高的5 000 个基因构建基因共表达模块。数据集处理后进行离群值检测，分析结果示该数据集中无明显离群值，故直接进行下一步分析，将基因聚类模块与临床分组表型进行分析（图1）。

图1 基于欧几里得距离的样本聚类树图Fig 1 The sample cluster tree diagram based on Euclidean distance

2.2 软阈值计算 利用WGCNA 构建共表达网络需先确定软阈值（soft-thresholding power），本研究设定1~20，对每个软阈值计算无标度拓扑拟合指数R2（scale free topology fitting index R2）和平均连通性（mean connectivity）。软阈值对应的R2＞0.8 且平均连通性接近零，说明此网络符合无尺度条件。因此，本研究选取软阈值β = 9（图2）。

图2 软阈值参数的确定Fig 2 Determination of soft threshold parameters

2.3 构建基因共表达模块 使用动态混合切割构建层次聚类树，树上的每一片叶子代表一个基因，具有相似表达数据的基因靠在一起，形成树的一个分支，代表一个基因模块，生成了20 个模块（图3）。

图3 模块构建Fig 3 Construction of modules

2.4 共表达模块之间的相关性分析 将层次聚类得到的基因模块进行相关性分析，观察模块与模块之间的关联，颜色越深则表示模块之间相关性越高（图4）。

图4 特征基因树图和特征基因邻接图Fig 4 Characteristic gene tree diagram and characteristic gene adjacency diagram

2.5 共表达模块的关联分析 根据拓扑重叠程度，选择全部的基因制作热图（图5）。

图5 共表达网络的热图可视化Fig 5 Thermal graph visualization of co-expression network

2.6 计算BA 疾病相关的模块相关性 对于每一个模块，计算基因表达与不同疾病组之间的相关性。与BA 相关性最大的模块是黄色模块，相关系数为0.69（P＜0.05），同时黄色模块与非BA 的胆汁淤积性疾病及正常组是负相关关系（图6），故将黄色模块识别为关键模块。

图6 模块基因与临床表型相关性图Fig 6 Correlation between modular genes and clinical phenotypes

2.7 识别关键模块内的枢纽基因 分析黄色模块内基因连通性（MM）及基因显著性（GS）之间的相关性，二者相关性良好，呈明显线性相关（图7）。以模块身份MM＞0.8 且基因显著性GS＞0.65 作为筛选条件，共筛选出16 个枢纽基因，包括基质金属蛋白酶7（matrix metallopeptidase 7，MMP7）、重组人分泌型磷蛋白1（secreted phosphoprotein 1，SPP1）、多功能蛋白聚糖（versican，VCAN）、含V-Set 域T 细胞激活抑制因子1（V-set domain containing T cell activation inhibitor，VTCN1）、重组人光蛋白聚酶（lumican，LUM）、G-蛋白信号转导调节子4（regulator of G protein signaling 4，RGS4）、上皮细胞黏附分子（epithelial cell adhesion molecule，EPCAM）、超长链脂肪酸延伸酶7（elongation of very long fatty acid elongase 7，ELOVL7）、MET 转录调节因子MACC1（MET transcriptional regulator MACC1，MACC1）、层黏连蛋白γ2（laminin subunit gamma 2，LAMC2）、染色体12开放阅读框75（chromosome 12 open reading frame 75，C12orf75）、IV 型胶原蛋白基因α1（collagen type IV alpha 1 chain，COL4A1）、波形纤维蛋白（vimentin，VIM）、跨膜蛋白（transmembraneprotein156，TMEM156）、胸腺素β10（thymosin beta 10，TMSB10）、肌球蛋白重链9（myosin heavy chain 9，MYH9）（表1）。

表1 黄色模块中的枢纽基因Tab1 Hub gene in yellow module

图7 黄色模块基因相关性散点图Fig 7 Scatter plot of module eigengenes in the yellow module

2.8 模块的富集分析 黄色模块GO 富集结果表明，生物学过程（biological process，BP）主要富集在细胞外基质组织和细胞外结构组织，细胞成分（cellular component，CC）主要富集于含胶原的细胞外基质，分子功能（molecular function，MF）主要富集于细胞黏附分子。黄色模块KEGG 富集结果示，BA 的发生可能与色氨酸代谢、甘油磷脂代谢、过氧化物酶体增殖物激活受体（peroxisome proliferator-activated receptor，PPAR）信号通路等密切相关。见图8。

图8 GO 和KEGG 富集气泡图Fig 8 Bubble diagrams of GO and KEGG enrichment

3 讨论

BA 是发生在婴幼儿期的一种严重的先天性疾病，由于肝门附近的胆道系统狭窄、闭锁或缺如，结合胆红素排入肠道受阻，最终导致肝脏胆汁性硬化、肝衰竭，甚至死亡，生后1 个月内行HPE 治疗肝脏存活率明显增高，但诊断该疾病需排除非BA 的胆汁淤积性疾病，如果不能区分BA 与非BA 的胆汁淤积性疾病可能导致手术治疗的延迟[1-2]。既往的研究通过差异分析的方法鉴定BA 的特征性基因，可能忽略一些具有重要功能但差异表达不明显的基因，也可能混杂与非BA 的胆汁淤积性疾病的相关基因。WGCNA 是一种无监督的层次聚类方法，能够识别与表型相关的“基因模块”，即可以整合在基础生物学途径中的高度互连的基因。与基因差异分析相比，WGCNA 划分的基因模块具有明显的生物学意义。结合临床信息分析，可以确定与疾病发病机制相关的重要模块和潜在枢纽基因[12]。为进一步分析BA 的发生、发展机制，并将非BA 的胆汁淤积性疾病加以鉴别，本研究采用WGCNA 分析BA、非BA 的胆汁淤积性疾病和无肝脏疾病肝脏转录组数据。

本研究构建了20 个基因共表达模块，鉴定了1 个与BA 高度相关的基因模块即黄色模块，其与非BA 的胆汁淤积性疾病和无肝脏疾病无明显相关性。本研究对关键模块进行GO 及KEGG 富集，发现与BA 显著相关的基因主要富集在细胞外基质（extracellular matrix，ECM）、细胞间黏附分子、色氨酸代谢、甘油磷脂代谢、PPAR 信号通路。ECM 含有大量信号分子，参与调节细胞的迁移、增殖、凋亡或分化，细胞不断通过合成、降解、重新组装和化学改造来重建和改造ECM，以保持组织平衡，不良的ECM 可能加剧疾病进展[13]。PPAR 受体激动剂能减少肾小球硬化和肾小管间质损伤[14]。

本研究根据模块身份和基因重要性识别出16 个BA 相关的枢纽基因，包括MMP7、SPP1、VCAN、VTCN1、LUM、RGS4、EpCAM、ELOVL7、MACC1、LAMC2、C12orf75、COL4A1、VIM、TMEM156、TMSB10和MYH9。有研究显示MMP7 是一种肝星形细胞分泌的蛋白酶，与BA 的肝纤维化与组织重塑相关[15]。血清蛋白组学研究提示MMP7 测定具有高灵敏度和特异性，可区别BA 与非BA 的胆汁淤积性疾病，可作为BA的可靠的生物标志[16-18]。这些提示MMP7 直接参与了BA 发病，但MMP7 对于诊断BA 的截断值尚未统一（1.43~52.85 ng/mL）[18]，因此MMP7 的临床应用仍需要更大规模的研究。SPP1 在胆汁异常患者中被过度表达，但它的灭活并没有影响表型的发展[19]，这表明SPP1 可能涉及肝脏对胆汁损伤的一般生理反应，但不一定与BA 的发病机制有关。在BA 患者及BA 动物模型的肝脏中均检测到未成熟的EpCAM，且EpCAM 阳性的不成熟胆道上皮细胞在短时间内迅速增生，EpCAM 可能是BA 患儿肝纤维化进展迅速的主要原因[20]。文献报道ELOVL7 能产生对于构建正常功能的人巨细胞病毒包膜必需的饱和超长链脂肪酸[21]，这与BA 的病因可能是病毒感染的说法一致。已有研究显示16 个枢纽基因中的部分基因在BA 发病中可能发挥重要作用，表明本研究采用的鉴定枢纽基因的方法精确可靠。VCAN、VTCN1、LUM、ELOVL7、MACC1、LAMC2、C12orf75、COL4A1、VIM、TMEM156、TMSB10 及MYH9 在BA发病机制中的研究，枢纽基因及富集分析所示功能是否能减缓甚至逆转BA 的进展，有待进一步深入研究。

综上所述，基于WGCNA 识别了与BA 发病相关的1 个关键模块及16 个枢纽基因，新发现了12 个可能和BA 的发病密切相关的基因，值得进一步实验验证。