单细胞现在研究到什么地步了？

wolf

想了解一下单细胞最前沿的研究

原文地址：https://www.zhihu.com/question/424463187

长长的路

R与Python，在生物信息学领域的博弈异常激烈。许多生信分析，两个阵营都发展出了自己的方法，比如单细胞数据分析，R有Seurat，Python就有Scanpy。这些层出不穷的方法不断地吸引着吃瓜群众的眼球，同时也让人患上了选择困难症。
到底谁优谁劣？一时竟难分高下。今天我们就以吃瓜的心态来看看两个阵营在单细胞数据分析方面的博弈。
R语言单细胞分析最佳实践
首先出场的是 Orchestrating Single-Cell Analysis with Bioconductor （Bioconductor OSCA）。由 Bioconductor 出品，旨在传授基于 R 的 Bioconductor 生态系统分析单细胞 RNA-Seq 的常见工作流程。其电子书网址是：https://bioconductor.org/books/release/OSCA/。




一篇同名论文[Amezquitaet al., 2020] 概述了使用 Bioconductor 进行单细胞分析。但在线的电子书含有大量示例代码，更详细。不过，OSCA 不包括其他单细胞组学，例如 scATAC-seq。空间转录组学在补充的 Orchestrating Spatially-Resolved Transcriptomics Analysis with Bioconductor （Bioconductor OSTA） 一书 （https://lmweber.org/OSTA-book/） 中进行了介绍。由于这些书籍是为 Bioconductor 生态系统设计的，因此它们仅使用 Bioconductor 上可用的工具。这些并不一定能得出书中所描述的最佳分析结果。Bioconductor 书籍对具有基本 R 和更强生物学背景的人特别有用，他们有兴趣学习如何使用 Bioconductor 分析单细胞和空间转录组学数据。
Python单细胞分析最佳实践
Python也不遑多让，好事者也推出了基于 Python 的单细胞分析最佳实践。电子书网址：
https://www.sc-best-practices.org/preamble.html，该书详细介绍了单细胞数据分析。有理论讲解和大量示例代码，并且一直处于更新当中。




作者团队先是在 2019 年发表了一篇论文：Current best practices in single‐cell RNA‐seq analysis: a tutorial，介绍了单细胞 RNA-Seq 分析的最佳实践。与Bioconductor OSCA 相比，这篇论文和示例分析没有受到其展示的工具的偏见，并且在涵盖主题的广度方面内容更加完整。
2023年，同一团队发表了新的论文：Best practices for single-cell analysis across modalities。提出了最新的最佳实践。其示例代码以 Python 为主，可以说是与 Bioconductor OSCA 针锋相对了。大家该如何选择呢？




海纳百川
正所谓海纳百川，有容乃大。我们在选择工具的时候，可以有偏见，但不能固执。Galaxy生信云平台（http://UseGalaxy.CN）就汇聚了大量生物信息学工具，不管是 R，Python 抑或是其他语言开发的，只要足够优秀，平台都做了整合，供大家自由选择。单细胞方面，就整合了基于 Seurat 和 Scanpy 的最佳分析实践，欢迎大家使用。




参考文献

• Orchestrating single-cell analysis with bioconductor. https://bioconductor.org/books/release/OSCA/, 2022.
  
• Orchestrating single-cell analysis with bioconductor. Nature Methods， 2022
  
• A single-cell survey of the small intestinal epithelium. Nature, 2017
  
• Orchestrating high-throughput genomic analysis with bioconductor. Nature Methods， 2015
  
• Current best practices in single-cell term`rna`-seq analysis: a tutorial. Molecular Systems Biology， 2019
  

同花顺

