文獻名：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）將原始測序數(shù)據(jù)轉換為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數(shù)據(jù)集在主成分和UMAP坐標方面的低維表示，我們采用了迭代潛在語義索引(iterative latent semantic indexing)方法（Granja et al.，2019）。該方法還確定了22個細胞簇和一組共有657930個cluster peaks。簡言之，在初始迭代中，根據(jù)20000個most accessible的5kb-tiling regions域確定了集群。在此，首先使用term frequency-inverse document frequency（TF-IDF）變換對計數(shù)進行歸一化，并基于這些歸一化計數(shù)計算singular values。使用Louvain聚類（在Seurat軟件包中實現(xiàn)，分辨率參數(shù)=0.6）根據(jù)前25個singular values確定初始聚類，排除第一個singular values，因為它與read深度的相關系數(shù)超過0.5。然后使用MACS2（v2.1.1）在每個cluster的所有cell的aggregated insertion sites上執(zhí)行Peak calling。通過從每組重疊峰中選擇得分最高的峰，獲得一組一致的、長度均勻的non-overlapping peaks。在第二次迭代中，其TF IDF歸一化計數(shù)在初始聚類中表現(xiàn)出最高可變性的50000個峰值為使用前50個derived singular values的精細聚類提供了基礎。在最后一次迭代中，經(jīng)過refined clusters中50000個最可變的峰被確定為最終峰集，并再次計算singular values。UMAP坐標和ATAC簇是根據(jù)這些最終singular values的前10個確定的。使用“uwot”R軟件包中實現(xiàn)的UMAP生成二維表示（v.0.1.8; parameter settings: ‘min.dist = 0.6’, ‘n.neighbors = 50’, ‘cosine’ distance metric)。

ChromVAR（v.1.6）用于使用JASPAR 2018數(shù)據(jù)庫中的位置權重矩陣獲得TFaccessibility profiles。使用“ChrAccR”將Gene activity scores計算為TSS相關峰的aggregated accessibility。為此，使用寬度參數(shù)sigma=10000 bp的徑向基函數(shù)（radial basis function RBF）分配的權重，將TSS 100000 bp內(nèi)的peak counts相加，將最小漸近權重(minimum asymptotic weight)設置為0.25。對于每個基因，所得分數(shù)通過權重之和標準化。對于可視化和下游分析，單個細胞的計數(shù)已重新調(diào)整為10000計數(shù)，并已進行l(wèi)og2標準化。為了增強二維UMAP空間中的可視化效果，利用奇異值空間中確定的細胞鄰域，使用MAGIC diffusion algorithm（van Dijk等人，2018）平滑了(smoothed )基因活性分數(shù)。由于此類imputation方法與 risk of oversmoothing相關（，我們限制了MAGIC在數(shù)據(jù)可視化中的應用。

