Nature | AI引导设计CRE顺式调控元件调节细胞类型特异性基因表达

2024年10月23日，美国巴港杰克逊实验室的 Ryan Tewhey 教授团队在Nature上发表了题为 Machine-guided design of cell-type-targeting cis-regulatory elements 的研究论文。研究通过使用来自于大规模平行报告基因测定（MPRA）技术获得高通量基因组实验数据来训练AI模型 Malinois，旨在从任何序列中预测细胞类型特异性的CRE活性。在此基础上，该团队提出了DNA活性的计算优化平台（CODA），一种从头设计新型合成CRE的方法，可以在三种转化细胞系中驱动细胞类型特异性转基因表达，并且在体内测试时在类似组织中实现了特异性。合成序列表现出与靶细胞类型中的活性以及同时降低脱靶细胞活性相关的独特基序词汇。

顺式调节元件（CRE）是调节基因表达的关键非编码DNA序列。它们控制基因的时间和空间表达，它为精准医疗开辟的新的可能性。目前有许多技术可以用于识别具有调节活性的CRE，例如全基因组关联研究（GWAS）。GWAS首先从队列中收集DNA和表型信息，然后进行全基因组测序以识别整个基因组中的变异位点。然后进行关联统计检验，通常使用线性回归来检验基因突变和表型之间的相关性。最后，结果在独立的验证队列中复制，以确认可靠的变异-表型关联（如图1所示）。GWAS可以筛选与常见疾病和性状相关的大规模遗传变异，包括非编码变异。然而，从众多变异位点中识别具有调控活性的变异是一个重大挑战。

荧光素酶报告基因测定是一种常用的方法，用于检测转录因子与目标基因启动子DNA之间的相互作用。首先，将要研究的基因调控元件插入报告基因的上游，构建报告质粒。然后，将转录因子质粒和报告基因质粒共同转染到细胞中。如果转录因子激活目标启动子，荧光素酶基因将被表达。荧光素酶的表达水平与转录因子活性的强度成正比。通过测量荧光强度可以评估荧光素酶活性，从而判断转录因子是否与目标启动子片段相互作用。双荧光素酶报告系统增加了一个额外的报告基因系统作为内对照，以消除实验过程中产生的组间误差。然而这种方法耗时且通量较低。

因此，研究者们开发了一种高通量分析方法，称为大规模平行报告基因测定（MPRA），用于研究成千上万个基因组调控元件的活性。MPRA在传统的荧光素酶报告基因中引入了一个独特的条形码。通过RNA-seq对转染DNA和相应mRNA的条形码进行测序，基因mRNA与DNA条形码读取比率分析CREs的活性。该方法可用于大规模分析，以验证GWAS的结果。

进一步考虑到致病性变异可能具有细胞特异性效应，科学家开发了单细胞大规模平行报告基因测定方法（scMPRA）。在MPRA的基础上，这种方法运行每个CREs同时标记有CRS特异性条形码cBC和细胞身份条形码rBC。通过量化这些条形码，可以获得CRS在单细胞中的表达，并与细胞身份相关联。这使得研究细胞类型特异性的基因调控成可能。

然而，自然进化生成的序列仅代表了一小部分，可能不足以用以用于治疗应用。那么，我们能否学会读取和编写这些调控元件的代码呢？

本研究通过使用来自MPRA的高通量基因组实验数据训练AI模型，提出了一个潜在的解决方案。该方法合成了成千上万种新的CREs，这些元件能够准确控制不同细胞类型中的基因表达，并且比天然CREs表现出更强的细胞类型特异性。

研究人员首先将MPRA与卷积神经网络结合，利用MPRA提供训练AI模型所需的大规模CRE序列数据集创建了Malinois模型，旨在从任何序列中预测细胞类型特异性的CRE活性。使用的数据来自数十万个人体DNA序列，包括血液（K562）、肝脏（HepG2）和大脑（SK-N-SH）细胞中的CRE活性。Malinois在K562、HepG2和SK-N-SH细胞中预测结果与实测活性结果高度相关（皮尔逊相关系数r=0.88-0.89； 斯皮尔曼等级相关系数p=0.81-0.83；所有P值均小于10^-300），并且在细胞特异性估计上与实验结果相当。

