Basic Information 英文标题: An integrated transcriptomic cell atlas of human neural organoids 中文标题:人类神经类器官的整合转录组细胞图谱 发表日期:20 November 2024 文章类型:Article 所属期刊:Nature 文章作者:Zhisong He | Barbara Treutlein 文章链接:https://www.nature.com/articles/s41586-024-08172-8 Abstract Para_01 人类神经类器官是通过体外多能干细胞生成的,是研究人类大脑发育、进化和疾病的有效工具。 然而,目前尚不清楚现有方案覆盖了人类大脑的哪些部分,也难以定量评估类器官的变化性和保真度。 在这里,我们将36个单细胞转录组数据集整合到一个综合的人类神经类器官细胞图谱中,该图谱包含超过170万个细胞。 将这些数据映射到发育中的人类大脑参考图谱上,可以显示在体外生成的主要细胞类型和状态,并估计主要细胞与类器官对应物之间的转录组相似性跨不同方案。 我们提供了一个编程接口来浏览图谱和查询新数据集,并展示了图谱在注解类器官细胞类型和评估新的类器官方案方面的强大功能。 最后,我们证明了该图谱可以用作多样化的对照队列,用于注解和比较神经疾病的类器官模型,识别可能构成病理机制基础的基因和途径。 人类神经类器官细胞图谱将有助于评估类器官的保真度,表征受扰动和疾病状态,并促进方案开发。 Main Para_01 人类神经类器官是在体外生长的自组织三维人脑组织,正成为研究人类大脑发育、进化和疾病的机制的强大工具。 它们可以使用外部模式因子(例如形态发生素)来引导其向特定脑区的发展或驱动特定细胞类型的出现(指导协议)。 相反,非指导协议依赖于类器官的自我模式能力来产生多样化的细胞类型和状态。 Para_02 单细胞 RNA 测序(scRNA-seq)是一种强大的技术,用于表征复杂组织中的细胞类型异质性,并揭示了在神经类器官中可以发育的各种前体、神经元和胶质细胞类型的显著异质性,以及某些神经谱系的分化轨迹。 这些数据还使得人类神经类器官细胞与原代人脑细胞之间的比较成为可能,大多数分析显示分子特征具有强烈的相似性。 也报告了显著的差异,包括与培养基成分相关的差异基因表达和与糖酵解相关的代谢特征的扰动。 尽管如此,对类器官组织的分析支持早期大脑发育的有效再现,scRNA-seq 方法已被应用于研究神经细胞命运决定的分子基础、灵长类动物的进化差异以及神经障碍的病理变化。 然而,目前尚不清楚现有的方案能生成发育中的中枢神经系统哪些部分,哪些部分仍然缺失。 系统地量化神经类器官细胞与其原代对应物的转录组保真度也一直具有挑战性。 Para_03 在这项研究中,我们通过将涵盖多种人类神经类器官协议的36个单细胞RNA测序(scRNA-seq)数据集整合成一个综合转录组细胞图谱来应对这些挑战。 我们建立了一个分析管道,允许对类器官图谱与发育中的人脑参考图谱进行全面和定量的比较。 我们协调了初级系统和类器官系统中细胞群体的注释,估计了不同神经类器官协议生成不同脑区的能力和精确度,并确定了在神经类器官中代表性不足的主要细胞群体。 我们评估了神经类器官中神经元的转录组保真度,并确定了先前描述的细胞压力是区分体外神经元与原代神经元代谢状态的一个普遍因素,而不会强烈影响神经元细胞类型的核心身份。 我们将一个神经类器官形态发生筛选的数据映射到整合图谱上,以评估区域特异性和新状态的产生。 我们还收集了11个单细胞RNA测序数据集,模拟10种不同的神经疾病,并将整合数据映射到神经类器官图谱上进行细胞类型注释和差异表达(DE)分析。 最后,我们展示了该图谱可以通过将新数据投影到当前图谱上来扩展。 总的来说,我们的工作提供了一个丰富的资源和一个新的框架,用于评估神经类器官的保真度、表征受扰和疾病状态以及简化协议开发。 Data curation, harmonization and integration Para_01 为了构建一个转录组人类神经类器官细胞图谱(HNOCA),我们收集了来自36个数据集的单细胞RNA测序数据和详细、统一的技术和生物元数据,其中包括34个已发表的数据集和两个尚未发表的数据集(补充表1),经过一致的预处理和质量控制后,共包含177万个细胞(图1a)。 HNOCA代表了使用26种不同的神经类器官分化方案生成的细胞类型和状态,包括三种非指导性和23种指导性方案,在7到450天的时间点上进行(图1b)。 为了消除批次效应,我们实施了一个三步集成管道。 首先,我们将HNOCA投影到发育中的人类大脑的单细胞图谱上,使用参考相似度谱(RSS)。 然后,我们开发了snapseed(方法部分)来进行基于初步标记的分层细胞类型注释。 最后,我们根据snapseed注释使用scPoli进行标签感知的数据集成。 使用先前建立的基准测试管道评估不同集成方法的表现显示,scPoli在这些数据集上表现最佳(扩展数据图1)。 我们在scPoli表示的基础上进行了聚类,并根据经典标记基因表达、类器官样本年龄以及自动生成的细胞类型标签对聚类进行注释。 均匀流形近似和投影(UMAP)嵌入突出了三种神经元分化轨迹,分别对应于背侧端脑、腹侧端脑和非端脑群体,以及从祖细胞到星形胶质细胞和少突胶质细胞前体等胶质细胞类型的轨迹(图1c-e和扩展数据图2)。 来自非指导性和指导性协议的细胞分布在所有轨迹上(图1f)。 Fig. 1: Integrated HNOCA.
a, HNOCA构建流程概述。 b, 包含在HNOCA中的生物样本元数据。 c–f, 集成HNOCA的UMAP图,按第2级细胞类型注释(c)、选定标记基因表达谱(d)、样本年龄(e)和分化协议类型(f)着色。 g, 分配给HNOCA中不同细胞类型的细胞比例。每个堆叠条形图代表一个生物样本,按数据集分组,并按样本年龄递增顺序排列。顶部条形图显示36个数据集、类器官分化协议、协议类型。底部条形图显示样本年龄。 h, 集成HNOCA的UMAP图,按实时信息转换矩阵之间的顶级扩散成分(DC1)着色。流箭头可视化了细胞状态向更成熟细胞的推断流动。 i, 沿皮质伪时间的标记基因表达谱。 j, 非端脑神经元的UMAP图,按聚类着色和标记。 k, 显示不同非端脑神经元聚类中选定基因相对表达的热图。彩色点表示如j所示的聚类身份。Cb,小脑;ChP,脉络丛;CP,脉络丛;Hy,下丘脑;max.,最大;MB,中脑;MH,延髓;min.,最小;Oligo,少突胶质细胞;OPC,少突胶质前体细胞;PSC,多能干细胞;telen.,端脑;Th,丘脑;vTelen,腹侧端脑。 Para_02 为了阐明细胞状态和类型的动态变化和转变,我们基于神经最优传输(neural optimal transport)使用 moscot 重建了 HNOCA 细胞的真实年龄指导下的伪时间(Fig. 1h)。 聚焦于背侧端脑神经轨迹,我们观察到标志基因如 SOX2(神经前体细胞 (NPCs))、BCL11B(深层皮质神经元)和 SATB2(上层皮质神经元)的一致伪时间表达谱型(Fig. 1i)。 为进一步解析非端脑神经元中的异质性,我们对这一群体进行了亚聚类,揭示了许多由不同标志基因表达特征的神经元群体(Fig. 1j,k)。 HNOCA projection to a human developing brain atlas Para_01 为了评估我们的细胞类型注释,并更精确地注释异质性的非端脑神经元群体,我们将 HNOCA 与最近发表的人类大脑发育的单细胞转录组图谱进行了比较(图 2a)。 我们对主要参考图谱应用了 scVI 和 scANVI,并使用 scArches 将 HNOCA 投影到相同的潜在空间。 共享的潜在空间使我们能够重建 HNOCA 细胞和主要参考图谱之间的加权二分 k-近邻(wkNN)图,该图用于将‘CellClass’和‘Subregion’标签以及神经母细胞和神经元的神经递质转运蛋白(NTT)信息转移到 HNOCA。 转移的标签与我们分配的标签高度一致(扩展数据图 3),并使我们能够精炼 HNOCA 非端脑 NPCs 和神经元的区域注释,以及非端脑神经元的 NTT 注释(图 2b),最终形成了最终的分层 HNOCA 细胞类型注释(扩展数据图 3)。 Fig. 2: Projection of HNOCA to primary developing human brain cell atlases assists organoid neural cell type annotation and estimation of primary cell type representation.
a, 人类发育大脑细胞图谱27的UMAP,按NTT亚型(左)、区域(中)和注释的细胞类别(右)着色。 b, HNOCA的UMAP,按映射的神经元NTT亚型(左)和NPC、中间前体细胞(IP)及神经元的区域标签着色。 c, 热图显示不同年龄类器官中的细胞比例与来自不同主要发育阶段的细胞匹配的情况。 d, 不同数据集中代表不同区域(端脑、间脑、中脑和后脑)的神经细胞百分比。x轴显示数据集,按总比例(条形高度)降序排列。基于非指导分化协议的数据集在下方用点标记。每个面板底部的条形图显示了类器官协议类型。 e, 人类发育大脑细胞图谱27的UMAP,按HNOCA数据集中细胞群体的最大存在评分着色。低评分表示细胞状态在HNOCA数据集中代表性不足。 f, 人类参考图谱27中不同细胞类别的最大存在评分分布。Eryt.,红细胞;Imm.,免疫;Vas.,血管;G-blast,胶质母细胞;F-blast,成纤维细胞;NC,神经嵴;Plac.,板状;RG,放射状胶质;IPC,中间前体细胞;N-blast,神经母细胞;N,神经元。 g, 盒形图显示了不同主要参考细胞簇中最大存在评分的分布。底部侧边栏显示了神经元与非神经元、细胞类别、主要参考区域的信息。 h, 人类发育大脑图谱的UMAP,显示在HNOCA中代表性不足的主要神经细胞类型或状态(红色)。Hippo,海马;HyTh,下丘脑;d,背侧;v,腹侧;CB,小脑。 Para_02 我们还试图将类器官细胞与人类大脑发育的第一孕期之后的阶段进行比较。 专注于背侧端脑,我们将HNOCA神经前体细胞和神经元的转录组谱型与涵盖第一孕期到青春期的人类皮层发育初级图谱中的细胞进行了比较。 我们观察到从第一孕期观察到的细胞状态向第二孕期皮层中观察到的更成熟状态的转变(图2c),并且没有检测到与后期阶段的显著匹配。 我们使用两个初级图谱对其他脑区进行了扩展比较,分别代表了第一孕期和第二孕期。 我们确认,在较老的类器官中,与其他脑区的第二孕期细胞状态的相似性增加(扩展数据图3)。 Para_03 我们评估了每种神经类器官方案生成不同脑区神经细胞的能力(图2d,扩展数据图3和4及补充表2)。 未引导的神经类器官数据集包含来自所有脑区的细胞,但各数据集中的比例各不相同,这表明未引导方案能够生成许多脑区,但具有高度变异性。 相比之下,从引导性类器官方案衍生的数据集则富含目标脑区的细胞,但通常显示出邻近脑区细胞的比例增加。 例如,几个从中脑类器官方案衍生的数据集也显示了高比例的后脑神经元,这表明形态发生素引导的不精确性。 Para_04 为了全面评估由 HNOCA 代表的类器官协议生成初级脑细胞类型的效果,我们估计了每个 HNOCA 数据集中每种初级细胞类型的出现分数(方法)。 较高的出现分数表明该类型细胞在 HNOCA 数据集中出现的频率高且可能性大。 通过对每个类器官数据集的分数进行归一化处理(扩展数据图 5 和补充表 3),我们得到了一个指标,用以描述每种初级细胞类型在至少一个 HNOCA 数据集中被代表的程度(图 2d)。 这一分析证实了 HNOCA 中缺乏红细胞、免疫细胞和血管内皮细胞,这些细胞均来源于非神经外胚层的生殖层(图 2e)。 正如预期的那样,端脑细胞类型在 HNOCA 中表现最为强烈。 相比之下,丘脑、中脑和小脑的细胞类型表现最少,包括丘脑网状核 GABA 能神经元、背侧中脑 m1 来源的 GABA 能神经元和 m1/m2 来源的谷氨酸能神经元,以及小脑浦肯野细胞(图 2f,g)。 值得注意的是,尽管这些细胞类型在 HNOCA 数据集中的数量少于原始图谱中的数量,但某些类器官协议可以生成其中一些代表性不足的细胞类型(例如,在后脑类器官协议中生成的浦肯野细胞)。 Transcriptomic fidelity organoid cell types Para_01 我们接下来的目标是了解由不同分化协议产生的类器官之间的转录组相似性和差异性,以及类器官与原代脑组织之间的转录组相似性和差异性。 我们鉴定了差异表达基因(DEGs),通过比较HNOCA中的神经细胞类型与其原代对应物(图3a和补充表4)。 我们发现,对于大多数神经细胞类型,超过三分之一(平均34.4%,标准差12.1%)的DEGs在至少一半的协议中是共享的(协议通用DEGs),这表明许多类器官与原代细胞之间的转录组差异独立于类器官协议(图3b)。 我们使用一个额外的人类皮层单细胞RNA测序数据集验证了我们的结果(扩展数据图6和补充表5)。 接下来,我们评估了区域神经细胞类型之间共享的差异转录程序,并确定了总共920个普遍存在的、协议通用的DEGs(uDEGs),这些DEGs在至少16种神经细胞类型中的14种中差异表达(图3c)。 这些uDEGs在神经元类型和协议之间显示出一致的倍数变化(r > 0.8)(图3d),代表了无论协议或神经元细胞类型如何,类器官中的神经元与原代组织中的神经元之间的一致分子差异。 在所有920个uDEGs中,363个基因持续上调,673个基因持续下调,只有59个基因(6%)在亚型或协议间表现出不一致的差异表达(图3e)。 Fig. 3: Transcriptomic comparison between organoid neurons and their primary counterpart reveals universal cell stress in organoids.
a, 在 HNOCA 中不同协议产生的神经细胞类型与其原始对应物之间的差异表达分析示意图27。 b, 不同神经细胞类型中表达基因的比例,在生成相应亚型的某些协议比例中显示出差异表达。左上角,谷氨酸能神经元;右下角,GABA 能神经元。颜色表示大脑区域。 c, 协议常见差异表达基因(DE 在至少一半协议中)的数量,按基因在其中差异表达的神经细胞类型的数量分组。 d, 不同神经元亚型*协议(即每种不同协议生成的神经细胞类型)中普遍存在的 DE 基因的表达对数变化(logFC)相关性的分布。 e, 每个类别的 DE 基因数量。 f, 下调(上部,蓝色)和上调(下部,红色)普遍存在的 DE 基因的基因本体富集分析。方块的大小与 -log 转换的调整后 P 值相关。 g,h, 初级神经细胞类型(上部,深色)和类器官对应物(下部,浅色)中线粒体 ATP 合成耦合电子传输模块得分(g)、经典糖酵解模块得分(h,左侧)和分子特征数据库标志性糖酵解模块得分(h,右侧)的分布。P 值,双侧 Wilcoxon 测试的显著性。 i, 热图显示了三个模块得分的两两相关性(corr.)。 j, 背侧端脑兴奋性神经元(dTelen VGLUT-N)的标志性糖酵解得分,分为三个初级发育中的人类大脑和至少包含 20 个 dTelen VGLUT-N 的 27 个类器官数据集。下图显示了可能与细胞应激相关的分化协议的选定特征。协议和出版物索引如扩展数据图 1 所示。Mat. media,成熟培养基。 k, HNOCA 中神经细胞类型和人类发育大脑图谱27 中神经细胞类型之间的基因表达谱的斯皮尔曼相关性,跨越可变转录因子(TFs)。数据集顺序与补充表 1 相同。 Para_02 通过对 uDEGs 进行基因本体富集分析,我们发现下调的 uDEGs 富集在神经发育过程,包括神经元细胞间粘附和突触组织(图 3f)。 上调的 uDEGs 富集在许多与代谢相关的术语中,包括线粒体 ATP 合成耦合电子传递(简称电子传递)和经典糖酵解(图 3f)。 先前的研究已将能量相关通路的富集与当前培养条件限制引起的代谢变化联系起来。 此外,分子特征数据库中的糖酵解基因集已被用于定义神经类器官的代谢状态。 对 HNOCA 和主要参考图谱中的线粒体电子传递、经典糖酵解和标志糖酵解基因集进行评分,我们发现所有三个术语在类器官和原代细胞之间显示出显著的分离(图 3g,h)。 使用参考文献 3 和 27 的数据集作为代表性例子,我们确定了所有神经细胞类型中糖酵解评分的类似分布,并且总体上类器官细胞的评分增加(扩展数据图 7)。 专注于背侧端脑神经元,我们比较了不同类器官分化协议下的糖酵解评分分布,并识别出几种与代谢细胞应激相关的协议特征。 例如,使用成熟培养基、切片或切割类器官以及在较小程度上摇动或旋转类器官会导致总体较低的糖酵解评分(图 3h)。 类器官和原代参考细胞类型在不同分化协议下的平均糖酵解评分和转录组相似性呈负相关。 当仅考虑可变转录因子时,这种相关性显著降低,表明类器官中的代谢变化对神经细胞类型的核分子身份影响有限(扩展数据图 7)。 这一观察结果与先前研究一致,这些研究显示神经类器官中细胞的特定代谢状态与原发组织相比不会影响神经元的命运决定和成熟。 Para_03 接下来,我们专注于 366 个可变转录因子的表达,以计算 HNOCA 数据集和主要参考图谱之间相应神经元细胞类型的关联性。 我们发现,无论是指导性还是非指导性的类器官分化协议,生成的神经元细胞类型与相应的原始参考细胞类型的相似度相当。 然而,我们观察到转录组相似度存在脑区依赖性差异。 例如,大多数脑区背侧的类器官神经元在类器官数据集中与其原始对应物的相似度高于大多数脑区腹侧衍生的细胞类型(图 3i)。 Para_04 为了识别除代谢状态以外降低类器官保真度的分子特征,我们将来自四个不同的人类大脑发育图谱的背侧端脑谷氨酸能神经元整合为一个综合的主要参考,并确定了神经元亚型和成熟状态的异质性(扩展数据图8)。 将HNOCA中的背侧端脑神经元投影到主要图谱上揭示了神经类器官中的相应异质性。 在差异表达分析中将代谢状态、成熟状态和细胞亚型作为协变量显著减少了差异表达基因的数量,支持了这些是区分类器官和原代脑细胞的主要因素的观点(扩展数据图8和补充表6)。 我们观察到丰富的生物过程包括突触囊泡循环和高电压门控钙通道活性的负调控(扩展数据图8),这表明类器官在这些过程中存在缺陷。 值得注意的是,这些差异在不同的类器官协议中都有观察到,并突显了体外与原代对应物之间一致的转录组学差异。 HNOCA facilitates organoid protocol evaluation Para_01 HNOCA 及我们建立的分析管道提供了一个框架,用于查询未包含在 HNOCA 中的新神经类器官单细胞 RNA 测序(scRNA-seq)数据集。 为了展示这一应用,我们从最近发表的一项多重神经类器官形态发生筛选中检索了 scRNA-seq 数据,并将其投影到 HNOCA 和主要参考文献的潜在空间(图 4a,扩展数据图 9 和补充表 7)。 我们转移了区域标签,并发现与提供的区域注释高度一致,但在前脑、中脑和后脑等广泛脑区内部具有更高的分辨率(图 4b)。 因此,我们的转移注释允许更全面地评估不同形态发生条件对产生不同脑区神经元的影响(图 4c)。 我们进一步计算了每个筛选条件下参考细胞的存在分数,并将不同筛选条件的数据与 36 个 HNOCA 数据集进行了比较。 使用平均存在分数进行层次聚类显示许多筛选条件具有不同的存在分数谱型(图 4d),表明区域细胞类型组成与 HNOCA 数据集不同。 接下来,我们总结了整个形态发生筛选数据的最大存在分数(图 4e),并与 HNOCA 数据的存在分数进行了比较,以确定筛选中增加存在的主要参考细胞类型(图 4f)。 该分析突出了在某些筛选条件下显著增加丰度的几个参考细胞簇(图 4g),例如腹侧端脑中的 LHX6/ACKR3/MPPED1 三重阳性 GABA 能神经元和腹侧中脑中的多巴胺能神经元。 总之,将形态发生筛选查询数据投影到 HNOCA 和主要参考文献上,使得对形态发生筛选数据的精炼注释成为可能,同时也对新分化协议生成以前在神经类器官中代表性不足或缺乏的神经元细胞类型的值进行了全面和定量的评估。 Fig. 4: Projection of neural organoid morphogen screen scRNA-seq data to HNOCA and human developing brain atlas allows cell type annotation and organoid protocol evaluation.
a, 投射神经器官形态发生筛选44单细胞RNA测序数据到HNOCA和人类发育大脑参考图谱27。UMAP显示了筛选条件组(左,使用形态发生素SAG(声波刺猬信号激动剂)、CHIR、BMP和FGF)和筛选数据的区域注释(右)。 b, 筛选数据的区域注释(行)与从主要参考数据转移的区域标签的比较。 c, 基于参考投影分配给不同区域的细胞比例。每个堆叠条形图代表一个筛选条件。 d, 基于主要参考数据中聚类的平均存在分数,对HNOCA数据集和筛选数据中的条件进行聚类。热图显示了主要参考数据中每个聚类的平均存在分数(列)。 e, 主要参考数据的UMAP按解剖区域着色(右)和筛选条件下最大存在分数(左)。 f, 相对于HNOCA数据集,筛选条件下细胞聚类覆盖率的增加,负值修剪为零。灰色水平线显示了定义筛选数据中获得的聚类的阈值(0.3)。 g, 主要参考数据的UMAP,获得的聚类以蓝色调突出显示。虚线圆圈突出了在端脑和中脑中覆盖率最高的两个聚类。 h, 在主要参考数据(上)和筛选数据集(下)中,g中标记的两个聚类的共表达分数。DA,多巴胺能。 HNOCA facilitates disease model interpretation Para_01 接下来,我们测试了整合的 HNOCA 是否可以作为评估神经疾病类器官模型的对照队列。 我们从10种神经类疾病类器官模型及其各自的对照组(小头症、肌萎缩侧索硬化症、阿尔茨海默病、自闭症、脆性X综合征(FXS)、精神分裂症、神经异位症、皮特-霍普金斯综合症、强直性肌营养不良和胶质母细胞瘤)中收集了11个单细胞RNA测序(scRNA-seq)数据集(图5a,扩展数据图10和补充表8)。 我们将数据投影到HNOCA和主要参考图谱上,以转移注释(图5b-f)。 我们发现,在大多数研究中,疾病模型类器官与其各自的研究特定对照类器官之间存在细胞类型和脑区组成的差异(图5g,h)。 这些差异可能代表了疾病表型,但也可能是细胞系变异的结果。因此,正确标注疾病和对照类器官的细胞类型和区域组成对于识别疾病表型至关重要,特别是在分析给定细胞类型中的疾病相关转录组变化时。 Fig. 5: The HNOCA as a control cohort to facilitate cell type annotation and transcriptomic comparison for neural organoid disease-modelling data.
a,疾病建模神经类器官图谱构建概述,以及向主要图谱27和HNOCA进行下游分析的投影。 b-f,整合的疾病建模神经类器官图谱的UMAP图,按预测的细胞类型注释(b)、NPCs、中间前体细胞和神经元的预测区域身份(c)、出版物(d)、疾病状态(e)和标记基因表达(f)着色。 g, h,分配给不同细胞类别(g)和区域(h)的细胞比例。 每个堆叠条形代表一个生物样本。 侧边条形显示疾病状态和出版物。 i,为疾病建模神经类器官图谱中的每个细胞重建匹配的HNOCA元细胞的示意图。 j,疾病建模神经类器官图谱的UMAP图,按与匹配的HNOCA元细胞的转录组相似性着色。 k,小提琴图表示按出版物划分的估计转录组相似性的分布。 左,对照细胞的分布;右,疾病细胞的分布。 l,热图显示GBM-2019数据集中AQP4+群体与其匹配的HNOCA元细胞之间顶级差异表达基因(DEGs)的表达。 行显示表达减少和增加最强烈的十个DEGs。 列显示疾病建模样本中AQP4+群体的平均表达(第一面板),每一样本的匹配HNOCA元细胞(第二面板),所有预测的对照星形胶质细胞和HNOCA中所有星形胶质细胞的平均表达。 m,火山图显示FXS-2021数据集中背侧端脑细胞与其匹配的HNOCA元细胞之间的差异表达分析(DE)。 DEGs用红色(在FXS中增加)和蓝色(在FXS中减少)着色。 圈出的点显示在SFARI数据库中注释的DEGs。 顶部栏显示了SFARI基因在增加(红色)和减少(蓝色)DEGs中的对数转换优势比。 GBM,胶质母细胞瘤。 Para_02 我们开发了一种基于wkNN的策略,为每个疾病模型类器官单细胞RNA测序数据集中的每个细胞生成匹配的HNOCA元细胞(图5i),并量化了它们的转录组相似性(图5j)。 胶质母细胞瘤类器官数据集显示,与原发对应物相比,其相似性显著低于其他疾病模型(图5k)。 为了评估这些转录组差异,我们在胶质母细胞瘤和匹配的对照元细胞之间进行了差异表达分析。 专注于AQP4+群体(扩展数据图10),我们确定了1,951个胶质母细胞瘤细胞相对于匹配的HNOCA元细胞的差异表达基因(补充表9),并发现RBM25、CALD1、HNRNPU和SPARC等基因的表达增加(图5l),所有这些基因都被报道与胶质母细胞瘤相关。 Para_03 接下来,我们专注于FXS58的类器官模型,在该模型中,对照组类器官中的NPCs和神经元具有非端脑身份,而疾病模型类器官主要包含端脑细胞(图5h和扩展数据图10)。 整合的HNOCA提供了对FXS新皮层神经元进行差异表达分析的机会,与匹配的HNOCA元细胞相比,鉴定了444个差异表达基因(DEGs)。 在FXS细胞中高表达的DEGs(122个基因)富含在西蒙斯基金会自闭症研究倡议(SFARI)数据库中标注的自闭症相关基因。 其中一个基因CHD2,在原始出版物中被报道为FXS的关键调节因子,其蛋白质水平增加,但在批量RNA测序实验中未能检测到其信使RNA(mRNA)水平的变化。 我们还检测到FMR1的表达降低,其功能丧失突变导致FXS。 Extending the HNOCA through data projection Para_01 新的单细胞 RNA 测序(scRNA-seq)数据集不断产生,重要的是要持续扩展和更新人类神经类器官单细胞图谱(HNOCA),以包含这些额外的数据。 因此,我们建立了一个计算工具包,用于将新的 scRNA-seq 数据投影到 HNOCA(图 6a)。 我们通过将来自另外六项研究的 scRNA-seq 数据纳入 HNOCA(扩展版 HNOCA;图 6b 和补充表 10),使用查询到参考的映射方法来演示该工具包的使用。 我们使用基于加权 k 近邻(wkNN)的标签转移来统一细胞类型注释,并将这些细胞置于由 HNOCA 表征的现有类器官单细胞转录组景观的背景下(图 6c-e)。 使用我们的方法将更多数据集映射到 HNOCA,通过增加对现有神经类器官协议和在类器官中产生的神经细胞类型的覆盖范围,增强了图谱。 Fig. 6: Extending the HNOCA by means of projection of extra datasets.
a, 社区进一步扩展 HNOCA 的单细胞 RNA 测序数据投影示意图。 b, UMAP 显示当前扩展的 HNOCA 的数据集组成。 c, UMAP 显示五个扩展数据集中细胞的预测细胞类型注释。NE,神经上皮;NC-D,神经嵴衍生物;MC,间充质细胞;EC,内皮细胞。 d,点图显示了扩展的 HNOCA 数据集中预测细胞类型中选定的细胞类型和区域标记的表达情况。 e,点图显示了扩展数据集中细胞类型的组成以及与匹配的 HNOCA 元细胞的平均相似性。 f,示意图展示了分析管道和不同的界面,以帮助社区分析神经类器官的单细胞 RNA 测序数据。 Para_02 为了使研究人员能够在自己的分析中使用 HNOCA,我们提供了多种探索和与图谱互动的选项(图 6f)。 HNOCA 可以通过一个在线门户进行浏览,实现基因表达的可视化和标记基因的发现。 我们还通过一个在线界面(http://www.archmap.bio/)提供 HNOCA,用于新数据集的交互映射,支持标签转移、存在分数计算和细胞状态的代谢评分。 最后,我们开发了 HNOCA-tools,这是一个 Python 包,实现了本文介绍的所有核心分析方法,如注释、参考映射、标签转移和差异表达测试方法。 Discussion Para_01 在这项研究中,我们通过整合由全球15个不同实验室使用26种不同的分化协议以及多样化的单细胞RNA测序技术生成的36个单细胞RNA测序数据集中的180万个细胞,构建了人类神经类器官的大规模整合细胞图谱,即HNOCA。 由此产生的图谱揭示了在现有协议条件下培养的神经类器官中可以发展的神经元、胶质和非神经细胞类型的高复杂性。 将HNOCA数据映射到各种人类发育大脑细胞参考图谱,允许对神经类器官协议进行全面评估,以生成不同脑区的细胞类型。 我们发现,前3个月培养的类器官最匹配第一孕期的主要数据,而大约3个月或更长时间培养的类器官则最匹配第二孕期的主要细胞状态。 我们没有观察到与较老发育阶段相匹配的显著神经元成熟和多样化特征,这表明当前神经类器官协议中存在神经元成熟的限制。 Para_02 我们在类器官神经元类型和它们的主要对应物之间进行了差异表达分析,以评估转录组保真度,并确定与糖酵解途径相关的代谢变化是区分类器官和主要细胞状态的主要因素,这与之前的报道一致。 尽管代谢压力对整体转录组保真度有负面影响,但区域细胞类型的分子身份得以保持,这一点通过与主要对应物高度一致的转录因子共表达模式得到了证明。 Para_03 我们展示了查询数据的映射,即将最近发表的单细胞转录组神经类器官形态发生筛选数据映射到HNOCA和主要参考文献,这使得细胞类型注释更加精确,并且能够与现有的神经类器官数据集进行组成比较。 我们强大的框架将有助于人类神经类器官的scRNA-seq数据的定量和比较分析,以及新神经类器官协议的基准测试。 Para_04 与之前的报告一致,我们发现无指导协议生成的神经细胞具有较高的脑区变异性,这在研究神经发育期间更广泛的命运决定时是有用的。 指导性协议导致目标脑区的显著富集。 我们还注意到,一些针对中脑的指导性协议显示出相对较低的特异性,并生成来自邻近脑区的神经细胞。 这个问题可能是由于类器官中神经干细胞对相同形态发生素信号的差异响应,或者是由于缺乏对定义中枢神经系统深层区域细胞所需的时间、浓度和形态发生素组合的完全理解。 Para_05 集成的 HNOCA 也是分析疾病建模神经类器官数据的优秀资源。 它有助于细胞类型注释,并提供了一个大规模的单细胞转录组对照队列用于比较。 例如,在许多研究中,我们观察到对照样本和疾病模型样本之间在细胞类型和区域组成上存在差异。 同时,HNOCA 提供了机会,可以在多线多协议的大规模对照队列背景下识别特定疾病的分子特征。 Para_06 我们展示了如何通过将额外的神经类器官单细胞转录组数据投影到图谱上来扩展和更新HNOCA。 此外,我们开发了一套计算工具HNOCA-tools,这将使其他研究人员能够重现我们在研究中应用的分析框架。 共同地,我们设想HNOCA将保持更新,并继续反映体外生成的人类神经细胞状态的景观,作为神经类器官社区的动态资源,使评估类器官保真度、表征扰动和疾病状态以及开发新协议成为可能。 Methods Metadata curation and harmonization of human neural organoid scRNA-seq datasets 人类神经类器官单细胞 RNA 测序数据集的元数据整理和 harmonization
Para_01 我们纳入了来自25篇出版物的33个人类神经类器官数据,加上三个未发表的数据集纳入我们的图谱(补充表1)。 我们通过 sfaira78 框架(GitHub 开发分支,2023年4月18日)整理了本研究中使用的所有神经类器官数据集。 为此,我们从每个包含出版物的数据可用性部分提供的位置或在未发表数据的情况下直接从作者处获取单细胞 RNA 测序计数矩阵和相关元数据。 我们根据 sfaira 标准(https://sfaira.readthedocs.io/en/latest/adding_datasets.html)统一了元数据,并手动整理了一个额外的元数据列 organoid_age_days,描述了类器官在采集前培养的天数。 Para_02 我们接下来移除了已发布数据集中不适用的子集:患病样本或表达疾病相关突变的样本,融合类器官,原代胎儿数据,激素处理的样本,在神经诱导前收集的数据和 share-seq 数据。 我们将所有剩余数据集统一到一个共同的特征空间,使用 ensembl79 版本 104 中的‘蛋白编码’或‘长链非编码 RNA’生物类型的任何基因,并用零填充任何数据集中缺失的基因。 平均而言,每个组成数据集报告了全部基因空间(36,842 个基因)的 50%。 然后,我们将所有剩余数据集合并创建了一个单一的 AnnData 对象。 Preprocessing of the HNOCA scRNA-seq data HNOCA 单细胞 RNA 测序数据的预处理
Para_01 所有处理和分析均使用 scanpy(v.1.9.3)进行,除非另有说明。 对于 HNOCA 的质量控制和过滤,我们移除了表达少于 200 个基因的任何细胞。 接下来,根据两个质量控制指标:表达的基因数量和线粒体计数百分比,移除异常值细胞。 为了基于每个质量控制指标定义异常值细胞,首先对所有细胞的值应用 z 变换。 任何 z 变换后的指标小于 -1.96 或大于 1.96 的细胞被定义为异常值。 对于任何使用 10X Genomics 的 v.3 化学试剂收集的数据集,在过滤后包含超过 500 个细胞的情况下,我们将高斯分布拟合到表示每个细胞表达基因数量的直方图上。 如果检测到双峰分布,我们将移除任何表达基因数量少于由分布两个峰值之间的谷值定义的细胞。 然后,我们通过将原始读取计数除以从 BioMart 获得的每个基因的最大基因长度来归一化所有 Smart-seq2 数据的原始读取计数。 接下来,我们将这些归一化的读取计数乘以数据集中所有基因的中位基因长度,并将这些长度归一化的计数等同于下游分析中使用唯一分子标识符获得的数据集的原始计数。 Para_02 下一步,我们通过首先将每个细胞的计数除以该细胞的总计数,然后乘以1,000,000,最后对每个计数+1取自然对数,生成了一个对数标准化表达矩阵。 我们使用scanpy的highly_variable_genes函数(flavor = ‘seurat_v3’,batch_key = ‘bio_sample’)以批次感知的方式计算了3,000个高变特征。 在这里,bio_sample表示数据集中原始元数据提供的生物样本。 平均而言,3,000个高变基因中有72%被报告在每个构成的HNOCA数据集中。 我们使用这3,000个特征通过主成分分析(PCA)计算了数据的50维表示,进而使用该表示计算了一个k近邻(kNN)图(n_neighbors = 30,metric = ‘cosine’)。 使用邻近图,我们通过UMAP计算了数据的二维表示,并使用Leiden聚类方法对未集成的数据进行了粗略(分辨率1)和精细(分辨率80)的聚类。 Hierarchical auto-annotation with snapseed 使用 Snapseed 进行分层自动标注
Para_01 Snapseed 是一种可扩展的自动注释策略,它根据提供的细胞类型层次结构及相应的细胞类型标记物对细胞进行注释。 该方法基于细胞聚类中标记基因表达的富集(优选高分辨率聚类),并且不一定需要数据整合。 Para_02 在这项研究中,我们使用 snapseed 获得标签感知集成的初始注释。 首先,我们构建了一个包括祖细胞、神经元和非神经类型在内的细胞类型层次结构,每种类型由一组标记基因定义(补充数据 1)。 接下来,我们通过最近发布的发育中人脑细胞图谱中的细胞簇平均表达谱来表示数据(RSS3)。 然后,我们在 RSS 空间中构建了一个 kNN 图(k = 30),并使用 Leiden 算法(分辨率 80)对数据集进行聚类。 对于这两个步骤,我们使用了通过 scanpy 提供的图形处理单元(GPU)加速的 RAPIDS 实现。 Para_03 对于层次结构中给定层级的所有细胞类型标记基因,我们计算了接收者操作特征曲线下的面积(AUROC)以及跨聚类的检出率。 对于每种细胞类型,通过将其标记基因的最大 AUROC 与最大检出率相乘来计算得分。 每个聚类随后被分配给得分最高的细胞类型。 此过程递归地应用于层次结构的所有层级。 使用未整合数据的精细(分辨率 80)聚类执行相同的过程,以获得用于下游作为整合方法基准测试的地面真实输入的未整合数据集的细胞类型标签。 Para_04 这种自动注释策略在 snapseed Python 包中实现,并在 GitHub 上可用(https://github.com/devsystemslab/snapseed)。 Snapseed 是一个轻量级包,用于支持图谱级别的数据集中的可扩展基于标记的注释,其中手动注释不可行。 该包实现了三个主要功能:annotate() 用于使用定义的标记基因对细胞类型列表进行非层次注释,annotate_hierarchy() 用于注释更复杂的手动定义的细胞类型层次结构,find_markers() 用于快速发现特定集群的特征。 所有功能都基于使用 JAX(https://github.com/google/jax)实现的 GPU 加速 AUROC 分数。 Label-aware data integration with scPoli 使用scPoli进行标签感知的数据整合
Para_01 我们使用 scArches 包中的 scPoli 模型对 HNOCA 的类器官数据集进行了整合。 我们将用于整合的批次协变量定义为数据集标识符(注释列‘id’)、生物学重复的注释(注释列‘bio_sample’)以及技术重复的注释(注释列‘tech_sample’)的组合。 这导致了 396 个单独的批次。 批次协变量在模型中表示为一个大小为五的学习向量。 我们使用基于 RSS 的 snapseed 细胞类型注释的前三个级别作为 scPoli 原型损失的细胞类型标签输入。 我们选择了一层 scPoli 编码器和解码器的隐藏层大小为 1,024,潜在嵌入维度为十。 我们使用 100 作为‘alpha_epoch_anneal’参数的值。 我们没有使用未标记的原型预训练。 我们总共训练了七轮模型,其中五轮是预训练轮。 Benchmark of data integration methods 数据集成方法的基准测试
Para_01 为了定量比较几种工具的类器官图谱整合结果,我们使用了 GPU 加速的 scib-metrics Python 包(版本 v.0.3.3),并使用了总体性能最高的嵌入进行所有下游分析。 我们比较了以下数据潜在表示的数据整合性能:未整合的 PCA,RSS3 整合,scVI(默认参数,但使用两层、30 维潜在空间和负二项似然)整合,scANVI(默认参数)整合,使用 snapseed 级别 1、2 或 3 注释作为细胞类型标签输入,scPoli 整合,使用 snapseed 级别 1、2 或 3 注释或同时使用所有三个注释级别作为细胞类型标签输入,以及使用 aggrecell 算法聚合的元细胞的 scPoli 整合(首先用作‘伪细胞’),使用 snapseed 级别 1 或 3 注释作为 scPoli 的细胞类型标签输入。 我们使用了以下评分来确定整合质量(每个评分在参考文献 46 中都有描述):Leiden 归一化互信息评分,Leiden 调整后的 Rand 指数,每种细胞类型标签的平均轮廓宽度,孤立标签评分(基于平均轮廓宽度评分)和细胞类型局部逆 Simpson 指数,以量化生物变异性的保持情况。 为了量化批次效应的去除,我们使用了每批标签的平均轮廓宽度,整合局部逆 Simpson 指数,kNN 批次效应测试评分和图连通性。 然后根据单独归一化(归一化到 [0,1] 范围内)的指标的总分对整合方法进行了排名。 在进行基准测试之前,我们迭代地从数据集中移除了任何具有与其他细胞相同潜在表示的细胞,直到数据集中的潜在表示不再包含任何重复行。 此过程共移除了 3,293 个重复细胞(占整个数据集的 0.002%),这是为了使基准测试算法能够无错误地完成。 我们在整合中使用了在未整合的 PCA 嵌入上计算的 snapseed 级别 3 注释作为真实的细胞类型标签。 Pseudotime inference 伪时间推断
Para_01 为了推断分化状态的全局排序,我们试图基于神经最优传输在 scPoli 隐空间中推断一个实时信息伪时间。 我们首先将类器官的年龄按天数分为七个区间((0, 15],(15, 30],(30, 60],(60, 90],(90, 120],(120, 150],(150, 450])。 接下来,我们使用 moscot 解决了一个时间神经问题。 为了根据预期增殖率对边缘分布进行评分,我们使用 score_genes_for_marginals() 方法获得了每个细胞的增殖和凋亡分数。 然后计算了边缘权重。 Para_02 最优传输问题使用了以下参数:迭代次数 = 25,000,计算瓦瑟斯坦基线 = False,批处理大小 = 1,024,耐心等待次数 = 100,预训练 = True,训练集大小 = 1。 为了计算每个年龄区间 i 中细胞的位移向量,我们使用了对应于 [i, i + 1] 传输图的子问题,除了最后一个年龄区间,我们使用了 [i − 1, i] 的子问题。 位移向量是通过从运输后的分布中减去原始细胞分布获得的。 使用 CellRank 提供的速度核,我们从位移向量计算了一个转移矩阵,并将其作为计算扩散图的输入。 负扩散成分 1 的排名被用作伪时间顺序。 Preprocessing of the human developing brain cell atlas scRNA-seq data 人类发育大脑细胞图谱单细胞 RNA 测序数据的预处理
Para_01 主要图谱27的细胞测序(scRNA-seq)数据经过 Cell Ranger 处理后,从其 GitHub 页面提供的链接(https://storage.googleapis.com/linnarsson-lab-human/human_dev_GRCh38-3.0.0.h5ad)获取。 为进一步的质量控制,过滤掉检测到基因少于300个的细胞。 转录本计数通过该细胞的总计数进行归一化,乘以10,000的比例因子,然后进行自然对数转换。 特征集与类器官图谱中检测到的所有基因相交,并使用 scanpy 函数 highly_variable_genes 选择变异最大的2,000个基因,其中‘供体’作为批次键。 创建了一个额外的‘neuron_ntt_label’列,用于表示从细胞簇元数据的‘AutoAnnotation’列中派生出的自动分类神经递质转运蛋白亚型标签(https://github.com/linnarsson-lab/developing-human-brain/files/9755350/table_S2.xlsx)。 Reference mapping of the organoid atlas to the primary atlas 类器官图谱与原发图谱的参考映射
Para_01 为了将我们的类器官图谱与人类初级发育大脑的数据进行比较,我们使用了 scArches 将其投影到上述的人类初级大脑单细胞 RNA 测序图谱上。 我们首先在初级图谱上预训练了一个 scVI 模型,其中‘供体’作为批次键。 该模型的构建参数如下:n_latent = 20, n_layers = 2, n_hidden = 256, use_layer_norm = ‘both’, use_batch_norm = ‘none’, encode_covariates = True, dropout_rate = 0.2,并且以 1,024 的批量大小训练了最多 500 个 epoch,同时设置了提前停止标准。 接下来,使用‘亚区域’和‘细胞类别’作为细胞类型标签,通过 scANVI 对模型进行了微调,以 1,024 的批量大小训练了最多 100 个 epoch,同时设置了提前停止标准和每标签样本数为 100。 为了将类器官图谱投影到初级图谱上,我们使用了由 scvi-tools 提供的 scArches 实现。 查询模型以 1,024 的批量大小进行了微调,最多 100 个 epoch,同时设置了提前停止标准和权重衰减为 0.0。 Bipartite weighted kNN graph reconstruction 二部加权k近邻图重构
Para_01 将主要参考文献和查询(HNOCA)数据投影到同一潜在空间后,通过识别每个查询细胞在参考数据中的100个最近邻来构建一个无权二分k近邻图,使用Python中的PyNNDescent或RAPIDS-cuML(https://github.com/rapidsai/cuml),具体取决于GPU加速的可用性。 同样,通过识别参考数据中每个参考细胞的100个最近邻,也构建了一个参考k近邻图。 对于参考-查询二分图中的每条边,两个连接细胞的参考邻居之间的相似性,分别定义为A和B,由Jaccard指数表示: Para_02 然后将 Jaccard 指数的平方分配为边的权重,以获得参考数据集和查询数据集之间的加权二分 kNN 图。 wkNN-based primary developing brain atlas label transfer to HNOCA cells 基于wkNN的主要发育大脑图谱标签转移到HNOCA细胞
Para_01 给定主要参考文献和查询(HNOCA)之间的wkNN估计,可以通过多数投票将参考文献的任何分类元数据标签转移到查询细胞上。 简而言之,对于每个类别,其支持度是通过计算链接到该类别的参考细胞的边的权重总和来确定的。 支持度最大的类别将被分配给查询细胞。 Para_02 为了获得非端脑神经前体细胞和神经元的最终区域标签,以及非端脑神经元的NTT标签,在转移过程中增加了约束。 对于区域标签,只有非端脑区域,即间脑、下丘脑、丘脑、中脑、中脑背侧、中脑腹侧、后脑、小脑、桥和延髓,被认为是有效的类别进行转移。 标签转移过程仅应用于HNOCA中的非端脑神经前体细胞和神经元。 在进行多数投票之前,HNOCA中所有非端脑神经前体细胞和神经元的有效类别的支持分数通过随机游走重启过程进行了平滑(重启概率alpha为85%)。 接下来,应用了考虑结构层次的分层标签转移。 首先,将考虑的区域分为间脑、中脑和后脑,每个结构的支持分数是其子结构分数的总和。 使用多数投票将每个细胞分配到这三个结构之一。 然后,对分配的结构下的子结构(例如,间脑下的下丘脑和丘脑)进行第二次多数投票。 Para_03 对于NTT标签,我们首先在参考文献中基于提供的NTT标签识别出有效的区域-NTT标签对,这些标签来自参考神经母细胞和神经元簇及其最常见的区域。 在这里,最常见的区域以分层方式重新估计到上述最精细的分辨率。 接下来,在转移NTT标签时,对于每个具有区域标签转移的非端脑神经元,仅考虑该区域认为有效的NTT标签进行多数投票。 Stage-matching analysis 阶段匹配分析
Para_01 为了将 HNOCA 中的端脑神经前体细胞和神经元与发育阶段相匹配,我们使用了最近发布的皮质发育图谱作为参考。 处理过的单核 RNA 测序数据从其数据门户(https://cell.ucsf.edu/snMultiome/)获得。 鉴于提供的元数据中的‘类’、‘亚类’和‘类型’标签作为注释,以及‘个体’作为批次标签,应用了 scPoli 进行标签感知的数据整合。 接下来,根据不同的发育阶段对数据进行了分割。 对于每个阶段,基于 scPoli 隐含表示应用了 Louvain 聚类(分辨率,5)。 所有阶段的聚类被合并,并根据变异系数识别高变基因,具体方法见此页面:https://pklab.med.harvard.edu/scw2014/subpop_tutorial.html。 最后,HNOCA 端脑神经前体细胞和神经元与每个聚类在已识别的高变基因上进行相关性分析。最佳相关的聚类的阶段标签被分配给查询的 HNOCA 细胞。 Para_02 为了将分析扩展到其他神经元细胞类型,还引入了第二孕期多区域人脑图谱。 处理后的计数矩阵和元数据从 NeMO 数据门户获得(https://data.nemoarchive.org/biccn/grant/u01_devhu/kriegstein/transcriptome/scell/10x_v2/human/processed/counts/)。 以提供的元数据中的‘cell_type’标签作为注释,‘individual’作为批次标签,运行 scPoli 进行标签感知的数据整合。 对 scPoli 隐含表示应用 Louvain 聚类以识别聚类(分辨率,20)。 同样,基于我们先前生成的 scANVI 隐含表示,以分辨率为 20 的 Louvain 聚类也应用于第一孕期多区域人脑图谱。 计算所有聚类的平均表达谱,并使用与上述相同的过程识别高度可变基因,适用于两个主要图谱组合的聚类。 接下来,HNOCA 中的每个 NPC 和神经元都与这些聚类的平均表达谱相关联。 确定了最佳相关的第一孕期和第二孕期聚类以及相关性。 两种相关性之间的差异用作指标,表明 HNOCA 中 NPC 和神经元的阶段匹配偏好。 Presence scores and max presence scores of cells in the primary developing brain atlas 初级发育大脑图谱中细胞的存在分数和最大存在分数
Para_01 给定一个参考数据集和一个查询数据集,存在分数是分配给参考数据集中每个细胞的分数,描述了该参考细胞的细胞类型或状态出现在查询数据中的频率或可能性。 在这项研究中,我们计算了每个 HNOCA 数据集中主要图谱细胞的存在分数,以量化我们在每个 HNOCA 数据集中看到由每个主要细胞代表的细胞类型或状态的频率。 Para_02 具体来说,对于每个 HNOCA 数据集,我们首先将 wkNN 图仅限于该数据集中的 HNOCA 细胞。 接下来,计算初级图谱中每个细胞的原始加权度,作为与该细胞相连的剩余边权重之和。 然后应用带有重启的随机游走过程,在初级图谱的 kNN 图上平滑原始得分。 简而言之,我们首先将初级图谱的 kNN 图表示为其邻接矩阵(A),然后进行行归一化,将其转换为转移概率矩阵(P)。 将原始得分表示为向量 s0,在每次迭代 t 中,我们生成 st 作为 Para_03 此过程执行了100次,以获得平滑的存在分数,随后进行了对数转换。 低于第5百分位或高于第95百分位的分数被裁剪。 裁剪后的分数被归一化到[0,1]范围内,作为HNOCA数据集中的最终存在分数。 Para_04 鉴于每个 HNOCA 数据集中的最终存在分数,整个 HNOCA 数据中的最大存在分数随后被轻松计算为初级图谱中每个细胞的所有存在分数的最大值。 较大的(接近于一)最大存在分数表明该细胞类型或状态至少在一个 HNOCA 数据集中出现频率较高,而较小的(接近于零)最大存在分数则表明在所有 HNOCA 数据集中代表性不足。 Cell type composition comparison among morphogen usage using scCODA 使用 scCODA 比较形态发生素使用中的细胞类型组成
Para_01 为了测试从不同类器官分化方案中引入某些形态发生素时细胞类型组成的改变,我们使用了 pertpy90 实现的 scCODA 算法。 scCODA 是一种用于检测单细胞 RNA 测序数据组成变化的贝叶斯模型。 为此,我们从每个分化方案中提取了关于添加的形态发生素的信息,并根据它们在神经分化中的作用将它们分组为 15 个广泛的分子组(补充表 1)。 这些分子组被用作模型中的协变量。 从主要图谱转移过来的区域标签被用作分析中的标签(细胞类型标识符)。 对于没有区域身份的细胞类型,使用了图 1c 中呈现的细胞类型标签。 由于多能干细胞和神经上皮细胞主要存在于早期类器官阶段,因此将它们从分析中移除。 我们使用生物样本作为样本标识符。 我们使用默认参数顺序运行了 scCODA,采用无转弯采样(run_nuts 函数),并将每种细胞类型依次作为参考。 我们使用多数投票系统来确定在超过一半的迭代中可信地发生变化的细胞类型。 Cell type composition comparison among morphogen usage using regularized linear regression 使用正则化线性回归比较形态发生素使用中的细胞类型组成
Para_01 为了补充使用 scCODA 进行的组成分析,我们设计了一种替代方法,使用正则化线性回归测试差异组成。 我们将区域组成矩阵作为响应变量 Y,并将分子使用情况作为自变量 X,拟合了一个广义线性模型: Para_02 该模型使用 Lasso 正则化(alpha = 1)和高斯噪声及恒等链接函数进行拟合。 正则化参数 lambda 通过 glmnet92 R 包中的 cv.glmnet() 函数实现的交叉验证自动确定。 所有非零系数 β 均被视为富集和耗竭的指示。 DE analysis between HNOCA neural cell types and their primary counterparts and functional enrichment analysis HNOCA神经细胞类型与其原代对应物之间的差异表达分析及功能富集分析
Para_01 为了研究类器官和原代细胞之间的转录组差异,我们使用最终一级注释将 HNOCA 限定为标记为‘神经元’的细胞。 我们还从人类发育大脑图谱中选取了在神经元_ntt_label 注释列中被分配有效标签的细胞。 我们添加了两个额外的数据集,分别是来自胎儿皮质细胞的数据集,来源于参考文献 39 和参考文献 28。 对于参考文献 39 的数据,我们将数据限定为标记为‘胎儿’的细胞,并使用 RSEM93 根据 STAR94 映射结果估计每个细胞中每百万读取的转录本数量。 然后,我们使用 scanpy 计算了主成分分析(PCA)、k近邻图(kNN graph)、UMAP 和 Leiden 聚类(分辨率 0.2)。 随后,我们选择了 STMN2 和 NEUROD6 表达最高的聚类作为皮质神经元聚类,并仅使用这些细胞。 对于参考文献 28 的数据,我们将数据集限定为在其出版物补充表 5(‘皮质注释’)中标注为‘神经元’的细胞,并计算了 PCA、邻域图和 UMAP 以可视化数据集。 我们发现只有来自个体 CS14_3、CS20、CS22 和 CS20 的样本含有可检测的 STMN2 和 NEUROD6 表达,因此我们将数据集进一步限定为仅来自这些个体的细胞。 Para_02 为了计算 HNOCA 细胞与其主要对应细胞之间的差异表达(DE),我们首先通过汇总每种区域神经细胞类型中的细胞数量来创建伪批量样本(注释列:HNOCA 的‘batch’;人类发育大脑图谱的‘SampleID’;参考文献 39 的‘sample’和参考文献 28 的‘individual’)使用 decoupler95 的 Python 实现(版本 1.4.0),同时丢弃少于十个细胞或总计数少于 1,000 的任何样本。 然后我们将特征空间缩小到所有数据集特征的交集,并移除表达少于 200 个基因的任何细胞。 我们进一步移除了在 HNOCA 中少于 1% 的神经元表达的任何基因以及位于 X 和 Y 染色体上的任何基因。 剩余的 11,636 个基因中,平均每个构成的 HNOCA 数据集中报告了 99% 的基因。 对于每种区域神经细胞类型,我们从伪批量数据中移除了与少于两个总样本或少于 100 个总细胞的类器官分化测定相关的任何样本。 接下来,我们使用 edgeR 来迭代计算每种区域神经细胞类型下,每个类器官分化协议与相应的主要细胞之间的差异表达基因,同时校正类器官年龄(天数)、每个伪批量样本的细胞数量、每个伪批量样本检测到的基因数量的中位数和标准差。 我们使用上述提到的人类发育大脑图谱(参考文献 27)、参考文献 28 和参考文献 39 的数据作为‘背侧端脑神经元 NT-VGLUT’细胞类型的 DE 比较的主要数据,而对所有其他细胞类型,我们使用人类发育大脑图谱作为胎儿数据集。 我们使用 edgeR 的基因特异性负二项广义线性模型与准似然 F 检验进行分析。 我们认为如果一个基因的错误发现率(Benjamini–Hochberg 校正)P 值小于 0.05 且绝对 log2 折叠变化大于 0.5,则该基因显著差异表达。 我们使用 GSEApy Python 包对我们的 DE 结果进行功能富集分析,使用的是‘GO_Biological_Process_2021’基因集。 Para_03 为了评估不同主要数据集对差异表达(DE)结果的影响,我们计算了使用参考文献6的协议生成的HNOCA子集中背侧端脑神经元NT-VGLUT与Braun等人27的主要数据集中相应细胞类型以及参考文献28的数据之间的DE。 为了避免技术效应影响这一分析,我们在比较中仅使用了通过10X Genomics 3' v.2协议生成的细胞。 我们如上所述生成了伪批量样本,并在DE比较中校正了类器官年龄(天数)和每个伪批量样本中的细胞数量。 我们使用了上述描述的相同基于edgeR的程序和截止值。 我们使用scipy的fcluster方法根据两个主要数据集中基因的对数倍变化来聚类基因。 我们将聚类分组以表示一致上调、一致下调以及三个不同调控不一致的基因组。 我们如上所述计算了每个基因组的功能富集。 Para_04 为了评估不同类器官数据集对基于协议的差异表达分析的影响,我们计算了每个类器官出版物中的背侧端脑神经元 NT-VGLUT(如果一个出版物使用了多个协议,则进一步按协议划分)与参考文献27数据集中相应细胞类型之间的差异表达。我们计算了伪批量样本,并使用与基于协议的差异表达分析相同的程序和截止值进行了差异表达分析。 , Transcriptomic similarity between HNOCA neural cell types and their primary counterparts in the human developing brain atlas HNOCA神经细胞类型与其在人类发育大脑图谱中的原代对应物之间的转录组相似性
Para_01 为了估计 HNOCA 和人类发育大脑图谱27中神经元之间的转录组相似性,我们首先总结了主要参考文献中每种神经细胞类型的平均表达量,以及 HNOCA 中每个数据集的情况。 对于每个 HNOCA 数据集,仅考虑至少包含 20 个细胞的神经细胞类型。 使用基于卡方检验的方差比测试(在具有伽玛分布的身份链接的广义线性模型上),以跨神经细胞类型的转录计数变异系数作为响应变量,跨神经细胞类型的平均转录计数的倒数作为自变量,识别出主要参考文献中各神经细胞类型中的高变异性基因。 调整后的 Benjamini–Hochberg P 值小于 0.01 的基因被认为是高变异性基因。 然后计算主要图谱中一种神经细胞类型与其在每个 HNOCA 数据集中对应物之间的相似性,方法是计算识别出的高变异性基因上的斯皮尔曼相关系数。 Para_02 为了估计由转录因子共表达定义的核心转录组身份的相似性,高度可变基因被子集化为仅包含转录因子以计算斯皮尔曼相关性。 转录因子列表从 AnimalTFDB v.4.0 数据库中获取。 Para_03 为了在数据集中识别代谢压力细胞,我们使用了 scanpy score_genes 函数(默认参数)来对从 enrichR GO_Biological_Process_2021 数据库获得的‘经典糖酵解’基因集进行评分,评分对象涵盖了来自 HNOCA 和相关文献的所有神经元细胞。 Para_04 为了估计糖酵解评分与全转录组相似性之间的相关性,以及糖酵解评分与核心转录组身份相似性之间相关性的差异显著性,我们生成了100个高变异基因的子集,每个子集的大小与高变异转录因子相同。 基于这些子集计算了转录组相似性,然后将其与糖酵解评分进行了关联。 Heterogeneity of the telencephalic trajectories 端脑轨迹的异质性
Para_01 为了表征 HNOCA 中端脑 NPCs 和神经元的异质性,我们首先将人类新皮层发育图谱中的细胞类型标签(如给定元数据中指示的‘类型’标签)转移到 HNOCA 的端脑 NPCs、中间前体细胞和神经元上,基于转录组相关性。 简而言之,我们将上述获得的每个主要图谱聚类分配给一个细胞类型,作为该聚类中最丰富的细胞类型。 然后,最佳相关的初级聚类的标签被转移到每个查询细胞。 鉴于转移的标签,结合图 1c 中显示的二级细胞类型注释,作为注释标签,scPoli 被应用于 HNOCA 的端脑子集进行数据整合。 Para_02 为了评估不同的整合策略在多大程度上恢复了神经元亚细胞类型的异质性,我们生成了四种不同的聚类标签:(1) 使用原始 scPoli 隐式表示的 Louvain 聚类(分辨率,2);(2) 使用更新后的 scPoli 表示的 Louvain 聚类(分辨率,2);(3) 使用 HNOCA 端脑亚集的 PCA 的 Louvain 聚类(分辨率,2)(基于端脑亚集中 3,000 个高变异基因的缩放表达,flavor = ‘seurat’),以及 (4) 每个样本单独进行的 Louvain 聚类(分辨率,1)(每个样本识别出 3,000 个高变异基因,flavor = ‘seurat’,随后进行数据缩放和 PCA)。 接下来,对于至少包含 500 个背侧端脑神经元的每个样本,计算上述四种聚类标签与作为金标准的转移细胞类型标签之间的调整互信息分数,这些分数是在被注释为第 2 级注释的背侧端脑神经元之间计算的。 Para_03 为了创建一个全面的背侧端脑神经元初级图谱,用于神经类器官和初级组织之间的差异表达分析,我们从四个不同的初级图谱中选取了背侧端脑神经元或新皮质神经元。 对于参考文献28,选择了作者定义的五个聚类(60、57、79、45、65)中的细胞,这些细胞高表达MAP2、DCX和NEUROD6。 对于参考文献29,基于NEUROD6和FOXG1的高表达,选择了具有以下‘clusterv2 - final’标签的细胞:‘Neuron_28’、‘Neuron_34’、‘GW19_2_29NeuronNeuron’、‘Neuron_30’、‘Neuron_66Neuron’、‘GW18_2_42NeuronNeuron’、‘Neuron_33’、‘Neuron_39Neuron’、‘Neuron_35’、‘Neuron_63Neuron’、‘Neuron_9’、‘Neuron_11’、‘Neuron_20’、‘Neuron_22’、‘Neuron_5Neuron’、‘Neuron_21’、‘Neuron_18’、‘Neuron_101Neuron’、‘Neuron_17’、‘Neuron_19’、‘Neuron_16’、‘Neuron_50Neuron’、‘Neuron_12’、‘Neuron_13’、‘Neuron_68Neuron’、‘Neuron_100Neuron’、‘Neuron_25’、‘Neuron_27’、‘Neuron_53Neuron’、‘Neuron_23’、‘Neuron_26’、‘Neuron_24’、‘Neuron_102Neuron’、‘Neuron_72Neuron’、‘Neuron_15’、‘Neuron_29’和‘Neuron_35Neuron’。 对于参考文献27,选择了解剖自背侧端脑且被注释为仅具有VGLUT NTT标签的神经元细胞。 对于参考文献30,选择了被注释为兴奋性神经元的细胞。 如前所述,Wang等人的初级图谱中经过整理的聚类也被子集化为那些具有兴奋性神经元标签的聚类。 所选的图谱中的背侧端脑神经元子集被合并到联合新皮质神经元图谱中。 Para_04 接下来,联合新皮层神经元图谱中的细胞与Wang等人图谱中每个兴奋性神经元簇的平均表达谱相关联。 将最佳相关的簇标签分配给联合新皮层神经元图谱中的每个细胞,从而使图谱中所有细胞的簇标签协调一致。 然后使用scPoli进行标签感知的数据整合。 基于scPoli的潜在表示,对联合新皮层神经元图谱执行Louvain聚类(分辨率,1)。 通过高度可变基因定义的联合新皮层神经元图谱中各簇的平均转录组谱,以最大相关方式将此聚类标签转移到HNOCA中的背侧端脑神经元。 Reference mapping of the neural organoid morphogen screen scRNA-seq data to the human developing brain atlas and HNOCA 神经类器官形态发生筛选的单细胞 RNA 测序数据与人类发育中大脑图谱及 HNOCA 的参考映射
Para_01 我们使用 scArches 将来自神经类器官形态发生筛选的 scRNA-seq 数据映射到人类发育大脑图谱的 scANVI 模型和 HNOCA 的 scPoli 模型。 在这两种情况下,筛选数据的‘数据集’字段被用作批次协变量,这表明属于以下三个类别之一:‘类器官筛选’、‘二次类器官筛选’或‘21孕周胎儿纹状体’。 对于映射到主要参考,我们使用了 scvi-tools 实现的 scArches,不使用细胞类型注释,并且使用权重衰减为 0 的参数训练模型 500 个周期,其余参数保持默认。 对于映射到 HNOCA,我们通过 scPoli 使用 scArches 并训练模型 500 个周期,不进行未标记原型训练。 Retrieval and harmonization of disease-modelling human neural organoid scRNA-seq datasets 疾病建模人脑类器官单细胞RNA测序数据集的检索与整合
Para_01 我们纳入了11个神经类器官的单细胞RNA测序(scRNA-seq)数据集,这些数据集旨在模拟包括小头症、肌萎缩侧索硬化症、阿尔茨海默病、自闭症、脆性X综合征、精神分裂症、神经异位、皮特-霍普金斯综合症、强直性肌营养不良和胶质母细胞瘤在内的10种不同的神经疾病。 对于提供了处理后数据的十个数据集,我们直接从Gene Expression Omnibus或ArrayExpress下载了计数矩阵和元数据。 对于只有FASTQ文件可用的数据集,我们下载了FASTQ文件,并使用Cell Ranger(v.4.0)将读取映射到从Cell Ranger网站获取的人类参考基因组和转录组(GRCh38 v.3.0.0),以进行基因表达量化。 所有数据集使用Python中的anndata连接在一起(join = ‘inner’)。 对于每个数据集,样本根据其疾病状态被分组为‘疾病’或‘对照’,其中‘疾病’代表来自患者细胞系、携带疾病相关等位基因的突变细胞系、CRISPR筛选中携带靶向引导RNA(gRNAs)的细胞和肿瘤衍生的类器官的数据,而‘对照’代表来自健康细胞系、突变校正细胞系和CRISPR筛选中仅携带非靶向gRNAs的细胞的数据。 Projection and label transfer-based annotation of the disease-modelling dataset 基于投影和标签转移的疾病模型数据集标注
Para_01 为了将疾病模型图谱与整合的HNOCA进行比较,我们使用scArches将它投影到HNOCA以及第一孕期的人类大脑单细胞RNA测序图谱。 对于向初级图谱的投影,使用了上述提到的将HNOCA映射到图谱的相同实现方法。 对于向HNOCA的投影,查询模型基于用HNOCA数据预训练的scPoli模型,并以16,384的批量大小进行了最多30轮的微调,其中包含20轮预训练。 基于投影到HNOCA的潜在表示,使用scanpy(默认参数)创建了疾病模型图谱的最近邻图,并使用scanpy(默认参数)创建了UMAP嵌入。 Para_02 接下来,对于 HNOCA 和疾病建模图谱,细胞由 HNOCA-scPoli 和初级 scANVI 模型的连接表示。 然后,通过识别每个疾病建模图谱细胞在 HNOCA 中的 50 个最近邻,如上所述重建了二分 wkNN 图。 基于二分 wkNN,应用多数投票标签转移方法将四级层次细胞类型注释和区域身份转移到疾病建模图谱中。 Reconstruction of matched HNOCA metacells 重建匹配的HNOCA元细胞
错误!!! - 待补充
Para_02 在这里,Ni 表示查询细胞 ci 的所有 HNOCA 最近邻,wik 表示查询细胞 i 和参考细胞 k 之间的边权重,ekj 表示基因 j 在参考细胞 k 中的表达水平。 Para_03 给定匹配的 HNOCA 元细胞转录组谱型,查询细胞与其在 HNOCA 中匹配的细胞状态之间的相似性则通过计算查询细胞转录组谱型与匹配的 HNOCA 元细胞转录组谱型之间的斯皮尔曼相关性来确定。 Re-analysis of GBM-2019 and FXS-2021 datasets GBM-2019和FXS-2021数据集的重新分析
Para_01 为了分析胶质母细胞瘤类器官数据集(GBM-2019),从综合疾病建模图谱中提取了该出版物中的细胞。 使用 Scanpy,默认参数下识别出高变异基因。 然后对高变异基因的对数标准化表达值进行了细胞间的缩放,并使用前 20 个主成分进行截断 PCA 以用于后续分析。 接下来,应用了和谐算法的 Python 实现版本 Harmonypy,以整合来自不同样本的细胞。 基于和谐整合的嵌入,重建了邻域图。 UMAP 嵌入和 Louvain 聚类(分辨率,0.5)是根据最近邻图创建的。 在识别出的 12 个聚类中,选择 AQP4 表达最高的两个聚类,即聚类-7 和聚类-0,进行后续差异表达分析。 Para_02 为了分析 FXS 数据集(FXS-2021),从出版物中提取了细胞子集,这些细胞来自整合的疾病建模图谱。 应用了与 GBM-2019 数据集相同的高变基因识别、数据缩放和主成分分析(PCA)程序。 接下来,基于前 20 个主成分直接创建了最近邻图。 然后,基于重建的最近邻图创建了 UMAP 嵌入和 Louvain 聚类(分辨率,1)。 在 30 个聚类中,选择了表达 EMX1 和 FOXG1 的 cluster-17 和 cluster-23,根据从 HNOCA 转移的标签,这两个聚类主要被预测为背侧端脑神经前体细胞和神经元,用于后续的差异表达分析。 F-test-based DE analysis for paired transcriptome 基于F检验的配对转录组差异表达分析
Para_01 为了比较两组配对细胞的表达水平,首先根据对数标准化表达值计算每个细胞对的每基因表达差异。 接下来,对于每个要测试差异表达的基因,其在计算出的每个细胞对表达差异上的方差(σ²)与表达差异平方和(对于基因i为di)除以细胞对数量进行比较: Data availability Para_01 所有精心策划的个体 HNOCA 数据集都可以通过 sfaira Python 工具轻松访问。 整合的 HNOCA 数据可在 Zenodo (https://doi.org/10.5281/zenodo.11203684) 和 CellxGene Discover Census (https://cellxgene.cziscience.com/collections/de379e5f-52d0-498c-9801-0f850823c847) 获取。 扩展的 HNOCA 社区版图谱也可通过 CellxGene Discover Census (上述相同网址) 获得。 两个版本的 HNOCA 均可通过 ArchMap 网络界面 (https://www.archmap.bio/) 进行参考映射。 HNOCA-tools 包提供了用于注释、参考映射和核心下游分析步骤的 Python 接口,并可在 https://github.com/devsystemslab/HNOCA-tools 获取。 有关可用工具的信息和 HNOCA-tools 的文档可在 https://devsystemslab.github.io/HNOCA-tools 查看。 用于重现分析的 Jupyter 笔记本和脚本可在 https://github.com/theislab/neural_organoid_atlas 获取。