多组 | Nature | 多组学与深度学习解析人类发育中的调控语法

Basic Information

• 英文标题：Multiomics and deep learning dissect regulatory syntax in human development
• 中文标题：多组学与深度学习解析人类发育中的调控语法
• 发表日期：08 April 2026
• 文章类型：Article
• 所属期刊：Nature
• 文章作者：Betty B. Liu | William J. Greenleaf
• 文章链接：https://www.nature.com/articles/s41586-026-10326-9

Abstract

Para_01

1. 转录因子通过以序列特异性方式结合调控DNA来在发育过程中建立细胞身份，通常促进局部染色质可及性并调控基因表达。
2. 绘制可及性染色质图谱可为转录调控提供关键见解，但目前可用于人类发育研究的数据集仅限于整体组织、单一器官或单一模态。
3. 本文介绍了人类发育多组学图谱，这是一个涵盖12个器官、203种细胞类型、817740个胎儿细胞的单细胞染色质可及性和基因表达图谱，包含超过100万个候选顺式调控元件，其中许多表现出器官特异性的体内增强子活性。
4. 通过训练深度学习模型，利用局部DNA序列预测染色质可及性，揭示了一整套影响可及性的基序语言，包括具有不同句法约束的复合基序，这些基序被预测可介导转录因子协同作用。
5. 我们识别出要求精确基序间距和方向的‘硬性’句法规则、允许灵活基序排列的‘软性’规则，以及普遍存在的抑制可及性的基序。
6. 基于模型对遗传变异的解读表明，具有正负效应的基序遭到破坏时，会与基因表达的变化产生一致的影响。
7. 我们的研究阐明了基序句法如何调控细胞类型特异性的染色质可及性，并为解析人类发育过程中顺式调控逻辑和解读遗传变异提供了基础性资源。

Inline Recommendations

Similar content being viewed by others

其他人正在查看的相似内容

1. 小鼠器官发生过程中染色质可及性的单细胞图谱
2. 成年小鼠大脑的单细胞DNA甲基组和三维多组学图谱
3. 染色质可及性图谱定义了跨组织的基质细胞身份

Main

Para_01

1. 在人类发育过程中，细胞类型的多样性通过转录因子的差异表达和活性产生，这些转录因子整合细胞内在和外在信号，指导基因调控。
2. 转录因子结合顺式调控元件中的特定DNA序列，通常诱导局部染色质可及性，并改变邻近基因的表达。
3. 然而，我们尚缺乏对驱动人类发育过程中染色质状态变化的转录因子模体的全面认识，这限制了我们理解转录因子结合位点的组织结构——即语法——如何参与调控。

Para_02

1. 利用DNase I超敏感位点测序（DNase-seq）和基于测序的转座酶可及性染色质检测（ATAC-seq）来绘制染色质可及性图谱，已经能够通过人类组织中的序列模体推断转录因子活性。
2. 然而，群体测量会掩盖细胞异质性，而大多数单细胞图谱主要集中于单一器官或单一组学模式。
3. 需要一种跨器官、多模态的视角，以捕捉顺式调控的细胞环境特异性，并将染色质状态与转录程序联系起来。

Para_03

1. 染色质可及性通常源于转录因子结合的协同作用，这种协同作用要么通过转录因子之间的直接相互作用实现，要么通过与核小体的竞争实现。
2. 这些机制分别对基序语法施加了硬（固定）或软（灵活）的限制。
3. 然而，这些规则在人类发育过程中的普遍性在很大程度上尚不清楚。
4. 此外，与复杂疾病相关的遗传变异富集于非编码基因组区域，但我们预测这些变异在特定细胞类型中破坏调控活性的能力仍然有限。

Para_04

1. 最近的深度学习模型通过局部DNA序列预测碱基分辨率的染色质可及性谱型，这些模型能够学习影响可及性的因果序列特征。
2. 除了从头发现具有预测性的基序和转录因子足迹外，这些模型还能通过预测DNA序列变化对可及性的定量影响，在计算机中探究调控序列的语法结构以及非编码遗传变异。
3. 因此，深度学习模型为解析转录因子结合如何影响染色质可及性的规律，以及将序列变异与顺式调控的破坏联系起来提供了一个强大的框架。

Para_05

1. 本文介绍了人类发育多组学图谱（HDMA），这是一个多组学、多器官的单细胞图谱，涵盖了12个人类胎儿器官中的染色质可及性和基因表达谱。
2. 我们绘制了超过一百万个可及的调控元件，并证实了它们在体内解析器官特异性和细胞类型特异性增强子活性的能力。
3. 我们训练并解析了深度学习模型，以预测细胞类型分辨的染色质可及性，定义了驱动可及性的调控序列基序词汇表，并推断出每种细胞类型在整个基因组中的预测性基序实例。
4. 对基序语法的分析揭示了严格的和宽松的语法规则。
5. 最后，我们优先筛选了可能在发育过程中干扰调控功能的疾病相关变异。
6. HDMA为解码顺式调控语法、将序列变异与基因调控联系起来，以及理解DNA序列如何影响人类发育过程中的转录调控提供了基础资源。

A multiomic human development atlas

Para_01

1. 我们利用SHARE-seq21这一可扩展的分裂-混合组合条形码平台，对受孕后第10至23周跨越12个器官的817,740个胎儿细胞同时进行了染色质可及性和基因表达的分析（图1a，扩展数据图1和图2a以及补充说明1）。
2. 所获得的高维度多组学数据（HDMA）来自相同细胞的真实多组学测量结果，在转录起始位点（TSS）富集程度、唯一分子标识符（UMI）数量以及检测到的基因数量方面，相比之前发布的多器官胎儿‘单组学’图谱均有提升（图1b，扩展数据图2b、c以及补充表1）。

Fig. 1: A multiomic, multi-organ atlas of human development.

- 图片说明

◉ a，HDMA数据集概览及按器官划分的样本量分布。◉ b，RNA测序和ATAC-seq技术的样本元数据和质量控制指标，样本来自受孕后第10至23周。每个器官的生物学独立样本数量如a所示。PCD，受孕后天数。◉ c，203种层次聚类细胞类型的标志基因表达和转录因子基序chromVAR偏差z分数。TF，转录因子。◉ d，内皮细胞不同层次细胞注释的示例。AG，肾上腺；BR，脑；EY，眼；HT，心脏；LI，肝脏；LU，肺；MU，肌肉；SK，皮肤；SP，脾脏；ST，胃/食管。

Para_02

1. 我们采用迭代方法，结合经典的区室标记、已知的细胞类型标记以及来自单细胞转录组的从头标记基因，注释了203个不同的细胞簇（L1簇）（方法）。
2. 通过在簇特异性可及染色质峰内显著富集的转录因子基序，进一步验证了细胞身份（补充说明2）。

Para_03

1. 对顶级标记基因的伪批量表达谱进行层次聚类显示，不同器官共有的细胞类型（如内皮细胞、成纤维细胞和免疫细胞）通常聚集在一起（图1c）。
2. 根据共有标记基因的表达情况，我们手动将相关的L1聚类归为134个合并后的细胞聚类（L2聚类），并进一步归为37个广泛的细胞类别（L3聚类）（图1d以及补充表2和3）。

Para_04

1. 通过全局 chromVAR22 偏差分析，鉴定了在同一细胞类型的不同 L1 簇间共有的转录因子基序，这些基序的可及性增加
2. 与标记基因相比，转录因子基序的可及性变化通常特异性较低，许多基序在不同细胞类型中表现出多效性的可及性模式
3. 例如，尽管 CUX2 标记基因的表达仅限于神经元，但在神经元中排名最高的 NEUROD1 基序却在神经元、成纤维细胞以及神经嵴来源的细胞类型中均显示出可及性增加，这反映了转录因子在发育过程中存在的基序简并性和依赖于环境的活性

Accessibility landscape of development

Para_01

1. 我们通过整合203种L1细胞类型中独立识别的峰，对人类发育过程中可及的调控图谱进行了编目。通过峰的迭代重叠，定义了一组全局性的1,032,273个染色质可及的顺式调控元件（caCREs），每个元件跨度为500个碱基对， collectively覆盖了约17%的基因组。大多数caCRE的最高标准化ATAC-seq信号出现在肝脏、眼睛、心脏以及胃和食管（胃/食管）细胞类型中（扩展数据图2d）。

Para_02

1. 与ENCODE v4数据库的比较显示，85.1%的HDMA caCREs与先前定义的候选顺式调控元件（cCREs）重叠（扩展数据图2e）。
2. 值得注意的是，我们检测到了所有ENCODE v4位点中的48.7%，其中包括56.2%的CTCF结合元件。
3. 在那14.9%未与ENCODE cCREs重叠的caCREs中，脑组织和眼部来源的元件所占比例较高，这表明单细胞分析能够揭示出在ENCODE基于群体细胞数据集中缺失的细胞类型特异性调控元件（扩展数据图2d）。

Para_03

1. 我们使用活性-接触（ABC）模型将caCREs与基因进行关联，该模型应用于每个L1细胞簇内。
2. 对于每个广泛的L3簇，我们汇总了其组成L1簇的ABC关联，并进行筛选，保留靶向具有可及启动子基因的元件。
3. 在每个L3簇中，我们将高度关联基因（HLGs）定义为ABC关联caCRE数量最多的前1%基因（补充表4）。
4. 基因本体（GO）分析显示，HLGs富集于细胞类型特异性过程，例如免疫细胞中的‘B细胞增殖’、角质形成细胞中的‘角质形成细胞分化’以及色素上皮细胞中的‘相机型眼形态发生’（图2a）。
5. 广泛的转录调控相关术语，包括‘RNA聚合酶II介导的转录正/负调控’，在37个L3簇中的25个中均显著富集。
6. 其他广泛富集的发育相关术语，如‘胚胎骨骼系统形态发生’和‘前后模式特化’，是由已知调控发育程序的HOX基因驱动的。
7. GO分子功能分析证实，在37个L3簇中的33个中，HLGs在转录因子结合活性方面显著富集（扩展数据图2f）。
8. 例如，内皮细胞HLGs包括经典调控因子（NR2F2、ELF4、SOX18、NFIC和FOXC1）、共调控因子（BCL9L），以及与内皮功能相关的基因（EGFL7、PLXND1和PLEC）。
9. 内皮HLGs中富集的顶级GO术语包括通用转录因子活性术语以及内皮特异性术语，如PDGFR结合（图2b,c）。

Fig. 2: Connecting DNA regulatory elements to genes and resolving organ and cell-type specificity in experimentally validated enhancers.

- 图片说明

◉ a，每个广泛类别细胞簇（L3细胞簇）中HLG富集的前5个GO生物学过程术语。Pos代表阳性；neg代表阴性；reg代表调控。◉ b,d，具有最多ABC关联caCRE的内皮细胞（b）和全局（d）HLG。◉ c,e，内皮细胞HLG（c）和全局HLG（e）的GO分子功能术语富集情况，统计显著性由单侧Fisher精确检验确定。◉ f，在与VISTA增强子重叠的HDMA caCRE中标准化并进行z分数处理的ATAC-seq信号。◉ g，VISTA胚胎mm101切片的明场（上图）和H&E染色图像（下图）。蓝色来自X-Gal染色，指示增强子活跃的位置。该检测基于一只携带该增强子的完整小鼠胚胎进行切片和组织学分析。在两个切片上重复成像，结果相似。A代表前部；D代表背侧；P代表后部；V代表腹侧。◉ h，在肝脏细胞类型中VISTA增强子mm101的可及性及其与ABC关联的α-珠蛋白基因（HBA2）的表达。RPKM表示每百万比对读段每千碱基的读段数。

