本方案整合了微生物代谢物靶点预测、直肠黏膜转录组学、蛋白质-蛋白质相互作用与通路富集分析、分子对接、分子动力学模拟以及分子力学/泊松-玻尔兹曼表面积(MM-PBSA)结合自由能估算,以生成一份排序后的候选代谢物相关宿主基因短名单,并对蛋白质-配体复合物按结构优先级排序,为后续实验验证提供假设生成依据。
方法文章
本方案整合了微生物代谢物靶点预测、直肠黏膜转录组学、蛋白质-蛋白质相互作用与通路富集分析、分子对接、分子动力学模拟以及分子力学/泊松-玻尔兹曼表面积(MM-PBSA)结合自由能估算,以生成一份排序后的候选代谢物相关宿主基因短名单,并对蛋白质-配体复合物按结构优先级排序,为后续实验验证提供假设生成依据。
目前尚无标准化的计算流程,用于从公开的化学、基因组和结构数据库中系统性地优先筛选与微生物代谢物相关的宿主基因及蛋白-配体复合物。本文描述了一个八阶段工作流程,该流程接受用户自定义的一组肠道微生物群衍生代谢物,输出一份排序后的候选代谢物相关宿主基因短名单、富集的生物学通路以及结构上优先排序的蛋白-配体复合物,供后续实验验证。该流程整合了以下步骤:(i)化学信息学代谢物分析;(ii)基于多数据库的候选靶点预测,使用蛋白-化学物相互作用数据库、基于配体的靶点预测工具及分子对接程序;(iii)对公开转录组数据进行差异表达基因分析;(iv)预测靶基因与差异表达基因的交集;(v)蛋白质-蛋白质相互作用网络构建与通路富集分析;(vi)分子对接分析,使用分子对接程序;(vii)采用分子动力学引擎结合用于分子动力学模拟的蛋白力场进行200 ns分子动力学模拟;(viii)MM-PBSA结合自由能估算。以实例演示,本研究处理了九种代表短链脂肪酸、胆汁酸、色氨酸衍生代谢物及尿石素A的肠道微生物群衍生或微生物群修饰代谢物,使用公开的肠易激综合征便秘型(IBS-C)直肠黏膜转录组数据集GSE36701。该工作流程在该数据集中筛选出17个独特且差异表达的预测代谢物相关基因。通过分子对接、分子动力学模拟和MM-PBSA分析,结构上优先排序出五个代谢物-蛋白复合物:lithocholic acid-VDR、lithocholic acid-NR1H4/FXR、ursodeoxycholic acid-NR1H4/FXR、tryptamine-HTR2A(在显式1-Palmitoyl-2-oleoyl-sn-甘油-3-磷脂酰胆碱(POPC)脂双层)以及尿石素A-CASP3。该方案设计为可适用于其他代谢物组合、疾病转录组数据集及靶点类别;所有输出结果均为假设生成型的计算预测,需经过独立的转录组学重复验证、蛋白水平验证以及功能性配体-反应实验,方可得出因果关系或治疗意义的结论。
伴便秘的肠易激综合征(IBS-C)是一种常见的功能性胃肠道疾病,其特征为反复发作的腹痛、排便习惯改变、腹胀和便秘,全球患病率估计约为总人口的10–15%1,2。目前的药物治疗包括促分泌剂、促动力药和解痉药,可在部分患者中改善个别症状;然而,治疗反应存在异质性,且持久缓解较为罕见,这反映了该疾病复杂的多因素致病机制1,3,4。因此,需要更全面地了解肠道微生物信号在黏膜水平的转导机制,以提出可验证的新治疗靶点假说。
肠道微生物群通过产生和生物转化多种化学结构不同的代谢物,如短链脂肪酸(SCFAs)、次级胆汁酸、色氨酸衍生物以及鞣花素等多酚类衍生物,参与维持下胃肠道的稳态5,6,7,8。这些分子通过一系列广泛但尚未完全明确的分子靶点与宿主细胞进行通讯,其作用范围远超经典的代谢物感应膜受体,还包括核受体、胞质酶、组蛋白修饰蛋白、肽类激素前体以及细胞内信号蛋白9。已有研究记录到肠易激综合征(IBS)患者肠道微生物群落组成及代谢物谱的改变,这为探讨与微生物代谢物响应相关的宿主基因在IBS-C患者直肠黏膜中是否存在转录水平的扰动提供了生物学依据10。
该九种代谢物组合被预先定义,旨在提供一组紧凑、化学多样性丰富且生物学意义明确的由肠道微生物群产生或经微生物群修饰的小分子。筛选基于五个标准:代表参与宿主-微生物群信号传导的主要微生物代谢物类别;具有已知或合理的远端肠道黏膜暴露性;具备明确的 PubChem 标识符和标准化学结构;分子大小和结构适合基于配体的靶点预测与分子对接;以及在便秘型肠易激综合征(IBS-C)中对上皮、神经免疫、肠内分泌、核受体或运动功能相关信号通路具有先前合理性。所选组合包括作为短链脂肪酸(SCFAs)的丁酸和丙酸;作为胆汁酸的鹅脱氧胆酸、石胆酸和熊脱氧胆酸;作为色氨酸衍生代谢物的色胺、吲哚-3-丙酸和吲哚-3-乳酸;以及作为肠道微生物群衍生多酚类代谢物的尿石素 A5,6,7,8,9,10。
以往大多数计算和实验研究都是孤立地考察单个代谢物-受体或代谢物-酶的相互作用,这种方法无法捕捉微生物代谢物信号在宿主通路中广泛分布且具有汇聚性的本质9,11。多个分析阶段的整合提供了相互增强的筛选能力,这是单一阶段无法独立实现的。针对经过审编的数据库进行计算靶点预测,可为每种代谢物获得一组广泛的潜在宿主蛋白候选。与疾病相关转录组数据的交集分析可显著缩小该候选集,仅保留那些在疾病背景下转录水平发生变化的候选分子。随后,通路富集分析和蛋白质-蛋白质相互作用网络分析将缩小后的候选列表映射到已知的生物学功能模块。分子对接提供了对每个候选复合物结合口袋互补性的初步计算评估,而进一步的200纳秒分子动力学(MD)模拟结合MM-PBSA结合自由能分解,则为结构优先级排序提供了时间分辨的热力学维度信息,这是仅靠对接评分无法获得的。若各步骤独立进行,缺乏系统性整合与逐级筛选,则所得候选列表将过于宽泛,难以进行实验验证,也无法识别出通路间的汇聚性结构特征。
在本实验方案的框架内,“代谢物相关基因”(MAG)指人类基因,其编码的蛋白质产物已被至少一个经过 curated 的计算预测数据库提名作为一种或多种肠道微生物群衍生代谢物的潜在分子靶点,并且该基因的转录本在用于演示本工作流程的疾病相关转录组数据集中表现出差异表达。这一操作性定义明确包括膜受体、核受体、胞质酶、信号转导蛋白、肽类激素前体以及其他细胞内蛋白。MAG 的指定并不代表代谢物与靶蛋白结合、形成蛋白-配体复合物、激活受体、改变蛋白丰度或导致疾病的实验证据,而是通过计算方法得出的、用于生成假设的提名,需进一步实验验证。
本方案描述了一个完整的八阶段计算工作流程(图1),提供了足够的操作细节,以支持独立重复、适应其他代谢物面板或疾病数据集,并扩展至其他宿主-微生物群相互作用的研究场景。该工作流程明确界定为一种假设生成和结构优先排序的框架,仅基于公开可用的组学和结构资源进行操作,不单独从计算结果中推断代谢物浓度变化、受体激活状态、蛋白质表达改变、下游信号活性或临床意义。本文以九种由肠道微生物群产生或经微生物群修饰的代谢物,以及公开的肠易激综合征便秘型(IBS-C)直肠黏膜转录组数据集GSE36701为例,演示该方案的具体应用,旨在识别代谢物-基因关联(MAGs),并对代谢物-蛋白质复合物进行优先排序,以供后续实验验证。
分析仅使用了来自GSE36701的公开可用的去标识转录组数据,以及公开可获取的化学、蛋白质和结构数据库。这些数据库的访问时间介于2026年1月至2026年5月之间。任何更晚的访问日期均记录在单独的材料表中。
1. 研究设计、硬件和软件要求
2. 代谢物筛选与化学信息学表征
3. 候选人类靶点预测
4. 转录组数据集与差异基因表达分析
5. 靶基因与差异表达基因重叠分析及统计评估
6. 蛋白质-蛋白质相互作用网络分析与通路富集
7. 分子对接
8. 分子动力学模拟
9. MM-PBSA 结合自由能估算
候选代谢物相关靶点
这九种代谢物在化学-蛋白质相互作用靶点预测和分子对接程序中产生了异质性的预测靶点集合。丙酸、色胺、胆汁酸和尿石素A均得到多个与胃肠道信号传导相关的已知靶点。预测的靶点谱包括经典的膜受体、核受体、细胞内酶、信号转导蛋白以及与肽类激素相关的蛋白。因此,后续结果被描述为代谢物相关基因(MAGs),而不仅限于受体发现(Table 1)。
与已报道的代谢物-蛋白质相互作用进行基准比较
为了将靶点预测结果与现有的实验知识进行基准比较,将预测的代谢物相关靶点关系分为三个证据等级:(i)实验支持的直接或近类别的代谢物-蛋白质相互作用,即已有文献报道该代谢物或结构相近的内源性代谢物能够结合、激活、抑制或功能调控编码蛋白;(ii)通路或靶点类别支持的相互作用,即预测的靶点属于已知的代谢物响应通路或受体家族,但针对该特定代谢物-蛋白质对的直接实验证据有限;(iii)仅基于计算关联,即在所查阅文献中未发现该对代谢物与蛋白质之间存在直接实验相互作用的证据。该基准比较旨在为预测的代谢物-靶点关联(MAGs)提供背景信息,而非用于验证其正确性。
多个预测结果重现了先前报道的生物学现象。丙酸-FFAR2被视为具有实验支持,因为FFAR2/GPR43是经典的短链脂肪酸受体。丁酸-HDAC3被归类为具有实验或类别支持,因为丁酸是公认的组蛋白去乙酰化酶抑制剂,且预测的相互作用涉及HDAC家族成员。涉及NR1H4/FXR和VDR的胆汁酸相关预测被认为得到了已确立的胆汁酸核受体生物学的支持,特别是对于疏水性胆汁酸(如LCA);而熊去氧胆酸(UDCA)相关的FXR预测则需谨慎解读,因为UDCA通常是一种较弱或依赖于情境的FXR配体。色胺相关的HTR1B、HTR2A、HTR2B和HTR6预测被归类为由5-羟色胺通路支持,而非已确认的直接受体特异性相互作用,因为色胺是一种由微生物色氨酸衍生的单胺类物质,而5-羟色胺受体已被证实为调控胃肠道运动和分泌的关键因子。尿石素A-CASP3被认为具有通路支持,依据是已有文献报道尿石素A与凋亡/半胱天冬酶相关反应之间的关联,但尚无CASP3直接结合的证据。吲哚-3-乳酸-KYAT1和吲哚-3-丙酸-KYAT1则被保留为仅基于计算的假设,因为现有广泛文献支持微生物来源的吲哚衍生物在宿主信号传导中的作用,但尚无这两种特定代谢物与KYAT1直接结合的证据7,8,38,39,40。
因此,表1区分了计算靶点提名与先前实验或通路支持的水平。对于每个靶点,还提供了预测来源(化学物-蛋白质相互作用靶点预测、分子对接程序,或两者均有)、化学物-蛋白质相互作用靶点预测的综合相互作用评分,以及通过分子对接程序鉴定靶点时的对接概率。那些缺乏直接先前实验证据的预测靶点被描述为候选的代谢物相关基因,需进行独立的蛋白质水平和配体响应验证。
预测靶基因与便秘型肠易激综合征(IBS-C)差异表达基因的重叠
联合预测靶基因列表与基因水平差异表达结果的交集分析,鉴定出17个独特的预测代谢物相关基因,这些基因在便秘型肠易激综合征(IBS-C)患者与健康志愿者的比较中表现出显著的差异表达。所有17个基因均呈下调表达。该基因集合包括膜受体和核受体(CASR、FFAR2、GPR68、HTR1B、HTR2A、HTR2B、HTR6、NR1H4、TBXA2R、VDR)以及非受体蛋白(CASP3、GCG、GNAQ、GPHN、HDAC3、KYAT1、MLN)(表1,图2A,B)。
全部17个MAGs的错误发现率(FDR)均低于0.05;其中16个达到更严格的FDR < 0.001标准,仅余下的一个基因(HTR1B)在FDR < 0.05水平上显著。在17个靶基因中,有7个(CASP3、GCG、GNAQ、GPHN、GPR68、HDAC3、TBXA2R)同时满足FDR < 0.001且绝对log2倍数变化超过1.0(logFC范围为−1.34至−1.10),表明该亚组基因表现出强烈且一致的下调。其余靶基因则表现为中等程度但具有统计学显著性的下调(|logFC|范围为0.45至0.97)。鉴于该数据集具有全基因组表达特征(见下文统计评估),对该一致描述性模式的解读需持谨慎态度。
靶基因-差异表达基因重叠的统计评估
为了正式评估17个基因重叠的统计学显著性,采用单尾Fisher精确检验,以17个预测靶基因作为查询集,以GSE36701中检测到的全部18,296个唯一基因折叠条目作为基因组背景。在此背景中,有17,296个基因(94.5%)在FDR < 0.05水平上差异表达,反映出在IBS-C直肠黏膜比较中近乎普遍的转录抑制现象。所有17个预测靶基因均属于差异表达基因(实际重叠为17/17,占100%)。鉴于背景中差异表达基因为94.5%,任意随机选择的17个基因集合的预期重叠数为16.1个基因。Fisher精确检验得到 p = 0.384,连续性校正后的优势比为2.03(95%置信区间0.12–33.73),在α = 0.05水平上无统计学显著性(图3A–C)。
该结果表明,观察到的17/17重叠并未超过在此数据集全基因组表达谱下随机预期的重叠程度。因此,这些发现应被解读为一种描述性的方向性模式,即所有17个预测靶基因在便秘型肠易激综合征(IBS-C)直肠黏膜组织中均一致且显著地下调,而非在基因组背景下具有统计学富集或独立验证的证据。正式的富集检验需要在差异表达谱更具选择性的转录组数据集中进行重复验证,其中显著差异表达的基因应远少于全部基因的一半。需要强调的是,这17个重叠基因的统一性下调是一种描述性观察结果,而非经过独立验证的统计学结果,因为该数据集的差异表达背景本身以整体下调为主,重叠基因之间存在共同的下调方向是预期之中的现象,且未经过正式的方向性检验。因此,这种一致的方向性不应被解释为协调性、代谢物特异性调控的独立统计学证据。
特异性代谢物模式
丙酸在重叠基因数量上最多,包括 CASR、FFAR2、GCG、GNAQ、GPHN、GPR68、MLN 和 TBXA2R,提示短链脂肪酸响应性及与 Gq 相关的信号通路可能参与其中。丁酸与 HDAC3 存在重叠,这与丁酸相关的组蛋白去乙酰化酶生物学功能一致,但仅凭 mRNA 下调尚不足以确立丁酸响应性的改变。胆汁酸相关重叠基因包括核受体 VDR 和 NR1H4,二者均为肠道中胆汁酸信号传导的已知效应分子38,39。色胺与 HTR1B、HTR2A、HTR2B 和 HTR6 存在重叠,提示5-羟色胺能信号通路可能作为一个候选模块,该系统在胃肠道运动和分泌中具有明确作用40。吲哚-3-乳酸和吲哚-3-丙酸与 KYAT1 存在重叠,尿石素 A 与 CASP3 存在重叠。
通路富集
对17个重叠基因的功能富集分析揭示了与G蛋白偶联受体(GPCR)下游信号传导、Gαq信号传导、GPCR配体结合、5-羟色胺能突触、神经活性配体-受体相互作用、钙信号转导、cAMP信号通路以及肽类激素分泌相关的通路。这些结果与基因集的组成一致,支持其生物学上的协调性,但反映的是所提交基因的功能注释,而非通路水平活性的独立证据。
蛋白质-蛋白质相互作用网络结构
蛋白质-蛋白质相互作用网络的构建和通路富集分析在三个互补网络中进行了解读。在整合的17基因元网络(网络1)中,最显著的、有注释支持的结构是以GNAQ为核心的GPCR/Gαq信号通路组分,将GNAQ与受体相关基因(包括TBXA2R、CASR、HTR2A和HTR2B)连接起来。5-羟色胺受体之间的连接较为有限,其中HTR2A与HTR2B之间的连接最为显著,而在所选置信度阈值下,其他多个基因则保持孤立或连接较弱。丙酸特异性网络(网络2)表现出更受限的拓扑结构,其中GNAQ保留了与CASR和TBXA2R之间有注释支持的连接,而FFAR2、GPR68、GCG、GPHN和MLN则处于孤立或弱连接状态。色胺/5-羟色胺网络(网络3)包含HTR1B、HTR2A、HTR2B和HTR6;在此子集中,HTR2A与HTR2B显示出主要的、有注释支持的连接,而HTR1B和HTR6在所选阈值下未直接相连(图4A–C)。
分子对接
对五种选定的代谢物-蛋白质复合物进行了分子对接。胆汁酸-核受体组合的Vina评分优于尿石素A-CASP3和色胺-HTR2A。LCA-VDR的评分最佳,为−10.0 kcal/mol,其次是LCA-NR1H4/FXR(−9.9 kcal/mol)和UDCA-NR1H4/FXR(−9.4 kcal/mol)。尿石素A-CASP3和色胺-HTR2A的评分较低,但仍合理,分别为−7.1 kcal/mol(表2)。
对于LCA-VDR复合物(PDB ID: 1DB1),预测的结合构象得到了一个经典氢键的支持,即LCA羧酸根氧原子与Ser278之间的氢键(4.29 Å),同时还存在广泛的疏水相互作用,涉及Leu230、Val234、Trp286、Val300、His305、Tyr295、Leu233和His397;此外还与Met272、Leu313、Ile271、Ile268、Leu309、Phe422、Val418、Ala231、Ala303、Cys288、Ser275和Phe150存在额外的范德华接触。排名最高的构象具有−10.0 kcal/mol的Vina评分、2055 Å3的空腔体积,以及坐标为(10, 19, 33)的网格中心(表3,图5A,B)。
对于LCA-NR1H4/FXR复合物(PDB ID: 3DCT),其对接得分为−9.9 kcal/mol,预测存在涉及His294和Ile335的氢键,与His294的π-σ相互作用,以及涉及Met290、Met328、Ala291、Leu287、Ile352和His447的疏水性烷基或π-烷基接触,并有进一步的范德华相互作用,支持类固醇骨架在FXR结合口袋中的容纳(表4,图6A,B)。
UDCA-NR1H4/FXR 复合物(PDB ID: 3DCT)的预测构象显示,其与 His447 形成一个常规氢键(3.66 Å),与 Gly322 形成另一个氢键(3.46 Å),与 Val325 存在 π-阴离子相互作用(4.96 Å),并与 Trp469 形成碳-氢键(4.51 Å)。相互作用图谱还发现与 Arg395(3.89 Å)和 Gln396(3.40 Å)存在不利的供体-供体接触,表明与同一受体结合口袋中 LCA 相比,UDCA 较低的 Vina 评分可能是由于局部几何结构或静电环境的不利所致(表 5,图 7A,B)。
在尿石素A-CASP3复合物(PDB ID: 2DKO)中,预测的结合模式显示与Gln161(3.78和4.19 Å)、Ser120(3.95 Å)以及Arg207(3.05和3.77 Å)形成常规氢键,并通过与Arg207的π-阳离子相互作用、与Cys163的π-供体氢键,以及涉及Arg64、Ala162、His121、Ser205和Trp206的额外π-烷基和范德华接触进一步稳定(表6,图8A,B)。
对于色胺-HTR2A复合物(PDB编号:6A93),预测的结合构象通过色胺的质子化胺基与Asp155之间形成的静电盐桥而稳定,Asp155是跨膜螺旋3中保守的天冬氨酸残基(Ballesteros-Weinstein编号中的D3.32),该残基在5-羟色胺及其相关受体中可锚定胺类配体的质子化胺基41,42,43;此外,还存在与Thr160和Ser159的氢键作用、与Phe340和Trp336的芳香环相互作用,以及与Val156和Ile163的π-烷基相互作用。与Tyr370、Phe339、Ser242、Phe243、Phe332和Leu123的额外范德华接触进一步支持了正构结合位点的结合模式(表7,图9A、B)。
对接方案验证
为评估对接方案的可靠性,进行了两项互补的对照实验。在自对接(阳性)对照中,将共结晶配体从其参考X射线结构中提取,并重新对接至其天然结合位点。维生素D类似物VDX在VDR/1DB1中的最高排序预测构象与其晶体结构位置的偏差为0.87 Å,而FXR/3DCT中共结晶配体WAY-362450的偏差为1.79 Å;这两个值均低于常规接受阈值2.0 Å,表明该对接方案在这些受体系统中具有几何学上的有效性(图10A,B)。在交叉对接(阴性)对照中,将石胆酸对接至半胱天冬酶-3(2DKO)——一种其并非已知配体的半胱氨酸蛋白酶——所得预测得分(−8.3 kcal/mol)比其在对应靶标VDR上的得分(−10.0 kcal/mol)低1.7 kcal/mol,符合预期的结合位点选择性。色胺对接至VDR的预测得分为−6.4 kcal/mol,而其在对应靶标HTR2A上的得分为−7.1 kcal/mol,差异为0.7 kcal/mol,该差异处于代谢物配体与靶蛋白分子对接得分的报道不确定性范围内,因此表明该小分子配体的预测选择性较弱(图10C)。综上所述,这些对照实验表明,在所测试条件下,该对接方案能够重现已知的结合几何构型,并能区分对应与非对应配对关系,但其计算预测结果仍不能替代实验测得的亲和力数据(表8)。
分子动力学模拟
对五种优先选择的复合物进行了200 ns生产轨迹的分子动力学模拟。四个可溶性核受体复合物在显式水溶剂中进行模拟,而色胺-HTR2A复合物则在显式POPC脂双层中进行模拟,以为该G蛋白偶联受体提供生理上合适的膜环境。分析评估了对接构象在时间依赖条件下的动态稳定性,并允许比较不同复合物之间的相对结构行为(表9)。
LCA-VDR/1DB1 复合物的 RMSD 曲线显示,在前 10 ns 内经历了一个短暂的平衡期,随后进入稳定平台期,波动范围主要在 0.20–0.28 nm 之间图11A)。RMSF 值较低,主链波动也较小 < 大多数残基为0.15 nm(图11B). 氢键分析显示存在一个持续的2–5个氢键的网络,偶尔可增加至7个(图11C). 回转半径(Rg)保持在1.25–1.75 nm范围内,溶剂可及表面积(SASA)维持在约130 nm2 (图11D, E).
urolithin A-CASP3/2DKO 复合物表现出更强的动态活性。RMSD 初始上升,随后在 0.4 至 0.7 nm 之间波动,并在约 165 ns 时出现短暂的高偏差事件(图 12A)。RMSF 分析显示残基水平具有较高的移动性,其中残基 175 附近的柔性环区域波动最大(图 12B)。氢键分析表明,在前 30–40 ns 存在一个包含约 2–5 个氢键的广泛网络,随后主要维持在 0 至 2 个间歇性氢键(图 12C)。相应的回转半径和溶剂可及表面积(SASA)曲线如 图 12D、E 所示。
对于NR1H4/FXR(3DCT)胆汁酸体系,主链RMSD值在大部分轨迹过程中保持在相对狭窄的范围内(图13A),而RMSF值显示核心区域的移动性较低,柔性区域的波动较高(图13B)。LCA-3DCT复合物在整个轨迹过程中维持了大约3至4个持续存在的氢键,而UDCA-3DCT复合物则表现出更大的氢键波动性,并在约125 ns后氢键数量减少。LCA和UDCA结合体系的回转半径曲线分别如图13C、D所示,相应的溶剂可及表面积(SASA)曲线如图13E、F所示。
色胺-HTR2A复合物的膜分子动力学
在包含258个脂质分子的显式POPC脂双层、显式三站点水模型以及0.15 M NaCl环境中,对色胺-HTR2A/6A93复合物进行了200 ns的模拟,总体系规模约为100,925个原子33,44,45。在整个轨迹过程中,受体稳定地嵌入脂双层中(图14)。主链RMSD在前100 ns内从约0.10 nm上升至接近0.15–0.20 nm的稳定平台期,并在此后保持稳定,所有值均低于0.25 nm,表明受体在膜环境中维持了稳定的构象,未发生全局去折叠(图15A)。残基分辨率的RMSF显示跨膜螺旋核心区域波动较小,而环区和末端区域表现出预期的较高灵活性,符合典型GPCR的柔性特征(图15B)。回转半径紧密限制在约2.06至2.12 nm之间,SASA在狭窄范围内波动且无持续漂移趋势,两者均证实了跨膜螺旋束的紧凑结构得以保持(图15C,D)。
在整个轨迹过程中,蛋白质与配体之间的氢键得以维持(图15E),氢键数量出现显著波动,范围为1至3个。为了特异性评估关键离子相互作用的持续性,监测了色胺质子化的铵基氮原子与Asp155(D3.32)羧基氧原子之间的最小距离在整个轨迹过程中的变化。该距离始终紧密分布在平均0.270 nm处(最小值0.247 nm,最大值0.424 nm),盐桥接触(< 0.4 nm)在99.9%的模拟时间内得以维持,仅出现两次短暂的瞬时偏离,未发生持续性解离事件(图16)。这些结果表明,在整个膜模拟过程中,保守的Asp155离子相互作用足以稳定色胺在HTR2A正构结合口袋中的构象。
MM-PBSA 结合自由能及残基分解
采用MM-PBSA分析为五个复合物(表10)增加了额外的能量优先级评估层。对于四个水相复合物,按残基分解分析确定了每种预测结合模式中的主要能量贡献残基。在LCA-VDR/1DB1复合物中,配体和Gln317表现出有利的贡献,而Trp286则表现出不利的贡献。在尿石素A-CASP3/2DKO复合物中,Arg64和Arg207显示出强烈的负向残基贡献,表明存在显著的极性或静电稳定作用;然而,相应的轨迹仍表现出高度动态性,说明仅有利的残基水平能量特性并不足以确保复合物的持续稳定性。对于3DCT体系,LCA的结合主要由Arg331驱动,而UDCA的结合则涉及一个更广泛分布的能量网络,包括Glu326、Asp394、Arg395、Arg441和Asp470。在四个水相体系中,MM-PBSA分解结果支持对基于LCA的复合物进行相对优先级排序。
对于膜嵌入的色胺-HTR2A/6A93复合物,对从双层膜轨迹中提取的蛋白质-配体子系统进行了MM-PBSA分析46,47。配体和Asp155(D3.32)表现出有利的贡献,其中Asp155在残基水平上的稳定作用最为显著,这与对接分析和轨迹距离分析中发现的盐桥相互作用一致。在周围正构结合口袋的残基(Ser86、Phe87、Phe133、Phe140、Phe141、Val156、Ser159、Thr160、Ile163、Val167、Tyr171)中,Trp137表现出最大的不利残基贡献,这些残基共同构成了结合口袋内衬的芳香性和极性接触网络。这些数值代表了用于结构优先排序的相对计算估计值,并非实验测得的结合亲和力。