不管是自测还是公共数据库，测序等组学数据是大家手头上最容易拿到的数据。先前大家通过生信分析可以批量水文章，现在组学数据和生信分析的意义正在发生改变，更重要的价值是如何从组学数据中发掘出“亮点”。但是，组学数据的特点是又多又杂，组学数据的解读和梳理、判断就是第一步，而后续的工作才是如何基于组学数据讲好故事。
今天我就来抛砖引玉分享一下，当拿到组学数据的时候，如何从千头万绪的结果中挑选到线头并把故事讲好，步骤大致分为四步：
1质控判断：从数据质量、样本、分组等整体情况看结果是否可信和可靠。
比如常见的2组3样本（3vs3 ）细胞的RNA测序结果，如果没有特殊情况下，组内差异（比如三个Con对照、三个药物组之间）一般小于组间差异的（药物组vs Con），这一点可以从样本PCA聚类等结果看出来，如果三个Con之间的差异要比组间差异还大，下一步筛选差异基因的时候用P值<0.05大概率很难筛选到很多差异基因（P值太大不符合<0.05的情况）。
当然，与均一性比较好的细胞相比，临床样本的异质性要大很多，特别是样本量比较大的情况下，出现离群样本很正常，这才是临床真实的情况，那离群样本如何处理呢？有的团队做法是舍弃这个样本，让数据看起来更好；有的则保留这个样本，但在讨论中说明可能的原因。不管是自测还是公共数据库，测序等组学数据是大家手头上最容易拿到的数据。先前大家通过生信分析可以批量水文章，现在组学数据和生信分析的意义正在发生改变，更重要的价值是如何从组学数据中发掘出“亮点”。但是，组学数据的特点是又多又杂，组学数据的解读和梳理、判断就是第一步，而后续的工作才是如何基于组学数据讲好故事。
今天我就来抛砖引玉分享一下，当拿到组学数据的时候，如何从千头万绪的结果中挑选到线头并把故事讲好，步骤大致分为四步：
1质控判断：从数据质量、样本、分组等整体情况看结果是否可信和可靠。
比如常见的2组3样本（3vs3 ）细胞的RNA测序结果，如果没有特殊情况下，组内差异（比如三个Con对照、三个药物组之间）一般小于组间差异的（药物组vs Con），这一点可以从样本PCA聚类等结果看出来，如果三个Con之间的差异要比组间差异还大，下一步筛选差异基因的时候用P值<0.05大概率很难筛选到很多差异基因（P值太大不符合<0.05的情况）。
当然，与均一性比较好的细胞相比，临床样本的异质性要大很多，特别是样本量比较大的情况下，出现离群样本很正常，这才是临床真实的情况，那离群样本如何处理呢？有的团队做法是舍弃这个样本，让数据看起来更好；有的则保留这个样本，但在讨论中说明可能的原因。
总体来说，这个步骤是质控步骤，也是后续分析的基石。另外一个例子就是单细胞RNA-Seq分析完成后需要看整个UMAP（或者TSNE）是否“干净”，如果在Cluster周围有很多散点、分叉等，后续分析得到的结论有时候会似是而非，所以很多团队会在质控这一步骤不断调整，直到Cluster清楚干净。这一步非常重要，需要先排除样本弄反、分组弄错、样本收集储存等过程的问题，然后才是对数据进行分析和解读。
2结论判断：初步看一下结果是否大致符合预期。
我们在看数据之前可以先思考一下：我关注的这个基因、药物或者疾病，如果差异基因、功能富集分析的通路和生物学过程出来后，我期待看到什么结果？
这个怎么判断呢？有两个参考：a是数据库和文献中关于基因、药物或者疾病的报道，比如A基因可能与细胞周期有关，那大致可以预期看到A基因沉默后的差异基因和通路与细胞周期有关；b. 是前期开展的功能实验，比如我们看到A基因沉默后细胞凋亡比较明显，那总要看到一些跟凋亡或者细胞死亡有关的通路吧。如果我们看到这个基因A在沉默后有对应的基因、通路富集（活化或者抑制），就可以在文章里面说明：如我们预期的（As expected）, RNA-Seq测序（或者其它组学结果）显示A 基因沉默后导致了细胞周期（或者凋亡）通路或者基因的富集和改变，这与先前……的报道一致（或者与我们功能实验的结果一致，一般是Consistent with ……）。
当然，另外一种情况就是：结果越看越奇怪。比如本来实验是A基因沉默（并且用qPCR和WB验证过了）的细胞做的测序，怎么这个沉默组里面A基因表达比对照组还高表达呢，A本来是促进细胞死亡的，怎么A沉默后细胞死亡有关的通路活性都升高了，是不是分组弄反了？！
另外，怎么这个数据呈现出来的差异基因和通路，与我关注的疾病没什么关系啊？明明做的是炎症和损伤的疾病，怎么连细胞因子都没差异啊？差异的基因跟疾病的病理特点都对不上啊？这个时候就要核实一下是不是分组、数据弄反了，或者样本的问题了。
3思考和查询：哪些结果是“超出预期”的？
第三步是在第二步的基础上进一步看数据：除了符合预期的结果外，哪些是没有想到的？就像文章里面经常出现的词：Surprisingly, Notably……，我们发现……。看到这里我们就知道要引出文章的重要创新点了。
而对大家来说，这个步骤就是我们创新点真正开始的地方：从数据分析引出主要发现。比如我们发现某个药物处理细胞后，除了符合预期的细胞铁死亡、细胞周期改变外，富集最显著的居然是调控细胞骨架的通路，假如我们的知识背景是细胞骨架只与细胞运动有关，就会认为是新发现；当然查文献或者看到推文（这个“平平无奇”的信号通路，居然能与这三个顶级热点联系起来，你上车了没？）后，就会想到双硫死亡，这样药物作用就有了创新点：药物不仅影响细胞铁死亡，还能通过调控细胞骨架诱导双硫死亡。
当然，有个观念很重要：做研究和讲故事的顺序和逻辑是不同的。做研究探索的时候，先有初步假设，然后做关键实验（我们可以理解为实验检查点：通过几个关键实验初步判断假设是否成立），再根据关键实验调整假说，先把大框架完成（也就整个工作量的30%），最后再差缺补漏补充更多实验数据；而讲故事的时候，则基本拿到完整数据了，需要考虑的是：研究的亮点是什么？怎么把数据和结论包装好，让编辑和审稿人更感兴趣？
回到上面整个问题，超出预期的发现可能会有很多：不过不用都列出来，一般一篇文章里面考虑2-3个新发现就够了，之所以考虑2-3个点是因为不能保证1个新发现就能验证出来；当然你可以挑选9个新发现，然后通过进一步分析再筛选到5、6个，最后再做实验验证，最后写成几篇文章就行了。需要说明的是：这里2-3个“新发现”需要结合文献报道来确认；另外，大家还是需要对KEGG通路等要有一些基本的了解，实在不行就挨个查文献，这里确实需要对信号转导、通路和生物学过程有一些背景知识。
4基于新发现的进一步分析和验证。
完成了第三步，相当于从一团乱麻中梳理出来几个“线头（新发现）”，下面其实也分两步：a) 分析层面的进一步验证；b) 实验层面的验证。
分析层面的验证其实不仅是用更多的数据集支持这个结论，还包括搞清楚新发现的大致情况，比如到底有多少基因发生了改变，其中最重要的基因是什么，以及A通路里面富集到的基因与其他（B）通路富集的基因是不是很多都重叠？比如A通路富集到5个基因（都是一个家族的），其中4个基因在另外一条不相关的B和C通路中也被富集到，那这个发现就往后放（原因是这些基因的功能指向性不强，既能A，又能B和C，可能只是偶然被富集到）。
另外，如果差异基因比较少的情况下做的功能富集（比如100个基因），虽然能富集到某个通路上只有3、4个基因，这个结果的参考价值也要打个问号，这时候需要考虑的是：为什么疾病与对照组相比只有100个显著差异基因？
一般情况下，A基因沉默或过表达的差异基因会比疾病vs对照少一些，毕竟只是一个基因的变化造成的扰动，但是疾病如果与对照只有100个显著差异基因，就要思考一下原因了：除非这100个RNA的变化就能反映疾病差异，或者疾病与对照的差异要在蛋白层面（而非RNA层面）显示出来。
最后一个实验验证就不展开说了，大家筛选到基因、通路或者功能以后如何做都比较清楚，干扰过表达、抑制剂激动剂，最后做功能机制等等。
总体来说，这个步骤是质控步骤，也是后续分析的基石。另外一个例子就是单细胞RNA-Seq分析完成后需要看整个UMAP（或者TSNE）是否“干净”，如果在Cluster周围有很多散点、分叉等，后续分析得到的结论有时候会似是而非，所以很多团队会在质控这一步骤不断调整，直到Cluster清楚干净。这一步非常重要，需要先排除样本弄反、分组弄错、样本收集储存等过程的问题，然后才是对数据进行分析和解读。
2结论判断：初步看一下结果是否大致符合预期。
我们在看数据之前可以先思考一下：我关注的这个基因、药物或者疾病，如果差异基因、功能富集分析的通路和生物学过程出来后，我期待看到什么结果？
这个怎么判断呢？有两个参考：a是数据库和文献中关于基因、药物或者疾病的报道，比如A基因可能与细胞周期有关，那大致可以预期看到A基因沉默后的差异基因和通路与细胞周期有关；b. 是前期开展的功能实验，比如我们看到A基因沉默后细胞凋亡比较明显，那总要看到一些跟凋亡或者细胞死亡有关的通路吧。如果我们看到这个基因A在沉默后有对应的基因、通路富集（活化或者抑制），就可以在文章里面说明：如我们预期的（As expected）, RNA-Seq测序（或者其它组学结果）显示A 基因沉默后导致了细胞周期（或者凋亡）通路或者基因的富集和改变，这与先前……的报道一致（或者与我们功能实验的结果一致，一般是Consistent with ……）。
当然，另外一种情况就是：结果越看越奇怪。比如本来实验是A基因沉默（并且用qPCR和WB验证过了）的细胞做的测序，怎么这个沉默组里面A基因表达比对照组还高表达呢，A本来是促进细胞死亡的，怎么A沉默后细胞死亡有关的通路活性都升高了，是不是分组弄反了？！
另外，怎么这个数据呈现出来的差异基因和通路，与我关注的疾病没什么关系啊？明明做的是炎症和损伤的疾病，怎么连细胞因子都没差异啊？差异的基因跟疾病的病理特点都对不上啊？这个时候就要核实一下是不是分组、数据弄反了，或者样本的问题了。
3思考和查询：哪些结果是“超出预期”的？
第三步是在第二步的基础上进一步看数据：除了符合预期的结果外，哪些是没有想到的？就像文章里面经常出现的词：Surprisingly, Notably……，我们发现……。看到这里我们就知道要引出文章的重要创新点了。
而对大家来说，这个步骤就是我们创新点真正开始的地方：从数据分析引出主要发现。比如我们发现某个药物处理细胞后，除了符合预期的细胞铁死亡、细胞周期改变外，富集最显著的居然是调控细胞骨架的通路，假如我们的知识背景是细胞骨架只与细胞运动有关，就会认为是新发现；当然查文献或者看到推文（这个“平平无奇”的信号通路，居然能与这三个顶级热点联系起来，你上车了没？）后，就会想到双硫死亡，这样药物作用就有了创新点：药物不仅影响细胞铁死亡，还能通过调控细胞骨架诱导双硫死亡。
当然，有个观念很重要：做研究和讲故事的顺序和逻辑是不同的。做研究探索的时候，先有初步假设，然后做关键实验（我们可以理解为实验检查点：通过几个关键实验初步判断假设是否成立），再根据关键实验调整假说，先把大框架完成（也就整个工作量的30%），最后再差缺补漏补充更多实验数据；而讲故事的时候，则基本拿到完整数据了，需要考虑的是：研究的亮点是什么？怎么把数据和结论包装好，让编辑和审稿人更感兴趣？
回到上面整个问题，超出预期的发现可能会有很多：不过不用都列出来，一般一篇文章里面考虑2-3个新发现就够了，之所以考虑2-3个点是因为不能保证1个新发现就能验证出来；当然你可以挑选9个新发现，然后通过进一步分析再筛选到5、6个，最后再做实验验证，最后写成几篇文章就行了。需要说明的是：这里2-3个“新发现”需要结合文献报道来确认；另外，大家还是需要对KEGG通路等要有一些基本的了解，实在不行就挨个查文献，这里确实需要对信号转导、通路和生物学过程有一些背景知识。
4基于新发现的进一步分析和验证。
完成了第三步，相当于从一团乱麻中梳理出来几个“线头（新发现）”，下面其实也分两步：a) 分析层面的进一步验证；b) 实验层面的验证。
分析层面的验证其实不仅是用更多的数据集支持这个结论，还包括搞清楚新发现的大致情况，比如到底有多少基因发生了改变，其中最重要的基因是什么，以及A通路里面富集到的基因与其他（B）通路富集的基因是不是很多都重叠？比如A通路富集到5个基因（都是一个家族的），其中4个基因在另外一条不相关的B和C通路中也被富集到，那这个发现就往后放（原因是这些基因的功能指向性不强，既能A，又能B和C，可能只是偶然被富集到）。
另外，如果差异基因比较少的情况下做的功能富集（比如100个基因），虽然能富集到某个通路上只有3、4个基因，这个结果的参考价值也要打个问号，这时候需要考虑的是：为什么疾病与对照组相比只有100个显著差异基因？
一般情况下，A基因沉默或过表达的差异基因会比疾病vs对照少一些，毕竟只是一个基因的变化造成的扰动，但是疾病如果与对照只有100个显著差异基因，就要思考一下原因了：除非这100个RNA的变化就能反映疾病差异，或者疾病与对照的差异要在蛋白层面（而非RNA层面）显示出来。
最后一个实验验证就不展开说了，大家筛选到基因、通路或者功能以后如何做都比较清楚，干扰过表达、抑制剂激动剂，最后做功能机制等等。

