Using R, Seurat, and CellChat to Analyze a Single-Cell Transcriptomics Dataset of Mouse Skin Wound Healing

In our lab, we use emerging tools such as single-cell and spatial transcriptomics with systems biology and bioinformatics approaches to investigate the spatial temporal cellular dynamics of differential healing outcomes. In recent years, we have seen a rapid adoption of single-cell transcriptomics to the study of wound healing in humans and in model organisms, such as mice. The full analysis of single-cell data sets is prohibitive for bench scientists with little to no experience of bioinformatics. This means that all too often, single-cell datasets are underutilized by scientists in the field of wound healing. This is the first comprehensive protocol that assumes no prior experience with bioinformatics, which takes a user all the way from dataset download to the output of relevant analyses in the context of wound healing research. our protocol should serve as a template for wound healing researchers to more fully analyze their own single-cell data sets and be able to extract new insights from publicly available data sets.

[Shalyn] To begin, navigate to the dataset files from the gene expression omnibus repository using the accession number GSE204777. Click on the first dataset titled GSM6190913. Scroll to the bottom of the GSM6190913 page and download the three listed files using either the FTP or HTML links. Using the computer's file explorer, move the downloaded files into a directory named b1, ensuring it is located within the working directory. Retrieve the directory path information for the single-cell sequencing files that were downloaded. Now, load the single-cell sequencing files into the working environment. Then, separate the gene expression and multiplexing HTO data from the working dataset. Create a Seurat object using the gene expression data while filtering out genes detected in fewer than five cells and cells with fewer than 200 genes. For datasets lacking HTO data, create the Seurat object with the same filtering parameters, and switch to the gene expression assay. Calculate the mitochondrial gene percentage in each cell and assign this value as a metadata variable. Visualize the distribution of gene counts, total RNA, and mitochondrial gene percentage across all cells. Remove cells with mitochondrial content exceeding 25% using a threshold and visualize the updated distributions after filtering out these low-quality cells. Detect likely doublets using the SC doublet finder method. Run the SC doublet finder pipeline using the commands and assign the resulting doublet scores as a new metadata variable. Now, visualize the distribution of doublet scores across all cells. Remove all cells with doublet scores above 0.25 and save the cleaned Seurat object as an RDS file in the working directory. Perform data normalization, scaling, and principle component analysis. Visualize the variance contribution across the first 50 principle components. Cluster the cells using the first 13 principle components and a clustering resolution of 0.1. Perform uniform manifold approximation and projection, or UMAP reduction in neighbor analysis, using the first 13 principle components and set the seed number to 123. Now, visualize cell clustering on a UMAP plot, followed by wound time and space annotations on a UMAP plot. Then, generate a table associating cell clusters with wound time and space annotations. Determine major cell type identities after calculating differentially expressed genes between all clusters and assign the resulting DEG lists to a variable before saving them as a delimited text file in the working directory. Now, open the dataset cluster markers file in a spreadsheet application. Use the Text Import Wizard to set the comma as a delimiter and format gene name columns as text to prevent automatic conversion of gene names into dates. In a spreadsheet, rank the average log two FC column from largest to smallest to order the rows by decreasing log twofold change values, followed by the cluster column from smallest to largest. To order the rows by increasing Seurat cluster numbers, filter the average log two FC column to include only values greater than or equal to 2.5, and then filter the PCT 1 column to include values greater than or equal to 0.4. Next, filter the PCT 2 column to include values less than or equal to 0.2. Finally, filter the P-value adjusted column to include values less than or equal to 0.01. Now, open the Enrichr web-based enrichment analysis tool. For each cluster, copy the list of differentially expressed genes into a separate Enrichr window and click analyze. Then, click the cell types tab above the analysis output and focus on the top five enrichments within the three curated cell marker databases. Based on the top enrichments from the Enrichr analysis, assign likely identities to the eight clusters. Combine clusters two and six into a single annotation labeled fibroblasts, and assign these annotations as a new metadata variable named cell types. Then, visualize the annotated cell types on a UMAP plot and display the localization of the bolded top cluster marker genes on a series of UMAP plots. Visualize the top cluster marker differentially expressed genes on a dot plot grouped by original cluster numbers, and then the top marker genes again in a dot plot, this time grouped by annotated cell types. To prepare for time series analysis, remove spatial annotation and simplify the dataset. Reassign the wound time and space metadata into a new variable called DPW for Days Post Wounding. Visualize the new DPW time course groupings on a UMAP plot and generate tables showing the number of cells of each type within each DPW group. Next, convert the cell counts to proportions to assess relative changes in cell type composition during healing and visualize the proportion of each DPW category within each cell type. Finally, visualize the proportion of each cell type within each DPW group and save the final Seurat object containing all annotations and filters as an RDS file into the working directory. All cells in the dataset clustered distinctly into major color-coded cell types on the UMAP plot, confirming successful annotation based on enriched cell type signatures. High expression of top cluster marker genes was localized within their respective cell type clusters on the UMAP plots. Dot plot visualization confirmed the highest expression levels of cluster marker genes were restricted to their annotated major cell types. The stacked bar plot revealed that neutrophils and macrophages were dominant at day one post-wounding, while fibroblasts, epithelial, and endothelial cells became more prevalent at later time points, reflecting the known cellular cascade of wound healing.