图 1:针对便秘型肠易激综合征(IBS-C)中与代谢物相关的宿主基因优先排序的计算工作流程。 该八阶段工作流程的示意图,整合了代谢物筛选、靶点预测、转录组差异表达分析、重叠分析、网络与通路富集分析、分子对接、分子动力学模拟以及MM-PBSA结合自由能分析。请点击此处查看此图的放大版本。

图2:IBS-C 黏膜中差异表达与代谢物靶基因重叠分析。(A)GSE36701 数据集中基因水平差异表达的火山图。蓝色点表示显著下调的基因;红色点表示显著上调的基因;灰色点表示无显著性差异的基因。标注了部分重叠的与代谢物相关的基因。(B)预测的330个独特代谢物靶基因与GSE36701中下调基因之间的重叠维恩图;共有17个基因重叠。请点击此处查看该图的放大版本。

图3:针对GSE36701数据集对17个预测的代谢物靶基因进行统计学评估。(A)17个基因中每个基因的log2倍数变化,按显著性层级着色。(B)背景基因与预测靶基因的差异表达率,采用Fisher精确检验。(C)Fisher精确检验所用的2×2列联表。所有17个靶基因均显著下调;该重叠结果被解释为一种描述性的方向性模式,而非统计学富集。请点击此处查看该图的放大版本。

