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

Summary

在这里，我们提出了一个分步的可视化工作流程，用于使用 R 分析小鼠皮肤伤口愈合的单细胞时程转录组学数据集。该协议包括一个标准管道，用于使用 Seurat 进行数据集下载、质量控制、可视化和细胞类型注释，以及使用 CellChat 进行细胞间相互作用分析。

Abstract

伤口愈合的过程受到不同细胞类型之间跨空间和时间的复杂相互作用的调节。通过对复杂环境中的单个细胞进行分析，单细胞转录组学方法能够研究伤口愈合过程中涉及的细胞异质性、细胞通讯网络和细胞间相互作用。然而，许多单细胞分析工具都是在计算机编码环境中运行的，由于明显缺乏生物信息学专业知识，伤口愈合科学家更广泛地使用它们受到阻碍。因此，提出了一个分步工作流程，展示了如何使用名为RStudio的图形编码环境对颞小鼠切除皮肤伤口愈合数据集进行基本的单细胞分析。这种可视化和引导式协议将使没有生物信息学背景的科学家能够下载以前发布的伤口愈合数据集，执行关键的质量控制步骤，运行标准的单细胞分析工作流程，包括使用 Seurat 的数据集可视化和细胞类型注释，运行细胞亚型分析，运行模块评分分析，使用 CellChat 运行细胞-细胞相互作用分析，并使用 Seurat 对多个数据集进行综合分析。为协议中的每个步骤提供叙述性解释，并呈现每行代码的图形结果，以安全地指导用户完成工作流程。这种对单细胞分析管道的可视化介绍的目标是使更多的伤口愈合科学家能够直接在自己的实验室中使用生物信息学工具，以便对他们自己的单细胞数据集进行更深入的分析，以及对以前发表的单细胞数据集进行更广泛的重新分析。

Introduction

伤口愈合是哺乳动物生物学中最复杂的过程之一，它涉及三个愈合阶段：炎症、增殖和决议 ^1,2。这些愈合阶段对数十种细胞类型及其数百种分子产物在伤口修复空间和时间上的协调作用进行了广泛分类³。数十年的组织学和分子研究基于整个愈合时间过程的伤口组织取样，已经阐明了^{组织修复的}总体细胞模式3，特别是在切除皮肤伤口愈合的可重复小鼠模型中^4,5,6。直到最近二十年，才有可能更充分地理解伤口愈合的复杂性，首先是对大块组织 7^、8、9 和细胞^{10、11、12、13、14} 尺度的伤口进行高通量转录组学分析.最近，几项研究在单细胞水平上对皮肤伤口进行了转录分析，识别了新的伤口细胞亚型，并展示了它们在愈合过程中如何相互作用 15,16,17,18,19,20。胡等人使用一种创新的空间单细胞 RNA 测序方法，在距伤口中心几个径向距离处的整个愈合时间过程中对皮肤伤口进行分析，这揭示了跨空间和时间的新细胞间和分子“运动”²⁰。这些研究正在以前所未有的细节揭示伤口愈合的复杂性，并且它们开始描绘出一幅巨大的细胞和分子异质性的图景。

生物信息学分析方法的最新重大进展使得对伤口愈合研究领域生成的复杂多组学数据集进行生物学理解成为可能。像 Seurat 这样的单细胞分析包提供了用于对数据集进行稳健分析和集成的工具，包括对伤口等复杂组织中的细胞类型进行分类²¹。对于单细胞数据的下游解释，CellChat等工具用于识别假定的细胞-细胞相互作用程序，这些程序可以解释细胞如何协调以修复伤口²²。虽然这些工具有据可查且被引用，但它们必须在计算机编码环境中运行，例如 R，这是一种统计和图形编程语言，最常用于基因组学和转录组学的生物信息学领域。虽然伤口愈合领域的生物学家和临床医生越来越多地使用单细胞方法来研究组织修复，但很少有人接受过直接在自己的实验室中使用 Seurat 和 CellChat 等工具所需的生物信息学培训。使用这些生物信息学工具的这种障碍不仅阻止科学家在没有生物信息学家帮助的情况下更深入地分析自己的数据集，而且还阻止科学家可靠地重新分析其他小组已经发表的大量单细胞数据。

因此，这里提出了一个分步工作流程，使没有生物信息学背景的科学家能够分析先前发表和公开可用的单细胞伤口愈合数据集²⁰。该协议使用称为 RStudio 的流行且免费的图形 R 编码环境，它演示了如何导航该环境以运行规定的代码行，从而能够使用 Seurat 和 CellChat 对复杂的单细胞数据集进行基本分析。在该协议中，提出了与伤口愈合研究相关的七种主要方法，包括：1）编码环境安装，2）下载数据集和关键质量控制步骤，3）单细胞分析工作流程，包括可视化和细胞类型注释，4）细胞亚型分析，5）模块评分分析，6）细胞-细胞相互作用分析，以及7）多个数据集的综合分析。在每种方法中，都提供了实际代码供用户与协议并排运行，并显示每行代码的实际图形结果以指导用户完成工作流程。对RStudio和基本单细胞分析工作流程的指导和可视化介绍的主要目标是使更多的伤口愈合科学家能够直接使用这些强大的工具，以便在研究领域取得更快的进展。