Para_04

1. 为了识别在不同细胞类型中由caCREs保守调控的基因，我们对caCREs进行了筛选，以获得在所有L3聚类中 consistently 被识别出的caCRE–基因连接。
2. 这些全局高连通性基因（gHLGs），对应于具有最稳定ABC连接caCREs的前1%基因，也在转录因子活性方面显著富集（图2d,e），并包含普遍存在的调控因子（KLF9、KLF4和BCL6）。
3. 这81个gHLGs构成了一组在人类发育过程中普遍存在且高度受调控的转录调控因子集合，表明转录因子是发育过程中最受调控的基因之一，这与我们之前在发育中的人脑中的观察结果一致（补充表4）。

Resolving regulatory element specificity

Para_01

1. 为了验证HDMA-derived增强子，我们检测了与VISTA数据库中通过小鼠胚胎第11.5天X-gal报告染色实验验证的增强子重叠的caCREs。
2. 在具有hg38坐标的3,272个VISTA增强子中，有851个与HDMA中检测的器官重叠，其中包括与1,032,273个caCREs相交的766个。
3. 对于注释为在脑、心脏和眼睛中具有活性的VISTA增强子，我们在对应器官的HDMA细胞类型簇中观察到可及性的显著富集（单侧Wilcoxon秩和检验，脑P = 10⁻⁴³⁷，心脏P = 10⁻⁴⁵²，眼睛P = 9.11 × 10⁻⁴⁸；AUROC概率分别为脑0.72、心脏0.75、眼睛0.59）（图2f），以及基因表达的显著富集（单侧Wilcoxon秩和检验，脑P = 1.04 × 10⁻³⁸，心脏P = 1.23 × 10⁻¹⁶⁶，眼睛P = 2.77 × 10⁻²³；AUROC概率分别为脑0.57、心脏0.67、眼睛0.57）（扩展数据图3a）。

Para_02

1. 令人惊讶的是，我们数据中几个先前被注释为心脏特异性的增强子仅在肝细胞类型中表现出强烈的可及性（图2f，虚线黑框）。
2. 鉴于小鼠肝脏天然的颜色及其位于心脏下方的位置，这些情况可能在此前仅通过视觉检查时被遗漏。
3. 专家对我们推定在肝脏中活跃的顶级候选增强子进行评审后，确认其中六个增强子（mm2143、mm69、mm291、mm101、mm257和mm18）在肝脏和心脏中均具有活性。
4. 对X-gal染色胚胎进行切片和成像验证了这六个候选增强子在肝脏中的活性（图2g和扩展数据图3b–f）。

Para_03

1. 一个候选元件mm101与已知的α-珠蛋白超级增强子位点重叠，并通过ABC分析与α-珠蛋白基因簇相连（图2h）。
2. α-珠蛋白是胎儿血红蛋白的一个亚基，仅存在于成红细胞中。
3. 由于胎儿肝脏在发育过程中是造血的场所，我们推测mm101增强子的活性将局限于成红细胞。
4. 事实上，mm101位点的染色质可及性以及HBA2基因的表达均特异性地在肝成红细胞和增殖中的成红细胞中升高（图2h）。
5. 组织学分析证实，在肝窦状隙内存在具有球形形态以及苏木精和伊红（H&E）染色所显示的特征性深紫色的肝成红细胞，这些细胞表现出X-gal染色阳性（图2g）。

The transcription factor motif lexicon

Para_01

1. 为了识别能够预测染色质可及性的顺式调控序列特征，我们训练了深度卷积神经网络（ChromBPNet），利用局部DNA序列来建模ATAC-seq信号的形态和强度。
2. 对于每种细胞类型，我们训练模型以预测在峰区和背景区域1000-bp窗口内测序读段的总数及其碱基分辨率分布，输入为2114 bp的局部序列上下文信息（图3a及方法部分）。
3. 我们对203种细胞类型分别训练了五折交叉验证模型，并保留了其中189种细胞类型的模型（总计n = 945个模型），这些模型在测试染色体的保留峰区中，预测值与实际观察到的对数读段数之间的中位皮尔逊相关系数达到0.78（图3b、c，扩展数据图4以及扩展数据图5a、b，补充表5）。

Fig. 3: Repertoire of transcription factor motifs in human development.

- 图片说明

◉ a，工作流程概述。对于每种细胞类型，训练一个 ChromBPNet 模型，以从局部DNA序列预测染色质可及性。通过解释模型得到代表对可及性贡献的单碱基评分。高贡献区域在细胞类型内部及之间聚类形成基序，这些基序用于注释每种细胞类型中的高贡献区域。CWM 表示贡献权重矩阵。◉ b，在心肌细胞中 SRF 基因座的示例，显示了观察到的和预测的染色质可及性、推断的核小体占据情况、单碱基贡献评分以及注释的基序实例。◉ c，上图：每个模型在峰值区域内预测值与观测值之间的皮尔逊相关系数分布（五个折叠的平均值）。下图：某一细胞类型、一个模型折叠在峰值区域内的预测值与观测值的对数计数对比。◉ d，独特的新发基序的总结与分类，以及来自部分类别的基序示例。◉ e，基础基序的汇总，每一行对应一大类基础基序。从左到右依次为：以CWM形式表示的新发基序、最相似的已知基序（位置权重矩阵PWM）、对可及性的贡献方向、跨细胞类型的基因组实例总数、实例重叠各种基因组特征的比例、各实例（按细胞类型）到最近转录起始位点TSS的中位距离分布、来自各组织隔室中细胞类型的实例比例、来自各器官中细胞类型的实例比例，以及每个大类基序内的基序变体数量。对于每个大类基序，括号中标注的是检测到匹配信号的细胞类型数量，并对应于每个箱线图分布中包含的细胞类型数。箭头指示文中提到的基序：普遍存在且以启动子为主的基序（深灰色）、普遍存在但位于远端的基序（浅灰色）以及组织特异性的远端基序（黑色）。◉ CWM 表示贡献权重矩阵。

Para_02

1. 我们研究了哪些顺式调控序列特征会影响染色质的可及性。我们预计这些序列对应于转录因子结合位点，这些因子会保护其遮蔽的核苷酸免受Tn5切割，但会促进其周围区域的局部可及性。
2. 我们使用DeepLIFT方法在每种细胞类型中分别估算了每个核苷酸对可及性的预测重要性，通过TF-MoDISco将特征聚合成基序，并在不同细胞类型间对基序进行聚类，从而构建出一个包含508个从头基序的词典，这些基序被预测可在189种细胞类型中调控染色质的可及性（图3a,d,e及方法部分）。

Para_03

1. 我们根据可能结合这些从头发现的基序的转录因子对其进行注释，通常在转录因子家族水平上进行注释，因为同一家族内的基序具有相似性（方法和补充表6）。
2. 一些基序与现有数据库中的已知转录因子结合基序（基础基序）或带有额外预测侧翼核苷酸的已知基序（带侧翼的基础基序）相匹配（图3d）。
3. 该基序图谱还识别出经典基序的变体，包括核心CTCF基序和上游CTCF基序，反映了不同CTCF锌指结构域与DNA的相互作用（扩展数据图5c）。
4. 超过一半的从头发现的基序为复合基序，由两个同类型位点（同源复合体）或异类型位点（异源复合体）组成。
5. 其中一组包含38个“未解析”基序，它们与已知数据库中的基序无显著相似性。
6. 大多数基序（n = 493，占97%）具有正向贡献分数，表明它们促进染色质可及性，而一小部分基序具有负向分数，提示其降低可及性（图3e和补充表6）。
7. 基于它们对可及性的方向性影响，我们将这些基序分别称为‘正向’和‘负向’基序。

Para_04

1. 接下来，我们使用Fi-NeMo在每种细胞类型可及的峰区中识别出每个基序的预测性实例（图3a，方法）。
2. 平均每种细胞类型包含839,544个基序实例，每个峰区含有2至8个实例，其中65%的实例位于峰顶点150 bp范围内（扩展数据图5d–j）。
3. 基序实例的数量随着每种细胞类型峰值总数的增加而增长，表明更深度的测序可能检测到更多的基序实例（扩展数据图5k）。

Para_05

1. 例如，我们重新研究了在胎肝成红血细胞中具有活性的mm101增强子。模型解析显示该增强子内存在两个得分较高的GATA基序实例（扩展数据图5l），这与GATA1在红细胞生成中的关键作用一致，并支持其在此背景下对α-珠蛋白基因激活的贡献。

Para_06

1. 接下来，我们将预测的基序实例按基因组背景和细胞类型特异性进行分层，并计算其到最近转录起始位点（TSS）以及通过NucleoATAC推断出的核小体dyad的距离（图3e）。
2. 该分析揭示了一组普遍存在、以启动子为主、靠近TSS的基序，包括NRF1、NFY、YY1或YY2（YY1/2）、TBP、SP或KLF（SP/KLF）、ETS和BANP。
3. 这些基序通常富含CG，与启动子区域富含CpG的特征一致。
4. 已知NRF1和BANP是DNA甲基化敏感的转录因子，而其他一些基序已被证实参与转录起始。
5. 相比之下，器官和细胞类型特异性的基序主要位于远端或内含子区域。
6. 这些基序包括SPI和RUNX（免疫细胞）、POU（眼和脑）以及HNF4（肝）（图3e）。
7. 在内皮细胞等出现在多个器官中的细胞类型中，独立训练的模型学习到了一致的序列贡献（扩展数据图6a），表明调控词汇的发现具有鲁棒性。
8. 我们还鉴定出两个在胃、肺和肝中有限使用的未解析回文元件（扩展数据图6b）。

Para_07

1. 为了理解上下文依赖的组合基序逻辑，我们在每种细胞类型中测试了所有基础基序对，以检测其基序实例是否存在显著的共现情况（扩展数据图7）
2. CTCF与其他基序的共现情况较为有限，而以启动子为主的基序则在不同组织中普遍存在相互共富集现象
3. 我们发现了若干依赖组织的基序配对关系，例如在眼部细胞类型中，ONECUT基序与HD和PAX基序共同富集，而在肝脏中则与FOX以及HNF1A或HNF1B基序共同富集

Inferring distinct modes of transcription factor synergy

Para_01

1. 调控元件上转录因子的组合结合可通过协同相互作用增强特异性和功能。
2. 在DNA介导的协同作用中，具有直接蛋白质-蛋白质相互作用的转录因子结合特定的复合识别位点，并呈固定排列；或者，在此类位点上的DNA可稳定转录因子之间较弱的相互作用。
3. 相比之下，核小体介导的协同作用源于转录因子与核小体之间的主动或被动竞争。
4. 我们推断，DNA介导的协同作用在蛋白质-蛋白质相互作用的长度尺度内（<20 bp）对结合位点施加了固定的间距和方向限制（硬语法），而核小体介导的协同作用则兼容结合位点的灵活组织（软语法）以及更长距离的位点间隔（20–150 bp）。

Para_02

1. 为了系统地识别每个从头鉴定出的复合基序的协同转录因子相互作用及其句法约束，我们使用 ChromBPNet 模型和 tangermeme 软件包实现了一种体外边际化分析框架，该方法全面评估了两个组成基序在不同间距或方向排列下对染色质可及性的联合效应，相较于它们的独立效应（图4a、b及方法部分）
2. 我们在每种细胞类型中重复了这一分析，但将重点放在表现出最多可预测复合基序实例的细胞类型上（方法部分、补充说明3以及补充表7）
3. 我们将协同效应定义为显著偏离独立基序效应的对数相加模型，并进一步根据基序间距和效应强度，将表现出协同作用的复合基序分类为具有严格或宽松句法特征，即其联合效应超过相加预期的情况（方法部分）