继续前进

以2024年的文献做剖面，单细胞数据积累到了要AI出手的地步。

Cole Trapnell 最近的综述里说，

Trapnell, C. Revealing gene function with statistical inference at single-cell resolution. Nat Rev Genet (2024). https://doi.org/10.1038/s41576-024-00750-w



如果我们要使用模型的内部结构来指导实验室的机理实验，我们必须要求的不仅仅是统计方法和人工智能（AI）工具的准确性。他们不仅要学习生物学，还要能够用我们理解的术语将其传达给我们。因此，人工智能社区正在大力努力使人工智能模型的输出更容易解释

可以看出现在数据量积累和维度的复杂程度已经到了非常有挑战性的阶段，工具的革新可能是新发现的关键。
顶刊新工具

Nature Methods 在力推单细胞的AI工具，这里汇总了一下2024年上半年的工具包及github链接：

1 Cui, H., Wang, C., Maan, H. et al. Nat Methods (2024).  【多大】

scGPT: toward building a foundation model for single-cell multi-omics using generative AI.
生成式预训练模型在语言和计算机视觉等各个领域都取得了显著的成功。具体来说，大规模多样化数据集和预训练转换器的结合已成为开发基础模型的一种有前途的方法。本研究在语言和细胞生物学之间建立了相似之处（与文本由单词组成类似：细胞由基因定义），本研究探讨了基础模型在推进细胞生物学和遗传研究方面的适用性。利用新兴的单细胞测序数据，本研究构建了一个单细胞生物学的基础模型，即scGPT，该模型基于一个跨越超过3300万个细胞的存储库的生成式预训练转换器。本研究结果表明，scGPT有效地提炼了有关基因和细胞的关键生物学见解。通过对迁移学习的进一步调整，可以优化scGPT，从而在各种下游应用中实现卓越的性能。这包括细胞类型注释、多批次整合、多组学整合、扰动反应预测和基因网络推断等任务。

https://github.com/bowang-lab/scGPT
2 Li, M.M., Huang, Y., Sumathipala, M. et al. Contextual AI models for single-cell protein biology. Nat Methods (2024). 【MIT;NASA】

Contextual AI models for single-cell protein biology了解蛋白质功能和开发分子疗法需要破译蛋白质起作用的细胞类型以及蛋白质之间的相互作用。然而，对于现有算法来说，跨生物环境对蛋白质相互作用进行建模仍然具有挑战性。在这里，我们介绍PINNACLE，这是一种几何深度学习方法，可以生成上下文感知的蛋白质表示。利用多器官单细胞图谱，PINNACLE 在上下文中的蛋白质相互作用网络上学习，从 24 个组织的 156 个细胞类型上下文中生成 394,760 个蛋白质表示。PINNACLE的嵌入空间反映了细胞和组织结构，从而实现了组织层次结构的零样本检索（zero-shot retrieval）。预训练的蛋白质表征可以适用于下游任务：增强基于 3D 结构的表征，以解决免疫肿瘤蛋白相互作用，并研究药物对细胞类型的影响。PINNACLE 在用于类风湿性关节炎和炎症性肠病的治疗靶点方面优于最先进的模型，并确定了比无上下文模型具有更高预测能力的细胞类型背景。PINNACLE能够根据其所处的环境调整其输出，这为生物学中大规模的特定环境预测铺平了道路。

https://github.com/mims-harvard/PINNACLE.
3 Hao, M., Gong, J., Zeng, X. et al. Large-scale foundation model on single-cell transcriptomics. Nat Methods (2024). 【清华;百图生科biomap】

Large-scale foundation model on single-cell transcriptomics
大型预训练模型已成为基础模型，从而在自然语言处理和相关领域取得突破。开发用于破译细胞“语言”和促进生物医学研究的基础模型是有希望的，但也具有挑战性。本文开发了一个大型预训练模型 scFoundation，具有涵盖约20,000个基因的1亿个参数，在超过5000万个人类单细胞转录组图谱上进行了预训练。scFoundation 在可训练参数的大小、基因的维度和训练数据量方面是一个大规模模型。其非对称transformer-like架构和预训练任务设计能够有效地捕获各种细胞类型和状态下基因之间的复杂上下文关系。实验表明，它作为基础模型的优点，在基因表达增强、组织药物反应预测、单细胞药物反应分类、单细胞扰动预测、细胞类型注释和基因模块推断等各种单细胞分析任务中取得了最先进的性能。
Large-scale foundation model on single-cell transcriptomics

4 Hou, W., Ji, Z. Assessing GPT-4 for cell type annotation in single-cell RNA-seq analysis. Nat Methods (2024). 【Duke】

Assessing GPT-4 for cell type annotation in single-cell RNA-seq analysis本文证明了大型语言模型 GPT-4 可以在单细胞 RNA 测序分析中使用标记基因信息准确注释细胞类型。当对数百种组织和细胞类型进行评估时，GPT-4 生成的细胞类型注释与手动注释表现出很强的一致性。此功能可以大大减少细胞类型注释所需的工作量和专业知识。此外，我们还开发了一个 R 软件包 GPTCelltype，用于 GPT-4 的自动细胞类型注释。

GPTCelltype_Paper这几篇文章都是很年轻的科学家在推进：



上面四篇文章的（部分）PI

作为十几年来生信方法的风向标，Nature Methods 的收稿暗示了，单细胞数据的AI驱动，很大概率会成为一个新的研究热点和应用趋势。

应用方面，肿瘤领域的临床应用从来都是跑到最前面的：
临床方面

• Sinha S, Vegesna R, Mukherjee S, Kammula AV, Dhruba SR, Wu W, Kerr DL, Nair NU, Jones MG, Yosef N, Stroganov OV, Grishagin I, Aldape KD, Blakely CM, Jiang P, Thomas CJ, Benes CH, Bivona TG, Schäffer AA, Ruppin E. PERCEPTION predicts patient response and resistance to treatment using single-cell transcriptomics of their tumors. Nat Cancer. 2024 Jun;5(6):938-952. doi: 10.1038/s43018-024-00756-7. Epub 2024 Apr 18. PMID: 38637658.
  

