文献分享-2023-scPlant：植物单细胞转录组数据分析的框架

Abstract

近年来，随着单细胞转录组测序技术在植物学研究中的广泛应用，我们能够更好地了解植物生长发育和对外界刺激的响应机制。然而，目前存在一些挑战：现有的单细胞转录组数据分析工具主要面向特定的任务设计，缺乏端到端的分析工具，并且不同工具之间的整合较为困难。在这里，我们提出了一个通用的框架scPlant，用于通过用户提供的输入数据探索植物单细胞图谱。scPlant pipeline实现了从基础数据处理到细胞类型注释与反卷积、轨迹推断、跨物种数据整合、细胞类型特异性基因调控网络构建等多种分析任务的功能。此外，多种可视化工具绑定在一个内置的Shiny程序中，能够动态探索单细胞转录组数据。

INTRODUCTION

（发展）文章首先回顾了拟南芥、水稻、玉米等对根系、花部分生组织、外部刺激响应等在单细胞分辨率下的工作，使剖析植物组织的复杂细胞结构成为可能。虽然单细胞转录组技术应用于植物上的技术难点已经被解决，许多植物的组织和器官可以进行常规单细胞分析，但是从这些一次实验中单个细胞成百上千的转录组数据分析获得新的生物学发现需要数据分析方面的专业技能和综合分析框架。（问题）然而，现有的分析工具是针对特定分析任务，对于没有生物信息学经验的生物学家来，说从繁多的分析工具中对新的数据建立分析管道仍是有挑战的。尽管绘制单细胞转录组图谱已有guidelines和best practices的提出，但因为缺乏植物的相关知识，还不能立即应用于植物的分析中。因此建立一个最新的workflow对植物图谱研究是十分有必要的。

（解决）在这里，我们描述了一个端到端的计算框架，称为scPlant，用于探索植物的单细胞转录组数据。scPlant的管线结合了常用分析工具和相关知识，从而允许用户执行各种分析任务。scPlant提供了基本单细胞转录组数据处理的核心模块，以及用于高级分析的多个附加模块，包括自动注释细胞类型（在可能的情况下），反卷积（对于bulk数据），轨迹推断，跨物种数据整合，细胞类型特异性基因调控网络构建。此外，scPlant还内置了由shiny搭建的多个可视化工具箱，可即时对数据进行探索。scPlant可在以下网址免费获得：

https://github.com/compbioNJU/scPlant。

RESULTS AND DISCUSSION

Overview of scPlant

scPlant整合了流行分析工具用于植物单细胞转录组数据分析，如Seurat，SCENIC等（下表）。

scPlant的输入仅需要一个或多个矩阵，这个矩阵可以从scRNAseq或snRNA-seq实验获得。scPlant的示意图如图1A所示。

scPlant框架由一个核心模块和几个高级附加模块组成（图1B）。

核心模块用于对单细胞转录组数据进行预处理，包括质量控制（例如，去除潜在噪声），变异校正（例如，批处理效应），标准化，降维，细胞聚类，差异表达分析和整合。使用核心模块之后（参见最后一节scPlant的实践指南），可以使用无监督聚类方法（如Louvain聚类算法），然后使用umap或tsne将高维表达矩阵嵌入低维空间。为了注释cluster，差异表达分析用于将感兴趣的cluster与所有其他细胞进行比较来找到每个cluster的标记基因集。因此，映射细胞cluster和鉴定marker基因集是在进行单细胞分析的首要任务。原则上，运行核心模块是所有高级模块的先决条件。在以下部分中，我们将演示如何使用已发布的scRNA-seq数据集使用scPlant的高级模块来获取生物学见解。

Case study 1: Cell-type annotation and deconvolution

通过上述核心模块，单细胞数据已经细胞进行分群并获得了相关的marker基因。每个cluster的细胞类型鉴定可以通过已知基因表达每个簇的细胞身份可以通过覆盖已知标记基因的表达共定位来确定。参考marker基因数据库(Chen et al., 2021; Jin et al., 2022)和单细胞图谱(Seyfferth et al., 2021)等资源对于细胞身份鉴定有非常重要的价值。在细胞类型注释方面，我们在scPlant中开发了可以基于现有marker基因数据库/GO条目/或已经建立的细胞图谱进行自动注释的功能（图2A）。