Fig. 4: Systematic inference of transcription factor cooperativity.

- 图片说明

◉ a，一个或两个不同基序的独特取向。◉ b，在计算机上进行边缘化的流程。◉ c–g，对包含HD位点和BHLH位点（E-box）的复合基序进行协同作用分析。◉ c，当基序以0–9 bp间距插入时，针对一种排列方式的预测谱图。谱图显示为五个fold的均值±标准差。◉ d，复合基序的CWM（卷积权重矩阵）。◉ e，在5 bp间距下的预测谱图，以及对插入基序的编辑序列进行模型解释后的平均贡献分数。高亮部分表示插入的核苷酸。◉ f，以对数计数形式量化的两个基序在可及性上的联合效应，按间距和方向分别展示。每个点代表一个fold，条形表示均值。◉ g，在n = 100条序列上的预测效应。P值通过单侧Wilcoxon符号秩检验获得。◉ h，推断出具有严格语法偏好的所有复合基序中，基序之间的距离（上图）以及基序中心之间的距离（下图）。◉ i，所有测试的基序对在其最优排列下的预测联合效应，与其独立效应之和的比较。◉ j，每种结果语法类别中的复合基序数量。◉ k，三种代表性复合基序的预测联合效应与独立效应之和的比较。效应曲线显示为五个fold的均值±标准差。◉ l，IKZF–RUNX复合基序在每种细胞类型中通过计算机边缘化得到的预测平均可及性谱图（左），平均效应（中），以及可能结合各组成基序的转录因子在每种细胞类型中的表达情况（右）。误差线表示模型不同fold之间的标准差。展示了预测可及性最高的40种细胞类型。◉ m，与l相同，但针对FOX–HNF4复合基序。◉ NK，自然杀伤细胞；TEC，胸腺上皮细胞。

Para_03

1. 我们鉴定出多个已知的受语法约束的复合元件。
2. 一个与近期报道的Coordinator元件相匹配的从头基序——由E-box和TAAT同源异形盒基序组成的复合体——在5 bp头尾相连的排列方式下表现出强烈的协同效应（扩展数据图8a–e）。
3. 值得注意的是，这种排列方式精确匹配了经过实验验证的Coordinator复合物的结合几何结构，其中TWIST1–TCF4异源二聚体与ALX4协同结合。
4. X射线晶体学研究表明，这种严格的基序间距对于稳定TWIST1骨架与ALX4侧链之间的相互作用是必需的。
5. 我们主要在皮肤细胞类型中检测到该基序，这与其在神经嵴来源的间充质中的功能一致。
6. 我们还观察到另一种由E-box和一种替代性同源异形盒位点（TGATTGAT）组成的复合基序表现出类似的受语法约束的协同效应，该基序主要存在于肌肉和胸腺细胞类型中（图4c–g）。

Para_04

1. 在测试的所有138个从头组合基序中，我们鉴定出67个具有显著协同作用的基序（图4h–j及补充表7）。
2. 硬性语法基序在基序中心间距约为7 bp和11 bp时表现出偏好，这与蛋白质结合在相邻或螺旋转角偏移的DNA沟槽中的情况一致，并在特定排列下显示出预测联合效应的显著峰值（图4h,k，扩展数据图8f,g）。
3. 相比之下，软性语法基序表现出更宽泛的距离偏好，其联合效应较温和，并在距离超过约100 bp后逐渐衰减至独立效应（扩展数据图8f,g）。
4. 在硬性语法基序中，我们发现了与已知二聚体转录因子相互作用一致的复合元件，包括SOX同源二聚体、p53同源二聚体以及GATA–TAL异源二聚体（补充说明3）。
5. 值得注意的是，预测的p53同源复合基序的最佳间距和方向与其两个NCATGN结合位点以头尾相接方式排列时经典的4 bp间隔相匹配。
6. 我们重点列出了一些已有报道显示存在协同作用但对其间距和方向限制此前研究较少的基序对，包括RUNX–RUNX、FOX–NR和ETS–NR，这些可能代表人类胎儿发育中的新语法特征。
7. 在所有细胞类型中对所有基序对进行筛选提供了隐式的阴性对照；我们的策略能够从绝大多数无协同作用的基序组合中区分出一小部分（8.9%）主要为细胞类型特异性的协同基序对（补充说明3）。

Para_05

1. 在大多数细胞谱系中检测到了协同基序语法。大约60%的细胞类型显示出至少一种硬语法和一种软语法的复合基序模式。在每种细胞类型中，具有硬语法和软语法的复合基序相较于无协同作用的基序和非复合基序表现出更高的总贡献分数，且硬语法基序尤为普遍，表明其对染色质可及性具有更强的预测影响力

Para_06

1. 接下来，我们通过将成对的基序以最佳排列方式插入背景序列，并预测所有细胞类型中的可及性，来检测协同效应是否具有细胞类型特异性（补充说明3）。
2. 例如，IKZF–RUNX复合基序仅在部分来源于胸腺和脾脏的免疫细胞类型中驱动了可及性（图4l）。
3. 这一模式与这些相同细胞中IKZF和RUNX家族转录因子的表达一致，表明协同效应也受到其组成因子可用性的限制。
4. 类似地，FOX–HNF4复合基序仅在肝细胞中表现出预测的可及性，这与这些细胞中特定FOX和HNF4家族成员的受限表达相关（图4m）。
5. 尽管由于单细胞RNA测序对低水平表达的转录因子基因检测能力有限，难以进行全面分析，但这些实例表明，基序协同作用共同依赖于精确的结合位点语法以及细胞类型特异性的转录因子表达。

Ubiquitous motifs reduce accessibility

Para_01

1. 大多数从头预测的基序具有正向贡献分数，但一小部分（15个基序）被预测具有负向效应（图5a、b）。
2. 负向基序出人意料地广泛存在于可及区域中，每个峰内有2至5个具有预测性的实例，并且在多个器官中占所有预测性基序实例的三分之一以上（图5b和扩展数据图5e–g）。
3. 许多此类基序与已知的具有抑制活性的转录因子家族相匹配，包括ZEB或SNAIL（ZEB/SNAIL）、HIC、BCL11A、NFY以及YY1/2，尽管它们对可及性的普遍负向效应此前未见报道（图5c），而其他一些则未能匹配到已知数据库中的条目（扩展数据图6b）。

Fig. 5: Ubiquitous motifs associated with reduced accessibility and directional eQTL effects.

- 图片说明

◉ a，肌细胞中MYOD1位点的观察和预测染色质可及性、推断的核小体占位、贡献分数以及基序实例注释，显示了一个负向贡献的ZEB/SNAIL基序示例。◉ b，按器官分类的从头基序分解，以及每类器官中正负基序的基序实例分解。◉ c，每类负向基序的从头基序（以CWM形式表示），以及外部数据库中相似度最高的已知PWM。◉ d，左侧和中间：热图显示从推断的核小体二联中心位置（左）和峰顶（中）每10个碱基对区间的基序实例计数，每个热图中每个基序的计数按距离区间进行z-score标准化；右侧：基因组实例重叠各种基因组特征的比例。每一行代表一个广泛的基序大类，仅展示在所有区间中实例总数超过50,000的类别。◉ e，确定不同类别胎儿基序在具有上调或下调eQTL变异中的组织特异性富集的工作流程。◉ f，eQTL变异分析结果。对于每个广泛的基序类别，点表示来自每个器官的每个独特基序的富集分数，颜色根据下调eQTL变异（左）和上调eQTL变异（右）的−log10 FDR值着色。堆叠条形图显示具有富集的基序实例所来源的器官比例。仅展示统计学上显著的富集结果。

Para_02

1. 值得注意的是，NFY和YY1/2基序具有双重作用，不同的基序实例可分别驱动正向或负向效应（图5c,d）。
2. 在正向基序实例上聚合的足迹显示出典型的特征，但在负向实例上则表现出相对于侧翼区域可及性降低的现象（扩展数据图9a）。
3. 因此，负向基序实例代表了在染色质可及区域内局部可及性降低的位点。
4. 与正向基序以及启动子元件基线相比，负向基序实例在不同细胞类型间的符合程度较低，提示存在细胞类型特异性的调控机制（扩展数据图9b,c）。
5. 在同一基序的不同细胞类型中，其正向与负向实例之间的重叠极为罕见，表明不同的基因组位点分别对可及性产生正向或负向影响（扩展数据图9d）。
6. 负向基序在核小体二联体附近富集，而在峰顶附近则减少；相比之下，正向基序倾向于避开二联体，并在峰顶附近聚集（图5d和扩展数据图5j）。
7. 接着，我们对每种基序的1000个实例进行了计算机模拟敲除分析，发现敲除正向基序会导致预测可及性下降，而敲除负向基序则导致预测可及性轻微上升（扩展数据图9e）。
8. 因此，具有相同序列特征的非编码区域可能对染色质可及性产生相反的影响。

Para_03

1. 为了评估基序对下游基因表达的影响，我们计算了来自基因型-组织表达（GTEx）数据集中与相应胎儿器官中预测的正向和负向基序实例重叠的表达数量性状位点（eQTL）中精细定位变异的富集情况（图5e及方法部分）。
2. 破坏正向基序的变异应会降低基因表达，而破坏负向基序的变异则应会增加表达。
3. 对于每个胎儿器官，我们量化了来自相关GTEx组织的上调和下调eQTL基因变异对在每个基序上的富集程度，并与随机打乱的背景集进行了比较。
4. 尽管只有一小部分eQTL变异与我们预测的基序实例重叠（扩展数据图9f,g），但这些富集结果与我们的假设一致。
5. 与正向基序重叠的eQTL在下调变异中显著富集（双侧Fisher精确检验，P = 9.91 × 10−11，优势比（OR）= 0.047），而与负向基序重叠的eQTL则在上调变异中富集（P = 8.01 × 10−4，OR = ∞）（图5f和补充表8）。
6. 我们注意到存在组织特异性的模式，在心脏中观察到GATA基序变异的富集，在肝脏中观察到HNF4基序变异的富集，在甲状腺中观察到NKX基序变异的富集（图5f）。

Disease variants in regulatory elements

Para_01

1. 数千项全基因组关联研究（GWAS）已将基因组位点与复杂的表型联系起来，但确定致病变异及其所影响的相关细胞类型仍然具有挑战性。
2. 因此，我们探究了胎儿细胞类型中的可及染色质图谱是否富集了与人类疾病相关的遗传变异（图6a）。
3. 我们将HDMA中的候选顺式调控元件与CAUSALdb中的精细定位变异进行整合分析，后者汇总了跨越349种不同性状的1483项GWAS的可信集合。

Fig. 6: Disease causal variants overlap motifs in fetal-only peaks.

- 图片说明

◉ a，识别在胎儿而非成年峰集中的疾病相关因果变异体的工作流程。◉ b，g-chromVAR富集结果中的前几位命中项，仅显示在相似MeSH术语中具有最高平均z分数的疾病特征，以及按L3细胞类型注释进一步分组的具有最高总z分数的L2胎儿细胞类型。◉ c，哮喘相关变异位点rs113892147是胎儿肺巨噬细胞中一个仅在胎儿中出现的阳性NRF1基序中的命中位点，预测其会降低染色质可及性。轨道从上到下分别为：编码基因、GWAS单核苷酸多态性（SNP）的原始GWAS P值、胎儿肺巨噬细胞簇的观测伪批量可及性、成人肺巨噬细胞中单核ATAC-seq（snATAC-seq）得到的峰、非效应等位基因和效应等位基因的预测可及性、PhyloP保守性评分、非效应和效应等位基因对预测可及性贡献的每碱基评分，以及非效应等位基因的基序实例。◉ d，冠状动脉疾病（CAD）相关变异位点rs12740374是肌肉内皮细胞中一个仅在胎儿中出现的与ZEB/SNAIL负性基序重叠的命中位点，预测其通过创建一个C/EBP位点来增加可及性（相应的基序作为插图显示在右下角）。轨道内容与c相同。