图4:重叠代谢物相关基因的复合蛋白质-蛋白质相互作用网络构建及通路富集蛋白质-蛋白质相互作用网络。(A)网络1:全部17个基因的整合元网络。(B)网络2:8个基因的丙酸特异性网络(CASR, FFAR2, GCG, GNAQ, GPHN, GPR68, MLN, TBXA2R)。(C)网络3:4个基因的色胺/血清素网络(HTR1B, HTR2A, HTR2B, HTR6)。所有网络均基于智人(Homo sapiens)构建,采用蛋白质-蛋白质相互作用网络构建及通路富集置信度≥0.700的条件。边表示有注释支持的功能关联。请点击此处查看该图的放大版本。

图5:胆石酸与维生素D受体(VDR)复合物的三维和二维结构示意图(PDB ID: 1DB1)。(A)三维表面与卡通图示,其中胆石酸以球状表示。(B)二维相互作用图,显示Ser278的氢键及周围的疏水和范德华接触。请点击此处查看该图的放大版本。

图6。石胆酸与NR1H4/FXR复合物的三维和二维结构示意图(PDB ID: 3DCT)。(A)三维表面及卡通模型示意图。(B)二维相互作用图,显示与His294和Ile335之间的氢键、π-σ相互作用以及周围的接触关系。请点击此处查看该图的放大版本。

