溃疡性结肠炎(ulcerative colitis,UC)是炎症性肠病的两种主要形式之一,它是一种慢性非特异性肠道疾病,主要表现为直肠和结肠的长期炎症和溃疡,其病变通常从直肠连续分布到近端结肠,且局限于肠黏膜和黏膜下层。UC病程表现为间歇性,活动期和非活动期交替出现,临床症状主要为腹痛、腹泻、便血和体重减轻等,并且由于病变持续反复的复发使得更加容易恶变为结肠癌。UC的发病机制目前尚未完全清楚,可能与免疫反应失调、肠道微生物紊乱、环境和遗传因素有关
[1]。近年,UC的发病率和患病率不断上升,并且由于复杂的发病机制和缺乏有效的治愈手段,这使得寻求更有效的临床诊断、监测和治疗方法日趋重要。
内质网应激(endoplasmic reticulum stress,ERS)是引起炎症性疾病发生、发展的重要因素。其生物特征是生理病理情况下的内质网中错误折叠蛋白质不断积累,从而使其蛋白折叠反应负荷增加,导致内质网稳态失调并引起内质网应激反应
[2]。研究表明在UC肠道病变过程中,持续的ERS会损害细胞功能,并转换为凋亡的适应机制,以清除不可逆的损伤细胞。过度的细胞凋亡会导致肠上皮细胞无法修复,从而破坏肠上皮的完整性并引发炎症反应
[3,4]。所以,探索ERS上游基因的调控网络对于调节UC肠道炎症病变具有重要意义。
MicroRNA(miRNA)是17~25个核苷酸的内源性单链非编码 RNA。它们靶向mRNA的3’非翻译区域,并根据互补碱基的程度抑制或降解靶基因。miRNA作为重要的基因表达调控因子,参与调控大多数人类基因的表达,在包括UC在内的许多自身免疫性疾病的发病机制中发挥着重要作用
[5]。近年来的研究表明,miRNA的异常表达在疾病发病机制的早期就发生了变化。此外,这种异常表达与UC和结肠炎相关癌症的疾病活动高度相关,几乎涉及UC所有关键的发病机制,包括调节肠道屏障功能、微生物菌群、免疫反应以及遗传易感性等
[6,7]。
孟德尔随机化(Mendelian randomization,MR)是一种使用遗传变异作为暴露的工具变量,来研究暴露因素与结果表型之间的因果关系。与传统的观察性研究相比,MR不太容易受到混淆因素和反向因果关系的影响
[8]。因此,MR分析可以被认为是一种替代的随机对照试验方法
[9]。大规模全基因组关联研究(genome‑wide association studies,GWAS)已经确定了许多与不同性状相关的遗传变异,这为MR提供了强大的数据来源
[10]。目前,关于miRNA与UC风险关联研究未见报道。由此,通过利用MR分析鉴定UC中miRNA的风险因果关系,并结合生信分析探索其靶向结合的内质网应激相关特征基因的潜在生物功能和免疫机制,为研究UC的有效诊断及靶向治疗提供一种新思路。
1 资料和方法
1.1 MR分析
本研究中的暴露数据来源于美国心肺研究院的一项microRNA表达数量性状位点的全基因组鉴定研究,并从中获取miRNA eQTL数据,该研究全面调查了5 000多个无血统关系的成熟人类miRNA的表达水平
[11]。另一方面,UC的遗传关联数据则来源于GWAS数据库(
https://gwas.mrcieu.ac.uk),GWAS ID: finn‑b‑ULCERNAS,包括212 507个样本和16 380 457个single nucleotide polymorphisms(SNP)位点。
研究用于MR分析的工具变量需满足3大基本假设:首先,遗传变异的SNP与暴露之间显著相关,同时确保每个遗传变异间彼此独立。设定
P<5e
-08;kb=1 000,
r2 <0.01去除存在连锁不平衡;用
F值排除弱工具变量偏倚,剔除
F<10的弱变量,计算公式为beta
2/S
x2,其中beta为遗传变异对暴露的效应值,S
x为遗传效应的标准误。其次,遗传变异仅通过暴露影响结果,而不通过其他生物学途径(即无水平多效效应)。当遗传变异在结局中
P>5e
-08,认为与结局无直接关联。最后,作为暴露的工具变量提取的遗传变异独立于与所选暴露和结果相关的混杂因素
[12]。
由于所有的遗传变异不一定都是有效的工具变量,因此采取多种方法进行MR分析。包括MR‑Egger回归,加权中位数,加权众数法,简单众数法和逆方差加权法(inverse variance weighted,IVW)5种方法
[13]。其中当遗传变异间不存在异质性和多效性时,IVW较其他分析方法有更好的统计性能,故以IVW的结果为主。IVW以结局方差的倒数为权重进行拟合,回归时不考虑截距项,当存在多效性时会引起偏倚,故此时用MR‑Egger为主。计算5种方法的beta或
OR值,方向一致则证明结果稳健。另外,异质性检验采用MR‑Egger和IVW法,当
P>0.05时则表示不存在异质性,相反当
P<0.05时,表示SNPs之间存在异质性;敏感性分析采用逐一排除法,以探讨单个SNPs对因果关联的影响。多效性分析则利用MR pleiotropy test函数,当
P<0.05时,表示存在多效性;
P>0.05则表示不存在多效性。
1.2 UC差异基因分析和风险miRNA的内质网应激相关基因鉴定
从GeneCards数据库(
https://www.genecards.org)获得ERS相关基因,通过关键词“内质网应激”和相关性评分≥10分进行筛选。
从GEO基因表达数据库(
https://www.ncbi.nlm.ni h.gov/geo/)中选取了6个关于UC患者及其对照组肠黏膜组织标本的高通量测序数据集,共149个UC患者和79个健康对照组。其中5个芯片数据集为GSE9452、GSE38713、GSE48958、GSE65114和GSE87466用于差异基因的分析,1个芯片数据集GSE13367作为最后的基因表达验证集。利用R软件的“limma”包对单个芯片数据集进行矫正,然后将矫正后的5个数据集利用R包“sva”进行多芯片数据合并。接下来使用Abel等
[14]提出的方法,用principal components analysis(PCA)分析消除数据集之间因实验操作、收集方法等差异引起的批次效应。然后通过利用“limma”包分析UC患者和健康对照组之间的差异表达基因(differentially expressed genes,DEGs),以“
P<0.05和差异倍数(Foldchange,FC)的|logFC|>0.585”作为筛选差异的阈值,得到具有显著的上调或下调的差异基因。随后将miRNA靶基因、内质网应激相关基因和UC差异表达基因取交集,得到miRNA靶向的差异内质网应激相关基因,最后通过R包对上述结果进行可视化。
1.3 Lasso回归和随机森林机器学习方法鉴定关键基因
将上述差异内质网基因整合到Lasso回归算法的分析中,该算法使用R包“glmnet”执行。另一个是随机森林机器学习方法,使用R包“randomForest”执行。接下来将两个方法得到的关键基因取交集,并通过receiver operating characteristic(ROC)曲线分析每个基因对样本的准确性,从而检测关键基因的诊断效能。最后利用R包对上述结果进行可视化。
1.4 关键基因的免疫浸润评估、GSVA分析以及表达验证
使用CIBESORT算法分析芯片数据集中22种免疫浸润细胞的比例,并通过R软件“tidyverse”包将5个关键基因与免疫浸润细胞进行Spearman相关性分析
[15]。另外,用R软件“GSVA”包对关键基因进行GSVA‑KEGG富集分析。最后,将关键基因在芯片数据集GSE13367中进行表达水平验证并进行可视化。
2 结果
2.1 MR分析结果
根据工具变量筛选标准,从原始数据(miRNA eQTL数据)中共收集了9 613个SNP作为工具变量,并剔除
F<10的弱工具变量偏倚。经MR分析,结果显示,hsa‑miR‑130b‑3p与UC存在因果关系(
表1)。5种分析方法所得的因果效应方向一致,其中IVW法的MR结果提示
P<0.05,
OR<1,表明hsa‑miR‑130b‑3p和UC存在负向因果关系(图
1A、
1B)。此外Cochran′s Q检验的
P>0.05,表明纳入的SNPs无明显异质性(
图1C)。对于水平多效性,MR‑Egger法截距值为-0.012(
P=0.755),表明纳入的SNPs无明显水平多效性,说明工具变量不通过暴露以外的途径影响结局。对于敏感性分析,在留一法检验显示去除任意SNP后结果仍然稳定(
图1D)。
2.2 UC差异基因分析和风险miRNA的内质网应激相关基因鉴定结果
首先对5个芯片数据集进行PCA分析,结果显示,批次矫正前不同实验数据集是分散的,说明存在批次效应(
图2A),而矫正后是随机分散且各个数据集中在一个范围内的,说明消除了批次效应对后续分析的影响(
图2B)。接下来对多芯片联合数据集进行DEGs分析,结果显示,UC患者与健康对照组之间的DEGs有1 060个,其中上调的有649个,下调的有411个(
图3)。接下来通过数据库对hsa‑miR‑130b‑3p的靶基因进行预测,共预测到1 286个靶基因。GeneCards数据库得到787个内质网应激相关基因。然后将UC差异基因、hsa‑miR‑130b‑3p靶基因和内质网应激相关基因取交集,最终筛选出5个hsa‑miR‑130b‑3p靶向差异的ERS基因(
图4A),其中在UC芯片数据集中genes encoding endothelin‑1 (
EDN1)、peroxisome proliferator‑activated receptor gamma (
PPARG)表达下调,而x‑box binding protein 1 (
XBP1)、mitogen‑activated protein kinase kinase kinase 5 (
MAP3K5)、leucine rich repeat kinase 2 (
LRRK2)表达上调(
图4B)。
2.3 机器学习对关键基因的鉴定结果
通过Lasso回归和随机森林机器学习两种方法鉴定UC中hsa‑miR‑130b‑3p靶向的差异ERS关键基因。结果显示,Lasso回归和随机森林算法同时鉴定出
XBP1、
PPARG、
LRRK2、
MAP3K5、
EDN1作为UC的关键基因(
图5)。随后对这5个关键基因构建ROC曲线,并计算曲线下面积(area under curve,AUC),其中AUC大于0.7表示良好的诊断性能。分析结果显示,5个关键基因的AUC均大于0.7(
图6A),说明这些特征基因去识别对照组和疾病组的样品准确性比较高。另外,联合使用5个关键基因预测疾病模型时AUC=0.961,95%
CI:0.934~0981(
图6B),表明联合5个关键基因预测UC与对照组时,相较于单个基因更具有敏感性和特异性。
2.4 关键基因GSVA分析
利用GSVA分析5个关键基因在高、低表达时的生物功能和通路。结果显示(
图7),
EDN1高表达时富集在ErbB信号通路、PPAR信号通路、甘油磷酸酯代谢等,低表达时富集在原发性免疫缺陷、肠道免疫网络、碱基切除修复等;
PPARG高表达时富集在氨基酸代谢和降解途径、PPAR信号通路、过氧化物酶体等,低表达时富集在肠道免疫网络、原发性免疫缺陷、细胞黏附分子、趋化因子信号通路等;
XBP1高表达时富集在白细胞迁移、趋化因子信号通路、大肠杆菌感染等,低表达时富集在糖代谢、氨基酸、磷脂代谢途径等;
LRRK2高表达时富集在趋化因子信号通路、细胞因子及其受体作用网络、细胞黏附分子等,低表达时富集在丙酮酸、磷酸、磷脂代谢等;
MAP3K5高表达时富集在霍乱弧菌感染、糖胺聚糖降解、亚麻酸代谢等;低表达时富集在氨基酸的合成、降解、代谢等过程。
2.5 免疫浸润评估
将关键基因进行免疫浸润分析,结果显示(
图8),5个关键基因均与T细胞亚群和B细胞亚群密切相关。另外,在中性粒细胞中
XBP1、
LRRK2、
MAP3K5呈现正相关,而
PPARG与中心粒细胞负相关;在巨噬细胞亚群中,
XBP1、
LRRK2、
MAP3K5与M2型巨噬细胞负相关且与M0型巨噬细胞正相关,而
PPARG则与之相反。除此之外,个别关键基因还与NK细胞、浆细胞、肥大细胞呈现关联。
2.6 关键基因的表达验证
为了评估预测结果的可靠性,利用独立检验数据集GSE13367再次验证5个关键基因在UC组织中的表达水平。结果显示,与正常对照组相比,
EDN1、
PPARG在UC组织中显著低表达;而
XBP1、
LRRK2、
MAP3K5显著高表达,且差异具有明显统计学意义(
图9),这与本课题组先前预测的结果完全一致。
3 讨论
UC作为一种难以治愈且反复发作的炎症性疾病,以慢性和异质性表现为特征,由环境、基因组、微生物和免疫因素相互作用诱发,这些相互作用导致了其发病机制的复杂性。为了准确定义UC的复杂交互机制,需要新的概念和工具来实现系统方法,从而揭示疾病网络中的关键参与者,精确定位炎症的驱动中心。通过利用先进的孟德尔随机分析方法和强大的生物信息学工具,旨为阐明UC中风险miRNA及其基因调控网络中的核心枢纽,从而能够捕获复杂的内在关联并增加对UC机制的了解,为未来精准靶向治疗提供新的思路。研究对miRNA暴露数据和全基因组UC遗传关联数据进行了MR分析,鉴定出hsa‑miR‑130b‑3p与UC存在显著的因果关系,并且hsa‑miR‑130b‑3p作为保护因素会减少UC的风险。研究表明,miR‑130b‑3p可在上皮细胞、内皮细胞和巨噬细胞等细胞类型中表达
[16,17],并且参与调控各类细胞生物反应,包括细胞凋亡、增殖、迁移、炎症和肿瘤等生理病理过程
[18,19]。在动物实验中发现,miR‑130b‑3p不仅可以通过靶向介导干扰素调节因子影响巨噬细胞的极化,还能调控脂质代谢、NF‑κB和PPAR等通路,从而改善UC小鼠组织中的炎症反应
[20,21]。此外,近期研究认为miR‑130b‑3p可作为一种新型的内源性抑制剂,通过TLR信号途径改善炎症反应,以用于炎症性疾病的分子靶向治疗
[22]。上述可见hsa‑miR‑130b‑3p与免疫炎症密切相关,并且对于UC疾病性状的具有重要影响。
本课题组为了进一步深入探究UC风险miRNA的下游基因调控机制和作用模式,通过对hsa‑miR‑130b‑3p的靶基因进行预测并与差异的ERS基因交集,再利用机器学习方法,最终鉴定出了5个ERS关键基因,分别为XBP1、PPARG、LRRK2、MAP3K5和EDN1,其中在UC芯片数据集中EDN1和PPARG显著下调,而XBP1、LRRK2和MAP3K5显著上调。同时,结合ROC模型预测分析,考虑联合5个关键基因的模型对于疾病的预测将加更具有敏感性和特异性。随后对5个关键基因进行GSVA富集分析和免疫浸润评估,以探索它们的生物功能通路和免疫调控网络。结果表明,下调的EDN1、PPARG和上调的XBP1、LRRK2、MAP3K5均在趋化因子、肠道免疫炎症调控等通路显著富集。这些分析揭示了UC发病机制的关键生理病理过程,表明了肠道免疫炎症调控网络对于疾病发生发展的重要影响。
XBP1已被广泛熟知,它是内质网应激通路中的关键基因,用以调控内质网跨膜启动信号以及内质网相关降解,并激活促炎信号和凋亡相关通路
[23]。许多UC的研究中表明XBP1可以通过介导肠上皮细胞中的内质网应激反应,从而调节细胞凋亡和炎症反应,并维持肠道内稳态
[4,24]。
LRRK2基因用以编码LRRK2蛋白激酶,它是一种具有激酶和鸟苷三磷酸酶活性的复合酶
[25]。在机体异常条件下,LRRK2错义突变导致LRRK2蛋白激酶过度活跃,上调的LRRK2可促进肠道炎症反应,并且LRRK2蛋白激酶能够诱导泛素化,导致肠道屏障功能损伤,从而增加炎症应激后肠道炎症的易感性
[26]。此外在动物实验中,LRRK2缺陷小鼠通过调节先天免疫反应表现出对实验性结肠炎的高度易感性
[27]。Apoptosis signal‑regulating kinase 1(ASK1,也称为 MAP3K5)作为p38丝裂原活化蛋白激酶和c‑Jun N末端激酶信号级联的关键激活因子,对UC发挥重要调控作用。在实验中ASK1抑制剂不仅作用肠道免疫调控,减轻炎症细胞浸润;还调节肠道紧密连接功能,维持肠道屏障稳定
[28,29]。EDN1是参与炎症过程的一种血管活性肽,已被证明具有调节其下游ERS介质的作用。除了其有效的血管收缩活性外,EDN1还参与细胞增殖、白细胞趋化性和新血管生成调节作用
[30]。PPARG是一类控制生殖、代谢、发育和免疫反应的核受体成员,它参与了多种自身免疫性疾病的发病机制,在调节巨噬细胞的活化和极化、树突状细胞的功能、介导T细胞的增殖和分化以及相关基质细胞方面发挥很大作用
[31]。相关研究报道PPARG在UC中表达下调
[32],另外PPARG激活剂可以介导氧化应激和内质网应激通路抑制细胞死亡
[33]。由此可见,这些关键基因对于UC的诊断、治疗及预后具有重要的潜在价值。
另一方面,通过将关键基因进行免疫细胞浸润分析,结果显示,这些基因与中性粒细胞、巨噬细胞亚群(M0、M2型巨噬细胞等)、和T细胞亚群(滤泡辅助性T细胞、调节性T细胞等)等免疫细胞密切相关。免疫细胞浸润是UC病理生物过程中不可或缺的关键因素,在UC患者的结肠组织和血液中发现大量聚集的活化中性粒细胞,它构成了细胞反应的第一线
[34]。最近研究指出,UC组织整体免疫细胞变化的特征是M0巨噬细胞和中性粒细胞的增加
[35]。在肠道发生炎症时,病原体通过刺激巨噬细胞的过度活跃,导致促炎因子的激活并加剧炎症反应。同时,M2型巨噬细胞的活性受到抑制,进一步阻碍了结肠组织的修复
[36,37]。研究还表明,在UC患者的结肠组织中发现滤泡辅助性T细胞(TFH)明显升高,且与疾病活动性正相关。当TFH细胞的数量或功能发生变化时,会使得生发中心异常反应,从而促进大量抗体和炎症因子的激活,并且调节性T细胞(Treg)的缺陷会进一步诱导过度的T细胞反应和炎症反应,TFH/Treg失衡导致UC肠道炎症病理的损伤
[38,39]。上述所有研究表明,关键基因密切相关的免疫细胞对于UC发病机制的重要意义。
最后,为了再次评估预测结果的可靠性,对5个关键基因进行独立数据集表达验证,结果显示与研究开始预测的结果完全一致。然而,本研究是通过对大数据的挖掘和分析,如若再进行相关体内外实验验证,这将使得本研究更加具有完整性。但是,在分析过程中,利用多个方法如ROC曲线、检验数据集分析等不断验证,已表明这些关键基因对于预测结果的可靠性和有效性。另外,由于MR分析结果与数据来源相关,采用不同来源研究的GWAS进行MR分析所显示UC的风险miRNA结果各不相同。为此,如果扩大后续数据样本量将可能有助于提高预测结果的全面性。
总而言之,本研究通过利用先进的MR分析方法鉴定出了UC中风险hsa‑miR‑130b‑3p。同时利用强大的生物信息学工具,对其靶向内质网应激相关关键基因进行分析,突出强调了免疫调控网络和肠道免疫细胞对于疾病发生发展的重要影响。此外研究还表明联合这5个关键基因模型对于UC疾病的诊断预测具有一定的敏感性和特异性。这些发现为UC的复杂分子机制提供了见解,并为新的诊断及靶向治疗提供了途径。
作者贡献度说明:
陈旭永:阅读文献、数据筛选、数据分析、文章绘图及撰写论文;杨小红:文章芯片、孟德尔随机化数据筛选;王柳丹、吴海东和郭一凡:参与查找、收集文献及提供相关修改意见;苗新普:文章结构设计、论文修改及内容校审。
所有作者声明不存在利益冲突关系。
国家自然科学基金资助项目(82160104)