Para_02

1. 我们使用g-chromVAR60量化了所有胎儿L2细胞类型中与性状相关的变异富集情况，并鉴定出131种具有显著富集的细胞类型（错误发现率（FDR）< 0.05）（补充表9）。
2. 为了聚焦于与疾病相关的信号，我们整理了一个包含68种细胞类型的子集，这些细胞类型在来自36个医学主题词（MeSH）疾病术语的79项研究中的变异表现出富集（补充表10）。

Para_03

1. 我们观察到，变异富集主要具有细胞类型特异性，而与器官无关：来自不同器官的相似细胞类型往往富集于相同的疾病性状（图6b）。例如，免疫细胞类型富集了与甲状腺疾病、皮炎和湿疹、鼻炎以及癌症相关的变异；肌细胞与心房颤动风险相关，而内皮细胞则富集了与高血压疾病相关的变异；特化的肝脏实质细胞，如肝细胞和胆管细胞，富集了与脂质代谢紊乱相关的变异。

Variant effect on local accessibility

Para_01

1. 为了研究遗传变异如何扰动调控，我们优先选取了高置信度的因果变异（SuSiE69 后验包含概率（PIP）≥ 0.8），这些变异与疾病富集的胎儿细胞类型中可及元件内的预测转录因子结合位点重叠（补充表11）。
2. 为了区分胎儿特异性效应，我们将这些变异与HDMA caCREs的重叠情况同来自ENCODE70的匹配成人细胞类型的ATAC-seq峰进行了比较（图6a和补充表12）。
3. 我们鉴定出28个仅在胎儿中出现的信号，以及80个在胎儿和成人峰集中均存在的信号。

Para_02

1. 使用 ChromBPNet 模型，我们预测了28个胎儿特异性变异对局部染色质可及性的影响（补充表13）。
2. 许多仅在胎儿中出现的变异被预测会通过破坏基序来改变可及性（扩展数据图10）。
3. 我们重点介绍了其中两个此类变异（图6c,d），其中一个涉及在肺泡巨噬细胞中与哮喘相关的正向NRF1基序（z得分=3.27，FDR=2.09×10−2），另一个是在肌肉内皮细胞中与冠状动脉疾病相关的负向ZEB/SNAIL基序（z得分=3.45，FDR=1.42×10−2）。

Para_03

1. 与哮喘相关的精细定位变异位点（rs113892147，SuSiE PIP = 0.819）重叠了一个在胎儿肺巨噬细胞可及区域内的保守NRF1结合基序（图6c），但在成人肺巨噬细胞样本中缺乏可及性。
2. 胎儿肺巨噬细胞ChromBPNet模型预测，相较于参考等位基因，效应A等位基因的引入会降低该区域的可及性，这与eQTL数据一致，后者将该等位基因与成人肺组织中ARHGEF25表达水平降低相关联（beta = −0.295，P = 6.9 × 10−7，Open Targets Platform）。
3. 尽管该发现与哮喘病理之间的直接联系尚不明确，但这一结果强调了胎儿巨噬细胞增强子可能在塑造哮喘风险中发挥的作用。

Para_04

1. 与冠状动脉疾病（CAD）相关的变异位点（rs12740374，SuSiE PIP评分=0.855）重叠于胎儿肌肉内皮细胞中的一个负性ZEB/SNAIL基序（图6d）。
2. 效应G等位基因破坏了较弱的负性ZEB/SNAIL基序，并创建了一个C/EBP结合位点，预测将导致染色质可及性增加。
3. 该位点（CELSR2–PSRC1–SORT1）在CAD风险中已有明确研究基础，且该变异可增强肝细胞中C/EBP的结合以及SORT1的表达。
4. 我们的研究结果将其潜在作用扩展至内皮细胞，这与近年来强调血管功能障碍在CAD发病机制中作用的新模型一致。
5. 值得注意的是，除肝脏外，rs12740374在肌肉组织中与邻近基因CELSR2（beta = 0.338，P = 1.2 × 10−36）和PSRC1（beta = 0.238，P = 4.6 × 10−8）表达升高相关，但在心室或心房组织中无此关联（Open Targets Platform），进一步凸显非心肌细胞在CAD中的潜在作用。

Discussion

Para_01

1. 在此，我们通过一个多组学的人类发育图谱定义了染色质可及性和转录组景观。
2. 我们在体内解析了器官特异性和细胞类型特异性的增强子活性，采用深度学习策略识别调控基序，揭示了转录因子协同作用的不同模式，并解析了疾病相关遗传变异对发育过程中染色质景观的影响。
3. 我们的工作通过揭示多个器官中单核苷酸分辨率、细胞类型特异性的调控逻辑，显著拓展了人类发育中调控图谱研究的分辨率和范围。

Para_02

1. 我们的深度学习框架与经典的基序发现方法相辅相成，例如来自DNase足迹分析或体外结合实验的基序富集。
2. 首先，与需要深度测序的足迹分析方法不同，ChromBPNet模型可以在单细胞簇的伪批量数据上进行训练和解释，即使测序覆盖度适中也能实现对原代组织中具有细胞环境特异性的稳健基序发现。
3. 我们的模型能够学习可能具有因果关系的序列特征，这些特征可直接预测染色质可及性，并支持在计算机中进行扰动分析，以系统地探究基序结构，从而提供机制层面的洞察。
4. 我们的从头基序分析结果与以往研究一致，表明有限的一组核心基序调控了大多数转录起始事件。

Para_03

1. 我们发现了67对在染色质可及性上表现出协同效应的基序对，其中48对对特定的间距和方向具有强烈偏好，表明存在严格的基序语法，类似于经典的IFNβ‘增强体’模型以及具有刚性结合结构的AP-1–IRF4（AICE）和ETS–IRF（EICE）复合元件。
2. 以往利用体外系统对转录因子协同作用进行的系统性研究也发现了类似的间距和方向偏好，但缺乏体内分辨率，且通常集中于结合位点融合或重叠的复合基序。
3. 我们分析所预测的具有协同作用且具组织特异性的复合基序为深入的机制研究提供了候选目标，并支持这样一种观点：密集排列且空间受限的转录因子相互作用是发育基因程序中相对常见的调控机制。

Para_04

1. 与此同时，具有松散语法的27个基序与间接协同作用的生物物理模型一致，这些模型包括转录因子与核小体竞争的质量作用模型，或招募能够驱逐核小体的染色质重塑复合物的模型，这与单分子足迹分析、CAP-SELEX（连续亲和纯化配体指数富集系统进化）实验以及ChIP-exo（染色质免疫沉淀结合外切酶消化）数据的深度学习分析结果相符。
2. 这种语法灵活性可能赋予更强的进化鲁棒性，使得调控元件在容忍基序间距、数量或方向发生改变的同时仍能保持其功能。

Para_05

1. 深度学习模型的解释识别出了一些负性基序，这些基序会降低染色质可及性，在核小体dyad区域和远端区域表现出明显的定位偏好，并且在上调的eQTL变异中显著富集。
2. 独立证据支持负性基序的调控作用：果蝇中一个ZEB2抑制子基序降低了增强子活性；一个新生变异形成了具有抑制作用的NR2F1基序，导致小鼠大脑中的报告基因表达下降；在发育中的肝脏中，富含ZEB1、TCF4和SNAI1等抑制性基序的区域预测到较低的增强子活性和基因表达水平，提示其可能参与调控的精细调节。
3. 在我们的数据集中，负性基序包括已知的抑制因子以及此前未被表征的、不同于沉默子元件的序列元件，表明存在更广泛的抑制信号参与可及性的调控。
4. 未来的研究需要鉴定出与这些负性基序结合的蛋白质，明确它们与染色质重塑复合物或核小体之间的相互作用，并阐明其如何调控染色质可及性和基因表达。

Para_06

1. 我们观察到YY1/2和NFY的基序可根据上下文发挥可及性的正向或负向调控作用，这与YY1已知的双重功能一致：在不同情况下既可作为激活因子也可作为抑制因子，同时也与NFY作为组蛋白折叠蛋白在活性调控元件处维持无核小体区域的作用相符。

Para_07

1. 尽管eQTL分析表明，正向基序的破坏通常会降低基因表达，这与激活输入的丧失一致，但有两个正向基序——RFX基序和一个SP/KLF类似位点——在上调eQTL中富集。
2. 这些特定的eQTL可能反映了亲和力优化的变异，导致转录因子结合能力增强以及基因表达升高。
3. 或者，这些转录因子可能在增加染色质可及性的同时抑制基因表达。
4. 事实上，RFX家族成员在不同组织中调控不同的基因集合，并在特定调控元件上作为抑制因子发挥作用。

Para_08

1. 最后，我们将与疾病相关的变异与特定的发育阶段细胞类型联系起来，并预测了它们对染色质可及性的影响。
2. 许多与成年期性状相关的变异位于胎儿特异性的可及区域，这表明即使受影响的染色质元件在成体组织中变得不可及，遗传变异仍可能通过干扰相关细胞谱系的发育而影响疾病发生。
3. 胎儿肺巨噬细胞中的哮喘相关变异提示，早期起作用的变异可能会损害起源于胎儿肝脏单核细胞并持续至成年的肺泡巨噬细胞的发育。
4. 或者，某些变异可能在生命后期导致胎儿调控元件的异常再激活，从而增加患病易感性。
5. 这些发现强调了发育背景在疾病病因中的重要性。

Para_09

1. 我们的研究存在局限性。
2. 尽管我们的图谱涵盖了广泛的细胞类型，但仍需要更深入的取样才能全面解析发育组织中的全部多样性。
3. 我们的数据并未在所有经过实验验证的VISTA增强子中重现器官特异性的增强子活性，这可能是因为这些增强子仅在小鼠胚胎第11.5天进行了检测，该时期在转录上最接近人类受孕后第5周，远早于我们所分析的样本时期。
4. 此外，VISTA集合仅代表一小部分强效的、组织特异性的增强子，因此从这些元件得出的观察结果可能无法推广到全部注释的候选顺式调控元件（caCREs）。
5. 基于染色质可及性训练的深度学习模型主要捕捉的是调控可及性的序列特异性DNA结合蛋白的影响，可能遗漏那些在可及性下游或独立于可及性发挥作用的调控因子。
6. 转录因子家族内部的基序简并性使得无法明确地将去新发现的基序归因于特定的单个转录因子。
7. 最后，尽管我们的变异效应预测提供了机制性假设，但仍需实验验证来确认这些预测对染色质状态和基因调控的影响。
8. 我们的工作代表了一个全面的、多组学的单细胞图谱，涵盖了人类胎儿发育过程，并提供了一个预测框架，用于解析跨细胞类型调控染色质可及性的顺式调控语法。
9. 我们提供了一个丰富的数据资源、训练好的模型、兼容基因组浏览器的轨道文件以及一个高效的预处理流程，以支持对人类发育过程中调控机制和疾病机制进行详细探索。

Methods

Ethics statement

伦理声明

Para_01

1. 去标识化的组织样本由斯坦福大学医学院在选择性终止妊娠手术中采集，并已获得用于组织研究的知情同意，且遵守相关的法律和机构伦理规定。
2. 未收集任何人口统计学信息。
3. 知情同意由医疗团队获取。
4. 相关组织样本的处理和分析是在斯坦福大学干细胞研究监督委员会（SCRO）批准的协议SCRO-796下进行的。