Protocol

注意： 在以下详细介绍七种生物信息学方法的工作流程中，协议的所有步骤都附有各自的代码块，这些代码块应按照列出的顺序直接在用户自己的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 菜单栏中选择 Session 并单击 Set Working Directory > Choose Directory 并 选择所需的文件夹来设置工作目录。
   1. 如果使用 Windows 计算机，请使用以下命令设置工作目录。将以下代码行中的 [Directory] 更改为实际的目录结构。请注意，R 中的目录分隔符是字符“/”
      setwd("C:/[Directory]")
   2. 如果使用 MacOS 计算机，以下命令还将设置工作目录。将以下代码行中的 [Directory] 更改为实际的目录结构。请注意，R 中的目录分隔符是字符“/”
      setwd("~/[Directory]")
   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 包的过程中，各种窗口出现和消失是正常的。如果出现一个窗口，要求编译包，请单击 YES。如果出现一个窗口，要求在安装包之前重新启动 R，请单击 NO。
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. 加载单细胞伤口愈合数据集并执行质量控制步骤

注意： 对于这个生物信息学工作流程，对先前发表的时空单细胞皮肤伤口愈合实验进行了重新分析²⁰.数据集文件存储在精选的 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 函数执行解复用。以下修拉小插图详细描述了此方法：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 毫米）分配（取自原始手稿），并将它们分配给名为 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. 作为重要的质量控制步骤，计算每个细胞中线粒体基因的百分比并将其分配为元数据变量。以下修拉小插图中详细描述了此方法：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 计数相关;这些是死亡或垂死的细胞。使用公平截止值从数据集中删除这些低质量的单元格。在这种情况下，使用先前发表的研究²⁰ 中描述的原始数据集的值，其中线粒体基因超过 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. 作为另一个重要的质量控制步骤，检测数据集中可能的双峰。这些细胞是在液滴测序过程中连接的，因此将导致不在单细胞水平的基因表达。为了解决这个问题，已经开发了几种工具。需要注意的是，每个工具都对单细胞数据做出了某些假设，因此用户在使用任何一种工具之前阅读所有相关文档非常重要。在此工作流中，实现名为 scDblFinder²³ 的方法。请注意，此方法使用施加固定预期双峰率的算法。使用以下命令运行 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
   修拉命令列表： https://satijalab.org/seurat/articles/essential_commands
   DefaultAssay（dataset） <- “RNA”
   dataset <- NormalizeData(object = dataset)
dataset <- FindVariableFeatures(object = dataset)
dataset <- ScaleData(object = dataset)
dataset <- RunPCA(object = dataset)
2. 可视化相对于前 50 个 PCA 维度的数据集变化量（补充图 5）。
   ElbowPlot(dataset, reduction = "pca", ndims = 50)
   大部分主要变化发生在前 13 个维度内。
3. 使用设置的 PCA 维度范围为 1-13 和设置的分辨率为 0.1 对数据集执行单元格聚类。
   dataset <- FindNeighbors(dataset, verbose = TRUE, dims = 1:13)
dataset <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
   注意： 该工作流程侧重于细胞类型的大规模差异。因此，它对 PCA 维度范围使用了一个相当保守的参数，其中前 13 个维度显示为代表数据集中的绝大多数变异。为了将细胞区分为更小、更稀有的亚型，用户可以使用更多的维度进行下游分析，因为这些稀有细胞亚型可能占较低水平的数据集变异。分辨率参数范围为 0 到 1，它决定了对数据集施加的分类分离的大小。该参数的设置取决于用户的研究问题。要将细胞聚类为许多小而罕见的亚型，请使用更高分辨率的下游分析。由于该工作流程旨在探索主要细胞类型之间更广泛的差异，因此它使用相当小的分辨率值 0.1，预计该值将细胞聚类为更少、更大的组。使用步骤 3.3 中提到的设置将区分 8 个唯一的单元簇。这些集群会自动分配给名为“seurat_clusters”的元数据变量。