研究人员随后开发了CODA（Computational Optimization of DNA Activity, DNA活性计算优化），这是一个用于设计具有编程功能的新型CRE的模块化平台。CODA遵循一个迭代循环：预测序列活性、使用目标函数量化序列和设计目标的契合度，然后更新序列以提高目标值。CODA通过计算MinGap来量化效果—即目标细胞类型中的预测MPRA活性与两个非目标细胞类型中的最大活性之间的差异。MinGap越大，表示非目标效应越小。除了使用MinGap来评估特异性外，还采用了一个径向坐标系来量化和可视化特异性。最具特异性的细胞沿着三个类型轴中的一个轴向外延伸，而各细胞类型中活性均匀的序列则趋向原点。此外，在CODA中，序列更新可以通过不同的序列设计算法进行控制。这些算法在设计细胞类型特异性元素方面的整体能力通过对超参数的选择具有较多的鲁棒性。然而，可以通过调整来平衡最大化目标与保持设计元素中k-mer多样性之间的权衡。

为了实证测试CODA的有效性，研究人员进了MPRA实验以测量合成序列的活性。对于每种细胞类型， 研究人员从CODA中的三个序列设计算法中各生成4,000个细胞类型特异性序列，共计36,000个合成候选序列。当最大化特异性时，Malinois倾向于强烈偏好某些序列基序，研究者们通过惩罚这些高频偏好基序的使用来鼓励CODA减少它们的出现。由CODA设计的合成序列中，大多数（94.1%）成功实现了细胞类型特异性。

在发现合成CRE比天然序列更具细胞类型特异性后，研究人员试图将序列内容与负责的转录因子（TF）联系起来。转录在一定程度上是通过个别转录因子与序列基序的结合以及转录因子之间的相互作用来调控。首先，研究人员使用Malinosis预测了每个序列在三种细胞类型中的核苷酸分辨率活性贡献分数，采用的事集成梯度法。接着使用TF-MoDiSco方法识别了66个由贡献分数指导的基序模式，从中提取出36个不冗余的核心基序（7-18bp），这些基序在MPRA测试库中富集，其中31个与已知的人类转录因子结合基序高度对齐。

进一步，研究人员检查了天然序列与合成序列在基序上的差异，使用基于贡献分数的基序扫描方法。这36个核心基序至少在和合成序列和天然序列中各出现一次，表明两类序列共享一个词汇。然而，基序的使用频率存在差异，合成序列中转录激活因子（GATA和HNF4A）和抑制因子（如MEIS2和GFI1B）的使用频率较倜然序列更高。词汇分析显示，合成序列由更多独特基序组成，且基序总数更多。此外惩罚和合成序列在非冗余基序使用上优于未惩罚的和合成序列。

对36个核心基序的调控活动贡献进行分析，28个基序对序列活性做出正向贡献，8个为负向贡献。包括GATA、HNF1B、HNF4A等已知激活因子在特定细胞类型中具有显著正向贡献，而GFI1B和MEIS2等抑制因子则对特定细胞类型有负向贡献。所有基序在随机背景中和替换为随机序列时表现出与其预测贡献一致的效应。

研究人员评估了合成的CREs的特异性是否能够超越用于设计的最初三种细胞系。为了确定在添加新的细胞系时特异性是否得以保留，研究人员为A549（肺上皮癌）和HCT116（结肠上皮癌）细胞训练了额外的模型，观察到合成CREs在其目标细胞类型中的最大预测活性保持不变，特别是使用Fast SeqProp生成的CREs，尽管最小间隙MinGap有所减少。

研究人员通过斑马鱼和小鼠的体内实验验证了几种合成CRE序列的有效性。在转基因斑马鱼中，CODA设计的HepG2细胞特异性CREs的活体成像结果显示，2/3的合成序列在发育中的肝脏中表现出强烈且一致的表达。为神经元设计的序列也显示出类似的成功。

此外，研究人员评估了在斑马鱼中表现出特异性的两个神经元序列的活性是否能够扩展到哺乳动物小鼠模型系统。结果表明，合成神将元CRE 1 （synN1）的表达仅限于发育中的大脑皮层，未在其他部位观察到表达。

总体而言，这项研究提出了一种全新的策略，用于从头设计新的合成CREs。首先，研究团队开发了Malinois，一个深度卷积神经网络，用于预测任意序列的细胞类型特异性CRE活性。在此基础上，团队引入了CODA，一种从零开始设计新型合成CREs的方法。随后，他们在三种不同的细胞系中测试并解释了功能性序列内容。最后，团队在体外评估了合成CREs的功能性和有效性。