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

摘要

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

摘要

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

引言

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

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

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

方案

NOTE: In the following workflows detailing seven bioinformatics methods, all steps of the protocols are accompanied by their respective blocks of code that should be run directly on the user's own RStudio interface in the order that they are listed. To make this protocol as user-friendly as possible, an R script file is included (Supplementary File 1: JoVE_Rscript.R), which may be loaded directly into the user's RStudio session, so that each line of code may simply be run. This avoids the user having to type or copy and paste the code from the protocol document, which could introduce errors. All protocol instructions are also included in the R script file in the form of comments, indicated by a hashtag '#' symbol at the beginning of each comment line.

1. Installing R, RStudio, and the required R packages for the single-cell analysis workflow

1. Download and install R (version 4.4.1) to the computer. Use the link that corresponds to the computer's operating system.
   1. If using a computer running Microsoft Windows, use this link: https://cran.rstudio.com/bin/windows/base/
   2. If using a computer running MacOS, use this link: https://cran.rstudio.com/bin/macosx/
2. Install the latest version of RStudio on the computer. Click on the following link and follow the instructions:
   https://www.rstudio.com/products/rstudio/download/#download
3. Install the Rtools (version 4.4), which will allow R to compile certain packages. Click on the following link and follow the instructions:
   1. If using Windows, use this link: https://cran.rstudio.com/bin/windows/Rtools/
   2. If using MacOS, use this link:
      https://mac.r-project.org/tools/
4. Set the local working directory; this is the folder on the computer where all files will be loaded in from and saved into. Set the working directory by selecting Session in the RStudio menu bar and clicking Set Working Directory > Choose Directory and selecting the desired folder.
   1. If using a Windows computer, use the following command to set the working directory. Change [Directory] in the following line of code to the actual directory structure. Be mindful that the directory delimiter in R is the character "/"
      setwd("C:/[Directory]")
   2. If using a MacOS computer, the following command will also set the working directory. Change [Directory] in the following line of code to the actual directory structure. Be mindful that the directory delimiter in R is the character "/"
      setwd("~/[Directory]")
   3. At any point during an R session, check the working directory using the following line of code:
      getwd()
   4. In RStudio, visually explore the working directory structure, including all files and folders contained therein, within the right-hand window in the Files tab. To navigate RStudio's Files explorer to the working directory, click the gear icon -> Go To Working Directory.
5. Install the following packages from the R package repository CRAN, which are necessary dependencies for the protocol. To install these packages, run the following commands.
   install.packages("devtools")
install.packages("readxl")
install.packages("openxlsx")
install.packages("tidyverse")
install.packages("scCustomize")
   NOTE: During installation of R packages, it is normal for various windows to appear and disappear. If a window appears asking to compile a package, click YES. If a window appears asking to restart R prior to installing the package, click NO.
6. Install the following packages from the curated R package repository Bioconductor
   (https://bioconductor.org/), which are necessary dependencies for the protocol. To install these packages, run the following commands. 
   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 the following packages, which are required for the workflow described in this manuscript.
   install.packages("Seurat")
devtools::install_github("jinworks/CellChat")
8. Load in each package to confirm that the installations were successful. In case any of the packages results in a "package not found" error, re-install it using the appropriate code above.
   library(readxl)
library(openxlsx)
library(tidyverse)
library(scCustomize)
library(edgeR)
library(scDblFinder)
library(Seurat)
library(CellChat)

2. Loading in a single-cell wound healing dataset and performing quality control steps

NOTE: For this bioinformatics workflow, a re-analysis of a previously published spatio-temporal single-cell skin wound healing experiment is conducted²⁰. The dataset files are stored in the curated NCBI Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/).

1. Navigate to the dataset files from GEO using the accession number GSE204777. Use the following link and review the study's experimental design: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204777
2. Within the GEO repository page, there are five individual batches of data obtained from five sequencing lanes. Click on the first dataset, entitled GSM6190913. The following is a direct link to the sample: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM6190913
3. Scroll down to the bottom of the page and download the following three files, using either the ftp or html links. Within the computer's file explorer, move these three files into a directory called b1. Make sure that the b1 folder is located within the working directory set in Step 1.4.
   File name: GSM6190913_b1_barcodes.tsv.gz / file size: 18.5 Mb
   File name: GSM6190913_b1_features.tsv.gz / file size: 254.1 Kb
   File name: GSM6190913_b1_matrix.mtx.gz / file size: 151.2 Mb