Sample collection and nuclei isolation

样品采集和细胞核分离

Para_01

1. 组织样本在冰上运送，并在处理前立即储存于液氮中。
2. 开发了一种适用于多种组织的细胞核分离方案，以高效分离稳定的细胞核用于后续文库构建。
3. 简言之，对于特定样本，将100–200 mg组织直接加入预冷的2 ml玻璃匀浆器（Kimble 885300-0002）中，内含1 ml细胞核提取缓冲液（250 mM蔗糖，25 mM KCl，5 mM MgCl2，20 mM HEPES-KOH，65 mM β-甘油，0.5% IGEPAL CA-630，1×蛋白酶抑制剂，1 mM DTT，0.2 mM亚精胺，0.5 mM精胺，60 U/ml RNasin Plus，2–5%正常山羊血清），并将匀浆器置于冰上。
4. 样本在冰上孵育10分钟。
5. 使用A型杵头匀浆20次，随后使用B型杵头匀浆20次。
6. 将样本转移至低吸附DNA管中。
7. 用额外300 µl细胞核提取缓冲液冲洗匀浆器中残留的细胞核并收集。
8. 样本在4 °C下以垂直旋转方式孵育5分钟。
9. 使用70 µm Flowmi滤器过滤样本。
10. 用额外的细胞核提取缓冲液调节总体积至1.2 ml。
11. 向样本中加入37%甲醛，使其终浓度为0.2%，并在室温下以垂直旋转方式孵育4分钟。
12. 用125 mM甘氨酸在室温下垂直旋转孵育8分钟以终止固定反应。
13. 向样本中加入细胞核提取缓冲液，使最终体积达到1.4 ml。
14. 制备碘克沙醇梯度以从匀浆液中富集细胞核。
15. 简言之，通过将60%碘克沙醇溶液稀释，并添加1 mM DTT、60 U/ml RNasin Plus和2–5%正常山羊血清，配制50%碘克沙醇溶液。
16. 将样本与适量碘克沙醇混合，使其最终碘克沙醇浓度为22%。
17. 在样本下方铺一层44%碘克沙醇溶液。
18. 然后，在样本与44%碘克沙醇层之间轻轻加入一层22%碘克沙醇溶液。
19. 样本在4 °C下以3,500g离心30分钟，关闭刹车。
20. 通过轻柔吹打分离出细胞核层，用于后续处理。

SHARE-seq library preparation

SHARE-seq文库制备

Para_01

1. 完整方案详见补充说明1，该方案改编自已发表的SHARE-seq方案。
2. 总共我们处理了来自23名个体的76个组织样本，这些样本跨越了12个组织处理和SHARE-seq文库制备批次，每个批次对应于特定器官的所有样本。

Library sequencing

文库测序

Para_01

1. 所有DNA文库均使用NovaSeq 6000平台，配合300循环S4 v1.5试剂盒及XP工作流程进行测序。
2. 采用96-99-8-96的配对末端测序配置（Read1-Index1-Index2-Read2）。
3. 我们使用Qubit和Tapestation对DNA文库进行定量，随后将文库混合至1.5 nM的浓度，最终上机浓度为300 pM。
4. 测序工作在斯坦福基因组技术中心完成。

VISTA embryo histology

VISTA胚胎组织学

Para_01

1. 我们从L. Pennachio处获得了在PBS中保存的X-Gal染色固定小鼠全胚胎，并将其转移至70%乙醇中储存。
2. 石蜡包埋由Histo-Tec实验室完成，采用无二甲苯脱水方案，因为二甲苯可能溶解X-Gal染色。
3. 具体而言，胚胎依次用80%、95%、100%、100%和100%乙醇进行脱水，每次20分钟，
4. 随后用50:50、80:20、90:10和100:0的石蜡：酒精混合液分别洗涤20分钟，以去除乙醇。
5. 之后对5微米切片采用标准流程进行包埋和H&E染色。

SHARE-seq data pre-processing

SHARE-seq 数据预处理

Para_01

1. 我们开发了一个高度并行化、快速且存储高效的预处理Snakemake（v7.15.1）流程，用于将测序仪的BCL文件转换为ATAC片段文件和RNA稀疏矩阵（扩展数据图1）。
2. 简而言之，原始BCL文件首先通过自定义脚本转换为FASTQ文件，该脚本按流动池的区域对bcl2fastq（v2.20.0.422，Illumina）转换过程进行并行化处理，解析读段周期，并根据Index2读段中的子文库条形码将原始FASTQ文件拆分为各个子文库。
3. 对于每个子文库，我们进一步将FASTQ文件随机分割为每份包含两千万条读段的小文件。

Para_02

1. 在每个ATAC-seq子文库的片段中，我们针对SHARE-seq条形码白名单进行条形码匹配，允许构成单细胞条形码的三轮8 bp条形码每轮存在1个碱基对的错配，随后使用fastp（v0.23.2）进行Nextera接头切除，通过Bowtie2（v2.5.0）进行基因组比对，并将输出的BAM文件转换为更节省存储空间的片段文件。
2. 然后，我们将一个子文库所有片段的片段文件合并，根据起始和终止坐标对每个细胞的片段进行去重，并基于第一轮细胞条形码将片段按样本解复用。
3. 最后，对于每个样本，我们将该样本在所有子文库中的解复用后片段文件合并，生成最终的ATAC-seq片段文件（.fragments.tsv.gz，.fragments.tsv.gz.tbi）。

Para_03

1. 在每个RNA亚文库的片段中，我们进行了条形码匹配，从Read2中解析10个碱基的UMI，并仅对Read1进行接头修剪，随后使用STAR（v2.5.4b）进行基因组比对，利用featureCounts（v2.0.1）进行基因注释，并将输出的BAM文件转换为更节省存储空间的TSV格式。
2. 然后，我们将同一亚文库所有片段的注释后TSV文件合并，根据第三轮条形码将其划分为12个条形码片段，针对每个条形码片段中的每个细胞、每个注释基因使用UMI-tools（v1.1.2）进行UMI去重。
3. 接着，基于第一轮细胞条形码，将去重后的TSV文件拆分为不同的样本，并将这些TSV文件转换为矩阵市场交换格式。
4. 最后，对于每个样本，我们将该样本在所有亚文库中的解复用后的矩阵市场交换文件合并，生成最终的RNA稀疏矩阵文件（*.matrix.mtx.gz, *.features.tsv.gz, *.barcodes.tsv.gz）。

Para_04

1. 平均而言，我们使用学术高性能计算集群可以在不到4小时内处理一个100亿读长的NovaSeq运行。
2. 该流程可以轻松适应其他基于分裂与池化（split-and-pool）的单细胞多组学数据的处理。
3. 所有文库均比对到hg38参考基因组。
4. 该流程可在https://github.com/GreenleafLab/shareseq-pipeline获取（稳定版本v1.0.0）。
5. 原始测序读长在公开存放之前已使用BAMboozle98进行匿名化处理，以保护供体隐私。
6. 匿名化代码可在shareseq-pipeline GitHub仓库的‘anonymize’分支上获得。

Para_05

1. 众所周知，在ATAC–seq实验中，Tn5转座酶会形成同源二聚体，两个Tn5分子之间存在9个碱基对的间隙，导致在每个可及位点的不同DNA链上产生相距9个碱基对的两次插入。
2. 当使用成对末端测序对DNA片段进行测序时，需要根据Tn5的插入偏移量调整起始和终止位置，以反映Tn5二聚体中点处可及位点的真实中心。
3. 为了校正Tn5的偏移，以往的ATAC–seq研究采用+4/−5偏移方法，即正链上的插入位点向前调整+4个碱基对，负链上的插入位点向后调整−5个碱基对。
4. 然而，这种方法实际上会导致来自同一次转座事件的两个片段在调整后的插入位点之间出现1个碱基对的错配（扩展数据图1b）。
5. 这种差异可能源于BAM和BED文件使用的端点排他性坐标系统，因为原始的+4/−5惯例仅在输出文件被解释为非标准的基于0且包含终点的基因组坐标系统时才是正确的。
6. 这种错配不会影响大多数将插入位点归入数百个碱基对尺度区间的下游ATAC–seq分析，但会影响对碱基对敏感的分析，例如转录因子足迹分析和基序分析。
7. 在此SHARE-seq预处理流程中，我们采用了+4/−4的偏移方式，从而获得一致的插入位点。关于使用这两种偏移方案校正读段后生成的示例基序，请参见参考文献17中的补充图3a。

SHARE-seq data QC and filtering

SHARE-seq 数据的质量控制与过滤

Para_01

1. 我们通过手动检查并设定以下指标的阈值，对每个样本进行了质量控制（QC）过滤：（1）ATAC-seq片段文件的TSS富集比和片段数量；（2）RNA稀疏矩阵的UMI数量、基因数量以及线粒体reads所占百分比；（3）RNA UMI与ATAC-seq片段数的比值，以去除某一检测模式下质量较低的细胞（扩展数据图2a）。所有样本过滤的阈值均汇总于补充表1中。未进行明确的批次效应校正，因为个体特异性效应通常与细胞类型组成的时间差异相互混淆，难以区分这些变异来源。

RNA normalization, ambient RNA removal, dimensionality reduction and clustering

RNA标准化、环境RNA去除、降维和聚类

Para_01

1. 我们使用R（v4.1.2）中的Seurat（v4.3.0）将每个器官过滤后的RNA稀疏矩阵处理为Seurat对象。
2. 我们采用迭代降维和聚类的工作流程，依次注释细胞类型并过滤掉额外的低质量聚类（扩展数据图2a）。
3. 在每次迭代中，我们首先对每个样本的RNA原始计数进行SCTransform v2和可变特征选择。
4. 然后使用Seurat中的SelectIntegrationFeatures函数，在所有样本中选择前3000个共有的可变特征，并排除线粒体基因、性染色体基因和细胞周期相关基因，以最大程度减少批次效应。
5. 我们将来自每个样本对象的原始RNA计数合并成一个单一矩阵，使用共有的特征进行SCTransform v2校正。
6. 然后使用celda（v1.6.1）包中的DecontX函数对SCT校正后的计数进行去污染处理，以去除每个细胞的环境RNA污染。
7. 去污染后的计数按样本拆分，标准化至每个细胞10,000个UMI，并进行对数归一化。
8. 类似于上述过程，我们从各样本的可变特征中再次选出一组前3000个共有的可变特征。
9. 基于这些共有特征，在合并的对象上进行主成分分析。
10. 随后使用Louvain算法在50个主成分的基础上进行细胞聚类，聚类分辨率为0.3，并通过均匀流形近似与投影（UMAP）进行嵌入可视化。
11. 接着我们检查每个聚类，移除那些与其他聚类相比UMI数量显著偏低、不同组织区室标志基因共表达水平异常高（在生物学上不可能，提示为双细胞）或缺乏明确细胞类型定义标志基因的低质量聚类。
12. 例如，同时高表达上皮和内皮区室标志基因的聚类被视为可疑双细胞而被剔除。
13. 在移除低质量聚类中的细胞后，我们从每个样本的RNA原始计数开始重复上述处理步骤。
14. 该过程反复执行，直到不再发现新的低质量聚类，通常需要一到三次迭代。
15. 最终通过此迭代质控流程保留下来的细胞聚类中的细胞被视为‘白名单’细胞。
16. 对于每个细胞类型聚类，我们通过Wilcoxon秩和检验（一类对所有其他类）识别标记基因，并筛选出log2倍变化大于1的基因作为有效标记基因。

Para_02

