方法文章

使用 R、Seurat 和 CellChat 分析小鼠皮肤伤口愈合的单细胞转录组数据集

5.3K 次观看

⸱

DOI：

10.3791/67266

⸱

2025年8月1日

 ,  ,  , 
1Department of Oral Biology, College of Dentistry, University of Illinois Chicago

本文内容

摘要

本文以小鼠皮肤伤口愈合的单细胞时间序列转录组数据集为例，利用 R 语言展示了一套逐步、可视化的分析工作流程。该方案包括使用 Seurat 进行数据集下载、质量控制、可视化和细胞类型注释的标准流程，以及使用 CellChat 进行细胞间相互作用分析。

摘要

伤口愈合过程受到不同细胞类型在时空上复杂相互作用的调控。通过对复杂微环境中单个细胞进行分析，单细胞转录组学方法能够研究参与伤口愈合过程的细胞异质性、细胞间通讯网络以及细胞-细胞相互作用。然而，许多单细胞分析工具需要在计算机编程环境中运行，而伤口愈合领域科研人员普遍缺乏生物信息学专业知识，这限制了这些工具的广泛应用。因此，本文提供了一套逐步操作流程，展示如何使用名为 RStudio 的图形化编程环境，对小鼠皮肤全层切口伤口愈合的时间序列数据集进行基础的单细胞分析。这一可视化且具有引导性的实验方案将帮助无生物信息学背景的科研人员下载已发表的伤口愈合数据集，执行关键的质量控制步骤，利用 Seurat 完成标准的单细胞分析流程（包括数据集可视化和细胞类型注释），开展细胞亚群分析、模块评分分析，使用 CellChat 进行细胞间相互作用分析，并利用 Seurat 实现多个数据集的整合分析。本方案为每一步操作提供了详细的叙述性说明，并展示了每一行代码对应的图形化结果，以安全、准确地引导用户完成整个分析流程。本可视化单细胞分析流程的目的是使更多伤口愈合领域的科研人员能够在自己的实验室中直接应用生物信息学工具，从而更深入地分析自身产生的单细胞数据集，并广泛地对已发表的单细胞数据集进行再分析。

引言

伤口愈合是哺乳动物生物学中最复杂的过程之一，包含三个阶段：炎症期、增殖期和修复完成期1,2。这些愈合阶段大致概括了数十种细胞类型及其数百种分子产物在空间和时间上的协调作用，贯穿整个伤口修复过程3。过去几十年中，基于在愈合时间进程中采集伤口组织样本所开展的组织学和分子研究，已阐明了组织修复的整体细胞模式3，尤其是在可重复的皮肤缺损小鼠模型中4,5,6。直到最近二十年，随着高通量转录组分析技术的发展，人们才得以更全面地认识伤口愈合的复杂性，最初应用于伤口整体组织水平7,8,9，随后扩展至单细胞水平10,11,12,13,14。最近，多项研究在单细胞层面描绘了皮肤伤口的转录图谱，鉴定出新的伤口相关细胞亚型，并揭示了它们在愈合过程中可能的相互作用方式15,16,17,18,19,20。Hu 等人采用一种创新的空间单细胞 RNA 测序方法，对愈合过程中从伤口中心不同径向距离处的皮肤伤口进行了分析，揭示了细胞间及分子层面在时空上的新“动态”变化20。此类研究以前所未有的细节解析了伤口愈合的复杂性，开始呈现出高度的细胞与分子异质性图景。

近年来，生物信息学分析方法的重大进展使得研究人员能够深入理解伤口愈合研究领域中产生的复杂多组学数据集。诸如 Seurat 等单细胞分析工具包提供了对数据集进行稳健分析与整合的手段，包括对伤口等复杂组织中的细胞类型进行分类21。在单细胞数据的下游解释方面，研究人员常使用 CellChat 等工具来识别潜在的细胞间相互作用程序，从而揭示细胞如何协调以实现伤口修复22。尽管这些工具具有完善的文档记录并被广泛引用，但它们必须在计算机编程环境中运行，例如 R 语言——一种统计与图形化编程语言，广泛应用于基因组学和转录组学等生物信息学领域。尽管伤口愈合领域的生物学家和临床医生正越来越多地采用单细胞技术研究组织修复，但极少有人具备在各自实验室中直接使用 Seurat 和 CellChat 等工具所需的生物信息学训练背景。这种技术门槛不仅阻碍了科研人员在无需生物信息学专家协助的情况下深入分析自身数据集，也限制了他们对其他研究团队已发表的大量单细胞数据进行可靠地重复分析的能力。

因此，本文提供了一个逐步操作的工作流程，旨在帮助没有生物信息学背景的科研人员分析一个已发表且公开可用的单细胞伤口愈合数据集20。该方案采用广泛使用且免费的图形化 R 编程环境 RStudio，演示了如何在该环境中运行指定代码行，利用 Seurat 和 CellChat 对复杂的单细胞数据集进行基础分析。本方案涵盖了伤口愈合研究中相关的七项主要方法，包括：1）编程环境的安装，2）数据集的下载及关键的质量控制步骤，3）单细胞分析工作流程（包括可视化和细胞类型注释），4）细胞亚型分析，5）模块评分分析，6）细胞间相互作用分析，以及 7）多个数据集的整合分析。在每种方法中，均提供了可实际运行的代码，供用户在执行方案时同步操作，同时展示了每一行代码生成的实际图形结果，以引导用户完成整个工作流程。本指南性且可视化地介绍 RStudio 及基础单细胞分析流程的主要目标，是使更多伤口愈合领域的科研人员能够直接使用这些强大工具，从而推动该研究领域更快发展。

方案

注意：在以下详细描述七种生物信息学方法的工作流程中，所有实验步骤均附有相应的代码块，用户应按照所列顺序直接在其自身的 RStudio 界面中运行这些代码。为了使本实验方案尽可能便于使用，本文提供了一个 R 脚本文件（补充文件 1：JoVE_Rscript.R），可直接加载到用户的 RStudio 会话中，以便逐行运行代码。此方式可避免用户手动输入或复制粘贴协议文档中的代码，从而减少出错的可能性。所有实验方案说明也以注释形式包含在 R 脚本文件中，每条注释行以井号“#”开头。