图7:熊去氧胆酸与NR1H4/FXR复合物的三维和二维结构示意图(PDB编号:3DCT)。(A)三维表面与卡通模型示意图。(B)二维相互作用图,显示与His447和Gly322形成氢键,与Val325形成π-阴离子相互作用,与Trp469形成碳-氢键,以及与Arg395和Gln396存在不利的供体-供体接触。请点击此处查看该图的放大版本。

图8:尿石素A与CASP3复合物的三维和二维结构示意图(PDB编号:2DKO)。(A)三维表面与卡通模型示意图。(B)二维相互作用图,显示与Gln161、Ser120和Arg207之间的氢键作用,与Arg207的π-阳离子相互作用,与Cys163的π-供体氢键,以及周围的接触作用。请点击此处查看该图的放大版本。

图9:色胺与HTR2A复合物的三维和二维结构示意图(PDB ID: 6A93)。(A)在三维分子可视化程序中生成的三维表面和卡通模型示意图。(B)利用分子可视化软件及二维相互作用图工具生成的二维相互作用图,展示了Asp155盐桥及其他结合位点的相互作用。请点击此处查看该图的放大版本。

图10:对接方案验证。(A,B)共结晶配体重新对接至VDR/1DB1(RMSD 0.87 Å)和FXR/3DCT(RMSD 1.79 Å);晶体结构构象与重新对接构象重叠,两者均低于2.0 Å的可接受阈值。(C)交叉对接选择性:石胆酸和色胺的同源与非同源Vina评分。请点击此处查看该图的放大版本。

图11。LCA-VDR/1DB1复合物在200 ns内的分子动力学轨迹分析。(A)RMSD变化曲线。(B)RMSF变化曲线。(C)氢键数量。(D)回转半径变化曲线。(E)溶剂可及表面积(SASA)变化曲线。请点击此处查看该图的放大版本。

图12:尿石素A-CASP3/2DKO复合物在200 ns内的分子动力学轨迹分析。(A)RMSD曲线,显示广泛的构象波动,并在约165 ns处出现短暂的高偏差事件。(B)RMSF曲线,显示在残基175附近存在显著的残基水平灵活性。(C)氢键数量。(D)回转半径曲线。(E)溶剂可及表面积(SASA)曲线。请点击此处查看该图的放大版本。

图13:NR1H4/FXR(3DCT)胆汁酸体系在200 ns内的分子动力学轨迹分析。(A)3DCT复合物的主链RMSD曲线。(B)主链RMSF曲线。(C)3DCT-LCA的回转半径曲线。(D)3DCT-UDCA的回转半径曲线。(E)3DCT-LCA的溶剂可及表面积(SASA)曲线。(F)3DCT-UDCA的溶剂可及表面积(SASA)曲线。请点击此处查看该图的放大版本。