我们实现了几种基于哺乳动物领域的自动注释方案，包括SingleR、CellAssign、Garnett、scCATCH等。我们手动注释了两个拟南芥scRNA-seq数据集(PRJNA517021和PRJNA706435)，水稻根(Zhang et al., 2021a)作为参考细胞图谱（图2B）。我们使用scPlant自动注释工具预测来自不同物种(PRJNA507252 和 CRA008788) (Ryu et al., 2019; Zhang et al., 2023)的新scRNA-seq/snRNA-seq数据集中的细胞类型，并比较了和手动注释的结果(ground truth; Figure 2C)。我们评估了基于scPlant marker数据库或参考细胞图谱的不同预测策略产生了希望的预测(Figure 2D和Supplemental Figure 1)。例如使用不同物种的参考注释或marker基因数据库，为复杂基因组（六倍体小麦）实现了高度准确的自动注释。注意，鉴于目前植物标记基因数据库远未完成，仍需要手动注释。然而，随着不同植物的参考数据越来越多，我们的自动注释为植物注释提供了一种可选方案。

特定细胞类型的生物学作用和功能特异性与细胞发挥作用的组织环境密切相关。作为细胞类型注释的替代方法，使用细胞类型特异性marker基因集构建的回归表达矩阵可以从bulk-RNA-seq数据中推断细胞类型组成及其对组织中全局基因表达变化的贡献，为此，我们在scplant管线中实现了细胞类型反卷积功能。解卷积功能对于比较不同条件下（如正常条件与应激反应）产生的bulk数据中细胞类型组成的变化特别有用，这样可以用于在进一步的功能研究中鉴定关键细胞类型或cluster。

为了说明这个概念，我们使用了从拟南芥根（Zhang et al。，2019a）获得的单细胞数据来卷积不同胁迫条件下产生的bulk-RNA-seq数据。显示维管组织和伴胞的组成 对各种应激源的响应有相当程度的改变（补充图2），突出了这些细胞类型在应激反应中的重要性。

Case study 2: Functional annotation of cell types

细胞类型注释的最终目标之一是基于所识别的簇的基因特征，用有意义的生物标记来表征每个细胞簇。为此，我们最近开发的CellFunTopic软件包中为细胞类型或cluster功能注释提供了一个新的模块。使用拟南芥根数据（Zhang et al。，2019a）进行案例研究，我们首先使用基因集富集分析（GSEA）来评估每个cluster的明确生物过程或通路（称为功能条目）（图3A和3B）。功能条目的参考注释可以从GO或KEGG数据库中获取。然后，我们使用主题模型（topic modeling）将可靠生物学意义的功能条目根据其主题贡献映射到相关细胞簇。

topic modeling：从文档集合中自动地找出一系列的主题（topics），每个文档集内可能有多个主题；主题：由众多词汇的概率分布（distribution）组成；
常用模型：LDA, pLSA， pLSI等，是一种无监督的学习过程；
Input：An unorganized collection of documents;
Output: An organized collection of topics

通过这种方式，我们能够用特定的功能条目注释每个细胞簇（图3C-3E）。

例如，光合细胞显示出与“叶绿体”，“光系统”和“响应”，“辐射相应”相关的富集功能。内皮层对应'对水的响应'和'水通道活性'。此外，scPlant能够对不同细胞cluster的功能相似性或特异性进行网络可视化（图3F）。它揭示了未注释的簇（unknown标记的18簇）与皮层紧密相连，与观察到的这个未知的cluster在细胞图谱上定位于皮层cluster相近一致（图2B）。

Case study 3: Inference of developmental trajectories