4. 使用前 13 个 PCA 维度执行 UMAP 降维和查找邻域分析。添加种子编号 123 以确保生成的数据投影的可重复性。
   dataset <- RunUMAP(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(dataset, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
   注意： UMAP 算法的随机性可能会生成略有不同的图，如图所示，显示了在运行 Windows 和 MacOS 的计算机上使用与上述相同的代码生成的替代 UMAP 图;请注意簇形状的细微差异。因此，用户必须在生成所有数据和图时保存它们并添加时间戳，并且仔细进行所有对簇的下游分析并牢记生物学理解，如下所述细胞类型注释。
6. 由于包括指代伤口愈合期间细胞来自何时何地的原始实验标签，因此在 UMAP 图上可视化细胞的伤口时间/空间注释（图 2）。
   DimPlot(dataset, group.by = "time_space", raster = FALSE, label = FALSE)
7. 生成与伤口时间/空间注释的细胞簇关联表。
   table(dataset$time_space, dataset$seurat_clusters)
8. 确定数据集中主要细胞类型的身份。为此，请计算所有簇之间的差异表达基因 （DEG）。获取集群的 DEG 列表，将它们分配给变量，并将输出保存为工作目录中的分隔文本文件。
   注意： 此步骤占用大量 CPU，可能需要很长时间，具体取决于用户的硬件。
   Idents(dataset) <- "seurat_clusters"
Cell_markers <- FindAllMarkers(dataset, 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. 通过复制文本文件的内容并使用文本导入向导将逗号分隔符和基因名称列的标识指定为文本，下载并打开电子表格（例如 Excel）中包含的 dataset_cluster_markers.txt 文件。指示基因名称为“文本”很重要，否则 Excel 会自动将某些基因名称转换为日期，例如 Sept7 更改为 September 7。
10. 在电子表格中，根据以下推荐参数筛选结果：
    1. 将 avg_log2FC 列从大到小排列，以便根据递减的对数₂ 倍变化 （log2FC） 排列所有行。
    2. 将 聚类 列从小到大排列，以根据增加的 Seurat 聚类数排列所有行。
    3. 筛选 avg_log2FC 列以查找大于或等于 2.5 的数字，以仅显示指示聚类中差异表达最多的基因 （DEG） 与其他聚类。
    4. 筛选 pct.1 列以查找大于或等于 0.4 的数字。该列是指指定簇中表达指定基因的细胞百分比（聚类%），将阈值设置为 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 基于 Web 的富集分析工具。
    使用链接： https://maayanlab.cloud/Enrichr/
12. 将每个集群的 DEG 列表复制到单独的 EnrichR 窗口中，然后单击分析。EnrichR 工具通过数百个精选数据库运行基因列表，并对每个类别中的每个富集术语进行排名。
13. 出于细胞类型注释的目的，单击上面的“细胞类型”选项卡，并关注左侧三个细胞标记精选数据库中的前 5 个富集（图 3）：
    细胞标记 2024 （http://bio-bigdata.hrbmu.edu.cn/CellMarker/）
    智人表 （https://tabula-sapiens-portal.ds.czbiohub.org/）
    PanglaoDB 增强型 （https://panglaodb.se/）
14. 根据这些数据库中 DEG 的丰富程度，确认 8 个集群的可能身份。请注意，有两个簇 （2， 6） 富集为成纤维细胞;因此，将这些聚类组合成单细胞类型注释。将单元格类型标识作为标签分配给名为 cell_types 的新元数据变量。
    Idents(dataset) <- "seurat_clusters"
dataset[["cell_types"]] <- Idents(dataset)
Idents(dataset) <- "cell_types"
dataset <- RenameIdents(dataset,
"0" = "Macrophage",
"1" = "Neutrophil",
"2" = "Fibroblast",
"3" = "Epithelial cell",
"4" = "Endothelial cell",
"5" = "T cell",
"6" = "Fibroblast",
"7" = "Smooth muscle cell"
)
levels(dataset)
dataset[["cell_types"]] <- Idents(dataset)
15. 将重命名的细胞簇可视化为 UMAP 图上的注释（图 4）。
    DimPlot(dataset, group.by = "cell_types", raster = FALSE, label = FALSE)
16. 在一系列 UMAP 图上可视化表 1 中顶部（粗体）簇标记基因的定位（图 5）。
    DefaultAssay(dataset) <- "RNA"
FeaturePlot(dataset, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4)
17. 在点图上可视化顶部聚类标记 DEG，按原始聚类编号分组（补充图 6）。
    DotPlot(dataset, group.by = "seurat_clusters", features = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "blue", mid = "grey", high = "red")
18. 在点图上可视化顶部簇标记 DEG，按带注释的细胞类型分组（图 6）。
    DotPlot(dataset, group.by = "cell_types", features = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "blue", mid = "grey", high = "red")
19. 要执行时间序列分析，首先简化数据集以删除空间分量。对于时程分析，使用名为“DPW”的新元数据变量将伤口时间/空间注释分组为伤口后总天数 （DPW）。
    Idents(dataset) <- "time_space"
dataset[["DPW"]] <- Idents(dataset)
Idents(dataset) <- "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) <- levels(dataset)
dataset <- RenameIdents(dataset, new.cluster.ids)
dataset[["DPW"]] <- Idents(dataset)
Idents(dataset) <- "DPW"
20. 在 UMAP 图上可视化新的伤口时间过程分组（补充图 7）。
    DimPlot(dataset, group.by = "DPW", raster = FALSE, label = FALSE)
21. 生成表格，显示每个 DPW 中出现每种类型的细胞数。
    ​table(dataset$DPW, dataset$cell_types)
    1. 可选步骤：要获取伤口时程组的 DEG 列表，请将它们分配给变量，并将输出另存为分隔文本文件。
       Idents(dataset) <- "DPW"
Cell_DPW_markers <- FindAllMarkers(dataset, 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）：
    pt1 <- 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 = "fill", width = 0.5) +
xlab("Sample") + RotatedAxis() +
ylab("Proportion") +
theme(legend.title = element_blank())
23. 可视化每个DPW中细胞类型的比例（ 图7）。
    pt2 <- 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 = "fill", width = 0.5) +
xlab("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
24. 将数据集 Seurat 对象作为 RDS 文件保存到工作目录中。
    saveRDS(dataset, "dataset_post_Method3.rds")

4. 使用修拉分析细胞亚型

注意： 单细胞分析的强大功能能够发现和分析上述主要细胞类型中的罕见亚型。这个例子侧重于成纤维细胞，成纤维细胞最初聚集成两个修拉簇，然后组合成一个类别。该协议的这一部分专门关注成纤维细胞，排除所有其他细胞类型，并探索它们在伤口愈合过程中的身份和时间特性。作为可选步骤，如果需要，将保存的 RDS 文件加载为 Seurat 对象。

dataset <- readRDS("dataset_post_Method3.rds")

1. 根据成纤维细胞身份对原始数据集进行子集化。
   标识（数据集） <- “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）将DEG过滤到顶级成纤维细胞亚型标记中。
9. 通过从 DEG 文本文件中复制三个簇中每个簇的前 5 个成纤维细胞亚型标记，将自定义基因列表定义为变量。
   FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
10. 通过调用点图的特征参数中的变量，可视化仅成纤维细胞数据集中列表中的基因（图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. 通过调用点图的特征参数中的变量，可视化原始单细胞数据集中列表中的基因（补充图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("Sample") +
ylab("Proportion") +
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("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
14. 将数据集 Seurat 对象作为 RDS 文件保存到工作目录中。
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

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

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

在这里，使用先前发表的研究² 中的基因列表，其中使用来自整个愈合连续体的批量 RNA 测序样本鉴定伤口愈合阶段特异性基因。基因列表被保存到一个制表符分隔的文本文件（补充文件 2：JoVE_PhaseSpecificGenes.txt）中，该文件现在可以下载到工作目录中，并用于生成识别三个主要愈合阶段的基因列表。

1. 通过读取 TEXT 文件将基因列表加载到变量中。
   PhaseSpecificGenes <- read_delim("JoVE_PhaseSpecificGenes.txt", delim="\t", col_names = T)
2. 将列分成单独的基因列表变量，并将基因更改为第一个字母大写的小鼠名称。
   PS_Inflammatory <- PhaseSpecificGenes[1]
PS_Inflammatory_ms <- lapply(PS_Inflammatory, str_to_sentence)
PS_Proliferative <- PhaseSpecificGenes[2]
PS_Proliferative_ms <- lapply(PS_Proliferative, str_to_sentence)
PS_Resolution <- PhaseSpecificGenes[3]
PS_Resolution_ms <- lapply(PS_Resolution, str_to_sentence)
   注意： （可选）如果需要，将保存的 RDS 文件作为 Seurat 对象加载。
   dataset <- readRDS("dataset_post_Method3.rds")
3. 使用基因列表作为模块，根据愈合的三个阶段对数据集中的每个细胞进行评分。
   dataset <- AddModuleScore(
object = dataset,
features = PS_Inflammatory_ms,
ctrl = 100,
name = 'Inflammatory'
)
dataset <- AddModuleScore(
object = dataset,
features = PS_Proliferative_ms,
ctrl = 100,
name = 'Proliferative'
dataset <- AddModuleScore(
object = dataset,
features = PS_Resolution_ms,
ctrl = 100,
name = 'Resolution'
)
4. 可视化每个细胞类别的聚合模块分数，包括 DPW 和主要细胞类型（图 10）。
   DotPlot(dataset, group.by="DPW", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
DotPlot(dataset, group.by="cell_types", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")

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

注意： 另一种有用且被广泛引用的单细胞数据集分析方法是推断细胞间相互作用。在此工作流程中，使用包 CellChat，它通过分析细胞群之间的差异配体-受体相互作用来推断细胞间通讯²²。最近，CellChat 的开发人员发布了一个详细的分步协议，用于其通用用途²⁴，这对于用户来说是一个极好的资源，因为他们完成以下工作流程并将其应用于他们的数据集。例如，以下工作流程比较了伤口后 1 天与 14 天伤口中所有主要细胞的相互作用 （DPW）。每个步骤都没有详细描述，因为所有步骤都已在官方 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)
   成纤维细胞显着增加了 D1 和 D14 DPW 之间的相互作用。
5. 显示所有重要的推断细胞间通讯途径的列表。
   cellchat_D1@netP$pathways
cellchat_D14@netP$pathways
   胶原蛋白途径是 D1 和 D14 DPW 上的重要途径之一。
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. 可视化胶原信号通路与成纤维细胞作为源细胞的相互作用（补充图17）。
   注意： cellchat 对象中的细胞类型按照它们在原始 Seurat 对象中分配的顺序列为 ID：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. 以成纤维细胞为源细胞，可视化胶原信号通路中每个配体受体对的贡献。
    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. 使用气泡图可视化与第 1 天相比，第 14 天以成纤维细胞为源细胞的胶原信号通路中单个配体-受体对的相对贡献（补充图 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 = "Increased signaling in D14", 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 = "Decreased signaling in D14", 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. 通过组合多个单细胞数据集进行综合分析的示例

注意： 单细胞数据集通常被分成多个文件，因为它们是分批或分组测序的。该工作流程显示了如何集成伤口愈合数据集 ²⁰ 的五批中的两批。当前的数据集集成方法由以下 Seurat 小插图描述：

scRNA-seq 集成简介：
https://satijalab.org/seurat/articles/integration_introduction

Seurat v5 中的综合分析：

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

注意： 有许多单细胞数据集集成方法，每种方法都有自己的优点和缺点。有关详细信息，请参阅积分方法的综合基准²⁵。对于用户来说，在依赖任何一种集成方法之前阅读所有相关文档非常重要。

1. 对另一批胡等人数据集²⁰ 重复方法 2 中的所有步骤。在以下协议中，使用批次 #3。包括补充 R 脚本文件，可用于处理批次 #3（补充文件 3：JoVE_Rscript_b3。请记住为数据集创建和使用新变量---在下面的代码中，对数据集批次 #3 使用“dataset_b3”。
   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 注释，然后对合并的数据集执行标准 Seurat 工作流程，如方法 3 中所述。
   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")

Results

从方法#2开始，该协议演练了在单细胞伤口愈合数据集上加载和执行质量控制步骤的步骤。创建修拉对象（步骤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 函数来识别数据集中可能的双峰，并为每个单元格分配一个名为“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）是一个简化的表格，显示了每个修拉簇的前 5 个上调和表达基因。然后，该协议描述了如何使用名为EnrichR的基于网络的功能富集分析工具根据顶级簇标记基因（步骤3.11-3.12）识别假定的细胞类型，以及这些步骤的代表性结果（图3）是EnrichR输出的裁剪屏幕截图，显示了八个细胞簇中每个细胞簇的最高富集细胞类型。然后，该协议根据其最丰富的细胞类型注释，为各自修拉簇中的所有细胞分配一个名为“cell_types”的新元数据标签（步骤 3.14）。在步骤3.15中，执行维度图功能，将重命名的细胞簇可视化为UMAP图上的细胞类型注释，该步骤的代表性结果（图4）表明数据集中的所有细胞都聚集在主要颜色编码的细胞类型周围。在步骤3.16中，使用特征图功能在一系列UMAP图上可视化顶部簇标记基因的定位（来自补充表2），并得出代表性结果（图5）是UMAP图的网格，显示了顶级细胞标记基因在其各自的主要细胞类型簇位置内的高表达。在步骤3.17和3.18中，执行点图功能以可视化细胞中顶部簇标记基因的相对表达水平，首先按其原始修拉簇数分组（步骤3.17），然后按带注释的细胞类型标记分组（步骤3.18）。这些步骤的代表性结果证实了顶级细胞标记基因仅在其各自的修拉簇中高水平表达（补充图6）并且仅在它们各自的主要细胞类型（图6).该协议的下一步将原始基于时空蛋白质的标记简化为严格的时间注释，根据细胞起源的伤后天数（DPW）来识别细胞。在步骤3.20中，执行维度图功能，将单元格可视化为UMAP图上的DPW注释，该步骤的代表性结果（补充图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）显示数据集中的成纤维细胞聚集在3种颜色编码的细胞亚型周围。根据成纤维细胞数据集的DPW注释（步骤4.6）可视化，产生UMAP图（补充图10），显示数据集中的成纤维细胞根据其DPW注释分布在整个数据集中。该协议接下来描述如何获取差异表达基因的列表并将其保存到文本文件中（步骤4.7），在Excel中打开数据表并执行各种过滤步骤，以获得每个细胞簇的排名靠前的簇标记（步骤4.8），并分配一个新变量，列出名为“FB_type_marker”的顶级成纤维细胞标记基因（步骤4.9）。在步骤4.10中，使用点图函数通过调用特征参数中的“FB_type_marker”变量来可视化仅成纤维细胞数据集中列表中的基因，该步骤的代表性结果（图9）是点图，确认成纤维细胞亚型标记仅在其各自的聚类类别（顶部）中高表达，但在整个DPW类别（底部）中均匀分布。在步骤4.11中，调用相同的特征变量来可视化整个伤口愈合数据集中的成纤维细胞标记基因，代表性结果（补充图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 数据集子集为受伤后的两个时间点，一个在炎症阶段（第 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）是弦图，显示了每个时间点所有细胞类型之间推断的胶原通路信号相互作用。正如预期的那样，这些结果表明成纤维细胞是胶原信号通路的主要来源细胞，尽管与 D14 相比，D1 时的信息流更多地局限于免疫细胞。为了在细胞间相互作用中将成纤维细胞作为源细胞，步骤 6.9 重复弦图函数添加源细胞参数，代表性结果（补充图17）是弦图，显示了推断的胶原蛋白通路信号传导在每个时间点与成纤维细胞作为源细胞的相互作用。在步骤6.10中，执行两个功能来可视化每个配体-受体对在胶原信号通路中以成纤维细胞为源细胞的贡献，一个使用气泡图（步骤6.10.1），另一个使用弦图（步骤6.10.2）。代表性结果显示了每个配体-受体对在D1（左）和D14（右）时间点以成纤维细胞为源细胞的胶原通路信号传导中的推断贡献，使用两个气泡图（补充图18） 和弦图 （补充图19).这些结果表明，在 D1 时，来自成纤维细胞的胶原蛋白途径仅限于中性粒细胞和巨噬细胞，以 Cd44 和 Sdc4 受体为主，但在 D14 中，其他细胞通过包括整合素在内的多种受体充当接收者。为了关注在成纤维细胞相互作用中显示出强大强度的 Col1a1-Cd44 配体-受体相互作用，设置一个参数（步骤 6.11），然后在步骤 6.12 中使用弦图函数来可视化所有细胞类型之间的这种特定配体-受体相互作用，并具有代表性的结果（补充图20）是弦图，显示了在D1（左）和D14（右）时间点推断出的所有细胞类型之间Col1a1-Cd44配体-受体相互作用。这些结果表明，虽然在 D1 中这种相互作用仅限于成纤维细胞作为源细胞，但在 D14 中，巨噬细胞和平滑肌细胞也充当源细胞。接下来，该协议描述了如何通过首先合并D1和D14 CellChat对象来执行差分CellChat分析（步骤6.13）。在步骤6.14中，执行比较相互作用功能以可视化两个伤口愈合时间点之间细胞-细胞相互作用的总数和相对强度，以及代表性结果（补充图21）是生成的条形图，显示了包括D1和D14伤口的细胞中推断的相互作用的总数（左）和强度（右），D14中的相互作用数量较多，而D1中相互作用的相对强度较高。在步骤6.15和6.16中，使用两个函数来可视化伤口从第1天过渡到第14天时每种细胞类型之间的不同细胞间相互作用强度，并具有各自的代表性结果，第一个是圆图（步骤6.15， 补充图22），第二个是热图（步骤 6.16， 补充图23），其中与 D1 相比，D14 中的相互作用增加以红色显示，减少的相互作用以蓝色显示。正如预期的那样，中性粒细胞和巨噬细胞介导的相互作用在 D1 中增加，成纤维细胞介导的相互作用在 D14 中增加。在步骤6.17中，使用排名函数创建一个图，该图与D1相比，在D14与D1相比，单个途径对D14时与成纤维细胞作为源细胞的细胞间相互作用的相对贡献进行排名，以及代表性结果（补充图24）显示了生成的秩图，其中D1在顶部以红色表示，D14在底部以蓝色表示，其中几个途径仅在D1或D14中表示，许多其他途径显示激活梯度。最后，在步骤6.18中，使用两个气泡图函数来显示D14与D1相比，以成纤维细胞为源细胞的胶原信号通路中单个配体-受体对的相对贡献，并具有相应的代表性结果（补充图25）显示 D14 中信号对增加（左）和减少（右）在 x 轴上的许多细胞间相互作用中与 D1 相比。正如预期的那样，与 D1 伤口相比，成纤维细胞在 D14 伤口中多个受体细胞之间的传出配体-受体对相互作用要增加得多，D1 伤口在炎症阶段的通讯更局限于中性粒细胞和巨噬细胞。

从方法 #7 开始，该协议演练了使用 Seurat 集成两个单细胞伤口愈合数据集的步骤。该协议首先描述了合并两批已发布的单细胞数据集并将标准 Seurat 工作流程应用于合并数据集的步骤（步骤 7.1-7.4）。在步骤7.5中，使用维度图功能根据合并但尚未集成的伤口愈合数据集的聚类和批号对UMAP图进行可视化。这一步的代表性结果（补充图26）是UMAP图，该图可视化了修拉簇（左）和批号（右）的分布，表明在数据集成之前，这两个数据集似乎没有任何显着的批次效应。然后，该协议使用 RPCA 方法和集成数据集的后续 Seurat 工作流程执行数据集成（步骤 7.7-7.8）。在步骤7.9中，使用维度图功能根据集成伤口愈合数据集的聚类和批号对UMAP图进行可视化。这一步的代表性结果（补充图27）是UMAP图，它可视化了修拉簇（左）和批号（右）的分布，表明现在两个批次在不同簇中的重叠更大。结果还表明，在数据整合后出现了一个额外的集群，这可能表明在控制数据批次的技术效应后，识别潜在重要细胞亚型的能力增强。

图 1：UMAP 图显示了数据集中的所有单元格围绕 8 个主要颜色编码的聚类组聚类。 从运行 Windows（左）和 MacOS（右）的计算机获得的结果。该图对应于步骤 3.5。 请点击此处查看此图的大图。

图 2：UMAP 图显示数据集中的所有单元格根据其时间/空间原点展开，根据时间/空间注释没有明显的聚类。 该图对应于步骤 3.6。 请点击此处查看此图的大图。

图 3：EnrichR 输出的裁剪屏幕截图，显示了每个细胞簇的顶级富集细胞类型。 该图对应于步骤 3.13。 请点击此处查看此图的大图。

图 4：UMAP 图显示了数据集中聚集在主要颜色编码细胞类型周围的所有细胞。 该图对应于步骤 3.15。 请点击此处查看此图的大图。

图 5：UMAP 图网格显示主要细胞类型簇中顶级细胞标记基因的高表达。 该图对应于步骤 3.16。 请点击此处查看此图的大图。

图 6：点图确认顶级细胞标记基因仅在其各自的主要细胞类型中高水平表达。 该图对应于步骤 3.18。 请点击此处查看此图的大图。

图 7：显示每个 DPW 类别中主要细胞类型的相对数量的比例图。 该图对应于步骤 3.23。 请点击此处查看此图的大图。

图 8：UMAP 图显示数据集中的成纤维细胞聚集在 3 种颜色编码的细胞亚型周围。 该图对应于步骤 4.5。 请点击此处查看此图的大图。

图9：点图确认成纤维细胞亚型标记仅在其各自的簇类别中高表达，但在整个DPW类别中均匀分布。 该图对应于步骤 4.10。 请点击此处查看此图的大图。

图 10：点图显示了每个 DPW 和每个主要细胞类型中主要愈合阶段模块在细胞中的平均表达。 该图对应于步骤 5.4。 请点击此处查看此图的大图。

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

补充图2：散点图显示有许多细胞具有较大的线粒体含量，这与低RNA计数相关---这些细胞是死亡或垂死的细胞。 该图对应于步骤 2.10。 请点击此处下载此图。

补充图 3：散点图显示，检测到的基因分布和每个细胞的线粒体 RNA 百分比现在更加正常，为稳健的下游分析扫清了道路。 该图对应于步骤 2.12。请点击此处下载此图。

补充图 4：小提琴图显示有许多细胞具有相对较高的双峰评分，并且 0.25 看起来是一个自然截止值，高于该值有一组可能的双峰。 该图对应于步骤 2.15。 请点击此处下载此图。

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

补充图6：点图确认顶级细胞标记基因仅在其各自的修拉簇中高水平表达。 该图对应于步骤 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天（红色）DPW时每种细胞类型之间的差异细胞间相互作用强度。 该图对应于步骤 6.15。 请点击此处下载此图。

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

补充图 24：排名图显示了第 1 天与第 14 天 DPW 时各个途径对成纤维细胞和其他细胞类型之间细胞间相互作用的相对贡献。 该图对应于步骤 6.17。 请点击此处下载此图。

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

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

补充图 27：UMAP 图显示数据集成后 Seurat 簇（左）和批号（右）的分布。 该图对应于步骤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 个上调和表达基因。请点击此处下载此表。

Discussion

在此协议中，RStudio 用于运行规定的代码行，从而能够使用 Seurat 对复杂的单细胞数据集进行基本分析。提出了几种与伤口愈合研究相关的方法，包括安装 R 编码环境、下载先前发布的单细胞伤口愈合数据集、执行关键质量控制步骤和标准单细胞分析工作流程，包括可视化、主要细胞类型注释、细胞亚型分析和使用 Seurat 进行综合分析，以及使用 CellChat 进行细胞间相互作用分析。

本文介绍的方法是使用 R 及其流行的开源科学包 Seurat²¹ 和 CellChat²² 进行单细胞分析的典型工作流程的简化小插曲。事实上，该工作流程只是使用复杂的单细胞伤口愈合数据集可以完成的分析类型的一个示例。对这种方法的可能修改几乎是无限的，唯一的限制是用户特定的科学探究。例如，用户可以根据他们可能希望对该数据集提出的研究问题更改一些关键参数，例如细胞类型和时间点。作者还希望用户感觉足够舒适，能够将此工作流程调整到他们自己感兴趣的单细胞数据集中;但是，使用此工作流程分析其他数据集时必须小心，因为每个实验都可能将技术和样品制备问题传播到数据本身。因此，在解释任何先前发布和重新分析的单细胞数据集的结果之前，用户必须阅读并理解所有实验细节。重要的是要记住，生物信息学工具是探索生物过程和假设生成的有力方法，并且任何对结果的关键生物学解释都必须在后续实验中得到验证。

在整个协议中，请注意，可能会对工作流程的特定区域进行重大修改以完成其他任务。但是，对工作流程的所有可能修改组合的细节超出了本手稿的范围。例如，用于细胞聚类的分辨率和用于 UMAP 分析的维度必然是主观的，这里介绍的工具既允许大规模分析（如此处对广泛定义的主要细胞类型所展示的那样）也允许非常具体的分析，这些分析可能需要将细胞亚聚类为更大数据集中更罕见的亚群。有关单细胞分析方法这方面的更多信息，以及有关单细胞分析管道中可能更改的所有其他参数的详细信息，作者建议用户访问修拉出版物^21,26 和网站 （https://satijalab.org/seurat/），该不断发展的工具的作者在其中提供了深入的解释、小插曲和教程。

本手稿介绍了单细胞转录组学文献中一些被引用和使用最多的工具，即 Seurat²¹ 和 CellChat²²，分别用于单细胞和细胞间相互作用分析。但是，存在其他工具以略有不同的方式执行类似的功能。对于单细胞数据集分析，有 Scran²⁷、Scater ²⁸ 和基于 Python 的 ScanPy²⁹，它们使用各种方法进行数据集集成²⁵。在该协议中，演示了细胞类型的手动注释，这依赖于用户解释簇细胞标记物富集的判断，但现在存在各种工具，可以自动分类细胞类型，例如 SingleR³⁰ 和 scGate³¹ 等。对于细胞-细胞通讯分析，CellChat 在该协议中得到了演示，但还有其他用于估计细胞-细胞通讯的工具，包括 CellPhoneDB³²、Cytotalk³³ 和在 LIANA（LIgand-receptor ANalysis 框架）共识框架³⁴ 中实施的其他配体受体数据库。所有生物信息学工具都是独一无二的，并具有自己的特点和可修改的参数。因此，用户在解释其使用产生的任何输出之前，仔细阅读每个工具的相关文档以了解其细微差别非常重要。最后，无论使用什么生物信息学工具，重要的是要记住，这些工具在不断发展，不同版本的软件包可能会产生不同的输出。

在 R 中，语法至关重要，错位的标点符号、引号、括号，甚至大写错误的字母都会导致错误。因此，用户在键入代码时注意细节，在复制代码行时要特别小心，以使其适应新的科学问题和数据集，这一点至关重要。为了解决可能遇到的具体错误，作者建议只需将错误消息复制粘贴到用户最喜欢的网络搜索引擎中，然后浏览生物信息学论坛（如 GitHub 和 Stack Overflow）的结果，因为最常见的错误可能已经由知识渊博的高级用户回答。在一些论坛中，最成功的答案会被其他用户“赞成”，他们认为解决方案最适合该问题。用户必须小心，不要简单地将他们在 Internet 上找到的代码行复制粘贴到他们自己的计算机中（特别是当解决方案要求在 R 编程语言之外更改系统设置时），因为此类程序可能是恶意的。一种令人兴奋的新兴解决编码错误的方法是使用强大的生成式大型语言 AI 模型，例如 OpenAI 的 ChatGPT、Microsoft 的 Copilot 或 Google 的 Gemini。事实证明，这些模型对于一般软件工程，特别是故障排除特别有用。为此，用户可以在向聊天机器人提供有关用户代码意图的简单提示后复制粘贴其代码的整行。通常需要注意的是，这些模型并非万无一失，用户可能必须尝试多个提示才能生成适合解决问题的答案。

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
笔记本电脑或台式电脑	不适用	不适用	运行 Windows 或 MacOS
R	不适用	版本 4.4.1	可从 https://cran.rstudio.com/ 免费下载
Rstudio 工作室	Posit 软件、PBC	版本 2024.09.0	可从 https://posit.co/download/rstudio-desktop/ 免费下载
办公Excel公司	Microsoft	任何版本	用于分析表数据
互联网浏览器	不适用	不适用	用于导航到网站
R 软件包	存储库	版本
开发工具	CRAN	2.4.5
阅读xl	CRAN	1.4.3
OpenXLSX	CRAN	4.2.7.1
整洁宇宙	CRAN	2.0.0
sc自定义	CRAN	2.1.2
生物经理	生物导体	1.30.25
NMF的	生物导体	0.28
复杂热图	生物导体	2.20.0
生物邻居	生物导体	1.22.0
单细胞实验	生物导体	1.26.0
循环	生物导体	0.4.16
边缘R	生物导体	4.2.1
scDbl查找器	生物导体	1.18.0
修拉	CRAN	5.1.0
手机聊天	Github的	2.1.2