Science
人脑中的单细胞DNA甲基化和3D基因组结构
引文信息
结构化摘要控制大脑内各种细胞类型的复杂基因调控机制对于理解其在健康和疾病中的功能至关重要。最近高通量表观基因组分析的激增预示着这些基因调控框架的突破性发现。特别是，DNA甲基化与三维（3D）染色质形成相结合，是基因调控的基础。本研究全面分析了跨越多个区域的成人人脑细胞中的DNA甲基化和染色质构象。
原理
脑细胞中的基因表达受 DNA 甲基化和染色质构象的调节，这些过程之间存在深刻的相关性。通过深入探究单细胞中的这些表观基因组特征，我们可以更好地了解人脑的基因调控复杂性。这项研究旨在仔细绘制成人脑细胞中的DNA甲基化景观和染色质结构。
结果
该研究揭示了基于表观基因组的脑细胞类型分类，阐明了基于其表观基因组印记的不同分类。发现神经元和非神经元之间染色质接触距离的明显差异。从区室、结构域和环中收集到对脑细胞 3D 基因组组织的见解，以及它们与染色质可及性和基因表达的关系。此外，在3D基因组特征中出现了独特的细胞类型特异性。该研究还揭示了细胞特异性DNA甲基化的模式及其对基因调控网络的总体影响。发现了皮质和基底神经节的区域差异。一项比较研究强调了人类和小鼠之间脑细胞类型和DMR的保守性。最后，单细胞甲基化条形码（scMCodes）的出现展示了在精确识别人脑细胞类型方面的巨大前景。
结论
这项全面的研究展示了人脑的单细胞DNA甲基化和3D基因组结构图谱。它阐明了整个大脑细胞的细胞类型特异性和多样化的表观遗传结构。这个表观基因组图谱有望成为一种宝贵的资源，推动脑细胞多样性、基因调控机制和新遗传工具的起源的进一步发现。
描述复杂细胞类型背后的基因调控程序对于理解健康和疾病中的大脑功能至关重要。在这里，我们通过在来自三个成年男性大脑的 46 个区域的 51.7 万个细胞（39.9 万个神经元和 11.8 万个非神经元）中以单细胞分辨率探测 DNA 甲基化和染色质构象，全面检查了人脑细胞表观基因组。我们鉴定了 188 种细胞类型并表征了它们的分子特征。综合分析显示，DNA 甲基化、染色质可及性、染色质组织和基因表达在细胞类型、皮质区域和基底神经节结构中发生了一致的变化。我们进一步开发了单细胞甲基化条形码，该条形码使用选定基因组位点的甲基化状态可靠地预测脑细胞类型。这种多模态表观基因组脑细胞图谱为成人大脑中细胞类型特异性基因调控的复杂性提供了新的见解。
高通量表观基因组分析已被用于阐明大脑中巨大细胞复杂性背后的基因调控程序 （1–3）。5′-甲基胞嘧啶 （5mCs） 是哺乳动物基因组中最常见的修饰碱基。脊椎动物基因组中的大多数 5mC 发生在胞嘧啶-鸟嘌呤二核苷酸 （CpGs） 处。CG 差异甲基化区 （DMR） 通常被认为是顺式调节元件 （CRE） 的指示 （4， 5）。然而，在脊椎动物神经元系统中，5mC 在非 CG（或 CH、H=A、C 或 T）环境中也被大量检测到 （6）。CG 和 CH 甲基化（mCG 和 mCH）在脑发育过程中都是高度动态的，并显示出细胞类型特异性 （1， 4， 7）。它们对基因调控和大脑功能也至关重要 （8）。此外，基因调控需要染色质折叠的适当三维 （3D） 构象，该构象被组织成活性 （A） 或抑制性 （B） 区室、拓扑关联结构域 （TAD） 和染色质环 （9）。这些 3D 结构促进了基因启动子与其调控元件之间的相互作用，提供了额外但关键的调控机制层。DNA 甲基化和染色质构象在调控基因表达方面相互作用和协调，并且这些过程高度相关 （3）。对脑细胞这些表观基因组特征的调查可以加深我们对人类大脑复杂性背后的基因调控的理解。在这里，我们使用单核表观基因组测序技术全面分析了来自皮质和皮质下区域的成人人脑细胞中的 DNA 甲基化和染色质构象。
基于表观基因组的脑细胞类型分类法
我们剖析了 46 个大脑区域，包括大脑皮层（CX，22 个区域）、基底前脑（BF，两个区域）、基底核（BN，11 个区域）、海马体（HIP，5 个区域）、丘脑（THM，2 个区域）、中脑（MB，1 个区域）、脑桥（PN，1 个区域）和小脑（CB，2 个区域）的大脑结构（图 1A、图 S1A 和表 S1 和 S2）。除了两个杏仁核区域（BM和CEN，各两个重复）外，大多数区域有来自三个成年男性供体的3个生物学重复（表S3）（图S1A和表S1）。荧光激活细胞核分选 （FANS） 用于分离每个样品中 90% 的 NeuN 阳性细胞和 10% 的 NeuN 阴性细胞（图 S1A）。然后，我们使用 snmC-seq3 （“mC”） （10） 在单细胞水平上分析所有 46 个大脑区域的 DNA 甲基化 （DNAm）。此外，我们使用 snm3C-seq （“m3C”） （3） 同时检查来自 CX、BF 和 BN 的 17 个大脑区域的单细胞 DNA 甲基化和染色质构象（图 1B 和图 S1A）。经过严格的质量控制，确认了 378,940 mC 和 145,070 m3C 细胞核适合进一步分析（图 S1B）。每个 mC 细胞平均产生 94 万个过滤读长，每个 m3C 细胞产生 ~220 万个读长，有 40.6 万个染色质接触。这种数据质量使我们能够可靠地测量基因组特征的DNAm（图S1C），识别可变的甲基化区域，并精确定位不同脑细胞类型的TAD和染色质环。
图 1.使用 snmC-seq3 和 snm3C-seq 对人脑细胞进行表观基因组分析。
（A） 涵盖的人脑结构和区域。（B） snmC-seq3 和 snm3C-seq 的分析模式示意图。（C）人脑核的迭代聚类和注释。来自整个mC数据集，来自抑制性/非端脑神经元细胞类和来自SubCtx-Cplx主要类型的细胞使用t-SNE依次可视化，由相应迭代中注释的细胞组着色。（D） 主要类型的稳健树状图和亚型编号、大脑结构和供体来源的元信息。调色板在本研究中共享。（E） 主要类型SubCtx-Cplx的兴奋性和抑制性标志物（SLC17A1和GAD1）的CH甲基化。（F） 人脑细胞由解剖区域着色。（G） 通过 snm3C-seq 分析的脑核的 2D 可视化。（H） 脑细胞类型中整体 CG 和 CH 甲基化的变化。（I） 全球DNA甲基化与MECP2和DNMT1基因表达的相关性。
通过对mC数据集进行迭代聚类（参见材料和方法），首先将细胞核分为三类：端脑兴奋性神经元、抑制性和/或非端脑神经元以及非神经元细胞（图1，C和D）。这些进一步分为40种主要类型和188种亚型（图1D;图S2，A至C;以及表S4和S5）。根据神经元细胞的 CH 低甲基化基因标记物和非神经元细胞的 CG 低甲基化标记物对细胞类型进行注释（参见材料和方法）。所有主要类型和亚型在供体之间都是保守的，尽管某些细胞类型的比例存在细微差异（图1D和图S2C）。稳健的树状图表明主要类型和亚型之间存在相似性（图1D，图S2C以及材料和方法）。端脑兴奋性和抑制性和/或非端脑神经元与非神经元细胞很好地分离，每种类型都形成一个特定的分支，但 CB 和浦肯野 （PKJ） 细胞除外，它们与非神经元细胞类型分组，可能是因为它们具有相似的整体 CG 和 CH 甲基化分数（图 1H 和图 S4A）。
非神经元主要类型均匀分布在大脑结构中，而神经元主要类型表现出相当大的空间特异性（图1，D和F）。大多数端脑兴奋性神经元按位置分组（图1D）。海马兴奋性神经元根据其亚结构[CA1、CA3和齿状回（DG）]进行分组。皮质兴奋性神经元按其皮质层（例如L2/3）和投射类型[例如，端脑内（IT）投射神经元（表 S4）]聚类。基底核兴奋性神经元，主要来自杏仁核，形成 Amy-Exc 组。端脑抑制神经元表现为 11 种主要类型，主要来自皮质区域（Pvalb、Pvalb-ChC、Sst、Lamp5、Lamp5-Lhx6、Sncg 和 Vip）和基底核或基底前脑（MSN-D1、D2、Foxp2 和 Chd7）。在丘脑中，确定了一种兴奋性和两种抑制性主要类型。一种抑制性主要类型THM-MB与一小群中脑细胞具有相似的DNA甲基化特征。另一种抑制性主要类型THM-Inh非常罕见（361个细胞或整个数据集的0.07%），由于夹层污染，可能起源于丘脑的甲壳核（图S2D）。来自 Pons 的脑桥核神经元构成了一种特定的主要类型 （PN）。小脑包含两种不同的主要类型：罕见细胞型PKJ（867个细胞或0.17%）和小脑颗粒细胞（CB）。最后，在基底核和中脑中发现的SubCtx-Cplx主要类型以其异质性而著称：其亚型由兴奋性和抑制性细胞组成（图1E），并具有神经递质受体、转运蛋白和神经肽基因的高度可变DNAm（图S2E）。
从单核DNAm图谱确定的细胞类型与来自同一人脑的单核转录组测序（snRNA-seq）和单核染色质可及性测序（snATAC-seq）数据相印证[参见材料和方法以及Siletti等人的配套手稿。（11）和Li等人。（12）]. 综合分析揭示了使用不同分子模式确定的细胞类型之间的强烈对应关系（图 S3A）。所有基于表观基因组的细胞亚型都与基于转录组的簇非常吻合（图 S3B），尽管基于转录组的簇来自 ~10 倍的细胞和 ~2 倍的大脑区域。
全球甲基化因主要类型而异：mCG 为 77.7% 至 85.5%，mCH 为 0.8% 至 10.7%。非神经元细胞和颗粒细胞（DG 和 CB）主要类型在 mCG 和 mCH 中的总体分数最低（图 1H 和图 S4A），与之前在小鼠中的研究一致 （1）。皮质抑制性神经元的mCG最高，而丘脑、中脑和脑桥的某些非端脑神经元表现出最高的mCH（图1H和图S4A）。细胞类型全局甲基化与DNAm读取器和修饰因子的基因表达相对应（图1I和图S4B）。MECP2 和 DNMT3A（主要的 mCH 读写器）的表达与整体 mCH 呈正相关 [Pearson 相关系数 （PCC） = 0.39 和 0.35]，与 mCG 呈弱相关（PCC = 0.17 和 0.08）（图 1I 和图 S4B）。DNA 甲基转移酶 DNMT1 的表达与不同细胞类型的 mCG 之间呈高度正相关 （PCC = 0.63）（图 1I），与其在成熟神经元中作为主要 mCG 维持者的作用相匹配 （13）。我们观察到 DNMT1 表达与 mCH 之间的相关性更高 （PCC = 0.72;图 1I），但被认为对 mCH 影响不大 （14）。这意味着 DNMT1 和 mCH 之间存在未知关系，或者是影响 DNMT1 表达和 mCH 的一些尚未发现的因素。
使用改进的 scHiCluster （15） 用于 m3C 细胞，我们能够仅通过染色质接触分离除 MSN-D1 和 D2 以外的所有主要类型（图 1G）。这也突出了大脑区域染色质构象的多样性（图S2F）。为了保证两个数据集之间注释的一致性，我们迭代地将mC和m3C单元共聚类，然后将单元类型注释从mC转移到m3C单元（表S5和材料和方法）。
神经元与非神经元接触距离的差异
为了研究不同尺度的细胞类型特异性基因组折叠，我们首先检查了基因组距离处每个细胞的接触比例。神经元在较短的距离（200 kb 至 2 Mb）上显示出相互作用的富集，而成熟的少突胶质细胞和非神经细胞在更远的接触（20 Mb 至 50 Mb）上富集。星形胶质细胞和少突胶质细胞祖细胞在两个范围内都表现出富集（图2，A至C，图S5，A和B）。在神经元细胞中，皮层兴奋性和皮层下神经元比皮层抑制细胞具有更多的短程相互作用（P < 1 × 10–300，Wilcoxon秩和检验;图2，A和B）。我们在小鼠（1）和不同技术[Dip-C（16）]的先前数据集中观察到了类似的模式，表明这些模式是保守的。在整个基因组中观察到神经元中短程接触的富集，包括神经元和非神经元基因位点（图 S5D）。较短和较长相互作用之间的比率与细胞的整体基因表达活性高度相关 （PCC = 0.87;图 S5E），并与细胞核的大小一致 [L5-ET >其他皮质兴奋性神经元>皮质抑制性神经元>非神经元 （17）]。这些结果表明，传统上与细胞周期阶段相关的接触距离谱 （18） 也可能因非分裂细胞中的细胞类型而异。
图 2.主要类型3D基因组结构的多样性。
（A） 每个细胞中与基因组距离的接触频率，每个细胞（列）内归一化的 Z 分数。细胞按主要类型分组，然后按细胞的中位数 log2 短/长比率排序。y 轴以 log2 比例分箱。（B） 主要类型的 Log2 空头/长头比率，顺序与 （A） 相同。（C） 四种主要类型的插补接触图。（D）热图表示（C）中距离归一化接触图的相关矩阵，折线图表示相关矩阵的PC1。（E） （D） 中两个矩阵的放大视图以及跨细胞 mCG 的相应相关矩阵。（F 和 G）FOXP2 位点（细胞类型 L4-IT 的标志物）（F）或 LAMP5 位点（Lamp5 和 Lamp5-Lhx6 的标记物）的兴奋性 IT 神经元的插补接触矩阵（热图）、边界概率（蓝线）、绝缘评分（橙线）、差分边界（线图中的红点）和差分环（热图中的青色点）（G）。灰色阴影代表基因体（TSS 到 TES）。（H） 使用结构域（顶部）或环（底部）作为特征的细胞 （n = 5707） 的 t-SNE 图，按主要类型着色。（I） 区室评分、边界概率或环相互作用强度与 ATAC 信号之间的 PCC、所有主要类型的 bin 的 mCG 和 mCH 分数（左）或仅顶部 DEG（右）。（J） 使用所有基因（左）或顶部 DEG（右）的不同类别重叠（x 轴）的区室评分、边界概率或环相互作用强度之间的 PCC，以及所有主要类型中的基因表达。（K 和 L） 相对于基因位于不同位置的所有箱中显着正相关或负相关的区室 （K） 或结构域边界 （L） 的比例，在顶部神经元 DEG 中平均。（M） 显著相关的环像素占所有环像素的比例（左），正相关和负相关环像素之间的比率（中），或位于相对于基因的不同位置的显著相关环像素（右）的平均PCC，顶部神经元DEG的平均值。（N） 在顶部神经元 DEG 中，具有显着正相关区室、与基因体重叠的域边界或基因体内的环像素或具有至少一个锚点的基因数量与基因的 TSS 或 TES 重叠。35个基因未包含在三个圆圈中的任何一个中。
接下来，我们研究了富集的长程或短程染色质相互作用与染色质区室或结构域之间的关系。我们以 100 kb 分辨率鉴定了每种主要类型中的染色体区室（图 2D），并以 25 kb 分辨率鉴定了结构域。非神经元中富集的长程相互作用主要是室内相互作用，尤其是B区室区域之间。神经元中的短程相互作用也在同一区室内富集（图 S5、F 和 G）。总的来说，我们观察到非神经元中筋膜室内相互作用的富集和筋膜间相互作用的消耗（图 S5、H 到 K 以及材料和方法），表明室强度更强。相比之下，神经元中短程相互作用的富集被发现是域内和域间（图 S5、L 和 M）。
脑细胞类型中的区室、结构域和环
我们假设，如果两个基因组位点在物理上接近，它们的甲基化状态将共变异。共甲基化系数矩阵描绘了单细胞基因组仓之间甲基化的相关性，显示出与染色质接触的区室结构相呼应的格子图案（图2E和图S6A）。这表明基因组被分离到局部共甲基化结构域中，这些共甲基化结构域构成了两组具有相反甲基化多样性的集合。在 scATAC-seq 数据中也观察到染色质可及性具有类似的共调节结构 （19），这加强了基因组区室化的证据。在探索DNA甲基化与3D基因组结构之间的联系时，我们观察到染色质相互作用的强度与其锚的平均甲基化分数之间的相关性也与染色体区室有关（图S6B），其中负相关更频繁地发生在活性区室中（P < 1 × 10–300;图S6C）。
然后，我们确定了单细胞中分辨率为25 kb的结构域，发现神经元的结构域（中位数为4813）比非神经元（中位数4308）多（P < 1 × 10–300），但平均大小较小，导致结构域覆盖的基因组比例相似（图 S6、D 和 E）。结构域的数量和大小与全局基因表达活性高度相关（图S6F）。基因组箱的边界概率被定义为它被识别为跨细胞域边界的频率，这反映了细胞类型伪散装接触图的绝缘分数（图 2、F 和 G）。
以 10 kb 分辨率描绘了 29 种主要类型（和 119 种细胞亚型）中每一种具有 ≥100 m3C 细胞的染色质环。我们检测到主要类型（亚型）的中位数为 524,935 （541,551） 个环像素，其中 45,140 （59,905） 个环峰（图 S6G）。其中，24.3%是远端DMR（有关DMR的系统描述，请参阅后面的部分）和基因启动子[转录起始位点（TSS）±2 kbp]之间的相互作用，远端DMR之间的相互作用为38.1%，启动子之间的相互作用为5.8%（图S6H）。
3D基因组特征的细胞类型特异性
使用区室评分、结构域边界概率或环路强度，我们能够区分细胞类型并确定其相似性的层次结构（图 2H;图 S7，A 到 C;以及材料和方法），表明这些 3D 结构的细胞类型特异性。特别是，染色体区室可以区分非神经元、兴奋性、抑制性和MSN神经元，但对于兴奋性或抑制性细胞类别中更精细的主要类型有困难（图S7、B和C）。相比之下，染色质结构域和环路在更精细的兴奋性和抑制性主要类型中区分得更好，并且环路表现最好（图 2H 和图 S7、B 和 C）。这强调了不同尺度的 3D 特征在跨细胞类型粒度的基因调控中的不同作用，突出了环可能比结构域更具体。这些分析的主要目的是在细胞类型特异性方面对比区室、结构域和环，而不是细胞类型聚类。染色质接触上细胞聚类的最新技术仍然基于基因组bin对，正如我们（15）和其他小组（20,21）所采用的那样（图S7，B和C，以及材料和方法）。
对所有（或神经元）主要类型的特定 3D 结构进行系统检查，确定了 1188 （1024） 个差分区室 （DC）、2050 （1720） 个差分域边界 （DB） 和 173,806 （148,395） 个差分环 （DL）（图 S8A 和表 S6）。染色质结构域通常被认为在细胞类型中是保守的 （22–25），而它们可以在细胞类型和发育中显示某些动态 （3， 26–28）。我们的数据进一步表明，即使在密切相关的细胞类型之间，染色质结构域也可能有所不同（图2，F和G）。DMR-DMR环比启动子-DMR或启动子-启动子环显示出更高的细胞类型特异性（图S8B）。在评估差异染色质环中的转录因子 （TF） 时，我们发现细胞类型特异性 TF（例如 NFIX 和 NHLH1）的基序在 DL 的锚点处更加富集，而 CTCF 是染色体结构的关键 TF，在管家环处高度富集（图 S8C）。这意味着CTCF的作用更多地在结构环中的作用，而不是在细胞类型特异性启动子-增强子相互作用中的作用。许多神经元TF（例如，NEUROG1和NEUROG2）在泛神经元环上富集，但在泛脑细胞环上富集（图S8C），这与它们的神经元特异性作用一致。
基因组组织与其他分子模式之间的关系
我们研究了不同的 3D 结构特征与其他表观基因组模式（mCG、mCH 和开放染色质）之间的联系。在神经元细胞类型中，mCG 和 mCH 都与区室评分、结构域边界概率和环路强度呈抗相关（图 2I 和图 S9、B 和 D 以及 S10、C 和 D）。相比之下，开放染色质信号与这些结构特征呈正相关，相关性相似或稍弱（绝对）（图2I和图S9、B和D，以及S10、C和D）。这些（反）相关性表明活性区室、强结构域和环相互作用之间的协调，以及与活性染色质状态相对应的开放染色质和甲基化耗竭。在不同细胞类型的差异结构特征（DC、DB 和 DL）之间，与 DC 和 DB 相比，DL 与 mCG、mCH 和开放染色质具有更强的（反）相关性（图 S10C），尤其是在跨细胞类型具有高变异性的环上（图 S10D）。所有细胞类型的相关性通常比单独的神经元弱（图2I和图S9和图10）。在 DNAm 和 3D 基因组结构之间观察到的反相关性可能是由于 DNAm 对驱动基因组折叠的因子（例如 CTCF）结合的影响、通过甲基化和基因组组织的高级结构形成或共享调节因子 [例如，小鼠皮层中的 Neurog2 （30)].需要进一步的发展或机制研究来解决因果关系（29,31）。
基因表达也与3D基因组结构相关，特别是对于细胞类型特异性基因（图2J）。我们鉴定了 1099 （1358） 个顶级差异表达基因 （DEGs），这些基因在神经元（所有）主要类型中成对排列。它们与与其基因体或启动子重叠的所有三个结构特征都表现出强烈的正相关（图 2J 和图 S11 至 13）。对于环，与锚定包涵的DEGs相比，锚定重叠的DEGs（在基因体或启动子上）的相互作用强度更相关（P < 1 × 10–300; 图2J）。bins的基因表达和/或结构信号变异性的增加与它们之间更高的正相关性有关，这证实了差异结构特征和差异基因表达之间的重叠（图S11，E和F，以及S12，E和F）。
我们进一步研究了 1099 个神经元 DEG 与其相关染色质结构之间的相对位置 [错误发现率 （FDR） < 0.01;参见材料和方法] 在周围区域 [TSS-5Mb 到转录末端位点 5Mb （TES 5Mb）]。相关区室主要位于基因体内（图2K），相关结构域边界在TSS和TES高度富集（图2L），表明基因体区室和结构域的动态与基因表达多样性相关。具有正相关的环在基因体内以及 TSS 和 TES 与基因体之间富集± 1 Mb 区域（图 2M）。具体而言，基因体内48%的环与基因表达相关，其中98%与基因表达呈正相关。相比之下，基因体外的环与表达相关比例要小得多。在DEGs的上游和下游区域以及上游或下游区域与基因体区域之间观察到较高比例的正相关环，表明基因结构的调控域。
在 1099 个 DEG 中，453 个 （41.2%） 的基因体被一个或多个具有正相关区室评分的基因组箱重叠，591 个 （53.8%） 与一个或多个相关的结构域边界重叠。共有 1037 个 （94.4%） DEGs 具有 TSS 或 TES 锚定的相关环，898 个 （81.8%） 在基因体内具有相关环。这些不同尺度的染色质结构动态总共覆盖了96.8%的DEGs（图2N），再次表明基因组结构与基因表达多样性之间存在很强的关联。总的来说，这些分析揭示了染色质结构的细胞类型特异性及其与其他表观基因组和转录组特征的关系，在人脑中以前所未有的细胞类型分辨率。
细胞类型特异性 DNA 甲基化模式和相关基因调控景观
为了描述细胞类型特异性甲基化谱，我们在 188 个脑细胞亚型中鉴定了 24,455 个 CH 和 13,096 个 CG 差异甲基化基因 （DMG）（图 S14A 和材料和方法）和 2,059,466 个 CG-DMR（图 3A 和材料和方法）。除了描述脑细胞身份的不同表观遗传特征外，这些甲基化模式还为理解脑细胞中的基因调控程序提供了重要的见解，基因体甲基化与基因表达呈负相关 （5， 7， 32），DMR 标记假定的 CRE （4， 5），TF 基序涉及候选细胞类型特异性调节因子 （32）。
图 3.脑细胞中的基因调控。
(A） 188 种细胞亚型的细胞类型特异性 DMR 的 mCG。（B） CH-低甲基化TF及其基序在CG低DMR中的富集。下图显示了TF PBX3在整个基因组中其潜在结合位点的平均甲基化分数。（C） 确定推定CRE的工作流程。 （D）推定CRE的甲基化与不同滤波的相应基因之间的相关性分布。（E） 推定的 CRE 数量和与不同过滤的开放染色质区域的重叠比例。（F） 热图显示推定的 CRE 的 mCG、mCH 和靶基因的表达，以及相应环的接触强度。（G） 主要类型 L2/3-IT 和 MSN-D1 中基因 SYT1 周围的基因体 mCH、DMR mCG 和 3D 染色质组织。（H） 热图显示了从主要人类细胞类型中鉴定出的与DMR中指定性状或疾病相关的变异的LDSC分析结果。星号表示 P 值的大小 （*P < 1 × 10–1， **P < 1 × 10–2， ***P < 1 × 10–3和 ****P < 1 × 10–4).
如果TF是低甲基化的DMG（图S14B和材料和方法），并且它们的基序在相同细胞类型的低甲基化DMR（hypo-DMR）上富集，我们将TF分配到特定细胞类型（参见材料和方法）。总共有 612 个 TF 被分配到主要的神经元类型和亚型，它们在塑造和维持细胞身份方面可能发挥重要作用。例如，TBR1 被分配到深层兴奋性神经元，特别是 L6-CT 和 L6b（图 3B），并且注意到它在皮质突出神经元的发育中起决定性作用 （33）。ZNF423 和 EBF2 均被分配到小脑细胞类型（图 3B）。它们对小脑发育都至关重要，而 EBF2 特别指导 PKJ 细胞的迁移 （34–36）。
分析亚型进一步突出了 TF 利用率的变化。例如，TF PBX3 被分配到纹状体中普遍存在的 MSN-D1 主要类型，仅在纹状体区室的亚型中低甲基化，而不是纹状体的基质区室（图 3B 和图 S14C）。这表明 PBX3 在纹状体中表达的偏好，证实了先前的观察结果 （37， 38）。对 PBX3 潜在结合位点（具有 PBX3 基序的低 DMR）的进一步检查显示，纹状体亚型的平均甲基化分数较低（图 3B），表明该 TF 在纹状体中具有区室特异性调节作用。
我们集成了 DMG、DMR 和差分环，以确定每种细胞类型的推定 CRE（图 3C）。如果一个基因的 TSS 在 DMR 的 5 Mb 范围内，则该基因与 DMR 相关。进一步细化仅保留了与环路或 DL 的两个锚点重叠的 DMR-DMG 对。计算不同细胞亚型的 DMR 的 mCG 分数和基因体的 mCH 分数之间的 Pearson 相关性以评估这种关联（图 S14D）。观察到增强的关联，特别是对于DL过滤的DMR（Fig. 3D和图S14D），其与开放染色质区域的重叠也增加（图3E）。我们在 1,122,919 个 DMR 和 12,327 个基因之间鉴定了 320 万个潜在的调节性 DMR 和基因对（表 S7）。这些DMR、DMG的甲基化分数及其相互作用（环）的强度是（反）相关的（图3F），它们可以共同协调特定的基因调控程序。例如，与MSN-D1神经元相比，编码突触结合蛋白（一种关键的突触囊泡蛋白）的基因SYT1在L2/3-IT神经元中表现出较低的远端DMR和SYT1基因体的甲基化分数，并且DMR和启动子之间的相互作用更强（图3G），导致SYT1在L2/3-IT神经元中的表达高于MSN-D1神经元（图3G）和图 S14E）。总体而言，CG 和 CH 甲基化与染色质构象的整合揭示了不同的细胞类型调控动力学。
全基因组关联研究 （GWAS） 已确定许多与脑部疾病相关的非编码位点，其中许多位于增强子区域 （39）。DMR 和环有助于将这些遗传变异定位到特定的细胞类型调控元件。使用连锁不平衡评分回归 （LDSC） （40），我们检测到 20 种脑部疾病或性状与人脑细胞中的 DMR 或环重叠 DMR 之间的关联（图 3H、图 S15A，以及材料和方法）。精神分裂症、双相情感障碍和神经质风险变异在皮质和 HIP 中兴奋性神经元的低 DMR 中显着富集，而阿尔茨海默病与小胶质细胞一致 [MGC;图3H（41）]。与基底神经节的 Foxp2 细胞类型相关的烟草使用障碍变异（图 3H），该区域与烟草成瘾有关 （42）。对疾病风险变异的进一步探索揭示了对基因调控的不同影响。尽管许多细胞类型与相同的疾病有关，但所涉及的风险变异可能是多种多样的。例如，与L6-CT相比，精神分裂症风险变异rs2789588与具有相似表观遗传特征的L2/3-IT和L6-CT神经元有关，而rs17194490仅与L2/3-IT有关，具有特异性DNA低甲基化，与相应基因的长程相互作用更强，基因表达更高（图S15B）。
皮质和基底神经节的区域异质性
除了细胞类型多样性之外，新皮层的基因表达 （43–45） 和 DNA 甲基化 （1） 都注意到跨区域共享细胞类型的异质性。我们广泛的表观基因组数据集进一步探索了更广泛的皮质区域和皮质下区域的基因调控异质性。为了从其他细胞类型异质性中辨别区域多样性，我们设计了一个工作流程来揭示单核DNA甲基化谱中的区域景观（图4A）。将这些图谱与大脑区域数据相结合，我们将细胞映射到“区域甲基化空间”（图4A以及材料和方法），其中更靠近的细胞具有来自相似大脑区域的甲基化邻居。在这个区域甲基化空间（图4A以及材料和方法）中，轨迹描绘了区域转换以及相关的甲基化组位移，从而增强了我们对区域DNA甲基化效应的掌握。
图 4.皮质和皮质下细胞的区域轴。
（A） 从单核 DNAm 确定区域轴的工作流程。（B） 区域空间中皮层神经元的 2D 可视化，按解剖位置着色。（C）皮层神经元之间的共同区域轴。散点图显示了区域指数在每个皮质区域的变化情况。（D） 示例皮质夹层位置示意图。（E） rDMR 的 mCG 区域梯度，以及 L2/3-IT 细胞中 mCH 和 rDMG 的表达。（F） NR2F1周围染色质构象的区域差异。蓝色和紫色数字分别表示每个域的相对域强度和启动子强度。（G） （F） 中标记的示例差分环路重叠 rDMR 的放大视图。在递减结构域（左）中，甲基化分数从 V1C 增加到 A46 到 LEC，而甲基化分数在递增结构域中减小（右）。（H） 基底神经节中的抑制性神经元在 DNA 甲基化中表现出侧-背-腹轴。（I） MSN-D1 的 2D t-SNE 可视化。突出显示了来自 NAC、CaB 和 Pu 的细胞。（J） MSN-D1中基因体mCH甲基化和LSAMP表达的区域差异。（K） MSN-D1 中 LSAMP 周围染色质构象的区域差异。
皮层兴奋性神经元在甲基化中表现出明显的区域多样性，尤其是LX-IT神经元（图4B）。皮质抑制神经元 （46） 的区域多样性较少研究，因为它们在转录组和表观基因组中不明显 （1， 47， 48）。我们的分析揭示了皮层抑制神经元之间的区域性差异，尽管不太明显（图4B）。通过单细胞轨迹分析构建每种皮质神经元细胞类型的区域轴 （49）。我们观察到皮层神经元之间轴上大脑区域的共同顺序，从大脑的后部区域 [例如，初级视觉皮层 （V1C）] 到前外侧区域 [例如，前额叶皮层 （A46） 和中部颞回 （MTG）]，然后到前内侧区域 [例如，前扣带皮层 （ACC） 和外侧内嗅皮层 （LEC）;图4，C和D]。只有 L6-CT 显示出这种后外侧-前前-中前趋势的特殊模式（图 4、B 和 C）。然而，这种共同趋势允许进一步分析皮层神经元的共识区域轴（图4C以及材料和方法）。
沿该轴的表观遗传改变表明 CX 的区域规范。例如，TF NR2F1（也称为 COUP-TFI）在大脑发育过程中具有梯度表达，这对于建立新皮层 （43） 和新皮层和内嗅皮层之间的边界 （50） 的尾喙区域特化至关重要。我们的数据显示，V1C （P） 和 LEC （MA） 的基因体甲基化程度较低，A46 （LA;图 4G 和图 4G。S16B），并伴有基因表达的反向趋势（图S16A）。与 NR2F1 相关的两个染色质结构域显示出相反方向变化的相互作用强度（图 4F）。在 V1C 中，与 LEC 相比，上游结构域与 NR2F1 的启动子相互作用更多，并且具有低 DMR。相比之下，下游结构域与 NR2F1 的启动子表现出更强的相互作用，并且具有 LEC 中低甲基化的 DMR 特征（图 4、F 和 G）。两个相邻基因 NR2F1-AS1 和 FAM172A 分别包含在这两个结构域中，显示出与结构域强度一致的表达趋势（图 S16A）。表观遗传学和转录的这种连贯变异意味着调控结构域的切换和替代 CRE 的使用以激活不同皮质区域的相同基因，这需要进一步研究。
对皮层神经元中区域差异表观遗传特征的系统检查，总共确定了 14,606 个（每种主要类型平均 2.9 万个）区域 DMG （rDMG）、88.54 万个（63.2 千个）区域 DMR （rDMR）、77.3 万个（71.2 千）个区域差异环 （rDL） 和 1495 （136） 个区域差异域边界 （rDB）（图 S16B 和材料和方法）。许多 rDMG 和 rDMR 沿后-外侧-前-中前轴显示出单调甲基化梯度（图 4E 和图 S16、C 和 D），而更复杂的模式（如 NR2F1）也存在。
基底神经节神经元也表现出显著的区域多样性。侧-背-腹轴在基底神经节中变得明显（图4H），伴有表观遗传变化。例如，从NAC到CaB再到Pu，LSAMP基因在mCH中增加（图4，I和J），染色质结构域及其周围环的强度降低（图4K）。我们在 MSN-D1 和 MSN-D2 细胞中测定了四种主要类型的基底神经节（MSN-D1、MSN-D2、Foxp2 和 Chd7;图 S16B）中的 6371 个 rDMG 和 39.88 万个 rDMR，以及 98,276 （50,271） 个 rDL 和 193 （99） 个 rDB（图 S16B）。大多数 rDMG 和 rDMR 与外侧-背-腹轴表现出强（反）相关性（图 S16，C 至 F），突出区域差异是细胞内基底神经节异质性的关键。背侧（CaB和Pu）和基底神经节腹侧部分（特别是其主要成分纹状体）之间的功能和神经连接的区别，之前已经注意到（51,52）。 我们的数据和分析提供了背腹差异的表观遗传基础，并完善了背侧基底神经节内的区域差异（图4H）。
相当数量（746 个中的 427 个）TF 基序在 rDMR 中富集（图 S16G 和材料和方法）。这些TF中约有47%在相应的细胞类型中表达（图S16H），表达（抗）与区域轴相关（例如，图S16I）。这些发现暗示了大脑中潜在的区域特异性调节机制，可能是潜在的功能多样性。
人与小鼠之间脑细胞类型和DMR的保存
在几个新皮质区域观察到灵长类动物和啮齿动物之间的脑细胞类型保守 （17， 53）。为了评估这种保守性是否适用于更广泛的大脑区域，我们比较了人和小鼠的单核 DNA 甲基化谱 （1），使用相应的区域，包括 CX、BF、BN 和 HIP（参见材料和方法）。整合分析表明，人脑中定义的三种主要类型与小鼠脑细胞存在差异（图5A）。小鼠 L4-IT 神经元仅与人类对应物的亚群对齐（图 5B），证实了人类 L4-IT 神经元中更大的异质性 （17）。人HIP-Misc1神经元与部分小鼠皮质IT神经元整合，HIP-Misc2神经元与任何小鼠细胞类型均不匹配。平行 snRNA-seq 数据集 （11） 验证了这两种人 HIP 细胞类型（图 5C 和图 S17B）。尽管不匹配的细胞类型需要进一步研究，但主要类型分类法在人类和小鼠之间更广泛的大脑区域通常是保守的（图5A和图S17A），而人类的整体CG和CH甲基化始终高于相应细胞类型的小鼠（图5D和图S17C）。
图 5.人和小鼠脑细胞甲基化组之间的跨物种比较。
（A） 使用 2D t-SNE 可视化的人和小鼠大脑之间的单细胞甲基化组的整合。（B） L4-IT、HIP-Misc1 和 HIP-Misc2 型中人脑和小鼠脑细胞类型之间的差异。（C） HIP-Misc1 和 HIP-Misc2 细胞类型中 TF TSHZ2 的 CH 低甲基化和基因表达。（D） 人与小鼠之间保守细胞类型的全局 mCH 和 mCG 的相关性。（E） 细胞型DMRs的跨物种匹配示意图。 （F）总体而言，~50%的DMRs在其他物种中具有直系同源序列，其中~25%是相互DMRs。 （G）DMR甲基化（红色）和随机洗牌背景（黑色）的跨物种相关性分布。（H） hcCnsvDMRs甲基化部分的例子。（I） 组蛋白修饰标记中 hcCnsvDMR 的富集。（J） 主要类型 Pvalb 中基因 INPP5J 周围 hcCnsvDMR 的浏览器视图。红色区域是在 （54） 中验证的细胞类型特异性远端增强子。
为了比较人类和小鼠大脑之间的基因调控，我们使用 liftOver 来匹配在单个物种中鉴定的主要类型的低 DMR（图 5E）。跨细胞类型中大约 40% 至 60% 的低 DMR 在其他物种中具有直系同源序列（我们称之为 OrthSeqs）。大约一半的OrthSeqs在其他物种（OrthDMRs）中也有直系同源物和低DMR（OrthDMRs）。大多数 （95%） 的 OrthDMR 是相互匹配的（CnsvDMRs;图5F和图S17D）。CnsvDMRs的甲基化分数在人和小鼠之间显示出跨细胞类型的显著相关性（图5，G和H），表明物种之间的功能保守性。
我们进一步选择了相关性最高的DMR（hcCnsvDMRs，图5G）。hcCnsvDMRs的功能富集分析表明，它们在与前脑发育相关的生物过程以及与树突和突触相关的细胞成分中富集（图S17，F和G，以及材料和方法）。与小鼠前脑中的组蛋白修饰进行比较 （4） 表明，这些 DMR 从异染色区 （H3K9me3） 中耗尽，并在增强子（H3K27ac 和 H3K4me1）、启动子 （H3K4me3） 和定位增强子 （H3K27me3;图 S17E 和材料和方法） 区域富集。根据染色质可及性将 hcCnsvDMR 进一步分类为开放或关闭状态 （2） 表明，开放 DMR 富含增强子和启动子。相比之下，封闭的DMRs在蓄势增强子中特别富集（图5I），这些增强子可能在发育过程中一直处于活跃状态。
物种之间的甲基化保守性暗示了通过比较表观遗传学发现增强子的策略。例如，Pvalb 神经元的特异性基因 INPP5J 具有许多远端和近端 hcCnsvDMR 与匹配的染色质可及区域重叠（图 5J），其中两个被验证为小鼠 Pvalb 神经元病毒靶向的特异性增强子（图 5J）（54）。
单细胞甲基化条形码可靠地预测人脑细胞身份
细胞基因组中的 DNA 甲基化变异包含代表过去和现在基因调控事件的分子“印记”（55）。我们在许多对脑细胞类型高度特异性的 CpG 位点上观察到不同的 DNA 甲基化模式（例如，图 S18B）。这促使我们设计了单细胞甲基化条形码 （scMCodes），以使用所选 CpG 位点的甲基化状态在单细胞水平上确定脑细胞类型（图 6A、图 S18A 以及材料和方法）。
图 6.脑细胞类型的 snMCodes。
（A） 派生 snMCodes 的工作流程。（B） 来自所有三个供体的 snMCodes。（C） snMCode 特征的细胞类型特异性示例。（D） 热图显示了 snMCodes 在预测细胞类型中的混淆矩阵。（E） 交叉供体测试中的细胞类型预测准确性。（F） snMCodes 以单细胞分辨率预测具有有限数量 CpG 位点的人类细胞类型。
我们首先选择CpG位点，以迭代方式区分脑细胞类型（参见材料和方法）。这些位点根据其跨细胞型甲基化模式进一步聚类为 39,000 个组。然后，我们通过交叉验证的机器学习模型评估了他们的细胞类型预测能力（图 S18A 和材料和方法）。选择800组共12,000个CpG位点作为scMCodes（图6，B和C，以及表S8和S9），以获得良好的预测能力（图6D），同时最小化特征数（图S18C）。这些scMCodes达到了~93%的准确度（图6D以及材料和方法）。
在这项研究的三个供体和一个外部个体中进行了交叉供体测试 （5）。结果显示预测准确率高（~92%至96%;图 6E），展示了 scMCode 方法的跨个体鲁棒性。单细胞测序的基因组覆盖率有限。平均而言，每个细胞中仅检测到~200个scMCodes的CpG位点（图6F），这强调了scMCode在使用几百个选定的甲基化位点确定人脑细胞类型方面的有效性。
讨论
对人脑中细胞多样性和独特基因调控机制的深刻理解对于阐明脑功能和制定脑部疾病的治疗方法至关重要。我们编制了一份全面的人类大脑单细胞DNA甲基化和3D基因组结构图谱，其中包含来自46个不同大脑区域的524,010个深度测序的细胞核，使我们能够识别188种表观遗传学上不同的细胞类型。在这项研究中，对大脑区域的广泛分析使我们能够识别皮质下区域特异性的细胞类型，并比较不同大脑区域中同一细胞类型的表观遗传多样性，这大大扩展了以前的工作（3,5,56）。 此外，我们在了解跨脑细胞类型和区域的 3D 基因组多样性方面取得了长足的进步，这得益于使用 snm3C-seq 的细胞分析增加了 30 倍。此外，还确定了 29 种细胞类型的结构域和环的特异性，使细胞类型分辨率大大超出了以前的研究 （28， 57–59）。
单核苷酸分辨率 DNAm 已被证明在预测表观遗传年龄 （60）、追踪细胞谱系 （61， 62） 和诊断危及生命的疾病 （63， 64） 方面很有价值。DNAm 中编码的复杂调控信息使我们能够提取一组 scMCodes，以实现可靠的细胞类型鉴定。鉴于游离 DNA 甲基化已被公认为癌症诊断的可靠工具 （58），并为脑部疾病提供了有前途的生物标志物 （59），我们的 scMCode 方法是脑部疾病非侵入性诊断的潜在变革性工具。它可以帮助查明病理性脑细胞类型并为治疗选择提供信息，标志着精准医学向前迈进了一大步。
然而，有几条途径值得进一步探索。首先，除了本研究中采样的46个大脑区域之外，人脑还有其他具有复杂细胞多样性的复杂结构，特别是在皮质下区域。除了这项研究可用的范围之外，更全面的大脑区域采样将提供对潜在基因调控复杂性的更深入的见解。其次，高质量组织的可用性限制了我们只有三名男性捐献者。虽然这满足了本研究的目的，即调查人类脑细胞类型并揭示其表观基因组模式，但扩大供体基础将进一步阐明脑细胞的个体变异以及对基因调控多样性的遗传影响。最后，我们的发现很大程度上源于分子模态相关性。验证这些关联对于描述调控元件的功能、绘制调控网络图谱以及利用假定的增强子进行细胞亚型研究至关重要。
总体而言，这种多模态人脑细胞图谱以基础表观基因组学的观点丰富了我们对脑细胞的理解。它不仅为探索细胞类型多样性、基因调控复杂性、区域变异和脑细胞内的进化保守提供了宝贵的资源，而且还为开发用于细胞类型特异性靶向的创新遗传工具提供了基本要素，例如假定的调控元件。
材料与方法人类死后组织标本筛选
这些研究旨在使用最新的单核甲基化组（本研究）、snRNA-seq 和 snATAC-seq 技术首次探索人类大脑的细胞、转录和表观基因组变异。这些方法在最高质量的组织上表现最佳，使用优化的方法制备，包括短死后间隔 （PMI） 靶向<12 小时、高度一致的组织擦拭和照片记录以确保解剖学上精确的取样、用过冷异戊烷冷冻以保持组织完整性，并在 –80°C 冰箱下真空下正确储存。除了低PMI外，还对RNA完整性（>7.0）、传染病、头部外伤、插管、神经病理学和死亡方式应用了严格的排除标准。
鉴于高度严格的排除标准，组织的可用性是一个挑战。符合这些标准的脑标本，并且可以为目前的研究获得整个大脑半球，这是非常罕见的，并且具有严重的男性偏见。在 2018 年至 2022 年期间，有 16 名捐赠者符合这些标准，只有 3 名女性捐赠者最终根据 QC 或其他排除标准被排除在外。通过所有排除和QC标准的三名捐赠者均为男性。
人体死后组织标本处理
去识别的成人死后人脑组织是在获得死者近亲的许可后获得的。组织收集是根据 2006 年美国统一解剖礼物法案进行的，该法案在《加州健康与安全法典》第 7150 节（2008 年 1 月 1 日生效）以及其他适用的州和联邦法律法规中有所描述。此外，西方机构审查委员会审查了组织收集程序，并确定它们不构成需要机构审查委员会审查的人类受试者研究。
考虑将 18 至 68 岁且没有已知神经精神或神经系统疾病史的男性供体纳入研究。使用供体血液样本对传染病（HIV、乙型肝炎和丙型肝炎）进行常规血清学筛查，将传染病检测呈阳性的供体排除在研究之外。对样本进行RNA质量筛选，并考虑将平均RNA完整性值≥7.0的样本纳入研究。死后脑标本的处理如前所述（17）（dx.doi.org/10.17504/protocols.io.bf4ajqse）。简而言之，以 1 厘米的间隔切割冠状脑板，拍照，在干冰冷却的异戊烷中冷冻，然后转移到真空密封袋中，在 –80°C 下储存，直到进一步使用。对于感兴趣大脑区域的解剖，神经解剖学家对组织板的照片进行了注释，以勾勒出解剖的目标区域。然后，将组织板从 –80°C 冰箱中取出并短暂转移到 –20°C，在那里保持 ~1 至 3 小时，以使组织平衡至 –20°C。 然后将组织转移到保持在–20°C的定制温控冷台中，并使用标准剃须刀片或手术刀去除感兴趣的区域。将组织块储存在 –80°C 的真空密封袋中，直至以后使用。
细胞核分离和 FANS
使用如前所述的标准方案进行细胞核分离 （dx.doi.org/10.17504/protocols.io.y6rfzd6）。如前所述，对 DAPI 和 NeuN 荧光强度进行设控 （17）。将 NeuN 阳性和 NeuN 阴性细胞核分选到单独的试管中，并在分选后以 90% NeuN 阳性和 10% NeuN 阴性细胞核的规定比例合并。将分选样品离心，冷冻在 1× 磷酸盐缓冲盐水 （PBS）、1% 牛血清白蛋白 （BSA）、10% 二甲基亚砜 （DMSO） 和 0.5% RNAsin Plus RNase 抑制剂 （Promega， N2611） 的溶液中，并储存在 –80°C 直至进一步处理。将预分选的细胞核沉淀解冻并重悬于Dulbecco的PBS（DPBS） 1%BSA中，离心，重悬于1ml DPBS中，并分选到384孔板中。制备来自供体 H19.30.001 和 H19.30.002 的细胞核，并将其分选到 384 孔板中。对于供体 H19.30.004，从 AIBS 接收的冷冻组织块按照前面描述的程序进行处理 （7）。如前所述，标记细胞核进行 NeuN 荧光，并将其分选到 384 孔板中 （1）。
文库制备和Illumina测序snmC-seq文库制备
snmC-seq3 文库是使用 snmC-seq2 的更新版本制备的。简而言之，样品经过亚硫酸氢盐转化，并用随机引物进行条形码。然后通过两次SPRI净化将样品合并，将16×384孔板压缩成1×96孔板。然后如前所述对混合的样品进行调整和扩增。接下来，通过另外两次SPRI清理来汇集和清理文库。最后，通过Qubit测定文库浓度并归一化以进行测序。使用带有S4流通池和150 bp双端模式的Illumina NovaSeq 6000仪器对从人脑组织生成的snmC-seq3和snm3C-seq（见下文）文库进行测序。
snm3C-seq 文库制备
对于来自供体 H19.30.001 和 H19.30.002 的一些样品，使用了预分选细胞核。将预分驴的细胞核沉淀解冻并重悬于DPBS 1%BSA中，离心并重悬于1ml DPBS中。对于供体H19.30.001和H19.30.002的剩余样品以及供体H19.30.004的所有样品，使用研钵和研杵粉碎冷冻组织。然后立即将所有样品与2%甲醛溶液交联5分钟，用0.2M甘氨酸淬灭5分钟，离心，用DPBS洗涤，并在–80°C下储存，直到准备进一步处理。接下来，使用适用于 snm3C-seq 的 Arima 试剂盒在 37°C 下调节并消化细胞核 1 小时，在 65°C 下消化 20 分钟以灭活酶，然后在室温下连接 15 分钟。最后，将细胞核重悬于1ml DPBS 1%BSA中，通过40μM过滤器过滤，并像snmC-seq3样品一样进行分选。
snmC/snm3C-seq 的详细实验方案已在前面 （dx.doi.org/10.5281/zenodo.8319891） 中描述。
供体特异性基因组gDNA文库制备方案
使用DNeasy Blood and Tissue Kit（Qiagen，Valencia，CA）从研磨的冷冻组织中提取基因组DNA。用 Covaris S2 （Covaris， Woburn， MA） 将 1 微克 DNA 片段化至 300 bp，然后进行末端修复 （Lucigen） 并添加 3′ A 碱基 （New England Biolabs）。将 Illumina（Illumina，SanDiego，CA）提供的胞嘧啶甲基化接头在 16°C 下用 T4 DNA 连接酶 （New England Biolabs） 与超声处理的 DNA 连接 16 小时。用 AMPure X P 微球（Beckman Coulter Genomics， Danvers， MA）进行两轮纯化分离接头连接的 DNA。通过四个循环的聚合酶链反应（PCR）富集衔接的DNA分子，反应组成如下：25μlKapa HiFi Hotstart（KapaBiosystems，Woburn，MA）和5μlTruSeq PCR引物混合物（Illumina）（最终50μl）。热循环参数为：95°C2 min，98°C30 s，98°C15 s，60°C30 s，72°C1 min4个循环，1个72°C5 min。使用AMPure X P微球纯化反应产物。接头连接纯化的PCR反应产生了用于NovaSeq 6000后续测序的文库。
来自供体基因组测序的变异检出
首先使用软件fastp（v0.20.1）对全基因组测序读数进行QC处理（65）。使用的命令行是“fastp -i input_PE_R1.fastq.gz -I input_PE_R2.fastq.gz -o output_PE_R1.fastq.gz -O output_PE_R2.fastq.gz -w 4”。然后使用软件 BWA （v0.7.17） （66） 将 QCed reads 映射到人类基因组组装 GRCh38 （hg38），并使用 BWA-MEM 算法，并通过软件 samtools （v1.10） （67） 以 bam 格式存储映射结果。具体来说，用于映射的命令行是“bwa mem -t 20 hg38-ref input_PE_R1.fastq.gz input_PE_R2.fastq.gz |samtools 视图 -Sb - > output.bam”。使用基因组分析工具包（GATK，v4.1.8.1）（68）的种系短变异发现工作流程分析定位的读数。简而言之，我们首先从映射的读段中删除重复的读段，然后通过基础质量评分重新校准步骤 （BQSR） 来生成可用于分析的读段。BQSR 步骤中使用的变异参考是 dbSNP138、Mills 和 1000 Genomes 金标准插入缺失，以及 1000 Genomes 1 期单核苷酸多态性 （SNP）。接下来，使用GATK的HaplotypeCaller从分析就绪读段中调用候选变异（SNPs InDels），并用变异质量评分重新校准步骤（VQSR）进一步过滤，分别确定高置信度SNP和InDels。用于重新校准 SNP 质量评分的变体参考是 Hapmap 3.3、OMNI 2.5、1000 Genomes phase 1 和 dbSNP138，InDels 是 Mills 和 1000 Genomes 金标准插入缺失和 dbSNP138。BQSR 和 VQSR 中使用的所有引用都是从 GATK 资源包 （https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0） 下载的。
供体特异性参考基因组
对于每个供体，我们使用 GATK 的 SelectVariants 函数选择高置信度纯合子 SNP，并使用 GATK 的函数 FastaAlternateReferenceMaker 将纯合子 SNP 代入 hg38 FASTA 文件，从而创建供体特异性参考基因组
常见的纯合子SNP
通过比较供体的纯合SNP，我们构建了三个供体之间共享的常见SNP列表。
映射和计数/特征矩阵生成
对于 snmC-seq3 和 snm3C-seq 数据集的序列读取映射，我们使用了自己的自定义管道（https://github.com/lhqing/cemba_data，版本 1.2.1.dev94 gc65e173）。该管道的主要步骤包括：（i） 将 FASTQ 文件解复用为单单元;（ii）读取级质量控制（QC）;（iii） 测绘;（iv） BAM 文件处理和 QC，以及 （v） 最终分子图谱生成。前面已经描述了这五个步骤的细节（10）。我们将所有reads映射到供体特异性基因组。定位后，我们计算了每个细胞中两组基因组特征的甲基胞嘧啶计数和总胞嘧啶计数。hg38基因组的非重叠染色体100-kbbin（由“bedtools makewindows -w 100000”生成）用于聚类分析，人类GENCODE v33定义的基因用于聚类注释和与数据集的集成。如前所述，特征的 CG 和 CH 甲基化水平均归一化 （1）。从每个特征集的归一化甲基化水平生成按特征划分的细胞矩阵。
质量控制措施
根据以下指标对测序的细胞进行过滤：（i） mCCC% < 0.06;（ii）全球mCG%>0.5;（iii）全球mCH%<0.15;（iv） 最终读数总数> 250,000 次;（v）映射率>0.5。对于用 snm3C-seq 分析的细胞，我们要求细胞具有> 50,000 个顺式接触，距离为 >2500 bp。
snmC-seq3数据的聚类和注释聚类分析
100 kb 基因组箱的 CG 和 CH 甲基化水平被用作聚类的输入特征。我们使用软件包 ALLCools （https://github.com/lhqing/ALLCools） 迭代进行聚类分析。在每次迭代中，首先通过删除平均总胞嘧啶碱基检出量为 <250 或 >3000 的箱来过滤 100 kb 的箱。那些与ENCODE黑名单（69）重叠的也被排除在聚类分析之外。然后使用支持向量回归 （SVR） 从 CG 和 CH 甲基化中分别选择前 5000 个高可变特征 （HVF）。然后，我们对每 5000 个特征应用主成分分析 （PCA） 以降低维度。通过双样本Kolmogorov-Smirnov检验，以调整后的P<0.1，为标准，为每种甲基化类型选择前n个主要成分（PC），直到第n个和（n 1）个PC的分布没有显著差异。我们进一步对每个顶级PC组进行了预聚类，并选择了在预聚类中富集的PC（调整后的P值≤0.05）。最后，将 CG 和 CH 甲基化 PC 中选定的 PC 连接起来进行进一步分析。我们在选定的 PC 上使用 Harmony （70） 来消除个体差异。协调特征被进一步输入到前面描述的共识聚类程序中 （1）。
双峰/碎片识别
一块板中每个细胞的读取数在 snmC-seq3 和 snm3C-seq 中均稳定。因此，我们采用了基于细胞相对读数的双峰/碎片检测策略。我们首先将每个细胞的读取数归一化为其板的平均读取。具有板相对读取数>1.2或<0.8的细胞被认为是双峰/碎片候选者。在每次聚类迭代后，如果聚类含有 >80% 的双峰/碎片候选者，则聚类将被标记为双峰/碎片，并从进一步分析中剔除。
单元格类型注释
这些簇根据其低甲基化基因手动注释为主要或亚型，这些基因要么是典型的脑细胞类型标记，要么是从当前数据集中从头确定的。与其他细胞类型相比，我们要求每种细胞类型在 CG 和 CH 甲基化中至少有 5 个 DMG。否则，它将与最近的集群合并。如果候选细胞类型的所有细胞都来自单个供体，则该候选细胞类型将被标记为异常值。
在主要类型水平上，在可能的情况下，使用文献中先前描述的已知脑细胞类型的命名法对细胞类型进行注释[例如，（13）];否则，根据细胞类型的区域组成或不同的标记基因对细胞簇进行注释。前一种方法的一个警告是，使用在啮齿动物或非人灵长类动物中定义的标记基因注释的细胞类型可能无法反映人类细胞类型中的相应基因活性。例如，SNCG基因在与小鼠Sncc细胞相对应的人类主要类型中低表达。然而，使用通用命名法有助于跨物种比较和现有知识转移。
snm3C-seq数据的聚类和注释
为了注释来自 snm3C-seq 的细胞，我们将它们与注释的 snmC-seq3 细胞相结合，并进行了类似于上述的迭代聚类分析。唯一的区别是，使用软件Scanorama（71）校正了个体和测序技术的批次效应。每次聚类迭代后，使用 K 最近邻 （KNN） 分类器将细胞类型注释从 snmC-seq3 细胞转移到 snm3C-seq 细胞。
细胞类型的稳健树状图
我们在不替换的情况下从每种细胞类型中重新采样一定数量的细胞，以计算具有 100 kb CG 和 CH 甲基化基因组特征的细胞类型的平均甲基化组谱。主要类型的重采样数为 800 次，亚型的重采样数为 500 次。然后使用平均剖面来计算成对相关距离。该过程重复 500 次以计算平均成对距离矩阵，然后通过具有平均连锁的分层聚类来构建最终的细胞型树状图。
确定 DMG
我们分别测定了用于 CG 和 CH 甲基化的细胞亚型之间的 DMG。为了避免因细胞类型细胞数量不平衡而引起的潜在偏差，我们将每种细胞亚型中的细胞样本减少到不超过 500 个。使用 Wilcoxon 秩和检验测试由人类 GENCODE v33 定义的所有蛋白质编码和长链非编码 RNA 基因 （lncRNA） 是否显着甲基化减少（或低甲基化）。使用 Benjamini-Hochberg （BH） 程序通过多测试校正调整 P 值。我们计算了候选基因的受试者工作特征下面积 （AUROC） 曲线。调整后的 P 值≤ 0.001 和 AUROC ≥ 0.8 的基因被认为是 CG 和 CH 甲基化中的成对 DMG。
确定 DMR
我们将单细胞DNA甲基化谱合并到细胞类型（主要类型/亚型）谱中，根据它们的簇注释，以供体聚集和供体分离的方式。在进一步分析之前，过滤掉这些甲基化谱的非常见纯合SNP CpG位点。然后，我们使用了软件 MethylPy [v1.4.2;（72）] 以确定供体聚集图谱的所有细胞类型的 mCG DMR。使用的命令行是“methylpy DMRfind–output-prefix OUTPUT_FILE_NAME–samples SAMPLE_NAMES–mc-type CGN–dmr-max-dist 250–sig-cutoff 0.01–allc-files MC_FILES”。如果连续的 DMR 距离在 250bp 以内，并且它们在 188 个亚型中的 mCG 分数的 Pearson 相关性大于 0.8，我们进一步合并了它们。我们通过评估供体聚集和分离图谱之间跨细胞类型甲基化模式的可重复性来进一步筛选每个 DMR。评价标准为不同细胞类型间mCG组分之间的Pearson相关系数为≥0.5，平均绝对误差为≤0.1。
然后根据其 mCG 分数与其稳健平均值的差异，将每个可重复的 DMR 分配为每种细胞类型中的低 DMR 或高 DMR。通过平均不同细胞类型第 25 个和第 75 个百分位数之间的 mCG 分数来计算每个 DMR 的稳健平均 m。mCG分数大于m 0.3的DMR被分配为每种细胞类型中的hyper-DMR，低于m-0.3的DMR 被分配为低DMR。仅包含一个 CG 位点或没有任何低 DMR 或高 DMR 分配的 DMR 被排除在进一步分析之外。
基序富集分析
共使用来自JASPAR2020 （73） 的 746 个 TF 结合图谱（基序）进行基序富集分析。首先将细胞类型特异性低 DMR 分割成 500 bp 的 bin，然后通过与基序的基因组位置相交来注释每个基序。基序基因组位置是从 http://expdata.cmmt.ubc.ca/JASPAR/downloads/UCSC_tracks/2020/hg38 下载的。为了测试主要细胞类型的基序富集，使用低DMR作为前景信号，并使用所有其他主要类型的低DMR作为背景。对于细胞亚型水平的测试，仅使用属于同一主要类型的其他亚型中的低DMR作为背景。采用单侧Fisher精确检验计算前景在背景下的富集P值。
不同单细胞数据集之间的集成用于人类单细胞 DNA 甲基化、表达和开放染色质的特征基质
从mC数据集中确定的主要和亚型水平的CG和CH细胞型标记基因作为整合分析的特征。与单细胞RNA测序（scRNA-seq）整合时[参见配套手稿Siletti等人。（11）] 或 snATAC-seq 数据集 [参见 Li 等人的配套手稿。（12）]，我们使用了相反的基因体甲基化值，因为它们通常与基因表达密切相关。scRNA-seq和snATAC-seq数据集均通过特征基因的平均总唯一分子标识符（UMI）计数进行归一化，然后通过log（x 1）进行转换。神经元细胞类型和非神经元细胞类型分别整合。神经元细胞类型采用CH甲基化，非神经元细胞类型采用CG甲基化。在整合非神经元细胞类型之前，应用了额外的过滤步骤，这要求每个细胞的特征基因的总UMI在scRNA-seq和snATAC-seq数据集中均为>3000。
用于人和小鼠单细胞 DNA 甲基化的特征基质
我们仅使用人和小鼠之间的同源基因进行整合分析。同源基因列表从小鼠基因组信息学 （MGI） 数据库 （https://www.informatics.jax.org/homology.shtml） 下载。同源基因是从与scRNA-seq和snATAC-seq数据集整合时使用的相同特征中选择的。来自 THM、MB、CB、PN 和内嗅皮层的人脑细胞被排除在外，因为公共小鼠数据集中不存在对应细胞 （1）。小鼠数据集的重新注释方式与人类数据集相同。CG甲基化分别用于整合神经元和非神经元细胞类型。
一种整合不同单细胞测序数据集的方法
在特征矩阵生成后，我们使用类似于 Seurat v3 的三步法将两个数据集 X 和 Y 投影到同一空间：（i） 使用规范相关分析来捕获数据集之间单元格之间的共享方差;（ii）在两个数据集之间找到五个相互最近邻的锚点;（iii）将两个数据集拉入同一空间。为了实现可扩展性，我们从每个数据集中随机选择了 20,000 个细胞 （X裁判和Y裁判） 作为拟合典型相关分析的参考，并转换其他像元 （Xqry（qry）和Yqry（qry）） 到同一个 CC 空间。具体来说，规范相关向量 （CCV）X裁判 和Y裁判（表示为U裁判和V裁判） 通过对其点积的奇异值分解计算，����������=���������哪里���������=我和���������=我.然后，CCVXqry（qry）和Yqry（qry）（表示为Uqry（qry）和Vqry（qry）） 的计算公式为����=����(���������)/�和����=����(����������).U 和 V 通过除以每行的 L2 范数进行归一化，并使用与 Seurat v3 相同的方法查找相互最近邻锚点并对锚点进行评分。X 和 Y 也被垂直组合，并且该组合矩阵的 PC 使用与 Seurat v3 相同的方法通过上一步生成的锚点集成在一起。此集成步骤将一个数据集（查询）的 PC 投影到另一个数据集（参考）的 PC，同时保持参考数据集的 PC 不变。生成的 PC 用于可视化和查找数据集之间的匹配聚类。
3D基因组分析
在单细胞和伪本体水平上分析了3D基因组特征。对于单细胞分析，我们使用 scHiCluster （15） 以 100 kb 的分辨率（焊盘 = 1）插补接触矩阵，25 kb 分辨率（焊盘 = 1）用于 10.05 Mb 以内的接触，以及 10 kb 分辨率（焊盘 = 2）用于 5.05 Mb 以内的接触。为了加快 10 kb 分辨率的插补，在每条染色体上每个 30 Mb 滑动窗口内以 10 Mb 的步长进行卷积和随机游走。仅使用滑动窗口中心 10 Mb 内的值作为最终结果。
对于伪体分析，我们通过取原始基质的总和或组内细胞上插补基质的平均值来合并每个组（主要类型、亚型或区域）的细胞，并且仅随机选择 1500 个细胞，如果该组包含 >1500 个细胞，则使用。我们还有一个小组合并了区室分析和嵌入比较中使用的所有细胞类型。该组总共包含 5707 个细胞，由随机选择的 200 个细胞生成，每种主要类型有 ≥100,000 个接触点，但 L5-ET 除外，仅鉴定出 107 个细胞，全部用于分析。方法的详细信息如下所述。
接触距离分布分析
我们根据接触的两个锚点之间的距离为每个单个细胞生成了接触的直方图。在log2距离尺度上平均分配，步长为0.125，范围从2500 bp到249 Mb（最长染色体的长度）。第 i 个 bin 是距离介于 2500 × 2 之间的触点数0.125我和 2500 × 20.125(我 1).在图 2B 和图 2 中。S5，短长比定义为第 51 个 （200 k） 至第 76 个 （2 M） 箱中的触点比例除以第 103 个 （20 M） 至第 114 个 （50 M） 箱中的触点比例。
隔间分析
使用每条染色体的 100 kb 分辨率的伪本体接触矩阵进行区室分析。我们首先使用了 5707 个小区的合并接触图，并过滤掉了覆盖率异常的 100 kb 箱。具体来说，bin i 在染色体 c 上的覆盖率（表示为RC，I） 定义为染色体 C 接触矩阵第 i 行的总和。我们只保留了覆盖率在 99% 之间的垃圾箱Rc中位数的两倍Rc减去第 99 个百分位数Rc.通过距离对接触矩阵进行归一化，并计算归一化矩阵的Pearson相关矩阵。这些合并的接触图用于拟合每条染色体的 PCA 模型。以第一台PC（PC1）作为区室评分，调整模型符号，确保CpG密度较高的区室得分为正。我们目视检查了合并矩阵的 PC1，以确保这些值对应于相关矩阵的格子图案，而不是染色体臂。对每种主要类型的接触图进行滤波，并以与上述相同的方式转换为相关矩阵，然后使用PCA模型进行转换。原始矩阵和插补矩阵均用于此分析。我们使用合并的原始矩阵来拟合 PCA 模型，并将每种细胞类型中原始矩阵的相关矩阵转换为原始区室评分，并将每种细胞类型中插补矩阵的相关矩阵转换为插补区室评分。通常，插补基质在较小的细胞群中效果更好，而当合并足够的细胞时，原始基质可提供更高的分辨率。
为了确定富集的长程相互作用是室间相互作用还是室内相互作用，我们通过接触锚点处区室分数的差异或总和对接触距离图进行分层。分数的差异反映了接触是室间还是室内，较大的差异代表室间。分数的总和区分了室内接触的接触是AA还是BB，大的正总和代表AA相互作用，而小的负总和代表BB相互作用。请注意，一般来说，长程细胞比短程细胞包含更多的室间相互作用，因此在计算原始接触计数时，具有更多长程相互作用的非神经元细胞也比神经元具有更多的室间接触。我们在图 S5 中报告的结果，从 F 到 K，使用了距离归一化接触计数，它反映了每个距离处室内或室间接触的相对比例。因此，我们的结果仅表明，非神经元细胞中较高比例的长距离接触是室内接触，这并不表明非神经元细胞中室内接触的总数更多。
鞍座图和隔室强度的计算方法与（74）中描述的相同。具体来说，在每条染色体中，我们根据区室分数对所有 100 kb 的 bin 进行排序，并将这些 bb 分组为 50 个等间隔组。对每组各分柱之间的距离归一化相互作用强度或每对分柱之间的 mCG 或 mCH 水平的 PCC 进行平均。轴按细胞类型的区室评分进行排序，因此 BB 相互作用位于左上角，AA 相互作用位于右下角。使用原始区室评分作为输入，在所有主要类型或神经元主要类型之间使用 dcHiC （75） 鉴定 DC。很大一部分 （>60%） 的 bin 被识别为差异，传统的 q 值阈值为 0.01。因此，我们只选择了Z分数变换马氏距离>1.960（标准正态分布的97.5个百分位数）的前DC。
域和差分域边界的识别
使用 scHiCluster 以 25 kb 分辨率得出结构域和绝缘分数。具体来说，在以 25 kb 分辨率的插补矩阵上使用 TopDom （76） 识别每个单个细胞内的结构域。使用伪体积插补矩阵（单个细胞的平均值）和 10 个箱的窗口大小计算每个箱的每个细胞组（大脑区域内的主要类型或主要类型）的绝缘分数。条柱的边界概率定义为将条柱称为域边界的条柱在组中总量中所占的比例。
在单个单元中识别的结构域数量与可能影响插补性能的短程读取数量相关。为了避免计算伪影，我们从每种细胞类型中选择细胞以匹配跨细胞类型短程接触的分布，并在图 S6、D 和 E 中观察到相同的趋势，这表明结构域数和大小之间的这些差异不能完全由神经元和非神经元之间的不同短/长比率来解释。
为了识别 n 个单元组之间的差分域边界，我们为每个 25 kb 的 bin 推导出了一个 n × 2 列联表，其中每行中的值表示组中将 bin 称为边界或不作为边界的单元数。我们计算了每个 bin 的卡方统计量和 P 值，并使用整个基因组的统计量峰值作为差分边界。峰值定义为 FDR 内卡方统计量的局部最大值 <1 × 10–3（Benjamini 和 Hochberg 程序）。如果两个峰彼此相距在五个区格以内，我们只保留具有较高卡方统计量的峰。我们还要求峰具有 Z 分数变换卡方统计量 >1.960（标准正态分布的 97.5 个百分位数），以及最大和最小边界概率之间的差异 >0.05。
染色质环和差分环的鉴定
使用 scHiCluster （15） 分别鉴定每个大脑区域内每个主要类型、亚型和主要类型的染色质环。我们只执行 50 kb 到 5 Mb 之间的循环调用，因为增加距离只会导致有效环路数量的有限增加。对于每个单个细胞，每条染色体的插补基质 Q细胞对数变换，并在每个对角线上对 Z 分数进行归一化（结果表示为 E细胞），并减去 ≥30 kb 和 ≤50 kb 之间的局部背景（结果表示为 T细胞），类似于 SnapHiC （77）。计算伪批量水平 t 统计量以量化细胞组中单个细胞中 E 和 T 与 0 的偏差，其中较大的偏差表示在全局 （E） 或局部 （T） 背景下的富集度较高。E细胞也在每个对角线上洗牌以生成 E沙狐球，然后是 T沙狐球，以估计 t 统计量的背景。可以通过比较观察到的细胞与洗牌细胞的 t 统计量来得出经验性 FDR。我们要求像素的平均 E > 0，环形和左下角背景的倍数变化>1.33，水平和垂直背景的倍数变化>1.2 （77），以及 FDR <0.01 与全局 （E） 和局部 （T） 背景相比。使用广度优先搜索算法从环路像素中选择环路峰峰，我们从具有最大 E 的环路像素开始，并将其与 20 kb（L0 距离）内具有较小 E 值的所有其他环路像素连接起来。在环路像素的每个连接分量中具有最大 E 值的环路像素被定义为环路顶峰。我们只在计算环路峰值时使用了峰的概念，在所有其他情况下，我们使用“环”来表示环像素。
为了比较不同细胞组之间环路的相互作用强度，我们采用方差分析（ANOVA）框架，使用Q计算至少一个细胞组中识别的每个环路的F统计量细胞（结果表示为 FQ） 或 T细胞（结果表示为 FT).然后，我们 Z 评分为 FQ和 FT在所有正在测试的循环中，并选择带有 F 的循环Q和 FT> 1.036（标准正态分布的第 85 个百分位数）作为微分环路。阈值是通过目视检查接触图以及相互作用与环锚 CG 甲基化的相关性来确定的。
在控制相互作用强度和局部背景下的富集后，对差分环和恒定环之间的基序富集进行了分析。我们首先生成了一个常量循环池，其 Z 得分为 FQ和 FT< 0.然后，我们将微分和常数环路分为 100 × 100 组，这些组基于 FQ和 FT.我们从每组中为微分环和常数环选择了相同数量的环，并将微分环或常数环中的基序富集与它们的并集进行了比较。
在 11 项比较中确定了差分环路，包括所有主要类型之间的比较;神经元主要类型之间;神经元、神经胶质细胞（ASC、ODC 和 OPC）、MGC、PC、EC 和 VLMC 之间;在胶质细胞主要类型之间;兴奋性神经元（L2/3-IT、L4-IT、L5-IT、L6-IT、L6-IT-Car3、L5/6-NP、L6-CT、L6b、L5-ET 和 Amy-Exc）、抑制性神经元（Lamp5、Lamp5-Lhx6、Sncg、Vip、Pvalb、Pvalb-ChC、Sst 和 Chd7）、脑核神经元（MSN-D1、MSN-D2 和 Foxp2）和 SubCtx-Cplx 之间;在兴奋性主要类型之间;抑制性主要类型之间;大脑核主要类型之间;IT专业类型（L2/3-IT、L4-IT、L5-IT和L6-IT）之间;在尾神经节隆起 （CGE） 衍生的主要类型（Lamp5、Lamp5-Lhx6、Sncg、Vip）之间;以及中棘神经元主要类型（MSN-D1 和 MSN-D2）之间。结果见表S6。图S8显示了一些比较的聚合峰分析。对于每个单个环路像素，选择从 –100 kb 到 100 kb 的插补接触图，并将最小值归一化为 0 到 1 的范围，并在文件夹变化为 Q 和 T 分别大于 1.2 和 1.5 的所有差分环路上取平均值，比较前景单元类型的平均值和背景单元类型的平均值。
基于不同3D基因组特征的单细胞包埋基于联系人
将分辨率为 100 kb、距离为 ≥100kb 和 ≤1 Mb 的插补触点用作奇异值分解降维的特征。前 30 个主成分通过每个单元的奇异值和 L2 范数进行归一化，然后用于图 1G 和图 S2F 中的 t 分布随机邻域嵌入 （t-SNE） 可视化。为了更好地可视化神经元细胞的异质性，我们将每个神经元细胞群下采样到1000个细胞，并与所有神经元一起拟合t-SNE，并投射回其他非神经元细胞。Higashi（20）和fastHigashi（21）用于图S7中嵌入的比较。由于内存限制（256 G），我们只以 500 kb 的分辨率而不是 100 kb 的分辨率运行这些工具（与他们的论文和 Github 页面中建议的相同）。供体信息被用作模型中的混杂因素，以避免由供体差异驱动的嵌入。
基于隔间
使用 CpG 密度法或特征向量法在原始接触图或由 scHiCluster 或 Higashi 插补的接触图上计算单细胞区室分数。CpG 方法使用基因组中每个 100 kb bins 的 CpG 密度，并为每个单个细胞计算每个 bin 的值，作为同一染色体上其他 100 kb bins 的 CpG 密度的平均值，由两个 bin 之间的相互作用强度加权。特征向量方法使用所有单细胞的插补矩阵的平均值来计算相关矩阵（如区室分析中所述）并拟合 PCA 模型以转换所有单细胞的相关矩阵。Higashi 还用于在生成单元嵌入后对接触图进行插补，如上面的“基于接触”部分所述，在嵌入空间上使用 0 个邻居或 5 个邻居来帮助插补。尽管使用 5 个邻居生成了具有更高细胞类型特异性的区室，但该方法强制平滑了细胞嵌入的信息，当嵌入可以很好地分离细胞类型时，这可能会人为地增加插补基质之间的差异。因此，这使得很难断言细胞类型的分离是由于单个细胞中区室的固有异质性还是由于嵌入的平滑，我们仍然在图 S7 中显示了 0 个邻居的结果。然后，我们对逐个单元室评分矩阵 X 进行奇异值分解 = 无人艇T推导嵌入 U 的单元。
基于循环
我们组合了在所有主要类型中识别的环路像素，以制作一个元循环列表，并生成了一个二进制逐个循环矩阵，其中每个元素都指示是否在环路像素处的单元格中检测到接触。对二元矩阵（表示为A）应用对数项频率的潜在语义分析来计算嵌入。具体来说，我们选择了 5 行以上有 1 的列，然后计算矩阵的列和 （�������=∑我=1#����一个我�） 并仅保留带有 Z 评分的箱���2������介于 –2 和 2 之间。通过对矩阵的行和进行除以归一化，对滤波矩阵进行归一化，生成项频率矩阵��，并进一步转换为用于奇异值分解 X 的 X = 无人艇T哪里�我�=���(��我�×100000 1)×���1 #�����������.
我们还生成了一个逐个单元的矩阵 B，其中每个元素表示每个单元中每个环像素处的插补接触。这是一个 5.7 k × 3.2 M 的密集矩阵，这限制了该方法扩展到我们数据集中所有单元格的能力。我们进行了特征值分解BB型T = 通常T并按特征值从大到小对特征向量进行排序，以推导出嵌入 U 的单元。
基于域边界
我们生成了一个二元单元 x 25 kb 的 bin 矩阵，其中每个元素都指示 bin 是否被标识为单元中的域边界。使用相同的LSI框架来获得单元嵌入。
聚类分析基准
我们在每个单元格中对所有嵌入的顶部维度应用了 L2 归一化。对于 t-SNE 可视化，我们使用了前 25 个维度，但在 Higashi 和 fastHigashi 中，我们使用了前 128 个维度。采用K-means进行聚类，使用前50个维度，除Higashi和fastHigashi外，使用128个维度。从 3 到 12 枚举 k，与聚类标签相比，调整后兰德指数 （ARI） 最高的结果如图 S7 所示。为了测试分离兴奋性细胞类型的能力，我们使用了 L2/3-IT、L4-IT、L5-IT、L6-IT、L6-IT-Car3、L5/6-NP、L6b、L6-CT 和 L5-ET。对于抑制性细胞类型，我们使用了 Lamp5-Lhx6、Lamp5、Sncg、Vip、Pvalb-ChC、Pvalb 和 Sst。
无法解析细胞类型可能是由于以下生物学原因，即不同细胞类型之间的区室差异很小，或者同一细胞类型内跨细胞的区室异质性很大，或者是技术原因，即算法在单细胞Hi-C数据上识别区室的能力有限。在其他分析的基础上，我们可以确定神经元细胞类型以及与基因表达密切相关的兴奋性或抑制性亚型之间的差异区室。这表明在更精细的细胞类型之间存在伪本体水平的区室差异，而这些差异在基于区室的单细胞嵌入中无法区分。因此，单细胞异质性或计算挑战可能是主要的决定因素。染色体成像数据的分析可以帮助进一步区分这两个因素，因为它们通常被认为是染色质结构的金标准，并且不需要插补算法来判断隔室。之前的一项研究得出结论，单个细胞的区室之间存在微小差异（78）。尽管这些差异相对于跨细胞类型的差异有多大仍然难以捉摸，并且需要来自复杂组织的染色体/基因组水平染色体成像数据才能解析。总之，这一结果是生物学特征特异性和计算局限性的综合作用，用于准确量化单细胞内的特征，结论是从迄今为止我们可以应用的最佳方法中得出的，并且可能会受到技术和算法改进的挑战。
差异 3D 基因组结构与其他模式的比较车厢
原始区室评分在细胞类型之间进行分位数归一化。对于每个 100 kb 的 bin，我们使用这个标准化分数来计算其 PCC，ATAC、mCG 和 mCH 信号在不同的细胞类型中位于相同的 100 kb bin 处。我们还计算了归一化区室评分与启动子 （TSS ± 2 kb） 或基因体 （TSS – 2 kb 至 TES 2 kb） 与 100 kb bin重叠的基因表达水平之间的 PCC。
域
边界概率在每个 25 kb 箱的起始位置定义。我们在边界的上游和下游 10 kb bin 处使用 ATAC、mCG 和 mCH 信号，并取两个 bin 的平均信号来计算具有边界概率的 PCC。我们还计算了边界概率与启动子或基因体与两个 10 kb bin重叠的基因表达水平之间的 PCC。
圈
为两个 10 kb 箱之间的每个环路像素定义交互强度。我们在环路的两个锚点箱上使用了 ATAC、mCG 和 mCH 信号，并取两个箱的平均信号来计算具有插补环路强度 （Q） 的 PCC。我们还计算了启动子或基因体与两个 10 kb bin重叠，或者基因体位于两个 10 kb bins 之间的基因的环强度与表达水平之间的 PCC。
请注意，在相关性分析中，区室和结构域或环路之间的差异可能是由于分辨率不同，因为ATAC使用100 kb分辨率，甲基化可能会稀释调节元件的信号。鉴于ATAC和甲基化信号的定量分析分辨率为10 kb，因此结构域和环路更具可比性。
DEGs和与3D基因组结构的比较
由于主要类型注释之间的差异，我们根据scRNA-seq数据和snmC-seq数据之间神经元细胞的整合，根据mC细胞为每个RNA细胞分配了一个主要类型标记。鉴于两个注释之间的明确对应关系，非神经元细胞根据其原始注释进行标记。我们从每种主要类型中随机选择了 1000 个 RNA 细胞，其中选择神经元细胞的概率与标记从 mC 细胞转移到该 RNA 细胞的置信度成正比。该程序从我们 3C 分析中使用的 29 种主要类型中总共提供了 29,000 个 RNA 细胞。对于每个簇对，P值是通过Wilcoxon秩和检验得出的，倍数变化计算为两个簇中细胞的平均表达水平之间的比率。绝对值为 log2 倍变化 >1 和 FDR（BH 程序）值 <0.01 的基因被认为是差异表达的。将具有最小 FDR（BH 程序）的前 100 个 DEG 用作聚类对之间的前 DEG，并将所有可能对的顶部结果连接起来并删除重复项以生成前 DEG 的最终列表。该分析确定了常染色体上神经元主要类型之间的 1099 个顶级 DEG，以及所有主要类型之间的 1358 个 DEG。
然后，我们计算了跨神经元主要类型的 3D 基因组结构和基因表达之间的 PCC。为了避免差异分析的临界选择导致的偏差，我们根据差异统计对箱和基因进行分组，并研究了分配给每组的箱和基因的相关性（图 S11、E 和 F，以及 S12、E 和 F）。对于每个 DEG 的表达水平，我们还计算了其与每个 100 kb bins 的分位数归一化区室分数、每个位置的边界概率和 25 kb 滑动间隔或 TSS – 5 Mb 至 TES 5 Mb 区域内环的强度的相关性基因（图 2，K 到 N）。我们对每种主要类型的 3D 基因组特征和基因进行了洗牌，以计算无效 PCC 并估计 FDR。对于基因 i 和 3D 基因组特征 j 之间的每个 PCC 值（表示为 x），左侧 FDR 计算为小于 x 的洗牌 PCC 比例与观察到的 PCC 小于 x 的比例之间的比率，右侧 FDR 计算为洗牌 PCC 比例大于 x 与观察到的 PCC 比例大于 x 之间的比率.计算左侧和右侧FDR对应<0.01的PCC阈值（表示为tl和tr），最终的PCC显著性阈值确定为±{max[abs（tl）， abs（tr）]}。
CRE预测
根据在细胞亚型之间确定的成对 CH-DMG，如果与其他亚型相比，在 187 对中至少有 40 对是低甲基化 DMG，我们将一个亚型中的基因分配为低甲基化基因。如果DMR在亚型中为CG低甲基化（参见上文“确定差异甲基化区域”部分）或其CG甲基化水平为<0.3，则将DMR分配给该亚型。如果 DMR 通过差异环连接到同一亚型中也是 DMG 的基因，则 DMR 被认为是候选 CRE。我们不要求连接DMR-DMG对的差分环路是在该对的子类型中检测到的环路。这种宽松标准的原因有三点。首先，差分环的强度与甲基化水平反相关（图2K）。如果 DMR-DMG 对与一个亚型中的差分环连接，则该环可能存在于另一个亚型中，该亚型具有相似的甲基化状态的 DMR-DMG 对。其次，CRE 通常是多效性的 （79）。如果 DMR-DMG 对的甲基化状态相似，则在 m3c 数据集涵盖的亚型中检测到的环可以在另一个未发现的亚型中重用。最后，由于计算方法的局限性或某些亚型的覆盖率不足，检测过程中可能会遗漏循环。DNA循环信息在亚型之间转移，在一定程度上可以应对这种情况。
脑部疾病风险变异与跨细胞类型的 DMR 之间的关联
我们获得了与神经系统疾病和智力控制特征相关的定量特征的GWAS汇总统计 （80）、教育程度 （81）、饮酒 （82）、阿尔茨海默病 （83）、双相情感障碍 （84）、注意力缺陷多动障碍 （85）、神经质 （86）、精神分裂症 （87）、肌萎缩侧索硬化症 （88）、烟草使用障碍 （89）、失眠 （90）、睡眠时间、冠状动脉疾病 （91）、身高、疲倦 （92）、1 型糖尿病 （93）、2 型糖尿病 （94）、过敏 （95）、出生长度 （96） 和出生体重 （97）。
我们以标准格式准备了联系不平衡得分回归的汇总统计量。接下来，我们使用 LiftOver 软件将主要型低 DMR 转换为人类基因组组装 GRCh37 （hg19） 坐标，并使用 1000 基因组计划第 3 阶段 SNP 进行注释 （98）。使用hypo-DMRs的超集作为背景。最后，我们使用了细胞类型特异性连锁不平衡评分回归[https://github.com/bulik/ldsc;（40）]来估计每个性状的每个注释的富集系数。
来自DNA甲基化谱的大脑区域轴一种细胞类型的 CG 甲基化和 CH 甲基化高度可变的 100 kb bin（聚类分析的特征相同）都用于计算 PCA 的较低维表示。首先，在PCA空间中构建单元的相邻图。然后，对于每个像元，通过平均该像元及其相邻像元的位置信息来计算区域身份向量。然后从区域身份向量构建成对的曼哈顿距离矩阵，以捕获大脑区域之间的关系。将原理坐标分析应用于该距离矩阵，以获得区域空间中的低维嵌入，并保持单元间的相对距离。因此，细胞从甲基化组空间转变为区域空间。在区域空间中，我们使用STREAM（49）中实现的弹性主图算法进行轨迹分析。手动调整参数epg_alpha、epg_mu和epg_lambda，以确保生成的轨迹能够很好地代表区域空间中像元的分布。每个像元根据其与轨迹的相对位置在[0,1]中被分配一个区域索引（或伪时间）范围。然后根据其区域指数将细胞沿轨迹分组为 20 个箱。可以计算每个 bin 的平均 DNA 甲基化图谱。
皮层和基底神经节的共识区域轴
首先计算每种细胞类型中每个皮质区域的细胞的平均区域指数。然后，通过对相应像元类型的平均区域指数进行平均来计算平均区域指数。最后，通过对平均指数进行排序来构建共识区域轴。
区域 DMG
我们使用一对一对休息策略来计算皮层和基底神经节主要类型中的区域特异性 CH-DMG （rDMG）。为了避免因不同区域细胞数量不平衡而引起的潜在偏差，我们将每个区域的细胞下采样到不超过 500 个。使用 Wilcoxon 秩和检验，测试蛋白质编码基因和 lncRNA 是否显着甲基化减少（或低甲基化）。使用 BH 程序通过多测试校正调整 P 值。调整后P值的基因≤1−10和 log2 倍变化 ≤ –0.1 被认为是 rDMG。
区域 DMR
来自同一大脑区域的细胞被合并为每种主要类型，以构建区域伪本体甲基化谱。然后，使用MethylPy软件的DMRfind功能来确定候选rDMR，其选项与确定细胞型DMRs相同。如果候选 rDMR 在测试区域内的 CG 甲基化变异≥ 0.6，则将其视为 rDMR。
TF的区域丰富基序
为简单起见，我们选择了甲基化水平单调变化（PCC ≥ 0.5 或 ≤ –0.5）的 rDMR，其区域轴在皮层和基底神经节中揭示（图 4、C 和 H），并针对细胞类型特异性低甲基化 DMR 进行 TF 基序富集分析。 调整后的 P <值为 1 × 10 的 TF 基序–50被认为是富集的。
人与小鼠之间保守DMR的富集分析hcCnsvDMRs的功能富集分析
使用基因组区域注释富集工具 （GREAT） （99） 计算 hcCnsvDMRS 的基因本体 （GO） 项富集。基因关联选择“基础 延伸”选项（上游 5.0 kb，下游 1.0 kb，最大延伸 100 kb），分析中包括“精选调控结构域”。
hcCnsvDMRs与小鼠前脑组蛋白修饰标记的比较
从 Encode 项目下载 P0 小鼠前脑组蛋白修饰标记的复制峰 （4）。特别是使用了H3K27ac（ENCFF044YBD）、H3K27me3（ENCFF461UUN）、H3K4me1（ENCFF467MYU）、H3K4me3（ENCFF066LGF）和H3K9me3（ENCFF997XJK）。使用软件基因组关联测试仪（GAT，v1.3.6）（100）计算组蛋白修饰标记中hcCnsvDMR的富集。通过将 hcCnsvDMR 与 P56 小鼠大脑中分析的 snATAC-seq 峰进行比较来确定 hcCnsvDMR 的可及性 （2）。
scMCode 结构候选 CpG 站点
我们为每种主要类型构建了伪体 mCG 图谱，然后迭代选择 CpG 位点以区分所有主要类型。在每次迭代中，根据以下标准选择CpG位点：（i）在所有其余主要类型中，它们要么几乎完全甲基化（mCG%≥80%），要么未甲基化（mCG%≤20%）;（ii）两种甲基化状态均出现在其余主要类型中;（iii）CpG站点应覆盖≥10%的其余主要类型≥80%。这些选定的 CpG 站点已添加到 CpG 站点池中，以便以后构建 scMCode。在每次迭代中，其余主要类型中所选 CpG 的甲基化水平如果为 ≥80% 或 ≤20%，则被二值化。计算二元化甲基化状态之间的成对距离，并保留与任何其他主要类型距离为 <20 的细胞类型以进行下一次 CpG 选择迭代。总共有 221,140 个 CpG 站点被选为 scMCode 建设的候选者。
scMCode 的 CpG 选址
各主要类型候选CpG位点的甲基化水平基于其DNAm组分进行三元离散化（在主要类型伪本体水平下，mCG%≤20%为-1，mCG%≥80%为1,20%<mCG%<80%为0）。根据这些主要类型的离散化 DNAm 状态，将 CpG 位点进一步分为 38,945 个特征。为了防止scMCode因细胞类型群体差异或供体个体变异而产生偏差，我们从每个供体的每种主要类型中随机选择300个细胞作为scMCode构建的数据集。在每个细胞中，每个特征的甲基化状态要么通过平均属于该特征的所有 CpG 位点的甲基化水平 （AverageCpG） 来计算，要么直接使用随机选择的属于该特征的单个 CpG 位点的甲基化水平 （RandomCpG） 来计算。在单细胞水平上，根据其DNAm分数对细胞特征矩阵进行三元离散化（离散化值为-1，mCG%<50%，mCG%为1>50%，mCG%=50%或未覆盖为0），然后用于训练随机森林模型来预测主要类型。采用4重交叉验证方案防止过拟合。最后，选取前 800 个最重要的特征来构建主要类型的 scMCode。我们观察到 AverageCpG 和 RandomCpG 在预测性能方面没有差异，这表明 scMCode 的鲁棒性。
通过 KNN 插补提高预测精度
我们直接使用逐个细胞特征矩阵实现了 ~88% 的单细胞预测准确率。鉴于单细胞数据的覆盖范围有限，可以进一步插补逐个细胞特征矩阵以提高预测准确性。在训练数据集中，我们随机选择了一半的细胞，并根据它们的主要类型将它们合并为伪细胞。这个过程重复了20次，并以同样的方式从这些伪细胞构建了一个伪细胞特征矩阵。KNN插补器就是在这个矩阵上建立的。首先用KNN插补器对测试数据集进行插补，然后馈送到随机森林模型进行预测。KNN插补步骤将预测准确率提高到~93%。
交叉供体试验
我们为每个供体推导了scMCodes，并相应地训练了KNN插补器和随机森林分类器。然后将单供体scMCodes和模型应用于其他供体的数据，以评估跨个体的稳健性。当训练和测试供体相同（图6E中的对角线）时，使用4倍交叉验证方案来防止过拟合并评估准确性。
确认
该出版物得到了美国国立卫生研究院（NIH）通过推进创新神经技术（BRAIN）计划进行的大脑研究 - 细胞普查网络（BICCN）的支持和协调。本出版物是人类细胞图谱（www.humancellatlas.org/publications）的一部分。我们感谢 J. R. Dixon、T. E. Bakken 和 N. L. Jorstad 的科学讨论和有益的评论;Ecker 小组的成员进行反馈和讨论;E. Melief 和 A. Schantz 提供出色的行政援助;L. Keene、K. Kern 和 A. Keen 寻求组织收集方面的帮助;最重要的是，脑组织捐献者及其家人，没有他们，这项工作就不可能完成。
资金：这项工作得到了美国国家心理健康研究所（NIMH U01MH121282对B.R.，M.M.B.和JRE的资助;UM1MH130994 J.R.E.、M.M.B. 和 B.R.;U01MH114812 E.L. 和 S.L.）以及 Nancy 和 Buster Alvord 捐赠基金（给 C.D.K.）。索尔克研究所的流式细胞术核心设施得到了美国国立卫生研究院国家癌症研究所的资助（癌症中心支持补助金 P30 014195和共享仪器补助金 S10-OD023689）。J.R.E是霍华德休斯医学研究所的研究员。
作者贡献：J.R.E. 监督了这项研究。W.T.， J.Z.， A.B.， R.G.C.， M.K.， J.A.， A.A.， J.R.N.， H.C.， N.D.J.， J.L.， J.K.O.， A.P.， N.E.， J.R.， J.L.， M.L.， N.C.， C.O.， C.V.， A.M.Y.， J.N.， N.D.， T.C.， N.S.， D.H.， R.H.、B.P.L. 和 C.D.K. 生成了数据。W.T.、J.Z. Q.Z. 和 J.X. 分析了数据。W.T.、J.Z.、K.S.、E.L.、B.R.、M.M.B. 和 J.R.E. 对数据进行了解释。W.T.、J.Z. 和 J.R.E. 撰写了手稿。所有作者都编辑并批准了手稿。
利益争夺：J.R.E. 是 Zymo Research Inc. 和 Ionis Pharmaceuticals 的科学顾问。B.R.是Arima Genomics Inc.的联合创始人和顾问，也是Epigenome Technologies的联合创始人。WT 和 J.R.E. 已提交与本文所述研究结果相关的临时专利（美国临时申请号 63/427,789，于 2022 年 11 月 23 日提交并分配给索尔克生物研究所）。其余作者声明没有竞争利益。
数据和代码可用性：本研究中分析的数据是通过 Brain Initiative 细胞普查网络 （BICCN：RRID：SCR_015820） 生成的，并存放在 NEMO 档案 （RRID：SCR_002001） 中，标识符为 nemo：dat-jx4eu3g，可在 https://assets.nemoarchive.org/dat-jx4eu3g 访问。原始数据和处理后的数据也以GSE215353号存入NCBI GEO/SRA。单细胞甲基化的浏览器可以在 https://cellxgene.cziscience.com/collections/fdebfda9-bb9a-4b4b-97e5-651097ea07b0 找到。数据可用性的摘要可在 http://neomorph.salk.edu/hba/ 找到。映射管道可在 https://hq-1.gitbook.io/mc/;分析代码可在 https://github.com/jksr/human-brain-atlas-code 和 https://zhoujt1994.github.io/scHiCluster/hba/intro.html 获得。
许可证信息：版权所有 © 2023 作者，保留部分权利;美国科学促进会（American Association for the Advancement of Science）的独家授权人。不要求美国政府原创作品。https://www.science.org/about/science-licenses-journal-article-reuse
