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

singular value:

一 scATAC processing

使用“cellranger atac mkfastq”（10x基因組學(xué)，v.1.2.0）將原始測序數(shù)據(jù)轉(zhuǎn)換為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深度的相關(guān)系數(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的精細聚類提供了基礎(chǔ)。在最后一次迭代中健民，經(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ù)庫中的位置權(quán)重矩陣獲得TFaccessibility profiles。使用“ChrAccR”將Gene activity scores計算為TSS相關(guān)峰的aggregated accessibility崇堵。為此型诚，使用寬度參數(shù)sigma=10000 bp的徑向基函數(shù)（radial basis function RBF）分配的權(quán)重，將TSS 100000 bp內(nèi)的peak counts相加鸳劳，將最小漸近權(quán)重(minimum asymptotic weight)設(shè)置為0.25狰贯。對于每個基因，所得分數(shù)通過權(quán)重之和標準化赏廓。對于可視化和下游分析涵紊，單個細胞的計數(shù)已重新調(diào)整為10000計數(shù)，并已進行l(wèi)og2標準化楚昭。為了增強二維UMAP空間中的可視化效果栖袋，利用奇異值空間中確定的細胞鄰域，使用MAGIC diffusion algorithm（van Dijk等人抚太，2018）平滑了(smoothed )基因活性分數(shù)。由于此類imputation方法與 risk of oversmoothing相關(guān)（昔案，我們限制了MAGIC在數(shù)據(jù)可視化中的應(yīng)用尿贫。

我們通過在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實施的典型相關(guān)分析（CCA）已分別應(yīng)用于每個妊娠時間點的匹配單細胞RNA和ATAC數(shù)據(jù)。為此捞稿，我們計算log-normalized and scaled gene activity scores又谋，作為scATAC-seq分析的細胞中基因表達的替代物拼缝。作為整合特征，我們使用每種模式中2000個最可變基因的結(jié)合作為Seurat的“FindTransferAnchors”功能的input彰亥，使用reduction method “cca”和參數(shù)“k.anchor=10”咧七。對于scRNA-seq分析的每個細胞和scATAC-seq分析的每個細胞，我們通過在聯(lián)合CCA L2空間中應(yīng)用最近鄰搜索任斋，在各自的其他模式中識別最近鄰細胞继阻。使用“FNN”R包確定最近鄰，使用帶有歐氏距離的“kd_tree”算法废酷。這些來自所有妊娠時間點的基于最近鄰的細胞匹配被連接起來瘟檩，以獲得跨兩種模式的數(shù)據(jù)集范圍的細胞匹配。

Linking gene regulatory elements and gene expression across all cell types
我們使用基于相關(guān)性的方法識別了peak-to-gene links澈蟆，該方法應(yīng)用于聚集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關(guān)聯(lián)喂击，并將lincRNA基因與相應(yīng)基因關(guān)聯(lián)剂癌，獲得候選峰基因?qū)Ατ诿總€候選峰值基因?qū)舶恚覀冇嬎憧杉靶院突虮磉_數(shù)據(jù)的CPM標準化計數(shù)的Pearson相關(guān)系數(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ù)獲得了一組相應(yīng)的76374條links裁奇。推斷的和multiome peak-gene links之間的重疊是通過為每個鏈接創(chuàng)建“‘GenomicInteraction”對象來計算的桐猬，peak作為第一個錨，基因啟動子作為第二個錨刽肠，然后應(yīng)用帶有參數(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”對象導(dǎo)入汽抚，并與我們的linkages overlap（“FindVerlaps”）。為了驗證伯病，both interactions 的兩個錨都需要重疊造烁。我們對所有可能的peak-gene links、顯著推斷的peak-gene links進行了分析午笛，并且惭蟋，由于這些正交數(shù)據(jù)類型不符合1:1，對于independent test药磺，我們也在overlap分析之前預(yù)篩選了與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和真實分布中相應(yīng)的比例。比例被用作繪制排列的抽樣概率卿叽。然后桥胞，根據(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。為了使這一分析對技術(shù)變異更具魯棒性询刹，我們將分析局限于背側(cè)前腦細胞中最可變的基因（1999個基因）谜嫉。具體而言，我們使用了URD包中的“findVariableGenes”函數(shù)凹联，參數(shù)為“diffCV.cutoff=0.15沐兰，mean.min=0.004”。對于每個可變基因蔽挠，我們計算ATAC細胞的基因活性得分向量與RNA數(shù)據(jù)中相應(yīng)最近鄰細胞的表達得分向量之間的Spearman相關(guān)系數(shù)住闯。我們還將這些相關(guān)性與每個基因的linked enhancers per gene進行了比較。從這個子集中象泵，我們將GPC定義為與至少10個CRE相關(guān)的 前10%基因活性表達相關(guān)性中的基因寞秃。

Calculation of motif synergy and correlation scores

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

Fuzzy c-means: clustering and re-projection approach

對于fuzzy clustering analysis，首先從膠質(zhì)細胞簇（單個細胞的10%）中隨機選擇1267個種子細胞，選擇的數(shù)量與簇大小成比例耳璧。通過將這些細胞與其scRNA PCA空間中的50個最近鄰結(jié)合成箫，對Pseudobulk數(shù)據(jù)集進行Pseudobulk。接下來旨枯，使用R軟件包“URD”中的功能“findVariableGenes”確定1957個可變表達基因蹬昌。通過對構(gòu)成每個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ù)計算所有模塊對之間的模塊連接性，并通過應(yīng)用基因共享的Jaccard指數(shù)的0.2閾值連接模塊败潦。通過應(yīng)用肘部法選擇該閾值本冲。為了可視化模塊之間的連接，使用R軟件包“UMAP”劫扒，使用中心矩陣（逐個模塊采樣）作為UMAP降維的基礎(chǔ)檬洞。

最后，重復(fù)該過程沟饥，將聚類參數(shù)（c添怔，m）和membership threshold across a range of values；從c=6到c=30贤旷，從m=1到2广料；以確保生成的嵌入結(jié)構(gòu)不會對聚類參數(shù)過于敏感。

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

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

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