1. 所有最终的聚类注释均包含在补充表2中。
2. 所有聚类标志物汇总于补充表3中，其中一部分在补充说明2的标志基因点图中进行了可视化展示。
3. 所有UMAP嵌入图均包含在补充说明2中。

ATAC–seq peak calling, motif enrichment and chromVAR

ATAC–seq峰识别、基序富集和chromVAR

Para_01

1. 我们使用 ArchR（v1.0.2）将过滤后的 ATAC 片段文件按器官处理为独立的 ArchR 项目。
2. 在通过迭代的 RNA 处理流程筛选出最终的白名单细胞条形码后，我们将聚类结果和细胞类型注释进行转移，并使用 Macs2（v2.2.7.1）对每个聚类调用峰区域。
3. 随后，利用 ArchR 的迭代重叠策略，将各聚类的峰区域合并为每个器官内统一的可重复峰集合，并构建了基于片段计数的细胞-峰矩阵。
4. 我们通过 Wilcoxon 秩和检验识别每个聚类的标志峰，并在标志峰内进行转录因子基序富集分析，显著性阈值设定为 FDR ≤ 0.1 且 log2（倍数变化）≥ 0.5。
5. 我们在每个器官的所有聚类中计算了 chromVAR 基序偏差。
6. 上述两项分析均使用了一个经过整理的 cisBP 基序集，包含 1141 个独特的人类转录因子基序。
7. 通过合并全部 12 个器官各自的 ArchR 项目，我们建立了一个全局 ArchR 项目，并通过迭代重叠来自所有器官各个聚类所调用的峰集合，构建了一个 HDMA 全局 caCRE 集合。

Linkage of regulatory elements to genes with modified ABC model

通过改进的ABC模型将调控元件与基因关联

Para_01

1. 我们采用ABC方法将caCRE与基因启动子进行关联。
2. 为了确保ABC增强子区域与我们HDMA全局caCRE集合中的一致，我们对ABC24进行了调整，使其能够接受自定义区域作为模型输入。
3. 在构建自定义区域集合时，我们使用bedtools将HDMA全局caCRE集合与hg38基因组的TSS集合合并。
4. 根据ABC官方文档推荐的单细胞ATAC分析流程（https://abc-enhancer-gene-prediction.readthedocs.io/en/latest/index.html），我们使用伪批量ATAC-seq信号（不包含H3K27ac染色质免疫沉淀测序（ChIP-seq））作为增强子活性，并利用基因组距离的幂律函数估算增强子与启动子之间的三维接触频率。
5. 我们在每个L1细胞簇上分别运行ABC模型。
6. 随后对结果进行筛选，保留ABC评分大于0.013的增强子-启动子连接，该阈值对应于基准CRISPR数据集中70%的召回率。
7. 我们修改后的ABC工作流程可在以下地址获取：https://github.com/GreenleafLab/ABC-Enhancer-Gene-Prediction-CustomRegions（commit b3d2156）。

VISTA enhancer analysis

VISTA增强子分析

Para_01

1. 我们筛选了结果以获得VISTA验证的增强子（于2024年1月24日访问），这些增强子来源于人类或具有人类序列同源物，并在HDMA所包含的任何器官中被注释为X-Gal阳性，且与HDMA全局caCREs重叠，并保留至少75%（375 bp）重叠度的VISTA增强子。
2. 如果多个caCRE与同一个VISTA增强子重叠，我们选择ATAC–seq信号最强的caCRE，依据是先前的观察表明VISTA增强子通常具有更小的核心元件，且增强子活性不依赖于增强子内所有调控元件。
3. ATAC–seq信号按L1簇进行缩放和对数归一化，然后对每个增强子在各簇间进行z-score标准化。
4. 对于每个器官，采用单侧Wilcoxon秩和检验计算在该器官中注释为阳性的VISTA增强子所重叠的caCRE中HDMA ATAC–seq信号富集的统计显著性。
5. 例如，为了检验脑部ATAC–seq信号富集的显著性，我们首先仅选取HDMA脑组织簇，然后使用单侧Wilcoxon秩和检验比较与VISTA脑阳性增强子重叠的caCRE和与VISTA脑阴性增强子重叠的caCRE中的ATAC–seq信号。
6. 效应量计算公式为W统计量除以（n1 × n2），其中n1为脑阳性组中caCRE的数量，n2为脑阴性组中caCRE的数量。
7. 该效应量代表脑阳性组中某个caCRE的ATAC信号高于脑阴性组中某个caCRE的概率，即AUROC概率。
8. 类似地，对于RNA数据，我们首先确定每个与VISTA增强子重叠的caCRE的最近基因，对原始RNA表达计数按L1簇进行缩放和对数归一化，然后对每个增强子在各簇间对表达值进行z-score标准化。
9. 对每个器官执行类似的Wilcoxon秩和检验，以评估VISTA阳性增强子中HDMA RNA信号富集的统计显著性。
10. 为了避免在Wilcoxon秩和检验中因机器精度导致数值下溢至零，P值在log10尺度上计算并相应报告。

Preparation of input regions for ChromBPNet models

为ChromBPNet模型准备输入区域

Para_01

1. 为了定义用于训练 ChromBPNet 模型的基因组区域，我们进行了第二轮峰值 calling 以获得一组宽松的可及区域。
2. 首先，通过将 SHARE-seq ATAC-seq 模态中属于每个 L1 细胞类型簇的所有细胞（来自所有样本）的片段进行拼接，生成了 203 个 L1 细胞类型簇各自的伪批量片段文件。
3. 对于每个伪批量数据，我们进一步生成了伪重复样本。
4. 针对每个片段，其起始位点和终止位点（对应 Tn5 插入位点）被随机分配到两个伪重复文件中的其中一个。
5. 同时，这两个伪重复文件也被合并为一个总的伪重复文件。
6. 使用 Macs2（v2.2.9.1）对这三个伪重复文件进行峰值识别，参数设置为：-p 0.01 –shift -75 –extsize 150 –nomodel -B –SPMR –keep-dup all –call-summits。
7. 仅保留那些在总伪重复文件中识别出的峰值，并且这些峰值必须同时在两个伪重复文件中也被识别出。
8. 排除与 GRCh38 ENCODE 黑名单（ENCODE 编号 ENCFF356LFX）重叠的峰值。
9. 将峰值坐标调整为以 Macs2 峰值峰顶为中心、长度为 1,000 bp 的区域。
10. 伪重复样本仅用于峰值识别，而伪批量片段文件则用于后续模型训练。

Para_02

1. 我们使用了 ChromBPNet 软件包（https://github.com/kundajelab/chrombpnet，commit a5c231），并遵循了所描述的工作流程。
2. 我们使用命令 chrombpnet prep nonpeaks 来定义与峰区域 GC 含量相匹配的背景区域。
3. 对于每种细胞类型，我们采用了五折交叉验证方案，其中每一折（标记为 0 到 4）包含不同的训练、验证和测试染色体集合，每个染色体至少出现在一个折的测试集中。
4. 我们使用了 ChromBPNet 提供的默认人类染色体划分方案。

ChromBPNet model training

ChromBPNet 模型训练

Para_01

1. ChromBPNet 模型是监督式卷积神经网络，经过训练用于利用峰值区域和背景区域中2,114个碱基对的一-hot编码DNA序列，预测输入区域中心1,000个碱基对窗口内的可及性谱型（以概率分布形式）和总的自然对数计数（以标量值形式）。
2. ChromBPNet 模型使用一个预训练的偏差模型，用以解释Tn5转座酶偏差未能捕获的残余可及性。
3. 我们训练了一个偏差模型，利用SHARE-seq设置下的心脏L1簇0伪批量数据（fold 0），通过chrombpnet偏差流程并设置偏差阈值因子为-b 0.4，来学习酶切偏差，该流程同时使用DeepLIFT进行模型解释。
4. 我们确认该偏差模型学习到了Tn5的基序特征，但未学习到转录因子的基序特征，并将此偏差模型用于后续训练每种203种L1细胞类型下五个fold的ChromBPNet模型（共1,015个模型），使用chrombpnet流程命令并结合ENCODE提供的GRCh38参考基因组（fasta: https://www.encodeproject.org/files/GRCh38_no_alt_analysis_set_GCA_000001405.15/，染色体大小来自ENCODE编号ENCFF667IGK）。

ChromBPNet model evaluation

ChromBPNet 模型评估

Para_01

1. 模型根据预测值与在保留的测试集染色体上的峰区中观察到的对数值之间的皮尔逊相关系数和斯皮尔曼相关系数，以及预测与观测图谱在峰区内的詹森-香农距离进行评估（补充表5和补充说明2）。
2. 对于4种细胞类型，若任一折叠中的斯皮尔曼相关系数低于0.5，通常对应于测序深度较低的伪批量样本，则排除这些模型。
3. 为了生成峰区跨折叠的平均预测可及性轨迹（表示每碱基的计数），对每个区域，将跨折叠的平均预测图谱logits通过softmax函数处理，转化为概率，然后乘以跨折叠的平均预测对数计数的指数值进行缩放。

ChromBPNet model interpretation

ChromBPNet 模型解释

Para_01

1. 我们进行了模型解释，以确定每个核苷酸对可及性的预测程度。
2. 我们运行了chrombpnet interpret命令，该命令使用DeepLIFT算法计算2,114个碱基对输入窗口中每个核苷酸相对于预测计数的贡献分数。
3. 针对所有峰区区域，每个模型折叠均计算出贡献分数，并进一步计算跨折叠的平均值。
4. 平均后的预测可及性谱图和贡献分数被转换为bigWig文件，并用于所有分析和图表绘制。

De novo motif discovery and assembly of the motif lexicon

从头基序发现及基序词典的组装

Para_01

1. 从头基序词典的构建需要三个步骤：（1）每种细胞类型中的从头基序发现；（2）在不同细胞类型间合并相似的基序；以及（3）基序注释。

