小鼠初级运动皮层的转录组和表观基因组细胞图集

　　以下方法适用于艾伦学院生成的以下转录组数据集：SCRNA SMART，SCRNA 10X V3 A，SCRNA 10X V2 A，SNRNA SMART，SNRNA 10X V3 A和SNRNA 10X V2 A. 　　所有程序均根据艾伦脑科学研究所的机构动物护理和使用委员会协议进行。为小鼠提供食物和水，并在每天不超过五只成年小鼠的日常12小时/夜周期中维持水和水。环境温度设置为72°F，相对湿度设置为40％。所有房间都在12/12-H的光线/黑暗周期上。在这项研究中，我们使用SNAP25-IRES2-CRE MICE41（MGI：J：220523）富含神经元，越过AI14（参考42）（MGI：J：220523），这些（MGI：J：220523）维持在C57BL/6J背景上（RRID：RRID：IMSR_JAX：IMSR_JAX：000664）。小鼠在产后53-59天安乐死。从男性和女性（Scrna Smart，SnRNA智能，SCRNA 10x V3 A和SnRNA 10x V2 A），仅雄性（SCRNA 10x V2 A）或仅女性（SNRNA 10x V3 A）收集组织。 　　我们通过调整先前描述的程序3,43来隔离单细胞。将大脑解剖，浸入人造脑脊液（ACSF）3中，嵌入2％琼脂糖中，并将其切成250-μm（智能seq）或350-μm（10x基因组学）冠状动脉剖面（精确仪器）。Allen Mouse Brain公共坐标框架版本3（CCFV3; RRID：SCR_002978）44本体学用于定义解剖的MOP（扩展数据图1B）。 　　对于Smart-Seq，将MOP从切片中进行了微切除，并用1 mg/ml pronase（p6911-1g，Sigma）分解为单个细胞，并如前所述处理3。对于10倍基因组学，在ACSF中用30 U/ml Papain（Pap2，Worthington）在30°C中消化组织块。通过将木瓜蛋白酶溶液与猝灭缓冲液（ACSF为1％FBS和0.2％BSA）交换三次，从而淬灭酶消化。将淬灭缓冲液中的组织块通过火抛光的移液器进行了三杆，直径约20次，直径为600μm。允许溶液沉降，并将含有单个细胞的上清液转移到新管中。将新鲜的淬灭缓冲液添加到沉降的组织块中，并使用300μm和150-μm的火抛光移液器重复进行三次和上清液转移。单细胞悬浮液通过70μm的滤光片进入15 ml圆锥管，底部有500μl高BSA缓冲液（ACSF为1％FBS和1％BSA），在摇摆桶中心的1分钟内，在100g的100g中离心10分钟。将上清液丢弃，并将细胞颗粒重悬于淬火缓冲液中。 　　使用130-μm喷嘴通过荧光激活的细胞分选（FACS； BD ARIA II； RRID：SCR_018091）收集所有细胞。通过将悬浮液通过70μm滤波器传递并添加DAPI（最终浓度为2 ng/ml）来制备细胞进行分类。排序策略如前所述3，大多数细胞使用TDTomato阳性标签收集。对于Smart-Seq，将单个细胞分为八孔PCR条的单个孔，这些井中含有来自Smart-Seq V4超低输入RNA试剂盒中的裂解缓冲液，用于测序（634894，Takara），用RNase抑制剂（0.17 U/μL），立即在干冰上冷冻，并在-80°C上储存在干冰上，并立即在-80°C上释放。对于10倍基因组学，将30,000个细胞在10分钟内分为含有500μl淬火缓冲液的管。每个等分试样的30,000个排序细胞的每个等分试样在200μl高BSA缓冲液的顶部轻轻地分层，并在摇摆的桶离心机中立即以230克离心10分钟。除去上清液，并将35μl缓冲液留在后面，其中重悬于细胞颗粒。定量细胞浓度并立即加载到10倍基因组铬控制器上。 　　这些方法适用于Broad Institute生成的SNRNA 10x V3 B数据集。 　　小鼠被安装在12小时的轻/黑暗时间表中，并在抵达后2周适应其住房环境。环境温度设置为70±2°F，相对湿度设置为40±10％。所有房间都在12/12-H的光线/黑暗周期上。马萨诸塞州理工学院涉及动物的所有程序均根据美国国立卫生研究院指南，以根据协议编号1115-111-18进行护理和使用实验动物的护理和使用指南，并由马萨诸塞州动物护理技术委员会批准。根据美国国立卫生研究院指南，针对实验动物的护理和使用规定，根据规程编号0120-09-16，都根据美国国立卫生研究院指南进行了所有涉及动物的程序。从雄性和雌性小鼠收集样品。 　　在60天龄的时候，通过在流动3％异氟烷的燃气室中给药1分钟，C57BL/6J小鼠通过给药来麻醉。通过检查负尾捏反应来确认麻醉。将小鼠移至解剖托盘中，并在过程的持续时间内通过流动3％异氟烷的鼻锥延长麻醉。使用含有110 mM NaCl，10 mM HEPES，25 mm葡萄糖，75 mM蔗糖，7.5 mM MGCL2和2.5 mM KCl的冰冷的pH 7.4 HEPES缓冲液进行心病灌注，以清除大脑和其他器官的血液。将大脑立即去除并在液氮蒸气中冷冻3分钟，并移至-80°C以进行长期储存。可以在协议中获得详细的协议。IO15。 　　冷冻小鼠的大脑被小脑牢固地安装在带有OCT嵌入化合物的低温恒温器上，以使包括拖把在内的整个前半部分被暴露并热不受干扰。使用眼托（P-715，P-715，羽毛安全剃须刀）在低温恒温器中手动对MOP的500μM前跨跨度（扩展数据）进行解剖（扩展数据图1C），并将其预先冷却至-20°C，并穿上了4×手术的Loupes。使用预冷的镊子将每个切除的组织解剖置于预冷的0.25-mL PCR管中，并存储在-80°C下。为了评估解剖精度，在每条500μm前后解剖连接处进行10μm冠状切片，并在NISSL染色后成像。根据根据McCarroll实验室慷慨提供的协议，使用柔和的，基于洗涤剂的解离从这些冷冻组织中提取核，该核可以使用柔和的，洗涤剂的解离，并加载到10倍Chromium V3系统中。逆转录和图书馆的生成是根据制造商的协议进行的。 　　与其他SNRNA-SEQ方法相比，该10x V3 SNRNA-SEQ方案导致恢复的基因数量更高。我们认为，这有三个原因，并且由这些结果结合起来所赋予的福利总和。 　　首先，将小鼠的大脑灌注于模拟ACSF的溶液中，然后在液氮蒸气上快速冷冻，以使RNA完整性高度保存。最初的脑组织的生物分析仪RIN得分通常为9.8。在解剖前的储存在1 ml OCT化合物中的水合下水槽的存在下，将大脑的储存在-80°C下，预冻结到5 ml储存管的底部。这样可以防止升华和随后的干燥依赖性RNA碎片化。 　　其次，我们进行了快速的样本处理。我们有一组训练有素的技术人员（如上所述）处理小鼠大脑，然后在一个连续方案中执行分离，FACS和10X处理（如下），而无需停顿。例如，每只鼠标都被灌注并准备在几分钟（10分钟）内进行解剖，我们将样本量限制在六只小鼠中，因此没有样品在等待整个过程中移动。 　　第三，冷冻组织SNRNA-SEQ方案结合了我们认为对质量很重要的两个主要特征，因为它们可以防止核“泄漏”有价值的信号，并同时用外源RNA信号污染了条形码的核混合物。第一个特征是离心水平非常低，我们发现这既会导致信号丢失和外源信号增加。第二个特征是根据McCarroll Laboratory Protocol46，包括赋形剂试剂Basf Kollidon VA-64。 　　以下方法适用于Salk Institute和San Diego分校生成的SNMC-Seq2和SNATAC-SEQ数据集。 　　根据美国国立卫生研究院指南，在第18-00006号协议下，由美国国立卫生研究院指南进行了涉及动物的程序，并由机构动物护理和使用委员会批准。雄性C57BL/6J小鼠在8周龄时从杰克逊实验室购买，并在12小时的黑暗/光周期中保持在12小时的黑暗/光周期，温度（20–22°C）和湿度（30-70％），而食物Ad Adbitum则在解剖前1周。 　　将大脑从56至63天大的小鼠中提取，并在冰冷的解剖介质中立即切成0.6毫米的冠状切片。根据艾伦脑的参考地图集，将拖把从沿前轴沿前轴的切片进行了解剖（扩展数据图1D）。在解剖过程中将切片保存在冰冷的解剖培养基中，并立即冷冻在干冰中，以进行随后的聚集和核产生。对于细胞核分离，为每种生物学重复汇集了从15至23只小鼠的摩托车的解剖区域，并为每个区域处理两个重复。如先前的研究9,10所述，通过流式细胞仪分离核。简而言之，如先前所述的9，通过在蔗糖缓冲液中的均质化产生核，并将产生的核颗粒分为两个等分试样。一个等分试样进行了蔗糖梯度纯化和Neun标记（SNMC-SEQ2），第二个等分试样直接进行标记（SNATAC-SEQ）。 　　先前针对SNMC-Seq2（参考文献14）描述了Bisulfite转换和库准备的详细方法，该协议可在协议上获得。IO47。使用Illumina Novaseq 6000仪器（RRID：SCR_016387）对SNMC-SEQ2库进行测序，该仪表板和150 bp配对端模式。 　　如前所述28,48进行组合条形码SNATAC-SEQ。将分离的脑核用摆动的铲斗离心机（在4°C下为500g持续5分钟； 5920r，eppendorf）。将核沉淀重悬于1 mL核核透透缓冲液中（5％BSA，0.2％igepal-Ca630，1 mM二硫代硫代醇并完整，无EDTA的无EDTA蛋白酶抑制剂鸡尾酒鸡尾酒鸡尾酒鸡尾酒（Roche），在PBS中进行500G，在4°C; 5920 rid; 5920 rid; 5920 rid;SCR_018092）。将核重悬于500μl高盐标记缓冲液中（36.3 mm三乙酸乙酸盐（pH 7.8），72.6 mm钾，乙酸钾，11 mm mg-乙酸盐和17.6％DMF），并使用血细胞计器计数。将浓度调节至每9μl的4,500个核，并将4,500个核分配到96孔板的每个孔中。为了进行标记，使用Benchsmart 96（Mettler Toledo; RRID：SCR_018093）添加1μL条形码TN5转座体48，混合了五次，并在37°C下用摇动在37°C下孵育60分钟（500 r.p.p.m.）。为了抑制TN5反应，将10μl的40 mM EDTA与台式96（Mettler Toledo）一起添加到每个井中，并在37°C下在37°C下孵育15分钟，并摇动（500 R.P.M.）。接下来，使用基准96（Mettler Toledo）添加20μL2×排序缓冲液（PBS中2％BSA和2 mM EDTA）。将所有井聚合到FACS管中，并用3μMDRAQ7（细胞信号传导）染色。使用SH800（Sony），将40个核分类为八个96孔板（总共768孔），其中含有10.5μlEB（25 pmol Primer i7，25 pmol引物i5和200 ng bsa（Sigma））。排序板的准备和所有下游移液步骤均在Biomek i7自动化工作站（Beckman Coulter; RRID：SCR_018094）上执行。添加1μl0.2％SDS后，将样品在摇动（500 R.P.M.）下在55°C下孵育7分钟。将Triton-X（12.5％;1μL）添加到每个孔中以淬灭SDS。接下来，添加12.5μlNEBNEXT高保真2×PCR主混合（NEB），并将样品放大PCR（72°C持续5分钟5分钟，98°C 30 s （98°C 10 s，63°C 30 s和72°C或60 s）×12个循环，保持在12°C）。PCR之后，将所有井合并。根据Minelute PCR纯化套件手册（Qiagen）使用真空歧管（Qiavac 24 Plus，Qiagen）纯化库，并使用Spri珠（0.55×和1.5×; Beckmann Coulter）进行尺寸选择。图书馆再净化一次，用尖峰珠（1.5倍，贝克曼·库尔特）净化。使用量子荧光表（Life Technologies; RRID：SCR_018095）对文库进行定量，并使用挂接（高灵敏度D1000，Agilent）验证核小体模式。该库是使用自定义测序引物，25％Spike-In库和以下读取长度进行测序的HISEQ2500 Sequencer（Illumina; rrid：scr_016383）：50 + 43 + 37 + 50（read1 + index1 + index1 + index2 + read2）11。 　　对于Smart-Seq处理，我们执行了具有正面和负面对照的程序，如前所述3。用于测序的Smart-Seq V4超低输入RNA试剂盒（634894，Takara）用于反向转录poly（a）RNA并扩增全长cDNA。将样品在八孔条中放大18个周期，一次以12-24条的速度进行。所有样品均通过Nextera XT DNA文库制备（FC-131-1096，Illumina）进行了使用Nextera XT索引KIT V2（FC-131-2001，Illumina）和一个自定义索引集（Integrated DNA Technologies）。Nextera XT DNA库制剂是根据制造商的说明进行的，并进行了修改，以将所有试剂的体积减少到原始协议的0.4×或0.5倍。 　　对于10x V2处理，我们使用了铬单细胞3'试剂盒V2（120237，10x基因组学）。我们遵循制造商的说明，以进行细胞捕获，条形码，逆转录，cDNA扩增和图书馆构造。我们针对每个单元格的测序深度为60,000个读取。 　　对于10x V3处理，我们使用了铬单细胞3'试剂盒V3（1000075，10x基因组学）。我们遵循制造商的说明，以进行细胞捕获，条形码，逆转录，cDNA扩增和图书馆构造。我们针对每个单元格的测序深度为120,000个读取。 　　如前所述3进行了智能seq V4库的处理。简而言之，在Illumina Hiseq2500平台上对库进行了测序（配对末端，读取长度为50 bp），并使用RefSeQ notation GFF gff FFF将2016年1月18日从NCBI取回的RefSeq notation grcm38.p3（mm10）对齐grcm38.p3（mm10）。（https://www.ncbi.nlm.nih.gov/genome/annotation_euk/all/）。使用Star v2.5.349进行序列比对。使用星级选项“ BAMRemaiduplicates”对PCR重复分解物进行掩盖和删除。仅使用独特对齐的读数进行基因定量。使用R基因组学软件包（RRID：SCR_018096）50计算基因计数，并在“ IntersectionNotempty”模式中分别在“ IntersectionNotempty”模式下函数分别计算出基因计数。对于Smart-Seq V4数据集，我们仅使用外显子区域进行基因定量。符合以下任何标准中的任何一个细胞已被删除： <100,000 total reads, <1,000 detected genes (CPM > 0), <75% of reads aligned to the genome or CG dinucleotide odds ratio > 0.5. Cells were classified into broad classes of excitatory, inhibitory and non-neuronal based on known markers, and cells with ambiguous identities were removed as doublets3. 　　10x v2 and 10x v3 libraries were sequenced on Illumina NovaSeq 6000 (RRID: SCR_016387), and sequencing reads were aligned to the mouse pre-mRNA reference transcriptome (mm10) using the 10x Genomics CellRanger pipeline (version 3.0.0; RRID: SCR_017344) with default parameters. Cells were classified into broad classes of excitatory, inhibitory and non-neuronal based on known markers. Low-quality cells that fit the following criteria were filtered from clustering analysis. Different filtering criteria were used for neurons and non-neuronal cells as neurons are bigger than non-neuronal cells and contain more transcripts. For scRNA datasets, we excluded neurons with fewer than 2,000 detected genes and non-neuronal cells with fewer than 1,000 detected genes; for snRNA datasets, we excluded neurons with fewer than 1,000 detected genes and non-neuronal cells with fewer than 500 detected genes. Doublets were identified using a modified version of the DoubletFinder algorithm51 and removed when the doublet score was greater than 0.3. 　　Paired-end sequencing reads were demultiplexed and aligned to the mm10 reference genome using bwa52. After alignment, we converted paired-end reads into fragments and for each fragment, we checked the following attributes: (1) mapping quality score MAPQ; (2) whether two ends are appropriately paired according to the alignment flag information; and (3) fragment length. We only keep the properly paired fragments whose MAPQ (–min-mapq) is greater than 30 with fragment length less than 1,000 bp (–max-flen). Because the reads have been sorted based on the names, fragments belonging to the same cell (or barcode) are naturally grouped together, which allows for removing PCR duplicates. After alignment and filtration, we used Snaptools (https://github.com/r3fang/SnapTools; RRID: SCR_018097) to generate a snap-format file that contains metadata, cell-by-bin count matrices of various resolutions and cell-by-peak count matrices. 　　The method for calculating enrichment at the transcription start site (TSS) was adapted from a previously described method53. TSS positions were obtained from the GENCODE database (RRID: SCR_014966). Briefly, Tn5-corrected insertions were aggregated ±2,000 bp relative (TSS strand-corrected) to each unique TSS genome-wide. Then, this profile was normalized to the mean accessibility ±1,900–2,000 bp from the TSS and smoothed every 11 bp. The maximum of the smoothed profile was taken as the TSS enrichment. We excluded any single cells that had fewer than 1,000 unique fragments or a TSS enrichment of less than 10 for any sample sets. 　　After filtering out low-quality nuclei, we used Scrublet (RRID: SCR_018098)54 to remove potential doublets for every sample set. Cell-by-peak count matrices were used as input, with default parameters. 　　We implemented a versatile mapping pipeline (cemba-data.rtfd.io) for all the single-cell methylome-based technologies developed by our group9,14,30. The main steps of this pipeline included: (1) demultiplexing FASTQ files into single-cell files; (2) reads-level quality control; (3) mapping; (4) BAM file processing and quality control; and (5) final molecular profile generation. The details of the five steps for snmC-seq2 have been previously described14. We mapped all the reads onto the mouse mm10 genome. After mapping, we calculated the methyl-cytosine counts and the total cytosine counts in two sets of genome regions for each cell: the non-overlapping 100-kb bins tiling the mm10 genome, which was used for methylation-based clustering analysis, and gene body regions ± 2 kb, which was used for cluster annotation and cross-modality integration. 　　We filtered the cells based on five quality metrics: (1) the rate of bisulfite non-conversion as estimated by the rate of methylation at CCC positions (mCCC) < 0.03 (the mCCC rate reliably estimates the upper bound of the bisulfite non-conversion rate9); (2) the overall mCG rate > 0.5; (3) the overall mCH rate < 0.2; (4) the total final reads (combining R1 and R2) > 500,000; and (5) the total mapping rate (using Bismark55) > 0.5。 　　先前描述了SNMC-SEQ2数据的聚类步骤30。简而言之，我们根据非重叠100 kb bins矩阵的β-二比分布计算后MCH和MCG速率。然后，我们选择了前3,000个高度可变特征来执行主组件分析（PCA），并分别找到MCH和MCG的主体组件。我们将两种甲基化类型的主要成分加在一起，以构建K-Neart邻居（KNN）图，并反复运行Leiden社区检测算法56以获得共识聚类结果。聚类的停止标准考虑了标记基因的数量，基于群集分配的可再现监督模型的准确性和最小簇大小。我们在两次迭代中进行了聚类，以获取主要类型和细粒类型，以与其他模式进行进一步整合。 　　有关库大小的估计，请参见图1E的扩展数据。对于每个数据集，我们估计了基于重复序列读取速率的唯一RNA或DNA片段（F）的总库大小。唯一映射读取的数量是，其中s是测序读取的总数。使用此方程，我们使用s，nunique的中值对f进行数值求解。 　　有关转录组分析，请参见图1。使用R套件SCRATTCH.HICAT3（RRID：SCR_018099；可在https：//github.com/alleninstitute/scrattutter/scrattch.hicat）中独立执行每个SCRNA-SEQ和SNRNA-SEQ和SNRNA-SEQ数据集的聚类。该软件包通过依次进行更精细的分裂来支持迭代聚类，同时确保所有成对的簇，即使在最佳级别，也可以通过严格的差分基因表达标准分离。对于SCRNA 10X数据集，我们使用Q1。th= 0.4，q.diff.th = 0.7，de.score.th = 150 and.cell.cells = 10。对于SnRNA 10x数据集，我们使用Q1.Th = 0.3，q.diff.th = 0.7，de.score.score.score.score.score.score.score.th = 100和min.cells = 10 sm s scrsn = 10。q.diff.th = 0.7, de.score.th = 150 and min.cells = 4. For the snRNA SMART dataset, we used q1.th = 0.4, q.diff.th = 0.7, de.score.th = 100 and min.cells = 4. We further performed consensus clustering by repeating iterative clustering on a subsample of 80% of cells, resampled 100 times, followed by final clustering based在共聚类概率矩阵上。使用此过程，我们可以微调群集边界并评估群集不确定性。 　　接下来，我们删除了低质量和双线驱动的簇。我们在每个子类中的每对簇之间进行了差异基因表达分析。如果任何集群的基因上调≤2（倍数变化> 2，fdr） < 0.01, with additional dataset-specific parameters listed in the previous paragraph) than another cluster, and had a substantially lower average number of detected genes per cell, we flagged the cluster as low quality and removed it from further analysis. Next, if the upregulated genes between any two clusters within a subclass were predominantly marker genes for a different subclass, and one of the clusters had a significantly higher average of genes detected per cell and UMI count, we flagged the cluster as a potential doublet cluster and removed it from further analysis. These criteria led to the exclusion of 8.3% of all cells, the vast majority of which came from the two 10x v3 datasets (scRNA 10x v3 A and snRNA 10x v3 B). While the 10x v3 platform boosts the gene detection for good cells, it does the same to damaged cells or debris, leading to an increased number of clusters that were excluded for these datasets. 　　To provide a consensus cell-type taxonomy across all transcriptomic datasets, we developed an integrative clustering analysis across multiple data modalities. This procedure is available via the harmonize function of the scrattch.hicat package. Unlike Seurat/CCA57, which aims to find aligned common reduced dimensions across multiple datasets, this method directly builds a common adjacency graph using the cells from all datasets, and then applies the Louvain community detection algorithm58. We extended the cluster merging algorithm in the scrattch.hicat package to ensure that all clusters can be separated by conserved differentially expressed genes across platforms. The i_harmonize function, similar to the iter_clust function in the single-dataset clustering pipeline, applies integrative clustering across datasets iteratively while ensuring that all the clusters at each iteration are separable by conserved differentially expressed genes. 　　To build a common adjacency matrix incorporating samples from all the datasets, we first chose a subset of datasets that we used as ‘reference datasets’. For this study, we used the 10x v2 single-cell dataset from the Allen Institute (scRNA 10x v2 A) and the 10x v3 single-nucleus dataset from the Broad Institute (snRNA 10x v3 B) as the reference datasets, as both are large datasets that provide comprehensive cell-type coverage and relatively sensitive gene detection. 　　The key steps of the pipeline are outlined: (1) perform single-dataset clustering (Methods described above). (2) Select the anchor cells for each reference dataset. For each reference dataset (scRNA 10x v2 A or snRNA 10x v3 B), we randomly sampled up to anchor cells per cluster to normalize coverage for each cell type. This is the only step that uses the dataset-specific clustering information. (3) Select highly variable genes. Highly variable gene selection and dimensionality reduction by PCA were performed using the scrattch.hicat package. We removed principal components with a Pearson correlation coefficient of more than 0.7 with log2(Ngenes). This step was implemented to mitigate the effect of cell or nucleus quality on gene expression variability, and to select only biologically relevant principal components. For each remaining principal component, Z-scores were calculated for gene loadings. The top 100 genes with an absolute Z-score greater than 2 were selected as highly variable genes. The highly variable genes from each reference dataset were combined. (4) Compute KNNs. For each cell in each query dataset, we computed its KNNs (k = 15) among anchor cells in each reference dataset (scRNA 10x v2 A or snRNA 10x v3 B), based on the highly variable genes selected above. The RANN package was used to compute KNN based on the Euclidean distance when the query and reference dataset was the same. To compute nearest neighbours across datasets, we used correlation as a similarity metric. (5) Compute the Jaccard similarity. For every pair of cells from all datasets, we computed their Jaccard similarity, defined as the ratio of the number of shared KNNs (among all anchors cells from all the reference datasets) divided by the number of combined KNNs. (6) Perform Louvain clustering. (7) Merge clusters. To ensure that every pair of clusters are separable by conserved differentially expressed genes across all datasets, for each cluster, we first identified the top three most similar clusters. For each pair of such closely related clusters, we computed the differentially expressed genes in each dataset. We focus on the conserved differentially expressed genes that are significant in at least one dataset, while also having more than twofold change in the same direction in all but one datasets. We then computed the overall statistical significance based on such conserved differentially expressed genes for each dataset independently. If any of the datasets passed our differentially expressed gene criteria described in the ‘clustering’ section, the pair of clusters remained separated; otherwise they were merged. Differentially expressed genes were recomputed for the merged clusters, and the process was repeated until all clusters were separable by the conserved differentially expressed genes criteria. If one cluster had fewer than the minimal number of cells in a dataset (4 cells for SMART-Seq and 10 cells for 10x), then this dataset was not used for differentially expressed gene computation for all pairs involving the given cluster. This step allows detection of unique clusters absent in some platforms. (8) Iterative clustering. Repeat steps 1–6 for cells within each cluster to gain finer-resolution clusters until no more clusters can be found. (9) Final compilation and merging of clusters. Concatenate all the clusters from all of the iterative clustering steps and perform the final merging as described in step 6. 　　For each pair of clusters, we computed the conserved differentially expressed genes, that is, those which are significantly differentially expressed in at least one dataset, with a twofold or more change in expression in the same direction among 70% of datasets. To allow computation of differentially expressed genes involving cell types only present in a subset of datasets, only the datasets with enough cells (based on min.cells parameter) for both cell types under comparison were used. We selected the top 50 genes in each direction. After pooling genes from all pairwise comparisons, we identified a total of 3,792 marker genes (Supplementary Table 6). 　　To facilitate direct comparison, we projected gene expression of all datasets to the space of a given reference dataset. To do that, we leveraged the KNN matrices computed during the iterative joint clustering step to adjust the expression values for systematic differences between datasets. During each iteration of the joint clustering, for cells in each dataset, we used the average gene expression of their KNNs among the anchor cells from the reference dataset as the adjusted expression in the reference space. At the top-level clustering, we imputed the expression for all genes. For each subsequent iteration, we only imputed the expression of the high-variance genes and the conserved differentially expressed genes for the clusters defined in that iteration. We used this iterative approach for imputation because the nearest neighbours based on the genes chosen at the top level may not reflect the distinction between the finer types, and the imputed values for the differentially expressed genes that define the finer types consequently are not accurate based on these nearest neighbours. Therefore, we deferred imputation of the differentially expressed genes between the finer types to the iteration when these types were defined. This method is provided in the impute_knn_global function in the scrattch.hicat package3. We imputed the gene expression matrix for both reference datasets used in the integrative clustering. 　　We first computed the average adjusted expression of marker genes for each cluster. This average was computed using each of the two reference datasets (scRNA 10x v2 A and snRNA 10x v3 B). Then, the two matrices were concatenated. We constructed a hierarchy (tree) using the build_dend_harmonize function in the scrattch.hicat package3. 　　We performed PCA based on imputed gene expression matrices of 3,792 marker genes using the 10x single-nucleus dataset from the Broad Institute as the reference, and selected the top 50 principal components (93% variance explained). We removed principal components with Pearson correlation coefficient > 0.6与log2（Ngenes）减少与检测基因数量相关的偏差。UMAP用于将细胞嵌入两个维度，其中nn.nn.nn.ne.nn.nn.nn.nn.nn.nn.nn.nn.nn.nn.nn.nn.nn。 　　有关Metaneighbor分析，请参见图1G。为了量化在七个转录组数据集中簇的可复制性，我们应用了无监督的metaneighbor（RRID：SCR_016727）的修改版本24。Metaneighbor使用邻居投票算法和跨数据库验证方案来量化多个数据集的群集相似性。它需要一组非均衡数据集，一组群集标签和一组高度可变的基因。我们使用了所有通过七个单细胞转录组数据集的质量控制标准以及通过独立聚类获得的标签（补充表5）的所有单元格（补充表5）。我们使用Metaneighbor中的变量程序选择了310个高度可变基因，这些基因被检测到所有数据集中的高度可变。 　　我们在两步过程中定义了可复制的簇：首先，我们量化了跨数据集的簇之间的相似性，然后我们提取了高度相似簇的组或“元群集”。我们使用Metaneighborus函数获得簇之间的初始相似性矩阵。默认情况下，群集相似性被量化为一个单VS-ALL AUROC：给定训练群集（在一个数据集中），我们询问与测试数据集中的所有其他单元格相比，测试群集（在另一个数据集中）的相似单元格是如何用于训练单元的。为了使群集变得更加严格，我们将单VS-ALL AUROC矩阵转换为一个最佳的AUROC矩阵：而不是从测试数据集中对所有单元格中的测试单元进行排名，而是从最佳匹配群集中将它们与单元格进行了比较。这种修改确保只有最佳匹配的AUROC> 0.5，从而有助于识别倒数最佳命中。为了解释性和计算效率，我们采用了以下惯例：通过将其与第二个最佳匹配集群进行比较，获得了群集的最佳AUROC，该集群的第二好的AUROC是通过计算1 - AUROC的最佳匹配群集群的AUROC，以及所有其他簇的AUROC，都可以获得0的AUROC，因为我们只能找到最佳的匹配项。为了提取元群集，我们将单VS最好的AUROC解释为图形，如果节点为簇，并且边缘连接节点，则如果它们是互惠的最佳命中。我们将元群集定义为此图中的连接组件。我们可以通过要求最佳命中率超过AUROC阈值来获得更健壮的元群集。在实践中，我们注意到，一VS最好的AUROC> 0.7在元群数和可重复性强度之间提供了良好的平衡。 　　为了可伸缩，我们通过以下方式修改了metaneighbor。在Metaneighborus函数中，我们删除了单元格相似性网络的等级标准化（通过将参数fast_version设置为true）和邻居投票的节点度归一化，从而可以对邻居投票程序进行分析简化。将变量程序应用于超过该大小的数据集的50,000个单元的随机子集。 　　Metaneighbor分析进一步使我们能够检查计算聚类程序的一致性（图1H）。我们运行了三个广泛使用的单细胞分析软件包25,26,27，以生成每个数据集的细粒聚类。这些集群分析未进行优化或手动策划；取而代之的是，我们使用“现成”计算程序来测试相对直接和自动化的分析结果的鲁棒性。因此，预计这些簇在生物学上的意义和鲁棒程度不如更自定义的程序，例如我们的参考聚类融合了差异表达的分析以验证细胞类型的生物学现实。使用三个现成的群集分析，我们通过迭代合并成对的群集对，从而创建了一系列日益粗糙的聚类序列，以最大程度地提高计算方法（ARI-Merering）的一致性。最后，在每个分辨率的每个级别上，我们都使用Metaneighbor来计算数据集中高度可复制（AUROC> 0.7）的簇数。该分析的结果表明，超过30-50个簇的数据的精细分区的可复制性有限。 　　有关SNMC-SEQ2的聚类分析，请参见图6a，b。我们将两种甲基化类型（CG和CH）的主要成分串联在一起，并使用它们构造了KNN图，其次是莱顿社区检测56。我们多次重复集群分析以获得共识聚类结果。聚类的停止标准考虑了标记基因的数量，基于群集分配的可再现监督模型的准确性和最小簇大小。我们在两个迭代中进行了聚类，以获得主要类型和细粒细胞类型，以与其他模态进行比较。 　　使用T-Distrib的随机邻居嵌入59（T-SNE; PERPLEXITY = 30）的二维嵌入，使用Scanpy Package60的实现基于顶部主组件进行计算。 　　有关SNATAC-SEQ的聚类分析，请参见扩展数据图6C，d。我们使用Snapatac Pipeline48来鉴定5-KB分辨率中的二核细胞序列矩阵作为输入。通过检查沿标记基因体的染色质可及性，将细胞簇注释到细胞类型。然后，在内侧神经节隆起（MGE）衍生和尾神经节杰出（CGE）衍生的抑制性GABA能中间神经元上进行了另一轮聚类，以识别子细胞类型。 　　有关多模式集成，请参见图2。 　　我们使用Liger（RRID：SCR_018100）来整合单细胞转录组和表观基因组数据，如前所述29，并进行了一种修改。我们在Liger软件包中使用优化功能在所有数据集上执行联合分解，除了甲基化（七个RNA数据集和一个ATAC数据集）以推断共享（W）和数据集特异性（VI）Metagene因子和细胞因子负载（HI）。然后，我们使用所得的W使用Liger包装中的索词函数来计算甲基化数据的细胞因子负载（HI）。我们发现，这种策略比共同分解所有八个数据集的集成更好，这可能是因为甲基化和所有其他数据集之间的数据集的逆关系和大小不平衡使人学习共享元语的学习变得复杂。我们的分析仅使用每个数据生成组注释的细胞作为传递质量控制。我们没有执行任何数据插补或平滑，而是使用Liger套件中的归一化和缩小的函数将每个数据集的原始单元格计数矩阵进行标准化和缩放。接下来，我们使用默认设置的QuantileAtignsnf函数来执行所有八个数据集的细胞因子矩阵（HI）的分数归一化。最后，我们在标准化的细胞因子矩阵（HI）上进行了卢旺聚类，以获得关节簇。我们进行了两轮整合和关节聚类。在第一轮中，我们分别整合了跨数据集和所有神经胶质的所有神经元。然后，我们针对四个神经元子类别进行了第二轮整合和聚类：兴奋性IT神经元，兴奋性非IT神经元，MGE中间神经元和CGE中间神经元。我们使用k = 40个因子进行非神经元分析，首轮神经元分析的k = 30，所有第二轮分析的k = 20。 　　SingleCellFusion30旨在鲁棒地整合DNA甲基化，ATAC-SEQ和/或RNA-SEQ数据。我们迭代地应用单次融合以整合来自八个数据集（补充表1）并共同调用细胞簇的所有神经元。为了整合宽阔的细胞类型，我们进行了三轮整合。对于上一轮生成的每个细胞簇，仅通过仅在该群集中的单元格上重新填充SCF，将其进一步分为较小的簇。在第一轮中，我们从8个数据集中的所有神经元上运行SCF，并获得了10个广泛的神经元簇。第二轮和第三轮分别产生29个簇和56个细粒簇（补充表3）。 　　该过程包括四个主要步骤：预处理，模式内平滑，跨模式插补以及聚类和可视化。（1）对于预处理步骤，我们为每个数据集定义了一个基因特征矩阵。基于液滴的RNA-seq特征（10倍）是log10（CPM + 1）归一化的；全长RNA-SEQ（SMART-SEQ）功能为Log10（TPM + 1）归一化。SNATAC-SEQ数据由基因体内的读数表示，通过Log10（RPM + 1）归一化，其中CPM代表映射的每百万个读取（计数标准化），TPM代表TPM代表每百万个读取映射的成绩单（长度归一化）和RPM读取的读数分别为每百万个读取（分别读取范围），分别绘制了范围（均值范围）。DNA甲基化数据由平均基因体MCH水平表示，每个细胞的全局（全基因组）平均MCH水平归一化。对于每个数据集，我们仅使用高质量的细胞（通过质量控制）和高度可变的基因（n = 4,000–6,300）进行进一步分析。为了选择高度可变的基因，对于RNA-seq和atac-seq数据集，我们首先删除了在不到1％的细胞中表达的基因。然后，根据细胞（CPM）的平均表达，我们将其余基因分为10个垃圾箱。对于每个垃圾箱，除表达最多的垃圾箱外，我们选择了表达分散最多的基因的前30％（方差/平均值）作为高度可变的基因。对于DNA甲基化数据集，我们首先选择了在超过95％的细胞中具有20多个胞质覆盖率的基因，然后根据其平均标准化MCH水平将其余基因分为10个垃圾箱 - 每个细胞的全球MCH标准化的原始MCH水平。对于每个垃圾箱，我们选择了最大的基因的最高基因作为高度可变的基因。（2）对于模式内平滑步骤，以减少特征矩阵的稀疏性和噪音， 我们使用数据扩散共享了具有相似曲线的单元格之间的信息。该过程改编自参考。61，并在参考文献中详细描述。30.在这里，我们完全遵循参考。30使用[NDIM = 50，K = 30，Ka = 5]，RNA-Seq数据集的[P = 0.7]，DNA甲基化数据集的[P = 0.9]，ATAC-SEQ数据集的[P = 0.9]和[P = 0.1]。（3）对于限制K-Partners（RKP）步骤的跨模式归类，为了整合所有八个数据集，我们将SCRNA 10x V2归为其他所有七个数据集中的单元格的基因特征。插补是在SCRNA 10x V2 A数据集和每个其他数据集之间成对完成的。对于每个成对的插补，我们遵循参考文献中描述的过程。30具有20个RKP和松弛参数3 [K = 20，Z = 3]。我们没有在低维空间中使用欧几里得距离，而是在两个数据集中都使用了（翻转的）Spearman相关系数，这是两个不同模态的细胞之间的距离度量。（4）对于聚类和可视化步骤，我们从单元格矩阵开始，其中包括八个数据集中的所有单元格，并且特征是SCRNA 10X V2 A数据集的高度可变基因。我们将特征的维度降低到了前50个主要组件中。接下来，我们在主组件矩阵上执行了UMAP嵌入40（n_neighbours = 60，min_dist = 0.5）。最后，我们在最终主成分矩阵（Euclidean距离，k = 30，分辨率= 0.1）上产生的KNN图（对称性，未加权）上进行了莱顿聚类。 　　对于扩展数据图7E，我们使用估算的SCRNA 10X V2 A基因特征（Log10（CPM + 1））为从单链灌注积分生成的八个不同数据集创建了簇质心的嵌入。群集由单个数据集群集以及与单链灌注的关节聚类来定义。通过平均估算的SCRNA 10x V2 A基因谱跨细胞计算簇质心。在获得基因群矩阵后，我们应用了PCA将其降低到50个特征尺寸，然后将umap嵌入给min_dist = 0.7和n_neighbours = 10。 　　对于图2E，为了比较跨数据方式的分子信号，将所有信号归一化为[0，1]。这是通过首先通过数据集特异性归一化获得分子信号来实现的（步骤1），然后进行线性转换（步骤2）。在步骤1中，对于Smart-Seq数据集，我们显示Log10（TPM + 1）;对于10倍RNA-seq数据集，我们显示Log10（CPM + 1）;对于ATAC-SEQ数据集，我们显示Log10（RPM + 1）标准化基因体数；对于DNA甲基化，我们显示了通过每个细胞的全局MCH水平标准化的基因MCH。对于步骤2，我们应用了线性转换以将信号的范围映射到[0，1]。对于DNA甲基化以外的数据集，我们应用了以下公式： 　　细胞的数据集特异性基因级信号在哪里，并且分别定义为所有细胞中最底层的2个百分位数和前两个百分位数。对于DNA甲基化数据集，我们应用了以下公式： 　　信号仍然映射到[0，1]但翻转 - 图上的高信号意味着低水平的DNA甲基化。我们这样做是为了使DNA甲基化信号与基因表达和开放染色质信号对齐，因为DNA甲基化是基因表达的抑制标志物，并且与之负相关。在这些公式中，分别定义为所有细胞中最底层的2个百分位数和前50个百分位数。 　　对于图2D，对于每个基因，细胞级信号的归一化方式与图2E步骤1中所述的方式相同。群集级信号是簇中细胞跨细胞的平均细胞级信号。在以这种方式获得基因矩阵后，对于非DNA甲基化数据集，通过每个簇的最大值（列）进一步将矩阵进一步标准化；对于DNA甲基化数据集，未进行进一步的归一化，因为它们已经通过细胞标准化了。 　　对于扩展数据，图7G，H，使用基因表达水平（log10（CPM + 1）（CPM + 1）；通过log10（cpm + 1）；通过scrna 10x v2 v2 v2 v2 a dataaset a date a Dataaset As emapity a scrna singecellfusion测量或估算的基因表达水平（CPM + 1）；使用基因表达水平（CPM + 1）；测量或估算的每个细胞表达水平，显示了每种细胞类型（单细胞输液）和每个数据集之间的细胞质心之间成对的Spearman相关系数。将质心轮廓计算为来自同一细胞类型和相同数据集的细胞跨细胞的平均值。行和色谱柱排序是相同的，它是通过平均链接和欧几里得距离的上述质心特征上的层次聚类生成的。扩展数据图7G显示了宽水关节聚类（10个子类；单链灌注L0）之间的相关性（补充表8）；扩展数据图7H显示了四个示例Broad级子类（MGE，CGE，L2/3 IT和L4/5 IT）的高级关节聚类（总共56个簇；并非全部显示； SingleCellFusion L2）（补充表8）（补充表8）。 　　对于扩展数据中的协议度量图7C，我们计算了Liger Paper29中所述的数据集协议指标。简而言之，我们使用非阴性矩阵分解（NMF；对于Liger）或PCA（用于单链灌注）进行了降低，并为每个数据集构建了KNN图。然后，我们使用Liger或SinglecellFusion的联合潜在空间构建了KNN图，并计算了关节空间中最近的邻居最近的邻居仍然是最近的邻居。该度量标准评估了联合潜在空间保持每个数据集的结构的能力。一项接近0的协议指标表明，从单个数据集中对结构的保存差，而一项接近1的协议则可以保留结构。 　　对于扩展数据中的对齐度量图7D，我们计算了Liger29和Seurat57论文中所述的数据集对齐指标，只是我们首先降低了样采样的单元格，以使群集比例和总细胞总数在所有数据集中都是相同的。接下来，我们使用Liger或SinglecellFusion的联合潜在空间构建了KNN图，并计算了每个数据集周围最近的邻居的哪一部分。然后，我们将指标归一化为0（无比对）和1（数据集的完美混合）。该度量标准评估了联合潜在空间对数据集的对齐程度。请注意，最大化一致性和最大化协议是竞争目标。例如，根据球形高斯分布将所有数据集的细胞随机混合通过随机混合细胞，可以轻率地最大化对齐。相反，通过简单地将非重叠潜在表示形式分配给所有数据集，就可以使人们最大程度地提高协议。但是，方法必须平衡这些竞争目标，以在对齐和协议指标上得分很高。 　　对于扩展数据图7F，要获得群集级基因信号，我们首先以与图2E的步骤1相同的方式获得了标准化的细胞级信号，然后以簇中的细胞跨细胞的平均细胞级信号进行操作。 　　根据图2中的细胞–细胞整合，要具有每种细胞类型的足够全基因组覆盖率，我们将共群体进一步合并为较高的水平，以增加每个群集的覆盖率，我们称其为表观基因组群集水平。 　　对于在SNMC-SEQ2数据中调用的DMR，我们将单细胞Allc文件合并到每个群集的伪式式级别中，然后使用Methylpy62 DMRFind函数来计算所有群集中的MCG DMR。在分析之前，将添加每个配对CPG站点的基本调用。简而言之，甲基纤维函数使用置换拟合优点的置换基均方根测试，在所有样品中同时鉴定差异甲基化位点，然后将250 bp的差异甲基化位点合并为DMR。然后通过检查预期计数观察到的计数的残基，将HYPO-DMR和超DMR分配给每个样品。我们还要求群集之间的MCG速率的最大差异大于0.3，从而过滤了DMR。 　　我们根据encode atac-seq管道（https://www.encodeproject.org/atac-seq/）调用峰值。对于每个细胞群，我们将所有正确配对的读数组合在一起，以生成一个伪膨化的ATAC-SEQ数据集，以进行单个生物学重复。此外，我们产生了两种伪复制物，每种伪造物都包括每个生物复制的一半读数。我们针对这四个数据集中的每个数据集以及两个生物学重复的数据池独立调用峰值。使用MACS263（RRID：SCR_013291）在TN5校正的单基插入上执行峰值调用：–Shift -75- shift-extsize 150 – Nomodel-call-call-call-summits – spmr-spmr – spmr – spmr – keps-spmr – keep-weep-weep-weep-dup all-q 0.01。我们将两侧的峰顶峰值延长至501 bp的最终宽度，用于合并和下游分析。为了生成可重复的峰列表，我们保持了峰值的峰值，即在汇总的数据集中检测到（1），并在两个单个生物学重复次数中重叠的50％或更多的峰值长度，或（2）在汇集的数据集中检测到（2），并且在两种峰值中均具有50％或更多的峰值长度。 　　To account for differences in performance of MACS2 based on read depth and/or the number of nuclei in individual clusters, we converted MACS2 peak scores (−log10(q value)) to score per million (SPM)64 and kept peaks with SPM > 2. We only kept reproducible peaks on chromosomes 1–19 and both sex chromosomes, and filtered ENCODE mm10 blacklist regions65（http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/mm10-mouse/mm10.blacklist.bed.gz）。最后，由于SNATAC-SEQ数据相对稀疏，因此我们仅选择了在每个群集中很大一部分细胞中被识别为开放染色质的元素。为此，我们定义了一组背景区域，通过随机从基因组中选择区域，同时从顺式调节元素的编码注册表（https://screen.encreen.encodeproject.orgg/）中排除了可访问的站点，从而匹配了每种单元类型的峰值区域。我们计算了每个细胞类型的核的比例，该细胞类型的ATAC片段映射到背景区域。接下来，我们拟合了零膨胀的beta模型，并在经验上确定了FDR <0.01的显着性阈值，以滤除潜在的假阳性峰。在下游分析中，至少一个集群中的FDR <0.01的峰区域包括在内。 　　我们使用“ BedTools相交”与“ -wa -u”参数来计算DMR和ATAC峰重叠66（RRID：SCR_006646）。 　　为了研究表观基因组数据中细胞数的调节元件识别效率，我们使用两种最丰富的细胞类型进行了饱和分析：L2/3 IT和L6 CT兴奋性神经元。分配给这两种单元格类型的总读取与散装seq相当。当我们有足够的细胞时，我们在每个群集中进行了三次，在每个群集中未替换了不同数量的细胞，并在可能的情况下分别使用细胞。在最后一组中，我们将每个单元格类型的所有单元格都用作最大参考。对于甲基组数据，我们在每个单元格组中称为L2/3 IT和L6 CT之间的DMR。每个细胞类型组都调用峰值。 　　我们使用Reptile35算法进行了增强剂预测。爬行动物是一种基于随机的监督方法，它结合了不同的表观基因组概况来源，并与基础水平DNA甲基化数据一起学习，然后区分增强子和基因组背景的表观基因组特征。我们使用CG甲基化，每个表观基因组簇和小鼠胚胎干细胞的染色质可及性，以与以前的研究相似的方式训练了该模型35,67。该模型首先是在小鼠胚胎干细胞数据上训练的，然后预测了我们称为每种细胞类型DMR的增强子得分的定量分数。这些阳性是在小鼠胚胎干细胞中前5,000个EP300峰的山顶上的2-KB区域。负面因素包括随机选择的5,000个启动子和30,000个2-KB基因组垃圾箱。垃圾箱与任何阳性或启动子没有重叠67。小鼠胚胎干细胞的甲基化和染色质可及性概况来自小鼠编码Project67。MCG Rate BigWig文件是使用软件AllCool（https://lhqing.github.io/allcools）从单元格 - 合并的Allc文件生成的。对于每种细胞类型的染色质可及性，我们合并了来自SNATAC-SEQ单元的所有片段，这些片段在集成分析中分配给该单元格类型，并使用“ DeepTools bamcoverage”来生成CPM均衡的Bigwig文件。所有Bigwig文件的垃圾箱尺寸均为50 bp。 　　我们使用了Jaspar 2020核心脊椎动物数据库68的724个基序位置权重矩阵（PWM），每个基序都能够分配相应的小鼠转录因子基因。对于每组爬行动物预测的增强剂，我们将区域长度标准化为中心±250 bp，并使用Meme Suite 69的FIMO工具来扫描每个增强器中的基序，而log grd率p值<10-6作为基序的阈值。为了计算基序富集，我们将成年非神经小鼠组织DMRS70用作背景区域。我们从背景设置的区域中减去了增强子，然后使用相同的方法在背景区域中扫描了图案。然后，我们使用Fisher的精确测试来找到富含该区域的主题和Benjamini -Hochberg程序来纠正多个测试。通过TFCLASS71分类对具有显着基序富集的转录因子进行分组。同一组中的基因共享非常相似的基序。 　　有关集群验证分析，请参见图4。 　　有关集群数量的下采样分析，请参见图4a – d。 　　预处理的方式与“与单链灌注的计算集成”部分中所述进行。预处理后，我们为每个数据集获得了一个基因特征矩阵。仅包括通过质量控制的神经元细胞（补充表1）和每个数据集的高度可变基因。 　　聚类（图4A）需要三个步骤。我们首先通过PCA缩小了特征尺寸（n = 50）。然后，我们使用主组件空间中的欧几里得距离之间在细胞之间构建了一个KNN图（k = 30）。我们最终使用固定分辨率参数（r = 6）应用了莱顿聚类算法。对于每个数据集，我们报告簇数，这是从完整数据集中随机降采样的单元格数的函数。错误条显示了S.E.M.（n = 10）的下采样。 　　该分析（图4C）旨在通过测试哪些聚类粒度最能保留细胞的基因级特征来估计数据集的最佳簇数。对于给定的数据集，一个基因矩阵，我们首先将基因特征随机分为两个集，分别用于聚类和验证。为了避免任何潜在的连锁，分裂是通过将染色体分为两组来完成的，从而使来自同一染色体的基因始终在同一组中。然后，我们仅使用不同的聚类分辨率进行了所有群集特征，在所有单元格上进行了所有单元格上的莱顿聚类（如图4A相关的方法所述）。聚类后​​，数据集中的每个单元格都会收到群集标签。接下来，我们将这些细胞随机分为训练和测试集。使用训练集的单元格，我们训练了一个监督模型，以预测基于群集分配的验证集基因特征。该模型是通过最小化模型预测和数据之间的MSE来训练的。这相当于预测细胞的基因特征为其簇质心。最后，我们通过计算测试集中的单元格的MSE来评估模型性能。这等同于估计测试集中单个单元格与使用训练集计算的群集质心之间的平方平方距离。作为簇数的函数（通过改变莱顿聚类中的分辨率参数），我们观察到MSE的U形曲线。曲线的最小点表示最合理的聚类分辨率。将此方案应用于每个数据集的每个数据集和不同的下采样级别，我们在图4C中报告了每个数据集的簇数与单元格数的函数。为了鲁棒，重复基因特征的随机分裂n = 5次；每次k = 5倍的交叉验证，重复n = 5次细胞的随机分裂。 　　扩展了图4C中使用的数据集聚类交叉验证方案，我们开发了一种交叉数据集交叉验证方法（图4D），通过将先前描述的dataset内部交叉验证方法与关节聚类方法相结合：singlecellfusion。首先，类似于数据集的交叉验证，我们将基因特征随机分为聚类和所有数据集的验证集。然后，我们通过在所有细胞和一半的基因特征（聚类特征集）上应用单链灌注来生成跨数据模式的集成聚类。聚类后​​，我们估算了上述每个数据集的验证功能集中的聚类MSE。将此方案应用于不同的细胞的不同下采样级别，我们在图4D中报告了簇数，这是每个数据集中的单元格数的函数。 　　有关综合分析，请参见图4E，f。我们收集了使用前面描述的四种集成聚类方法获得的群集（Conos，Liger，图1和单链球化的RNA共有聚类），以及从RNA数据集的独立群集中获得的“子类”水平。每个集成方法以两个粒度水平返回簇。我们命名了聚类L1的更粗级和聚类L2群集的较高水平。当我们希望研究神经元群集层次结构的一致性时，我们将分析重点放在转录组数据的神经元簇上。 　　为了量化可复制性，我们使用了Metaneighbor的相同修改版本，相同的数据集和与上面定义的相同的可变基因（请参见“ Metaneighbor Analysis”）。我们使用单VS最佳AUROC来获得群集相似性得分，然后计算出每个集成群集的平均AUROC分数（在存在群集的每对数据集上平均）。对于每种方法，我们都报告了整个集成簇的中位AUROC作为最终可重复性得分。为了量化聚类结果的总体相似性，我们计算了调整后的RAND索引。必要时，我们将调整后的兰德指数计算限制为标记的单元格的交点（每对方法都会重新计算该交集）。 　　为了评估不同平台支持不同细胞亚群的程度，我们评估了在有和没有平台内比较的情况下恢复相应细胞的能力的差异。使用conos37（图4G，H）进行细胞的聚类，使用Walktrap群落检测来识别分层细胞群体。分层簇的稳定性如下估计：进行了20个随机细胞亚采样弹，每个数据集中的每个细胞中的每个样品中的每个样品中的95％采样，并重复WallTrap层次层次群集程序。对于原始行走绑带树中的每个节点，我们通过找到最佳匹配的子树来评估稳定性相对于亚地采样的集合的最小特异性和灵敏度。为了评估仅基于跨平台比较恢复亚群的能力，我们删除了关节图（CONOS生成）中的平台内边缘（连接由同一平台生成的数据集）。这样，仅根据映射到另一个平台来检测亚群。修改的方法促进了在不同平台中常见的细胞群体的分组，因为它删除了关节图中的平台特定信息。 　　为了评估给定细胞类型的不同平台检测到的表达曲线的相似性（图4H），我们使用Jensen – Shannon差异来评估四个RNA-Seq平台之间基因表达模式的总体相似性（SCRNA 10X V3 A，SNRNA 10X V3 A，SNRNA 10X V3 A，SCRNA SMART SMART和SNRNA SMART）。具体而言，每个平台从每种单元格中对1,000个单元进行采样。如果细胞类型的细胞数量少于1,000个细胞，则进行替换的采样。省略了任何特定平台中不到1％（少于300个单元）的细胞类型。然后将每个基因检测到的分子在每个平台中的每个细胞类型中汇总在所有采样细胞中。通过每种细胞类型或平台的分子总数进行标准化，并计算了Jensen -Shannon差异。 　　有关研究设计的更多信息可在与本文有关的自然研究报告摘要中获得。