非酒精性脂肪性肝病(NAFLD) 作为慢性肝脏损伤的常见病症, 通常与血脂、血糖指标升高以及脂质代谢功能紊乱等代谢疾病的发生密切相关^[1]。随着生产生活方式的变化, 肝脏的代谢负荷增加, 正常生理功能严重受损, 导致NAFLD在我国及全球患病率剧增^{[2, 3]}。NAFLD依据疾病谱可分为非酒精性单纯性脂肪肝(NAFL) 和非酒精性脂肪性肝炎(NASH)。随着病情持续恶化, NASH又分化为可逆性肝纤维化、不可逆性肝硬化甚至原发性肝细胞癌(HCC) 等不同病理阶段^[4]。而肝纤维化作为重要病理阶段之一, 其特征为肝脏损伤、炎症因子浸润、气球样变以及不同分期的纤维样变。纤维化会导致肝脏损伤进程更快, 且不可逆地向HCC发展, 如果对其进行有效防治则可阻止终末期肝病的发生^[5]。遗憾的是, 非酒精性脂肪性肝纤维化的发病机制尚未明确, 目前公认的机制为肝脏“二次打击”学说^[6], 因此通过生物信息技术从基因芯片进行挖掘可以从分子层面为该病机制研究提供新思路。

本文利用基因表达数据库(GEO) 检索非酒精性脂肪性肝纤维化基因芯片, 筛选肝纤维化模型组与对照组差异表达基因, 进而对差异基因进行GO功能分析和KEGG通路分析, 同时构建差异基因蛋白互作网络, 明确关键基因。最后将差异基因转换为对应探针列表, 导入Connectivity Map (CMap) 筛选潜在化合物, 以期为非酒精性脂肪性肝纤维化的临床治疗提供更多选择。最终通过上述分析, 结合生物信息学分析和基因芯片探究非酒精性脂肪性肝纤维化的分子机制及潜在治疗药物提供新研究思路。

