Basic Information
- 英文标题:Spatially resolved mapping of cells associated with human complex traits
- 中文标题:与人类复杂性状相关的细胞空间分辨图谱
- 发表日期:19 March 2025
- 文章类型:Article
- 所属期刊:Nature
- 文章作者:Liyang Song | Jian Yang
- 文章链接:https://www.nature.com/articles/s41586-025-08757-x
Abstract
Para_01
- 描绘与疾病相关的细胞的空间分布对于理解疾病病理学至关重要。
- 在此,我们提出了针对复杂性状的基因信息驱动的细胞空间定位方法(gsMap),该方法整合了空间转录组学数据和全基因组关联研究的汇总统计数据,以空间解析的方式将细胞映射到人类复杂性状,包括疾病。
- 利用覆盖25个器官的胚胎空间转录组学数据集,我们通过模拟和验证各器官中已知的与性状相关联的细胞或区域对gsMap进行了基准测试。
- 将gsMap应用于脑部空间转录组学数据后,我们发现与精神分裂症相关的谷氨酸能神经元的空间分布更接近于认知性状,而非情绪性状(如抑郁)。
- 与精神分裂症相关的谷氨酸能神经元分布在背侧海马附近,并且钙信号和调控基因的表达上调。
- 而与抑郁症相关的谷氨酸能神经元分布在深层内侧前额叶皮层附近,其神经可塑性和精神药物靶基因的表达上调。
- 我们的研究提供了一种用于性状相关细胞的空间解析映射的方法,并展示了通过这些图谱获得生物学见解(例如性状相关细胞及其特征基因的空间分布)的价值。
Main
Para_01
- 组织内细胞的组成和空间组织对其功能至关重要,也可以作为其健康状态的指标。
- 空间转录组学(ST)的进步使得在细胞原始空间位置上对基因表达水平进行分析成为可能,这为研究细胞空间组织及揭示相关生物学机制提供了有前景的方向。
- 近年来,越来越多的研究利用ST技术探索不同组织中的细胞空间组织。
- 然而,在识别与复杂性状或疾病最相关的细胞并绘制它们在组织内的空间分布方面,仍然存在显著的知识空白。
Para_02
- 为了识别与性状相关的细胞或细胞类型,先前的研究提出了基于遗传学信息的策略,这些策略整合了复杂性状(包括疾病)的全基因组关联研究 (GWAS) 数据与单细胞 RNA 测序 (scRNA-seq) 数据。
- 尽管这些方法可以精确定位与性状相关的细胞,但由于 scRNA-seq 数据中缺乏细胞空间位置信息,它们在绘制这些细胞的空间分布时面临挑战。
- 从原则上讲,这些基于 scRNA-seq 的方法可以应用于空间转录组学 (ST) 数据。
- 然而,由于缺乏对细胞空间坐标的建模以及 ST 数据中高水平的技术噪声,这些方法在空间感知的性状相关细胞映射方面能力有限。
- 虽然先前的研究将空间域与复杂性状相关联,但这些分析并未达到单细胞分辨率,限制了其描绘性状相关细胞空间分布的能力。
- 因此,需要开发新方法,能够将 ST 数据整合到 GWAS 中,以实现性状相关细胞的精细空间解析映射。
Para_03
- 在本研究中,我们介绍了 gsMap,这是一种整合高分辨率空间转录组数据和全基因组关联研究摘要统计数据的方法,用于对与性状相关的细胞进行空间定位映射。
- 利用覆盖 25 个器官的胚胎空间转录组数据集,我们通过模拟的全基因组关联研究数据评估了 gsMap 的特异性,并通过重现不同器官中的细胞与一系列复杂性状之间已知的关联验证了该方法的敏感性。
- 将 gsMap 应用于脑部空间转录组数据集后,我们生成了涵盖 30 种与人类脑部相关的复杂性状的大规模性状–脑细胞关联图谱。
Overview of gsMap
Para_01
- gsMap 的基本概念涉及评估位于或靠近在空间转录组(ST)数据中某一点高度表达的基因的遗传变异(主要是单核苷酸多态性,SNPs),是否在与某一特定性状相关的遗传关联中富集(方法和补充说明,第1-3节)。
- 每个 ST 切片需要包括全转录组的基因表达谱和单个点的空间坐标。
- 在这里,"点"指的是高分辨率 ST 平台中的单个细胞(例如 Stereo-seq 细胞模式)或传统 ST 平台中的细胞簇(如 10X Visium)。
- gsMap 包含三个步骤。
- 首先,gsMap 利用点的同质性来解决 ST 数据中的稀疏性和技术噪声问题。
- gsMap 使用图神经网络(GNN)根据基因表达模式和空间位置识别每个焦点周围的同质点。
- 通过聚合这些同质点的信息,计算每个点的基因特异性评分(GSS),以表示某个基因在该点的表达水平的相对排名(图1a)。
- 其次,gsMap 将每个点的 GSS 分配给每个基因转录区域上游 50 kb 和下游 50 kb 窗口内的 SNPs,并结合使用表观基因组数据建立的 SNP 到基因映射,生成每个点的独特 SNP GSS 注释集。
- 将每个点视为一个 SNP 注释集,gsMap 使用分层连锁不平衡得分回归(S-LDSC)评估具有较高 GSS 的 SNP 是否不成比例地解释了该性状更大比例的遗传力,条件是基于基线注释(图1b)。
- 富集 P 值用于衡量某一点与性状之间关联的统计显著性。
- 最后,为了量化特定空间区域与性状关联的显著性,gsMap 使用柯西组合检验来聚合该空间区域内单个点的 P 值(图1c)。
- 本质上,gsMap 可被视为一种在细胞分辨率下将复杂性状进行遗传信息映射到 ST 数据的方法。
Fig. 1: Schematics of the gsMap method.
- 图片说明
◉ gsMap 首先使用图神经网络 (GNN) 学习嵌入表示,这些嵌入整合了基因表达水平、空间坐标,并可选地整合细胞类型注释先验信息。◉ 随后,gsMap 根据嵌入中的余弦相似性为每个位点识别同质位点,以形成一个微区域。◉ 每个位点依次被视为焦点位点,并通过将基因在微区域内的平均排名除以其在整个空间转录组 (ST) 切片中的平均排名,计算每个基因在每个焦点位点的特异性得分。◉ D,位点的微区域;F,基因表达特异性;G,位点空间图;R,排名;X,基因表达矩阵;Z,嵌入表示。◉ 然后,基于每个位点的 GSS(基因空间特异性得分)与其到转录起始位点 (TSS) 的距离以及 SNP 到基因的链接图,将这些 GSS 映射到单核苷酸多态性 (SNP),从而为每个位点生成一组独特的 SNP 注释。◉ 对于每个位点的 SNP GSS 注释,gsMap 使用分层连锁不平衡回归 (S-LDSC) 来评估具有更高 GSS 的 SNP 是否对目标性状的遗传力有显著富集。◉ LD,连锁不平衡。◉ 为了量化某个空间区域与性状关联的显著性,gsMap 使用柯西组合检验来聚合该空间区域内位点的 P 值。
Validation of gsMap
Para_01
- 作为原理验证,我们试图通过验证不同组织和性状之间已知的关联来验证 gsMap。
- 对于这一验证分析,我们使用了从胚胎晚期阶段收集的空间转录组数据,因为它涵盖了广泛的组织谱系。
- 由于缺乏人类在胚胎晚期阶段的空间转录组数据,我们使用了 Chen 等人提供的小鼠数据集,其中包括胚胎第 16.5 天(E16.5)发育阶段的 25 种不同器官(图 2a 和方法部分)。
- 该策略假设小鼠的基因表达谱与人类相似,这一假设得到了先前证据的广泛支持,表明小鼠和人类之间超过 80% 的基因在组织水平表达上具有相关性,并且在细胞类型水平上有超过 75% 的相关性(方法部分和补充图 1)。
- Chen 等人的小鼠空间转录组数据集是通过 Stereo-seq 技术生成的,包含‘bin50’分辨率的数据,其中每个点代表几个细胞,并且还包含单细胞分辨率的数据。
- 为了确保稳健性,我们首先使用了 bin50 分辨率的数据。
- 我们在分析中纳入了 110 种复杂性状的公开可用 GWAS 汇总统计数据(平均样本量 n = 385,000;补充表 1 和补充说明第 4 节)。
- 我们的 gsMap 成功重现了不同器官与性状之间的已知关联。
- 例如,智商(IQ)被映射到大脑中的点,平均红细胞血红蛋白浓度(MCHC)与肝脏中的点相关联,而身高则被映射到多种组织,软骨及其原基显示出最高的相关性(图 2b)。
- 通过柯西组合检验获得的汇总组织水平关联 P 值证实了这些观察结果(图 2c 和补充表 2)。
- 敏感性分析表明,这些结果并非由每个点上具有最高 GSS 的少数基因驱动(补充图 2)。
- 我们还将 gsMap 应用于 Carnegie 第 8 阶段(CS8)完整人类胚胎的空间转录组数据。
- 尽管在此发育阶段器官尚未完全形成,gsMap 仍成功将性状映射到其各自的器官前体。
- 例如,智商(IQ)被映射到外胚层(P = 1.5 × 10⁻⁶),即大脑的前体,而 MCHC 被映射到内胚层(P = 1.5 × 10⁻⁹),即肝脏的前体(扩展数据图 1 和补充表 3)。
Fig. 2: Mapping human complex traits to mouse embryo.
- 图片说明
◉ a, 小鼠 E16.5 胚胎 ST 数据,分辨率为 bin50,点根据组织类型着色。GI,胃肠道。◉ b, 使用小鼠 E16.5 胚胎 ST 数据的 gsMap 结果,针对 IQ、MCHC 和身高的分析,颜色表示关联的显著性。基因组膨胀因子 (λGC) 是 GWAS 统计功效的指标。◉ c, 组织-性状关联,行代表组织类型,列代表性状。颜色表示关联的显著性。性状缩写见补充表 1。◉ d, 小鼠 E16.5 胚胎 ST 数据,分辨率为单细胞,点根据细胞类型着色。◉ e, 左图,小鼠 E16.5 胚胎中上皮细胞的分布,白色矩形显示六个空间区域。棕色点代表上皮细胞。右图,MPB 的 gsMap 结果,点根据关联的显著性着色。◉ f, 六个空间区域中 MPB 的 gsMap 结果。◉ g, 面部 (L2) 区域中上皮细胞的密度。从蓝色到黄色的渐变颜色表示细胞密度从低到高。黑点代表上皮细胞。◉ h, 上皮细胞密度与面部区域上皮细胞与 MPB 关联显著性之间的相关性。上皮细胞密度估计为 100 个网格内上皮细胞的平滑相对计数,并使用插值法估计每个网格中上皮细胞-MPB 关联的显著性。x 轴表示上皮细胞密度,y 轴表示上皮细胞-MPB 关联的显著性。数据点代表网格,根据关联的显著性着色。黑线表示回归线;P 值来自双侧 t 检验 (n = 100)。
Para_02
- 为了验证所识别的性状与位点之间关联的可靠性,我们基于来自英国生物样本库的10万名无亲缘关系的欧洲血统个体的真实基因型数据,模拟了四种不同水平的遗传率和多基因性的零假设情景(方法)。
- 在这些情景中,因果变异要么随机分布,要么集中在特定的基因组区域,但并未富集于具有空间依赖性表达的基因中。
- 在所有零假设情景下,gsMap 的错误发现率(FDR)得到了良好控制,因为它未在 FDR < 0.05 的情况下识别出任何显示显著关联的位点(扩展数据图 2 和补充图 3 和 4)。
Para_03
- 为了展示 gsMap 在将细胞空间定位到某一性状上的能力,我们将其应用于小鼠胚胎单细胞分辨率的空间转录组数据(图 2d 和方法部分)。
- 以男性型秃发(MPB)为例,我们发现 MPB 与面部区域的上皮细胞之间存在关联(图 2e)。
- 值得注意的是,与 MPB 相关的上皮细胞并非随机分布,而是倾向于在空间上聚集在一起(图 2f)。
- 细胞与 MPB 关联的显著性水平与上皮细胞密度之间存在强相关性(r = 0.78,P = 2.2 × 10⁻¹³,图 2g、h)。
- 上皮细胞紧密的空间排列表明它们可能形成毛囊,这一结论得到了毛囊相关标记基因 KRT15、KRT5 和 KRT17 表达谱的支持(补充图 5)。
Para_04
- 如上所述,将 GWAS 与 scRNA-seq 数据整合的现有方法(例如单细胞疾病相关性评分(scDRS))可以重新用于 ST 数据,即使它们不包含空间坐标信息。
- 我们使用模拟和真实的 GWAS 数据集(补充图 6 和 7 以及补充说明第 6 节)比较了 gsMap 和 scDRS。
Mapping human traits to mouse brain
Para_01
- 我们首先将 gsMap 应用于将人类复杂性状映射到大脑区域,这主要得益于大量脑空间转录组数据的可用性,尽管大多数并非来自人类样本。
- 为了确保统计功效,我们在分析中纳入了 30 种不同的脑相关性状,这些性状具有大规模全基因组关联研究样本(平均样本量 n = 315,000),涵盖了认知、情绪和行为。
- 由于缺乏人类全脑空间转录组数据,我们开始整合全基因组关联研究的汇总统计数据与小鼠脑空间转录组数据,并主要集中于进化上保守的大脑区域,以确保结果的广泛外推性。
- 我们处理了 Chen 等人提供的成年小鼠半脑冠状切片空间转录组数据,该数据包括来自 14 个脑区和 13 种细胞类型的 50,140 个细胞(图 3a)。
- 我们首先通过使用柯西组合检验汇总每个脑区单个细胞的 P 值来评估整个脑区与某一性状的关联(图 3b)。
- 我们观察到,皮层在大多数性状中显示出最高的相关性,例如智商(IQ)(P = 9.8 × 10⁻¹⁵)、精神分裂症(SCZ)(P = 1.7 × 10⁻¹⁷)和抑郁症(P = 1.0 × 10⁻¹⁸)。
- 次高相关的脑区因性状而异:智商被映射到海马 CA1 区(P = 7.2 × 10⁻¹³),精神分裂症映射到海马 CA1 区(P = 4.4 × 10⁻¹⁴),而抑郁症则映射到中脑(P = 1.1 × 10⁻¹⁵)。
Fig. 3: Mapping human complex traits to adult mouse brain.
- 图片说明
◉ a,成年小鼠脑空间转录组数据,点根据脑区(左)或细胞类型(右)着色。◉ b,脑区与性状关联,颜色表示关联的显著性。行代表脑区,列代表性状。右侧条形图显示每个脑区中谷氨酸神经元的比例,顶部条形图表示λGC,作为GWAS统计功效的指标。◉ 脑区和细胞类型的缩写列于补充说明第5节。◉ c,不同脑区中的细胞类型比例。◉ d,脑区与精神分裂症(SCZ)关联显著性(x轴)与其谷氨酸神经元比例(y轴)的相关性。灰色线为回归线;P值来自双侧t检验(n = 14)。◉ e,性状关联显著性与不同脑区中某细胞类型比例的相关性(n = 30)。x轴显示不同细胞类型,y轴显示皮尔逊相关系数。点代表性状,并按细胞类型着色。◉ 每个箱线图中,中心线表示中位数,凹槽代表95%置信区间,箱子表示四分位距,须延伸至1.5倍四分位距。◉ AD,阿尔茨海默病;PD,帕金森病。◉ f,中脑和CA区域中智商(IQ)、精神分裂症(SCZ)和抑郁症的gsMap结果。每个点代表一个单独的谷氨酸神经元,颜色表示其与性状关联的显著性。◉ g,大脑皮层(左)、CA(中)和中脑(右)的优势比(ORs)。误差线表示95%置信区间,中心表示OR值。P值来自双侧卡方检验(n = 19,598)。柱状图确保小范围误差线的可见性(未显示单个数据点)。◉ CoPe,认知表现。
Para_02
- 由于大脑区域包含多种细胞类型,我们通过评估特定细胞类型的较高比例是否与大脑区域-性状关联的显著性增加相关,来探索它们对这种关联的贡献。
- 我们发现,在所有细胞类型中,平均而言,谷氨酸能神经元(Glu-neurons)在大脑区域-性状关联中的贡献最强(中位数 r = 0.73)。
- 例如,大脑区域与精神分裂症(SCZ)关联显著性与其谷氨酸能神经元比例之间的相关性为 0.80(P = 6.3 × 10⁻⁴)(图 3c-e)。
- 使用来自小鼠 E16.5 胚胎的脑矢状切面空间转录组数据(ST 数据),这一关系得到了重复验证(扩展数据图 3)。
- 然而,有两个异常值,即帕金森病和阿尔茨海默病。
- 我们观察到,多巴胺能神经元(DA-neurons)比例较高的大脑区域与帕金森病的关联更强(r = 0.81,P = 4.9 × 10⁻⁴;图 3e)。
- 具体而言,已知多巴胺能神经元丰度最高的黑质/腹侧被盖区(SN/VTA),显示出与帕金森病的最高相关性(P = 5.2 × 10⁻⁵;扩展数据图 4)。
- 这一发现与先前的神经学研究一致,表明帕金森病患者的 SN/VTA 区域存在显著损伤。
Para_03
- 在注意到谷氨酸能神经元在大多数分析性状中扮演重要角色后,我们进一步探索了与这些性状相关的谷氨酸能神经元是否在分布模式上因性状不同而有所差异。
- 考虑到直接比较不同性状间的 gsMap 结果可能会受到全基因组关联研究统计功效差异的干扰,我们使用了优势比(OR),其计算方法是特定区域中与性状相关和不相关的谷氨酸能神经元比例除以其他所有区域的比例(见方法部分)。
- 我们的分析集中于精神分裂症(SCZ)、认知性状(如智商)和情绪性状(如抑郁),因为谷氨酸能神经元对这些性状与脑区关联的贡献显著(图3e),且其全基因组关联数据具有足够的统计效力(λGC > 1.3)。
- 我们发现,除了大脑皮层外,海马体CA区域的谷氨酸能神经元与认知性状密切相关(例如,智商的优势比为1.5,P值为4.3 × 10^-25)。
- 而与情绪性状相关的谷氨酸能神经元则倾向于分布在中脑区域(例如,抑郁的优势比为6.7,P值为7.9 × 10^-221;图3f、g)。
- 值得注意的是,尽管精神分裂症包含认知和情绪症状,但我们发现与其相关的谷氨酸能神经元表现出与认知性状相似的空间分布模式,并且在海马体CA区域特别富集(优势比为4.5,P值为7.6 × 10^-158;图3g)。
Para_04
- 我们的研究结果共同揭示了以下几点:
- 在大脑区域中,皮层与大多数性状表现出最强的关联,但帕金森病仅与黑质/腹侧被盖区显著相关;
- 谷氨酸能神经元对性状与脑区之间的关联贡献显著;
- 分布在不同脑区的谷氨酸能神经元显示出与性状的不同关联。
- 分布在海马CA区的谷氨酸能神经元与认知性状和精神分裂症表现出强关联,
- 而分布在中脑的谷氨酸能神经元则与情绪障碍(如抑郁症)表现出强关联。
Trait-associated neurons in mouse hippocampus
Para_01
- 在展示了与性状相关的谷氨酸能神经元(Glu-neurons)比例在不同脑区中存在差异后,我们接下来试图探讨这些神经元是否在单个脑区内表现出特定的空间分布模式。
- 我们首先聚焦于海马体的阿蒙角1区(CA1),因为CA1区域的空间结构相对简单,并且先前的研究表明,CA1中的锥体神经元的电生理特性在其空间位置上表现出变化。
- CA1区域内的细胞沿三个空间轴密集排列:近端-远端(P-D)轴、背侧-腹侧(D-V)轴和深层-浅层(D-S)轴(图4a)。
- 由于本研究使用的脑空间转录组数据对应于冠状切面,因此无法追踪Glu-neurons沿P-D轴的空间分布。
- 因此,我们的研究集中在Glu-neurons沿D-V轴和D-S轴的空间分布上。
Fig. 4: Trait-associated Glu-neurons in mouse hippocampus.
- 图片说明
◉ a,小鼠海马体三个空间轴的示意图。◉ b,SCZ在CA1区域的gsMap结果,颜色表示关联的显著性。◉ c,Glu神经元与SCZ关联显著性与其沿CA1背-腹侧轴空间位置的相关性,细胞根据关联显著性着色。黑色回归线的P值来自双侧t检验(n = 2,824)。◉ d,相关分析的z统计量。x轴显示Glu神经元基因表达水平与其沿背-腹侧轴空间位置相关性的z统计量,y轴显示其与SCZ关联的显著性。FDR < 0.01的基因以橙色突出显示。◉ e,基因表达水平与Glu神经元沿背-腹侧轴空间位置相关(n = 241)或与SCZ关联显著性(n = 100)的基因的GO术语富集分析。◉ f,Glu神经元与性状关联显著性及其沿CA1背-腹侧轴空间位置相关性的z统计量。虚线表示FDR = 0.01的阈值。◉ g,顶部为在CA1深层、浅层和中脑高表达基因的重叠情况。底部为在中脑和CA1浅层高表达基因与在中脑和CA1深层高表达基因对比的比值比(OR)。误差条表示95%置信区间,中心线表示OR(n = 259)。◉ h,在中脑和CA1浅层均高表达基因(n = 55)的GO术语富集分析。◉ i,CA1 Glu神经元沿背-腹侧轴的位置(x轴)与基因表达与中脑Glu神经元余弦相似性(y轴)的相关性,细胞根据其沿背-腹侧轴位置着色。黑色回归线的P值来自双侧t检验(n = 2,824)。e、g和h中的P值来自双侧卡方检验。
Para_02
- 首先,沿着背腹轴(D–V轴),我们观察到靠近CA1背侧的谷氨酸能神经元(Glu-neurons)与特征之间的关联更强(补充图8),这一趋势在精神分裂症(SCZ)中尤为明显(r = −0.51,P = 3.5 × 10−190;图4b,c)。
- 我们通过使用小鼠脑多重错误鲁棒荧光原位杂交(MERFISH)空间转录组数据集复制了这一沿D–V轴的关联强度下降现象,该数据集包括来自多个个体的冠状面和矢状面切片(扩展数据图5)。
- 为了研究潜在机制,我们进行了两项表达关联分析,以识别其表达水平与Glu-neurons沿D–V轴位置相关的基因,以及其表达水平与Glu-neurons对SCZ相关性相关的基因(方法)。
- 这两项分析的z统计量显示出强相关性(r = 0.78,P < 5 × 10−324),其中HPCA(编码钙结合蛋白)、PPP3CA(编码钙调磷酸酶亚基)和ATP2B1(编码钙转运蛋白)在两项分析中均被列为优先排名的前三个基因(图4d)。
- 随后的基因本体(GO)术语富集分析显示,这些鉴定出的基因在与钙和离子转运相关的通路中显著富集(图4e)。
Para_03
- D-S轴的空间尺度较窄,但我们仍观察到沿该轴的谷氨酸能神经元(Glu-neurons)存在异质性。
- 利用gsMap,我们发现靠近CA1浅层侧的谷氨酸能神经元与抑郁的关联更强(r = 0.18,P = 2.7 × 10⁻²⁰;图4f)。
- 这一发现与最近的一项研究一致,该研究表明操纵CA1浅层神经元而非深层神经元可以改善小鼠的抑郁样行为。
- 受上述结果的启发,即中脑中的谷氨酸能神经元与抑郁高度相关(图3g),我们进一步研究了CA1浅层侧和中脑中谷氨酸能神经元之间的关联。
- 我们发现,在中脑高表达(FDR < 0.05)的基因与在CA1浅层侧高表达的基因相比,与CA1深层侧高表达基因的重叠更多(55对3个基因,OR = 14.75,P = 2.25 × 10⁻⁹;图4g)。
- GO富集分析表明,在中脑和CA1浅层侧均高表达的基因在轴突发育和神经元包鞘相关的通路中显著富集(图4h)。
- CA1和中脑中谷氨酸能神经元的转录组相似性随着CA1 D-S轴增加而增强(r = 0.41,P = 2.6 × 10⁻¹¹⁷;图4i)。
- 这些结果表明,中脑中的谷氨酸能神经元与CA1浅层侧的谷氨酸能神经元在功能上存在重叠。
Para_04
- 总之,我们的研究结果揭示了海马体CA1区域中与性状相关的谷氨酸能神经元(Glu-neurons)在空间上的分布模式。
- 靠近CA1背侧的谷氨酸能神经元显示出与精神分裂症(SCZ)更强的关联,并且钙信号传导和调节相关基因的表达增加。
- 靠近CA1浅层侧的谷氨酸能神经元则表现出与抑郁症更强的关联,并且轴突发育相关基因的表达增加。
Mapping human traits to macaque cortex
Para_01
- 上述针对成年小鼠空间转录组数据的分析集中于海马体和中脑,这些区域在不同哺乳动物物种中从空间结构和生物学功能上来看是进化保守的。
- 为了将我们的分析扩展到进化上更高级的大脑区域——大脑皮层,我们遇到了挑战,这是由于人和小鼠皮层之间存在功能和结构上的差异。
- 目前可用的人类皮层空间转录组数据集仅覆盖了背外侧前额叶皮层的一小部分,因此其覆盖范围和分辨率有限。
- 因此,我们使用了陈等人提供的猕猴大脑皮层空间转录组数据集,该数据集更接近人类皮层,并提供了来自三只成年猕猴的143个皮层区域的单细胞分辨率(补充表4和表5)。
Para_02
- 猕猴ST数据中连续皮层切片的可用性为验证gsMap的稳健性提供了额外的机会。
- gsMap结果显示相邻ST切片之间具有显著的一致性(例如,SCZ的中位r = 0.92)。
- 接下来,我们比较了来自不同生物学重复实验的空间匹配ST切片之间的SCZ gsMap结果。
- 这些比较也显示出高度一致性,使用猕猴皮层ST数据集时中位r值为0.84,使用小鼠胚胎ST数据集时中位r值为0.83。
- 此外,尽管缺乏单细胞分辨率的人类全皮层数据,我们检查了人类DLPFC和猕猴前额叶皮层(PFC)之间gsMap结果的一致性。
- 由于人类DLPFC采样位置的不确定性,我们选择使用猕猴PFC,这使得从猕猴中准确选择匹配的皮层区域变得困难。
- 基于猕猴PFC ST数据(Stereo-seq)的gsMap结果与基于人类DLPFC ST数据(10X Visium)的结果高度相关(r = 0.51,P = 8.6 × 10−8),尽管存在包括ST平台差异和皮层定义差异在内的不一致。
- 鉴于观察到的一致性,我们将30种与人类大脑相关的特征映射到整个猕猴左大脑皮层的ST切片上,达到了空间分辨的单细胞分辨率。
- 这些映射结果可在https://yanglab.westlake.edu.cn/gsmap获取。
Para_03
- 为了系统分析脑相关性状与皮层细胞之间的空间关联图谱,我们首先估计了成对性状之间关于它们与整个大脑皮层细胞关联的相关性,以下称为细胞相关性相关(CRC)。
- 我们观察到精神分裂症(SCZ)与其他认知性状(例如认知表现)聚为一类,而情绪性状则形成另一类(扩展数据图6c)。
- 在调整全基因组关联研究(GWAS)统计功效后,这一结果仍然一致,并且与我们使用小鼠脑空间转录组(ST)数据的研究结果一致(图3和图4以及补充图12和13),表明认知性和情绪性状在相关脑区上存在差异。
- 值得注意的是,两种性状之间的CRC可能与其遗传相关性有很大差异(扩展数据图6d)。
- 例如,尽管糖化血红蛋白(HbA1c)与高密度脂蛋白(HDL)之间的遗传相关性为负(rG = −0.18,P = 1.6 × 10−11),但它们的CRC为正(CRC = 0.79,P < 5.0 × 10−324),这是由于它们与分布在肝脏中的细胞有共同的关联(补充图14)。
- 精神分裂症(SCZ)是一种复杂的精神性疾病,涉及认知、情绪和行为症状的组合,通常不被单纯归类为认知或情绪障碍。
- 然而,在我们的研究中,我们发现SCZ与认知性状(例如,认知表现与SCZ之间的CRC = 0.72)在皮层细胞中的CRC比与情绪性状(例如,抑郁与SCZ之间的CRC = 0.55)更强(扩展数据图6e)。
- 这一结果表明,SCZ与认知表现之间广泛共享细胞和皮层区域,提示与认知相关的脑区损伤可能是SCZ患者的主要病理变化。
Trait-associated neurons in macaque PFC
Para_01
- 利用上述生成的空间解析皮层细胞-性状关联图,我们研究了成年猕猴前额叶皮层(PFC)中与性状相关的神经元的空间分布,考虑到先前的研究已表明前额叶皮层在情感和认知中具有重要作用。
- 其他皮质叶中的神经元同样至关重要;然而,由于篇幅限制,我们重点关注前额叶皮层,并在我们的交互式网站上展示了其他皮质叶的结果。
Fig. 5: Trait-associated Glu-neurons in macaque PFC.
- 图片说明
◉ a,猕猴前额叶皮层(PFC)空间转录组数据,细胞按皮层区域着色。EBZ表示切片位置。◉ b,谷氨酸能神经元的比值(ORs),底部注释表示空间转录组切片。◉ c,抑郁症的gsMap结果,细胞按关联显著性着色。◉ d,14r区域中谷氨酸能神经元亚型的比例。◉ e,14r及其他区域中谷氨酸能神经元亚型的gsMap结果。每个图中显示双侧Wilcoxon秩和检验的P值(n = 88,857)。中央线表示中位数,凹槽代表95%置信区间,箱子表示四分位距,触须延伸至1.5倍四分位距。◉ f,L2(L2/3)谷氨酸能神经元沿L-M轴的位置(x轴)与其与抑郁症关联显著性(y轴)之间的相关性。点表示单个谷氨酸能神经元,按亚型着色。黑色线为回归线,阴影区域表示95%置信区间。◉ g,谷氨酸能神经元特性与其沿L-M轴位置之间关联显著性的相关性z分数。行对应谷氨酸能神经元亚型,列表示特性。◉ h,14r外侧(左,n = 1,344)和内侧(右,n = 1,340)谷氨酸能神经元的比值(ORs)。柱状图确保小范围误差条可见(未显示单独数据点)。虚线表示14r的OR值。◉ i,与小鼠中脑或CA1(S)相比,猕猴PFC 14r内侧区域(n = 143)高表达基因的富集倍数,相对于其他PFC区域(n = 601)。◉ j,14r内侧区域(n = 143)高表达基因的GO术语富集分析。◉ k,在精神药物靶点中,14r内侧区域(n = 143)高表达基因的富集倍数,相对于所有其他检测到的基因(左,n = 13,313)或其他PFC区域高表达基因(右,n = 601)。h、i、k中的误差条表示95%置信区间,中心表示OR值。P值来自双侧卡方检验。
Para_02
- 我们首先探讨了与性状相关的谷氨酸能神经元(Glu-neurons)在前额叶皮层(PFC)不同区域的分布情况。
- 我们发现,位于PFC 14r区域(直回)的谷氨酸能神经元与抑郁症表现出强烈的关联(OR = 5.3,P = 5.2 × 10^-123;图5b、c),并且当使用相邻PFC ST切片的数据时,这些结果仍然保持一致(补充图15)。
- 这种显著的关联在大多数谷氨酸能神经元亚型(层)中都很明显,但罕见的被标注为L3/4/5的谷氨酸能神经元除外(图5d、e和补充图16)。
- 此外,我们发现PFC 14r区域中性状与谷氨酸能神经元之间的关联与其局部空间分布相关。
- 沿着14r区域的外侧-内侧(L-M)轴,我们观察到谷氨酸能神经元与情绪障碍的相关性逐渐增加,特别是在浅层皮质层(如L2和L3;图5f、g)。
- 将14r区域分为外侧和内侧两部分后,抑郁症在外侧(OR = 4.1,P = 2.9 × 10^-78)和内侧(OR = 6.7,P = 3.6 × 10^-169)均表现出与谷氨酸能神经元的显著关联,其中内侧更为显著(图5h)。
- 为了研究潜在的生物学机制,我们将PFC 14r区域内侧谷氨酸能神经元的基因表达水平与其他PFC区域进行比较,并在14r内侧鉴定出143个高表达基因(FDR < 0.05;补充图17)。
- 基因本体(GO)富集分析显示,这些基因富集于突触组织和细胞连接组装相关的通路中(图5j)。
- 回顾上述发现,即分布在中脑和CA1浅层的谷氨酸能神经元与抑郁症高度相关(图3和图4),我们还评估了这些基因与小鼠脑ST数据中中脑或CA1浅层高表达基因的重叠情况。
- 与高表达于其他PFC区域的基因相比,这种重叠显著更高(OR = 6.0,P = 4.4 × 10^-11;图5i)。
- 此外,高表达于中脑和CA1浅层的基因富集于轴突发生和细胞黏附相关的通路(图4i),这些通路与在PFC 14r内侧鉴定出的突触组织和细胞连接组装通路密切相关。
- 这些通路共同促进了神经可塑性,这是神经元适应环境刺激的基本机制,而这一机制可能在抑郁症患者中受到破坏。
Para_03
- 接下来,为了展示由 gsMap 优先排序的 PFC 14r 区域的临床价值,我们进行了富集分析,以探讨在 PFC 14r 区域高表达的基因是否在已批准或上市的精神疾病药物靶基因中富集。
- 我们从 Drug Repurposing Hub 和 DrugBank 数据库中收集了 417 个精神疾病药物靶基因。
- 与在其他 PFC 区域高表达的基因(OR = 5.2,P = 1.6 × 10⁻¹⁰)或猕猴 ST 数据中捕获的所有其他基因(OR = 5.1,P = 3.7 × 10⁻⁴;补充图 18)相比,在 PFC 14r 区域高表达的基因在精神疾病药物靶基因中显示出显著富集。
- 对于在 PFC 14r 区域内侧高表达的基因,这种富集更为显著(OR = 16.0,P = 2.2 × 10⁻¹⁰;图 5k)。
- 这些基因还与通过常用复杂性状基因优先排序方法鉴定的基因存在显著重叠(补充图 19 和方法部分)。
- 一致地,我们观察到分布在 PFC 14r 内侧的细胞相较于分布在其他 PFC 区域的细胞表现出最高的药物模块评分(扩展数据图 9)。
- 我们进一步探索了那些靶基因在 PFC 14r 内侧高表达基因中富集的药物,并将其与所有其他基因进行比较(补充说明,第 7 节和补充表 6)。
Para_04
- 我们的研究结果共同揭示了前额叶皮层(PFC)中谷氨酸能神经元(Glu-neurons)与特定性状之间存在空间模式化的关联。
- 我们发现,靠近PFC 14r区(直回)内侧的谷氨酸能神经元与抑郁症密切相关,并且在该区域高表达的基因富集于神经可塑性相关通路以及精神疾病药物的靶点。
Discussion
Para_01
- 在此我们介绍了 gsMap,它整合了细胞基因表达谱、细胞空间坐标、SNP 到基因的链接图以及 GWAS 汇总数据,以将细胞在空间上映射到人类复杂性状。
- gsMap 方法已实现为 Python 包,并可免费获取于 https://github.com/JianYang-Lab/gsMap。
- 通过使用真实的空间转录组学(ST)和真实的及模拟的 GWAS 数据集进行广泛的基准分析,我们证明了 gsMap 在空间感知下识别性状与细胞之间的关联时具有准确性、鲁棒性和强大功能(补充图 20–25)。
- 通过将 gsMap 应用于高分辨率脑部 ST 数据集,我们生成了详细描述 30 种复杂性状与空间解析单细胞分辨率关联的性状-脑图。
- 这些图覆盖了进化保守的脑区(如海马体)和高级脑区(如大脑皮层),以及与认知、情感和行为相关的各种复杂性状,且我们的部分发现具有临床价值(补充说明,第 8 节)。
- 我们开发了一个交互式网络工具(https://yanglab.westlake.edu.cn/gsmap)来可视化和下载这些性状-脑关联图。
- 尽管我们在本研究中主要将 gsMap 作为一种遗传分析工具,但它可以集成到 ST 数据分析工具包中,用于探索单个细胞与人类疾病之间的关联,从而提供对空间模式化人类细胞疾病相关性的见解以改进诊断和治疗,同时还可揭示动物细胞与人类疾病的关联以实现更好的疾病建模。
- 此外,我们在 gsMap 中提供了一个选项,用于联合分析多个 ST 切片(方法),并证明了单独分析和联合分析结果之间的一致性(补充图 26)。
- 因此,我们建议使用联合分析模式以更好地解释从多个技术或生物学重复获得的 ST 数据中的 gsMap 结果。
Para_02
- 大脑中的细胞通常被归类为常规的细胞类型(例如,谷氨酸神经元)以研究它们与疾病的相关性。
- 我们的研究结果表明,由于细胞类型内部存在异质性,仅将细胞归类为这些常规细胞类型不足以理解它们在疾病中的作用。
- 例如,大脑中相同常规细胞类型的细胞在与性状或疾病的关联方面表现出广泛的异质性,而这种异质性与其空间分布相关(图3-5)。
- 因此,利用空间背景注释来捕捉细胞类型内部的异质性,不仅可以绘制出与性状相关的细胞在细胞类型内的空间分布图,还可以增强检测与性状相关的组织或细胞类型的能力(补充图20和21)。
- 我们的发现结合先前的形态学、电生理学和RNA测序研究,表明空间模式化的细胞类型内部异质性可能是哺乳动物大脑细胞的基本组织规则(补充说明,第9节)。
Para_03
- 另一个问题是 gsMap 是否依赖于空间上下文依赖的表达数量性状位点(eQTL)效应的假设。
- 上下文依赖的 eQTL 效应通常基于在每个上下文中标准化的基因表达水平进行检测,研究结果一致表明大多数 eQTL 并非上下文依赖。
- 然而,即使在标准化表达水平不存在上下文依赖的 eQTL 效应的情况下,由于不同细胞群体中原始表达水平方差的差异,eQTL 变异对复杂性状的影响可能因通过不同细胞群介导而有所不同(扩展数据图 10 和补充说明,第 10 节)。
- 探索空间上下文依赖的基因表达遗传调控需要群体规模的空间转录组(ST)数据,而这些数据目前尚不可用。
- 一种替代方法是解卷积批量样本中空间域的丰度,类似于细胞状态丰度的解卷积,然后通过交互模型(空间域与基因型交互)检测基因表达的空间依赖遗传效应。
Para_04
- 我们注意到本研究存在若干局限性,这些局限性可作为未来研究的潜在方向。
- 首先,由于目前可用的人类空间转录组数据有限,我们依赖于从小鼠和猕猴获得的脑空间转录组数据,以空间映射与人类复杂性状相关的脑细胞。
- 尽管我们的研究结果在不同的空间转录组数据集中得到了验证,并且在小鼠、猕猴和人类数据中展示了gsMap结果的一致性,但使用非人类的空间转录组数据不可避免地会因人类与这些动物模型之间的固有差异而导致统计功效的降低。
- 尽管来自大规模研究的GWAS汇总数据可能部分抵消这种功效损失,鉴于空间组学技术的快速发展以及相关成本的降低,我们预计未来研究将利用人类组织生成广泛且高质量的空间组学数据。
- 这些数据集将成为进一步验证我们结果和发现无法通过非人类空间转录组数据实现的新发现的宝贵资源。
- 其次,gsMap识别的性状与细胞之间的关联并不意味着存在因果关系。
- 由于细胞间的基因表达谱存在相关性,那些在生物学上并未对某一性状作出贡献的细胞也可能被识别为与该性状相关。
- 例如,我们使用小鼠胚胎空间转录组数据检测到背根神经节细胞与精神分裂症之间的关联,这可能是由背根神经节细胞与脑神经元之间基因表达值的相关性引起的。
- 第三,本研究专注于单独分析每个空间转录组切片或联合分析多个解剖学上同质的切片,而将多个解剖学上异质的空间转录组切片的联合分析留待进一步探索。
- 第四,将基因调控网络(例如,通路)纳入GSS计算可能会提高gsMap的功效。
- 最后,本研究的主要目标是应用gsMap探索性状相关细胞在脑中的空间分布,而对其他组织的探索则留待进一步研究。
Methods
Inclusion and ethics
包容与伦理
Para_01
- 本研究已获得西湖大学伦理委员会的批准(批准编号:20200722YJ001)。
gsMap method
gsMap 方法
Design principles
设计原则
Para_01
- 为了有效说明 gsMap,我们首先总结其设计原则。
- gsMap 利用 S-LDSC 框架来评估位于或靠近 ST 数据中特定位置特异性表达基因的遗传变异(主要是 SNPs)是否在与目标性状相关的遗传关联中富集。
- 为了精确估计单个位置的基因表达特异性,gsMap 会从同质位置聚合信息,这是关键的一步,因为 ST 数据中单个位置的基因表达谱通常具有稀疏性和高技术噪声。
- 仅使用空间坐标不足以识别同质位置,因为空间上相邻的位置不一定属于同一细胞类型。
- 仅使用基因表达谱也可能由于技术噪声导致对同质位置的识别存在偏差。
- 为了解决这些限制,gsMap 使用图神经网络 (GNN) 学习嵌入表示,该嵌入整合了空间坐标和基因表达谱,然后根据嵌入矩阵中的相似性为每个焦点位置识别同质位置。
- gsMap 随后通过聚合其同质位置的信息来估计每个焦点位置的基因表达特异性评分 (GSSs)。
- 尽管可以将现有工具生成的 ST 数据的位置嵌入应用于 gsMap,但为了增强我们软件工具的可用性并提高嵌入质量,我们开发了一个内置的 GNN 模型(见下文),该模型在可用时结合位置注释,并使用图注意力层 (GAT) 在识别焦点位置的同质位置过程中防止过度平滑。
Data input
数据输入
Para_01
- gsMap需要以下输入:(1) GWAS汇总统计;(2) 基于测序的空间转录组数据,包括全转录组基因表达谱和单个点的空间坐标;(3) LD参考数据;以及(4) 可选的SNP到基因的链接图谱。
- 我们注意到,不同平台的空间转录组数据分辨率差异显著,从测序芯片上的一个点代表一组细胞(例如,10X Visium),到亚细胞分辨率的点(例如,Stereo-seq)。
- 为了确保gsMap在处理亚细胞分辨率的空间转录组数据时的鲁棒性,需要进行细胞分割分析,该分析将测序芯片上的原始点合并为单个细胞。
- 为简化起见,我们继续使用‘点’来表示空间转录组数据中的数据点,其中在一个高分辨率空间转录组平台中,一个点代表一个单独的细胞(经过细胞分割分析后),而在传统空间转录组平台中,一个点代表多个细胞。
Processing gene expression data
处理基因表达数据
错误!!! - 待补充
Building a spatial graph of the spots
构建点的空间图
错误!!! - 待补充
Learning embedding matrix
学习嵌入矩阵
错误!!! - 待补充
Para_02
- gsMap 中的图注意力自编码器的损失函数包括两部分:(1) 用于重建基因表达矩阵的均方误差;以及 (2) 用于预测位点细胞类型的交叉熵损失。
Para_03
- 这里,({x"}{ig"}) 表示基因 g 在位置 i 的归一化表达值,而 ({\widehat{x"}}{ig"}) 是其重构值。
- ({p"}_{ik"}) 表示位置 i 属于细胞类型 k 的概率,如果位置 i 属于细胞类型 k,则编码为 1,否则为 0。
- ({q"}_{ik"}) 表示预测概率,c 是细胞类型的数量,γ 是平衡重构损失和交叉熵损失的超参数。
- 在实践中,γ 在大多数情况下被设置为 0.5。
- 在没有注释细胞类型的情况下,γ 被设置为 1。
- 值得注意的是,细胞类型注释被包括进来以提高嵌入的准确性,但它们并不是严格必要的,因为它们主要可以通过基因表达谱捕获(补充图 22)。
- 在训练过程中,使用 Adam 优化器来最小化损失函数,并使用指数线性单元(ELU)作为激活函数。
- 权重衰减设置为 10−4,最大迭代次数设置为 1,000。
- 当 (|L"}^{t+1}-{L"}^{t"}< {10^{-4}) 时,认为迭代已收敛。
Identification of homogeneous spots
同质斑点的识别
Para_01
- 嵌入矩阵 Z 整合了来自基因表达值、点的空间位置和细胞类型先验的信息。
- 然后,我们基于这些点在此潜在空间中的余弦相似性,为每个焦点识别出同质的点。
Para_02
- 其中 ({{\bf{Z"}}}_{i"}\in {{\mathbb{R"}}}^{m\times 1}) 是点 i 的嵌入向量。
- 对于每个焦点位置,我们选择与其余点中余弦相似度最高的 d 个点。
- 这一过程将每个单独的点与 d 个空间上接近且具有相似转录组谱的点对齐,这些点被称为微域 (D)。
- 在本研究中,对于由 10X Visium 平台生成的空间转录组数据,我们将 d 设置为 20;而对于由 Stereo-seq 平台生成的空间转录组数据,我们将 d 设置为 50。
Estimation of GSS
GSS的估计
Para_01
- 我们根据基因在单个位点的表达水平对其进行排序,较高的表达值获得较高的排名。
- 我们选择使用基因表达排名而不是基因表达值,因为排名在不同技术与生物学重复中对伪迹更具鲁棒性(补充图23),正如之前所建议的那样。
- 对于每个基因,通过计算其在焦点位点微域内的表达排名的几何平均值,并除以其在空间转录组数据所有位点中的表达排名的几何平均值,来评估其在每个焦点(个体)位点内的表达特异性。
- 对于基因g,({F"}_{ig"}) 表示其在位点i的表达特异性,计算公式如下:
Para_02
- 其中 ({R"}_{ig"}) 表示基因 g 在位置 i 的表达排名,(D_i) 表示位置 i 的微域(即同质位置的集合),d 表示位置的数量。
- 如果 ({F"}_{ig"} < 0),表明该基因在焦点位置内没有特异性表达,我们将其设置为 0。
- 此外,我们将每个基因在焦点位置微域中的表达比例与其在所有位置中的比例进行比较。
- 如果该比例小于 1,表明较大的 ({F"}{ig"}) 可能是由于少数异常高表达排名的位置导致的,我们同样将 ({F"}{ig"}) 设置为 0。
- 为了对齐 GWAS 汇总统计数据与估计的基因特异性之间的尺度,我们进行了指数投影,以进一步区分高表达特异性的基因。
Para_03
- 其中,({S"}_{ig"}) 表示基因 g 在位置 i 的最终特异性评分。
Mapping GSSs to SNPs
将 GSS 映射到 SNP
Para_01
- 我们将每个基因的特异性评分映射到其转录区域上下游 50 kb 窗口范围内的相应 SNPs。
- 通过敏感性分析,我们已经证明 gsMap 对不同窗口大小具有鲁棒性。
- 我们在 gsMap 软件工具中加入了一个选项,用于基于从表观基因组数据(例如 Roadmap 和 Activity-by-Contact 模型)衍生的 SNP 到基因链接图谱,将 100 kb 窗口外的 SNPs 映射到基因。
- 这一过程为每个位点生成了一组独特的 SNP 注释。
- 在本研究使用的 ST 数据集中,平均而言,每个位点有超过 150,000 个 SNPs 和 3,000 个基因具有非零 GSS。
Linking genomic annotations with GWAS data
将基因组注释与GWAS数据关联
Para_01
- 将每个位点视为一组SNP注释,gsMap使用S-LDSC框架评估具有较高GSS的SNP是否在目标性状的遗传力中富集,条件是基于基线SNP注释。
- gsMap中的S-LDSC可以被视为GWAS χ2统计量与使用来自单个位点的SNP注释计算的分层LD分数之间的线性回归分析。
- 关于S-LDSC的详细解释可以在补充说明的第2节中找到。
- 本研究中使用的LD参考数据来自千人基因组计划第三阶段(1KGP3)。
- 根据之前的研究,我们使用块刀切法来估计S-LDSC中回归系数的标准误差。
- P值通过单侧z检验计算,用于评估回归系数是否显著大于0。
- 较小的P值表明焦点位点与目标性状的相关性更强。
Estimating the strength of enrichment for a spatial region
估计某空间区域富集强度
Para_01
- 为了评估特定空间区域与目标性状之间的相关性,gsMap 使用柯西组合检验来聚合该空间区域内各个点的 P 值:
Para_02
- 其中 ({P"}_{i"}) 表示点 i 属于空间区域的 P 值。
- 该空间区域的聚合 P 值近似为:
Para_03
- 柯西组合检验通过整合空间域内所有位点的信号,与之前的域水平方法表现出高度一致性(中位数 r = 0.82),同时提供了更高的统计功效(补充图 21)。
Running time
运行时间
Para_01
- 我们优化了 gsMap 代码,以确保其在处理大规模 ST 数据时的效率。
- gsMap 在 3 个 CPU 小时内分析具有 12 万个位点的大规模 ST 数据集,而即使使用现有的单一位点 GSS,利用 S-LDSC 软件完成相同任务则需要超过 20,000 个 CPU 小时(每个位点约 10 分钟)。
- 关于 gsMap 每个步骤运行时间的摘要可在补充说明的第 3 节中找到。
Diagnostic tools
诊断工具
Para_01
- 我们在 gsMap 软件工具中加入了一个模块,该模块使用户能够绘制选定基因的表达水平和 GSS 的空间分布图,同时生成曼哈顿图以突出显示映射到这些基因的 SNP 关联信号,从而辅助数据的诊断分析。
Joint analysis mode
联合分析模式
Para_01
- 我们开发了一种联合分析模式,用于分析具有多个技术重复和/或生物学重复的ST数据,生成统一的结果并简化解释。
- 简而言之,每个位点的GSS是通过计算基因在同质位点中的平均排名,并通过其在所有位点(包括来自不同生物学重复的位点)中的平均排名进行标准化得出的。
- 该步骤实现了跨生物学重复的GSS注释的协调统一。
- 然后,我们进行遗传力富集分析(例如,使用S-LDSC)以将每个位点与性状关联起来,并应用柯西组合检验来整合具有相同注释的所有位点的P值。
Quantifying the association of a spatial region with a trait relative to all other regions
量化一个空间区域相对于所有其他区域与某特性之间的关联
Para_01
- 我们使用了一个比值比(OR),其计算方法是特定区域中与性状相关联的位点数与非关联位点数的比率,除以所有其他区域的相应比率,以此来比较具有不同GWAS统计功效的性状在特定空间区域的gsMap结果。
- 该指标量化了某一区域与某一性状之间的关联强度,相对于所有其他区域而言。
- 因此,它弥补了由于GWAS统计功效差异导致的gsMap结果的不同,从而能够更有意义地对跨性状的gsMap结果进行比较。
- OR的显著性通过卡方检验进行评估。
Simulations
模拟
Para_01
- 这项模拟研究使用了来自英国生物银行的真实基因型数据,涉及 10 万名欧洲血统的无关个体。
- 我们使用了 HapMap3 的单核苷酸多态性(SNP),并过滤掉了次要等位基因频率(MAF)< 0.01 或哈迪-温伯格平衡检验 P 值 < 10^-6 的 SNP,最终保留了总计 1,195,548 个 SNP。
- 我们使用 GCTA(V1.94.1)基于一组选定的因果变异的真实基因型数据生成定量性状。
- 我们模拟了四种零假设场景,其中因果变异分别:(1)随机分布在全基因组中;(2)富集在候选顺式调控元件(cCREs)中;(3)富集在连锁不平衡(LD)区块内;以及(4)富集在非空间可变基因(NSVGs)的顺式区域中。
- 这里的"零假设"意味着模拟的因果变异并未富集在具有上下文依赖表达的基因附近或周围,也就是说,模拟的性状预计不会与任何特定的位点群相关联。
Para_02
- 在情景1中,我们模拟了具有不同多基因性水平(即因果单核苷酸多态性(SNP)的比例)和遗传力(即表型变异中归因于因果SNP的比例)的表型。
- 因果SNP的数量从10,000到500,000不等,遗传力范围从0.1到0.6。
- 在情景2到情景4中,我们将因果SNP的数量固定为100,000,遗传力固定为0.3。
- 对于情景4,我们使用了两种策略来选择非空间变量基因(NSVGs)。
- 首先,我们选择了由gsMap计算得出的最大GSS小于1.5的基因。
- 其次,为了选择完全独立于gsMap分析过程的基因,我们使用SPARK-X来测试空间依赖性表达,并选择了P值大于0.85的基因。
- 每次模拟重复进行了三次。
- 我们使用PLINK(V1.90)将SNP与模拟表型关联起来,并将从SNP衍生出的前十个主成分作为协变量进行拟合。
GWAS summary statistics
GWAS 汇总统计
Para_01
- 我们从公共领域收集了涵盖广泛复杂性状和疾病的全基因组关联研究(GWAS)汇总统计量,这些性状和疾病跨越八大类:自身免疫、精神疾病、生殖、行为、代谢、血液学、人体测量学和癌症。
- 为了确保足够的GWAS统计功效,我们仅纳入那些LDSC遗传力估计的卡方统计量超过25的性状,这是基于之前研究的建议。
- 总之,我们分析了来自英国生物样本库(UK Biobank)及其他公开可用来源的110种复杂性状(包括疾病)的GWAS汇总统计量(平均样本量n = 385,000,补充表1)。
- 由于主要组织相容性复合体(MHC)区域的复杂性,我们将其从所有分析中排除。
Spatial transcriptomics datasets
空间转录组学数据集
Para_01
- 本研究纳入了六个空间转录组学数据集:小鼠胚胎数据(Stereo-seq)、人类胚胎数据(Stereo-seq)、小鼠脑数据(Stereo-seq 和 MERFISH)、猕猴皮层数据(Stereo-seq)以及人类 DLPFC 数据(10X Visium)。
- 为了将小鼠或猕猴基因与人类基因对齐,我们使用了 biomaRt(V3.18)R 包进行同源基因转换。
- 经过同源转换后,小鼠数据集的平均基因数量为 16,330,猕猴数据集的平均基因数量为 13,536。
- 按照标准分析流程,我们利用 Scanpy(V1.9.6)Python 包处理每个空间转录组学数据集。
- 每个数据集的详细信息总结如下。
Mouse embryonic ST data
小鼠胚胎ST数据
Para_01
- 我们分析了来自小鼠胚胎数据集的54个冠状切片,这些切片来源于8只C57BL/6J小鼠,每个切片平均包含81,125个点,涵盖了从胚胎阶段E9.5到E16.5的时间范围。
- 本研究包含了bin50分辨率(53个切片)和单细胞分辨率(1个切片)的数据。
- 我们获取了h5ad文件的访问权限,每个文件包括基因表达计数矩阵、细胞类型注释以及切片中点的空间坐标。
- 我们基于已知的标记基因验证了细胞类型注释的准确性。
Human embryonic ST data
人类胚胎ST数据
Para_01
- 我们分析了来自人类胚胎数据集的62个横切面,这些数据源自一个中国血统的CS8阶段人类胚胎,共包含38,562个位点。
- 我们可以访问h5ad文件,每个文件包含基因表达计数矩阵、细胞类型注释和位点空间坐标。
Mouse brain ST data
小鼠脑空间转录组数据
Para_01
- 我们分析了来自小鼠脑数据集的两个切片:一个来自成年C57BL/6J小鼠脑的冠状切片(包含50,140个点)和一个来自E16.5胚胎期C57BL/6J小鼠脑的矢状切片(包含65,303个点)。
- 这两个切片均达到了单细胞分辨率。
- 我们获取了h5ad文件,并使用已知的标记基因验证了细胞类型注释。
Mouse brain MERFISH ST data
小鼠脑MERFISH空间转录组数据
Para_01
- 我们分析了一个小鼠脑 MERFISH 空间转录组数据集,该数据集包含两个切片:一个来自成年 C57BL/6J-1 小鼠脑的冠状切片,包含 41,181 个点,以及一个来自成年 C57BL/6J-3 小鼠脑的矢状切片,包含 105,934 个点。
- 该数据集是通过基于图像的平台生成的,其中对约 1,100 个基因进行了成像。
- 其余基因由作者使用配对的单细胞 RNA 测序数据进行了推算。
- 我们获取了推算后的空间转录组数据集,其中包括 15,768 个基因的表达水平、细胞空间坐标、细胞类型注释和脑区注释。
Macaque cortical ST data
猕猴皮层ST数据
Para_01
- 我们分析了来自猕猴大脑皮层数据集的162个冠状切片,这些切片来源于三只成年雄性食蟹猴,每个切片平均包含266,654个点。
- 所有切片均为单细胞分辨率。
- 我们获取了sctransform的基因表达矩阵和元数据文件。
- 对于每个切片,我们对齐了点的空间坐标、细胞类型注释、皮层区域注释以及切片切割位置(EBZ),并将上述所有信息整合到一个h5ad文件中。
Human DLPFC ST data
人类DLPFC(背外侧前额叶皮层)ST数据
Para_01
- 我们分析了来自成人DLPFC数据集的八个冠状切片,这些数据来源于两名欧洲血统的捐赠者,每个切片平均包含3,973个位点。
- 这些数据是使用10X Visium生成的,每个位点包含几十个细胞。
- 我们可以访问h5ad文件,其中每个文件包括基因表达计数矩阵、位点空间坐标和皮层层次注释。
Single-cell RNA-seq datasets
单细胞 RNA 测序数据集
Para_01
- 在这项研究中,我们使用了六个单细胞RNA测序数据集来估算小鼠和人类之间的基因表达相关性。
- 这些数据集包含900万个细胞,涵盖了来自小鼠和人类的25种主要组织。
- 细胞类型注释由单细胞RNA测序数据集的作者提供,我们将其重新聚类为八大主要细胞类别:内皮细胞、上皮细胞、胶质细胞、神经元、免疫细胞、肌肉细胞、干细胞和其他细胞。
Comparison of the gsMap results from the human and macaque datasets
人类和猕猴数据集的 gsMap 结果比较
Para_01
- 我们将 gsMap 应用于八个人类 DPLFC ST 切片,并观察到这些切片之间结果高度一致。
- 然后我们为每个皮层计算了一个优势比 (OR),代表某一皮层与某种性状之间的关联强度,相对于所有其他皮层。
- 为了确保人类和猕猴结果之间的稳健比较,我们将八个 DPLFC ST 人类切片中每个皮层的中位 OR 值与九个 PFC ST 猕猴切片中的中位 OR 值进行比较。
- 考虑到人类和猕猴数据集中仅有五个匹配的皮层,我们汇总了来自 22 个脑相关性状的 OR 值,在比较分析中共得到 110 个数据点。
Integration of GWAS summary statistics with ST data using scDRS
使用 scDRS 将 GWAS 汇总统计与 ST 数据整合
Para_01
- scDRS9 是一种可以将 GWAS 汇总统计与单细胞 RNA 测序数据相结合以识别与性状相关细胞的方法。
- 简而言之,scDRS 计算性状富集分数,以检查一个细胞是否在一组与性状相关的基因中表现出过高的表达水平。
- 这些基因是通过基于基因的关联测试(例如 MAGMA76)从 GWAS 汇总统计数据中获得的。
- 为了评估性状富集分数的统计显著性,通过将该分数与具有匹配表达特征的对照基因计算出的分数进行比较,计算 P 值。
- 尽管 scDRS 最初是为单细胞 RNA 测序数据开发的,但原则上可以通过将空间转录组数据中的每个点视为一个细胞来应用于空间转录组数据。
- 按照 scDRS 的标准分析协议(V1.03),我们首先使用 MAGMA 从 GWAS 汇总统计数据生成基于基因的测试统计量,并利用从 1KGP3 获得的参考连锁不平衡数据,如同 gsMap 分析中所做的那样。
- 接下来,我们使用 scDRS 中的 ‘munge-gs’ 命令并将 ‘n-max’ 参数设置为 1,000,以生成一个包含从 MAGMA 分析中鉴定出的前 1,000 个性状相关基因的基因权重文件。
- 最后,我们使用 scDRS 中的 ‘compute-score’ 命令并将 ‘n-ctrl’ 参数设置为 1,000,通过蒙特卡罗采样生成 1,000 组对照基因,以获得每个点的性状富集 P 值。
Expression association analysis
表达关联分析
Para_01
- 我们进行了两项表达关联分析,以识别其在Glu神经元中的表达水平与其沿CA1背腹轴位置相关的基因,以及其在Glu神经元中的表达水平与精神分裂症相关性相关的基因。
Para_02
- 表达关联分析-1(针对基因 g):
Para_03
- 其中,({\bf{g"}}\in {{\mathbb{R"}}}^{k\times 1}) 表示基因在 k 个 Glu 神经元中的标准化表达水平,({\bf{x"}}\in {{\mathbb{R"}}}^{k\times 1}) 表示 k 个 Glu 神经元沿 D-V 轴的空间位置,β 为系数,({\bf{c"}}\in {{\mathbb{R"}}}^{k\times 1}) 表示位点文库大小向量(每个位点所有 UMI 计数的总和),用于校正因位点文库大小导致的基因表达变异,e 表示残差。
Para_04
- 表达关联分析-2(针对基因g):
Para_05
- 其中,({\bf{z"}}\in {{\mathbb{R"}}}^{k\times 1}) 表示 k 个谷氨酸神经元与目标性状(例如精神分裂症)之间关联的 -log10 gsMap P 值,γ 为系数,其余参数定义如上。
Gene prioritization for complex traits
复杂性状的基因优先排序
Para_01
- 我们将 gsMap 优先排序的基因与其他基因优先排序方法(包括 COLOC、FUSION、SMR 和 PoPS)所识别的基因进行了比较,这些方法针对精神疾病(包括抑郁症、双相情感障碍和精神分裂症),其中有许多药物靶点基因可用于验证。
- 对于 FUSION(V1.0.0),我们使用了来自 GTEx 的所有 eQTL 数据集;对于 COLOC(V5.2.3)和 SMR(V1.3.1),我们使用了来自 GTEx 和 MetaBrain 的所有 eQTL 数据集;而对于 PoPS(V0.2),我们使用了 PoPS 作者提供的特征矩阵。
Psychiatric drug targets
精神科药物靶点
Para_01
- 我们从DrugBank数据库和药物重定位中心数据库中收集了药物靶点基因。
- 对于DrugBank数据库,我们选择了用于治疗根据国际疾病分类(ICD)代码F00至F99分类的精神疾病的药物。
- 从药物重定位中心数据库中,我们选择了归类于精神疾病领域的药物。
- 总共,我们从那些已获批或正在进行临床试验以治疗精神疾病的药物中鉴定出417个靶点基因。
- 基于鉴定出的药物靶点基因,我们使用费舍尔精确检验评估在猕猴前额叶皮层14r区域高表达的基因是否在这些药物靶点基因中富集,相较于所有检测到的基因或其他前额叶皮层区域高表达的基因。
- 使用Seurat(V4.4.0)R包中的AddModuleScore函数,按照默认设置计算了这些417个药物靶点基因在每个位点的模块得分。
Data availability
Para_01
- 本研究使用了110个性状的GWAS汇总统计信息,详见补充表1。
- 小鼠胚胎和脑ST数据可在国家基因库数据库(CNGBdb)中获取,登录号为CNP0001543。
- 人类胚胎ST数据可在国家基因组数据中心(NGDC)基因组序列档案库(GSA)中获取,登录号为HRA005567。
- 小鼠脑MERFISH ST数据可通过https://cellxgene.cziscience.com/collections/0cca8620-8dee-45d0-aef5-23f032a5cf09获取。
- 猕猴皮质ST数据可在CNGBdb中获取,登录号为CNP0002035。
- 人类DLPFC ST数据可通过https://research.libd.org/globus/获取。
- 来自人类和小鼠的scRNA-seq数据可通过https://cellxgene.cziscience.com/获取。
- 参考LD数据(由1KGP3生成)可通过ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/2013050获取。
- LDSC基线功能注释数据可通过https://data.broadinstitute.org/alkesgroup/LDSCORE获取。
- DrugBank数据库可通过https://go.drugbank.com/releases/latest获取。
- 药物再利用中心数据库(2020年3月24日版本)可通过https://repo-hub.broadinstitute.org/repurposing#download-data获取。
- 人类、小鼠和猕猴参考基因组数据可通过https://nov2020.archive.ensembl.org/index.html获取。
- 不同性状和ST数据集的gsMap结果可通过https://yanglab.westlake.edu.cn/gsmap/home可视化并下载。
- 本论文提供了源数据。
Code availability
Para_01
- gsMap 的源代码可以在 GitHub (https://github.com/JianYang-Lab/gsMap) 和 Zenodo (https://doi.org/10.5281/zenodo.14744887) 上获取。