PERCEPTION predicts patient response and resistance to treatment using single-cell transcriptomics of their tumors - PubMed (nih.gov)为个体癌症患者量身定制最佳治疗仍然是一项重大挑战。为了解决这个问题，我们开发了PERCEPTION（PERsonalized Single-Cell Expression-Based Planning for Treatments In ONcology），这是一种精准的肿瘤学计算管道。我们的方法使用来自大规模细胞系药物筛选的公开匹配的批量和单细胞表达谱。这些配置文件有助于建立基于患者单细胞肿瘤转录组学的治疗反应模型。PERCEPTION在预测培养和患者肿瘤来源的原代细胞对靶向治疗的反应方面取得了成功，并在多发性骨髓瘤和乳腺癌的两项临床试验中取得了成功。它还捕获了接受酪氨酸激酶抑制剂治疗的肺癌患者的耐药性发展。在所有临床队列中，PERCEPTION 的表现优于已发表的最先进的基于单细胞和基于bulk的预测因子。我们的工作展示了使用肿瘤的单细胞表达谱对患者进行分层，将鼓励在临床环境中采用单细胞组学分析，增强基于单细胞组学的精准肿瘤学工具。
PERCEPTION2.  Zargari A, Topacio BR, Mashhadi N, Shariati SA. Enhanced cell segmentation with limited training datasets using cycle generative adversarial networks. iScience. 2024 Apr 12;27(5):109740. doi: 10.1016/j.isci.2024.109740. PMID: 38706861; PMCID: PMC11068845.
Enhanced cell segmentation with limited training datasets using cycle generative adversarial networks - PubMed深度学习正在改变生物图像分析，但由于缺乏大型、多样化的注释数据集，其在单细胞分割中的应用受到限制。我们通过引入基于 CycleGAN 的架构 cGAN-Seg 来解决这个问题，该架构使用有限的注释数据集增强了细胞分割模型的训练。在训练过程中，cGAN-Seg 生成带注释的合成相差或荧光图像，其形态细节和细微差别与真实图像非常相似。这增加了分割模型的可变性，增强了合成样本的真实性，从而提高了预测准确性和泛化性。实验结果表明，与传统训练技术相比，cGAN-Seg显著提高了广泛使用的分割模型的性能。我们的方法有可能加速显微镜图像分析基础模型的开发，表明其在通过高效的训练方法推进生物图像分析方面具有重要意义。（这个方向的AI是进步最快的，算是和单细胞也相关）
GitHub - abzargar/cGAN-Seg3.  immunity ： AI and immunology VOICES | VOLUME 57, ISSUE 6, P1177-1181, JUNE 11, 2024   DOI:https://doi.org/10.1016/j.immuni.2024.05.016
人工智能正迅速成为日常生活许多方面的一部分，其影响遍及所有研究领域。我们请研究人员分享他们对人工智能如何改变免疫学研究的看法，前进的必要条件，潜力和陷阱，以及随着该领域进入新时代，什么将保持不变。（前瞻性讨论）
Dominic Grün的观点：（免疫学中）从细胞到组织：最近用于分析单细胞基因组学数据的新型人工智能工具激增，彻底改变了我们整合大量异质数据集的能力，以构建组织图谱并整合染色质可及性和基因表达等数据模式。目前，我们正在目睹单细胞分辨率空间数据的出现，这些数据揭示了细胞相互作用如何影响组织结构和由组织生态位驱动的细胞状态反应。未来的一个关键挑战将是超越以细胞为中心的模型，并利用生成式人工智能开发全组织模型，使我们能够以单细胞分辨率预测组织对特定干预措施的反应。这些模型将更深入地了解组织环境中细胞扰动的原因和后果。这一目标要求在开发可解释的人工智能模型方面取得进一步进展。尽管利用变分自编码器或转换器模型的现有方法可以产生多维数据的信息表示，作为多模态整合的基础，但对基因调控的机理见解受到设计的限制。克服这些局限性的一个有希望的途径是设计人工智能代理，这些代理可以被训练以主动学习受生成任务约束的因果作用，例如重建组织中细胞状态的耦合集合。然后，这些药物可以探索单个细胞中明确的干预措施对整个组织的影响，并大大扩展我们的范围，以更深入地了解免疫细胞与其组织环境在健康和疾病中的相互作用。这将使突破性应用成为可能，例如免疫疗法对肿瘤微环境影响的in silico测试。
4. New study uses AI technology to understand cell dysfunction in ALS
New study uses AI technology to understand cell dysfunction in ALS » College of Public Health & Health Professions » College of Public Health and Health Professions » University of Florida (ufl.edu)来自佛罗里达大学公共卫生与卫生专业学院和英国谢菲尔德大学的一组研究人员获得了运动神经元疾病协会将单细胞分析与称为深度学习的人工智能方法相结合，以详细和系统的方式分析受 ALS（也称为 Lou Gehrig 病）影响的不同类型的中枢神经系统细胞。
“我们希望发现有关这些细胞如何相互作用以及它们的功能障碍如何与运动神经元退化相关的重要信息，”该研究的联合首席研究员说Sai Zhang，博士，UF PHHP的助理教授流行病学系.“这些知识对于寻找治疗ALS和开发靶向疗法的新方法至关重要。

数据集的整合

数据越多，AI越强
1，Pan, L., Parini, P., Tremmel, R. et al. Single Cell Atlas: a single-cell multi-omics human cell encyclopedia. Genome Biol 25, 104 (2024). https://doi.org/10.1186/s13059-024-03246-2

http://www.singlecellatlas.org.单细胞测序数据集是生物学和医学的关键，它以前所未有的分辨率揭示了对异质细胞群的见解。在这里，我们通过深入表征来自 125 个健康成人和胎儿组织的五个单细胞组学、空间转录组学和两个体组学的数据集，构建了人体组织的单细胞多组学图谱。我们构建了其基于补体网络的平台，即单细胞图谱（SCA，http://www.singlecellatlas.org），以实现对人类胎儿和成人组织中深度多组学特征的广泛交互式数据探索。图谱资源和数据库查询旨在为各种组学研究提供一站式、全面且省时的资源。
https://github.com/eudoraleer/sca文献挖掘，单细胞中的基因知识关联
2. bioRxiv  [Preprint]. 2024 Mar 5:2023.10.16.562533. doi: 10.1101/2023.10.16.562533. （预印）
最近在利用大规模基因表达数据开发单细胞生物学基础模型方面取得了重大进展。Geneformer 和 scGPT 等模型从数百万个细胞的基因表达谱中隐式学习基因和细胞功能，这需要大量的数据管理和资源密集型训练。在这里，我们通过利用基于文献的 ChatGPT 基因嵌入来探索一个更简单的替代方案。我们的提案 GenePT 使用 GPT-3.5 对单个基因的 NCBI 文本描述来生成基因嵌入。从那里，GenePT以两种方式产生单细胞包埋：（i）通过平均基因包埋，按每个基因的表达水平加权;或 （ii） 使用按表达水平排序的基因名称，为每个细胞创建一个嵌入句子。无需数据集管理和额外的预训练，GenePT 高效且易于使用。在用于评估最新单细胞基础模型的许多下游任务中，例如，对基因特性和细胞类型进行分类，GenePT的性能与Geneformer和其他模型相当，而且通常更好。GenePT表明，文献的大型语言模型嵌入是生物基础模型的一条简单而有效的途径。
https://github.com/yiqunchen/GenePT3. Yuan, Q., Duren, Z. Inferring gene regulatory networks from single-cell multiome data using atlas-scale external data. Nat Biotechnol (2024). https://doi.org/10.1038/s41587-024-02182-7
Inferring gene regulatory networks from single-cell multiome data using atlas-scale external data现有的基因调控网络 （GRN） 推断方法仅依赖于基因表达数据或较低分辨率的批量数据。尽管最近整合了染色质可及性和RNA测序数据，但从有限的独立数据点学习复杂的机制仍然是一个艰巨的挑战。在这里，我们介绍了 LINGER（用于基因调控的终身神经网络），这是一种机器学习方法，用于从单细胞配对基因表达和染色质可及性数据中推断 GRN。LINGER将不同细胞环境中的图谱规模外部批量数据和转录因子基序的先验知识整合为流形正则化。与现有方法相比，LINGER的准确度相对提高了四到七倍，并揭示了全基因组关联研究的复杂调控格局，从而能够增强对疾病相关变异和基因的解释。根据参考单细胞多组数据进行GRN推断，LINGER能够仅从批量或单细胞基因表达数据中估计转录因子活性，利用丰富的可用基因表达数据来识别病例对照研究中的驱动调节因子。
https://github.com/Durenlab/LINGER业界动态