4. Get the directory information for the single-cell sequencing files downloaded in step 2.3.
   b1_data_dir <- file.path(getwd(), "b1")
5. Load in the single-cell sequencing files. The gene.column parameter specifies the gene/feature naming used. In this case, use gene.column = 2 for gene symbols (gene.column = 1 is for Ensembl gene names).
   b1 <- Read10X_GEO(data_dir = b1_data_dir, gene.column = 2)
   ​NOTE: Most single-cell datasets do not have additional multiplexed data, so the 10x files generated using this step would not have multiple layers. For the current multiplexed dataset, proceed with step 2.6. For datasets without multiplexing data, skip to step 2.7.
6. Demultiplex the single-cell dataset using spatio-temporal barcodes.
   1. For the working dataset, separate the gene expression and HTO (multiplexing) data.
      dataset_rna_counts <- b1$GSM6190913_b1_$`Gene Expression`
​dataset_barcodes <- b1$GSM6190913_b1_$`Antibody Capture`
   2. Create a Seurat object using the gene expression data, while immediately filtering out genes expressed in less than 5 cells and cells with less than 200 genes detected.
      dataset <- CreateSeuratObject(counts = dataset_rna_counts, min.cells = 5, min.features = 200)
   3. Create a gene expression dataset as a layer and generate a list of cells and barcodes common to both assays.
      dataset_rna_counts2 <- LayerData(dataset, data = "RNA")
​dataset_joint.barcodes <- intersect(colnames(dataset_rna_counts2), colnames(dataset_barcodes))
   4. Subset gene expression and HTO counts by joint cell barcodes.
      dataset_rna_counts2 <- dataset_rna_counts2[, dataset_joint.barcodes]
​dataset_barcodes2 <- as.matrix(dataset_barcodes[, dataset_joint.barcodes])
   5. Confirm that the HTO have the expected barcode names.
      ​rownames(dataset_barcodes2)
   6. Create a new assay to store barcode information and add this assay to the previously created Seurat object.
      dataset_barcode_assay <- CreateAssayObject(counts = dataset_barcodes2)
​dataset[["barcodes"]] <- dataset_barcode_assay
   7. Validate that the object now contains multiple assays.
      DefaultAssay(dataset)
   8. Normalize barcode data and perform demultiplexing via the HTODemux function. This method is described in detail in the following Seurat vignette: https://satijalab.org/seurat/articles/hashing_vignette
      ​dataset <- HTODemux(dataset, assay = "barcodes", positive.quantile = 0.99)
   9. Group cells based on global classification results and remove cells with no barcode classification.
      Idents(dataset) <- "barcodes_classification.global"
​dataset <- subset(dataset, idents = "Negative", invert = TRUE)
   10. Group cells based on the maximum HTO signal.
       Idents(dataset) <- "barcodes_maxID"
   11. Visualize the distribution of detected genes in cells per their multiplexed barcodes (Supplementary Figure 1).
       VlnPlot(dataset, features = "nFeature_RNA", pt.size = 0.1, log = TRUE)
   12. Rename the barcodes to their actual wound time (days post-wounding) and space (2-8 mm) assignments (taken from the original manuscript) and assign them to a new metadata variable called 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. For datasets without multiplexing data: Create a Seurat object, while immediately filtering out genes expressed in less than 5 cells and cells with less than 200 genes detected.
   dataset <- CreateSeuratObject(counts = b1, min.cells = 5, min.features = 200)
8. Switch to analyzing the gene expression assay of the dataset.
   DefaultAssay(dataset) <- "RNA"
Idents(dataset) <- "data"
9. As an important quality control step, calculate the percentage of mitochondrial genes in each cell and assign it as a metadata variable. This method is described in detail in the following Seurat vignette: https://satijalab.org/seurat/articles/pbmc3k_tutorial#qc-and-selecting-cells-for-further-analysis
   dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^mt-")
10. Visualize the distribution of detected genes, numbers of RNA and the mitochondrial percentage in all cells (Supplementary Figure 2).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
11. There are a number of cells with large mitochondrial content, which correlates with low RNA counts; these are dead or dying cells. Remove these low-quality cells from the dataset using a fair cutoff. In this case, use the values of the original dataset described in the previously published study²⁰, in which cells with more than 25% mitochondrial genes are removed.
    dataset <- subset(dataset, subset = nFeature_RNA > 200 & percent.mt < 25)
