顺式调控元件（CREs）控制基因表达，并且其结构和功能是动态变化的，反映了不同效应蛋白组成随时间的变化。然而，测量基因组中CREs处效应蛋白组织的方法有限，阻碍了将CRE结构与其在细胞命运和疾病中的功能联系起来的努力。在此研究中，我们开发了PRINT，这是一种计算方法，可以从批量和单细胞染色质可及性数据中识别DNA-蛋白质相互作用的足迹。利用这些多尺度足迹，我们创建了seq2PRINT框架，该框架使用深度学习来精确推断转录因子和核小体结合，并解释CREs的调控逻辑。我们将seq2PRINT应用于人类骨髓的单细胞染色质可及性数据，观察到造血过程中围绕先驱因子中心的CREs的逐步建立和扩展。此外，我们还发现了小鼠造血干细胞中与年龄相关的CRE结构变化，包括核小体足迹的广泛减少和新识别的Ets复合基序的增加。总体而言，我们建立了一种从染色质可及性数据中获得DNA结合蛋白动态丰富见解的方法，并揭示了分化和衰老过程中调控元件的结构。
为了克服Tn5转座酶序列偏差对足迹检测的显著干扰，我们训练了一个卷积神经网络，基于来自细菌人工染色体（BACs）的去蛋白化DNA上的Tn5插入数据（扩展数据图1a-d和补充表1）。该模型在高GC含量区域表现出色，特别是在提取的人类基因组DNA数据上的表现（R=0.92），优于Tn5偏差校正的ChromBPNet（扩展数据图1e-k）。我们为人类基因组和常见模式生物提供了预计算的Tn5偏差预测，以及一个预训练的深度学习模型作为领域的资源（“数据可用性”）。
PRINT可以在多种蛋白质大小范围内识别DNA结合蛋白的足迹。我们开发了一种统计方法，量化相对于给定位置估计背景分散的观察到的Tn5插入缺失的显著性，以得出足迹评分（方法）。与先前的足迹检测方法相比，这种方法在去蛋白化DNA上减少了数量级的假阳性检测（扩展数据图1l-v）。受早期使用MNase或ATAC-seq片段大小推断不同大小DNA结合蛋白的方法启发，我们在4-200 bp的窗口范围内计算了足迹评分。我们在体外验证了这种方法，在去蛋白化DNA上孵育了纯化的MYC/MAX或CEBPA。只有在存在纯化转录因子的情况下，才在转录因子基序位点检测到强足迹，而背景信号非常低，而一种已建立的ATAC-seq足迹检测方法并未区分前景和背景（图1c-f和扩展数据图2a-d）。此外，我们发现在较高浓度（100 vs 50 nM）的MYC/MAX下，低亲和力位点的足迹增加，表明足迹评分对特定位点的转录因子占有率敏感（图1g）。
我们发现PRINT可以检测哺乳动物细胞中的足迹。我们观察到了与核小体和特定转录因子相对应的不同足迹模式（图1h-i）。转录因子结合模式可以分为四个代表性类别（扩展数据图2e）。与之前使用DNase I的研究一致，我们发现不同转录因子之间的足迹强度有所不同，包括一些未留下可检测足迹的转录因子，这可能是由于弱或瞬时结合DNA所致。也可以检测到抑制性转录因子的足迹（扩展数据图2f）。我们通过与ChIP-exo数据进行基准测试验证了足迹，发现在转录因子结合位点上的一致性，以及ChIP-exo可能的假阴性（扩展数据图2g-h）。总之，我们的结果显示，多尺度足迹分析与PRINT可以稳健地检测许多不同的DNA结合蛋白。
接下来，我们试图使用多尺度足迹来预测特定蛋白质与DNA的结合。我们设计了预测转录因子和核小体结合的模型（图2a、方法和补充说明）。该模型仅以局部DNA序列作为输入，就能预测核小体和转录因子足迹（图2b）。我们在HepG2细胞的ATAC-seq数据中观察到，预测的多尺度足迹与观察到的足迹之间整体相关性为0.75，对读深的子采样具有鲁棒性（扩展数据图4b）。
我们提取了seq2PRINT学习到的关键序列特征，发现由此产生的逐碱基DNA序列归因分数使我们能够剖析CRE内的转录因子结合架构。在一个示例位点中，相对于整个CRE计算的归因分数突出了跨越该区域的短序列，这些序列与转录因子基序位置重叠，而针对特定足迹对象计算的分数则突出了特定基序（图2b,c）。对于足迹α和γ，序列模型识别出支撑该足迹的基序，而对于足迹β，则检测到足迹下方的NFE2L2基序及其邻近的NFYB基序，表明附近转录因子之间可能存在结合协调。足迹δ可能代表一个核小体，并由附近的转录因子NRF1和NFYB基序预测，显示出更长距离的依赖关系，并揭示了最相关的核小体定位因素。后两个例子进一步表明，一些缺乏明显足迹的转录因子（例如NFYB）可以通过seq2PRINT检测到，可能是由于对相邻元素的影响，而且这种方法可用于建模CRE内的DNA结合蛋白之间的相互作用。
seq2PRINT可用于预测全基因组范围内的转录因子结合。我们使用来自seq2PRINT的序列归因分数生成了一个训练用于预测ChIP-seq数据的转录因子结合评分（方法）。转录因子结合评分能够以高精度预测转录因子结合，并优于先前的方法（图2d）。我们甚至能够预测那些弱或无直接足迹的转录因子，而其他方法对此类转录因子的表现特别低（图2e，扩展数据图2e和4c-f，补充表2）。为了研究模型如何检测这些转录因子，我们通过打乱其基序序列模拟了转录因子结合的丧失。模型预测了包括核小体在内的跨尺度邻近足迹的变化，这些变化与实验确定的足迹变化（包括CTCF降解诱导的CTCF耗尽）非常吻合（聚合R=0.93，每位点中位数R=0.66；图2f）。类似的結果也在地塞米松诱导的糖皮质激素受体重定位和DNA结合以及IFN诱导的Stat2结合DNA中获得（扩展数据图4g-j）。我们广泛应用了这种方法，并识别出对核小体有强影响的转录因子（例如JUN和YY1）或在同一位置的转录因子（例如ZKSCAN1）（图2g）。预测的足迹评分变化在转录因子家族中高度相似，基本独立于结合位点的相似性（扩展数据图4k）。
seq2PRINT归因分数识别出预测足迹的关键DNA序列模式，使得能够从头识别基序。使用训练好的模型，我们聚类和对齐局部序列归因分数，并识别出106个从头基序。这些从头基序以无偏的方式恢复了已知基序，以及如SOX二聚体等复合基序（图2h）。一些从头基序与已知基序数据库没有显著匹配，但与强转录因子或核小体足迹相关（图2h，底部）。我们提供了本研究中发现的所有从头基序及其相关的多尺度足迹的完整列表（补充数据1和2）。
多尺度足迹分析和seq2PRINT解析了造血过程中转录因子结合的动态变化。我们使用同时高通量ATAC和RNA表达测序（SHARE-seq）生成了来自七名人类供体的874,480个骨髓单核细胞的联合scATAC-seq和scRNA-seq数据集（图3a和扩展数据图5a-g）。为了实现足迹分析，我们将单细胞合并成1,000个伪批量，代表所有主要细胞类型和发展过渡（扩展数据图5h-m和方法）。将深度学习序列模型应用于所有数据表示22亿次读取，并使用大型模型的低秩适应（LoRA）对每个伪批量进行微调，实现了与单独训练每个伪批量相比约80倍的速度提升和预测准确性提高（扩展数据图6a-c）。
转录因子结合预测显示同一CRE在不同细胞类型中结合不同组的转录因子。在髓系主调节因子SPI1（PU.1）启动子处，seq2PRINT结合评分在髓系细胞中对SPI1和AP-1位点较高，而在红细胞中仅预测强烈的GATA1结合，这与这些转录因子已知的调控关系一致（图3b）。值得注意的是，这些不同细胞类型的转录因子结合模式和其他位点上的模式不能仅通过测量启动子的整体可及性来区分（扩展数据图6d-g）。我们使用主成分分析量化了每个CRE处转录因子结合模式的复杂性，发现复杂的CRE（多于一个成分）在远端CRE中相对于启动子高度富集（69.8%对17.2%；图3c），突显了增强子的细胞类型特异性利用。
对红细胞分化轨迹的转录因子结合分析表明，CRE的建立是逐步进行的。这在血红蛋白位点控制区的HS3增强子中尤为明显。在造血干细胞和共同髓系祖细胞中，核小体无规则排列，低转录因子结合被预测。随着细胞沿红细胞发育的进展，按伪时间排序，β-珠蛋白（HBB）的表达增加，核小体足迹首先在CRE边缘出现，并在GATA/TAL基序位点有强转录因子结合评分。然后CRE逐渐扩展，在边缘添加了KLF1/NFE2因子结合（图3d和扩展数据图6h）。HBB启动子也表现出相同的顺序转录因子结合模式，且在该位点的PRINT/seq2PRINT结合预测与先前发布的大量平行报告分析数据非常吻合（扩展数据图6i-m）。
我们发现这种模式在整个基因组中围绕中央转录因子向外扩展的红细胞相关CRE中普遍存在。GATA和TAL因子的转录因子结合评分在红细胞伪时间早期增加，而整体CRE开放和NFE2、KLF1、NR2F1和AP-1因子的结合则在分化后期发生（图3e和扩展数据图6n）。有趣的是，在可及性获得之前结合并靶向基因表达的转录因子也更靠近CRE中心，而后期结合的转录因子基序位于侧翼区域（图3f,g和补充表3）。这些观察结果主要独立于CRE的复杂性（扩展数据图6o-r）。这一观察结果与之前的研究一致，表明CRE内转录因子基序的刻板排列和染色质重塑剂在CRE边缘的富集。
接下来，我们使用seq2PRINT研究了衰老过程中CRE组织的变化。生物学衰老是一个多因素过程，影响广泛的组织生理，包括造血干细胞（HSC）的功能和增殖显著变化。由于之前的研究表明，衰老伴随着广泛的表观遗传变化，我们认为seq2PRINT非常适合检测衰老过程中转录因子活性和CRE结构的差异。我们通过荧光激活细胞分选（FACS）从年轻（11周龄，n=10）和老年（24个月龄，n=5）雄性小鼠的骨髓中分离出造血祖细胞（Lin−）和HSC（Lin−Sca1+cKit+CD48−CD150+），并使用10x Multiome平台获得了涵盖14,640个HSC和33,585个造血祖细胞的48,225个细胞的联合ATAC-RNA谱型（图4a,b和扩展数据图7a-g）。年轻和老年小鼠的HSC分别聚集在一起，并与之前报道的衰老特征和标记基因相匹配（图4c,d，扩展数据图8a-e和补充表4）。为了进一步探索HSC异质性，我们收集了先前描述的基因集，涵盖了谱系偏向和衰老，以及从我们的数据中学习到的基因程序（方法）。通过对这些基因集进行评分和聚类，我们识别出了五个HSC亚群，这些亚群按年龄和谱系潜力区分（图4e,f和扩展数据图8f-m）。这五个HSC细胞状态代表了与先前发现一致的主要类别，并且可以进一步分区，表明额外的衰老相关变异轴（如氧化磷酸化调控、未折叠蛋白反应）可能与衰老免疫系统病理学有关（补充表5）。
我们发现，随着衰老，CRE发生了广泛的转录因子结合变化。我们对所有HSC伪批量应用了seq2PRINT，并比较了年轻和老年HSC，发现NF-I、Runx、Ets、Gata和AP-1（例如Fos、Jun）的活性显著增加（图4g和补充表6）。考虑到比较所有年轻和老年HSC可能会受到衰老过程中HSC亚群组成的混淆，我们分别对巨核细胞（Mk）偏向和多线祖细胞进行了亚群特异性比较。亚群特异性比较结果总体相似，表明大多数变化在亚群间共享（扩展数据图9a-d）。然而，某些Ets基序，特别是Spi1/Spib/Spic和Elf，在多线祖细胞中表现出与衰老相关的下调，表明衰老过程中亚群特异性转录因子变化（扩展数据图9c,d）。
为了全面、无偏地检查seq2PRINT学到的序列基序，我们识别了差异活动的从头基序（图4g，扩展数据图9e和补充表6），包括许多在衰老过程中丢失的富含CG的从头基序，可能与DNA结合因子对甲基化变化的识别有关（扩展数据图9f）。我们发现，含有Ets同源或异源二聚体的复合基序增加（图4g-h和扩展数据图9g-h）。此外，我们发现，在衰老过程中，Ets/Runx复合基序在老年多线祖细胞中特别富集（扩展数据图9h,i）。支持这一发现，之前的工作提出了Ets/Runx共现对HSC维持和髓系命运的作用。机制上，Ets与DNA的结合受到自抑制结构域的负调控，该结构域在与Runx或第二个Ets家族转录因子共结合后释放，并且在两个实验确定的结构中显示了物理相互作用（PDB-4L0Z和PDB-2NNY；图4h）。为了测试seq2PRINT预测的Ets复合基序是否代表类似的直接相互作用，我们使用AlphaFold3预测结构，发现所有seq2PRINT识别的Runx/Ets构型无论方向如何都显示出类似的物理相互作用。支持这些预测的有效性，基于已知Ets/Ets和Ets/Runx基序构型的AlphaFold3结构与实验测量结构高度一致（预测结构10和PDB-4L0Z之间的均方根偏差为0.825 Å），并且溶剂不可接近的碱基与每个转录因子的核心基序匹配（图4h，扩展数据图9j，补充说明和补充表7）。seq2PRINT因此揭示了衰老过程中转录因子的重新布线，并通过从头基序分析和AlphaFold3预测了涉及衰老和HSC自我更新的广泛Ets和Runx共结合构型。我们还提供了跨伪批量的TF评分，以使探索与基因程序相关的TF成为可能（补充表5）。
在多个系统中，各种研究表明，衰老和细胞衰老伴随有核小体的全局丢失，尽管关于这种丢失是全局性的还是局限于特定的TF相关位点仍存在争议。为了进一步探讨HSC的表观遗传衰退，我们使用PRINT测量了衰老过程中CRE内的核小体占有率。我们观察到，在老年HSC中，CRE内的核小体足迹广泛减少（图4i和扩展数据图9k-q）。我们注意到，38%的CRE具有差异核小体足迹，但没有显著的可及性变化，突显了足迹分析在解析变化方面的效用。核小体足迹的减少并不一定意味着基因组中核小体数量的减少，因为这也可能是由于核小体定位/相位的丧失，且分析主要限于CRE。为了确定这些位点上是否有特定的转录因子发生变化，我们使用seq2PRINT基序并发现下调的TF包括Yy1、Nrf1和Ctcf，而上调的TF包括Nfyb、Atf、Jun和Klf（图4j,k 和 扩展数据图9m,n,p,q 和 10a-f）。值得注意的是，Yy1、Nrf1和Ctcf也被预测在多种细胞系中调控核小体位置（图2c,g）。Yy1已被证明是HSC自我更新和静止的调控因子。总体而言，我们发现衰老过程中核小体足迹广泛改变，转录因子结合重新布线，提示衰老过程中CRE结构失调（图4l）。
讨论
我们的研究结果展示了细胞分化和衰老过程中CRE处转录因子结合和核小体重定位的复杂动态变化。之前的足迹研究表明，转录因子结合主要由CRE的整体开放或关闭决定，而不是CRE内转录因子的差异结合。相比之下，我们的研究发现不同细胞类型中的CRE被不同组的转录因子占据。这在SPI1/PU.1和Wasl启动子分析中得到了体现（图3b和4k），尽管整体可及性水平相似，但在多个核小体和转录因子配置中观察到，显示出标准染色质可及性分析会遗漏的基因调控机制。直接支持这一模型的研究表明，转录因子在发育过程中会发生切换。在连续的造血分化轨迹中，我们发现CREs围绕中心转录因子逐步扩展，侧翼转录因子在发育的后期阶段结合，表明增强子的建立是一个类比而非数字的过程。
更广泛地说，PRINT生成了给定细胞群体中所有DNA结合蛋白的同时图像。使用seq2PRINT建模足迹不仅可以解释转录因子结合，还可以从头识别转录因子基序和协同结合，表明转录因子具有专门的功能，如重塑或稳定邻近核小体。这些属性使seq2PRINT区别于ChromBPNet，后者旨在预测可及性谱型而非解释支撑足迹的序列特征。使用LoRA，我们减少了足迹预测的计算负担，使seq2PRINT能够扩展到单细胞表观基因组学数据。我们预计这种方法将使高分辨率足迹分析与基因组结构和局部基因表达等多种表观基因组学数据类型相结合。同样，以碱基对分辨率识别转录因子结合并将足迹归因于特定序列，也可以为之前因峰基分析而被掩盖的致病遗传变异赋予新功能。
这些方法存在局限性。我们的方法不包括转座片段长度，这可能是分析的一个附加特征。此外，由于深度学习需要多个示例来识别模式，罕见或非典型结合事件可能会被忽略。此外，单细胞ATAC-seq的足迹分析需要伪批量处理以产生足够的Tn5插入事件来识别足迹。我们设想，PRINT可以通过其他方法加以补充，如基于甲基转移酶的单分子足迹分析和DNA序列突变分析，以进一步分析特定CRE的结构和功能。然而，由于ATAC-seq已被广泛采用和广泛用于单细胞分析，seq2PRINT可能独特地使人们能够回顾性和前瞻性地研究健康和患病人类组织中广泛的CRE足迹图谱。综上所述，我们的方法从未维染色质可及性数据中提取了丰富的多维特征空间，揭示了CRE在高基因组和细胞状态分辨率下的动态结构。
(全文结束)