1. 赛诺菲首席执行官：人工智能有望开启一个伟大的药物发现时代，从根本上改变医学，但前提是我们允许它实现 （2024/02/19）
Sanofi CEO: AI promises a great era of drug discovery that could fundamentally change medicine–but only if we allow it to deliver on its promise 2. Deepcell 和 NVIDIA 合作推进单细胞研究中的 AI （2024/1/9）
Deepcell 和 NVIDIA 合作推进单细胞研究中的 AI - PharmaTimes3.  bridgeinformatics 从单细胞测序到生物学故事：单细胞测序和人工智能如何彻底改变生物学 （2024/6/14）
From Single Cell to A Story: How Single-Cell Sequencing & AI Are Revolutionizing Biology - Bridge Informatics4.  sonraianalytics 人工智能驱动的单细胞分析和自动细胞注释
人工智能驱动的单细胞分析 - Sonrai Analytics5. 微流控技术的AI改进
Enhancing single-cell biology through advanced AI-powered microfluidics - PubMed6. developer nvidia  
    RAPIDS加速单细胞的多组学研究。

Accelerating Single Cell Genomic Analysis using RAPIDS | NVIDIA Technical BlogNVIDIA Parabricks v4.3 的加速功能
Boost Multi-Omics Analysis with GPU-Acceleration and Generative AI | NVIDIA Technical Blog 总结

最后请允许我借用Nature Biotechnology的编辑总结：

生物学的悄然革命正在进行中，因为单细胞工具（+AI）推动了下一波药物发现，并有望将疗法与个体相匹配。
Can single-cell biology realize the promise of precision medicine? | Nature Biotechnology很喜欢quiet revolution这个词，在“寒武纪式”的行业周期里，有一种“静水流深，至生命之本”的美感。

卡卡

作者按
本章节主要讲解了单细胞数据的细胞组成（比例）的分析方法，讲解了为什么直接分析比例不行的原因与例子，并介绍了“有标签”，“有标签有层次”，“无标签”三个维度的细胞组成分析方法。本教程首发于单细胞最好的中文教程，未经授权许可，禁止转载。
全文字数|预计阅读时间: 5000|10min
——Starlitnightly（星夜）

1. 背景

在单细胞分析中，我们除了关注基因表达模式受不同条件所影响导致的改变，我们还会关注细胞组成（例如细胞类型的比例）也会在不同条件下发生变化。例如药物处理，外源感染，细胞癌变等，这些刺激因素将会导致细胞类型的变化。由于细胞组成是一种整体性的变化，因此我们需要大量的细胞以及样本数量，才能证实组成的变化发生。



细胞组成变化

在本章中，我们将详细探讨细胞组成的变化的分析手段，并阐明清楚为什么不能直接对细胞比例进行分析？分析细胞比例的局限性在哪？这会是很有意思的一章教程。
import omicverse as ov
import scanpy as sc
import scvelo as scv
import pertpy as pt

ov.utils.ov_plot_set()

____            _     _    __                  
  / __ \____ ___  (_)___| |  / /__  _____________