计算建模可以用于推断不同细胞类型之间的谱系关系。这些模型可用于识别控制发育的关键调控基因和信号通路。在植物系统中，单细胞发育轨迹的推断已被用于研究特定组织和器官的形成，如叶和根。为了进行这种分析，scPlant结合了各种广泛使用的伪时间推理工具，如Monocle2（Qiu et al。，2017），Monocle3（Cao et al。，2019），CytoTRACE（Gulati et al。，2020），和SlingShot（Street et al。，2018）。我们使用拟南芥根的单细胞转录组学数据（Zhang et al。，2019a）推断发育轨迹， 细胞类型手动注释过程如上所述（图4A-4D）。SlingShot通过构建最小生成树连接cluster，然后拟合每个分支的曲线（图4A），CytoTRACE从scRNA-seq数据预测细胞的分化状态（图4B）。SlingShot和CytoTRACE都不需要先验信息来推断细胞轨迹。然而，CytoTRACE预测了数据中轨迹错误的root（木质部细胞）。因此，可以提供先验信息作为轨迹可能的起点。

Monocle2和Monocle3允许指定root（比如本例中的分生组织细胞），这可能有助于工具发现正确的轨迹（图4C和4D）。

在测试数据集中，CytoTRACE, Monocle2和Monocle3构建的细胞发育轨迹的排序各不相同（图4E）。

然而，这些方法预测的pseudotimes总体显示出正相关（图4F）。

因此，方法的选择应主要基于轨迹的预期拓扑和数据集的维度(Saelens et al., 2019)。

Case study 4: Comparative analysis of single-cell atlases across plant species

随着越来越多不同植物的的单细胞图谱绘制，比较跨物种单细胞的基因表达是很有吸引力的工作，这为研究同一细胞类型中的保守基因表达，细胞类型谱系分析，和推断细胞身份的进化起源提供了机会( Ryu et al , 2021)。对于这类分析，我们在scPlant框架中提供了one-to-one的同源基因作为锚点的工具，对匹配的器官/组织单细胞数据进行跨物种整合。值得注意的是，scPlant主要使用CCA 和RPCA策略进行数据整合(Stuart et al., 2019)。虽然这些方法已经在不同的应用中得到了广泛的使用，但在植物单细胞数据中测试其他先进的整合方法的性能(Luecken et al., 2021)是很有吸引力的。

作为一个案例研究，我们通过分析已发表的拟南芥( Zhang et al , 2019a)、水稻( Zhang et al , 2021a)和玉米( Ortiz-Ramèrez et al , 2021)(图5A ~ 5D)的scRNA - seq数据，试图在根中生成一个跨物种的整合细胞图谱。

我们利用已知的marker基因(图5B , 5E , 5F)对每个物种的单细胞图谱进行手动注释和整合。

我们发现整合前后大多数细胞身份是一致的，并且整合后物种特异性的细胞类型仍然保留(图5C )，说明了整合分析的鲁棒性。例如，分生组织细胞在每个物种( Zhang et al . , 2019a , 2021a ; Ortiz-Ramèrez et al , 2021)中都被单独注释(图5B )，这些细胞在整合的细胞图谱中聚在一起，表达具有高度保守性和特异性的标记基因(图5F和5G)。进一步研究在同一细胞类型中跨物种表达的保守或特有的基因，以及基因调控网络会是有趣的。

Case study 5: Construction of cell-type-specific gene regulatory networks

基因调控网络（GRN）由转录因子（TF）与其靶基因（TG）之间的相互作用组成，协调细胞特异性基因表达模式，进而决定细胞功能。为了使用单细胞数据构建细胞类型特异性GRN，我们将单细胞GRN分析工具SCENIC(Aibar et al., 2017)整合到scPlant管道中，并准备了支持SCENIC进行在植物数据分析所需的数据库。使用上述分析中注释的细胞图谱（Zhang et al。，2019a），我们可以使用SCENIC基于共表达分析和TF基序富集来鉴定细胞类型特异性GRN和regulons（图6）。我们随scPlant提供了各种可视化工具使用不同方法来展示单细胞GRN结果。