微阵列分析数据  在GEO (https://www.ncbi.nlm.nih.gov/geo/) 数据库中输入检索词“nonalcoholic liver fibrosis”, 下载GSE109345基因表达谱。GSE109345基于Affymetrix GPL19057平台(Affymetrix Illumina NextSeq 500), 由van Koppen A等提交。GSE109345数据集包含78个样本数据, 其中有48个纤维化样本和24个正常对照样本(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109345)。

数据处理及差异表达基因的筛选  将下载的样本数据进行分组并命名, 分别为纤维化组和对照组。将分组结果导入BioJupies (https://amp.pharm.mssm.edu/biojupies/) 进行处理及筛选差异基因。主成分分析(PCA) 是一种统计方法, 通常用于在RNA-seq数据集中探索生物样本的相似性^[7]。将经过处理数据进行下一步处理以识别差异基因(DEGs)。|log₂ fold change|≥2且校正后P < 0.05作为差异基因的筛选条件。

GO功能富集和KEGG通路分析  GO分析是功能富集研究中常用的方法, 基因功能包括生物过程(BP)、分子功能(MF) 和细胞成分(CC)。KEGG是一个在生物信息学分析中广泛使用的数据库, 它存储大量关于基因组、生物途径和疾病的数据。使用DAVID (https://david.ncifcrf.gov/) 对本研究中的DEGS进行GO注释分析和KEGG通路分析。其中GO分析中筛选条件为P < 0.05, 基因计数 > 5; KEGG通路分析P < 0.05。

PPI网络构建  互作基因检索数据库(https://string-db.org/) 可以将差异基因之间的关系进行可视化处理。将筛选得到的差异基因提交到String数据库, 选择研究物种为“Homo sapiens”, 对潜在蛋白的相互作用进行分析, 设置互作分值为高可信(high confidence: 0.700), 隐藏蛋白网络中没有互作关系的节点, 构建疾病-蛋白互作网络, 导出TSV格式文件。

CMap分析  CMap是药物基因组学研究领域的重要数据库, 收录1 300多种化合物, 该数据库是基于相关疾病基因表达转录组分析的基础, 揭示疾病基因与潜在化合物之间的关系^{[8, 9]}。将上述筛选差异基因经R语言处理后转换为对应探针列表, 上传CMap数据库(https://portals.broadinstitute.org/cmap/) 筛选潜在化合物。

PCA分析显示数据的前三个主要成分(PC) 的交互式三维散点图。每个点代表一个RNA-seq样品, 具有相似基因表达谱的样品在三维空间中更接近(图 1)。因此, 对照组与模型组的样品表达相互独立, 进行下一步分析。将对照组和模型组进行差异基因表达分析, 共筛选得到109个差异表达基因, 包括70个上调基因和39个下调基因, 对结果进行可视化分析, 绘制聚类图及火山图(图 2)。

Figure 1

Figure 1

Figure 1 Principal component analysis result graph

Figure 2

Figure 2

Figure 2 Differential gene clustering heat map and volcano map. A: Differential gene heat map between normal group and liver fibrosis patients; B: Volcanic map of different genes between normal group and liver fibrosis patients

运用DAVID对本研究中差异基因进行GO、KEGG分析。通过对109个差异表达基因进行分析, 得到结果如下(图 3)。GO富集分析显示, 差异基因主要富集的生物过程(BP) 为蛋白激酶B信号转导的正调控(positive regulation of protein kinase B signaling)、基因表达的正调控(positive regulation of gene expression)、对二苯乙烯类化合物的反应(response to stilbenoid)、环氧酶P450途径(epoxygenase P450 pathway) 和细胞对肿瘤坏死因子的反应(cellular response to tumor necrosis factor) 等; 对差异基因的组分分析(CC) 主要集中在胞外区(extracellular region)、细胞外隙(extracellular space)、细胞表面(cell surface)、细胞外基质(extracellular matrix) 和细胞器膜(organelle membrane) 等; 差异基因的分子功能(MF) 有转运活性(transporter activity)、小分子结合(small molecule binding)、信息素结合(pheromone binding)、铁离子结合(iron ion binding) 和血红素结合(heme binding); KEGG通路分析结果显示, 差异表达基因主要富集的通路有视黄醇代谢通路(retinol metabolism)、类固醇激素合成通路(steroid hormone biosynthesis) 和花生四烯酸代谢通路(arachidonic acid metabolism), 各通路富集基因及数目如表所示(表 1)。

Figure 3

Figure 3

Figure 3 GO enrichment analysis results of differential genes

表1(Table 1)

Table 1 KEGG pathways analysis Term Description Gene Count mmu00830 Retinol metabolism CYP2B9, CYP26A1, CYP4A12B, CYP2B13, CYP2A4, CYP2B10, CYP2C40 7 mmu00140 Steroid hormone biosynthesis CYP2B9, CYP2B13, HSD3B5, SULT1E1, CYP2B10, CYP2C40 6 mmu00590 Arachidonic acid metabolism CYP2B9, CYP4A12B, CYP2B13, CYP2B10, CYP2C40 5

Table 1 KEGG pathways analysis

将差异基因输入String平台后, 经处理后共得到104个节点, 173条连线, 剔除无互作关系的节点后绘制PPI图(图 4)。依据图中所示, 与非酒精性脂肪性肝纤维化发生发展较为密切的基因主要为金属蛋白酶组织抑制因子1 (TIMP1)、CC趋化因子配体2 (CCL2)、信号调节蛋白B1 (SIRPB1) 和细胞色素P450 (CYP) 等。

Figure 4

Figure 4

Figure 4 Differentially expressed gene protein interaction network

CMap中的mean值、富集得分范围为: -1~1, 分数正值表示该药物与疾病发展呈正相关, 负值表示该药物能够抑制疾病发展, 绝对值越大则表明药物与疾病相关性越高^[9]。将处理后的差异基因探针列表提交至CMap数据库, 筛选出关联强度居前5位的负相关化合物, 吡格列酮(pioglitazone)、槟榔碱(arecoline)、米多君(midodrine)、曲格列酮(troglitazone) 和小白菊内酯(parthenolide) 等化合物对非酒精性脂肪性肝纤维化差异表达基因具有抑制作用, 转换后探针ID列表及筛选得到的化合物如下(表 2、3)。

表2(Table 2)

Table 2 Differential gene corresponding gene probe ID list Down Up Gene symbol ID Gene symbol ID CADM4 214267_s_at A4GNT 221131_at SRGAP3 209794_at ADAM8 205180_s_at CLIC3 219529_at ANXA2 210427_x_at CYP26A1 206424_at VCAN 215646_s_at GRM8 216255_s_at APOA4 206894_at SQLE 213562_s_at ARHGAP22 206298_at SLC13A2 207567_at CAPG 201850_at PITX3 208277_at CCL2 216598_s_at CCR2 207794_at CD63 200663_at CHST11 219634_at CLEC4E 219859_at CLEC7A 221698_s_at COL1A1 217430_x_at CORIN 220356_at DBP 209782_s_at DCLK1 205399_at DEFB1 210397_at EGR2 205249_at FGF21 221433_at GPNMB 201141_at HAO2 220801_s_at IL1RN 216245_at ITGAX 210184_at LCN2 212531_at LGALS3 208949_s_at LY6D 206276_at MMP12 204580_at MMP13 205959_at NRG1 206237_s_at NTRK2 221796_at PLS1 205190_at RGS16 209325_s_at SERPINE1 202627_s_at SLC35F2 218826_at SMPD3 219695_at SPRR1A 214549_x_at SULT1E1 219934_s_at SYBU 218692_at THBS1 201109_s_at TIMP1 201666_at TTC39A 217101_at

Table 2 Differential gene corresponding gene probe ID list

表3(Table 3)

Table 3 The top 5 potential treatments for liver fibrosis in the CMap list CMap name Mean n Enrichment P Specificity Percent non-null Pioglitazone -0.420 11 -0.571 0.000 74 0.042 4 63 Midodrine -0.298 5 -0.784 0.000 86 0.030 1 60 Arecoline -0.259 4 -0.793 0.003 68 0.012 4 50 Troglitazone -0.267 16 -0.379 0.014 62 0.013 9 50 Parthenolide -0.506 4 -0.656 0.032 35 0.220 7 75

Table 3 The top 5 potential treatments for liver fibrosis in the CMap list

随着工业化进程以及生活方式的转变, 作为主要代谢性疾患之一的NAFLD患病率呈逐年上升趋势, 现已成为全球第一大慢性肝病。同时, 受到儿童肥胖率的升高、2型糖尿病患病率的攀升以及人口老龄化等因素的影响, 预计未来NAFLD的患病率以及进展期肝病患者的比例将持续增加^[10]。研究表明, 晚期纤维化的NAFLD患者死亡率明显高于早期患者, 因此对于非酒精性脂肪性肝纤维化的机制研究及临床防治显得极为重要^[11]。随着基因芯片技术的发展, 一方面利用生物信息技术能够精确阐明疾病发病机制, 另一方面能够借助生物信息的数据优势作为药物研发的新途径。

本研究利用GEO数据库丰富的数据信息, 以“NASH”、“nonalcoholic fatty liver fibrosis”作为关键词检索并筛选芯片数据, 获取GSE109345数据集数据。将数据分为纤维化组和对照组, 利用生物信息学分析筛选得到差异表达基因109个, 其中上调表达70个, 下调表达39个。进一步对得到差异基因进行功能富集分析(GO) 及通路分析(KEGG), 明确差异基因所参与及富集的功能和通路, 同时构建差异基因蛋白互作网络, 将其进行可视化分析, 以明确与非酒精性脂肪性肝纤维化相关的主要基因。最后利用CMap筛选潜在化合物, 为非酒精性脂肪性肝纤维化的临床防治提供候选治疗药物。

研究发现, 慢性肝脏疾病均伴随代谢紊乱症, 非酒精性脂肪性肝纤维化的发生与脂质代谢紊乱、胰岛素抵抗(insulin resistance, IR) 和慢性氧应激等密切相关。相关临床及实验研究表明, 表 3中的潜在药物通过改善代谢紊乱达到治疗肝纤维化的效果。

IR能够导致肝脏脂肪过量沉积, 诱发游离脂肪酸(FFA) 过量释放, 导致肝星状细胞(HSC) 活化和细胞外基质(ECM) 在肝脏堆积, 加速肝脏纤维因子的合成^[12]。吡格列酮是临床治疗糖尿病的常用药, 能明显改善IR, 实验证明应用吡格列酮干预可明显提高IGT-OLETF大鼠对外源性胰岛素的反应性, 改善其高胰岛素血症, 使胰岛素敏感指数基本恢复正常。同时甘油三酯(TG) 和FFA水平降低, 提示吡格列酮能够降低肝脂含量, 减轻肝脏脂质沉积, 改善脂肪组织正常代谢功能^{[13, 14]}。为进一步明确吡格列酮治疗纤维化的具体机制, 通过构建吡格列酮干预非酒精性脂肪性肝纤维化模型, 结果表明吡格列酮通过调控Toll样受体4 (TLR4)/细胞核因子κB (NF-κB) 及转化生长因子-β1 (TGF-β1)/Smad3蛋白(Smad3) 信号通路, 抑制TLR4、IκB激酶β (IKK-β)、NF-κB等基因的表达, 进而抑制HSC活化、肝细胞炎症、纤维化相关信号通路, 延缓或阻止非酒精性脂肪性肝纤维化的进展^{[15, 16]}。

米多君是α1肾上腺素能受体激动剂, 研究发现在常规治疗基础上加用米多君可提高肝硬化患者的生存率, 有助于控制肝硬化造成的腹水^[17]。肝纤维化作为肝硬化进展过程中的重要病理阶段, 若能在肝纤维化早期进行干预, 能有效阻止终末期肝病的发生。目前尚未有米多君治疗肝纤维化的临床报道, 但CMap结果提示米多君具有干预纤维化进展的可能。

槟榔碱是槟榔所含的主要生物碱, 研究发现其对胰岛B细胞分泌功能具有直接保护作用, 能够促进胰岛素分泌, 维持肝细胞正常生理功能。文献报道槟榔碱能通过“外源性化合物感受器”组成性雄甾烷受体(CAR)、孕烷x受体(PXR) 改善2型糖尿病模型大鼠肝糖代谢紊乱, 具体机制为槟榔碱通过CAR、PXR直接作用于叉头状转录因子O1 (FoxO1)、过氧化物酶体增殖物激活受体γ共激活因子-1α (PGC-1α), 调节肝糖代谢关键酶磷酸烯醇式丙酮酸羧激酶基因(PEPCK) 和葡萄糖-6-磷酸酶(G6Pase) 基因表达; 同时通过NF-κB降低2型糖尿病大鼠肝炎性因子肿瘤坏死因子-α (TNF-α) 和白细胞介素-6 (IL-6) 表达, 改善肝胰岛素敏感性^{[18, 19]}。

曲格列酮是一种噻唑烷二酮类胰岛素增敏剂, 临床主要用于治疗2型糖尿病和胰岛素抵抗。研究提示, 曲格列酮毒性作用机制可能与其诱导活性氧簇(ROS) 生成有关, 过量ROS不仅可以改变基因和蛋白表达, 影响细胞的功能与表型, 还可以攻击细胞内的核酸、蛋白质和脂质等重要成分, 抑制线粒体生物合成, 进一步加剧线粒体功能障碍, 最终引发细胞的凋亡^{[20, 21]}。

小白菊内酯是从中草药小白菊中提取到的一种倍半萜内酯类化合物, 具有抗炎和抗肿瘤等药理作用。在高脂饮食诱导大鼠NAFLD模型中, 应用小白菊内酯干预后胆固醇(TC)、低密度脂蛋白(LDL) 水平降低, 高密度脂蛋白(HLD) 水平升高, 提示其可以降低高血糖, 改善胰岛素敏感性, 调节肝脏血脂异常, 增强细胞的抗氧化能力, 延缓NAFLD进展^[22]。此外, 小白菊内酯能够与NF-κB通路结合发挥抗炎作用, 值得注意的是该通路与HSC的活化、ECM的合成密切相关。因此, 小白菊内酯可能通过与NF-κB通路结合发挥抗纤维化的药理作用^[23]。

综上所述, 本实验运用生物信息学结合基因芯片, 探究非酒精性脂肪性肝纤维化的差异表达基因, 并对其进行初步分析, 得到差异基因主要参与的信号通路及其功能, 进一步挖掘其潜在治疗药物。为未来肝纤维化的临床治疗和机制研究提供理论依据, 但仍需要更加严谨的实验对其加以验证。

作者贡献: 本文设计思路由杨婧、王锐提供; 赵耀伟负责靶点筛选、通路富集、网络构建及CMap药物筛选; 牛捷负责数据可视化处理。

利益冲突: 本文不存在利益冲突。