图14:嵌入显式POPC脂双层中的色胺-HTR2A复合物。 受体以跨越脂双层的卡通形式展示,POPC脂质以线条表示,其中磷酸头部基团被突出显示,色胺位于正构结合口袋内。膜的上下方显示有水分子。请点击此处查看该图的放大版本。

图15:色胺-HTR2A/6A93复合物在显式POPC脂双层中200 ns的分子动力学轨迹分析。(A)主链RMSD曲线。(B)残基级RMSF曲线。(C)回转半径曲线。(D)溶剂可及表面积(SASA)曲线。(E)蛋白-配体氢键数量。请点击此处查看该图的放大版本。

图16:色胺与Asp155(D3.32)离子相互作用在200 ns膜模拟轨迹中的持续性。 色胺铵氮原子与Asp155羧酸根氧原子之间的最小距离随时间变化情况如图所示;虚线表示0.4 nm的盐桥接触阈值。在整个模拟过程中,该接触保持了99.9%的时间。请点击此处查看该图的放大版本。
| 基因符号 | 来源代谢物 | 功能类别 | log2FC | FDR(校正P值) | 显著性等级 |
| GCG | 丙酸 | 肽类激素相关蛋白 | −1.342 | 1.97e−7 | FDR <0.001 & |logFC > 1 |
| HDAC3 | 丁酸 | 酶 | −1.234 | 2.44e−6 | FDR <0.001 & |logFC| > 1 |
| CASP3 | 尿石素A | 酶 | −1.198 | 6.66e−7 | FDR <0.001 & |logFC| > 1 |
| GPR68 | 丙酸 | 膜受体 | −1.137 | 4.35e−6 | FDR <0.001 & |logFC| > 1 |
| GNAQ | 丙酸 | 细胞内信号蛋白 | −1.122 | 1.05e−6 | FDR <0.001 & |logFC| > 1 |
| GPHN | 丙酸 | 其他细胞内蛋白 | −1.109 | 1.13e−6 | FDR <0.001 & |logFC| > 1 |
| TBXA2R | 丙酸 | 膜受体 | −1.104 | 4.04e−7 | FDR <0.001 & |logFC| > 1 |
| HTR6 | 色胺 | 膜受体 | −0.967 | 2.17e−5 | FDR <0.001 |
| VDR | 石胆酸 | 核受体 | −0.942 | 5.73e−7 | FDR <0.001 |
| HTR2A | 色胺 | 膜受体 | −0.937 | 4.99e−6 | FDR <0.001 |
| FFAR2 | 丙酸 | 膜受体 | −0.889 | 1.44e−4 | FDR <0.001 |
| NR1H4 | 石胆酸 / 熊去氧胆酸 | 核受体 | −0.861 | 3.68e−6 | FDR <0.001 |
| HTR2B | 色胺 | 膜受体 | −0.702 | 1.29e−4 | FDR <0.001 |
| MLN | 丙酸 | 肽类激素相关蛋白 | −0.605 | 7.39e−5 | FDR <0.001 |
| KYAT1 | 吲哚-3-乳酸 / 吲哚-3-丙酸 | 酶 | −0.530 | 3.61e−4 | FDR <0.001 |
| CASR | 丙酸 | 膜受体 | −0.483 | 4.05e−4 | FDR <0.001 |
| HTR1B | 色胺 | 膜受体 | −0.455 | 3.18e−2 | FDR <0.05 |
表1:IBS-C直肠黏膜数据集中差异表达基因与预测的代谢物相关靶基因的重叠情况。 所列重叠基因均为下调基因。表1作为独立的电子表格提交,其中列出了每个靶基因对应的来源代谢物、功能类别、靶点预测来源(化学-蛋白质相互作用靶点预测、分子对接程序,或两者均有)、化学-蛋白质相互作用靶点预测的综合相互作用评分,以及在可用情况下的分子对接程序预测概率、预测等级、log2倍数变化、FDR值及表达显著性等级。数据来源:基因表达值来自基因合并后的GSE36701差异表达表(每个基因取FDR最低的探针)。靶点预测来源及置信度值根据化学-蛋白质相互作用靶点预测和分子对接程序的输出结果整合而成,筛选阈值为化学-蛋白质相互作用靶点预测综合相互作用评分≥0.700,分子对接程序预测概率≥0.70。化学-蛋白质相互作用靶点预测评分为0–1范围内的综合评分;STP表示分子对接程序预测概率。等级1(Tier 1)= 化学-蛋白质相互作用靶点预测的严格支持;等级1+(Tier 1+)= 化学-蛋白质相互作用靶点预测的严格支持并得到分子对接程序的交叉验证。
| 复合物 | 蛋白质 (PDB ID) | 配体 | Vina 评分 (kcal/mol) | 空腔体积 (A^3) | 网格中心 X,Y,Z (A) | 搜索框 (A) |
| LCA-VDR | VDR (1DB1) | 石胆酸 | −10.0 | 2055 | 10, 19, 33 | 25 x 25 x 25 |
| LCA-NR1H4/FXR | NR1H4/FXR (3DCT) | 石胆酸 | −9.9 | 3395 | 137, 31, 78 | 25 x 25 x 25 |
| UDCA-NR1H4/FXR | NR1H4/FXR (3DCT) | 熊去氧胆酸 | −9.4 | 3395 | 137, 31, 78 | 25 x 25 x 25 |
| 尿石素 A-CASP3 | CASP3 (2DKO) | 尿石素 A | −7.1 | 233 | 37, 34, 32 | 25 x 25 x 25 |
| 色胺-HTR2A | HTR2A (6A93) | 色胺 | −7.1 | 3238 | 12, −1, 61 | 25 x 25 x 25 |
表2:分子对接结果:五种优先级最高的蛋白-配体复合物中,代谢物配体与靶标蛋白对接的最高排名得分及结合腔参数。 结合腔大小以 Å3 表示。数据来源:Docking_Validation/Results/Docking_Validation_Results.xlsx,工作表“Original_Docking_Scores”。对代谢物配体与靶标蛋白进行分子对接;对接搜索强度(exhaustiveness)= 8,随机种子(seed)= 42(固定),所有复合物的生成构象数(num_modes)= 9;报告最高排名(模式1)的对接构象。
| 相互作用类型 | 残基 | 距离 (Å) | 备注 |
| 传统氢键 | Ser278 | 4.29 | LCA 羧基氧 |
| 疏水/π-烷基接触 | Leu230, Val234, Trp286, Val300, His305, Tyr295, Leu233, His397 | - | |
| 范德华接触 | Met272, Leu313, Ile271, Ile268, Leu309, Phe422, Val418, Ala231, Ala303, Cys288, Ser275, Phe150 | - |
表3:鹅脱氧胆酸与VDR(PDB ID: 1DB1)对接生成的结合模式。 来源:分子可视化及二维相互作用图工具“2D配体-残基相互作用图”,如论文结果部分(分子对接)所述。“-”表示该相互作用接触的距离值未单独列出。
| 相互作用类型 | 残基 | 距离 (Å) | 备注 |
| 氢键 | His294 | - | |
| 氢键 | Ile335 | - | |
| Pi-σ 相互作用 | His294 | - | |
| 烷基 / Pi-烷基(疏水相互作用) | Met290, Met328, Ala291, Leu287, Ile352, His447 | - | |
| 范德华接触 | 其他结合口袋残基(来源中未单独列出) | - | 支持类固醇骨架的容纳 |
表4:石胆酸与NR1H4/FXR(PDB ID: 3DCT)对接生成的结合模式。
来源:分子可视化和二维相互作用图工具2D配体-残基相互作用图,如论文结果部分(分子对接)所述。“-”表示该接触的距离值未单独报告。
| 相互作用类型 | 残基 | 距离 (Å) | 备注 |
| 传统氢键 | His447 | 3.66 | |
| 氢键 | Gly322 | 3.46 | |
| π-阴离子相互作用 | Val325 | 4.96 | |
| 碳-氢键 | Trp469 | 4.51 | |
| 不利的供体-供体接触 | Arg395 | 3.89 | |
| 不利的供体-供体接触 | Gln396 | 3.40 |
表5:熊去氧胆酸与NR1H4/FXR(PDB ID: 3DCT)对接生成的结合模式。 来源:分子可视化及二维相互作用图工具生成的2D配体-残基相互作用图,如论文结果部分(分子对接)所述。“-”表示该相互作用接触的距离值未单独报告。
| 相互作用类型 | 残基 | 距离 (Å) | 备注 |
| 常规氢键 | Gln161 | 3.78 | |
| 常规氢键 | Gln161 | 4.19 | 第二个接触点 |
| 常规氢键 | Ser120 | 3.95 | |
| 常规氢键 | Arg207 | 3.05 | |
| 常规氢键 | Arg207 | 3.77 | 第二个接触点 |
| π-阳离子相互作用 | Arg207 | - | |
| π-供体氢键 | Cys163 | - | |
| π-烷基 / 范德华接触 | Arg64, Ala162, His121, Ser205, Trp206 | - |
表6:鞣花素A与CASP3(PDB ID: 2DKO)对接生成的结合模式。 来源:分子可视化及二维相互作用图工具“2D配体-残基相互作用图”,如论文结果部分(分子对接)所述。“-”表示该相互作用接触的距离值未单独列出。
| 相互作用类型 | 残基 | 距离 (Å) | 备注 |
| 静电盐桥 | Asp155 (D3.32) | - | 色胺的质子化胺基 |
| 氢键 | Thr160 | - | |
| 氢键 | Ser159 | - | |
| 芳香环接触 | Phe340, Trp336 | - | |
| π-烷基相互作用 | Val156, Ile163 | - | |
| 范德华接触 | Tyr370, Phe339, Ser242, Phe243, Phe332, Leu123 | - |
表7:色胺与HTR2A(PDB ID: 6A93)对接生成的结合模式。 来源:分子可视化及二维相互作用图工具“2D配体-残基相互作用图”,如文稿结果部分(分子对接)所述。“-”表示该相互作用接触的距离值未单独报告。
| (A) 重新对接验证(阳性对照) | ||||||
| PDB ID | 蛋白质 | 共结晶配体 | Vina 评分 (kcal/mol) | RMSD (Å) | 阈值 (Å) | 结果 |
| 1DB1 | VDR | VDX (维生素D类似物) | −13.0 | 0.87 | 2.0 | PASS |
| 3DCT | FXR | WAY-362450 (064) | −11.9 | 1.79 | 2.0 | PASS |
| (B) 交叉对接验证(阴性对照) | ||||||
| 配体 | 同源靶标 (PDB) | 同源评分 (kcal/mol) | 非同源靶标 (PDB) | 非同源评分 (kcal/mol) | 差值 (kcal/mol) | 选择性 |
| Lithocholic acid | VDR (1DB1) | −10.0 | CASP3 (2DKO) | −8.3 | 1.7 | 已确认 |
| Tryptamine | HTR2A (6A93) | −7.1 | VDR (1DB1) | −6.4 | 0.7 | 中等(在Vina不确定性范围内 ±0.5–1.0) |
表8:对接协议验证结果:重新对接的RMSD值(阳性对照)和交叉对接评分(阴性对照)。 数据来源:Docking_Validation/Results/Docking_Validation_Results.xlsx 和 Docking_Validation/Logs/*.log(代谢物配体与靶标蛋白的分子对接,搜索穷尽性 = 8,随机种子 = 42,对接盒子尺寸为 25 Å × 25 Å × 25 Å)。RMSD 值通过重原子按原子名称匹配计算(无叠加)。
| 复合物 | RMSD (nm),均值 ± 标准差(范围) | Rg (nm),均值 ± 标准差(范围) | SASA (nm^2),均值 ± 标准差(范围) | 氢键数,均值 ± 标准差(范围) | RMSF (nm),均值(最大值) |
| LCA-VDR/1DB1 | 0.230 ± 0.025 (0.167–0.296) | 1.889 ± 0.009 (1.863–1.919) | 130.4 ± 2.3 (122.3–137.4) | 1.9 ± 0.9 (0–7) | 0.093(最大值 0.600,位于残基 120) |
| LCA-NR1H4/FXR/3DCT | 0.190 ± 0.020 (0.135–0.281) | 1.824 ± 0.008 (1.804–1.849) | 129.7 ± 2.3 (121.9–138.1) | 3.8 ± 0.7 (1–6) | 0.113(最大值 0.298) |
| UDCA-NR1H4/FXR/3DCT | 0.190 ± 0.020 (0.135–0.281) | 1.834 ± 0.013 (1.809–1.921) | 131.0 ± 3.4 (121.6–143.5) | 1.1 ± 1.1 (0–5) | 0.113(最大值 0.298) |
| Urolithin A-CASP3/2DKO | 0.521 ± 0.058 (0.244–0.755) | 1.892 ± 0.024 (1.839–1.984) | 134.9 ± 3.0 (126.4–146.4) | 0.6 ± 0.7 (0–3) | 1.172(最大值 2.532,位于残基 175) |
| Tryptamine-HTR2A/6A93 (membrane) | 0.177 ± 0.017 (0.131–0.227) | 2.089 ± 0.007 (2.070–2.116) | 165.1 ± 2.7 (156.–172.7) | 1.7 ± 0.7 (0–4) | 0.090(最大值 0.319) |
表9:五种优先级蛋白-配体复合物(包括膜嵌入的色胺-HTR2A系统)200纳秒分子动力学模拟行为的汇总。 数据来源:分子动力学轨迹分析工具(.xvg文件)——gmx rms、gmx gyrate、gmx sasa、gmx hbond、gmx rmsf——根据方案步骤8.8,对每次200纳秒生产阶段模拟的最后150纳秒(50–200纳秒)进行计算。RMSD/Rg基于主链拟合;SASA探针半径为0.14 nm;氢键供体-受体截断距离为0.35 nm / 30°。LCA-3DCT和UDCA-3DCT共享一条蛋白主链轨迹(RMSD、RMSF),但具有配体特异性的Rg/SASA/氢键数据。
| 色胺-HTR2A/6A93(膜)— 按残基定量分解 | ||
| 残基 | 总ddG贡献值(kcal/mol),均值 ± 标准差 | 作用方向 |
| Asp155 (D3.32) | −89.94 ± 6.81 | 稳定作用(主导) |
| 色胺(配体) | −13.01 ± 6.22 | 稳定作用 |
| Tyr171 | 13.62 ± 4.54 | 去稳定作用 |
| Val167 | 23.32 ± 3.96 | 去稳定作用 |
| Val156 | 20.03 ± 3.81 | 去稳定作用 |
| Thr160 | 4.86 ± 3.64 | 去稳定作用 |
| Ser159 | 24.16 ± 3.48 | 去稳定作用 |
| Ser86 | 24.48 ± 3.65 | 去稳定作用 |
| Phe87 | 35.18 ± 4.04 | 去稳定作用 |
| Phe133 | 32.80 ± 3.70 | 去稳定作用 |
| Phe140 | 30.63 ± 3.84 | 去稳定作用 |
| Phe141 | 35.25 ± 3.55 | 去稳定作用 |
| Ile163 | 27.64 ± 3.71 | 去稳定作用 |
| Trp137 | 53.77 ± 4.32 | 去稳定作用(最不利) |
| 其余四个复合物 — 按残基分解鉴定出的残基(定性) | ||
| 复合物 | 残基 | 作用方向 |
| LCA-VDR/1DB1 | 配体(LCA) | 有利作用 |
| LCA-VDR/1DB1 | Gln317 | 有利作用 |
| LCA-VDR/1DB1 | Trp286 | 不利作用 |
| LCA-NR1H4/FXR/3DCT | Arg331 | 有利作用(主导) |
| UDCA-NR1H4/FXR/3DCT | Glu326 | 混合/分布网络 |
| UDCA-NR1H4/FXR/3DCT | Asp394 | 混合/分布网络 |
| UDCA-NR1H4/FXR/3DCT | Arg395 | 混合/分布网络 |
| UDCA-NR1H4/FXR/3DCT | Arg441 | 混合/分布网络 |
| UDCA-NR1H4/FXR/3DCT | Asp470 | 混合/分布网络 |
| 尿石素A-CASP3/2DKO | Arg64 | 强有利作用(极性/静电) |
| 尿石素A-CASP3/2DKO | Arg207 | 强有利作用(极性/静电) |
表10:每残基MM-PBSA分解简要摘要:五个优先考虑的蛋白-配体复合物中具有稳定或去稳定作用的残基(绝对贡献值 ≥ 0.5 kcal mol⁻1),包括嵌入膜中的色胺-HTR2A体系。 数据来源:Membrane Simulation/03_MMPBSA/results/FINAL_DECOMP_MMPBSA.dat(分子力学/连续溶剂结合能计算工具广义玻恩模型(GB)的每残基分解结果,“复合物:总能量分解”)。残基编号已从CHARMM-GUI构建系统的内部编号(偏移量+68)转换为本论文其他部分所使用的原始PDB编号6A93。
来源:先前的分子动力学模拟数据/1DB1,2KD0, LCA & UDCA_3DCT}/mmpbsa_*/分解_NORMAL_GB_复合物_TDC*.svg 及论文结果(MM-PBSA 结合自由能及按残基分解结果)。这四个复合物在项目目录中没有数值化的按残基 .dat/.csv 输出文件(仅有渲染后的 SVG 图像,其文本为矢量路径形式,无法机器提取);仅在论文文本中报告了残基身份及其贡献的有利或不利方向。这四个复合物的精确 kcal/mol 贡献值在源代码仓库中不可获得。
这项探索性计算研究展示了一种整合的、可重复的工作流程,用于优先筛选与微生物代谢物相关的宿主基因及蛋白-配体复合物,并在此应用于一个公开的便秘型肠易激综合征(IBS-C)直肠黏膜转录组数据集。通过该工作流程,一组预测的与微生物代谢物相关的基因与数据集中持续下调的基因重叠,这些基因富集于G蛋白偶联受体(GPCR)、5-羟色胺能、钙信号传导、神经活性配体-受体以及核受体相关通路中,这些系统在微生物群-宿主通信中的作用日益受到关注48,49。这些发现应严格视为假设生成性的结果:本分析并未测量微生物代谢物浓度、受体蛋白丰度、配体结合、受体激活、下游信号传导、肠道运动、分泌功能、疼痛反应或临床结局。最有力的可支持结论是,所鉴定出的基因和通路是值得进一步实验验证的候选对象,而非已被证实的疾病机制。
相较于以往仅孤立研究单一代谢物-受体配对的工作,本方案的关键意义在于将靶点预测、公共转录组学数据、网络分析、分子对接(含验证对照)、分子动力学模拟以及MM-PBSA方法整合为一个连续的优先级排序流程。每个阶段均对前一阶段产生的候选结果进行筛选和生物学背景化,这种逐级过滤机制使得最终的候选列表在实验上具有可操作性。本文鉴定出的以GNAQ为核心的GPCR模块以及与5-羟色胺受体相关的模块在生物学上具有合理性:Gq信号通路在磷脂酶C激活、肌醇1,4,5-三磷酸生成、钙离子动员、分泌过程及肠内分泌功能中发挥重要作用;短链脂肪酸和色氨酸代谢产物介导的信号通路在黏膜稳态中具有明确功能,而5-羟色胺能信号在胃肠道运动、分泌、内脏敏感性及肠-脑通讯中亦有公认作用10,18,50,51。
本研究的一个关键方法学特征是对膜嵌入型受体 HTR2A 的处理。由于溶液相模拟无法重现调控 G 蛋白偶联受体构象行为的脂质环境,因此在明确的 POPC 双分子层中对色胺-HTR2A 复合物进行了模拟。在此膜环境中,受体在整个 200 ns 的轨迹过程中保持结构稳定,且色胺的铵基与 Asp155(D3.32)之间的盐桥在几乎整个模拟过程中均得以维持。三种独立证据——对接构象、轨迹中持续的接触距离以及每残基 MM-PBSA 分析的主要贡献——均指向相同的保守 D3.32 相互作用,这为预测的色胺结合模式提供了内在一致性,该模式重现了胺类配体在5-羟色胺受体上的典型结合几何构型。
在重复本工作流程时,需考虑某些方法学问题。经典结构错误或泛测定干扰化合物(PAINS)标记的化合物会通过靶点预测和分子对接传播,因此需要准确选择代谢物并进行化学信息学的审慎整理。应始终应用置信度标准(化学-蛋白质相互作用靶点预测:≥0.700;分子对接程序:≥0.70;蛋白质-蛋白质相互作用网络构建与通路富集:≥0.700),以减少由噪声驱动产生的靶点集合。预测的靶点应按功能类别进行分组,以避免将所有与代谢物相关的基因误判为受体。PDB结构的精确预处理、配体能量最小化以及在已知结合残基周围合理设置网格是对接过程中的关键环节;本文引入的自对接和交叉对接对照可为对接方法的正确性提供客观评估。分子动力学的可重复性范围由力场参数化、适当的溶剂化或膜结构构建、分阶段平衡以及充分的生产阶段采样共同决定。
常见的调整和故障排除步骤包括:如果目标预测未返回任何结果,可放宽阈值;对于具有多个探针的基因,需在探针水平检查方向一致性;将孤立的蛋白质-蛋白质相互作用网络构建和通路富集节点解释为依赖于阈值的结果,而非生物学上无关的信号。对于膜受体,应使用明确的脂质双分子层模拟而非水相模拟,本文所述的HTR2A方法即为范例。当需要进行残基级别的能量分解时,必须使用具备分解功能的计算引擎,并将报告的残基编号与天然受体的编号对应,以避免歧义。我们建议将通路富集结果视为候选基因列表的组织性背景,而非通路水平的验证依据。从机制上看,只要基因列表中包含多个5-羟色胺受体基因,无论蛋白质水平是否存在共调控,都可能导致G蛋白偶联受体(GPCR)、5-羟色胺能或钙信号通路术语的富集。在分析嵌入膜中的HTR2A体系的RMSD、Rg和RMSF值时,应考虑脂质双分子层的影响:轨迹后期Rg的降低可能反映的是双分子层驱动的跨膜结构域构象适应,而非整体去折叠;配体-蛋白质间持续存在的氢键应结合整体RMSD的稳定性一并解读。
本研究的局限性较大,限制了结果的解释。研究基于一个相对较小的公共数据集,且在主要的公共转录组数据库(基于网络的差异基因表达分析工具和 ArrayExpress)中未发现与本研究设计和平台相当的独立肠易激综合征便秘型(IBS-C)直肠黏膜转录组数据集,可用于分析时的重复验证队列。缺乏独立的转录组重复验证是一项主要局限性,本文中的任何结论均不应被解读为对单数据集发现结果的外部验证。该数据集显示出近乎普遍的差异表达(约 94.5% 的基因具有显著性,其中绝大多数为下调),这一特性使得传统的富集统计方法失去信息价值,无法得出关于靶基因下调相对于基因组背景特异性的结论;因此,重叠结果仅作为描述性的方向性模式报告,而非统计学上的富集。 bulk 黏膜转录组学无法区分真实的基因调控与细胞组成变化。mRNA 表达水平不能决定蛋白质丰度或功能反应。靶点预测数据库存在注释偏差,而分子对接、分子动力学(MD)和 MM-PBSA 的结果依赖于力场的选择、配体参数化、起始位置、模拟时间以及采样充分性。部分网络服务器和软件包组件的精确版本/构建标识未知,包括 CHARMM 兼容配体参数化服务、CHARMM-GUI 和统计计算环境
软件包构建以及分子力学/连续溶剂结合能计算工具的子版本无法从存档的项目记录中完全恢复,应注明在单独的材料表中提供。MM-PBSA 值为相对估算值,不包含显式的构型熵项,不应被解释为实验测得的亲和力。本研究缺乏代谢组学数据,无法确定在便秘型肠易激综合征(IBS-C)中配体的可及性是否发生改变,也无法判断所观察到的表达变化是病因、结果、代偿性反应还是无关的相关性。
该方法的未来应用应包括独立的转录组学重复验证、定量聚合酶链式反应(qPCR)和蛋白质水平验证,通过单细胞或空间转录组学进行细胞类型定位,对相关代谢物类别进行代谢组学分析,以及在患者来源的结肠类器官、黏膜外植体或类似模型中开展功能性配体-反应检测。与腹泻型肠易激综合征、混合型肠易激综合征、炎症性肠病及非肠易激综合征型便秘队列1,2的比较将有助于确立疾病的特异性。在结构研究方面,重复分子动力学(MD)轨迹模拟、采用不同初始构象进行敏感性分析,并完整记录拓扑结构、轨迹以及MM-PBSA输入和输出文件的存放信息,将进一步增强研究的可重复性。实验性的配体-反应检测仍是确定优先结合复合物是否具有功能相关性的必要手段;目前的结果尚不支持任何临床或治疗方面的主张。
作者声明无利益冲突。
本研究未获得外部资金支持。谨此感谢GSE36701数据集以及STITCH、SwissTargetPrediction、SwissADME、STRING、RCSB蛋白质数据库、基因表达综合数据库、CHARMM-GUI和膜蛋白取向(OPM)资源的公开获取,同时感谢AutoDock Vina、GROMACS、CHARMM36m、CGenFF、gmx_MMPBSA、Open Babel、PyMOL和Discovery Studio Visualizer软件的公开使用。
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| AutoDock Vina | Scripps Research / 开源 | v1.2.7;https://vina.scripps.edu/ | 代谢物配体与靶蛋白的分子对接。 |
| CGenFF/ParamChem | SilcsBio / 马里兰大学 | v4.6;https://cgenff.com/ | 用于分子动力学模拟的配体力场参数化。 |
| CHARMM36m 力场 | CHARMM 开发团队 / 开源 | CHARMM36m;https://www.charmm.org/charmm/resources/charmm-force-fields/ | 用于分子动力学模拟的蛋白质力场。 |
| CHARMM-GUI 膜构建器 | CHARMM-GUI / 利哈伊大学 | 网络服务器;具体版本不可追溯;https://www.charmm-gui.org/?doc=input/membrane | 显式POPC膜系统的构建与平衡设置。 |
| Discovery Studio Visualizer | BIOVIA (Dassault Systèmes) | 2021;https://discover.3ds.com/discovery-studio-visualizer-download | 二维配体-残基相互作用分析。 |
| GEO2R | NCBI 基因表达综合数据库 | 网络工具;访问时间:2026年1月至5月;https://www.ncbi.nlm.nih.gov/geo/geo2r/ | GSE36701 的差异表达分析。 |
| GeneCards | 魏茨曼科学研究所 | 网络数据库;访问时间:2026年1月至5月;https://www.genecards.org/ | 在靶点标准化过程中进行基因符号和基因信息验证。 |
| gmx_MMPBSA | 开源(Valdés-Tresanco 等) | 1.5.x;https://valdes-tresanco-ms.github.io/gmx_MMPBSA/ | MM-PBSA 结合自由能估算及残基分解分析。 |
| GROMACS | GROMACS 开发团队 / 开源 | 2024.2;https://www.gromacs.org/ | 分子动力学模拟引擎。 |
| GSE36701 转录组数据集 | NCBI 基因表达综合数据库 | GSE36701;https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36701 | 公开的肠易激综合征便秘型(IBS-C)直肠黏膜表达数据集。 |
| Open Babel | 开源 | 3.2.0;https://openbabel.org/ | 化学文件格式转换、三维配体生成及配体预处理。 |
| OPM 数据库 | 密歇根大学 | 网络数据库;访问时间:2026年1月至5月;https://opm.phar.umich.edu/ | 用于对齐 HTR2A 的膜蛋白取向坐标数据库。 |
| ParmEd | ParmEd 开发团队 / 开源 | 4.x;https://parmed.github.io/ParmEd/html/index.html | 氢原子质量重分配及分子模拟拓扑结构处理。 |
| PyMOL | Schrödinger / 开源 | 2.x;https://www.pymol.org/ | 三维结构可视化及受体-配体图像绘制。 |
| RCSB 蛋白质数据库 | RCSB PDB | 网络数据库;访问时间:2026年1月至5月;https://www.rcsb.org/ | 实验蛋白质结构及PDB坐标的来源。 |
| STITCH | STITCH 联盟(EMBL) | v5.0;https://stitch.embl.de/ | 化学-蛋白质相互作用靶点预测。 |
| STRING | STRING 联盟 / ELIXIR | v12.0;https://version-12-0.string-db.org/ | 蛋白质-蛋白质相互作用网络构建与通路富集分析。 |
| SwissADME | SIB 瑞士生物信息学研究所 / 洛桑大学 | 网络工具;访问时间:2026年1月至5月;https://www.swissadme.ch/ | 化学信息学描述符、药代动力学预测及PAINS评估。 |
| SwissTargetPrediction | SIB 瑞士生物信息学研究所 / 洛桑大学 | 网络工具;访问时间:2026年1月至5月;https://www.swisstargetprediction.ch/ | 基于配体的人类蛋白质靶点预测。 |
| UniProt ID 映射 | UniProt 联盟 | 网络服务;访问时间:2026年1月至5月;https://www.uniprot.org/id-mapping | 将蛋白质标识符映射为标准化的HGNC批准基因符号。 |
| NVIDIA RTX 3080 | NVIDIA 公司 | RTX 3080;≥8 GB 显存;论文中未指定 CUDA/驱动版本 | 支持 CUDA 的图形处理器,用于分子动力学模拟。 |
| CUDA 兼容 GPU | NVIDIA 公司 | 论文中未指定 CUDA 工具包版本;≥8 GB 显存 | 支持 CUDA 且显存 ≥8 GB 的 GPU;工作站还需配备 ≥32 GB 内存和 6 核 CPU。 |
| Ubuntu Linux | Canonical Ltd. / 开源 | 22.04 LTS | 64 位 Linux 操作系统。 |
| Python 3.9 | Python 软件基金会 | 3.9 | 用于工作流脚本编写和数据分析的通用编程环境。 |
| 基因表达综合数据库(GEO) | NCBI / 美国国家医学图书馆 | 公共网络存储库;论文中未指定软件版本 | 公共功能基因组学数据存储库。 |
| AutoDockTools/MGLTools | 分子图形学实验室,Scripps Research | 1.5.7 | 分子结构与对接输入准备工具包。 |
| GROMACS 分析工具 | GROMACS 开发团队 / 开源 | 2024.2 | 分子动力学轨迹分析工具。 |
| CHARMM-GUI 六步协议 | CHARMM-GUI / 利哈伊大学 | 网络协议;具体版本不可追溯 | 基于网络的多阶段分子系统准备与平衡工作流程。 |
| cgenff_charmm2gmx_py3.py | 开源转换脚本;论文中未注明来源 | 论文中未指定版本 | 力场拓扑结构转换脚本。 |
| Python | Python 软件基金会 | 3.9 | 通用编程环境。 |
| SciPy | SciPy 社区 / 开源 | 论文中未指定版本 | 科学计算库。 |
| scipy.stats.fisher_exact | SciPy 社区 / 开源 | 论文中未指定 SciPy 版本 | Fisher’s 精确检验实现。 |
| R | 统计计算基金会 | 4.3.x | 统计计算环境。 |
| Bioconductor | Bioconductor 项目 / 开源 | 3.18 | 生物信息学软件框架。 |
| limma | Bioconductor 项目 / 开源 | 论文中未指定版本 | 差异基因表达分析软件包。 |
| Benjamini–Hochberg 方法 | 统计方法 | 不适用(统计过程) | 错误发现率校正方法。 |
| NVIDIA RTX 3080 | NVIDIA 公司 | RTX 3080;≥8 GB 显存;论文中未指定 CUDA/驱动版本 | 至少配备 8 GB 显存的图形处理器。 |
| CUDA 兼容 GPU | NVIDIA 公司 | 论文中未指定 CUDA 工具包版本;≥8 GB 显存 | 支持通用并行计算的图形处理器。 |
| Ubuntu Linux 22.04 LTS | Canonical Ltd. / 开源 | 22.04 LTS | 64 位 Linux 操作系统。 |
| TIP3P | CHARMM 力场开发团队 / 开源 | TIP3P;无适用软件版本 | 三站点显式水模型。 |
| MM/PBSA | gmx_MMPBSA 开发团队 / 开源 | gmx_MMPBSA 1.5.x | 分子力学/泊松–玻尔兹曼表面面积结合能方法。 |
申请许可以重复使用本 JoVE 文章的文本或图表
申请许可