单细胞多组学研究揭示基因调控功能

资源利用单细胞多组学研究基因调控的功能图形摘要亮点d静息和免疫刺激的人血d数据集成识别跨细胞d图R揭示了染色质可及性d基因调控模型识别疾病相关的调控相互作用作者维奈角法比亚娜·卡尔萨杜阿尔特，严虎，...，Andrew S.放大图片创作者：Michael D.埃斯基罗斯特罗对应jason_buenrostro@harvard.edu简言之用于测量染色质可及性（ATAC-seq）和基因表达（RNA-seq）的单细胞方法正在迅速发展，但整合数据和推断基因调控关系的工具仍然有限。在这里，我们生成了静息和刺激的人血细胞的多组学数据，并提出了一个新的计算框架，用于构建基因调控网络（GRNs）。具体来说，我们描述了基因调控的功能推断（图1），一个工作流程：（1）将scATAC-seq与scRNA-seq配对，（2）将顺式调控元件连接到靶基因，(3)确定TF-基因关系。Kartha等人，2022，细胞基因组学2，1001662022年9月14日？2022作者。https://doi.org/10.1016/j.xgen.2022.100166会会~开放获取资源利用单细胞多组学研究基因调控的功能推断维奈角Kartha，1，2Fabiana M.Duarte，1，2Yan Hu，1，2Sai Ma，1，2Jennifer G.Chew，3Caleb A.Lareau，4Andrew Earl，1，2Zach D.Burkett，3Andrew S.Kohlway，3Ronald Lebofsky，3和Jason D.1，2，5，*1干细胞和再生生物学系，哈佛大学，剑桥，MA 02138，美国2基因调控观察站，麻省理工学院和哈佛大学布罗德研究所，剑桥，MA 02142，美国3数字生物学组，Bio-Rad，普莱森顿，CA 94588，美国4病理学系，斯坦福大学医学院，斯坦福，CA 94305，美国5引线触点* 通讯地址：https://doi.org/10.1016/j.xgen.2022.100166jason_buenrostro@harvard.edu总结细胞在对环境刺激做出反应时需要协调控制基因表达在这里，我们将scATAC-seq和单细胞RNA测序（scRNA-seq）应用于静息和刺激的人血细胞。总的来说，我们产生了91，000个单细胞谱，使我们能够探测跨细胞类型、刺激和时间的免疫应答的顺式调节景观。推进工具来整合多组学数据，我们开发了基因调控的功能推断（FigR），这是一种计算配对scA-TAC-seq与scRNA-seq细胞的框架，将远端顺式调控元件连接到基因，并推断基因调控网络（GRN）以识别候选转录因子（TF）调控因子。利用这些配对的多组学数据，我们定义了免疫刺激的调节染色质（DORC）的结构域，并发现细胞以分钟的时间尺度改变染色质可及性和基因表达刺激GRN的构建在疾病相关的DORC处洗脱TF活性总的来说，FigR能够阐明单细胞数据中的调控相互作用，为了解组织内细胞的功能提供了新的机会。介绍真核细胞已经进化出精致的控制能力，可以持续感知和响应外部环境的提示。1-3-6这种高度环境特异性和通常异质性的基因激活促进适当的抗病毒或炎症反应，包括免疫的标志之一随着时间的推移，我们对免疫的理解也在不断发展;例如，研究表明，染色质可能会引发细胞的免疫反应，7，8导致分泌状态，9或通过产生关键信号分子进一步协调周围细胞的激活10单细胞基因组学方法极大地促进了我们对免疫细胞的细胞多样性的理解。例如，表征小鼠10和人14免疫细胞中的时间和刺激依赖性转录特征的单细胞RNA测序（scRNA-seq）已经鉴定了随时间活化或抑制的不同转录程序，并突出了响应免疫刺激剂的15同时，一些先前的研究已经应用染色质可及性和基因表达测定来定义顺式调节图谱跨越静息12、16、17和刺激14、16免疫细胞类型。最近，2019冠状病毒病（COVID-19）大流行促使人们使用单细胞ATAC-seq和RNA-seq工具来表征对感染的这些不同的努力试图阐明免疫细胞功能的表观遗传控制;即，定义细胞内基因调控网络（GRN）的细胞回路。虽然这些努力已经导致了巨大的见解的转录控制的免疫细胞，这些研究是有限的，由现有的计算工具建模单细胞之间的基因调控动力学。从单细胞数据构建GRNs的进展20，21促进了新的机会，以揭示刺激后细胞功能和适应机制。然而，仅利用共表达的大多数方法20、22、23在以下能力方面受到限制：（1）确定关键的顺式调节元件和（2）阐明主TF调节因子对基因表达的功能。大量先前的工作已经证明，表观基因组学数据可以大大提高功能GRNs的确定。24-在另一个实例中，使用共表达的scRNA-seqGRN方法依赖于以下假设：CellGenomics 2，100166，September 14，2022 <$2022作者。1这是一篇基于CC BY-NC-ND许可证的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。会开放获取资源2Cell Genomics2，100166，2022TF激活基因;然而，广泛的功能实验表明，TF可以沉默染色质以抑制靶基因。36-在这里，我们创建了一个用于构建免疫细胞GRNs的示例数据集。为此，我们将多种刺激剂与染色质可及性和基因表达单细胞分析相结合，以表征和评估与人外周血单核细胞（PBMC）中免疫细胞刺激相关的顺式调节景观的动力学然后，我们建立了基因调控的功能推断（图R），这是一种用于独立或同时分析单细胞ATAC-seq（scATAC-seq）和scRNA-seq的可遗传方法，（1）计算配对scATAC-seq和scRNA-seq数据集（需要时），（2）推断顺式调控相互作用，以及（3）定义TF-基因GRN。利用这些综合数据，我们确定染色质可及性的变化预示着单核细胞免疫刺激后基因表达的变化。最后，我们强调了如何使用这种方法来识别关键的转录因子及其与靶基因的关系，包括刺激反应和疾病相关的调控蛋白结构域（DORC）。我们的工作强调了血液刺激结合高通量单细胞多组学的使用，以及使用FigR作为模型来推导免疫细胞激活所需的关键转录调控模块来开发增强子GRNs的进展。结果静息和刺激PBMC的组合高通量单细胞表观基因组为了表征与宿主对人血液中刺激物的响应相关的染色质可及性和转录景观，我们在刺激物暴露的不同时间点对静息和刺激的人PBMC进行基于液滴的scATAC-seq和scRNA-seq（图1A;STAR方法）。具体而言，将来自健康供体的细胞（n = 3或4;表S1）暴露于已知可引发抗病毒样或核心炎症反应的刺激物1或6 h，包括脂多糖（LPS;细菌细胞膜的一种组分）、佛波肉豆蔻酸酯（PMA）加离子霉素（一种激活核因子κB [NF-κB]信号传导的强效酯）39或干扰素γ（IFN-γ;一种内源性产生的免疫调节细胞因子），同时在单细胞分析之前的每个时间点使用DMSO对照（STAR方法）。选择这些刺激物是因为它们已显示出诱导具有独特转录动力学的不同时间和10，14，39-总的来说，我们产生了超过150亿个读段，产生了高覆盖度的单细胞调控图谱，其包含跨越所有条件的67，581个scATAC-seq和23，754个scRNA-seq细胞（图1B），其中每个细胞平均8，865.2（±SD = 4，837）个比对的独特核片段，并且对于scATAC-seq概况的细胞（图S1 A和S1 B），峰中读段的平均分数（FRiP）为0.6（±SD = 0.05），并且对于scRNA-seq概况的细胞（图S1C和S1 D;STAR方法），平均3，021个UMI（±SD =425.77）。聚集scATAC-seq和scRNA-seq细胞（STAR方法）产生离散的细胞簇，主要代表单核细胞、T（CD 4/CD 8）和B淋巴细胞以及自然杀伤（NK）细胞，每个簇和条件下来自所有供体的细胞均匀分布（图1C、1D和S1E-S1重要的是，这些更广泛的簇中的每一个包括通过刺激条件的细胞的亚簇（图IE和S1H）。为了正式注释scRNA-seq细胞的细胞类型，我们首先使用先前描述的计算方法（在Seurat中）42跨批次（此处定义为每个处理条件）比对细胞，使得能够跨条件共聚类和注释scRNA-seq细胞。使用该方法的细胞聚类产生不同的分组（图S2A和S2 B），其富集细胞类型和刺激特异性基因表达标记物并用于注释细胞类型（图S2C）。骨髓细胞在基因启动子（scATAC-seq）和基因表达水平（scRNA-seq）周围的可及性峰的检查证实了刺激特异性和时间特异性变化（图1F、1G和S2D）。重要的是，所有主要细胞类型在所用处理条件下以相对均匀的比例捕获（图S2E），从而能够实现独立测定的染色质可及性和下游基因表达谱的多组学整合。精确整合单细胞染色质可及性和基因表达谱我们推断，来自配对背景的数据可能能够确定GRNs，促进对免疫细胞刺激的关键调控过程的解释支持scATAC和scRNA-seq数据29、42、43的整合的当前框架依赖于识别“锚"细胞，即在共同的低维空间中代表共享生物状态的细胞，然后从一个数据集中找到另一个数据集中的代表性细胞。虽然对于匹配相应细胞类型的细胞（即，注释级配对），这些方法通常（1）导致高的一对多细胞条形码匹配率，导致下游的细胞使用总体较低，或（2）不能充分解决数据集之间的细胞类型不平衡。为了解决这一挑战，我们开发了一种方法（scOpt-Match），其使用受约束的最佳细胞映射方法鉴定scATAC-seq和scRNA-seq数据之间的细胞对（图2A）。对于这种方法，我们首先使用典型相关分析（CCA）创建scATAC-seq和scRNA-seq细胞的共享共嵌入，类似于先前在Seurat中描述的功能。42接下来，我们通过首先对整个单元空间进行子聚类并构建单元k-最近邻来Cell Genomics2，100166，2022年9月14日3会开放获取资源图1.静息和刺激的人血细胞的高通量单细胞表观基因组和转录谱分析(A) 刺激实验的示意性高亮用DMSO对照、脂多糖（LPS）、干扰素γ（IFN-γ）或佛波醇肉豆蔻酸酯乙酸酯（PMA）加离子霉素刺激人外周血单个核细胞（PBMCs）1或6 h，6 h处理条件下有或无高尔基体抑制剂（GI）然后将细胞分裂并使用scATAC-seq和scRNA-seq对所考虑的每个条件和时间点进行分析。(B) 通过scATAC和scRNA-seq的质量控制过滤的每个条件下分析的细胞总数(C) 基于潜在语义索引（LSI）降维的scATAC-seq细胞的均匀流形近似和投影（UMAP），其中细胞由处理条件着色。(D) 基于主成分分析（PCA）降维的scRNA-seq细胞的UMAP，其中细胞通过处理条件着色(E) scATAC-seq细胞（上）和scRNA-seq细胞（下）的UMAP，突出显示对照（6小时）和PMA（1和6小时）条件下的各个条件。(F) scATAC-seq单核细胞在基因IFITM 3和HES 4周围的聚集可及性概况。(G) 基于scRNA-seq单核细胞的刺激特异性基因标志物的插补scRNA-seq计数的单细胞表达水平的分布，如（F）中所示。(kNN)共包埋空间中ATAC和RNA细胞之间的图，在给定的kNN子图内对来自两种测定的细胞进行采样（STAR方法）。在下采样以匹配测定之间的细胞数量时（即，scATAC或scRNA）在给定的亚基中，图中，使用受约束的全局匹配算法44将ATAC-RNA细胞之间的子图测地线距离用作成本函数来配对细胞类似于旅行推销员问题，这确保了产生的ATAC-RNA细胞4Cell Genomics2，100166，2022会开放获取资源图2.基于稀疏kNN的ATAC-RNA细胞配对允许scATAC-seq和scRNA-seq数据的最佳配对和整合(A) 示意性突出显示scOptMatch(B) 描绘来自人PBMC的特定免疫细胞类型的实验珠富集的示意图(C) 当使用greedy（左）与scOptMatch方法进行PBMC分离物数据集配对时配对的RNA细胞条形码的实例数的分布。(D) 使用两种不同配对策略配对的总scATAC和scRNA-seq细胞的百分比。(E) PBMC分离物细胞类型之间scATAC-scRNA-seq配对的准确性热图，通过与相应scRNA-seq细胞类型正确配对的scATAC-seq细胞的百分比着色。(F) scRNA-seq刺激的细胞的UMAP显示在图ID中，其中细胞在刺激条件下对齐以实现细胞类型注释，通过注释的细胞类型着色。（图例接下页）Cell Genomics2，100166，2022年9月14日5会开放获取资源对于所有可能对的组合中的总测地距离，对被最小化。重要的是，只有在一定距离内的ATAC-RNA细胞（测地线kNN）被考虑用于配对作为先验，相对于考虑所有可能的配对（STAR方法），进一步加速了计算时间为了创建基准scOptMatch的参考数据集，我们分离了PBMC内的细胞类型，并对scRNA-seq和scATAC-seq进行了分析（以单独的方式）。13完整的数据反映了对应于五个PBMC亚群的scATAC-seq（n = 17，920个细胞）和scRNA-seq（n = 8，089个细胞）数据（图2B;STAR方法）。使用这些数据，我们使用（1）上述最佳匹配（scOptMatch）或（2）“贪婪”最佳匹配方法（在CCA空间中为每个ATAC细胞选择最接近的RNA细胞）确定ATAC-RNA细胞对。如所预期的，我们发现scOptMatch导致在所 有 细 胞 中 从 两 个 数 据 集 配 对 的 细 胞 的 数 量 显 著 更 大（92.06%scATAC和98.4%scRNA;图 2C和S3F），这是与贪婪方法（22.2%scRNA）相比更少的ATAC-RNA细胞多重映射实例（图2D和S3G）的结果。重要的是，scOpt-Match方法还能准确地 绘 制 相 同 参 比 细 胞 类 型 的 细 胞 图 （ 图 2E ） . 为 了 确 认scOpMatch的配对性能，我们将其应用于先前生成的初级运动皮层细胞（n = 84，178个细胞）的SNARE-Seq 2数据，45代表具有用于实验配对的染色质可及性和每个细胞的RNA表达谱的真实标记的较不离散的细胞群体（图S3H）。实施与未配对数据相同的整合工作流程导致具有共享细胞类型簇注释的细胞之间的平均ATAC-RNA作图率为90.05%（±22.3 SD）（n = 10个簇组;图S3I ） ， 尽 管 对 于 完 全 相 同 的 细 胞 几 乎 没 有 ATAC-RNA 匹 配（<1%），其中错配主要发生在非常罕见的细胞类型（例如，L5ET-2，n = 29个细胞）或密切相关的簇注释（例如，Sngc/Vip细胞）。因此，scOpt-Match可以使用独立生成的或多模式的谱在多个测定中配对相似细胞类型的细胞。如果scOptMatch能够在测定之间实现具有相似注释的细胞的近似一对一配对，我们推断它可以促进从匹配的细胞背景产生的scATAC-seq和scRNA-seq数据的整合分析。因此，我们试图将其应用于配对我们的刺激多组学数据集。使用scOptMatch将scATAC和scRNA-seq细胞按条件配对（图S3重要的是，这种细胞配对通过简单地使用从scRNA-seq基因表达标记物定义的注释进 一 步 实 现 了 scATAC-seq 的 细 胞 类 型 注 释 （ 图 2F 、 2G 和S3M）。通过细胞类型和条件聚集单细胞并过滤以获得足够的计数，得到139个假批量（平均总共194万RNA和230万ATAC聚集体计数）。利用这种高深度的资源，我们发现，染色质的可及性在远端峰是高度细胞类型的spe-特异性，甚至比基因表达更特异，而启动子可及性在细胞类型和刺激条件下相对不变（图2H），验证了先前的报道。17，46总的来说，这些数据的高质量和远端染色质可及性的精致细胞类型特异性促使进一步分析GRN潜在的刺激反应。我们推断，这种用于细胞配对的scOptMatch方法，使得scATAC与scRNA谱能够近似均匀地配对，将建立整合的数据集，并且可以用于下游分析，类似于在同一细胞内伴随的可接近性和RNA表达谱。使用整合的单细胞数据鉴定跨刺激的远端峰-基因相互作用我们接下来试图将顺式调节峰的变化与基因表达相关联，作为优先考虑作为免疫应答GRN的一部分的特征的手段。为此，我们将计算方法作为FigR框架的一部分，以确定显著的远端峰与基因表达的相互作用，如先前在多模态数据上进行的那样。28具体而言，我们使用计算配对的细胞（每次测定n = 62，219个细胞）将每个基因转录起始位点（TSS）周围固定窗口（100 kb）内发现的峰的可及性与该基因的表达相关，通过这种方式，我们在全基因组范围内鉴定了总共34，370个独特的染色质可及性峰，显示出与基因表达的显著相关性（排列p % 0.05），跨越总共11，304个基因。基于显著相关峰的总数对基因进行优先排序，我们鉴定了与高密度峰-基因相互作用相关的基因子集，我们最近将其描述为DORC28（图3B; nR 7显著峰-基因关联，n = 1，128个基因，n = 12，583个峰;数据S1）。DORC相关基因的列表包括许多已知的与先天性和适应性免疫应答途径相关的免疫应答的介质，10、40、47、48，这也通过基因集富集分析（GSEA）证实（图S4A;数据S1）。值得注意的是，在这些基因中，我们看到大部分远端顺式调节缔合（距离基因TSS>5 kb;图 3C、S4B和S4 C）。通过使用每个DORC的总相关峰可及性信号对细胞进行评分（称为DORC可及性评分），我们确定了单个细胞（图3D、S4D和S4 E）或每个DORC、刺激条件和细胞类型的假批量（图3E）中染色质可及性和基因表达在与匹配的对照条件（DMSO对照）比较后，我们观察到PMA处理对DORC可及性和表达的最大影响，如在大多数细胞类型中所见，以及IFN γ或LPS处理的更适度的影响，如主要在单核细胞中所见（图3F）。值得注意的是，我们发现，刺激诱导细胞转录组的较大变化，(G) 未对齐的scRNA-seq细胞（图1D中所示）的UMAP由注释的细胞类型（左）着色，并且scATAC-seq刺激的细胞（图1C中所示）由成对的scRNA-seq细胞注释（右）着色，使得能够对刺激的scATAC和scRNA-seq分析的细胞进行下游数据整合。(H) 与基因启动子（左）、启动子远端（中）和成对基因表达（右）相关的聚合单细胞染色质可及性谱的成对皮尔逊相关性，按细胞类型和条件聚合。6Cell Genomics2，100166，2022会开放获取资源图3.整合的多组学分析确定了与单细胞刺激反应相关的关键调控模块(A) 使用计算配对的scATAC-seq和scRNA-seq刺激数据集鉴定显著的染色质可及性峰-基因关联的顺式调节分析的示意图(B) 基于所有细胞类型和刺激条件下显著基因-峰相关性的数量的最高命中(C) 使用（A）中概述的方法确定的突出显示DORCTRAF 1(D) TRAF1的DORC可及性评分（左）和配对RNA表达（右）的UMAP。(E) 所有DORC中每种条件每种细胞类型的聚合DORC可及性评分和细胞RNA表达的成对Pearson相关性，通过DORC评分相关性使用分层聚类进行聚类。(F) 基于每种刺激条件与其相应对照条件下DORC间总评分的Pearson相关系数显示的总体DORC可及性（上图）和基因表达（下图）变化，按条件/细胞类型注释显示。(G) 热图显示了在各种条件和细胞类型中前10个差异DORC的联合的单细胞DORC可及性的平均差异（n = 53个基因）。细胞类型颜色条表示在该试验的所有条件下具有最显著变化的细胞组。与染色质可接近性相比，细胞类型的免疫调节作用显著，但细胞类型同时改变了染色质和表达，从而诱导免疫基因的激活（图3E和3F）。有趣的是，我们还发现GI的添加强烈减弱了对PMA（CD8、NK和B细胞）和LPS（CD8）的免疫应答，这可能是抑制旁分泌信号传导的结果，并且响应于IFN（单核细胞），这可能是抑制自分泌信号传导的结果。令人惊讶的是，只有少数细胞类型，包括B和CD8 T淋巴细胞，在PMA暴露6小时后，当同时用GI处理时，表现出这种可接近性和基因表达变化的抑制反应，表明在大多数情况下，DORC是内在调节的。DORC间的单细胞差异检测确定了许多免疫应答（图3G和S4F;数据S2和S3）。这包括共享的LPS和IFN诱导的基因（MX 1，IFIT 3，OAS 3和OASL）和PMA诱导的与细胞凋亡和存活相关的基因（NR 4A 1/2/3，EGR 1，REL和TRAF1）。有趣的是，我们还观察到主要配体编码基因（IL1A，IL1B和CCL3）和免疫抑制剂（CD274[也称为PDL 1]，NFKBIA和TNIP 1）在这些顶级差异DORC中。值得注意的是，我们的顺式调节分析回收了DORC，其中大部分（~79%）包括先前注释为在不同细胞背景中与超级增强子区域连接的基因31（STAR方法;图S4 G和S4 H），其余（n = 238个基因）包括几种刺激应答基因（IFIT 1、MX 1、OAS 13、IL 13、IL 3RA和IL 27 RA）和细胞类型标志物（CD 14、NKG 7、GZMK和CD 8B）。我们的方法来Cell Genomics2，100166，2022年9月14日7会开放获取资源图4.染色质和基因表达动力学与刺激反应时间的关系(A) scATAC细胞的UMAP由每个刺激条件下估计的NN刺激（刺激）时间着色。(B) 突出显示针对三种不同刺激条件确定的单核细胞约束的DORC基因的重叠的UpSet图(C) 热图突出显示了DORC基因（n = 38）的平滑归一化DORC可及性、RNA表达和残留（DORC-RNA）水平，这些基因被鉴定为与对照（1 h）和刺激（1 h/6h）单核细胞（n = 1，776个细胞）中的LPS NN刺激时间相关。(D) 对照（1 h）和刺激（1 h/6 h）单核细胞（PMA +对照n = 2，002个细胞，IFNγ+对照n = 2，601个细胞）的染色质（DORC）与DORCFOSB（左）和IFIT 3（右）相对于平滑PMA和LPS NN刺激时间的基因表达（RNA）虚线表示与从DORC可及性或RNA表达水平的滑动平均值获得的值的LOESS拟合（n = 100个细胞/滑动窗口箱）。颜色条指示每个箱中最常见的（模式）细胞条件。(E) 与（D）相同，但TRAF 1相对于LPS刺激和对照（1 h）单核细胞。(F) （D）中所示的对照和LPS刺激的单核细胞中与TRAF1表达相关的单个顺式调节元件的平滑可及性评分，按LPS NN刺激时间排序。鉴定DORC揭示了广泛染色质控制下的基因受刺激的细胞的特征是染色质可及性景观的早期变化，这些变化引发基因表达之前，我们使用多模态数据来表明DORC可接近性预示着沿着发育轨迹的基因表达，并且该活性预测细胞状态转换。为此，我们试图使用配对的多组学数据来确定细胞是否通过其染色质可及性状态引发免疫反应。推导轨迹伪时间的方法通常需要定义单个根单元类型。因为我们确定了18种离散的细胞类型，排除了伪时间的使用，我们试图利用一种替代方法来定义轨迹。为此，我们计算了每次治疗的细胞最近邻（NN）刺激时间，其表示基于实验治疗标签的刺激暴露时间的加权平均值简而言之，我们对每个细胞取细胞NN（k= 50），并计算对照、1-h刺激和6-h刺激细胞的平均邻域，分别分配权重0、1和2（图4A和S5A;STAR方法）。这种连续的时间测量使我们能够研究沿着刺激轨迹的染色质可及性和基因表达动力学（图S5B）。使用这些刺激时间定义，我们试图确定染色质可及性是否在基因表达之前激活以“对于该分析，我们选择关注单核细胞群体，因为它响应于我们的炎症因子而被直接激活，如先前文献和我们对GI的观察所述（图3F和3G），以评估染色质和基因表达。8Cell Genomics2，100166，2022会开放获取资源-表达动力学与刺激时间有关将我们的峰-基因相关性方法严格限制为对照1 h、刺激1 h和刺激6 h单核细胞，我们鉴定了一组与单核细胞的LPS（n = 38个基因）、IFN γ（n = 33个基因）或PMA（n = 39个基因）刺激相关的DORC基因。 这些DORC包括在骨髓细胞中刺激时诱导的已知表达标志物。有趣的是，我们还发现这些单核细胞特异性DORC的一小部分在多种刺激中共享（图4B）。 通过对每个处理的响应于NN刺激（刺激）时间的单细胞DORC可接近性和RNA水平求平均（STAR方法），我们可视化了沿着对照（0 h）至6 h刺激时间轴的染色质可接近性和基因表达的变化（图4C、S5C和S5D）。计算染色质与RNA的差异（残差），我们主要观察到染色质可及性变化先于早期时间点的表达变化（高残差）在稍后的时间点，我们发现残差较低，反映了免疫刺激后RNA的积累这些观察结果分别由基因FOSB和IFIT 3用LPS和PMA处理定型（图4C、4D和S5E）。这些变化发生在相对较快的时间尺度上;例如，FOSB的染色质变化是在60分钟时间点内发生的早期事件。值得注意的是，DORC可接近性的这些变化由单独的顺式调节元件构成，其中一些元件变得快速可接近（即，启动子），而其他的则沿着刺激时间轴缓慢地变得可接近（一些远端调节元件），如LPS应答基因肿瘤坏死因子（TNF）受体相关因子1（TRAF1）所突出显示的（图4E和4F）。相反，我们注意到一些例外，包括PMA应答性热休克蛋白编码基因HSP 90 AA 1和HSPH 1，与DORC可及性的相应变化相比，其表现出早期表达增益（图S5F）。我们证明，使用计算配对的多组学数据，在与刺激样细胞状态相关的基因表达之前检测染色质可及性激活的能力。一种确定DORC活性候选TF调节因子的计算方法在FigR的核心，我们开发了一种计算方法来使用多组学数据定义免疫应答的GRN。在这个阶段，FigR使用配对的scATAC-seq和scRNA-seq数据，并具体测试TF基序在预定的顺式调节元件（即，DORC）以及TF表达与给定DORC基因的总体可及性水平（DORC评分）的相关性，以推断可能的TF激活子和阻遏子（图5A）。首先，对于给定的DORC基因，我们基于其DORC可及性kNN确定DORC顺式调节元件库这假设跨整个细胞空间共变的DORC然后，我们使用TF基序的参考数据库中相对于GC含量和全局峰可达性匹配的背景组的基序匹配频率对TF基序富集的显著性（Z检验）进行统计学检验。同时，我们计算TF RNA表达水平和DORC可及性评分之间的Spearman相关系数。最后，确定激活子和阻遏子，我们结合相对于所有TF的相对基序富集（Z检验P）和RNA表达抑制（Z检验P）的显著性估计，计算我们称为“调节概率”（-log 10尺度上的“调节得分”）的有符号概率得分，代表基序富集和RNA相关TF的交叉。为了能够使用这种方法来阐述新的调节剂，我们策划了一组扩展的独特的人（n = 1，141）和小鼠（n = 890）TF结合序列基序，其扩展了先前建立的数据库49（STAR方法）。为了证明FigR GRN方法的实用性，我们将其应用于配对刺激scATAC-seq和scRNA-seq数据，以揭示刺激反应的关键调节因子为此，我们首先测试所有刺激响应性DORC基因（n = 1，128）和参考TF基序（图5B;数据S4）。使用调节分数阈值（abs（调节分数）R1）过滤TF-DORC关联，然后我们可以查询给定DORC的推定TF调节子（图5C、5D、S6A和S6 B）以及可能由特定TF驱动的DORC集合（图S6C）。例如，FigR鉴定了MX1的已知活化剂，包括TF的IRF家族：IRF3、IRF7、IRF9和STAT2，所有这些都属于IFN信号传导途径。50我们通常根据所有DORC的平均调节评分（图5E）或通过正调节和负调节的DORC的分数（图56D）区分TF激活剂和TF抑制剂例如，我们看到SPI1（PU.1）、BACH1和BCL 11A作为顶级转录激活因子，其作用已在之前描述17、51、52，BCL 11B作为顶级转录阻遏物（图5E、S6C和S6D）。重要的是，BCL11B已被证明是T细胞成熟的关键阻遏物。图53[61]我们的方法估计，与DORC相关的35.6%的TF（绝对（调节评分）R1）具有抑制性关联（所有靶DORC的平均调节评分为0），这与以前的工作一致，反映了大部分TF作为阻遏物发挥作用的理解。最后，证明了FigR的GRN推断的无偏性质和扩展的TF基序数据库的利用的价值，我们鉴定了支配免疫细胞功能的新的我们强调ZEB2，它被认为是CD8+ T细胞功能的抑制因子。图62验证了其表达活性并鉴定了下游靶DORC（n =132;调节评分1.5），包括IL 7 R和TCF 7，其与T细胞中的免疫记忆相关。63我们还鉴定了激活剂ZNF 467，这是一种全新的调节剂，它调节86个基因（调节评分> 1.5）。在ZNF 467的直接靶点中，我们发现了CD14，它是单核细胞的关键标志物，在单核细胞活化和分化中对免疫功能至关重要重要的是，我们观察到这些TF的表达是刺激特异性和细胞类型特异性的（图S6F），表明FigR能够利用单细胞数据中的协方差来确定上下文特异性为了突出这种方法的广泛通用性，我们将我们的方法应用于来自替代组织的现有多模态数据，包括来自人类皮层的SNARE-Seq 2数据和来自小鼠皮肤组织的SHARE-seq数据SNARE-Seq2脑数据的顺式调节相关性（图S3H）鉴定了与DORC相关的基因子集（n = 432;图S7 A），包括Cell Genomics2，100166，2022年9月14日9会开放获取资源图5.FigR基因调控网络（GRN）工作流程的设计和应用(A) FigR GRN工作流程示意图。(B) 散点图显示所有DORC与TF的关联，用带符号的监管评分着色(C) MX1的候选TF调节子。突出显示的点是具有abs（调节评分）R1（-log 10量表）的TF，所有其他TF以灰色显示。(D) （C）中突出显示的TF的调节评分（签名，-log 10量表）(E) 所有DORC（n = 1，128）中每个TF（n = 870）的平均调节评分（带符号，-log 10量表），突出显示选择的TF激活剂（右偏）与TF表达剂（左偏）。(F) 涉及GWAS变体的DORC的所有显著TF-DORC富集的DORC调节评分的热图（左）（abs（调节评分）R1.5; n = 89个TF，n = 73个DORC）。还显示了所考虑的所有疾病中每个DORC的相应最小GWASP(G) SLE GWAS SNP相关DORC（橙色节点）及其相关TF（灰色节点）的TF-DORC网络可视化（F）。边缘按带符号的法规分数缩放10Cell Genomics2，100166，2022会开放获取资源已知的神经元（兴奋性：FEZF 2、RORB;抑制性：PVALB、LAMP 5、GAD 1）和非神经元（小胶质细胞：PAX 6、SLC 1A3;少突胶质细胞：MBP、FA 2 H）分化的标志物（图S7B和S7C）。使用这些皮质特异性DORC应用FigR接下来，使用预定义的峰-基因关联，我们的GRN推断方法恢复了我们先前发现与鼠毛囊分化相关的DORC的TF调节子（图S7H）。28这包括激活子Lef1、Hoxc13和Grhl1以及阻遏子Tcf12和Pou2f3（图S7I）。我们确定激活子Dlx365和阻遏子Zebl和Barx 266是顶级TF调节子（图S7J和S7 K）。为了评估FigR的在基因子集周围引入顺式调节元件的方法DORC）和使用scRNA-seq共表达来衍生GRN的替代方法，我们进行了使用我们的方法和SCENIC20确定的优先调节剂的比较，当应用于刺激PBMC数据时（STAR方法）。比较任一方法中由给定TF正调控的DORC基因的数量，我们观察到TF的子集，与SCENIC相比，通过FigR确定的靶点数量更大，反之亦然（图S8A）。有趣的是，查询GM 12878染色质免疫沉淀测序（ChIP-seq）数据以获得来自任一类别的顶级TF表明，SCENIC优先级的TF调节剂倾向于结合启动子区域，而图R优先级的那些与映射到远端增强子的调节元件重叠，多于基线所预期的（图S8B）。这表明FigR确实优先考虑与增强子相关性更高的TF。因此，我们表明，FigR可以利用不同的单细胞技术的实验配对的多模式数据，以获得GRN使用经验推导的峰-基因和TF-峰基序关联，以达到候选TF调节。接下来，我们想确定从FigR推断的刺激反应GRN是否可用于揭示疾病相关遗传变异及其非编码调控元件的调控机制。为了揭示疾病相关的细胞状态，我们对与GWAS SNP重叠峰相关的单细胞可及性进行评分。(GWAS，第10- 7页;图S8 C）。 我们观察到所测试的不同炎性疾病的染色质可及性的刺激特异性和细胞类型特异性富集（图S8 D），验证了先前的工作16，67，表明免疫刺激揭示了富集疾病GWAS变体的调节元件。例如，对于系统性红斑狼疮（SLE），我们观察到LPS和IFN γ刺激的B淋巴细胞和单核细胞中GWAS相关可及性的富集升高，对于过敏，我们观察到IFN γ和PMA刺激的CD 4/CD 8淋巴细胞中GWAS相关可及性的富集升高（图S8 E）。我们发现我们的免疫刺激揭示了与自身免疫和相关遗传变异相关的细胞状态及其相应的染色质可及性概况。接下来，我们推断，我们基于GRN的分析可以识别疾病相关遗传变异的相关机制。例如，已知调节因子NF-κB跨细胞类型起作用以促进炎性基因表达。68的确，我们发现NF-κB驱动大部分GWAS变体相关DORC的活性（图5F）。将该分析扩展到所有DORC（n = 77），我们发现了89个推定的TF驱动因子（abs（调节评分）R1.5），揭示了跨越一种或多种疾病的谱系年龄决定性以及刺激响应性TF的组合（图5F）。对SLE特异性DORC的子集（n = 15个DORC，n = 48个相关TF）的更仔细检查揭示了关键的调节关联，包括先前确定的SLE基因：BLK、IRF5、IRF8和NCF 2（图5G）。我们的方法可以优先考虑DORC及其推定的TF调节因子，以剖析与各种自身免疫性疾病有关的调节程序。我们包括推断的刺激-响应PBMC GRN ， 其 可 以 通 过 RShiny 网 络 应 用 程 序 （ https ：//buenrostrolab. shinyapps.io/stimFigR/）上提供。讨论在这里，我们生成了人体血液中免疫刺激的调控图谱。这项工作是通过高覆盖率的单细胞数据和支持多组学数据整合、顺式调节分析和基于单细胞概况构建增强子感知的GRN的新计算在这一过程中，我们克服了三个关键挑战：(1) 我们实现了一种方法，以更好地计算配对单细胞，（2）我们将远端顺式调节峰与靶基因相关联，（3）我们将TF与靶基因相关 联 。 重 要 的 是 ， FigR 使 用 独 立 或 伴 随 产 生 的 单 细 胞ATAC/RNA数据来推断GRN的能力将广泛地实现跨广泛范围的scATAC-seq和相关多组学技术的DORC-GRN分析与仅使用共表达或共可及性的静态测量的现有方法不同，使用我们提出的FigR框架的GRN构建通过跨单个细胞的这些特征的相关性来简化染色质和RNA动力学，提供了鉴定跨越细胞状态的基因调控关系的手段。为此，我们利用经验统计方法来计算TF-基因相互作用的概率，避免使用参数化机器学习方法。我们观察到，使用单细胞数据的表达相关顺式调节元件的这种确定固有地选择通常由超级增强子调节的基因，而不需要组蛋白修饰的独立分析，正如我们在使用多模态数据的先前工作中也证实的，28同时在单细胞中具有上下文特异性变异性。重要的是，我们还表明，FigR中的统计工具（峰值基因和TF基因）是可推广的，可以应用于真正的多模态数据集，分析来自同一细胞的染色质可及性和基因表达。我们认为，一种基