[ol]- 1) First, for de novo motif discovery on sequences driving chromatin accessibility, we used TF-MoDISco31, which, in brief, identifies seqlets, corresponding to short spans of contiguous high positive-importance or high negative-importance nucleotides, and clusters them into recurrent 30 bp patterns. We used the implementation in the TF-MoDISco-lite package (https://github.com/kundajelab/tfmodisco, v2.0.7) on the mean contribution scores for each cell type, sampling 1,000,000 seqlets for both positively and negatively contributing seqlets (parameter -n 1,000,000), and using the default behaviour to search for seqlets in the central 400 bp of input regions. Each de novo motif is represented by a 4 × 30 CWM, computed as the mean contribution score per position per nucleotide across its seqlets, and a 4 × 30 position probability matrix (PPM), computed by normalizing the nucleotide frequencies per position across its seqlets. We manually inspected CWMs learned in each cell type, and the ChromBPNet models and motifs for ten cell types which predominantly learned low-complexity motifs were excluded from downstream analysis. This resulted in 189 cell types (945 models total) passing quality control and used throughout our study, which collectively learned 6,362 motifs including both positively contributing and negatively contributing motifs. - 2) Next, we automatically consolidated the 6,362 motifs into a non-redundant set. We first derived clusters of motifs which were broadly similar. CWMs were trimmed by removing positions where the total contribution across nucleotides was less than 30% of the maximum total contribution among all positions17. PPMs were trimmed to the same positions as the trimmed CWMs, and converted to position frequency matrices (PFMs) by multiplying PPMs by the total number of seqlets associated with each motif. PFMs from all cell types in each organ were first leniently clustered, separately for positive and negative motifs, using the gimmemotifs package109 (v0.18.0, with command gimme cluster -t 0.8), which returns an average PFM for each motif cluster. Average PFMs from across organs were then subjected to a second round of clustering with gimmemotifs. Within each broad motif cluster, we then collapsed the constituent CWMs originating from individual cell types using the SimilarPatternCollapser functionality in TF-MoDISco, which merges together similar motifs using the same method it does for seqlets during initial motif discovery in step (1). This procedure resulted in 834 motifs. - 3) Finally, we performed motif quality control, and annotated and categorized each motif. For annotation, we used TOMTOM110 (v4.11.2) to compute similarities between the 834 de novo motif CWMs and a curated set of 5,193 known transcription factor binding site position weight matrices (PWMs) from JASPAR, CIS-BP, and other studies (obtained from https://resources.altius.org/~jvierstra/projects/motif-clustering-v2.1beta/), using the command tomtom -no-ssc -oc.–verbosity 1 -text -min-overlap 5 -mi 1 -dist pearson -evalue -thresh 10.0. For every motif, we manually inspected the most similar PWMs to assign a provisional label, typically at the transcription factor family level. We further collapsed highly similar motifs missed by the clustering approach, retaining the motif with the highest number of seqlets across cell types. We flagged 107 motifs which were low-complexity, noisy, or dominated by CG dinucleotides for exclusion. We categorized motifs as ‘base’ if they matched known PWMs; ‘base with flanks’ if they matched known PWMs but exhibited additional high-contribution nucleotides; and ‘homocomposite’ or ‘heterocomposite’ if the motifs clearly matched two similar or distinct known motifs, respectively. Motifs that did not resemble known PWMs were labelled as ‘unresolved’. After exclusion of low-quality motifs, this resulted in a set of 508 non-redundant motifs used for downstream analysis (Supplementary Table 6). Motif labels have the naming scheme ‘ | : / # ’. The unique ID is a value from 0 to 508 (one ID corresponding to the low-quality motifs was filtered out). The index is used to distinguish separate motifs which each match the same known motif, typically representing subtle variations in nucleotide preferences or flanks. Composite motifs used a similar naming scheme, with labels for constituent motifs separated by underscores. Motifs were also assigned a broad label (for example, ‘CTCF’ and ‘CTCFupstream’ motifs shared a broad ‘CTCF’ label), used throughout figures and analyses to aggregate results.

Identification of predictive motif instances

识别预测性基序实例

Para_01

1. 为了识别每种细胞类型中从头生成的基序的基因组实例，我们使用了Fi-NeMo软件包（https://github.com/kundajelab/finemo_gpu，v0.23，提交版本b81876d），该软件包以一种关注贡献值的方式进行基序扫描。
2. 简而言之，Fi-NeMo为每个峰区域拟合一个稀疏线性回归模型，以最小化每个区域中的贡献分数与基于修剪后的输入CWM加权组合重建的分数之间的差异。
3. 在系数矩阵中，某个位置上的非零系数表示该位置存在一个基序实例。
4. 对于每种细胞类型，Fi-NeMo的输出是一个覆盖所有峰区的BED文件，代表那些既在序列上匹配基序又具有相对较高绝对贡献分数的短区域。
5. 我们首先使用全部834个聚类后的CWM运行Fi-NeMo，参数设置为–alpha 0.8 –trim-threshold 0.3，其余使用默认参数（参数alpha现在称为lambda）。
6. 由于在线性建模方法中基序之间存在竞争关系，我们在排除复合基序后，使用缩减后的436个基序集合再次运行Fi-NeMo。
7. 我们进行了事后过滤以获得高质量的注释用于下游分析。
8. 为了评估特定基序实例调用的质量，Fi-NeMo计算输入CWM与通过平均该基序所有Fi-NeMo识别出的实例的贡献分数所得到的CWM之间的相关性。
9. 对于每种细胞类型，将两次Fi-NeMo运行的结果实例合并，并对实例CWM与输入CWM相关性低于0.9的基序进行标记。
10. 接着，过滤掉这些被标记基序的所有实例。
11. 最后，为了减少冗余，如果多个相同注释的实例重叠超过3个碱基对，则仅保留‘hit_correlation’值最高的实例，该实例表示其贡献分数与相应输入CWM具有最高相似性。
12. 此步骤最终生成了每种细胞类型的最终基序注释，包含其预测的基序实例。

Inference of nucleosome positions from ATAC modality

从ATAC模态推断核小体位置

Para_01

1. 我们使用 NucleoATAC34 从 SHARE-seq 的 ATAC 模态数据中推断核小体的位置和占据情况。
2. 简而言之，为了确定核小体占据度，NucleoATAC 将观察到的 ATAC 片段大小分布建模为核小体片段和无核小体片段的混合分布，并将特定位置上核小体片段的最大似然比例用作连续的占据信号。
3. 我们对原始代码进行了调整，使其以片段作为输入（https://github.com/sjessa/NucleoATAC，v0.4.1），并在用于训练 ChromBPNet 模型的相同峰值区域内，对每个聚类的伪批量片段文件运行 NucleoATAC 流程。
4. 后续分析使用了 nucleoatac occ 命令的输出结果：核小体二联体位置的调用用于基序实例的分析，而每个核苷酸位置上核小体片段的最大似然比例则用于将核小体占据情况可视化为基因组轨迹。

Annotation of motif instances

基序实例的注释

Para_01

1. 为了根据基因组特征注释预测的基序实例，如果实例位于转录起始位点（TSS）上游或下游2 kb范围内，则将其归类为位于启动子区；如果与外显子重叠，则归类为外显子区；如果与基因体重叠但不与外显子重叠，则归类为内含子区；其余情况归类为远端区域。
2. 基因组特征的定义基于Bioconductor TxDb.Hsapiens.UCSC.hg38.knownGene（v3.14.0）注释，该注释对应于GENCODE V38的UCSC knownGene轨道，并通过ArchR软件包中的createGeneAnnotation函数进行整合。
3. 基序实例到TSS的距离被计算为实例中心位置与上述注释中定义的TSS之间的距离。
4. 类似地，对于每种细胞类型，通过计算实例中心位置与由NucleoATAC推断出的最近核小体二联体或峰顶之间的距离，分别得到基序到核小体二联体或峰顶的距离。
5. 在分析中，我们从二联体或峰顶开始，在0到250 bp范围内以10 bp为区间统计了基序实例的数量。

Motif co-occurrence analysis

模体共现分析

Para_01

1. 为了识别共富集的基序，我们在每种细胞类型的 ChromBPNet 训练峰内对基序对的共现富集进行了检验。
2. 我们将分析限制在从头预测基序文库中的基础基序集合，并根据其广泛的注释对基序进行分组。
3. 对于每种细胞类型中的每一对基序，我们构建了一个 2 × 2 的列联表，记录包含两个、一个或不包含任何基序的峰的数量。
4. 然后我们进行单侧费舍尔精确检验，以计算富集的 P 值。
5. 为了校正所有基序对之间的多重假设检验，我们使用贝叶斯-霍赫伯格方法对 P 值进行了调整。
6. 我们将调整后 P 值小于 0.05 的基序对判定为具有显著共现。

In silico marginalizations to assess motif synergy

通过计算机模拟的边缘化分析来评估基序协同作用

Para_01

1. 我们使用训练好的ChromBPNet模型，按照先前的方法评估了基序协同作用。
2. 在这种计算机模拟的边缘化策略中，通过将一段短序列插入到一组背景性的、不可及的基因组区域库中（替换其中心核苷酸），并利用训练好的ChromBPNet模型进行前向传播，预测每个背景区域及其编辑后区域的可及性，从而量化该序列对染色质可及性的影响。
3. 通过对多个区域中编辑区域与背景区域之间预测的自然对数计数值差异取平均，我们估计出该序列的‘边际’效应。
4. 可以同时插入两个或多个序列，以估计它们在可及性上的联合效应。
5. 具体而言，对于两个基序A和B，我们将插入一个基序或同时插入两个基序后区域的预测对数计数值分别定义为 ({y"}{{\rm{A"}}})、({y"}{{\rm{B"}}}) 和 ({y"}{{\rm{j"}}})；未编辑背景区域的预测对数计数值为 ({y"}{0})。
6. 基序A和B的边际效应（以对数计数表示）分别为 ({\Delta }{{\rm{A"}}}={y"}{{\rm{A"}}}-{y"}{0}) 和 ({\Delta }{{\rm{B"}}}={y"}{{\rm{B"}}}-{y"}{0})，它们的联合效应为 ({\Delta }{{\rm{J"}}}={y"}{{\rm{J"}}}-{y"}_{0})。
7. 我们将基序A和B的独立效应定义为 ({\Delta }{{\rm{S"}}}={\Delta }{{\rm{A"}}}+{\Delta }_{{\rm{B"}}})，这对应于独立效应的对数相加模型，或在计数单位下的乘法模型。
8. 因此，协同效应可被定义为两个序列的联合效应 ({\Delta }{{\rm{J"}}}) 相对于其独立效应 ({\Delta }{{\rm{S"}}}) 的显著增强。

Para_02

1. 为了实施协同效应分析，我们使用了 tangermeme 软件包（https://github.com/jmschrei/tangermeme，v0.4.3）。
2. 我们首先对从头预测的复合基序集合进行筛选，使得每个复合基序均由一对独特的组成基序构成，最终获得138个复合基序。
3. 对于每个组成基序，我们识别出具有最多基序实例的基础从头预测基序，按上述方法修剪CWM，并将修剪后基序中每个位置上贡献分数最高的核苷酸定义为共识序列。
4. 所测试的基序及相关序列列于补充表7中。
5. 序列经过人工调整，以进一步去除无信息量的侧翼区域或更好地匹配复合基序，并进行去重处理，确保相同的序列对仅被测试一次。
6. 针对每种细胞类型，从GC含量匹配的背景区域中随机选取100个背景区域。
7. 对于两个相同基序的组合，存在三种独特的方向；而对于两个不同基序的组合，则存在四种独特的方向。
8. 我们将每种方向与基序间距的组合视为一种基序‘排列方式’。
9. 对于每个复合基序，将其组成基序序列A和B插入到100条背景序列中所有可能的方向和距离（从0到200 bp）位置，并使用全部五个ChromBPNet模型折叠（针对该复合基序预测实例最多的细胞类型）来预测背景序列和编辑序列的可及性，从而计算Δ_J（补充说明3）。
10. 类似地，将A和B单独插入以计算独立效应及独立效应之和Δ_A、Δ_B和Δ_S。
11. 每个序列和模型折叠的效果均被计算，并报告其平均效果（补充表7）。
12. 对于联合效应，报告的是在最优排列方式下（即具有最大平均效应Δ_J的方向与距离组合）基序对的效果。
13. 我们考虑每个基序对在其最优排列方式下的情况，并采用Wilcoxon符号秩检验来检验该排列下联合效应与独立效应之间的配对差异（Δ_J−Δ_S）是否显著大于0。
14. 多重检验校正采用Benjamini–Hochberg方法进行。
15. 将调整后P值小于0.001且（Δ_J−Δ_S）大于0.15的复合基序标注为具有协同作用。
16. 为了确认插入的序列确实驱动了预测的协同效应，我们如前所述使用DeepLIFT对编辑后的序列进行模型解释，并验证了预测可及性的序列是否与我们插入的序列一致。
17. 这一过程识别出18个复合基序，其预测效应由与插入序列不同的核苷酸驱动；对于这些基序，我们未进行协同作用分类，最终得到120个具有协同作用分类的基序（补充表7）。

Para_03

1. 为了定义具有语法偏好的协同基序（即协同作用仅限于或在特定结合位点排列下增强），对于每个复合基序，我们计算了所有排列下联合效应 Δ_J 的z分数。
2. 在任何排列中，若复合基序的效应超过均值四个标准差（z分数 > 4），则被标注为具有硬性语法。
3. 我们还考虑到基序可能存在较弱的长距离偏好，即软性语法。
4. 在组成基序相距20至150 bp的任何排列中，若复合基序的 (Δ_J - Δ_S) 大于0.15，则被标注为具有软性语法。

Para_04

1. 为了验证我们预测的特异性，我们对所有189种通过质量控制的细胞类型中的138个基序对，使用ChromBPNet模型重复进行了计算机模拟的边缘化实验（补充说明3）。
2. 对于每种细胞类型，将序列插入其相应的背景区域中，并在100条序列和5个模型折叠上计算平均效应。
3. 效应的计算方式如前所述，并采用Benjamini–Hochberg方法进行多重检验校正。
4. 为了评估预测协同效应的细胞类型特异性，我们针对部分复合基序，分析了在所有细胞类型中最优排列下的计算机模拟边缘化结果。

In silico motif ablations

计算机模拟的基序消融

Para_01

1. 类似地，对于选定的模体，我们使用tangermeme在计算机中进行了消融实验，并利用为一种细胞类型（Heart_c3，成纤维细胞）训练好的ChromBPNet模型。对于每个模体，我们根据Fi-NeMo命中相关性得分将其模体实例分为四分位数，并从每个四分位数中随机抽取250个模体实例（总共1000个序列）。为了在计算机中消融模体，我们通过将每个位置的碱基替换为具有最小绝对贡献值的碱基来“中和”每个实例，即最中性的碱基。这些假设的贡献分数是基于DeepLIFT方法计算得出的反事实估计值，用于评估在该序列背景下每个位置上其他三种替代碱基的贡献情况，是ChromBPNet模型解释工作流程中的一部分。

Ranking of motifs by prevalence and importance

根据普遍性和重要性对主题进行排序

Para_01

1. 为了评估基序的普遍性和重要性，我们聚焦于在细胞类型特异性TF-MoDISco基序发现步骤中每种细胞类型所学习到的基序，并仅考虑正向基序。
2. 为了计算基序的普遍性，针对每种细胞类型，我们获取了在该细胞类型中学到的基序集合，以及对应词典基序在该细胞类型中的基序实例数量。
3. 我们在每种细胞类型内根据基序实例的数量对基序进行排序，并通过将每个排名除以该细胞类型中的最高排名来对排名进行归一化，使得归一化后的排名落在[0, 1]范围内。
4. 类似地，为了计算基序的重要性，针对每种细胞类型，我们获取了在该细胞类型中学到的基序集合，并对相应修剪后的CWM在所有核苷酸和位置上的贡献分数进行求和。
5. 最后，为了比较在协同作用分析中定义的不同类别基序的普遍性和重要性，我们根据基序是否具有硬语法、软语法、无预测协同作用，或未被测试（意味着该基序不是复合基序，或如上所述在协同作用分析中被过滤掉）对其进行分组。
6. 同时具有硬语法和软语法的基序被归入硬语法基序组。
7. 然后，我们在每种细胞类型内计算每个组中基序的平均归一化排名。

Motif footprinting

模体足迹分析

Para_01

1. 对于选定的广泛性基序（CTCF、ZEB/SNAIL、NFY、NFY负向基序、YY1/2以及YY1/2负向基序），我们使用了基于深度学习得到的基序实例，并筛选出每种基序最常见的变异形式，然后根据Fi-NeMo命中相关性得分将这些实例分为四个四分位数。
2. 为了从SHARE-seq的ATAC模态数据中计算基序足迹的元图谱，我们使用了BPCells的footprint功能，该功能统计落在以基序实例为中心的500 bp窗口内每个位置上的插入片段（即片段的起始和终止）数量。
3. 随后，将每个位置上的插入片段计数除以窗口两侧外部各10%区域的平均插入片段计数，后者被用作局部背景可及性的衡量指标。

Enrichment of eQTL variants in motifs

在基序中富集eQTL变异

Para_01

1. 我们获取了组织特异性的GTEx v8 eQTL数据，并按器官来源将这些数据合并以匹配我们的胎儿器官，具体组织归类如下：肾上腺、脑（杏仁核、前扣带皮层BA24、尾状核基底节、小脑半球、小脑、大脑皮层、额叶皮层BA9、海马体、下丘脑、伏隔核基底节、壳核基底节、脊髓颈段C-1、黑质）、心脏（主动脉、冠状动脉、心耳附着物、左心室）、肺、肝、肌肉（胫动脉、骨骼肌）、皮肤（非日光暴露的耻骨上区、日光暴露的小腿下部）、脾、胃/食管（胃、胃食管交界处、食管黏膜、食管肌层）以及甲状腺。
2. 对于每个器官，我们在存在重复的情况下选择后验包含概率分数较高的变异位点，从而获得唯一的变异-基因对列表。

Para_02

1. 对于每个器官，使用每个变异-基因对的log2等位基因折叠变化效应大小（aFCs）来确定变异对基因表达方向的影响（上调或下调表达）。
2. 针对每个器官分别进行分析时，我们将每种细胞类型的全部基序实例合并，然后去除具有相同基序名称和基因组位置的重复条目。
3. 分别统计与匹配胎儿器官中阳性或阴性基序实例重叠的上调型和下调型变异的数量，以获得观察到的计数。
4. 随后将aFCs随机打乱并重新分配效应方向，再重复上述计数过程。
5. 我们对每个器官执行了100,000次打乱操作。
6. 富集得分定义为观察到的计数除以100,000次打乱后计数的平均值。
7. 富集P值计算为打乱后计数大于观察计数的打乱次数所占的比例。
8. 对P值采用Benjamini–Hochberg方法进行多重检验校正，如果某个基序的FDR小于0.05，并且其观察计数高于该类型eQTL变异打乱计数的95%置信区间，则认为该基序在上调或下调变异中显著富集。
9. 分别对阳性和阴性基序进行双侧Fisher精确检验，比较在上调或下调eQTL变异中显著富集与非显著富集的基序数量（补充表8）。

Enrichment of cell types with disease variants using g-chromVAR

使用g-chromVAR对具有疾病变异的细胞类型进行富集分析

Para_01

1. 我们使用g-chromVAR60来计算疾病变异在细胞类型特异性上的富集情况。
2. 由于CAUSALdb59中所有遗传变异均以hg19坐标表示，因此将来自胎儿细胞类型的ChromBPNet模型输入峰值集合从hg38提升至hg19坐标，并将其用作g-chromVAR的输入。
3. 所有坐标转换均使用rtracklayer（v1.54.0）R软件包112中实现的liftOver111工具，以及从UCSC基因组浏览器获取的hg38ToHg19.over.chain.gz或hg19ToHg38.over.chain.gz链文件完成。
4. 针对数据库中列出的每个变体，我们分别整理了CAUSALdb中全部13,710项研究的SuSiE得分。
5. CAUSALdb包含了一些在某些GWAS中缺失连锁不平衡信息的变异，这会影响精细定位过程中对因果变异的估计。
6. 这些变异在数据库中被赋予默认的PIP值-1，并从我们的分析中剔除。
7. 我们针对每种器官分别运行g-chromVAR（v0.3.2），以获得每项研究中可信变异对应的细胞类型富集的z分数和P值。
8. 我们使用SuSiE69生成的PIP值，并剔除那些研究中所有可信变异的PIP值总和小于5的研究，以排除统计效能不足的情况，最终保留了13,194项研究的结果。
9. 对所有P值采用Benjamini–Hochberg方法进行多重检验校正，若某细胞类型在特定研究中的变异富集FDR小于0.05，则认为其富集具有显著性。
10. 随后，我们手动提取与HDMA中各器官相关的疾病性状所对应的部分显著结果（相关研究列表见补充表10）。
11. 将潜在的因果变异（SuSiE PIP ≥ 0.8）从hg19坐标提升至hg38坐标，并与胎儿基序坐标进行重叠分析。
12. 对于数据集中每种胎儿细胞类型，我们随后手动识别出ENCODE中具有snATAC–seq数据的匹配成人组织和细胞类型，并从ENCODE70获取其snATAC–seq伪复制峰值集合（见补充表12）。
13. 接着，将与胎儿基序实例重叠的因果变异进一步与相应的成人峰值集合（峰值的−log10 P值 > 2）进行重叠分析（如数据可用）。
14. 仅当变异出现在少于两个匹配的成人峰值集合中时，才被视为‘仅胎儿特有’的命中结果。

Prediction of variant effect using ChromBPNet models

使用ChromBPNet模型预测变异效应

Para_01

1. 我们使用训练好的ChromBPNet模型预测并解释了特定非编码变异对染色质可及性的影响，方法与之前的工作一致。
2. 我们使用tangermeme软件包进行预测和模型解释。
3. 对于每个变异位点，我们使用相应细胞类型中模型训练所用的1,000个碱基对的峰区域，提取与该变异重叠峰区的参考基因组序列。
4. 将该序列（向两侧等长扩展至总长2,114个碱基对）输入到相应细胞类型的五个五折交叉验证训练的ChromBPNet模型中，以获得峰区内的预测可及性谱型和聚合的对数计数。
5. 对于每一轮交叉验证，为了将预测谱型的logits转换为可及性谱型，将谱型logits进行softmax处理，并通过指数化的预测对数计数进行缩放；同时使用DeepLIFT方法针对计数输出进行模型解释。
6. 接着，将效应等位基因替换到变异位置的序列中，并像参考序列一样获取预测结果和贡献分数。
7. 对于每个模型，我们计算了在变异位点中心、跨度为100个碱基对的窗口内，每个碱基预测读段计数差值的总和，并计算了跨折叠模型的平均效应得分。
8. 我们还计算了效应等位基因与非效应等位基因在峰区预测计数之间的log2倍数变化，其中log2倍数变化大于0表示预测效应等位基因会增加可及性。
9. 在图示中展示的是跨折叠模型的平均预测谱型和贡献分数。

Genome browser visualizations

基因组浏览器可视化

Para_01

1. 对于所有特定基因组位点的数据可视化，我们使用了 BPCells R 软件包以及我们代码仓库中包含的自定义脚本。
2. 所有基因组轨迹均已在线托管，可通过 WashU 基因组浏览器进行交互式可视化。

Statistics and reproducibility

统计与可重复性

Para_01

1. 没有使用统计方法预先确定样本量。
2. 分析中未排除任何数据。
3. 实验未进行随机化。
4. 研究者在实验过程中和结果评估时未对分组情况设盲。
5. 对于图中所有的箱线图，各组成部分含义如下：中心线表示中位数；箱子的上下边界表示上四分位数和下四分位数；须线表示距离不超过1.5倍四分位距的最小值和最大值。

Reporting summary

报告摘要

Para_01

1. 有关研究设计的更多信息，请参见与本文相关的《自然》系列期刊报告摘要。

Data availability

Para_01

1. 所有已处理的数据（包括片段文件、计数矩阵、细胞注释、全局caCRE注释、ChromBPNet模型、基序词典、基序实例以及基因组轨迹）均已存放在https://zenodo.org/communities/hdma。
2. 原始的匿名测序数据已存入序列读取档案数据库（Sequence Read Archive），登录号为BioProject PRJNA1402391。
3. 本研究产生的原始基因组数据的元信息已存放在Zenodo，DOI为https://doi.org/10.5281/zenodo.17259745。
4. 所有数据类型及其对应URL的说明见补充表14。
5. ENCODE v4 cCRE数据从公开数据库https://downloads.wenglab.org/Registry-V4/GRCh38-cCREs.bed下载。

Code availability

Para_01

1. 所有分析代码均可在 https://github.com/GreenleafLab/HDMA 获取，并已存档于 Zenodo，网址为 https://doi.org/10.5281/zenodo.17298234。