12. Visualize the distribution of detected genes, numbers of RNA, and the mitochondrial percentage in all cells after removing low-quality cells (Supplementary Figure 3).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
13. As an additional important quality control step, detect likely doublets in the dataset. These are cells that were joined during droplet sequencing and will thus result in gene expressions that are not at the single-cell level. To counter this problem, several tools have been developed. It is important to note that each tool makes certain assumptions about single-cell data, therefore it is important for the user to read all relevant documentation before using any one tool. In this workflow, implement a method called scDblFinder²³. Please note that this method uses an algorithm that imposes a fixed expected doublet rate. Use the following commands to run the scDblFinder pipeline.
    DefaultAssay(dataset) <- "RNA"
sce <- scDblFinder(GetAssayData(dataset, assay="RNA", slot="counts"), samples=Idents(dataset))
    NOTE: Multiplexed single-cell datasets such as this one may also be screened for doublets by removing cells that express multiple barcodes. This workflow did not demonstrate this method because most single-cell datasets do not have this unique feature. Instead, a more generalizable pipeline for screening for doublets is shown using the scDblFinder method.
14. Assign the doublet score to a new metadata variable.
    dataset$scDblFinder.score <- sce$scDblFinder.score
15. Visualize the distribution of the doublet score in all cells (Supplementary Figure 4).
    VlnPlot(dataset, features = "scDblFinder.score", raster=FALSE, pt.size=0.5)
16. Remove the cells above the 0.25 double score threshold. This threshold was chosen based on the violin plot generated above, which showed that most cells in the dataset could be assigned to either very high or very low double scores, and 0.25 is a reasonable cutoff for this dataset that would remove the vast majority of likely doublets without removing many unlikely doublets.
    dataset <- subset(dataset, scDblFinder.score < 0.25)
17. Save the dataset Seurat object as an RDS file into the working directory.
    saveRDS(dataset, "dataset_post_Method2.rds")

3. Analyzing a single-cell wound healing dataset using Seurat

NOTE: (Optional step) If starting the workflow here, load in the saved RDS file as a Seurat object.

dataset <- readRDS("dataset_post_Method2.rds")

1. Perform the standard Seurat workflow for single-cell dataset normalization, scaling, and Principal Component Analysis (PCA). This standard workflow is described in the following Seurat vignettes:
   Seurat - Guided Clustering Tutorial: https://satijalab.org/seurat/articles/pbmc3k_tutorial
   Seurat Command List: 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. Visualize the amount of dataset variation with respect to the first 50 PCA dimensions (Supplementary Figure 5).
   ElbowPlot(dataset, reduction = "pca", ndims = 50)
   Much of the major variation occurs within the first 13 dimensions.
3. Perform cell clustering of the dataset using a set PCA dimension range of 1-13 and a set resolution of 0.1.
   dataset <- FindNeighbors(dataset, verbose = TRUE, dims = 1:13)
dataset <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
   NOTE: This workflow focuses on large scale differences in cell types. Therefore, it uses a fairly conservative parameter for the PCA dimension range, wherein the first 13 dimensions are shown to represent the vast majority of variation within the dataset. For discrimination of cells into smaller and rarer subtypes, the user can use a higher number of dimensions for downstream analysis since those rare cell subtypes likely account for lower levels of dataset variation. The resolution parameter ranges from 0 to 1, and it determines the magnitude of categorical separation imposed on the dataset. The setting of this parameter is dependent on the user's research question. For clustering of cells into numerous small and rare subtypes, use higher resolutions for downstream analysis. Since this workflow aims to explore more broad differences between major cell types, it uses a fairly small resolution value of 0.1, which is expected to cluster cells into fewer, larger groups. Use of the settings mentioned in step 3.3 will discriminate 8 unique cell clusters. These clusters are automatically assigned to a metadata variable called "seurat_clusters".