1. 安装用于单细胞分析流程的 R、RStudio 及所需 R 软件包

  1. 下载并安装 R（版本 4.4.1）到计算机。使用与计算机操作系统对应的链接。
    1. 如果使用运行 Microsoft Windows 的计算机，请使用此链接： https://cran.rstudio.com/bin/windows/base/
    2. 如果使用运行 MacOS 的计算机，请使用此链接： https://cran.rstudio.com/bin/macosx/
  2. 在计算机上安装最新版本的 RStudio。点击以下链接并按照说明操作：
    https://www.rstudio.com/products/rstudio/download/#download
  3. 安装 Rtools（版本 4.4），这将允许 R 编译某些软件包。点击以下链接并按照说明操作：
    1. 如果使用 Windows，请使用此链接： https://cran.rstudio.com/bin/windows/Rtools/
    2. 如果使用 MacOS，请使用此链接：
      https://mac.r-project.org/tools/
  4. 设置本地工作目录；这是计算机中用于加载和保存所有文件的文件夹。通过选择以下方式设置工作目录： 会议 在 RStudio 菜单栏中并点击 设置工作目录 > 选择目录 并选择所需的文件夹。
    1. 如果使用 Windows 计算机，请使用以下命令设置工作目录。将以下代码行中的 [Directory] 更改为实际的目录结构。请注意，R 中的目录分隔符为字符 "/"
      setwd("C:/[目录]")
    2. 如果使用 MacOS 计算机，以下命令还将设置工作目录。请将以下代码行中的 [Directory] 更改为实际的目录结构。注意，R 中的目录分隔符为字符 "/"
      setwd("~/[目录]")
    3. 在 R 会话的任何时刻，均可使用以下代码行检查工作目录：
      getwd()
    4. 在 RStudio 中，通过右侧窗口直观浏览工作目录结构，包括其中包含的所有文件和文件夹 文件 选项卡。要导航至 RStudio 文件浏览器中的工作目录，单击 齿轮 图标 -> 转到工作目录.
  5. 从 R 软件包仓库 CRAN 安装以下软件包，这些是本实验方案所必需的依赖包。要安装这些软件包，请运行以下命令。
    install.packages("devtools")
    install.packages("readxl")
    install.packages("openxlsx")
    install.packages("tidyverse")
    install.packages("scCustomize")
    注意： 在安装 R 包期间，各种窗口的出现和消失是正常现象。如果出现要求编译包的窗口，请单击 是如果出现窗口提示在安装包之前需要重启 R，请单击 否.
  6. 从经过筛选的 R 软件包仓库 Bioconductor 安装以下软件包
    (https://bioconductor.org/），这些是本实验方案所必需的依赖包。要安装这些软件包，请运行以下命令。 
    if (! requireNamespace("BiocManager"，quietly = TRUE) )
      install.packages("BiocManager")
    BiocManager::install("NMF", update=F)
    BiocManager::install("ComplexHeatmap", update=F)
    BiocManager::install("BiocNeighbors", update=F)
    BiocManager::install("SingleCellExperiment", update=F)
    BiocManager::install("circlize", update=F)
    BiocManager::install("edgeR", update=F)
    BiocManager::install("scDblFinder", update=F)
  7. 安装本研究中所述工作流程所需的以下软件包。
    install.packages("Seurat")
    devtools::install_github("jinworks/CellChat")
  8. 加载每个包以确认安装成功。如果任何包出现 "包未找到" 错误，请使用上方的相应代码重新安装。
    library(readxl)
    library(openxlsx)
    library(tidyverse)
    library(scCustomize)
    library(edgeR)
    library(scDblFinder)
    library(Seurat)
    library(CellChat)

2. 加载单细胞伤口愈合数据集并执行质量控制步骤

注意：本生物信息学工作流程对先前发表的一项时空单细胞皮肤伤口愈合实验进行了重新分析20。数据集文件存储在经过 curated 的美国国家生物技术信息中心（NCBI）基因表达综合数据库（GEO）仓库中（https://www.ncbi.nlm.nih.gov/geo/）。

  1. 使用登录号 GSE204777 导航至 GEO 上的数据集文件。使用以下链接并查看该研究的实验设计：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204777
  2. 在 GEO 仓库页面中，有来自五个测序通道的五个独立数据批次。点击第一个数据集，标题为 GSM6190913。以下是该样本的直接链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM6190913
  3. 向下滚动至页面底部，使用 ftp 或 html 链接下载以下三个文件。在计算机的文件浏览器中，将这三个文件移至名为 b1 的目录中。确保 b1 文件夹位于步骤 1.4 中设置的工作目录内。
    文件名：GSM6190913_b1_barcodes.tsv.gz / 文件大小：18.5 Mb
    文件名：GSM6190913_b1_features.tsv.gz / 文件大小：254.1 Kb
    文件名：GSM6190913_b1_matrix.mtx.gz / 文件大小：151.2 Mb
  4. 获取在步骤 2.3 中下载的单细胞测序文件的目录信息。
    b1_data_dir <- file.path(getwd(), "b1")
  5. 加载单细胞测序文件。gene.column 参数指定所使用的基因/特征命名方式。此处使用 gene.column = 2 表示基因符号（gene.column = 1 用于 Ensembl 基因名称）。
    b1 <- Read10X_GEO(data_dir = b1_data_dir, gene.column = 2)
    ​注意：大多数单细胞数据集不包含额外的多重检测数据，因此使用此步骤生成的 10x 文件将不具有多个图层。对于当前的多重检测数据集，请继续执行步骤 2.6。对于无多重检测数据的数据集，请跳至步骤 2.7。
  6. 使用时空条形码对单细胞数据集进行去多重化。
    1. 对于工作数据集，分离基因表达和 HTO（多重检测）数据。
      dataset_rna_counts <- b1$GSM6190913_b1_$`Gene Expression`
      ​dataset_barcodes <- b1$GSM6190913_b1_$`Antibody Capture`
    2. 使用基因表达数据创建 Seurat 对象，同时立即过滤掉在少于 5 个细胞中表达的基因以及检测到少于 200 个基因的细胞。
      dataset <- CreateSeuratObject(counts = dataset_rna_counts, min.cells = 5, min.features = 200)
    3. 将基因表达数据创建为一个图层，并生成两个检测共有的细胞和条形码列表。
      dataset_rna_counts2 <- LayerData(dataset, data = "RNA")
      ​dataset_joint.barcodes <- intersect(colnames(dataset_rna_counts2), colnames(dataset_barcodes))
    4. 根据共有的细胞条形码对基因表达和 HTO 计数进行子集筛选。
      dataset_rna_counts2 <- dataset_rna_counts2[, dataset_joint.barcodes]
      ​dataset_barcodes2 <- as.matrix(dataset_barcodes[, dataset_joint.barcodes])
    5. 确认 HTO 具有预期的条形码名称。
      ​rownames(dataset_barcodes2)
    6. 创建一个新检测以存储条形码信息，并将此检测添加到先前创建的 Seurat 对象中。
      dataset_barcode_assay <- CreateAssayObject(counts = dataset_barcodes2)
      ​dataset[["barcodes"]] <- dataset_barcode_assay
    7. 验证对象现在是否包含多个检测。
      DefaultAssay(dataset)
    8. 对条形码数据进行归一化，并通过 HTODemux 函数执行去多重化。该方法在以下 Seurat 教程中有详细描述： https://satijalab.org/seurat/articles/hashing_vignette
      ​dataset <- HTODemux(dataset, assay = "barcodes", positive.quantile = 0.99)
    9. 根据全局分类结果对细胞进行分组，并移除无条形码分类的细胞。
      Idents(dataset) <- "barcodes_classification.global"
      ​dataset <- subset(dataset, idents = "Negative", invert = TRUE)
    10. 根据最大 HTO 信号对细胞进行分组。
      Idents(dataset) <- "barcodes_maxID"
    11. 可视化每个多重条形码对应的细胞中检测到的基因分布（补充图 1）。
      VlnPlot(dataset, features = "nFeature_RNA", pt.size = 0.1, log = TRUE)
    12. 将条形码重命名为其实际的伤口时间（伤后天数）和空间（2-8 mm）分配（取自原始论文），并将其赋值给一个名为 time_space 的新元数据变量。
      Idents(dataset) <- "barcodes_maxID"
      levels(dataset)
      dataset <- RenameIdents(dataset,
      "Barcode-1" ="D01_2mm",
      "Barcode-2" ="D01_4mm",
      "Barcode-3" ="D01_6mm",
      "Barcode-4" ="D01_8mm",
      "Barcode-5" ="D03_2mm",
      "Barcode-6" ="D03_4mm",
      "Barcode-7" ="D03_6mm",
      "Barcode-8" ="D03_8mm",
      "Barcode-9" ="D07_2mm",
      "Barcode-10" ="D07_4mm",
      "Barcode-11" ="D07_6mm",
      "Barcode-12" ="D07_8mm",
      "Barcode-13" ="D14_2mm",
      "Barcode-14" ="D14_4mm",
      "Barcode-15" ="D14_6mm",
      "Barcode-16" ="D14_8mm",
      "Barcode-17" ="UW")
      levels(dataset)
      dataset[["time_space"]] <- Idents(dataset)
  7. 对于无多重检测数据的数据集：创建 Seurat 对象，同时立即过滤掉在少于 5 个细胞中表达的基因以及检测到少于 200 个基因的细胞。
    dataset <- CreateSeuratObject(counts = b1, min.cells = 5, min.features = 200)
  8. 切换至分析数据集的基因表达检测。
    DefaultAssay(dataset) <- "RNA"
    Idents(dataset) <- "data"
  9. 作为一项重要的质控步骤，计算每个细胞中线粒体基因的百分比，并将其赋值为一个元数据变量。该方法在以下 Seurat 教程中有详细描述： https://satijalab.org/seurat/articles/pbmc3k_tutorial#qc-and-selecting-cells-for-further-analysis
    dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^mt-")
  10. 可视化所有细胞中检测到的基因数量、RNA 数量及线粒体基因百分比的分布（补充图 2）。
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
    plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
    plot1 + plot2
  11. 存在一些线粒体含量较高的细胞，这与低 RNA 计数相关；这些是死亡或濒死细胞。使用合理的阈值从数据集中移除这些低质量细胞。此处采用先前发表的研究20中描述的原始数据集的值，即移除线粒体基因比例超过 25% 的细胞。
    dataset <- subset(dataset, subset = nFeature_RNA > 200 & percent.mt < 25)
  12. 在移除低质量细胞后，可视化所有细胞中检测到的基因数量、RNA 数量及线粒体基因百分比的分布（补充图 3）。
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
    plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
    plot1 + plot2
  13. 作为另一项重要的质控步骤，检测数据集中可能的双细胞（doublets）。这些是在液滴测序过程中合并的细胞，因此其基因表达并非真正的单细胞水平。为解决此问题，已开发出多种工具。需注意的是，每种工具对单细胞数据均有一定假设，因此用户在使用任一工具前应仔细阅读所有相关文档。在此工作流程中，采用一种名为 scDblFinder23 的方法。请注意，该方法使用一种设定固定预期双细胞率的算法。使用以下命令运行 scDblFinder 流程。
    DefaultAssay(dataset) <- "RNA"
    sce <- scDblFinder(GetAssayData(dataset, assay="RNA", slot="counts"), samples=Idents(dataset))
    注意：此类多重单细胞数据集也可通过移除表达多个条形码的细胞来筛选双细胞。本工作流程未展示此方法，因为大多数单细胞数据集不具备此独特特征。相反，此处展示了使用 scDblFinder 方法进行更具普适性的双细胞筛选流程。
  14. 将双细胞评分赋值给一个新元数据变量。
    dataset$scDblFinder.score <- sce$scDblFinder.score
  15. 可视化所有细胞中双细胞评分的分布（补充图 4）。
    VlnPlot(dataset, features = "scDblFinder.score", raster=FALSE, pt.size=0.5)
  16. 移除双细胞评分高于 0.25 阈值的细胞。该阈值基于上述生成的小提琴图选择，图中显示数据集中大多数细胞的双细胞评分均处于极高或极低水平，0.25 是该数据集的一个合理截断值，可移除绝大多数可能的双细胞，同时避免移除大量非双细胞。
    dataset <- subset(dataset, scDblFinder.score < 0.25)
  17. 将数据集 Seurat 对象保存为工作目录中的 RDS 文件。
    saveRDS(dataset, "dataset_post_Method2.rds")

3. 使用 Seurat 分析单细胞伤口愈合数据集

注意：（可选步骤）如果从此处开始工作流程，请将已保存的 RDS 文件作为 Seurat 对象加载。

dataset <- readRDS("dataset_post_Method2.rds")

  1. 执行 Seurat 标准工作流程，对单细胞数据集进行归一化、缩放及主成分分析（PCA）。该标准工作流程详见以下 Seurat 示例文档：
    Seurat - 受指导的聚类分析教程： https://satijalab.org/seurat/articles/pbmc3k_tutorial
    Seurat 命令列表： https://satijalab.org/seurat/articles/essential_commands
    DefaultAssay(dataset) <- "RNA"
    数据集 <- NormalizeData(object = dataset)
    数据集 <- FindVariableFeatures(object = dataset)
    数据集 <- ScaleData(object = dataset)
    数据集 <- RunPCA(object = dataset)
  2. 通过前 50 个 PCA 维度可视化数据集变异程度补充图5).
    ElbowPlot(数据集, reduction = "pca"，ndims = 50
    大部分主要变异出现在前13个维度中。
  3. 使用设定的主成分分析（PCA）维度范围1-13和设定的分辨率为0.1，对该数据集进行细胞聚类。
    数据集 <- FindNeighbors（数据集，verbose = TRUE，dims = 1:13）
    数据集 <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
    注意： 该工作流程侧重于细胞类型的大规模差异，因此在主成分分析（PCA）的降维范围中采用了一个相对保守的参数，前13个维度能够代表数据集中的绝大部分变异。为了将细胞进一步区分为更小且更稀有的亚型，用户可在下游分析中使用更高数量的维度，因为这些稀有细胞亚型可能仅贡献数据集变异的较低比例。分辨率参数的取值范围为0至1，用于控制施加于数据集上的类别分离程度。该参数的设定取决于用户的具体研究问题。若需将细胞聚类为多个较小且稀有的亚型，应在下游分析中使用较高的分辨率值。由于本工作流程旨在探索主要细胞类型之间的广泛差异，因此采用了一个较小的分辨率值0.1，预期将细胞聚类为数量较少但规模较大的群组。使用第3.3步中所述的参数设置可区分出8个独特的细胞簇。这些细胞簇将被自动分配至一个名为 "seurat_clusters".
  4. 使用前13个主成分分析（PCA）维度进行UMAP降维和寻找邻近点分析。添加随机种子数123，以确保数据投影结果的可重复性。
    数据集 <- 运行UMAP(dataset, verbose = TRUE, dims = 1:13, seed.use = 123)
    注意： UMAP 算法具有随机性，会将随机性引入降维过程（参见 "稳定性和可重复性" 在 https://cran.r-project.org/web/packages/umap/vignettes/umap.html）。在使用一致的种子值时，可提供一种可重复的方法 "最低程度的可重复性" 由于算法的随机性，生成的图表可能与代表性图示及下游结果略有差异。测试发现，使用 Windows 系统和 MacOS 系统的计算机所得结果可能存在显著差异，这可能是由于这些操作系统对随机性的实现方式不同所致。
  5. 在 UMAP 图上可视化细胞的聚类情况（图1).
    DimPlot(数据集, group.by = "seurat_clusters"，raster = FALSE，label = TRUE）
    注意： UMAP算法的随机性可能导致生成略微不同的图，如图所示，在Windows和MacOS系统上使用相同代码生成的备用UMAP图中可见；请注意各簇形状的细微差异。因此，用户必须在生成数据和图表的同时进行保存并添加时间戳，并且在对簇进行后续分析时应谨慎操作，结合生物学背景进行理解，如下文细胞类型注释所述。
  6. 由于原始实验标签包含了细胞在伤口愈合过程中来源的位置和时间信息，可在UMAP图上可视化细胞的伤口时间/空间注释（图2).
    DimPlot(数据集, group.by = "时间_空间"，raster = FALSE，label = FALSE）
  7. 生成细胞簇与伤口时间/空间注释的关联表格。
    table(dataset$time_space, dataset$seurat_clusters)
  8. 确定数据集中主要细胞类型的分类身份。为此，计算所有簇之间的差异表达基因（DEG）。获取各簇的DEG列表，将其赋值给一个变量，并将输出结果保存为工作目录中的分隔文本文件。
    注意： 此步骤对 CPU 资源消耗较大，所需时间长短取决于用户的硬件配置。
    Idents(dataset) <- "seurat_clusters"
    细胞标志物 <- FindAllMarkers(数据集, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
    write.csv(Cell_markers, file = file.path(getwd(), "dataset_cluster_markers.txt"))
  9. 下载并打开所包含的 dataset_cluster_markers.txt 在电子表格（例如 Excel）中通过复制文本文件的内容并使用该 文本导入向导 指定逗号分隔符以及基因名称列的身份 文本. 表明基因名称为‘文本’ 很重要，否则 Excel 会自动将某些基因名称转换为日期，例如 Sept7 会变为 9月7日。
  10. 在电子表格中，根据以下推荐参数筛选结果：
    1. 排序 avg_log2FC 按从大到小的顺序排列该列，以使所有行按 log 值递减排序2 倍数变化（log2FC）。
    2. 排序 簇 将列按从小到大的顺序排列，以根据 Seurat 聚类编号递增的方式重新排列所有行。
    3. 筛选 avg_log2FC 用于显示大于或等于 2.5 的数值列，以展示指定簇与其他簇相比差异表达最显著的基因（DEGs）。
    4. 筛选 pct.1 用于表示数值大于或等于0.4的列。该列表示在指定簇中表达特定基因的细胞百分比（Cluster %），将阈值设为0.4意味着仅显示在该簇至少40%的细胞中表达的基因。
    5. 过滤 pct.2 用于表示小于或等于 0.2 的数值的列。该列表示未在指定簇中的细胞表达指定基因的百分比（非簇细胞百分比），将阈值设为 0.2 意味着仅显示在最多 20% 的非指定簇细胞中表达的基因。
    6. 筛选 p_val_adj 小于或等于 0.01 的数值对应的列。该列指校正后的 P 值或错误发现率（FDR），用于表示所指出的差异表达基因（DEG）的统计显著性，将阈值设为 0.01 意味着仅保留 FDR < 0.01 的结果被显示。
      注：补充表1 （JoVE_DEGs_cellMarkers.xlsx）包含在方案步骤 3.10 中使用的排序后的差异表达基因的完整输出结果。 补充表 2 显示本分析中每个簇的前 5 个基因，加粗的基因为后续可视化所使用。
  11. 为实现对细胞簇的无偏倚细胞类型注释，使用基于网页的富集分析工具 EnrichR。
    使用链接： https://maayanlab.cloud/Enrichr/
  12. 将每个簇的差异表达基因（DEG）列表分别复制到独立的 EnrichR 窗口中，然后点击分析。EnrichR 工具会将基因列表与数百个经过人工整理的数据库进行比对，并对每个类别中富集的术语进行排序。
  13. 为了进行细胞类型注释，点击上方的“细胞类型”标签页，并关注左侧三个经人工整理的细胞标志物数据库中排名前五的富集结果（图3):
    CellMarker 2024http://bio-bigdata.hrbmu.edu.cn/CellMarker/)
    人类图谱（Tabula Sapiens）https://tabula-sapiens-portal.ds.czbiohub.org/)
    PanglaoDB 增强版（https://panglaodb.se/)
  14. 根据这些数据库中差异表达基因（DEGs）的富集情况，确认8个聚类的可能身份。注意有两个聚类（2、6）富集为成纤维细胞，因此将这两个聚类合并为单一的细胞类型注释。将细胞类型身份作为标签分配给名为 cell_types 的新元数据变量。
    Idents(数据集) <- "seurat_clusters"
    数据集[["细胞类型"]] <- Idents(数据集)
    Idents(数据集) <- "细胞类型"
    数据集 <- 重命名Idents(数据集,
    "0" = "巨噬细胞",
    "1" = "中性粒细胞",
    "2" = "成纤维细胞",
    "3" = "上皮细胞",
    "4" = "内皮细胞",
    "5" = "T细胞",
    "6" = "成纤维细胞",
    "7" = "平滑肌细胞"
    )
    dataset 的水平
    数据集[["细胞类型"]] <- Idents（数据集）
  15. 将重命名的细胞簇作为注释可视化呈现在UMAP图上（图4).
    DimPlot(数据集, group.by = "细胞类型", raster = FALSE, label = FALSE
  16. 根据表1中加粗显示的前导聚类标记基因，在一系列UMAP图上可视化其定位情况（图5).
    DefaultAssay(dataset) <- "RNA"
    FeaturePlot(数据集, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4
  17. 在点图上按原始聚类编号分组展示各聚类前导标志基因的差异表达基因（DEGs）补充图6).
    DotPlot(数据集, group.by = "seurat_clusters"，特征 = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "蓝色"，中 = "灰色"，高 = "红色")
  18. 在点图上可视化按注释细胞类型分组的各簇前几位标志基因的差异表达基因（DEGs）图6).
    DotPlot(数据集, group.by = "细胞类型"，特征 = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "蓝色"，mid = "灰色"，高 = "红色")
  19. 为了进行时间序列分析，首先简化数据集以去除空间成分。对于时间过程分析，将伤口时间/空间注释按伤后总天数（DPW）进行分组，并创建一个名为 "伤后天数".
    Idents(数据集) <- "时间与空间"
    数据集[["天后（days post-wounding, DPW）"]] <- Idents（数据集）
    Idents(数据集) <- "DPW"
    new.cluster.ids <- c("D1", "D1", "D1", "D1", "D3", "D3", "D3", "D3", "D7", "D7", "D7", "D7", "D14", "D14", "D14", "D14", "UW")
    names(new.cluster.ids) <- 查看数据集的因子水平
    数据集 <- 重命名 dataset 中的细胞簇标识为 new.cluster.ids
    数据集[["DPW"]] <- Idents（数据集）
    Idents(数据集) <- "天（DPW）"
  20. 在 UMAP 图上可视化新的伤口时间序列分组（补充图7).
    DimPlot(数据集, group.by = "天后（days post-wounding, DPW）"，raster = FALSE，label = FALSE）
  21. 生成表格，显示每个损伤后时间点（DPW）中每种细胞类型的细胞数量。
    ​table(dataset$DPW, dataset$cell_types)
    1. 可选步骤：获取伤口时间序列各组的差异表达基因（DEG）列表，将其赋值给变量，并将输出结果保存为分隔符分隔的文本文件。
      Idents(数据集) <- "DPW"
      细胞_伤后天数标记基因 <- FindAllMarkers(数据集, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
      write.csv(Cell_DPW_markers, file = file.path(getwd(), "dataset_DPW_markers.txt"))
  22. 将各细胞类别的细胞数量转换为类别比例，以更好地理解愈合过程中细胞类型组成随时间的相对变化。可视化每种细胞类型中DPW的比例补充图8):
    第1部分 <- table(dataset$DPW, dataset$cell_types)
    pt1 <- as.data.frame(pt1)
    pt1$Var1 <- as.character(pt1$Var1)
    ggplot(pt1, aes(x = Var2, y = Freq, fill = Var1)) +
    theme_bw(base_size = 15) +
    geom_col(position = "填充"，width = 0.5) +
    xlab("样本"）+ 旋转坐标轴() +
    ylab("比例") +
    theme(legend.title = element_blank())
  23. 可视化各时间点（DPW）中细胞类型的比例（ 图7).
    第2部分 <- table(dataset$cell_types, dataset$DPW)
    pt2 <- as.data.frame(pt2)
    pt2$Var1 <- as.character(pt2$Var1)
    ggplot(pt2, aes(x = Var2, y = Freq, fill = Var1)) +
    theme_bw(base_size = 15) +
    geom_col(position = "填充"，width = 0.5) +
    xlab("样本") +
    ylab("比例") +
    theme(legend.title = element_blank())
  24. 将 Seurat 对象数据集保存为 RDS 文件至工作目录。
    saveRDS(dataset, "dataset_post_Method3.rds")

4. 使用 Seurat 分析细胞亚型

注意：单细胞分析的强大之处在于能够发现和分析上述主要细胞类型中存在的稀有亚型。本示例聚焦于成纤维细胞，这些细胞最初被聚类为两个 Seurat 聚类，随后被合并为一个类别。本部分实验方案专门针对成纤维细胞，排除所有其他细胞类型，以探究其在伤口愈合过程中的身份特征和时间动态特性。作为可选步骤，如有需要，可加载已保存的 RDS 文件作为 Seurat 对象。

dataset <- readRDS("dataset_post_Method3.rds")

  1. 根据成纤维细胞身份对原始数据集进行子集筛选。
    Idents(dataset) <- "cell_types"
    dataset_fibroblast <- subset(dataset, idents = "Fibroblast")
  2. 对该较小的数据集执行主成分分析（PCA），并可视化数据集在PCA各维度上的变异程度（补充图9）。
    dataset_fibroblast <- RunPCA(dataset_fibroblast, verbose = TRUE)
    ElbowPlot(dataset_fibroblast, reduction = "pca", ndims = 50)
    注意： 大部分主要变异集中在前9个维度内。
  3. 使用PCA维度范围1–9及设定分辨率为0.1，对数据集进行细胞聚类。
    dataset_fibroblast <- FindNeighbors(dataset_fibroblast, verbose = TRUE, dims = 1:9)
    dataset_fibroblast <- FindClusters(dataset_fibroblast, verbose = TRUE, resolution = 0.1)
    注意： 使用这些参数设置，算法可区分出3个独特的成纤维细胞聚类。这些聚类将自动分配至一个名为 seurat_clusters 的元数据变量中。
  4. 使用前9个PCA维度执行UMAP降维及邻域分析。添加随机种子数123，以确保数据投影结果的可重复性。
    dataset_fibroblast <- RunUMAP(dataset_fibroblast, verbose = TRUE, dims = 1:9, seed.use = 123)
  5. 在UMAP图上可视化细胞的聚类情况（图8）。
    DimPlot(dataset_fibroblast, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
  6. 在UMAP图上可视化细胞的伤口时间进程注释（补充图10）。
    DimPlot(dataset_fibroblast, group.by = "DPW", raster = FALSE, label = FALSE)
  7. 获取三个成纤维细胞亚型的差异表达基因（DEG）列表，并将其保存至工作目录中的文本文件。
    Idents(dataset_fibroblast) <- "seurat_clusters"
    fibroblast_markers <- FindAllMarkers(dataset_fibroblast, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
    write.csv(fibroblast_markers, file = file.path(getwd(), "fibroblast_cluster_markers.txt"))
  8. 在电子表格中参照上述步骤（步骤3.9–3.10），将差异表达基因筛选为各成纤维细胞亚型的前导标志基因。
  9. 通过从差异表达基因文本文件中复制每个聚类的前5个成纤维细胞亚型标志基因，定义一个自定义基因列表变量。
    FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
  10. 在仅含成纤维细胞的数据集中，通过在点图（dotplot）的 features 参数中调用该变量，可视化该基因列表的表达情况（图9）。
    DotPlot(dataset_fibroblast, group.by="seurat_clusters", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
    DotPlot(dataset_fibroblast, group.by="DPW", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
  11. 在原始单细胞数据集中，通过在点图的 features 参数中调用该变量，可视化该基因列表的表达情况（补充图11）。
    DotPlot(dataset, group.by="cell_types", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
  12. 可视化各DPW类别中成纤维细胞亚型的比例分布（补充图12）。
    pt3 <- table(dataset_fibroblast$seurat_clusters, dataset_fibroblast$DPW)
    pt3 <- as.data.frame(pt3)
    pt3$Var1 <- as.character(pt3$Var1)
    ggplot(pt3, aes(x = Var2, y = Freq, fill = Var1)) +
    theme_bw(base_size = 15) +
    geom_col(position = "fill", width = 0.5) +
    xlab("样本") +
    ylab("比例") +
    theme(legend.title = element_blank())
  13. 可视化各成纤维细胞亚型类别中DPW成纤维细胞的比例分布（补充图13）。
    pt4 <- table(dataset_fibroblast$DPW, dataset_fibroblast$seurat_clusters)
    pt4 <- as.data.frame(pt4)
    pt4$Var1 <- as.character(pt4$Var1)
    ggplot(pt4, aes(x = Var2, y = Freq, fill = Var1)) +
    theme_bw(base_size = 15) +
    geom_col(position = "fill", width = 0.5) +
    xlab("样本") +
    ylab("比例") +
    theme(legend.title = element_blank())
  14. 将处理后的Seurat数据集对象保存为RDS文件至工作目录。
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

5. 通过模块评分进行后续分析的示例

注意： 分析单细胞数据集的一种有用方法称为模块评分。在此工作流程中，可以根据先前的知识定义一个基因列表，然后计算模块评分，从而识别每个细胞内该基因列表的潜在富集情况。这些评分可以在细胞注释间进行平均，以揭示潜在的富集模式。

此处使用先前发表的研究2中的基因列表，该研究通过愈合连续过程中采集的批量RNA测序样本鉴定了伤口愈合各阶段特异性基因。这些基因列表已保存为制表符分隔的文本文件（补充文件 2：JoVE_PhaseSpecificGenes.txt），可下载至工作目录，用于生成能够识别三个主要愈合阶段的基因列表。

  1. 通过读取文本文件，将基因列表加载到变量中。
    阶段特异性基因 <- read_delim("JoVE_阶段特异性基因.txt"，分隔符="伤口愈合过程受到不同细胞类型在时空上复杂相互作用的调控。通过对复杂微环境中单个细胞进行分析，单细胞转录组学方法能够研究参与伤口愈合过程的细胞异质性、细胞通讯网络以及细胞间相互作用。然而，许多单细胞分析工具需在编程环境中运行，由于缺乏生物信息学专业知识，阻碍了伤口愈合研究者对这些工具的广泛应用。因此，本文提供了一套逐步操作流程，展示如何使用名为 RStudio 的图形化编程环境，对小鼠皮肤切除伤口愈合的时间序列数据集进行基础的单细胞分析。这一可视化、引导式操作指南将帮助无生物信息学背景的研究人员下载已发表的伤口愈合数据集，执行关键的质量控制步骤，使用 Seurat 完成标准的单细胞分析流程（包括数据集可视化和细胞类型注释），开展细胞亚群分析、模块评分分析，利用 CellChat 进行细胞间相互作用分析，并使用 Seurat 对多个数据集进行整合分析。本方案为每一步操作提供了文字说明，并展示了每一行代码对应的图形化结果，以安全引导用户完成整个分析流程。本单细胞分析流程的可视化导引旨在使更多伤口愈合领域的科研人员能够在其实验室直接使用生物信息学工具，从而促进对其自身单细胞数据集的深入分析，以及对已发表单细胞数据集的广泛再分析。", col_names = T)
  2. 将各列拆分为独立的基因列表变量，并将基因名称更改为首字母大写的鼠源基因名称。
    PS_炎症性 <- 时期特异性基因[1]
    PS_炎症性_ms <- lapply(PS_Inflammatory, str_to_sentence)
    PS_增殖期 <- PhaseSpecificGenes[2]
    PS_增殖期_小鼠 <- lapply(PS_Proliferative, str_to_sentence)
    PS_分辨率 <- PhaseSpecificGenes[3]
    PS_分辨率_ms <- lapply(PS_Resolution, str_to_sentence)
    注意： （可选）如有需要，将保存的 RDS 文件作为 Seurat 对象加载。
    数据集 <- readRDS("dataset_post_Method3.rds")
  3. 使用基因列表作为模块，根据愈合的三个阶段对数据集中的每个细胞进行评分。
    数据集 <- AddModuleScore(
    object = dataset,
    features = PS_Inflammatory_ms
    ctrl = 100,
    name = '炎症性'
    )
    数据集 <- AddModuleScore(
    object = dataset,
    features = PS_Proliferative_ms,
    ctrl = 100,
    name = '增殖期'
    数据集 <- AddModuleScore(
    object = dataset,
    features = PS_Resolution_ms,
    ctrl = 100,
    名称 = '分辨率'
    )
  4. 按细胞类别（包括DPW和主要细胞类型）可视化聚合的模块评分图10).
    DotPlot(数据集, group.by="天"，特征 = c("炎症性1","增殖期1", "分辨率1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "蓝色"，中 = "灰色"，高 = "红色")
    DotPlot(数据集, group.by="细胞类型"，特征 = c("炎症性1","增殖期1", "分辨率1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "蓝色"，中 = "灰色"，高 = "红色")

6. 通过 CellChat 进行后续分析的示例

注意：分析单细胞数据集的另一种有用且被广泛引用的方法是推断细胞间相互作用。在此工作流程中，使用 CellChat 软件包，该软件包通过分析不同细胞群之间配体-受体相互作用的差异来推断细胞间通讯22。最近，CellChat 的开发人员发布了一份详细的逐步操作协议24，供用户通用参考，这在用户执行以下工作流程并将其应用于自身数据集时是非常优秀的资源。例如，以下工作流程比较了伤口愈合后1天与14天（DPW）时所有主要细胞间的相互作用。由于所有步骤已在 CellChat 的官方出版物24以及相关教程中详细描述，此处不再详述每一步操作，教程链接如下：

使用
CellChat 进行细胞间通信的推断与分析：https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html

使用 CellChat 进行多个数据集的比较分析：
https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html

可选步骤：如有需要，将已保存的 RDS 文件作为 Seurat 对象载入：

dataset <- readRDS("dataset_post_Method3.rds")

  1. 根据 DPW 注释将原始数据集划分为两个数据集。
    Idents(dataset) <- "DPW"
    dataset_D1 <- subset(dataset, ident = "D1")
    dataset_D14 <- subset(dataset, ident = "D14")
  2. 定义用于执行 CellChat 的注释——本例中使用主要细胞类型。
    Idents(dataset_D1) <- "cell_types"
    Idents(dataset_D14) <- "cell_types"
  3. 创建 CellChat 对象并遵循典型的 CellChat 工作流程。有关每一步的详细参考，请参见上方链接的教程。
    cellchat_D1 <- createCellChat(dataset_D1, group.by = "ident", assay = "RNA")
    cellchat_D14 <- createCellChat(dataset_D14, group.by = "ident", assay = "RNA")
    CellChatDB <- CellChatDB.mouse
    CellChatDB.use <- CellChatDB
    cellchat_D1@DB <- CellChatDB.use
    cellchat_D14@DB <- CellChatDB.use
    cellchat_D1 <- subsetData(cellchat_D1)
    cellchat_D14 <- subsetData(cellchat_D14)
    future::plan("multisession", workers = 4)
    cellchat_D1 <- identifyOverExpressedGenes(cellchat_D1, do.fast = F)
    cellchat_D14 <- identifyOverExpressedGenes(cellchat_D14, do.fast = F)
    cellchat_D1 <- identifyOverExpressedInteractions(cellchat_D1)
    cellchat_D14 <- identifyOverExpressedInteractions(cellchat_D14)
    cellchat_D1 <- computeCommunProb(cellchat_D1, type = "triMean", population.size = TRUE)
    cellchat_D14 <- computeCommunProb(cellchat_D14, type = "triMean", population.size = TRUE)
    cellchat_D1 <- filterCommunication(cellchat_D1, min.cells = 10)
    cellchat_D14 <- filterCommunication(cellchat_D14, min.cells = 10)
    cellchat_D1 <- computeCommunProbPathway(cellchat_D1)
    cellchat_D14 <- computeCommunProbPathway(cellchat_D14)
    cellchat_D1 <- aggregateNet(cellchat_D1)
    cellchat_D14 <- aggregateNet(cellchat_D14)
    cellchat_D1 <- netAnalysis_computeCentrality(cellchat_D1, slot.name = "netP")
    cellchat_D14 <- netAnalysis_computeCentrality(cellchat_D14, slot.name = "netP")
  4. 可视化在每个伤口愈合时间点所有主要细胞类型的传入与传出相互作用强度（补充图14）。
    netAnalysis_signalingRole_scatter(cellchat_D1)
    netAnalysis_signalingRole_scatter(cellchat_D14)
    成纤维细胞在伤后第1天（D1）与第14天（D14）之间的相互作用显著增强。
  5. 显示所有显著推断出的细胞间通讯通路列表。
    cellchat_D1@netP$pathways
    cellchat_D14@netP$pathways
    胶原蛋白通路是伤后第1天和第14天均显著的重要通路之一。
  6. 聚焦于胶原信号通路及其与成纤维细胞的相互作用。
    pathways.show <- c("COLLAGEN")
  7. 使用环形图可视化不同细胞类型间胶原信号通路的相互作用（补充图15）。
    par(mfrow=c(1,2))
    netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "circle")
    netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "circle")
    par(mfrow=c(1,1))
  8. 使用弦图可视化不同细胞类型间胶原信号通路的相互作用（补充图16）。
    par(mfrow=c(1,2))
    strwidth <- function(x) {0.5}
    netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6)
    netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6)
    par(mfrow=c(1,1))
  9. 可视化以成纤维细胞为信号来源细胞时，COLLAGEN 信号通路的相互作用（补充图17）。
    注意： CellChat 对象中的细胞类型按其在原始 Seurat 对象中分配的顺序以编号列出：1 = 巨噬细胞，2 = 中性粒细胞，3 = 成纤维细胞，4 = 上皮细胞，5 = 内皮细胞，6 = T 细胞，7 = 平滑肌细胞。
    par(mfrow=c(1,2))
    strwidth <- function(x) {0.5}
    netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6, sources.use = 3)
    netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6, sources.use = 3)
    ​par(mfrow=c(1,1))
  10. 可视化以成纤维细胞为信号来源细胞时，COLLAGEN 信号通路中各配体-受体对的贡献。
    1. 使用气泡图（补充图18）：
      gg1 <- netVisual_bubble(cellchat_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, remove.isolate = FALSE)
      gg2 <- netVisual_bubble(cellchat_D14, sources.use = 3, targets.use = NULL, signaling = pathways.show, remove.isolate = FALSE)
      ​gg1 + gg2
    2. 使用弦图（补充图19）：
      par(mfrow=c(1,2))
      strwidth <- function(x) {0.4}
      netVisual_chord_gene(cellchat_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, show.legend= F)
      netVisual_chord_gene(cellchat_D14, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, legend.pos.x = 60)
      ​par(mfrow=c(1,1))
  11. 聚焦于 COLLAGEN 信号通路中的 Col1a1-Cd44 配体-受体相互作用。
    ​LR.show <- "COL1A1_CD44"
  12. 使用弦图可视化不同细胞类型间 Col1a1-Cd44 配体-受体相互作用（补充图20）。
    strwidth <- function(x) {0.5}
    netVisual_individual(cellchat_D1, signaling = pathways.show, pairLR.use = LR.show, layout = "chord")
    netVisual_individual(cellchat_D14, signaling = pathways.show, pairLR.use = LR.show, layout = "chord")
  13. 通过生成一个合并的 CellChat 对象，执行差异性 CellChat 分析。
    object.list_D14_v_D1 <- list(D1 = cellchat_D1, D14 = cellchat_D14)
    cellchat_D14_v_D1 <- mergeCellChat(object.list_D14_v_D1, add.names = names(object.list_D14_v_D1))
  14. 可视化伤口愈合时间点中细胞间相互作用的总数及相对强度（补充图21）。
    gg1 <- compareInteractions(cellchat_D14_v_D1, show.legend = F)
    gg2 <- compareInteractions(cellchat_D14_v_D1, show.legend = F, measure = "weight")
    gg1 + gg2
  15. 使用环形图可视化伤口从第1天过渡到第14天过程中，各细胞类型间细胞间相互作用强度的差异（补充图22）。
    netVisual_diffInteraction(cellchat_D14_v_D1, weight.scale = T, measure = "weight")
  16. 使用热图可视化伤口从第1天过渡到第14天过程中，各细胞类型间细胞间相互作用强度的差异（补充图23）。
    netVisual_heatmap(cellchat_D14_v_D1, measure = "weight")
  17. 使用排序图可视化以成纤维细胞为信号来源细胞时，各通路在第14天与第1天对细胞间相互作用的相对贡献（补充图24）。
    rankNet(cellchat_D14_v_D1, mode = "comparison", measure = "weight", sources.use = 3, targets.use = NULL, stacked = T, do.stat = TRUE)
  18. 使用气泡图可视化以成纤维细胞为信号来源细胞时，胶原信号通路中各配体-受体对在第14天相对于第1天的相对贡献（补充图25）。
    gg1 <- netVisual_bubble(cellchat_D14_v_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, comparison = c(1, 2), max.dataset = 2, title.name = "第14天信号增强", angle.x = 45, remove.isolate = F)
    gg2 <- netVisual_bubble(cellchat_D14_v_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, comparison = c(1, 2), max.dataset = 1, title.name = "第14天信号减弱", angle.x = 45, remove.isolate = F)
    gg1 + gg2
  19. 与 Seurat 对象类似，CellChat 对象可保存并打开为 RDS 文件。
    saveRDS(cellchat_D1, "cellchat_D1.rds")
    saveRDS(cellchat_D14, "cellchat_D14.rds")
    saveRDS(cellchat_D14_v_D1, "cellchat_D14_v_D1.rds")
  20. 可选步骤：也可从 RDS 文件中打开 CellChat 对象。
    cellchat_D1 <- readRDS("cellchat_D1.rds")
    cellchat_D14 <- readRDS("cellchat_D14.rds")
    cellchat_D14_v_D1 <- readRDS("cellchat_D14_v_D1.rds")

7. 通过整合多个单细胞数据集进行整合分析的示例

注意：单细胞数据集通常被分为多个文件，因为它们是分批次或分组测序的。本工作流程展示了如何整合伤口愈合数据集五个批次中的其中两个20。当前的数据集整合方法在以下 Seurat 示例文档中有详细描述：

单细胞RNA测序整合简介： 
https://satijalab.org/seurat/articles/integration_introduction

Seurat v5 中的整合分析：

https://satijalab.org/seurat/articles/seurat5_integration

注意：单细胞数据集整合方法有多种，每种方法都有其优势和局限性。详情请参见整合方法的全面基准评估25。用户在依赖任何一种整合方法之前，务必阅读所有相关文档。

  1. 对 Hu 等人数据集的另一批次重复方法 2 中的所有步骤20。在以下方案中，使用第 3 批次（batch #3）。补充的 R 脚本文件已提供，可用于处理第 3 批次（补充文件 3：JoVE_Rscript_b3.R）。请记住为数据集创建并使用一个新的变量——在下方代码中，使用“dataset_b3”表示第 3 批次数据集。
    1. 可选步骤：如有需要，从工作目录中保存的 RDS 文件打开两个数据集作为 Seurat 对象：
      dataset <- readRDS("dataset_post_Method2.rds")
      dataset_b3 <- readRDS("dataset_b3_post_Method2.rds")
  2. 为每个数据集分配一个名为“batch”的新变量，以便在后续分析中标记数据集来源。
    dataset@meta.data$batch <- "b1"
    dataset_b3@meta.data$batch <- "b3"
  3. 使用 Seurat 合并两个数据集，添加基于批次的细胞 ID 注释，然后对合并后的数据集执行方法 3 中所述的标准 Seurat 分析流程。
    dataset_merged <- merge(x = dataset, y = c(dataset_b3), add.cell.ids = c("b1", "b3"), merge.data = TRUE)
    DefaultAssay(dataset_merged) <- "RNA"
    dataset_merged <- NormalizeData(dataset_merged)
    dataset_merged <- FindVariableFeatures(dataset_merged)
    dataset_merged <- ScaleData(dataset_merged)
    dataset_merged <- RunPCA(dataset_merged)
    ElbowPlot(dataset_merged, reduction = "pca", ndims = 50)
  4. 在数据整合之前，对合并后的数据集进行聚类和 UMAP 分析。
    dataset_merged <- FindNeighbors(dataset_merged, dims = 1:10, reduction = "pca")
    dataset_merged <- FindClusters(dataset_merged, resolution = .1, cluster.name = "unintegrated_clusters")
    dataset_merged <- RunUMAP(dataset_merged, dims = 1:10, seed.use = 123, reduction = "pca", reduction.name = "umap.unintegrated")
  5. 根据聚类和批次编号可视化 UMAP 图（补充图 26）。
    DimPlot(dataset_merged, reduction = "umap.unintegrated", group.by = c("seurat_clusters", "batch"))
  6. 显示各聚类中按批次编号分布的细胞数量。
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    注意： 从 UMAP 图和表格来看，这两个数据集之间似乎不存在显著的批次效应。批次效应的证据表现为两个数据集之间聚类分布出现意外差异，这可能意味着数据集之间存在潜在的技术差异，掩盖了真实的生物学相似性。
  7. 使用 RPCA 方法进行 Seurat 数据整合。有关此方法及其他数据整合方法的更多信息，请参阅上方链接的 Seurat 教程文档。
    dataset_merged <- IntegrateLayers(
    object = dataset_merged, method = RPCAIntegration,
    orig.reduction = "pca", new.reduction = "integrated.rpca",
    verbose = TRUE
    )
  8. 在数据整合后，对合并后的数据集进行聚类和 UMAP 分析。
    FindNeighbors(dataset_merged, dims = 1:10, reduction = "integrated.rpca")
    dataset_merged <- FindClusters(dataset_merged, resolution = .1, cluster.name = "rpca_clusters")
    dataset_merged <- RunUMAP(dataset_merged, dims = 1:10, seed.use = 123, reduction = "integrated.rpca", reduction.name = "umap.rpca")
  9. 在整合后根据聚类和批次编号可视化 UMAP 图（补充图 27）。
    DimPlot(dataset_merged, reduction = "umap.rpca", group.by = c("seurat_clusters","batch"))
  10. 在整合后显示各聚类中按批次编号分布的细胞数量。
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    从整合后数据的 UMAP 图和表格可以看出，两个批次在不同聚类之间现在具有良好的重叠。有趣的是，在使用相同聚类参数的情况下，数据整合还额外识别出一个聚类。
  11. 在数据集整合完成且进行下游分析之前，必须将合并数据集的图层进行合并。
    dataset_merged <- JoinLayers(dataset_merged)
  12. 将Seurat 对象数据集保存为 RDS 文件至工作目录。
    saveRDS(dataset_merged, "merged_dataset_post_Method7.rds")

结果

从方法 #2 开始，本方案逐步介绍如何加载单细胞伤口愈合数据集并对其进行质量控制。在创建 Seurat 对象（步骤 2.6.2）后，通过一系列操作将数据集中的两个检测内容（RNA 和蛋白质；步骤 2.6.3–2.6.7）进行合并，并根据时空条形码对蛋白质检测数据进行解复用（步骤 2.6.8–2.6.9）。解复用函数为数据集中的每个细胞分配多个元数据标签，包括“barcodes_maxID”，用于标识每个细胞最可能对应的时空条形码（步骤 2.6.10）。在步骤 2.6.11 中，使用小提琴图函数可视化基于多重条形码的细胞中检测到的基因分布情况。该步骤的代表性结果（补充图 1）显示，每种条形码对应的检测基因数量分布较为均匀，这对于数据集的完整性以及伤口愈合时间点的下游分析至关重要。在为蛋白质条形码分配适当的标签后（步骤 2.6.12），本方案进一步展示了如何对数据集的 RNA 检测部分进行质量控制，首先计算每个细胞中线粒体基因所占百分比（步骤 2.9）。在步骤 2.10 中，使用特征散点图函数可视化所有细胞中检测基因数、RNA 分子数以及线粒体基因百分比的分布情况。该步骤的代表性结果（补充图 2）显示，部分细胞具有较高的线粒体基因含量，这与较低的 RNA 计数相关，提示这些细胞可能为死亡或濒死细胞。在去除 RNA 计数较低且线粒体基因含量较高的细胞后（步骤 2.11），于步骤 2.12 对筛选后的子集数据再次执行特征散点图函数，该步骤的代表性结果（补充图 3）显示，每个细胞中检测到的基因数及线粒体 RNA 百分比的分布现已趋于正常，为后续稳健的分析奠定了基础。接下来，本方案描述了如何使用 scDblFinder 函数识别数据集中可能的双细胞（doublet）事件，并为每个细胞分配一个名为“scDblFinder.score”的新元数据（步骤 2.13–2.14）。在步骤 2.15 中，使用小提琴图函数可视化数据集中双细胞评分的分布情况，该步骤的代表性结果（补充图 4）显示，部分细胞具有相对较高的双细胞评分，且 0.25 似乎为一个自然的分界值，高于该值的细胞群体很可能为双细胞。因此，后续步骤采用该阈值对数据集进行筛选，保留评分低于该阈值的细胞（步骤 2.16），从而完成对该单细胞数据集的质量控制流程。

从方法 #3 开始，本实验方案逐步介绍如何使用 Seurat 软件包及工作流程分析经质量控制的单细胞伤口愈合数据集。首先对 RNA 数据进行标准化和缩放处理，随后进行主成分分析（PCA）（步骤 3.1）。在步骤 3.2 中，使用肘部图函数可视化前 50 个 PCA 维度中数据集变异程度的变化，该步骤的代表性结果（补充图 5）显示，大部分主要变异集中在前 13 个维度内，这由图中曲线的拐点所确定。接着，方案展示了如何基于前 13 个 PCA 维度以及较低的聚类分辨率参数（0.1）来寻找细胞邻居并进行细胞聚类（步骤 3.3）以及 UMAP 降维分析（步骤 3.4），选择这些参数旨在识别伤口中最具普适性的主要细胞类型。在步骤 3.5 中，通过降维图函数将细胞聚类结果可视化于 UMAP 图中，该步骤的代表性结果（图 1）显示，数据集中所有细胞围绕 8 个主要的颜色编码 Seurat 聚类组分布，且在运行 Windows（左侧）与 MacOS（右侧）的计算机上获得的 UMAP 图略有差异。在步骤 3.6 中，再次使用降维图函数可视化细胞的时间/空间注释信息，该步骤的代表性结果（图 2）显示，数据集中所有细胞根据其来源的时间/空间信息分散分布，未见明显按时间/空间注释形成的聚类。接下来，方案描述了如何获取差异表达基因列表并将其保存为文本文件（步骤 3.8），在电子表格中打开数据表，并执行多种筛选步骤以获得每个细胞聚类的最高排名聚类标志基因（步骤 3.9–3.10.6）。这些步骤的代表性结果（补充表 1）为包含全部排序后差异表达基因的最终电子表格文件，而另一代表性结果（补充表 2）则为简化表格，展示每个 Seurat 聚类中前 5 个上调且高表达的基因。随后，方案描述了如何使用名为 EnrichR 的基于网页的功能富集分析工具，根据各聚类的顶级标志基因推断潜在的细胞类型（步骤 3.11–3.12），这些步骤的代表性结果（图 3）为 EnrichR 输出结果的截图，展示了八个细胞聚类中各自富集程度最高的细胞类型。随后，方案根据各 Seurat 聚类中最显著富集的细胞类型注释，为所有细胞分配一个名为“cell_types”的新元数据标签（步骤 3.14）。在步骤 3.15 中，使用降维图函数将重命名后的细胞聚类以细胞类型注释的形式可视化于 UMAP 图中，该步骤的代表性结果（图 4）显示，数据集中所有细胞围绕主要颜色编码的细胞类型聚类分布。在步骤 3.16 中，使用特征图函数将顶级聚类标志基因（来自补充表 2）在一系列 UMAP 图中进行定位可视化，其代表性结果（图 5）为一组 UMAP 图网格，显示各主要细胞类型聚类区域内其对应顶级标志基因的高表达情况。在步骤 3.17 和 3.18 中，使用点图函数可视化顶级聚类标志基因在细胞中的相对表达水平，首先按原始 Seurat 聚类编号分组（步骤 3.17），再按注释的细胞类型标签分组（步骤 3.18）。这些步骤的代表性结果证实，顶级细胞标志基因仅在其对应的 Seurat 聚类中高表达（补充图 6），也仅在其对应的各大细胞类型中高表达（图 6）。方案的下一步将原始基于空间-时间的蛋白标签简化为严格的时序注释，即根据细胞来源的伤口后天数（DPW）进行标识。在步骤 3.20 中，使用降维图函数将细胞以 DPW 注释形式可视化于 UMAP 图中，该步骤的代表性结果（补充图 7）展示了伤口时间序列注释在整个单细胞伤口愈合数据集中的分布情况。如预期所示，第 1 天（D1）的注释主要集中在中性粒细胞和巨噬细胞聚类中，而较晚的伤口愈合时间点则在其他细胞类型中更为丰富。方案后续步骤使用堆叠柱状图，首先可视化不同细胞类型中 DPW 的比例分布（步骤 3.22），然后可视化不同时间点中细胞类型的比例分布（步骤 3.23）。这些步骤的代表性结果为比例图，分别展示每种主要细胞类型类别中 DPW 细胞的相对数量（补充图 8）以及每个 DPW 类别中主要细胞类型的相对数量（图 7）。这些结果验证了已知的皮肤伤口愈合细胞级联过程：在炎症阶段早期，免疫细胞（中性粒细胞和巨噬细胞）占主导地位；在增殖阶段，其他细胞类型（上皮细胞和内皮细胞）开始出现，而成纤维细胞在伤口修复后期尤为显著。

从方法 #4 开始，本方案概述了使用 Seurat 对单细胞数据集中的某一主要细胞类型进行分析的步骤，以鉴定伤口愈合过程中潜在的细胞亚型。本方案聚焦于成纤维细胞，这些细胞最初被聚类为两个 Seurat 聚类，随后合并为单一类别，并描述了如何创建一个新的 Seurat 对象，该对象仅包含原始数据集中来自成纤维细胞的数据（步骤 4.1）。随后在该成纤维细胞特异性数据集上执行 Seurat 工作流程（步骤 4.2–4.4），其中步骤 4.2 生成一个拐点图（补充图 9），显示成纤维细胞数据集中的大部分主要变异发生在前 9 个 PCA 维度内。在步骤 4.5 中，使用降维图函数在 UMAP 图上可视化细胞的聚类情况，该步骤的代表性结果（图 8）显示数据集中的成纤维细胞围绕三种颜色编码的细胞亚型聚集。根据其 DPW 注释对成纤维细胞数据集进行可视化（步骤 4.6），得到一张 UMAP 图（补充图 10），显示数据集中的成纤维细胞按照其 DPW 注释广泛分布。接下来，方案描述了如何获取差异表达基因列表并将其保存为文本文件（步骤 4.7），在 Excel 中打开数据表并执行多种筛选步骤，以获得每个细胞聚类中排名靠前的聚类标志基因（步骤 4.8），并创建一个名为“FB_type_marker”的新变量，列出排名靠前的成纤维细胞标志基因（步骤 4.9）。在步骤 4.10 中，使用点图函数通过在 features 参数中调用“FB_type_marker”变量，在仅含成纤维细胞的数据集中可视化该基因列表，该步骤的代表性结果（图 9）为点图，证实成纤维细胞亚型标志基因仅在其对应的聚类类别中高表达（上图），但在 DPW 类别中分布较为均匀（下图）。在步骤 4.11 中，调用相同的 features 变量以在整体伤口愈合数据集中可视化成纤维细胞标志基因，其代表性结果（补充图 11）是一张点图，证实成纤维细胞亚型标志基因主要在原始成纤维细胞中高表达。最后，方案后续步骤使用堆叠柱状图，首先可视化三种成纤维细胞亚型中 DPW 比例的分布情况（步骤 4.12），然后可视化不同时间点成纤维细胞亚型的比例分布（步骤 4.13）。这些步骤的代表性结果为比例图，分别显示每种成纤维细胞亚型类别中 DPW 细胞的相对数量（补充图 12）以及每种 DPW 类别中成纤维细胞亚型的相对数量（补充图 13）。这些结果表明，在伤口愈合的时间进程中，成纤维细胞亚型的比例发生了显著变化：第一种成纤维细胞亚型（聚类 0）在早期伤口（D1 和 D3）中占主导地位，第二种亚型（聚类 1）在伤口消退期（D14）占主导地位，而第三种亚型（聚类 2）在伤口愈合的增殖期（D7）表达水平最高。

从方法步骤 #5 开始，本方案逐步介绍如何使用 Seurat 中的模块评分功能分析单细胞伤口愈合数据集。方案首先描述了如何使用制表符分隔的文本文件将基因集上传至 R 中的变量（步骤 5.1–5.2），随后将模块评分功能应用于与伤口愈合三个主要阶段相关的三个基因集（步骤 5.3）。在步骤 5.4 中，使用点图函数可视化两个不同元数据类别中的综合模块评分，该步骤的代表性结果（图 10）为点图，分别显示在受伤后天数类别（DPW，右侧）和主要细胞类型类别（左侧）中，各细胞内主要愈合阶段模块的平均表达水平。这些结果表明，以伪批量方式将基于批量测序的基因表达谱应用于单细胞表达数据集，是一种强大的比较生物信息学方法，可充分利用伤口愈合领域已发表的数据集。

从方法 #6 开始，本方案根据特定科学问题——即比较早期与晚期伤口来源的细胞——逐步介绍如何使用 CellChat 软件包及工作流程分析由 Seurat 生成的单细胞伤口愈合数据集。首先，将整体 Seurat 数据集按损伤后两个时间点进行子集划分：一个处于炎症期（伤后第 1 天（D1）），另一个处于伤口修复期（伤后第 14 天（D14））（步骤 6.1）。随后创建两个 CellChat 对象，并执行 CellChat 方案中所有典型功能，以计算在本方案方法 #3 中鉴定出的各细胞类型之间的所有潜在相互作用（步骤 6.2–6.3）。在步骤 6.4 中，执行信号散点图功能，以可视化各主要细胞类型在两个伤口愈合时间点的输入和输出相互作用强度。该步骤的代表性结果（补充图 14）为散点图，分别显示 D1（左）和 D14（右）时间点各主要细胞类型的输入相互作用强度（y 轴）与输出相互作用强度（x 轴）。这些结果表明，中性粒细胞和巨噬细胞等免疫细胞在炎症期具有最强的细胞间相互作用强度，而纤维细胞则在伤口修复期主导了细胞间相互作用，这与数十年来的伤口愈合研究结果一致。后续步骤聚焦于一条显著富集的信号通路——胶原蛋白通路（步骤 6.5–6.6）。在步骤 6.7 中，执行环形图功能，以可视化两个时间点间各细胞类型在胶原蛋白信号通路中的相互作用。该步骤的代表性结果（补充图 15）为环形图，显示 D1（左）和 D14（右）时间点所有细胞类型间推断出的胶原蛋白通路信号相互作用。在步骤 6.8 中，使用弦图功能对相同相互作用进行可视化，其代表性结果（补充图 16）为弦图，展示各时间点所有细胞类型间推断出的胶原蛋白通路信号相互作用。如预期所示，这些结果表明纤维细胞是胶原蛋白信号通路的主要信号源细胞，尽管在 D1 时信息流主要局限于免疫细胞，而在 D14 时更为广泛。为进一步聚焦纤维细胞作为信号源细胞在细胞间相互作用中的角色，步骤 6.9 重复执行弦图功能，并添加源细胞参数，其代表性结果（补充图 17）为弦图，展示以纤维细胞为源细胞时，各时间点所有细胞类型间推断出的胶原蛋白通路信号相互作用。在步骤 6.10 中，执行两个功能以可视化在以纤维细胞为源细胞的情况下，胶原蛋白信号通路中每一对配体-受体的贡献：一个使用气泡图（步骤 6.10.1），另一个使用弦图（步骤 6.10.2）。代表性结果通过气泡图（补充图 18）和弦图（补充图 19）展示 D1（左）和 D14（右）时间点以纤维细胞为源细胞时，胶原蛋白通路中每一对配体-受体的推断贡献。这些结果表明，在 D1 时，来自纤维细胞的胶原蛋白通路主要局限于中性粒细胞和巨噬细胞，且以 Cd44 和 Sdc4 受体为主；而在 D14 时，其他细胞通过包括整合素在内的多种受体作为信号接收细胞。为进一步聚焦 Col1a1-Cd44 配体-受体相互作用（该相互作用在纤维细胞间表现出较强强度），在步骤 6.11 中设置相应参数，并在步骤 6.12 中将其用于弦图功能，以可视化所有细胞类型间该特定配体-受体相互作用，其代表性结果（补充图 20）为弦图，展示 D1（左）和 D14（右）时间点所有细胞类型间推断出的 Col1a1-Cd44 配体-受体相互作用。这些结果表明，在 D1 时该相互作用的源细胞仅限于纤维细胞，而在 D14 时，巨噬细胞和平滑肌细胞也作为源细胞参与其中。接下来，本方案描述如何进行差异性 CellChat 分析：首先合并 D1 和 D14 的 CellChat 对象（步骤 6.13）。在步骤 6.14 中，执行“比较相互作用”功能，以可视化两个伤口愈合时间点间细胞间相互作用的总数及相对强度，其代表性结果（补充图 21）为柱状图，分别显示 D1 和 D14 伤口中细胞推断出的相互作用总数（左）和强度（右），其中 D14 的相互作用数量更高，而 D1 的相互作用相对强度更高。在步骤 6.15 和 6.16 中，使用两个功能分别可视化伤口从第 1 天过渡到第 14 天过程中各细胞类型间细胞间相互作用强度的差异，其相应代表性结果分别为环形图（步骤 6.15，补充图 22）和热图（步骤 6.16，补充图 23），其中 D14 相较于 D1 增强的相互作用以红色表示，减弱的以蓝色表示。如预期所示，中性粒细胞和巨噬细胞介导的相互作用在 D1 更强，而纤维细胞介导的相互作用在 D14 更强。在步骤 6.17 中，使用排序功能生成一个排序图，以比较 D14 与 D1 时间点以纤维细胞为源细胞时，各通路对细胞间相互作用的相对贡献。其代表性结果（补充图 24）显示排序图，D1 以红色位于上方，D14 以蓝色位于下方，多个通路仅在 D1 或 D14 中特异性表达，其余则呈现激活梯度。最后，在步骤 6.18 中，使用两个气泡图功能展示 D14 与 D1 时间点以纤维细胞为源细胞时，胶原蛋白信号通路中各配体-受体对的相对贡献，其相应代表性结果（补充图 25）显示在 x 轴上多个细胞间相互作用中，D14 相较于 D1 增强（左）和减弱（右）的信号配对。如预期所示，在 D14 伤口中，纤维细胞向多个接收细胞发出的配体-受体对相互作用显著增加，而在 D1 伤口中，信号交流在炎症期主要局限于中性粒细胞和巨噬细胞。

从方法步骤 #7 开始，本方案逐步介绍使用 Seurat 整合两个单细胞伤口愈合数据集的操作流程。首先描述合并已发表的两个单细胞数据集批次，并对合并后的数据集应用标准的 Seurat 分析流程（步骤 7.1–7.4）。在步骤 7.5 中，使用降维图函数可视化尚未整合的合并伤口愈合数据集按聚类和批次编号分布的 UMAP 图。该步骤的代表性结果（补充图 26）为 UMAP 图，分别显示 Seurat 聚类分布（左图）和批次编号分布（右图），表明在数据整合前，这两个数据集之间未见明显的批次效应。随后，方案采用 RPCA 方法进行数据整合，并继续对整合后的数据集执行后续的 Seurat 分析流程（步骤 7.7–7.8）。在步骤 7.9 中，再次使用降维图函数，按聚类和批次编号可视化整合后伤口愈合数据集的 UMAP 图。该步骤的代表性结果（补充图 27）为 UMAP 图，分别显示 Seurat 聚类分布（左图）和批次编号分布（右图），结果显示两个批次在不同聚类间的重叠程度进一步增加。结果还显示，在数据整合后出现了一个新的额外聚类，这可能提示在控制了数据批次的技术效应后，识别潜在重要细胞亚型的能力有所提升。

Windows 与 MacOS 上的 UMAP 聚类比较，seurat_clusters，数据可视化图表。
图 1：UMAP 图显示数据集中所有细胞围绕 8 个主要颜色编码的聚类组进行聚类。 结果分别来自运行 Windows（左）和 MacOS（右）的计算机。该图对应于步骤 3.5。请单击此处查看此图的放大版本。

UMAP结果示意图；降维数据可视化；不同颜色表示不同聚类。
图2：UMAP图显示数据集中所有细胞根据其时间/空间来源分布展开，未见依据时间/空间注释形成明显聚类。 本图对应步骤3.6。 请点击此处查看该图的放大版本。

细胞标志物分析图表，聚类显示差异表达基因在数据库中的分布情况。
图 3：EnrichR 分析结果的截屏（已裁剪），显示各细胞聚类中富集程度最高的细胞类型。 本图对应步骤 3.13。请点击此处查看该图的放大版本。

细胞类型UMAP图；巨噬细胞、中性粒细胞、成纤维细胞。降维分析。
图4：UMAP图显示数据集中所有细胞围绕主要颜色编码的细胞类型聚类。 本图对应步骤3.15。请点击此处查看该图的放大版本。

UMAP基因表达图谱；Arg1、Retnlg、Fgf7、Dsp、Tie1、Cd3g、Cdh4、Rgs5的分析。
图5：UMAP图网格展示主要细胞类型聚类中前导细胞标志物基因的高表达。 本图对应步骤3.16。请点击此处查看该图的放大版本。

细胞类型标志基因表达点图，显示各特征的平均表达水平和表达百分比。
图6：点图验证了各主要细胞类型的标志性基因仅在其对应细胞类型中高表达。 本图对应步骤3.18。请点击此处查看此图的放大版本。

细胞类型比例条形图；内皮细胞、上皮细胞、成纤维细胞、巨噬细胞、中性粒细胞、平滑肌细胞、T细胞。
图7：各DPW类别中主要细胞类型相对数量的比例图。 该图对应于步骤3.23。请点击此处查看此图的放大版本。

UMAP 聚类图；seurat_clusters；数据可视化；降维；基因表达分析
图 8：UMAP 图显示数据集中成纤维细胞围绕三种不同颜色编码的细胞亚型聚类。 该图对应步骤 4.5。请点击此处查看该图的放大版本。

基因表达点图；可视化百分比和平均表达水平；示意图；数据分析。
图9：点图证实成纤维细胞亚型标志物仅在其对应的簇类别中高表达，但在DPW类别中分布较为均匀。 本图对应步骤4.10。请点击此处查看此图的放大版本。

细胞表达水平气泡图；特征：炎症性、增殖性、消退性；数据：细胞身份、表达百分比、平均表达量。
图10：点图显示不同愈合阶段模块在每DPW及主要细胞类型中的细胞平均表达水平。 本图对应步骤5.4。请点击此处查看该图的放大版本。

补充图1：结果显示每个条形码检测到的基因分布较为均匀，这对于数据集的完整性和伤口愈合时间点的下游分析至关重要。该图对应于步骤2.6.11。请点击此处下载该图。

补充图 2：散点图显示存在大量线粒体含量较高的细胞，这些细胞与较低的 RNA 计数相关——这些是死亡或即将死亡的细胞。 该图对应于步骤 2.10。请点击此处下载该图。

补充图3：散点图显示，检测到的基因分布以及每个细胞的线粒体RNA百分比现更趋近于正态分布，为可靠的下游分析奠定了基础。 该图对应于步骤2.12。请点击此处下载该图。

补充图4：小提琴图显示存在一些双细胞（doublet）评分相对较高的细胞，且0.25看起来是一个自然的截断值，高于该值的细胞群体很可能是双细胞（doublet）。 该图对应于步骤2.15。请点击此处下载该图。

补充图5：肘部图显示大部分主要变异发生在前13个维度内。 该图对应于步骤3.2。请点击此处下载该图。

补充图6：点图验证了各细胞标志物基因仅在其对应的Seurat聚类中高表达。 该图对应步骤3.17。请点击此处下载该图。

补充图7：UMAP图显示伤口愈合数据集中伤口时间序列注释的定位。 该图对应于步骤3.20。请点击此处下载该图。

补充图8：显示每种主要细胞类型类别中DPW细胞相对数量的比例图。 该图对应于步骤3.22。请点击此处下载该图。

补充图9：肘部图显示，成纤维细胞数据集中的大部分主要变异发生在前9个维度内。 该图对应于步骤4.2。请点击此处下载该图。

补充图10：UMAP图显示根据DPW注释对数据集中的成纤维细胞进行分布展示。 该图对应步骤4.6。请点击此处下载该图。

补充图11：点图证实成纤维细胞亚型标志物主要在原始成纤维细胞聚类中高表达。 该图对应步骤4.11。请点击此处下载该图。

补充图12：显示每个DPW类别中成纤维细胞亚型相对数量的比例图。 该图对应于步骤4.12。请点击此处下载该图。

补充图13：显示在每个成纤维细胞亚型类别中跨DPW的成纤维细胞相对数量的比例图。 该图对应于步骤4.13。请点击此处下载该图。

补充图14：散点图显示了在第1天（D1，左侧）和第14天（D14，右侧）主要细胞类型传入（y轴）和传出（x轴）相互作用的强度。 该图对应于步骤6.4。请点击此处下载该图。

补充图15：圆形图显示了每个DPW类别中所有细胞类型之间推断的胶原蛋白通路信号相互作用。 该图对应于步骤6.7。请点击此处下载该图。

补充图16：弦图显示了每个DPW类别中所有细胞类型之间推断的胶原蛋白通路信号相互作用。 该图对应于步骤6.8。请点击此处下载该图。

补充图17：弦图显示了推断的胶原蛋白通路信号相互作用，其中成纤维细胞作为每种DPW类别的来源细胞。 该图对应于步骤6.9。请点击此处下载该图。

补充图18：气泡图显示了在每个DPW类别中以成纤维细胞为信号来源细胞时，胶原蛋白通路信号传导中每对配体-受体的推断贡献。 本图对应步骤6.10.1。请点击此处下载该图。

补充图19：弦图显示了在每个DPW类别中以成纤维细胞为信号来源细胞时，胶原蛋白通路信号传导中每对配体-受体的推断贡献。 该图对应步骤6.10.2。请点击此处下载该图。

补充图 20：显示每个 DPW 类别中所有细胞类型间推断的 Col1a1-Cd44 配体-受体相互作用的弦图。 该图对应于步骤 6.12。请点击此处下载该图。

补充图21：显示第1天和第14天伤口中推断出的相互作用数量（左）和强度（右）的柱状图。该图对应于步骤6.14。请点击此处下载该图。

补充图22：圆形图显示伤口从伤后第1天（蓝色）过渡到第14天（红色）过程中，各细胞类型之间细胞-细胞相互作用强度的差异。 该图对应步骤6.15。请点击此处下载该图。

补充图23：热图显示伤口从第1天（蓝色）到第14天（红色）愈合过程中各细胞类型之间细胞-细胞相互作用强度的差异。 该图对应于步骤6.16。请点击此处下载该图。

补充图24：显示成纤维细胞与其他细胞类型在损伤后第1天与第14天之间细胞间相互作用中各通路相对贡献度的排序图。 该图对应步骤6.17。请点击此处下载该图。

补充图25：气泡图显示在第1天与第14天DPW时，成纤维细胞作为信号来源细胞时，胶原信号通路中各配体-受体对的相对贡献。该图对应于步骤6.18。请点击此处下载该图。

补充图26：数据整合前Seurat聚类（左）和批次编号（右）的UMAP分布图。 该图对应步骤7.5。请点击此处下载该图。

补充图27：数据整合后Seurat聚类（左）和批次编号（右）的UMAP图分布。 该图对应步骤7.9。请点击此处下载该图。

补充文件 1：JoVE_Rscript.R： 主 R 代码脚本文件，包含本方案所有部分所述的全部步骤和说明。请点击此处下载该文件。

补充文件 2：JoVE_PhaseSpecificGenes.txt。 一个以制表符分隔的文本文件，其中包含在实验方案第 5.1 步中加载的基因列表。 请点击此处下载该文件。

补充文件 3：JoVE_Rscript_b3.R。 补充的 R 代码脚本文件，包含分析数据集第 3 批次所需的所有步骤和说明，用于本方案的步骤 7.1。 请点击此处下载该文件。

补充表 1：JoVE_DEGs_cellMarkers.xlsx。 该 Excel 文件包含本方案第 3.10 步所用的排序后的差异表达基因完整结果。请单击此处下载该表格。

补充表 2：各 Seurat 聚类中前 5 个上调且高表达的基因。 请点击此处下载该表格。

讨论

在本实验方案中，使用 RStudio 运行指定的代码行，以利用 Seurat 对复杂的单细胞数据集进行基础分析。介绍了多种与伤口愈合研究相关的方法，包括 R 编程环境的安装、下载已发表的单细胞伤口愈合数据集、执行关键的质量控制步骤以及标准的单细胞分析工作流程（包括可视化、主要细胞类型注释、细胞亚型分析和基于 Seurat 的整合分析），并使用 CellChat 进行细胞间相互作用分析。

本文介绍的方法是使用 R 及其流行的开源科学软件包 Seurat21 和 CellChat22 进行单细胞分析的典型工作流程的简化示例。实际上，该工作流程仅展示了利用复杂的单细胞伤口愈合数据集所能完成的分析类型之一。对该方法可能进行的修改几乎是无限的，唯一的限制在于用户特定的科学问题。例如，用户可根据希望从该数据集中提出的研究问题，调整某些关键参数，如细胞类型和时间点。作者还希望用户能够自如地将此工作流程适配于自身感兴趣的单细胞数据集；然而，在使用该工作流程分析其他数据集时必须谨慎，因为每个实验都可能将技术性问题和样本制备问题带入数据本身。因此，在解读任何已发表并被重新分析的单细胞数据集结果之前，用户必须仔细阅读并充分理解所有实验细节。需要牢记的是，生物信息学工具是探索生物学过程和提出假设的有力手段，但对结果所作的任何关键性生物学解释都必须通过后续实验加以验证。

在整个实验流程中，需注意可根据具体需求对工作流程的特定环节进行重大修改，以实现其他目标。然而，所有可能的流程修改组合的详细说明已超出本文的范围。例如，用于细胞聚类的分辨率以及UMAP分析所采用的维度本质上具有主观性，本文介绍的工具既支持大规模分析（如本研究中对广义主要细胞类型的分析），也支持更精细的特定分析，后者可能涉及将细胞进一步划分为更大数据集中更稀有的亚群。关于单细胞分析方法的这一方面，以及单细胞分析流程中其他可调整参数的详细信息，作者建议读者参考Seurat的相关出版物21,26及官方网站（https://satijalab.org/seurat/），该工具的开发者在网站上提供了深入的解释、示例分析（vignettes）和操作教程。

本文介绍了单细胞转录组学文献中引用和使用最广泛的几种工具，分别是用于单细胞分析的 Seurat21 和用于细胞间相互作用分析的 CellChat22。然而，也存在其他以略微不同方式实现类似功能的工具。在单细胞数据集分析方面，有 Scran27、Scater28 以及基于 Python 的 ScanPy29，这些工具采用多种方法进行数据集整合25。本实验方案展示了细胞类型的 manual 注释方法，该方法依赖用户对聚类细胞标志物富集结果的判断；但目前已有多种可实现细胞类型自动分类的工具，例如 SingleR30 和 scGate31 等。在细胞间通讯分析方面，本方案演示了 CellChat 的使用，但还存在其他用于推断细胞间通讯的工具，包括 CellPhoneDB32、Cytotalk33，以及集成在 LIANA（LIgand-receptor ANalysis framework）共识框架中的其他配体-受体数据库34。所有生物信息学工具均具有独特性，各自带有特定的特性与可调节参数。因此，用户必须仔细阅读每种工具的相关文档，以充分理解其细微差异，再对分析结果进行解释。最后，无论使用何种生物信息学工具，都必须牢记这些工具在持续发展，不同版本的软件包可能产生不同的分析结果。

在 R 中，语法至关重要，标点符号、引号、括号的位置错误，甚至字母大小写错误，都会导致报错。因此，用户在输入代码时必须注意细节，尤其是在复制代码行并将其修改以适应新的科学问题和数据集时更应格外谨慎。针对可能遇到的具体错误，作者建议将错误信息直接复制粘贴到用户常用的网络搜索引擎中，并浏览来自生物信息学论坛（如 GitHub 和 Stack Overflow）的搜索结果，因为大多数常见错误很可能已被经验丰富的高级用户解答。在某些论坛中，其他用户会通过“点赞”方式标记出他们认为最有效的解决方案。用户必须注意，切勿直接将从互联网上找到的代码行复制粘贴到自己的计算机中（特别是当解决方案要求更改 R 编程环境之外的系统设置时），因为这些程序可能存在恶意风险。一种新兴且有效的排错方法是使用强大的生成式大型语言人工智能模型，例如 OpenAI 的 ChatGPT、微软的 Copilot 或谷歌的 Gemini。这些模型已被证明在软件工程领域，尤其是问题排查方面具有显著帮助。使用时，用户可在向聊天机器人简要说明代码意图后，将其代码整行复制粘贴给模型。需要注意的是，这些模型并非万无一失，用户可能需要尝试多个提示词（prompt），才能获得适用于解决当前问题的正确答案。

披露

作者声明无利益冲突。

致谢

M.S. Wietecha 实验室获得了美国国立卫生研究院/国家普通医学科学研究所（NIH/NIGMS）R35-GM154921 项目基金、伤口愈合学会研究基金以及伊利诺伊大学芝加哥分校牙科学院口腔生物学系的资助。

材料

本文使用的材料清单
姓名公司目录编号评论
笔记本电脑或台式计算机N/AN/A运行 Windows 或 MacOS 
RN/A版本 4.4.1可从 https://cran.rstudio.com/ 免费下载
RstudioPosit 软件，PBC版本 2024.09.0可从 https://posit.co/download/rstudio-desktop/ 免费下载
Office Excel微软任意版本用于表格数据的分析
互联网浏览器N/AN/A用于访问网站
R 软件包知识库版本
开发工具CRAN2.4.5
readxlCRAN1.4.3
openxlsxCRAN4.2.7.1
tidyverseCRAN2.0.0
scCustomizeCRAN2.1.2
BiocManagerBioconductor1.30.25
NMFBioconductor0.28
ComplexHeatmapBioconductor2.20.0
BiocNeighborsBioconductor1.22.0
SingleCellExperimentBioconductor1.26.0
circlizeBioconductor0.4.16
edgeRBioconductor4.2.1
scDblFinderBioconductor1.18.0
SeuratCRAN5.1.0
CellChatGithub2.1.2

参考文献

  1. Eming, S. A., Martin, P., Tomic-Canic, M. Wound repair and regeneration: mechanisms, signaling, and translation. Sci Transl Med. 6 (265), 265sr6(2014).
  2. Wietecha, M. S., et al. Phase-specific signatures of wound fibroblasts and matrix patterns define cancer-associated fibroblast subtypes. Matrix Biol. 119, 19-56 (2023).
  3. Rodrigues, M., Kosaric, N., Bonham, C. A., Gurtner, G. C. Wound healing: a cellular perspective. Physiol Rev. 99 (1), 665-706 (2019).
  4. Chen, L., Mirza, R., Kwon, Y., DiPietro, L. A., Koh, T. J. The murine excisional wound model: contraction revisited. Wound Repair Regen. 23 (6), 874-877 (2015).
  5. Rhea, L., Dunnwald, M. Murine excisional wound healing model and histological morphometric wound analysis. J Vis Exp. (162), e61616(2020).
  6. Fischer, K. S., et al. Protocol for the splinted, human-like excisional wound model in mice. Bio-Protocol. 13 (3), e4606(2023).
  7. Iglesias-Bartolome, R., et al. Transcriptional signature primes human oral mucosa for rapid wound healing. Sci Transl Med. 10 (451), aap8798(2018).
  8. Chen, L., Arbieva, Z. H., Guo, S., Marucha, P. T., Mustoe, T. A., DiPietro, L. A. Positional differences in the wound transcriptome of skin and oral mucosa. BMC Genomics. 11, 471(2010).
  9. Leonardo, T. R., et al. Transcriptional changes in human palate and skin healing. Wound Repair. 31 (2), 156-170 (2023).
  10. Rognoni, E., et al. Fibroblast state switching orchestrates dermal maturation and wound healing. Mol Syst Biol. 14 (8), e8174(2018).
  11. Bergmeier, V., et al. Identification of a myofibroblast-specific expression signature in skin wounds. Matrix Biol. 65, 59-74 (2018).
  12. Plikus, M. V., et al. Regeneration of fat cells from myofibroblasts during wound healing. Science. 355 (6326), 748-752 (2017).
  13. Shook, B. A., et al. Myofibroblast proliferation and heterogeneity are supported by macrophages during skin repair. Science. 362 (6417), aar2971(2018).
  14. Rinkevich, Y., et al. Skin fibrosis. Identification and isolation of a dermal lineage with intrinsic fibrogenic potential. Science. 348 (6232), aaa2151(2015).
  15. Guerrero-Juarez, C. F., et al. Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds. Nat Commun. 10 (1), 650(2019).
  16. Gay, D., et al. Phagocytosis of Wnt inhibitor SFRP4 by late wound macrophages drives chronic Wnt activity for fibrotic skin healing. Sci Adv. 6 (12), eaay3704(2020).
  17. Haensel, D., et al. Defining epidermal basal cell states during skin homeostasis and wound healing using single-cell transcriptomics. Cell Rep. 30 (11), 3932-3947.e6 (2020).
  18. Phan, Q. M., Sinha, S., Biernaskie, J., Driskell, R. R. Single-cell transcriptomic analysis of small and large wounds reveals the distinct spatial organization of regenerative fibroblasts. Exp Dermatol. 30 (1), 92-101 (2021).
  19. Foster, D. S., et al. Integrated spatial multiomics reveals fibroblast fate during tissue repair. Proc Natl Acad Sci U S A. 118 (41), e2110025118(2021).
  20. Hu, K. H., et al. Transcriptional space-time mapping identifies concerted immune and stromal cell patterns and gene programs in wound healing and cancer. Cell Stem Cell. 30 (6), 885-903.e10 (2023).
  21. Hao, Y., et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol. 42 (2), 293-304 (2024).
  22. Jin, S., et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 12 (1), 1088(2021).
  23. Germain, P. -L., Lun, A., Garcia Meixide, C., Macnair, W., Robinson, M. D. Doublet identification in single-cell sequencing data using scDblFinder. F1000Research. 10, 979(2021).
  24. Jin, S., Plikus, M. V., Nie, Q. CellChat for systematic analysis of cell-cell communication from single-cell transcriptomics. Nat Protoc. 20 (1), 180-219 (2024).
  25. Luecken, M. D., et al. Benchmarking atlas-level data integration in single-cell genomics. Nat Methods. 19 (1), 41-50 (2022).
  26. Hao, Y., et al. Integrated analysis of multimodal single-cell data. Cell. 184 (13), 3573-3587.e29 (2021).
  27. Lun, A. T. L., McCarthy, D. J., Marioni, J. C. A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Research. 5, 2122(2016).
  28. McCarthy, D. J., Campbell, K. R., Lun, A. T. L., Wills, Q. F. Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics (Oxford, England). 33 (8), 1179-1186 (2017).
  29. Wolf, F. A., Angerer, P., Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19 (1), 15(2018).
  30. Aran, D., et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 20 (2), 163-172 (2019).
  31. Andreatta, M., Berenstein, A. J., Carmona, S. J. scGate: marker-based purification of cell types from heterogeneous single-cell RNA-seq datasets. Bioinformatics. 38 (9), 2642-2644 (2022).
  32. Efremova, M., Vento-Tormo, M., Teichmann, S. A., Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. 15 (4), 1484-1506 (2020).
  33. Hu, Y., Peng, T., Gao, L., Tan, K. CytoTalk: de novo construction of signal transduction networks using single-cell transcriptomic data. Sci Adv. 7 (16), eabf1356(2021).
  34. Dimitrov, D., et al. Comparison of methods and resources for cell-cell communication inference from single-cell RNA-Seq data. Nat Commun. 13 (1), 3224(2022).

重印与许可

标签

Seurat 分析CellChat 工作流程小鼠皮肤数据集细胞类型注释UMAP 可视化差异基因表达细胞间通讯生物信息学工作流程
使用 R、Seurat 和 CellChat 分析小鼠皮肤伤口愈合的单细胞转录组数据集
联系我们
向图书馆推荐

隐私政策

使用条款

政策

版权所有 © 2026 MyJoVE Corporation。保留所有权利。