我們通過在200bp基因組平鋪窗口中對簇pseudobulk samples的insertion counts求和來創(chuàng)建ATAC信號軌跡，并提供與WashU表觀基因組瀏覽器兼容的trackhub(http://epigenomegateway.wustl.edu)除了推測的CRE-gene links外，還包含這些profiles。

Matching of single-cell transcriptomes and epigenomes
Seurat實施的典型相關分析（CCA）已分別應用于每個妊娠時間點的匹配單細胞RNA和ATAC數(shù)據(jù)。為此，我們計算log-normalized and scaled gene activity scores，作為scATAC-seq分析的細胞中基因表達的替代物。作為整合特征，我們使用每種模式中2000個最可變基因的結合作為Seurat的“FindTransferAnchors”功能的input，使用reduction method “cca”和參數(shù)“k.anchor=10”。對于scRNA-seq分析的每個細胞和scATAC-seq分析的每個細胞，我們通過在聯(lián)合CCA L2空間中應用最近鄰搜索，在各自的其他模式中識別最近鄰細胞。使用“FNN”R包確定最近鄰，使用帶有歐氏距離的“kd_tree”算法。這些來自所有妊娠時間點的基于最近鄰的細胞匹配被連接起來，以獲得跨兩種模式的數(shù)據(jù)集范圍的細胞匹配。

Linking gene regulatory elements and gene expression across all cell types
我們使用基于相關性的方法識別了peak-to-gene links，該方法應用于聚集scATAC和scRNA計數(shù)的pseudobulk samples 。通過從整個scATAC seq數(shù)據(jù)集中隨機抽取200個細胞來定義這些pseudobulk samples 。將這200個種子細胞與其各自的99個最近鄰細胞在ATAC-PC空間中組合，使得每個pseudobulk samples 總共包含100個細胞。峰的pseudobulk ATAC insertion counts通過對各單細胞成員的峰插入計數(shù)求和獲得。通過選擇與CCA空間中的100個ATAC細胞相似的最近鄰的100個scRNA細胞，獲得匹配的RNA細胞。通過對各單細胞成員的基因計數(shù)求和獲得pseudobulk RNA基因計數(shù)。類似地，在多組數(shù)據(jù)集中，從ATAC模式中采集了100個細胞的200個pseudobulk samples，并在RNA空間聚集相同的細胞。每個匹配的pseudobulk samples分別用其或有RNA和ATAC細胞的多數(shù)簇和年齡分配進行注釋。

然后，我們通過將基因組距離在1到250kb之間的峰與蛋白質(zhì)編碼的TSS關聯(lián)，并將lincRNA基因與相應基因關聯(lián)，獲得候選峰基因對。對于每個候選峰值基因對，我們計算可及性和基因表達數(shù)據(jù)的CPM標準化計數(shù)的Pearson相關系數(shù)，并根據(jù)其t統(tǒng)計量計算這些系數(shù)的FDR調(diào)整P值。我們通過僅保留| PCC |>0.4和FDR調(diào)整的P值<0.05的配對，定義了一組64878個高置信peak-to-gene links。使用相同的方法，為多組數(shù)據(jù)獲得了一組相應的76374條links。推斷的和multiome peak-gene links之間的重疊是通過為每個鏈接創(chuàng)建“‘GenomicInteraction”對象來計算的，peak作為第一個錨，基因啟動子作為第二個錨，然后應用帶有參數(shù)“use.region = ‘both’”的函數(shù)“findOverlaps”

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

為了使用orthogonal分析驗證上述linkages，采用了兩種方法。首先，對于multiome and singleome linkages，使用“GenomicScore”軟件包中的“gscores”函數(shù)，計算linked and unlinked peaks的phastCons 100-way vertebrate conservation scores。使用 Wilcoxon rank-sum test比較linked and unliked peaks的得分。

其次，我們使用最近發(fā)布的鄰近連接輔助芯片測序（PLAC-seq）數(shù)據(jù)集對3D接觸進行了分析，該數(shù)據(jù)集針對大腦皮層4種分類發(fā)育人類細胞類型中的H3K4me3位點（Song等人，2020年）。數(shù)據(jù)集由來自FACS分類的中間神經(jīng)元、興奮性神經(jīng)元、放射狀膠質(zhì)細胞和來自分離的人類大腦皮層組織的中間祖細胞的啟動子捕獲3D接觸文庫組成。因此，這個orthogonal 3D contact dataset提供了兩個驗證軸：第一，增強子-啟動子linkages，第二，這些linkages的細胞類型特異性。

來自PLAC-seq數(shù)據(jù)的Interaction calls作為“GenomicInteractions”對象導入，并與我們的linkages overlap（“FindVerlaps”）。為了驗證，both interactions 的兩個錨都需要重疊。我們對所有可能的peak-gene links、顯著推斷的peak-gene links進行了分析，并且，由于這些正交數(shù)據(jù)類型不符合1:1，對于independent test，我們也在overlap分析之前預篩選了與PLAC-seq區(qū)域的任何一維相互作用的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和真實分布中相應的比例。比例被用作繪制排列的抽樣概率。然后，根據(jù)該length-matched null model計算PLAC-seq overlaps。

最后，我們推斷，如果inferred linkages是有效的，驗證的基因也會表現(xiàn)出細胞類型限制性表達模式，與3D contact的分類細胞類型一致。為了確定這一點，我們計算了RNA-seq數(shù)據(jù)中主要細胞類型中l(wèi)inked genes的表達。然后，我們根據(jù)PLAC-seq相互作用的細胞類型來源來劃分這些表達值。同樣，對于linkages的ATAC-seq峰值，我們計算了相同邊界上scATAC-seq數(shù)據(jù)的mean accessibility。

Identification of genes with predictive chromatin (GPCs)
GPC的定義主要基于單個細胞之間的high gene activity-expression correlations。為了使這一分析對技術變異更具魯棒性，我們將分析局限于背側前腦細胞中最可變的基因（1999個基因）。具體而言，我們使用了URD包中的“findVariableGenes”函數(shù)，參數(shù)為“diffCV.cutoff=0.15，mean.min=0.004”。對于每個可變基因，我們計算ATAC細胞的基因活性得分向量與RNA數(shù)據(jù)中相應最近鄰細胞的表達得分向量之間的Spearman相關系數(shù)。我們還將這些相關性與每個基因的linked enhancers per gene進行了比較。從這個子集中，我們將GPC定義為與至少10個CRE相關的 前10%基因活性表達相關性中的基因。

Calculation of motif synergy and correlation scores

我們使用“getAnnotationSynergy”chromVAR函數(shù)計算motif簇之間的pairwise synergy scores。這些分數(shù)定義為包含兩個不同motif簇的結合位點的CREs中染色質(zhì)活性的差異，以及隨機子樣本CREs中的可達性差異，該隨機子樣本CREs僅包含一個motif簇的結合位點（差異較大的基序簇）。這一定義基于這樣一種直覺，即與只有一個TF可以結合的基因座相比，兩個TF可以潛在結合的基因組基因座的可及性更高的動態(tài)性（方差）暗示了TFs的潛在共同依賴性。因此，正協(xié)同分數(shù)對應于潛在共結合所解釋的可及性可變性超過獨立基序發(fā)生所解釋的可變性的相互作用。為了區(qū)分motif accessibility中的co-dependence概念和簡單相關性，我們還使用chromVAR中的“getAnnotationCorrelation”函數(shù)計算了相關系數(shù)。該分數(shù)定義為分別從包含一個但不包含另一個基序簇的基序的CRE計算的aggregate motif activity scores（偏差分數(shù)）之間的相關性。

Fuzzy c-means: clustering and re-projection approach

對于fuzzy clustering analysis，首先從膠質(zhì)細胞簇（單個細胞的10%）中隨機選擇1267個種子細胞，選擇的數(shù)量與簇大小成比例。通過將這些細胞與其scRNA PCA空間中的50個最近鄰結合，對Pseudobulk數(shù)據(jù)集進行Pseudobulk。接下來，使用R軟件包“URD”中的功能“findVariableGenes”確定1957個可變表達基因。通過對構成每個Pseudobulk的各個單細胞成員的特征計數(shù)求和，形成pseudobulk counts matrix 。

使用R軟件包“e1071”中的函數(shù)“cmeans”對該pseudobulk matrix 進行pseudobulk matrix 聚類，參數(shù)c=14，m=1.25，產(chǎn)生了一個按模塊劃分的基因“membership matrix”和一個按模塊劃分的樣本“centers matrix”。為了確定下游分析的“fixed”或binarized module membership，我們將threshold membership score定義為將所有基因分配給一個簇的最大得分（閾值=0.06）。使用R包“clusterProfiler”中的函數(shù)“enrichGO”計算每個模塊的基因本體豐富度。使用Jaccard指數(shù)計算所有模塊對之間的模塊連接性，并通過應用基因共享的Jaccard指數(shù)的0.2閾值連接模塊。通過應用肘部法選擇該閾值。為了可視化模塊之間的連接，使用R軟件包“UMAP”，使用中心矩陣（逐個模塊采樣）作為UMAP降維的基礎。

最后，重復該過程，將聚類參數(shù)（c，m）和membership threshold across a range of values；從c=6到c=30，從m=1到2；以確保生成的嵌入結構不會對聚類參數(shù)過于敏感。

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

scATAC細胞的Pseudobulk samples是使用上述基因活性評分方法生成的。該矩陣被子集以匹配RNA模糊聚類分析的特征（基因）。在缺少特征的情況下，使用其中間基因活性估算值。為了將ATAC-seq細胞投射到RNA模糊聚類嵌入中，我們轉置了membership matrix，并用偽塊矩陣將其與基因活性相乘。最后，我們使用R“stats”中的“predict”函數(shù)，模糊聚類UMAP模型作為第一個參數(shù)，得到的轉置乘積矩陣作為第二個參數(shù)，以確定ATAC偽塊的UMAP坐標。

為了執(zhí)行投影操作，我們獲取樣本X基因活性分數(shù)矩陣，并將其乘以來自模糊C-均值聚類的特征loadings（genes X loadings）。生成的矩陣被輸入用于創(chuàng)建original manifold.的同一UMAP模型。這提供了流形上投影點和地標之間相似性的可視化。為了將此操作限制在GPC中，我們強制其他基因的基因活性分數(shù)為中位數(shù)。因此，同樣的樣本被預測兩次，一次考慮所有基因，另一次只考慮GPC。這兩點是圖6F中arrow visualization的基礎。最后，為了提供該轉化的基線，我們使用隨機和定義的對照基因集執(zhí)行了該操作。