例如，跨cluster的regulon activity模式可以在热图（Gu，2022）中展示（图6A和6B）。高度特异性的regulons可以在散点图（图6C）和网络图中展示（图6D）。我们发现光敏色素相互作用TF PIF5是光合细胞的top regulon TF，在光合作用和光保护中起调节作用。

最近，开发了植物特定的GRN推测工具MINI-EX(Ferrari et al., 2022)。我也将此工具整合到我们的scPlant管道中。为了比较MINI-EX和SCENIC的结果，我们对它们的输出结果以类似的方式进行可视化（补充图3）。

总体而言，MINI-EX和SCENIC预测的regulon在不同细胞类型中是一致的（图6E和6F）。

例如，两种方法都将bZIP44鉴定为伴胞的特异性调节子TF，并且ANAC030被鉴定为木质部细胞的调节子TF(Pawittra et al., 2020)。然而，我们注意到，就regulons富集评分而言，MINI-EX产生的信噪比高于SCENIC（图6E）。这种差异可能反映了每种工具识别regulons的不同策略。MINI-EX为特定细胞类型分配regulons时使用TF和TG表达水平来过滤。作为对比，SCENIC根据regulon（TFs+TGs）的总体表达计算regulon活性评分，导致在特定细胞类型中有低TF表达值的regulon。

我们还适配了PMET(Rich-Griffin et al., 2020)来预测细胞类型特异性标记基因的启动子区域内的TF结合基序对。为了以细胞类型特异性方式检测regulon的共相关，在PMET分析中仅使用来自上述SCENIC分析的regulon TF。例如，图6G中显示了基序在特定细胞cluster（光合细胞簇23）启动子内一致基因的基序富集对。一个高度显著的TCP家族基序对被鉴定为光合作用细胞cluster特异，包括已知参与光响应的TCP14和TCP15(Viola et al., 2023)。此外，ABI5和ABF基序对也被富集，它们已被证明抑制光合作用，促进叶绿素分解代谢和叶片衰老( Collin et al , 2021)。

Practice guidelines

为了使用scPlant对单细胞转录组数据进行分析，我们提供了运行上述案例研究的真实代码(Box 1)。当前的protocol需要2 - 3天的时间来建立工作流程并对多个已发表的单细胞数据进行完整的分析(Zhang et al., 2019a, 2021a; Ryu et al., 2019; Ortiz-Ramı ́rez et al., 2021)。该pipeline可以很容易地适应于分析各种植物的新数据集。该workflow面向对计算机软件和Linux命令行环境基本熟悉的植物生物学家。提供了在线指导和实践的notebook可以帮助用户使用推荐的工作流程( https://compbionju.github.io/scPlant/)。

最后，我们设计了一个内置的Shiny应用程序( scPlantApp )，它结合了R编程语言的可视化能力和用户友好的网页界面。(Chen et al.,  2018)，可以对于上述workflow所得的分析结果进行交互探索(图7 )。这种功能有利于对感兴趣的cluster、基因和GRNs进行深度挖掘。

scPlant features

ScPlant是专门为植物scRNAseq数据分析而设计的平台。其一些特征让它对于植物scRNA-seq的分析十分适用：

（1）植物热定基因组注释：scPlant平台使用的marker基因、调控基序数据库和功能注释是植物特定的，这有助于在植物scRNA-seq数据中对细胞类型进行准确鉴定和功能注释。此外，封装的GRN推断工具MINIEX是专门为植物单细胞数据分析而设计的。

（2）植物特异性参数：sc Plant平台可以使用针对植物scRNA - seq数据分析而特定的参数。例如，在植物中，每个细胞中检测到的基因数低于人和动物。因此，在使用scPlant平台时，可以使用不太严格的cut off对数据进行质量控制，以避免过滤掉基因低表达的细胞。此外，植物特有的细胞器，如叶绿体和线粒体，具有独特的基因表达谱，可以用于质量控制。

( 3 )可视化和探索：scPlant平台包含交互式可视化工具，使用户能够以用户友好的方式对其sc RNA-seq数据进行探索和分析。这包括可视化基因表达模式、识别共表达基因模块及探索基因网络和通路的工具。