文献名：Chromatin and gene-regulatory dynamics of the developing human cerebral cortex at single-cell resolution

singular value:

一 scATAC processing

使用“cellranger atac mkfastq”（10x基因组学，v.1.2.0）将原始测序数据转换为fastq格式。scATAC-seq reads与GRCh38（hg38）参考基因组，并使用“cellranger atac count”（10x Genomics, v.1.2.0）进行定量。

使用“chracr”R包（v.dev.0.9.11+）进一步处理Fragment data。我们筛选出sequencing fragments少于1000或超过50000的细胞。使用（Granja et al.，2019）中描述的方法计算TSS富集作为信噪比的度量，我们丢弃TSS富集小于4的细胞。性染色体和线粒体DNA上的片段被排除在下游分析之外。

为了获得单细胞ATAC数据集在主成分和UMAP坐标方面的低维表示，我们采用了迭代潜在语义索引(iterative latent semantic indexing)方法（Granja et al.，2019）。该方法还确定了22个细胞簇和一组共有657930个cluster peaks。简言之，在初始迭代中，根据20000个most accessible的5kb-tiling regions域确定了集群。在此，首先使用term frequency-inverse document frequency（TF-IDF）变换对计数进行归一化，并基于这些归一化计数计算singular values。使用Louvain聚类（在Seurat软件包中实现，分辨率参数=0.6）根据前25个singular values确定初始聚类，排除第一个singular values，因为它与read深度的相关系数超过0.5。然后使用MACS2（v2.1.1）在每个cluster的所有cell的aggregated insertion sites上执行Peak calling。通过从每组重叠峰中选择得分最高的峰，获得一组一致的、长度均匀的non-overlapping peaks。在第二次迭代中，其TF IDF归一化计数在初始聚类中表现出最高可变性的50000个峰值为使用前50个derived singular values的精细聚类提供了基础。在最后一次迭代中，经过refined clusters中50000个最可变的峰被确定为最终峰集，并再次计算singular values。UMAP坐标和ATAC簇是根据这些最终singular values的前10个确定的。使用“uwot”R软件包中实现的UMAP生成二维表示（v.0.1.8; parameter settings: ‘min.dist = 0.6’, ‘n.neighbors = 50’, ‘cosine’ distance metric)。

ChromVAR（v.1.6）用于使用JASPAR 2018数据库中的位置权重矩阵获得TFaccessibility profiles。使用“ChrAccR”将Gene activity scores计算为TSS相关峰的aggregated accessibility。为此，使用宽度参数sigma=10000 bp的径向基函数（radial basis function RBF）分配的权重，将TSS 100000 bp内的peak counts相加，将最小渐近权重(minimum asymptotic weight)设置为0.25。对于每个基因，所得分数通过权重之和标准化。对于可视化和下游分析，单个细胞的计数已重新调整为10000计数，并已进行log2标准化。为了增强二维UMAP空间中的可视化效果，利用奇异值空间中确定的细胞邻域，使用MAGIC diffusion algorithm（van Dijk等人，2018）平滑了(smoothed )基因活性分数。由于此类imputation方法与 risk of oversmoothing相关（，我们限制了MAGIC在数据可视化中的应用。