4. Perform UMAP dimensional reduction and finding neighbors analysis using the first 13 PCA dimensions. Add the seed number 123 to ensure the reproducibility of the resulting data projection.
   dataset <- RunUMAP(dataset, verbose = TRUE, dims = 1:13, seed.use = 123)
   NOTE: The UMAP algorithm is stochastic and introduces randomness into the dimensional reduction (see "Stability and Reproducibility" in https://cran.r-project.org/web/packages/umap/vignettes/umap.html). While using a consistent seed provides a "minimal level of reproducibility" to the algorithm, the resulting plot may still be slightly different from what is shown in the representative figures and downstream results. Testing has found that results will especially differ between computers running Windows and those running MacOS, likely due to different implementations of randomness in these operating systems.
5. Visualize the clustering of the cells on a UMAP plot (Figure 1).
   DimPlot(dataset, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
   NOTE: The randomness of the UMAP algorithm may generate slightly different plots, as shown in the figure showing alternate UMAP plots generated using the same code as above on computers running Windows and MacOS; note the slight differences in the shapes of the clusters. It is, therefore, imperative that the user saves and time-stamps all data and plots as they are generated, and that all downstream analyses on clusters be made carefully and with biological understanding in mind, as it is described below for cell type annotation.
6. Because the original experiment labels that refer to where and when the cells came from during wound healing are included, visualize the wound time/space annotation of the cells on a UMAP plot (Figure 2).
   DimPlot(dataset, group.by = "time_space", raster = FALSE, label = FALSE)
7. Generate a table of cell cluster association with the wound time/space annotation.
   table(dataset$time_space, dataset$seurat_clusters)
8. Determine the identities of the major cell types in the dataset. To do this, calculate the differentially expressed genes (DEG) between all clusters. Get DEG lists for clusters, assign them to a variable, and save the output as a delimited text file within the working directory.
   NOTE: This step is CPU intensive and can take a long time, depending on the user's hardware.
   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. Download and open the included dataset_cluster_markers.txt file in a spreadsheet (e.g., Excel) by copying the contents of the text file and using the Text Import Wizard to specify the comma delimiter and the identity of the gene name columns as Text. Indicating that the gene names are ‘Text’ is important, otherwise Excel will automatically convert certain gene names to dates, e.g. Sept7 changing to September 7..
10. In a spreadsheet, filter the results according to the following recommended parameters:
    1. Rank the avg_log2FC column from largest to smallest in order to arrange all rows according to decreasing log₂ fold changes (log2FC).
    2. Rank the cluster column from smallest to largest to arrange all rows according to increasing Seurat cluster numbers.
    3. Filter the avg_log2FC column for numbers greater than or equal to 2.5 to show only the most differentially expressed genes (DEGs) in the indicated cluster versus other clusters.
    4. Filter the pct.1 column for numbers greater than or equal to 0.4. This column refers to percentage of cells in the indicated cluster expressing the indicated gene (Cluster %) and setting the threshold at 0.4 means that only genes expressed in at least 40% of the cells in the indicated cluster are shown.
    5. Filter the pct.2 column for numbers less than or equal to 0.2. This column refers to the percentage of cells NOT in the indicated cluster expressing the indicated gene (Non-cluster %), and setting the threshold at 0.2 means that only genes expressed in at most 20% of the cells NOT in the indicated cluster are shown.
    6. Filter the p_val_adj column for numbers less than or equal to 0.01. This column refers to the adjusted P value or the false discovery rate (FDR), indicating the statistical strength of the indicated DEG, and setting the threshold at 0.01 means that only genes with an FDR < 0.01 are shown.
       NOTE: Supplementary Table 1 (JoVE_DEGs_cellMarkers.xlsx) contains the full output of ranked differentially-expressed genes used in step 3.10 of the protocol. Supplementary Table 2 shows the top 5 genes for each cluster in this analysis, with the bolded genes used for subsequent visualizations.
11. For unbiased cell-type annotation of clusters, use the EnrichR web-based enrichment analysis tool.
    Use the link: https://maayanlab.cloud/Enrichr/
12. Copy the lists of DEGs for each cluster into a separate EnrichR window, then click Analyze. The EnrichR tool runs the gene list through hundreds of curated databases and ranks every enriched term in each category.
13. For the purposes of cell type annotation, click the Cell Types tab above and focus on the top 5 enrichments in the three cell marker curated databases on the left-hand side (Figure 3):
    CellMarker 2024 (http://bio-bigdata.hrbmu.edu.cn/CellMarker/)
    Tabula Sapiens (https://tabula-sapiens-portal.ds.czbiohub.org/)
    PanglaoDB Augmented (https://panglaodb.se/)
14. Based on the enrichments of the DEGs in these databases, confirm the likely identity of the 8 clusters. Note that there are two clusters (2, 6) that enrich as fibroblasts; therefore, combine these clusters into single cell type annotations. Assign the cell type identities as labels to a new metadata variable called 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. Visualize the renamed cell clusters as annotations on a UMAP plot (Figure 4).
    DimPlot(dataset, group.by = "cell_types", raster = FALSE, label = FALSE)
16. Visualize the localization of the top (bolded) cluster marker genes from Table 1 on a series of UMAP plots (Figure 5).
    DefaultAssay(dataset) <- "RNA"
FeaturePlot(dataset, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4)
17. Visualize the top cluster marker DEGs on a dot plot, grouped by the original cluster numbers (Supplementary Figure 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. Visualize the top cluster marker DEGs on a dot plot, grouped by the annotated cell types (Figure 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. To perform time-series analyses, first simplify the dataset to remove the spatial component. For time-course analyses, group the wound time/space annotations into overall days post-wounding (DPW) with a new metadata variable called "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. Visualize the new wound time-course groupings on a UMAP plot (Supplementary Figure 7).
    DimPlot(dataset, group.by = "DPW", raster = FALSE, label = FALSE)
21. Generate tables showing how many cells of each type occur in each DPW.
    ​table(dataset$DPW, dataset$cell_types)
    1. Optional Step: To also get DEG lists for wound time-course groups, assign them to a variable, and save the output as a delimited text file.
       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. Convert the cell numbers into proportions per category to better understand the relative changes in cell type composition throughout the time course of healing. Visualize the proportion of DPW in each cell type (Supplementary Figure 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. Visualize the proportion of cell types in each DPW ( Figure 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. Save the dataset Seurat object as an RDS file into the working directory.
    saveRDS(dataset, "dataset_post_Method3.rds")

4. Analyzing cell subtypes using Seurat

NOTE: The power of single-cell analysis enables the discovery and analysis of rare subtypes within the major cell types that are analyzed above. This example focuses on fibroblasts, which initially clustered into two Seurat clusters before they were combined into a single category. This part of the protocol focuses on fibroblasts specifically, excluding all other cell types, and explores their identities and temporal properties during wound healing. As an optional step, if needed, load the saved RDS file as a Seurat object.

dataset <- readRDS("dataset_post_Method3.rds")

1. Subset the original dataset according to the fibroblast cell identity.
   Idents(dataset) <- "cell_types"
   dataset_fibroblast <- subset(dataset, idents = "Fibroblast")
2. Perform PCA on this smaller dataset and visualize the amount of dataset variation with respect to PCA dimensions (Supplementary Figure 9).
   dataset_fibroblast <- RunPCA(dataset_fibroblast, verbose = TRUE)
ElbowPlot(dataset_fibroblast, reduction = "pca", ndims = 50)
   NOTE: Much of the major variation occurs within the first 9 dimensions.
3. Perform cell clustering of the dataset using a set PCA dimension range of 1-9 and a set resolution of 0.1.
   dataset_fibroblast <- FindNeighbors(dataset_fibroblast, verbose = TRUE, dims = 1:9)
dataset_fibroblast <- FindClusters(dataset_fibroblast, verbose = TRUE, resolution = 0.1)
   NOTE: Using these settings, the algorithm discriminates 3 unique fibroblast clusters. These clusters are automatically assigned to a metadata variable called seurat_clusters.
4. Perform UMAP dimensional reduction and finding neighbors analysis using the first 9 PCA dimensions. Add the seed number 123 to ensure the reproducibility of the resulting data projection.
   dataset_fibroblast <- RunUMAP(dataset_fibroblast, verbose = TRUE, dims = 1:9, seed.use = 123)
5. Visualize the clustering of the cells on a UMAP plot (Figure 8).
   DimPlot(dataset_fibroblast, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
6. Visualize the wound time-course annotation of the cells on a UMAP plot (Supplementary Figure 10).
   DimPlot(dataset_fibroblast, group.by = "DPW", raster = FALSE, label = FALSE)
7. Get DEG lists for the three fibroblast subtypes and save them into a text file in the working directory.
   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. Follow similar steps as above (steps 3.9-3.10) in a spreadsheet to filter the DEGs into the top fibroblast subtype markers.
9. Define a custom gene list as a variable by copying the top 5 fibroblast subtype markers for each of the three clusters from the DEG text file.
   FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
10. Visualize the genes in the list in the fibroblast-only dataset by calling the variable in the features parameter of the dotplot (Figure 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. Visualize the genes in the list in the original single-cell dataset by calling the variable in the features parameter of the dotplot (Supplementary Figure 11).
    DotPlot(dataset, group.by="cell_types", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
12. Visualize the proportion of fibroblast subtypes in each DPW category (Supplementary Figure 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. Visualize the proportion of DPW fibroblast cells in each fibroblast subtype category (Supplementary Figure 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. Save the dataset Seurat object as an RDS file into the working directory.
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

5. Example of a follow-up analysis via module scoring

NOTE: One useful method for analyzing single-cell datasets is called module scoring. In this workflow, one can define a gene list according to previous knowledge and then calculate module scores, which can identify potential enrichments of the gene list within each cell. These scores can be averaged across cell annotations to reveal potential patterns of enrichments.

Here, use gene lists from a previously published study², where wound healing phase-specific genes were identified using bulk RNA-sequenced samples from across the healing continuum. The gene lists were saved into a tab-delimited text file (Supplementary File 2: JoVE_PhaseSpecificGenes.txt) that can now be downloaded into the working directory and used to generate gene lists that identify the three major healing phases.

1. Load in gene lists into a variable by reading the TEXT file.
   PhaseSpecificGenes <- read_delim("JoVE_PhaseSpecificGenes.txt", delim="\t", col_names = T)
2. Separate the columns into individual gene list variables and change the genes to mouse names that have their first letter capitalized.
   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)
   NOTE: (Optional) If needed, load the saved RDS file as a Seurat object.
   dataset <- readRDS("dataset_post_Method3.rds")
3. Use the gene lists as modules to score each cell in the dataset according to the three phases of healing.
   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. Visualize the aggregate module scores per cell category, including DPW and major cell types (Figure 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. Example of a follow-up analysis via CellChat

NOTE: Another useful and well-cited method for the analysis of single-cell datasets is inferring cell-cell interactions. In this workflow, use the package CellChat, which infers cell-cell communications by analyzing differential ligand-receptor interactions between cell groups²². Recently, the developers of CellChat published a detailed step-by-step protocol for its generalized use²⁴, and this is an excellent resource for users as they work through the following workflow and apply it to their datasets. As an example, the following workflow compares interactions of all major cells in wounds at 1 vs 14 days post-wounding (DPW). Each step is not described in great detail, since all steps have already been described in the official CellChat publication²⁴ as well as tutorials, linked here:

Inference and analysis of cell-cell communication using
CellChat: https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html

Comparison analysis of multiple datasets using CellChat:
https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html

Optional step: If needed, load in the saved RDS file as a Seurat object:

dataset <- readRDS("dataset_post_Method3.rds")

1. Subset the original dataset into two datasets according to the DPW annotation.
   Idents(dataset) <- "DPW"
dataset_D1 <- subset(dataset, ident = "D1")
dataset_D14 <- subset(dataset, ident = "D14")
2. Define the annotation by which CellChat will be performed --- in this case, use the major cell types.
   Idents(dataset_D1) <- "cell_types"
Idents(dataset_D14) <- "cell_types"
3. Create the CellChat objects and follow the typical CellChat workflow. See tutorials linked above as detailed references for each step.
   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. Visualize the incoming vs outgoing interaction strengths in all major cell types at each wound healing time point (Supplementary Figure 14).
   netAnalysis_signalingRole_scatter(cellchat_D1)
netAnalysis_signalingRole_scatter(cellchat_D14)
   Fibroblasts dramatically increase their interactions between D1 and D14 DPW.
5. Show the lists of all significant inferred cell-cell communication pathways.
   cellchat_D1@netP$pathways
cellchat_D14@netP$pathways
   The collagen pathway is one of the significant pathways on both D1 and D14 DPW.
6. Focus on the collagen signaling pathway and its interaction with fibroblasts.
   pathways.show <- c("COLLAGEN")
7. Visualize the collagen signaling pathway interactions between cell types using circle diagrams (Supplementary Figure 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. Visualize the collagen signaling pathway interactions between cell types using chord diagrams (Supplementary Figure 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. Visualize the COLLAGEN signaling pathway interactions with fibroblasts as source cells (Supplementary Figure 17).
   NOTE: The cell types in the cellchat objects are listed as IDs in the order that they were assigned in the original Seurat object: 1 = Macrophage, 2 = Neutrophil, 3 = Fibroblast, 4 = Epithelial cell, 5 = Endothelial cell, 6 = T cell, 7 = Smooth muscle cell.
   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. Visualize the contributions of each ligand-receptor pair in the COLLAGEN signaling pathway with fibroblasts as source cells.
    1. Using bubble plots (Supplementary Figure 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. Using chord diagrams (Supplementary Figure 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. Focus on the Col1a1-Cd44 ligand-receptor interaction within the COLLAGEN signaling pathway.
    ​LR.show <- "COL1A1_CD44"
12. Visualize the Col1a1-Cd44 ligand-receptor interactions between cell types using chord diagrams (Supplementary Figure 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. Perform differential CellChat analysis by generating a combined CellChat object.
    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. Visualize the total numbers and relative strengths of cell-cell interactions in the wound healing time points (Supplementary Figure 21).
    gg1 <- compareInteractions(cellchat_D14_v_D1, show.legend = F)
gg2 <- compareInteractions(cellchat_D14_v_D1, show.legend = F, measure = "weight")
gg1 + gg2
15. Visualize using a circle plot the differential cell-cell interaction strengths between each cell type as the wound transitions from day 1 to day 14 (Supplementary Figure 22).
    netVisual_diffInteraction(cellchat_D14_v_D1, weight.scale = T, measure = "weight")
16. Visualize using a heatmap the differential cell-cell interaction strengths between each cell type as the wound transitions from day 1 to day 14 (Supplementary Figure 23).
    netVisual_heatmap(cellchat_D14_v_D1, measure = "weight")
17. Visualize using a rank plot the relative contributions of individual pathways to cell-cell interactions with fibroblasts as source cells at day 14 vs day 1 ( Supplementary Figure 24).
    rankNet(cellchat_D14_v_D1, mode = "comparison", measure = "weight", sources.use = 3, targets.use = NULL, stacked = T, do.stat = TRUE)
18. Visualize using bubble plots the relative contributions of individual ligand-receptor pairs in the collagen signaling pathway with fibroblasts as source cells at day 14 compared to day 1 (Supplementary Figure 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. Just like Seurat objects, CellChat objects can be saved and opened as RDS files.
    saveRDS(cellchat_D1, "cellchat_D1.rds")
saveRDS(cellchat_D14, "cellchat_D14.rds")
saveRDS(cellchat_D14_v_D1, "cellchat_D14_v_D1.rds")
20. Optional step: CellChat objects can also be opened from RDS files.
    cellchat_D1 <- readRDS("cellchat_D1.rds")
cellchat_D14 <- readRDS("cellchat_D14.rds")
cellchat_D14_v_D1 <- readRDS("cellchat_D14_v_D1.rds")

7. Example of an integrative analysis by combining multiple single-cell datasets

NOTE: Single-cell datasets are often separated into multiple files because they were sequenced in batches or groups. This workflow shows how to integrate two of the five batches of the wound healing dataset ²⁰. The current methods for dataset integration are described by the following Seurat vignettes:

Introduction to scRNA-seq integration: 
https://satijalab.org/seurat/articles/integration_introduction

Integrative analysis in Seurat v5:

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

Note: There are numerous single-cell dataset integration methods, each with its own strengths and weaknesses. For details, see the comprehensive benchmark of integration methods²⁵. It is important for the user to read all relevant documentation before relying on any one integration method.

1. Repeat all steps in Method 2 for another batch of the Hu et al. dataset²⁰. In the following protocol, batch #3 is used. The supplementary R script file is included and can be used to process batch #3 (Supplementary File 3: JoVE_Rscript_b3.R). Remember to create and use a new variable for the dataset --- in the code below, use "dataset_b3" for dataset batch #3.
   1. Optional Step: If needed, open the two datasets as Seurat objects from their RDS files saved in the working directory:
      dataset <- readRDS("dataset_post_Method2.rds")
dataset_b3 <- readRDS("dataset_b3_post_Method2.rds")
2. Assign a new variable to each dataset named "batch" to label the dataset of origin in subsequent analyses.
   dataset@meta.data$batch <- "b1"
dataset_b3@meta.data$batch <- "b3"
3. Perform Seurat merging of the two datasets, adding batch-based cell ID annotations, and then perform the standard Seurat workflow for the merged dataset, as described in Method 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. Perform clustering and UMAP analysis on the combined dataset prior to data integration.
   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. Visualize the UMAP plot according to cluster and batch numbers (Supplementary Figure 26).
   DimPlot(dataset_merged, reduction = "umap.unintegrated", group.by = c("seurat_clusters", "batch"))
6. Show the distribution of cell numbers in each cluster according to batch number.
   table(dataset_merged$batch, dataset_merged$seurat_clusters)
   NOTE: From the UMAP plot and the table, there does not appear to be any significant batch effects for these two datasets. Evidence of batch effects would be manifested as unexpected discrepancies in cluster distribution between the two datasets, which could mean that there are potential technical differences between the datasets that are overriding actual biological similarities.
7. Perform Seurat data integration using the RPCA method. For more information on this and other data integration methods, please read the Seurat vignette linked above.
   dataset_merged <- IntegrateLayers(
object = dataset_merged, method = RPCAIntegration,
orig.reduction = "pca", new.reduction = "integrated.rpca",
verbose = TRUE
)
8. Perform clustering and UMAP analysis on the combined dataset after data integration.
   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. Visualize the UMAP plot according to cluster and batch numbers after integration (Supplementary Figure 27).
   DimPlot(dataset_merged, reduction = "umap.rpca", group.by = c("seurat_clusters","batch"))
10. Show distribution of cell numbers in each cluster according to batch number after integration.
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    From the UMAP plot and the table of the integrated data, there is now excellent overlap between the two batches across different clusters. Interestingly, the integration of the data resulted in the identification of an additional cluster using the same clustering parameters.
11. After dataset integration and prior to downstream analyses, the layers of the merged dataset must be joined.
    dataset_merged <- JoinLayers(dataset_merged)
12. Save the dataset Seurat object as an RDS file into the working directory.
    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 函数识别数据集中可能的双细胞（doublets）事件，并为每个细胞分配一个名为“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 中，使用特征图函数在一系列 UMAP 图上可视化各聚类中排名靠前的标志基因的定位情况，其代表性结果（图 5）为一组 UMAP 图网格，显示各主要细胞类型聚类区域内其对应标志基因的高表达水平。在步骤 3.17 和 3.18 中，使用点图函数可视化细胞中排名靠前的聚类标志基因的相对表达水平，首先按原始 Seurat 聚类编号分组（步骤 3.17），然后按注释的细胞类型标签分组（步骤 3.18）。这些步骤的代表性结果证实，排名靠前的细胞标志基因仅在其对应的 Seurat 聚类中高表达（补充图 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）显示数据集中的成纤维细胞围绕三种颜色编码的细胞亚型聚集。根据 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 生成的单细胞伤口愈合数据集。首先，将整体 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 聚类分布（左图）和批次编号分布（右图），显示两个批次在不同聚类间的重叠程度进一步增加。结果还显示，在数据整合后出现了一个新的聚类，这可能表明在控制了数据批次的技术效应后，识别潜在重要细胞亚型的能力有所提升。

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

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

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

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

材料

本文使用的材料清单

姓名	公司	目录编号	评论
笔记本电脑或台式计算机	N/A	N/A	运行 Windows 或 MacOS
R	N/A	版本 4.4.1	可从 https://cran.rstudio.com/ 免费下载
Rstudio	Posit 软件，PBC	版本 2024.09.0	可从 https://posit.co/download/rstudio-desktop/ 免费下载
Office Excel	微软	任意版本	用于表格数据的分析
互联网浏览器	N/A	N/A	用于访问网站
R 软件包	知识库	版本
开发工具	CRAN	2.4.5
readxl	CRAN	1.4.3
openxlsx	CRAN	4.2.7.1
tidyverse	CRAN	2.0.0
scCustomize	CRAN	2.1.2
BiocManager	Bioconductor	1.30.25
NMF	Bioconductor	0.28
ComplexHeatmap	Bioconductor	2.20.0
BiocNeighbors	Bioconductor	1.22.0
SingleCellExperiment	Bioconductor	1.26.0
circlize	Bioconductor	0.4.16
edgeR	Bioconductor	4.2.1
scDblFinder	Bioconductor	1.18.0
Seurat	CRAN	5.1.0
CellChat	Github	2.1.2