/ / / / __ `__ \/ / ___/ | / / _ \/ ___/ ___/ _ \
/ /_/ / / / / / / / /__ | |/ /  __/ /  (__  )  __/
\____/_/ /_/ /_/_/\___/ |___/\___/_/  /____/\___/                                             

Version: 1.5.6, Tutorials: https://omicverse.readthedocs.io/2.加载数据

本章将介绍这两种方法并将其应用于 Haber 数据集[ Haber等人。，2017 ]。该数据集包含来自小鼠小肠和类器官的 53,193 个单个上皮细胞。一些细胞还受到细菌或蠕虫感染，例如分别通过沙门氏菌和Heligmosomoides polygyrus感染。在本教程中，我们使用完整 Haber 数据集的子集，其中仅包括专门为此目的收集的对照细胞和受感染细胞。值得注意的是，我们排除了仅收集大单元格的附加数据集，以加快计算速度并降低复杂性。
数据下载地址：https://figshare.com/ndownloader/files/38169900
adata=ov.read('data/haber_count.h5ad')
adata
<hr/>adata.obs['condition'].cat.categories
<hr/>Index(['Control', 'Hpoly.Day3', 'Hpoly.Day10', 'Salmonella'], dtype='object')数据分 10 个批次进行收集。条件变量是对照组（Control）、沙门氏菌（Salmonella）、肠道寄生线虫感染3天（Hpoly.Day3）和肠道寄生线虫感染10天（Hpoly.Day10）。cell_label存放了细胞类型。
3. 组成差异

我们在分析单细胞数据中的细胞比例的时候，通常是描述细胞在不同分组中的占比来进行叙述的，在一般的分析教程中，例如药物治疗后，某类T细胞的比例上升。但是，这类T细胞真的上升了吗？这是一个值得去思考的问题。我们来看下面这个例子：
在一项癌症研究中，我们想要比较健康器官与患病器官的细胞类型组成。我们假设在这种器官中存在三类细胞，分别是A，B，C，其中：

• 健康器官中，三类细胞的细胞数量相同，都是2000个
  
• 患病器官中，A类细胞变成了4000个，但是B类和C类细胞仍是2000个细胞
  

import pandas as pd
import numpy as np
healthy_tissue = [2000, 2000, 2000]
diseased_tissue = [4000, 2000, 2000]
example_data_global = pd.DataFrame(
    data=np.array([healthy_tissue, diseased_tissue]),
    index=[1, 2],
    columns=["A", "B", "C"],
)
example_data_global["Disease status"] = ["Healthy", "Diseased"]
example_data_global
<hr/>

	A	B	C	Disease status
1	2000	2000	2000	Healthy
2	4000	2000	2000	Diseased

import seaborn as sns
import matplotlib.pyplot as plt
plot_data_global = example_data_global.melt(
    "Disease status", ["A", "B", "C"], "Cell type", "count"
)

fig, ax = plt.subplots(1, 2, figsize=(12, 6))
sns.barplot(
    data=plot_data_global, x="Disease status", y="count", hue="Cell type", ax=ax[0]
)
ax[0].set_title("Global abundances, by status")

sns.barplot(
    data=plot_data_global, x="Cell type", y="count", hue="Disease status", ax=ax[1]
)
ax[1].set_title("Global abundances, by cell type")

plt.show()
<hr/>​



模拟真实数据

​
对于分析的目的而言，我们在获得健康器官和患病器官后，想探究清楚，患病后哪一类细胞的比例会变化。很显而易见，A类细胞翻倍变化了。但是我们在分析的时候，真的如我们理想情况一般吗？
答案是否定的，我们都知道，单细胞测序技术有着通量限制，即使是技术迭代发展，测序仪单次所能通过的细胞也在5000-6000之间。这意味着，我们每个病人器官所获的的细胞总数是有限的。我们无法获得器官内的所有细胞，这个时候。我们会用所测序的细胞来作为这个“病人器官”的代表。一个新的问题也随之出现，我们所获的的细胞类型计数只是一种比例，其总和为1。在统计学上，我们把这种数据称之为“成分数据”，特征是一个样本中的所有特征（细胞类型）的相对丰度总和为1.
因此，我们前面的A、B、C三类细胞，我们并不能获得6000和8000个细胞，我们假设测序仪的限制为600，因此我们只能从中随机获得600个细胞。我们采用随机抽样模拟这一过程。
np.random.seed(112)
healthy_sample = np.random.multinomial(
    pvals=healthy_tissue / np.sum(healthy_tissue), n=600
)
diseased_sample = np.random.multinomial(
    pvals=diseased_tissue / np.sum(diseased_tissue), n=600
)
example_data_sample = pd.DataFrame(
    data=np.array([healthy_sample, diseased_sample]),
    index=[1, 2],
    columns=["A", "B", "C"],
)
example_data_sample["Disease status"] = ["Healthy", "Diseased"]
example_data_sample
<hr/>

	A	B	C	Disease status
1	193	201	206	Healthy
2	296	146	158	Diseased

<hr/>plot_data_sample = example_data_sample.melt(
    "Disease status", ["A", "B", "C"], "Cell type", "count"
)

fig, ax = plt.subplots(1, 2, figsize=(12, 6))
sns.barplot(
    data=plot_data_sample, x="Disease status", y="count", hue="Cell type", ax=ax[0]
)
ax[0].set_title("Sampled abundances, by status")

sns.barplot(
    data=plot_data_sample, x="Cell type", y="count", hue="Disease status", ax=ax[1]
)
ax[1].set_title("Sampled abundances, by cell type")
plt.show()
<hr/>​



模拟真实数据采样

​
我们可以很直观地感受到，虽然A类细胞的比例依然是最大的，但是B、C两类细胞似乎在病人器官中出现了下降的情况，但对应我们的真实情况而言，B、C两类细胞在病人与健康人中的比例应该是不变的。
读到这里，你可能会去思考，如何才能将B、C两类细胞的误差进行消除，使得我们分析能得到一个合理的结果。
可能你注意到，例如我们同时对B、C两类细胞在病人中的样本减少1000个细胞，A类细胞保持2000个那么我们在采样600个细胞的时候，所获得的测序结果与A类细胞变为4000个是一致的。因此，我们可以考虑固定数据的参考点，我们假设该参考点在所有样本中保持不变。例如我们保持B、C细胞不变，那么A类细胞的分析结果将与真实情况一致。参考点的选取可以是单一细胞类型、多种细胞类型的聚合（例如几何平均值）或一组正交基。
4. 基于已有的标记

scCODA属于一种通过预定义类（例如细胞类型）来统计得出细胞的成分变化。该算法受微生物组数据分析的启发。scCODA 提出了一种贝叶斯方法来解决单细胞分析中常见的低重复问题。它使用分层狄利克雷多项式模型对细胞类型计数进行建模，该模型通过对所有测量的细胞类型比例进行联合建模来解释细胞类型比例的不确定性和负相关偏差。为了确保唯一可识别的解决方案和易于解释，scCODA 中的参考被选择为特定的细胞类型。因此，scCODA 检测到的任何成分变化始终必须相对于所选参考进行查看。
但需要注意的是，scCODA 假设协变量和细胞丰度之间存在对数线性关系，这在使用连续协变量时可能并不总是反映潜在的生物过程。scCODA 的另一个限制是无法推断除成分效应之外的细胞成分之间的相关结构。此外，scCODA 仅模拟平均丰度的变化，但不检测响应变异性的变化
4.1 初始化scCODA模型。

我们使用 load 函数准备 MuData 对象以供后续处理，并根据输入数据创建成分分析数据集。我们将 cell_type_identifier 指定为cell_label，sample_identifier 指定为batch，covariate_obs 与condition我们的例子一样。
sccoda_model = pt.tl.Sccoda()
sccoda_data = sccoda_model.load(
    adata,
    type="cell_level",
    generate_sample_level=True,
    cell_type_identifier="cell_label",
    sample_identifier="batch",
    covariate_obs=["condition"],
)
sccoda_data
<hr/>MuData object with n_obs × n_vars = 9852 × 15223
  2 modalities
    rna:        9842 x 15215
      obs:        'batch', 'barcode', 'condition', 'cell_label'
    coda:        10 x 8
      obs:        'condition'
      var:        'n_cells'<hr/>为了概述不同条件下的细胞类型分布，我们可以使用 scCODA 的boxplots. 为了更好地了解数据的分布方式，红点显示实际的数据点。
pt.pl.coda.boxplots(
    sccoda_data,
    modality_key="coda",
    feature_name="condition",
    figsize=(12, 5),
    add_dots=True,
    args_swarmplot={"palette": ["red"]},
)
plt.show()
<hr/>​



细胞比例箱线图

​
箱线图突出显示了细胞类型分布的一些差异。明显值得注意的是沙门氏菌病症的肠上皮细胞（Enterocyte）比例很高。但其他细胞类型，例如转运扩增（TA）细胞，与对照组相比，沙门氏菌条件下的丰度也表现出明显差异，我们必须正确评估这些差异是否具有统计显着性。
当然，你除了通过箱线图，你也可以通过堆叠柱状图来可视化细胞类型的比例变化。
pt.pl.coda.stacked_barplot(
    sccoda_data, modality_key="coda", feature_name="condition", figsize=(2, 2)
)
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
plt.ylabel('Prop',fontsize=11)
plt.legend(bbox_to_anchor=(1.05,1),fontsize=10)
plt.show()
​



细胞比例堆叠柱状图

​
4.2 scCODA模型训练

我们在初始化scCODA模型后，还需要指定两个主要参数：

• 公式（formula）：斜变量（格式为R语言的样式）
  
• 参考细胞类型（reference_cell_type）：用作参考的细胞类型
  

在上述例子中，我们将条件（condition）指定为唯一的斜变量，该变量包括来4种类型，scCODA将会对这四种类型进行分别建模。此外，我们发现内分泌细胞（Endocrine）在4种类型中的变化很小，即在所有样品中具有几乎恒定的相对丰度。如果参考细胞类型不清楚，我们可以令scCODA自身去计算合适的参考细胞类型，这时候，公式变量的输入如下
formula = "covariate_1 + covariate_2"reference_cell_type"automatic"
sccoda_data = sccoda_model.prepare(
    sccoda_data,
    modality_key="coda",
    formula="condition",
    reference_cell_type="Endocrine",
)
sccoda_model.run_nuts(sccoda_data, modality_key="coda", rng_key=112)
<hr/>sample: 100%|██████████| 11000/11000 [00:53<00:00, 207.54it/s, 127 steps of size 2.16e-02. acc. prob=0.76]
sccoda_data["coda"].varm
<hr/>AxisArrays with keys: intercept_df, effect_df_condition[T.Hpoly.Day3], effect_df_condition[T.Hpoly.Day10], effect_df_condition[T.Salmonella]计算结果如下表所示，我们选择effect_df_condition[T.Salmonella]进行解释，这代表了在沙门氏菌的组别的细胞比例变化情况，其中各参数的意义如下：

• Final Parameter: 统计学模型上的参数，类似y=kx+b中的k
  
• 预期样本（Expected Sample）与对数变化（log2-fold change）：我们首先假设当前协变量（T.Salmonella）的值为1，其他协变量的值为0，然后计算该预期样本与截距部分中没有活动协变量的预期样本之间的对数倍数变化。由于我们的数据是成分数据，如果我们的协变量的效应为0，那么预期样本与截距样本相同，那么log2的倍数变化也为0
  

在沙门氏菌病例中，我们仅看到肠上皮细胞的可信增加，而所有其他细胞类型均未受到该疾病的影响。具有相同总细胞计数的对照样品和感染样品之间肠细胞的对数倍数变化约为 1.54。
sccoda_data["coda"].varm["effect_df_condition[T.Salmonella]"]
<hr/>

	Final Parameter	HDI 3%	HDI 97%	SD	Inclusion probability	Expected Sample	log2-fold change
Cell Type
Endocrine	0.000000	0.000	0.000	0.000	0.0000	32.610693	-0.527456
Enterocyte	1.545527	0.994	2.066	0.283	1.0000	383.050645	1.702268
Enterocyte.Progenitor	0.000000	-0.463	0.526	0.143	0.2900	126.297502	-0.527456
Goblet	0.000000	-0.337	1.010	0.289	0.4288	52.701307	-0.527456
Stem	0.000000	-0.725	0.333	0.180	0.3205	135.184465	-0.527456
TA	0.000000	-0.913	0.326	0.223	0.3546	79.015201	-0.527456
TA.Early	0.000000	-0.352	0.638	0.158	0.3119	152.267715	-0.527456
Tuft	0.000000	-1.272	0.597	0.359	0.4152	23.072473	-0.527456

4.3 调整错误发现率

scCODA 根据包含概率选择可信效应。可信效应和不可信效应之间的界限取决于所需的错误发现率 (FDR)。较小的 FDR 值会产生更保守的结果，但可能会遗漏一些效应，而较大的 FDR 值会选择更多的效应，但会产生更多的错误发现。
默认情况下，该值为 0.05。由于根据数据集的不同，FDR 可能会对结果产生重大影响，因此我们建议尝试高达 0.2 的不同 FDR，以获得最显着的效果。
sccoda_model.set_fdr(sccoda_data, 0.4)
为了获得每种细胞类型的成分变化的二元分类，我们credible_effects在结果对象上使用 scCODA 函数。标记为“True”的每种细胞类型或多或少都存在。倍数变化描述了细胞类型是否更多或更少存在。因此，我们将把它们与下面的二元分类一起绘制。
sccoda_model.credible_effects(sccoda_data, modality_key="coda")
<hr/>Covariate                 Cell Type            
condition[T.Hpoly.Day3]   Endocrine                False
                          Enterocyte                True
                          Enterocyte.Progenitor    False
                          Goblet                    True
                          Stem                      True
                          TA                        True
                          TA.Early                 False
                          Tuft                      True
condition[T.Hpoly.Day10]  Endocrine                False
                          Enterocyte                True
                          Enterocyte.Progenitor    False
                          Goblet                    True
                          Stem                     False
                          TA                       False
                          TA.Early                  True
                          Tuft                      True
condition[T.Salmonella]   Endocrine                False
                          Enterocyte                True
                          Enterocyte.Progenitor    False
                          Goblet                   False
                          Stem                     False
                          TA                       False
                          TA.Early                 False
                          Tuft                     False
Name: Final Parameter, dtype: bool
pt.pl.coda.effects_barplot(sccoda_data, modality_key="coda",
                           covariates="condition")
plt.show()
<hr/>​



细胞组成分析图

​
这些图很好地显示了条件对细胞类型的显着且可信的影响。这些效应在很大程度上与 Haber 论文中的发现一致，该论文使用非组合泊松回归模型得出了他们的发现：

• “感染沙门氏菌后，成熟肠上皮细胞（Entercyte）的频率大幅增加。” [哈伯等人。，2017年]
  
• “蠕虫感染后将导致杯状细胞（Goblet）和簇状细胞（Tuft）丰度增加。” [哈伯等人。，2017年]
  

熟悉原始出版物的读者可能想知道为什么 Haber 等人使用的模型。发现比 scCODA 更显着的效果，例如沙门氏菌感染时干细胞和转运扩增细胞的减少。这与scCODA的结果不完全一致，这是由于细胞比例数据是成分数据，一种细胞类型的相对丰度的增加，将导致其他细胞类型的相对丰度减少。所以如果只使用非组合泊松回归模型便会出现这种问题。
5. 基于已有的标记以及层次结构

除了每种细胞类型的丰度之外，典型的单细胞数据集还以基于树的分层排序的形式包含有关不同细胞相似性的信息。这些层次结构可以通过基因表达的聚类（通常用于发现属于同一细胞类型的细胞簇）自动确定，也可以通过生物信息层次结构（如细胞谱系）自动确定。 tascCODA是 scCODA 的扩展，它将分层信息和实验协变量数据集成到成分计数数据的生成模型中[ Ostner et al. ，2021 ]。这对于高分辨率的细胞图谱工作特别有益。
5.1 层次结构聚类

要使用 tascCODA，我们首先必须定义细胞类型的分层排序。一种可能的层次聚类使用八种细胞类型，并根据它们在 PCA 表示中的相似性（皮尔逊相关性）对它们进行排序sc.tl.dendrogram。由于这种结构在我们的数据中非常简单，因此不会给我们带来很多新的见解，因此我们希望有一个更复杂的聚类。最近获得此类簇的一种方法是schist包[ Morelli等人。，2021 ]，它使用嵌套随机块模型，以不同的分辨率级别对细胞群进行聚类。使用标准设置运行该方法需要一些时间（在我们的数据上约为 15 分钟），并为我们提供了将每个单元格分配给adata.obs. 首先，我们需要通过 PCA 嵌入定义单元格之间的距离度量
ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,)
ov.pp.scale(adata)
ov.pp.pca(adata,layer='scaled',n_pcs=50)
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=30,
               use_rep='scaled|original|X_pca', random_state=112)
# Calculate UMAP for visualization purposes
#sc.tl.umap(adata)
adata.obsm["X_mde"] = ov.utils.mde(adata.obsm["scaled|original|X_pca"])
ov.utils.cluster(adata,method='schist', random_seed=112)
adata

<hr/>Begin robust gene identification
After filtration, 15215/15215 genes are kept. Among 15215 genes, 15215 genes are robust.
End of robust gene identification.
Begin size normalization: shiftlog and HVGs selection pearson
normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
['Cck', 'Defa24', 'Fabp6', 'Gcg', 'Ghrl', 'Gip', 'Nts', 'Reg3b', 'Reg3g', 'Reg4', 'Sct', 'Spink4', 'Sst', 'Tff3', 'Zg16']
    finished (0:00:00)
extracting highly variable genes
--> added
    'highly_variable', boolean vector (adata.var)
    'highly_variable_rank', float vector (adata.var)
    'highly_variable_nbatches', int vector (adata.var)
    'highly_variable_intersection', boolean vector (adata.var)
    'means', float vector (adata.var)
    'variances', float vector (adata.var)
    'residual_variances', float vector (adata.var)
End of size normalization: shiftlog and HVGs selection pearson
computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:16)
<hr/>WARNING: Threading backend does not work with this version of graph-tool
Switching to loky backend
minimizing the nested Stochastic Block Model
        minimization step done (0:09:20)
        consensus step done (0:12:35)
    done (0:12:35)
    finished: and added
    'nsbm', the cluster labels (adata.obs, categorical) (0:12:35)
AnnData object with n_obs × n_vars = 9842 × 15215
    obs: 'batch', 'barcode', 'condition', 'cell_label', 'nsbm_level_0', 'nsbm_level_1', 'nsbm_level_2', 'nsbm_level_3', 'nsbm_level_4', 'nsbm_level_5', 'nsbm_level_6'
    var: 'robust', 'highly_variable_features'
    uns: 'scaled|original|pca_var_ratios', 'scaled|original|cum_sum_eigenvalues', 'neighbors', 'schist'
    obsm: 'scaled|original|X_pca', 'X_mde', 'CM_nsbm_level_0', 'CM_nsbm_level_1', 'CM_nsbm_level_2', 'CM_nsbm_level_3', 'CM_nsbm_level_4', 'CM_nsbm_level_5', 'CM_nsbm_level_6'
    varm: 'scaled|original|pca_loadings'
    layers: 'counts', 'scaled', 'lognorm'
    obsp: 'distances', 'connectivities'<hr/>ov.utils.embedding(adata,
                   basis='X_mde',
                    frameon='small',
                   color=["nsbm_level_1", "nsbm_level_2",
                           "nsbm_level_3","cell_label"],
                   ncols=2)
<hr/>​



细胞层次簇

​
adata.write_h5ad('tmp/tasccoda_schist.h5ad',compression='gzip')
5.2 初始化TAsccoda模型

对于TAsccoda模型而言，我们需要在树的层面上作为样本的输入，即从cell_level变成cluster_level，我们希望簇的规模能尽可能大一些，而不是过小，所以我们从nsbm_level_1开始。
tasccoda_model = pt.tl.Tasccoda()
tasccoda_data = tasccoda_model.load(
    adata,
    type="cell_level",
    cell_type_identifier="nsbm_level_1",
    sample_identifier="batch",
    covariate_obs=["condition"],
    levels_orig=["nsbm_level_4", "nsbm_level_3", "nsbm_level_2", "nsbm_level_1"],
    add_level_name=True,
)
tasccoda_data5.3 TAsccoda模型训练

tascCODA 中的模型设置和执行与 scCODA 类似，参考的自由参数和公式也相同。此外，我们可以通过参数phi和lambda_1参数来调整树聚合和模型选择pen_args（有关更多信息，请参阅[ Ostner et al. , 2021 ] ）。在这里，我们使用无偏设置phi=0和模型选择，该模型选择比默认的(lambda_1=1.7)稍微宽松一些。
值得注意的是，我们发现内分泌细胞组主要是簇18，那么我们可以用簇18来做reference_cell_type
tasccoda_model.prepare(
    tasccoda_data,
    modality_key="coda",
    reference_cell_type="18",
    formula="condition",
    pen_args={"phi": 0, "lambda_1": 3.5},
    tree_key="tree",
)
tasccoda_model.run_nuts(
    tasccoda_data, modality_key="coda", rng_key=112,
    num_samples=10000, num_warmup=1000
)
<hr/>Zero counts encountered in data! Added a pseudocount of 0.5.

sample: 100%|██████████| 11000/11000 [04:52<00:00, 37.61it/s, 63 steps of size 8.74e-02. acc. prob=0.91]我们发现，TAscCODA的acc达到了0.91，远高于scCODA的0.76，这表明TAscCODA有着更好的准确性与可解释性。
5.4 TAsccoda模型的解释性

tascCODA 的结果首先应该被解释为对树节点的影响。节点上的非零参数意味着该节点下所有细胞类型的聚合计数显着变化。我们可以轻松地将这三种疾病状态的树状图可视化。蓝色圆圈表示增加，红色圆圈表示减少
pt.pl.coda.draw_effects(
    tasccoda_data,
    modality_key="coda",
    tree="tree",
    covariate="condition[T.Salmonella]",
    show_leaf_effects=False,
    show_legend=False,
)
或者，对内部节点的影响也可以通过树转换到细胞类型级别，从而允许像 scCODA 一样计算对数倍数变化。为了可视化细胞类型的对数倍数变化，我们绘制了与 scCODA 相同的图。
pt.pl.coda.effects_barplot(tasccoda_data, modality_key="coda",
                           covariates="condition")
<hr/>​



细胞组成分析2

​
通过绘制 UMAP 嵌入上每个条件的效应大小，并将其与细胞类型分配进行比较，可以获得另一种富有洞察力的表示
kwargs = {"ncols": 3, "wspace": 0.25, "vcenter": 0, "vmax": 1.5, "vmin": -1.5}
tasccoda_data["rna"].obsm['X_umap']=tasccoda_data["rna"].obsm['X_mde']
pt.pl.coda.effects_umap(
    tasccoda_data,
    effect_name=[
        "effect_df_condition[T.Salmonella]",
        "effect_df_condition[T.Hpoly.Day3]",
        "effect_df_condition[T.Hpoly.Day10]",
    ],
    cluster_key="nsbm_level_1",
    **kwargs
)
ov.utils.embedding(
    tasccoda_data["rna"],
    basis='X_mde',
    color=["cell_label", "nsbm_level_1"],
    ncols=2, wspace=0.5,
    frameon='small'
)
<hr/>​




结果与 scCODA 的发现非常相似：

• 对于沙门氏菌感染，我们得到簇的聚集增加，大约代表细胞类型簇中的肠细胞。对于簇 12,18,7,9，这种增加甚至更强，如对叶水平的额外积极影响所示
  
• 对于蠕虫感染，3 天后我们没有得到可信的变化。10 天后，我们发现包含干细胞（Stem）和转运扩增细胞（TA）的细胞簇减少，肠上皮细胞（Enterocytes）和肠上皮细胞祖细胞（Enterocyte progenitors）的减少也不太明显，scCODA 也发现了这一点。
  

此时的发现与原始论文则较为一致。
6. 无明确的细胞类型标记（发育）

还有一种有意思的情况是发育过程中，细胞的状态演变及比例变化。此时细胞位于过渡态居多，如果我们继续使用scCODA或者是TAscCODA，我们根据明确的细胞注释作为参考是不太恰当的。在这里有一些根据KNN邻域图的方法，或许会对发育过程中细胞谱系的变化，有一定帮助。注意，探究发育过程中的细胞比例变化的意义是不大的。这是因为多种细胞会有共享状态。此时，我们会考虑，不同的药物处理组，是否会带来不同的发育轨迹，而不是细胞比例的相关分析。

• DA-seq 使用一系列 k 值，根据细胞邻域中两种生物状态的细胞的相对流行率，为每个细胞计算一个分数[ Zhao et al. ，2021 ]。这些分数用作逻辑分类器的输入，以预测每个细胞的生物状况。
  
• Milo 将细胞分配到 KNN 图上部分重叠的邻域，然后使用广义线性模型 (GLM) 对细胞计数进行差异丰度 (DA) 测试[ Dann等人，2017]。，2022 ]。
  
• MELD 使用基于图形的密度估计来计算在每种条件下观察每个细胞的相对似然估计[ Burkhardt等人。，2021 ]。
  

受限于章节内容长度，本小节内容不做详细展开，感兴趣的可以自行去使用Milo进行分析，教程地址：https://www.sc-best-practices.org/conditions/compositional.html
7. 总结

• 如果主要兴趣在于已知细胞类型或状态之间的组成变化，请使用 scCODA 或 tascCODA 来统计评估丰度变化。
  
• 如果数据没有明显聚类（例如在发育过程中），如果我们对可能出现在细胞类型之间的过渡状态或特定子集中的细胞丰度差异感兴趣，则应使用基于 KNN 的方法，例如 DA-Seq 或 MILO给定细胞类型的细胞。
  

8. 思考

• 从视觉上直观推断细胞比例变化是很困难的。为什么？
  
• 为什么需要将细胞类型丰度解释为比例而不是绝对计数？不这样做会有什么问题？
  
• 在什么情况下应该使用使用簇信息（例如细胞类型）的工具，在什么情况下应该使用不使用簇信息的工具？
  

队长是我

单细胞现在无外乎三部分：一是大家关注的科学问题，一般是一个课题最初的起点。目前首先在动物中，比如癌症，发育，进化，免疫等领域单细胞都是解析异质性的一大利器。当然植物的研究也在跟上---对解析度的要求越来越高---单细胞水平。因为实验科学更多是Hypothesis driven，所以有科学假设，能够在单细胞水平去验证假设，基本上所有科研工作者应该都愿意吧，所以该技术的发展在需求上在近年来应该还会是井喷状态。
第二是测序技术的开发，也就是回答你科学问题的手段。近几年单细胞井喷状态的一个很重要原因是商业化的测序技术，这样大家只要有经费，就可以测序，而不是仅仅几个实验室单独掌握这项技术。商业化的程序使技术平民化了和流行，现在华大也进军单细胞领域了，一些小公司未来也会不断进入这个领域。从测序本身来看，未来的趋势应该是1.多组学同时在单细胞水平表征，2.空间转录组 3.(超)高通量4.样本的通用性等。这些还在不断的开发，再测试和商业化。
第三就是测序数据的处理。获得的单细胞测序数据，维度高，数据稀疏，不同模态数据格式不一等，这给算法开发，特别是目前大数据时代，人工智能的介入提供了很好的机会。目前很多算法开发的组提供的分析工具更新换代速度很快，比如Seurat马上到4了，当然还有很多其他的分析工具。我记得之前有个asesome single cell统计，应该有100多个分析工具。希望这些工具流程的标准化，以及未来能有大一统的分析工具。
我是比较积极看待单细胞近几年的。目前该技术确实是井喷状态，但也不得不说，目前常规的一些单细胞文章确实也遇到瓶颈了，未来该技术也许也会变成pcr一样比较常规的技术。