我们通过在200bp基因组平铺窗口中对簇pseudobulk samples的insertion counts求和来创建ATAC信号轨迹，并提供与WashU表观基因组浏览器兼容的trackhub(http://epigenomegateway.wustl.edu)除了推测的CRE-gene links外，还包含这些profiles。

Matching of single-cell transcriptomes and epigenomes
Seurat实施的典型相关分析（CCA）已分别应用于每个妊娠时间点的匹配单细胞RNA和ATAC数据。为此，我们计算log-normalized and scaled gene activity scores，作为scATAC-seq分析的细胞中基因表达的替代物。作为整合特征，我们使用每种模式中2000个最可变基因的结合作为Seurat的“FindTransferAnchors”功能的input，使用reduction method “cca”和参数“k.anchor=10”。对于scRNA-seq分析的每个细胞和scATAC-seq分析的每个细胞，我们通过在联合CCA L2空间中应用最近邻搜索，在各自的其他模式中识别最近邻细胞。使用“FNN”R包确定最近邻，使用带有欧氏距离的“kd_tree”算法。这些来自所有妊娠时间点的基于最近邻的细胞匹配被连接起来，以获得跨两种模式的数据集范围的细胞匹配。

Linking gene regulatory elements and gene expression across all cell types
我们使用基于相关性的方法识别了peak-to-gene links，该方法应用于聚集scATAC和scRNA计数的pseudobulk samples 。通过从整个scATAC seq数据集中随机抽取200个细胞来定义这些pseudobulk samples 。将这200个种子细胞与其各自的99个最近邻细胞在ATAC-PC空间中组合，使得每个pseudobulk samples 总共包含100个细胞。峰的pseudobulk ATAC insertion counts通过对各单细胞成员的峰插入计数求和获得。通过选择与CCA空间中的100个ATAC细胞相似的最近邻的100个scRNA细胞，获得匹配的RNA细胞。通过对各单细胞成员的基因计数求和获得pseudobulk RNA基因计数。类似地，在多组数据集中，从ATAC模式中采集了100个细胞的200个pseudobulk samples，并在RNA空间聚集相同的细胞。每个匹配的pseudobulk samples分别用其或有RNA和ATAC细胞的多数簇和年龄分配进行注释。

然后，我们通过将基因组距离在1到250kb之间的峰与蛋白质编码的TSS关联，并将lincRNA基因与相应基因关联，获得候选峰基因对。对于每个候选峰值基因对，我们计算可及性和基因表达数据的CPM标准化计数的Pearson相关系数，并根据其t统计量计算这些系数的FDR调整P值。我们通过仅保留| PCC |>0.4和FDR调整的P值<0.05的配对，定义了一组64878个高置信peak-to-gene links。使用相同的方法，为多组数据获得了一组相应的76374条links。推断的和multiome peak-gene links之间的重叠是通过为每个链接创建“‘GenomicInteraction”对象来计算的，peak作为第一个锚，基因启动子作为第二个锚，然后应用带有参数“use.region = ‘both’”的函数“findOverlaps”

Validation of inferred peak-gene links using conservation and chromosome conformation capture data

为了使用orthogonal分析验证上述linkages，采用了两种方法。首先，对于multiome and singleome linkages，使用“GenomicScore”软件包中的“gscores”函数，计算linked and unlinked peaks的phastCons 100-way vertebrate conservation scores。使用 Wilcoxon rank-sum test比较linked and unliked peaks的得分。

其次，我们使用最近发布的邻近连接辅助芯片测序（PLAC-seq）数据集对3D接触进行了分析，该数据集针对大脑皮层4种分类发育人类细胞类型中的H3K4me3位点（Song等人，2020年）。数据集由来自FACS分类的中间神经元、兴奋性神经元、放射状胶质细胞和来自分离的人类大脑皮层组织的中间祖细胞的启动子捕获3D接触文库组成。因此，这个orthogonal 3D contact dataset提供了两个验证轴：第一，增强子-启动子linkages，第二，这些linkages的细胞类型特异性。

来自PLAC-seq数据的Interaction calls作为“GenomicInteractions”对象导入，并与我们的linkages overlap（“FindVerlaps”）。为了验证，both interactions 的两个锚都需要重叠。我们对所有可能的peak-gene links、显著推断的peak-gene links进行了分析，并且，由于这些正交数据类型不符合1:1，对于independent test，我们也在overlap分析之前预筛选了与PLAC-seq区域的任何一维相互作用的significant links。

为了解释significant links的skewed length distribution，我们还从所有可能链接的空间（每个链接10000个）生成了1000个长度匹配permutations。首先，对于significant peak-gene links，计算peak-promoter distance。将距离分为25个0-250kb的等分箱，并计算每个箱中peak-gene links的比例。接下来，我们将所有可能的peak-gene links分配给一个bin和真实分布中相应的比例。比例被用作绘制排列的抽样概率。然后，根据该length-matched null model计算PLAC-seq overlaps。

最后，我们推断，如果inferred linkages是有效的，验证的基因也会表现出细胞类型限制性表达模式，与3D contact的分类细胞类型一致。为了确定这一点，我们计算了RNA-seq数据中主要细胞类型中linked genes的表达。然后，我们根据PLAC-seq相互作用的细胞类型来源来划分这些表达值。同样，对于linkages的ATAC-seq峰值，我们计算了相同边界上scATAC-seq数据的mean accessibility。

Identification of genes with predictive chromatin (GPCs)
GPC的定义主要基于单个细胞之间的high gene activity-expression correlations。为了使这一分析对技术变异更具鲁棒性，我们将分析局限于背侧前脑细胞中最可变的基因（1999个基因）。具体而言，我们使用了URD包中的“findVariableGenes”函数，参数为“diffCV.cutoff=0.15，mean.min=0.004”。对于每个可变基因，我们计算ATAC细胞的基因活性得分向量与RNA数据中相应最近邻细胞的表达得分向量之间的Spearman相关系数。我们还将这些相关性与每个基因的linked enhancers per gene进行了比较。从这个子集中，我们将GPC定义为与至少10个CRE相关的 前10%基因活性表达相关性中的基因。

Calculation of motif synergy and correlation scores

我们使用“getAnnotationSynergy”chromVAR函数计算motif簇之间的pairwise synergy scores。这些分数定义为包含两个不同motif簇的结合位点的CREs中染色质活性的差异，以及随机子样本CREs中的可达性差异，该随机子样本CREs仅包含一个motif簇的结合位点（差异较大的基序簇）。这一定义基于这样一种直觉，即与只有一个TF可以结合的基因座相比，两个TF可以潜在结合的基因组基因座的可及性更高的动态性（方差）暗示了TFs的潜在共同依赖性。因此，正协同分数对应于潜在共结合所解释的可及性可变性超过独立基序发生所解释的可变性的相互作用。为了区分motif accessibility中的co-dependence概念和简单相关性，我们还使用chromVAR中的“getAnnotationCorrelation”函数计算了相关系数。该分数定义为分别从包含一个但不包含另一个基序簇的基序的CRE计算的aggregate motif activity scores（偏差分数）之间的相关性。

Fuzzy c-means: clustering and re-projection approach

对于fuzzy clustering analysis，首先从胶质细胞簇（单个细胞的10%）中随机选择1267个种子细胞，选择的数量与簇大小成比例。通过将这些细胞与其scRNA PCA空间中的50个最近邻结合，对Pseudobulk数据集进行Pseudobulk。接下来，使用R软件包“URD”中的功能“findVariableGenes”确定1957个可变表达基因。通过对构成每个Pseudobulk的各个单细胞成员的特征计数求和，形成pseudobulk counts matrix 。

使用R软件包“e1071”中的函数“cmeans”对该pseudobulk matrix 进行pseudobulk matrix 聚类，参数c=14，m=1.25，产生了一个按模块划分的基因“membership matrix”和一个按模块划分的样本“centers matrix”。为了确定下游分析的“fixed”或binarized module membership，我们将threshold membership score定义为将所有基因分配给一个簇的最大得分（阈值=0.06）。使用R包“clusterProfiler”中的函数“enrichGO”计算每个模块的基因本体丰富度。使用Jaccard指数计算所有模块对之间的模块连接性，并通过应用基因共享的Jaccard指数的0.2阈值连接模块。通过应用肘部法选择该阈值。为了可视化模块之间的连接，使用R软件包“UMAP”，使用中心矩阵（逐个模块采样）作为UMAP降维的基础。

最后，重复该过程，将聚类参数（c，m）和membership threshold across a range of values；从c=6到c=30，从m=1到2；以确保生成的嵌入结构不会对聚类参数过于敏感。

Projecting ATAC-seq data into fuzzy clustering space and GPC projection

scATAC细胞的Pseudobulk samples是使用上述基因活性评分方法生成的。该矩阵被子集以匹配RNA模糊聚类分析的特征（基因）。在缺少特征的情况下，使用其中间基因活性估算值。为了将ATAC-seq细胞投射到RNA模糊聚类嵌入中，我们转置了membership matrix，并用伪块矩阵将其与基因活性相乘。最后，我们使用R“stats”中的“predict”函数，模糊聚类UMAP模型作为第一个参数，得到的转置乘积矩阵作为第二个参数，以确定ATAC伪块的UMAP坐标。

为了执行投影操作，我们获取样本X基因活性分数矩阵，并将其乘以来自模糊C-均值聚类的特征loadings（genes X loadings）。生成的矩阵被输入用于创建original manifold.的同一UMAP模型。这提供了流形上投影点和地标之间相似性的可视化。为了将此操作限制在GPC中，我们强制其他基因的基因活性分数为中位数。因此，同样的样本被预测两次，一次考虑所有基因，另一次只考虑GPC。这两点是图6F中arrow visualization的基础。最后，为了提供该转化的基线，我们使用随机和定义的对照基因集执行了该操作。