Verwendung von R, Seurat und CellChat zur Analyse eines Einzelzell-Transkriptomik-Datensatzes der Wundheilung der Maushaut

Summary

Hier präsentieren wir einen visuellen Schritt-für-Schritt-Workflow zur Analyse eines Einzelzell-Zeitverlauf-Transkriptomik-Datensatzes der Wundheilung der Maushaut mit R. Das Protokoll umfasst eine Standard-Pipeline für das Herunterladen von Datensätzen, Qualitätskontrolle, Visualisierungen und Zelltyp-Annotationen mit Seurat sowie die Analyse von Zell-Zell-Interaktionen mit CellChat.

Abstract

Der Prozess der Wundheilung wird durch komplexe Wechselwirkungen zwischen verschiedenen Zelltypen über Raum und Zeit hinweg reguliert. Durch die Profilierung einzelner Zellen in ihrer komplexen Umgebung ermöglichen Einzelzell-Transkriptomik-Methoden die Untersuchung der zellulären Heterogenität, der Zellkommunikationsnetzwerke und der Zell-Zell-Interaktionen, die am Wundheilungsprozess beteiligt sind. Viele Einzelzellanalysewerkzeuge werden jedoch in einer Computercodierungsumgebung ausgeführt, und ihre breitere Verwendung durch Wundheilungswissenschaftler wird durch den offensichtlichen Mangel an bioinformatischem Fachwissen vereitelt. Daher wird ein Schritt-für-Schritt-Workflow vorgestellt, der zeigt, wie eine grafische Codierungsumgebung namens RStudio verwendet wird, um eine grundlegende Einzelzellanalyse eines temporalen Maus-Exzisionsdatensatzes zur Heilung von Hautwunden durchzuführen. Dieses visuelle und geführte Protokoll ermöglicht es Wissenschaftlern ohne bioinformatischen Hintergrund, einen zuvor veröffentlichten Wundheilungsdatensatz herunterzuladen, kritische Qualitätskontrollschritte durchzuführen, einen standardmäßigen Einzelzellanalyse-Workflow einschließlich Datensatzvisualisierungen und Zelltyp-Annotationen mit Seurat auszuführen, Zell-Subtyp-Analysen durchzuführen, Modul-Scoring-Analysen durchzuführen, Zell-Zell-Interaktionsanalysen mit CellChat durchzuführen und integrative Analysen mehrerer Datensätze mit Seurat durchzuführen. Für jeden Schritt im Protokoll werden narrative Erklärungen bereitgestellt und grafische Ergebnisse aus jeder Codezeile werden präsentiert, um den Benutzer sicher durch den Arbeitsablauf zu führen. Das Ziel dieser visuellen Einführung in eine Einzelzell-Analyse-Pipeline ist es, mehr Wundheilungswissenschaftlern die Möglichkeit zu geben, bioinformatische Werkzeuge direkt in ihren eigenen Laboren zu nutzen, um tiefere Analysen ihrer eigenen Einzelzell-Datensätze sowie eine umfassendere Re-Analyse von zuvor veröffentlichten Einzelzell-Datensätzen zu ermöglichen.

Introduction

Die Wundheilung ist einer der komplexesten Prozesse in der Säugetierbiologie und umfasst ein Spektrum von drei Heilungsphasen: Entzündungsphasen, Proliferativ und Auflösung ^1,2. Diese Heilungsphasen klassifizieren grob die koordinierten Aktionen von Dutzenden von Zelltypen und Hunderten ihrer molekularen Produkte über Raum und Zeit der Wundheilung³. Mehrere Jahrzehnte histologischer und molekularer Studien, die auf Wundgewebeproben über den zeitlichen Verlauf der Heilung basieren, haben die übergreifenden zellulären Muster der Gewebereparatur^{aufgeklärt 3}, insbesondere in reproduzierbaren Mausmodellen der exzisionellen Hautwundheilung ^4,5,6. Erst in den letzten zwei Jahrzehnten ist es möglich geworden, die Komplexität der Wundheilung besser zu verstehen, beginnend mit dem Aufkommen transkriptomischer Hochdurchsatzanalysen von Wunden auf der Skala von Massengeweben ^7,8,9 und Zellen 10,11,12,13,14. In jüngster Zeit haben mehrere Studien Hautwunden auf Einzelzellebene transkriptionell profiliert, neue Wundzell-Subtypen identifiziert und gezeigt, wie sie während der Heilung miteinander interagieren können 15,16,17,18,19,20. Hu et al. verwendeten einen innovativen räumlichen Einzelzell-RNA-Sequenzierungsansatz, um Hautwunden während des gesamten Heilungsverlaufs in mehreren radialen Abständen vom Wundzentrum zu profilieren, was neue interzelluläre und molekulare "Bewegungen" über Raum und Zeit aufzeigte²⁰. Solche Studien enträtseln die Komplexität der Wundheilung in noch nie dagewesener Detailgenauigkeit und zeichnen ein Bild von enormer zellulärer und molekularer Heterogenität.

Wichtige Fortschritte in der bioinformatischen Analysemethode in jüngster Zeit ermöglichen es, die komplexen Multi-Omics-Datensätze, die im Bereich der Wundheilungsforschung generiert werden, biologisch zu verstehen. Einzelzellanalysepakete wie Seurat bieten Werkzeuge für die robuste Analyse und Integration von Datensätzen, einschließlich der Klassifizierung von Zelltypen in komplexen Geweben wie Wunden²¹. Für die nachgelagerte Interpretation von Einzelzelldaten werden Werkzeuge wie CellChat verwendet, um mutmaßliche Zell-Zell-Interaktionsprogramme zu identifizieren, die erklären könnten, wie Zellen sich koordinieren, um Wunden zu reparieren²². Obwohl diese Tools gut dokumentiert und gut zitiert sind, müssen sie in einer Computercodierungsumgebung wie R ausgeführt werden, einer statistischen und grafischen Programmiersprache, die am häufigsten in den Bioinformatikbereichen Genomik und Transkriptomik verwendet wird. Während Biologen und Kliniker auf dem Gebiet der Wundheilung zunehmend Einzelzellansätze verwenden, um die Gewebereparatur zu untersuchen, verfügen nur wenige über die Bioinformatik-Ausbildung, die erforderlich ist, um Tools wie Seurat und CellChat direkt in ihren eigenen Labors zu verwenden. Eine solche Barriere bei der Nutzung dieser Bioinformatik-Tools hindert Wissenschaftler nicht nur daran, ihre eigenen Datensätze ohne die Hilfe von Bioinformatikern tiefer zu analysieren, sondern hindert Wissenschaftler auch daran, die Fülle an Einzelzelldaten, die bereits von anderen Gruppen veröffentlicht wurden, zuverlässig erneut zu analysieren.

Daher wird hier ein Schritt-für-Schritt-Workflow vorgestellt, der es Wissenschaftlern ohne bioinformatischen Hintergrund ermöglicht, einen zuvor veröffentlichten und öffentlich zugänglichen Einzelzell-Wundheilungsdatensatz^{zu analysieren 20}. Das Protokoll verwendet die weit verbreitete und kostenlose grafische R-Codierungsumgebung namens RStudio und zeigt, wie in dieser Umgebung navigiert wird, um vorgeschriebene Codezeilen auszuführen, die eine grundlegende Analyse eines komplexen Einzelzelldatensatzes mit Seurat und CellChat ermöglichen. In diesem Protokoll werden die sieben wichtigsten Methoden vorgestellt, die für die Wundheilungsforschung relevant sind, darunter: 1) Installation der Kodierungsumgebung, 2) Herunterladen des Datensatzes und kritischer Qualitätskontrollschritte, 3) Arbeitsabläufe für Einzelzellanalysen einschließlich Visualisierungen und Zelltyp-Annotationen, 4) Zell-Subtyp-Analysen, 5) Modul-Scoring-Analysen, 6) Zell-Zell-Interaktionsanalysen und 7) integrative Analysen mehrerer Datensätze. Innerhalb jeder Methode wird dem Benutzer der tatsächliche Code zur Verfügung gestellt, der parallel mit dem Protokoll ausgeführt werden kann, und tatsächliche grafische Ergebnisse aus jeder Codezeile werden angezeigt, um den Benutzer durch den Workflow zu führen. Das Hauptziel dieser geführten und visuellen Einführung in RStudio und einen grundlegenden Arbeitsablauf für die Einzelzellanalyse ist es, mehr Wundheilungswissenschaftlern die direkte Nutzung dieser leistungsstarken Werkzeuge zu ermöglichen, um schnellere Fortschritte auf dem Gebiet der Forschung zu ermöglichen.

Protocol

ANMERKUNG: In den folgenden Workflows, die sieben Bioinformatik-Methoden beschreiben, werden alle Schritte der Protokolle von ihren jeweiligen Codeblöcken begleitet, die direkt auf der benutzereigenen RStudio-Oberfläche in der Reihenfolge ausgeführt werden sollten, in der sie aufgelistet sind. Um dieses Protokoll so benutzerfreundlich wie möglich zu gestalten, ist eine R-Skriptdatei enthalten (Supplementary File 1: JoVE_Rscript.R), die direkt in die RStudio-Sitzung des Benutzers geladen werden kann, so dass jede Codezeile einfach ausgeführt werden kann. Dadurch wird vermieden, dass der Benutzer den Code aus dem Protokolldokument eingeben oder kopieren und einfügen muss, was zu Fehlern führen könnte. Alle Protokollanweisungen sind auch in der R-Skriptdatei in Form von Kommentaren enthalten, die durch ein Hashtag-Symbol "#" am Anfang jeder Kommentarzeile gekennzeichnet sind.

1. Installieren von R, RStudio und den erforderlichen R-Paketen für den Workflow für die Einzelzellanalyse

1. Laden Sie R (Version 4.4.1) herunter, und installieren Sie es auf dem Computer. Verwenden Sie den Link, der dem Betriebssystem des Computers entspricht.
   1. Wenn Sie einen Computer mit Microsoft Windows verwenden, verwenden Sie diesen Link: https://cran.rstudio.com/bin/windows/base/
   2. Wenn Sie einen Computer mit MacOS verwenden, verwenden Sie diesen Link: https://cran.rstudio.com/bin/macosx/
2. Installieren Sie die neueste Version von RStudio auf dem Computer. Klicken Sie auf den folgenden Link und folgen Sie den Anweisungen:
   https://www.rstudio.com/products/rstudio/download/#download
3. Installieren Sie die Rtools (Version 4.4), die es R ermöglichen, bestimmte Pakete zu kompilieren. Klicken Sie auf den folgenden Link und folgen Sie den Anweisungen:
   1. Wenn Sie Windows verwenden, verwenden Sie diesen Link: https://cran.rstudio.com/bin/windows/Rtools/
   2. Wenn Sie MacOS verwenden, verwenden Sie diesen Link:
      https://mac.r-project.org/tools/
4. Legen Sie das lokale Arbeitsverzeichnis fest. Dies ist der Ordner auf dem Computer, in den alle Dateien geladen und gespeichert werden. Legen Sie das Arbeitsverzeichnis fest, indem Sie in der RStudio-Menüleiste Sitzung auswählen und auf Arbeitsverzeichnis festlegen > Verzeichnis auswählen klicken und den gewünschten Ordner auswählen.
   1. Wenn Sie einen Windows-Computer verwenden, verwenden Sie den folgenden Befehl, um das Arbeitsverzeichnis festzulegen. Ändern Sie [Directory] in der folgenden Codezeile in die tatsächliche Verzeichnisstruktur. Beachten Sie, dass das Verzeichnistrennzeichen in R das Zeichen "/" ist.
      setwd("C:/[Directory]")
   2. Wenn Sie einen MacOS-Computer verwenden, wird mit dem folgenden Befehl auch das Arbeitsverzeichnis festgelegt. Ändern Sie [Directory] in der folgenden Codezeile in die tatsächliche Verzeichnisstruktur. Beachten Sie, dass das Verzeichnistrennzeichen in R das Zeichen "/" ist.
      setwd("~/[Directory]")
   3. Überprüfen Sie zu einem beliebigen Zeitpunkt während einer R-Sitzung das Arbeitsverzeichnis mit der folgenden Codezeile:
      getwd()
   4. Untersuchen Sie in RStudio visuell die Struktur des Arbeitsverzeichnisses, einschließlich aller darin enthaltenen Dateien und Ordner, im rechten Fenster auf der Registerkarte Dateien . Um den Datei-Explorer von RStudio zum Arbeitsverzeichnis zu navigieren, klicken Sie auf das Zahnradsymbol -> Zum Arbeitsverzeichnis wechseln.
5. Installieren Sie die folgenden Pakete aus dem R-Paketrepository CRAN, bei denen es sich um notwendige Abhängigkeiten für das Protokoll handelt. Führen Sie die folgenden Befehle aus, um diese Pakete zu installieren.
   install.packages("devtools")
install.packages("readxl")
install.packages("openxlsx")
install.packages("tidyverse")
install.packages("scCustomize")
   ANMERKUNG: Während der Installation von R-Paketen ist es normal, dass verschiedene Fenster angezeigt und ausgeblendet werden. Wenn ein Fenster angezeigt wird, in dem Sie aufgefordert werden, ein Paket zu kompilieren, klicken Sie auf JA. Wenn ein Fenster angezeigt wird, in dem Sie aufgefordert werden, R vor der Installation des Pakets neu zu starten, klicken Sie auf NEIN.
6. Installieren Sie die folgenden Pakete aus dem kuratierten R-Paketrepository Bioconductor:
   (https://bioconductor.org/), die notwendige Abhängigkeiten für das Protokoll sind. Führen Sie die folgenden Befehle aus, um diese Pakete zu installieren.
   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. Installieren Sie die folgenden Pakete, die für den in diesem Manuskript beschriebenen Workflow erforderlich sind.
   install.packages("Seurat")
devtools::install_github("jinworks/CellChat")
8. Laden Sie jedes Paket, um zu bestätigen, dass die Installationen erfolgreich waren. Falls eines der Pakete den Fehler "Paket nicht gefunden" auslöst, installieren Sie es mit dem entsprechenden Code oben neu.
   library(readxl)
library(openxlsx)
library(tidyverse)
library(scCustomize)
library(edgeR)
library(scDblFinder)
library(Seurat)
library(CellChat)

2. Laden eines einzelligen Wundheilungsdatensatzes und Durchführen von Qualitätskontrollschritten

ANMERKUNG: Für diesen bioinformatischen Workflow wird eine Re-Analyse eines zuvor veröffentlichten räumlich-temporalen Experiments zur Wundheilung von Einzelzellen der Haut durchgeführt²⁰. Die Datensatzdateien werden im kuratierten NCBI Gene Expression Omnibus (GEO) Repository (https://www.ncbi.nlm.nih.gov/geo/) gespeichert.

1. Navigieren Sie zu den Datensatzdateien aus GEO mit der Zugangsnummer GSE204777. Verwenden Sie den folgenden Link und überprüfen Sie das Versuchsdesign der Studie: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204777
2. Auf der GEO-Repository-Seite gibt es fünf einzelne Datenchargen, die aus fünf Sequenzierungsbahnen stammen. Klicken Sie auf den ersten Datensatz mit dem Titel GSM6190913. Im Folgenden finden Sie einen direkten Link zur Stichprobe: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM6190913
3. Scrollen Sie zum Ende der Seite und laden Sie die folgenden drei Dateien herunter, entweder über die FTP- oder HTML-Links . Verschieben Sie diese drei Dateien im Datei-Explorer des Computers in ein Verzeichnis mit dem Namen b1. Stellen Sie sicher, dass sich der Ordner b1 in dem Arbeitsverzeichnis befindet, das in Schritt 1.4 festgelegt wurde.
   Dateiname: GSM6190913_b1_barcodes.tsv.gz / Dateigröße: 18.5 Mb
   Dateiname: GSM6190913_b1_features.tsv.gz / Dateigröße: 254.1 Kb
   Dateiname: GSM6190913_b1_matrix.mtx.gz / Dateigröße: 151.2 Mb
4. Rufen Sie die Verzeichnisinformationen für die Einzelzell-Sequenzierungsdateien ab, die Sie in Schritt 2.3 heruntergeladen haben.
   b1_data_dir <- file.path(getwd(), "b1")
5. Laden Sie die Einzelzell-Sequenzierungsdateien. Der Parameter gene.column gibt die verwendete Gen-/Merkmalsbenennung an. Verwenden Sie in diesem Fall gene.column = 2 für Gensymbole (gene.column = 1 steht für Ensembl-Gennamen).
   b1 <- Read10X_GEO(data_dir = b1_data_dir, gene.column = 2)
   ANMERKUNG: Die meisten Einzelzellen-Datasets verfügen nicht über zusätzliche Multiplexdaten, sodass die 10x-Dateien, die mit diesem Schritt generiert werden, nicht über mehrere Layer verfügen. Für den aktuellen Multiplex-Datensatz fahren Sie mit Schritt 2.6 fort. Für Datensätze ohne Multiplexing von Daten fahren Sie mit Schritt 2.7 fort.
6. Demultiplexen des Einzelzell-Datasets mit räumlich-zeitlichen Barcodes.
   1. Trennen Sie für den Arbeitsdatensatz die Genexpressions- und HTO-Daten (Multiplexing).
      dataset_rna_counts <- b1$GSM6190913_b1_$`Gene Expression`
​dataset_barcodes <- b1$GSM6190913_b1_$`Antibody Capture`
   2. Erstellen Sie ein Seurat-Objekt unter Verwendung der Genexpressionsdaten, während Sie sofort Gene herausfiltern, die in weniger als 5 Zellen exprimiert werden, und Zellen mit weniger als 200 nachgewiesenen Genen.
      dataset <- CreateSeuratObject(counts = dataset_rna_counts, min.cells = 5, min.features = 200)
   3. Erstellen Sie einen Genexpressionsdatensatz als Schicht und generieren Sie eine Liste von Zellen und Barcodes, die beiden Assays gemeinsam sind.
      dataset_rna_counts2 <- LayerData(dataset, data = "RNA")
​dataset_joint.barcodes <- intersect(colnames(dataset_rna_counts2), colnames(dataset_barcodes))
   4. Teilmenge der Genexpression und HTO-Anzahl durch gemeinsame Zell-Barcodes.
      dataset_rna_counts2 <- dataset_rna_counts2[, dataset_joint.barcodes]
​dataset_barcodes2 <- as.matrix(dataset_barcodes[, dataset_joint.barcodes])
   5. Vergewissern Sie sich, dass die HTO über die erwarteten Barcodenamen verfügen.
      ​rownames(dataset_barcodes2)
   6. Erstellen Sie einen neuen Assay, um Barcode-Informationen zu speichern, und fügen Sie diesen Assay dem zuvor erstellten Seurat-Objekt hinzu.
      dataset_barcode_assay <- CreateAssayObject(counts = dataset_barcodes2)
​dataset[["barcodes"]] <- dataset_barcode_assay
   7. Überprüfen Sie, ob das Objekt jetzt mehrere Assays enthält.
      DefaultAssay(dataset)
   8. Normalisieren Sie Barcode-Daten und führen Sie ein Demultiplexing über die HTODemux-Funktion durch. Diese Methode wird in der folgenden Seurat Vignette ausführlich beschrieben: https://satijalab.org/seurat/articles/hashing_vignette
      ​dataset <- HTODemux(dataset, assay = "barcodes", positive.quantile = 0.99)
   9. Gruppieren Sie Zellen basierend auf globalen Klassifizierungsergebnissen und entfernen Sie Zellen ohne Barcode-Klassifizierung.
      Idents(dataset) <- "barcodes_classification.global"
​dataset <- subset(dataset, idents = "Negative", invert = TRUE)
   10. Gruppieren Sie Zellen basierend auf dem maximalen HTO-Signal.
       Idents(dataset) <- "barcodes_maxID"
   11. Visualisieren Sie die Verteilung der detektierten Gene in Zellen anhand ihrer gemultiplexten Barcodes (Ergänzende Abbildung 1).
       VlnPlot(dataset, features = "nFeature_RNA", pt.size = 0.1, log = TRUE)
   12. Benennen Sie die Barcodes in ihre tatsächliche Wundzeit (Tage nach der Verwundung) und ihre Raumzuordnungen (2-8 mm) (aus dem Originalmanuskript) um und weisen Sie sie einer neuen Metadatenvariablen namens time_space zu.
       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. Für Datensätze ohne Multiplexing: Erstellen Sie ein Seurat-Objekt, während Sie sofort Gene herausfiltern, die in weniger als 5 Zellen exprimiert werden, und Zellen mit weniger als 200 nachgewiesenen Genen.
   dataset <- CreateSeuratObject(counts = b1, min.cells = 5, min.features = 200)
8. Wechseln Sie zur Analyse des Genexpressionsassays des Datensatzes.
   DefaultAssay(dataset) <- "RNA"
Idents(dataset) <- "data"
9. Als wichtiger Schritt der Qualitätskontrolle berechnen Sie den Prozentsatz der mitochondrialen Gene in jeder Zelle und weisen Sie ihn als Metadatenvariable zu. Diese Methode wird in der folgenden Seurat Vignette ausführlich beschrieben: https://satijalab.org/seurat/articles/pbmc3k_tutorial#qc-and-selecting-cells-for-further-analysis
   dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^mt-")
10. Visualisieren Sie die Verteilung der detektierten Gene, die Anzahl der RNA und den mitochondrialen Anteil in allen Zellen (Ergänzende Abbildung 2).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
11. Es gibt eine Reihe von Zellen mit einem hohen mitochondrialen Gehalt, der mit niedrigen RNA-Zahlen korreliert; Dabei handelt es sich um tote oder sterbende Zellen. Entfernen Sie diese Zellen mit geringer Qualität aus dem Dataset unter Verwendung eines fairen Grenzwerts. Verwenden Sie in diesem Fall die Werte des ursprünglichen Datensatzes, der in der zuvor veröffentlichten Studie²⁰ beschrieben wurde, in der Zellen mit mehr als 25 % mitochondrialen Genen entfernt werden.
    dataset <- subset(dataset, subset = nFeature_RNA > 200 & percent.mt < 25)
12. Visualisieren Sie die Verteilung der detektierten Gene, die Anzahl der RNA und den mitochondrialen Prozentsatz in allen Zellen nach der Entnahme minderwertiger Zellen (Ergänzende Abbildung 3).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
13. Als weiterer wichtiger Schritt der Qualitätskontrolle ist die Erkennung wahrscheinlicher Dubletten im Datensatz. Dabei handelt es sich um Zellen, die bei der Tröpfchensequenzierung zusammengefügt wurden und somit zu Genexpressionen führen, die nicht auf Einzelzellebene liegen. Um diesem Problem entgegenzuwirken, wurden mehrere Tools entwickelt. Es ist wichtig zu beachten, dass jedes Tool bestimmte Annahmen über Einzelzellendaten trifft, daher ist es wichtig, dass der Benutzer alle relevanten Dokumentationen liest, bevor er ein Tool verwendet. Implementieren Sie in diesem Workflow eine Methode namens scDblFinder²³. Bitte beachten Sie, dass diese Methode einen Algorithmus verwendet, der eine feste erwartete Dublettenrate vorschreibt. Verwenden Sie die folgenden Befehle, um die scDblFinder-Pipeline auszuführen.
    DefaultAssay(dataset) <- "RNA"
sce <- scDblFinder(GetAssayData(dataset, assay="RNA", slot="counts"), samples=Idents(dataset))
    ANMERKUNG: Multiplex-Einzelzelldatensätze wie dieser können auch auf Dubletten untersucht werden, indem Zellen entfernt werden, die mehrere Barcodes ausdrücken. In diesem Workflow wurde diese Methode nicht veranschaulicht, da die meisten einzelligen Datasets nicht über diese einzigartige Funktion verfügen. Stattdessen wird eine generalisierbarere Pipeline für das Screening auf Dubletten mit der scDblFinder-Methode gezeigt.
14. Weisen Sie den Doublet-Score einer neuen Metadatenvariablen zu.
    dataset$scDblFinder.score <- sce$scDblFinder.score
15. Visualisieren Sie die Verteilung des Doublet-Scores in allen Zellen (Ergänzende Abbildung 4).
    VlnPlot(dataset, features = "scDblFinder.score", raster=FALSE, pt.size=0.5)
16. Entfernen Sie die Zellen über dem Schwellenwert von 0,25 mit doppelter Punktzahl. Dieser Schwellenwert wurde auf der Grundlage des oben generierten Violin-Diagramms ausgewählt, das zeigte, dass die meisten Zellen im Datensatz entweder sehr hohen oder sehr niedrigen Doppelwerten zugeordnet werden konnten, und 0,25 ist ein vernünftiger Grenzwert für diesen Datensatz, der die überwiegende Mehrheit der wahrscheinlichen Dubletten entfernen würde, ohne viele unwahrscheinliche Dubletten zu entfernen.
    dataset <- subset(dataset, scDblFinder.score < 0.25)
17. Speichern Sie das Seurat-Objekt des Datensatzes als RDS-Datei in das Arbeitsverzeichnis.
    saveRDS(dataset, "dataset_post_Method2.rds")

3. Analyse eines Einzelzell-Wundheilungsdatensatzes mit Seurat

HINWEIS: (Optionaler Schritt) Wenn Sie den Workflow hier starten, laden Sie die gespeicherte RDS-Datei als Seurat Objekt.

dataset <- readRDS("dataset_post_Method2.rds")

1. Führen Sie den standardmäßigen Seurat-Workflow für die Normalisierung, Skalierung und Principal Component Analysis (PCA) von Einzelzelldatensätzen durch. Dieser Standard-Workflow wird in den folgenden Seurat Vignetten beschrieben:
   Seurat - Geführtes Clustering-Tutorial: https://satijalab.org/seurat/articles/pbmc3k_tutorial
   Liste der Seurat-Befehle: https://satijalab.org/seurat/articles/essential_commands
   DefaultAssay(Datensatz) <- "RNA"
   dataset <- NormalizeData(object = dataset)
dataset <- FindVariableFeatures(object = dataset)
dataset <- ScaleData(object = dataset)
dataset <- RunPCA(object = dataset)
2. Visualisieren Sie das Ausmaß der Variation des Datensatzes in Bezug auf die ersten 50 PCA-Dimensionen (Ergänzende Abbildung 5).
   ElbowPlot(dataset, reduction = "pca", ndims = 50)
   Ein Großteil der großen Variation findet innerhalb der ersten 13 Dimensionen statt.
3. Führen Sie das Zell-Clustering des Datensatzes mit einem festgelegten PCA-Dimensionsbereich von 1 bis 13 und einer festgelegten Auflösung von 0,1 durch.
   dataset <- FindNeighbors(dataset, verbose = TRUE, dims = 1:13)
dataset <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
   ANMERKUNG: Dieser Workflow konzentriert sich auf große Unterschiede bei den Zelltypen. Daher wird ein relativ konservativer Parameter für den PCA-Dimensionsbereich verwendet, wobei die ersten 13 Dimensionen die überwiegende Mehrheit der Variation innerhalb des Datensatzes darstellen. Für die Unterscheidung von Zellen in kleinere und seltenere Subtypen kann der Benutzer eine höhere Anzahl von Dimensionen für die nachgelagerte Analyse verwenden, da diese seltenen Zell-Subtypen wahrscheinlich ein geringeres Maß an Variation des Datensatzes aufweisen. Der Auflösungsparameter reicht von 0 bis 1 und bestimmt die Größe der kategorischen Trennung, die dem Dataset auferlegt wird. Die Einstellung dieses Parameters hängt von der Forschungsfrage des Benutzers ab. Für die Clusterung von Zellen in zahlreiche kleine und seltene Subtypen verwenden Sie höhere Auflösungen für die nachgelagerte Analyse. Da dieser Workflow darauf abzielt, allgemeinere Unterschiede zwischen den Hauptzelltypen zu untersuchen, wird ein relativ kleiner Auflösungswert von 0,1 verwendet, der die Zellen in weniger, aber größere Gruppen gruppieren soll. Bei Verwendung der in Schritt 3.3 genannten Einstellungen werden 8 eindeutige Zellcluster unterschieden. Diese Cluster werden automatisch einer Metadatenvariablen namens "seurat_clusters" zugewiesen.
4. Führen Sie eine UMAP-Dimensionsreduktion und eine Analyse der Suche nach Nachbarn anhand der ersten 13 PCA-Dimensionen durch. Fügen Sie die Startnummer 123 hinzu, um die Reproduzierbarkeit der resultierenden Datenprojektion sicherzustellen.
   dataset <- RunUMAP(dataset, verbose = TRUE, dims = 1:13, seed.use = 123)
   ANMERKUNG: Der UMAP-Algorithmus ist stochastisch und führt Zufälligkeit in die Dimensionsreduktion ein (siehe "Stabilität und Reproduzierbarkeit" in https://cran.r-project.org/web/packages/umap/vignettes/umap.html). Während die Verwendung eines konsistenten Ausgangswerts dem Algorithmus ein "minimales Maß an Reproduzierbarkeit" verleiht, kann sich das resultierende Diagramm dennoch geringfügig von dem unterscheiden, was in den repräsentativen Abbildungen und nachgelagerten Ergebnissen gezeigt wird. Tests haben ergeben, dass sich die Ergebnisse insbesondere zwischen Computern mit Windows und MacOS unterscheiden, wahrscheinlich aufgrund unterschiedlicher Implementierungen von Zufälligkeit in diesen Betriebssystemen.
5. Visualisieren Sie die Cluster-Bildung der Zellen in einem UMAP-Diagramm (Abbildung 1).
   DimPlot(dataset, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
   ANMERKUNG: Die Zufälligkeit des UMAP-Algorithmus kann leicht unterschiedliche Diagramme erzeugen, wie in der Abbildung gezeigt, die alternative UMAP-Diagramme zeigt, die mit demselben Code wie oben auf Computern mit Windows und MacOS erstellt wurden. Beachten Sie die geringfügigen Unterschiede in den Formen der Cluster. Es ist daher zwingend erforderlich, dass der Benutzer alle Daten und Diagramme speichert und mit einem Zeitstempel versieht, sobald sie erstellt werden, und dass alle nachgelagerten Analysen von Clustern sorgfältig und unter Berücksichtigung des biologischen Verständnisses durchgeführt werden, wie es im Folgenden für die Zelltyp-Annotation beschrieben wird.
6. Da die ursprünglichen Experimentbeschriftungen, die sich darauf beziehen, woher und wann die Zellen während der Wundheilung stammten, enthalten sind, visualisieren Sie die Wundzeit-/Raum-Annotation der Zellen in einem UMAP-Diagramm (Abbildung 2).
   DimPlot(dataset, group.by = "time_space", raster = FALSE, label = FALSE)
7. Generieren Sie eine Tabelle der Zellcluster-Assoziation mit der gewickelten Zeit-/Raumanmerkung.
   table(dataset$time_space, dataset$seurat_clusters)
8. Bestimmen Sie die Identitäten der wichtigsten Zellentypen im Dataset. Berechnen Sie dazu die differentiell exprimierten Gene (DEG) zwischen allen Clustern. Rufen Sie DEG-Listen für Cluster ab, weisen Sie sie einer Variablen zu, und speichern Sie die Ausgabe als Textdatei mit Trennzeichen im Arbeitsverzeichnis.
   ANMERKUNG: Dieser Schritt ist CPU-intensiv und kann je nach Hardware des Benutzers viel Zeit in Anspruch nehmen.
   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. Laden Sie die enthaltene dataset_cluster_markers.txt Datei herunter und öffnen Sie sie in einer Tabellenkalkulation (z. B. Excel), indem Sie den Inhalt der Textdatei kopieren und den Textimport-Assistenten verwenden, um das Kommatrennzeichen und die Identität der Gennamenspalten als Text anzugeben. Die Angabe, dass es sich bei den Gennamen um 'Text' handelt, ist wichtig, da Excel sonst bestimmte Gennamen automatisch in Daten umwandelt, z.B. vom 7. September bis zum 7. September.
10. Filtern Sie die Ergebnisse in einer Tabelle nach den folgenden empfohlenen Parametern:
    1. Ordnen Sie die avg_log2FC Spalte von der größten zur kleinsten, um alle Zeilen nach abnehmenden 2-fachen Änderungen (log2FC) anzuordnen.
    2. Ordnen Sie die Clusterspalte von der kleinsten zur größten, um alle Zeilen entsprechend der steigenden Anzahl von Seurat-Clustern anzuordnen.
    3. Filtern Sie die avg_log2FC Spalte nach Zahlen größer oder gleich 2,5, um nur die am stärksten differentiell exprimierten Gene (DEGs) im angegebenen Cluster im Vergleich zu anderen Clustern anzuzeigen.
    4. Filtern Sie die Spalte pct.1 nach Zahlen, die größer oder gleich 0,4 sind. Diese Spalte bezieht sich auf den Prozentsatz der Zellen im angegebenen Cluster, die das angegebene Gen exprimieren (Cluster %), und die Festlegung des Schwellenwerts auf 0,4 bedeutet, dass nur Gene angezeigt werden, die in mindestens 40 % der Zellen im angegebenen Cluster exprimiert werden.
    5. Filtern Sie die Spalte pct.2 nach Zahlen kleiner oder gleich 0,2. Diese Spalte bezieht sich auf den Prozentsatz der Zellen, die NICHT im angegebenen Cluster das angegebene Gen exprimieren (Nicht-Cluster-%), und die Festlegung des Schwellenwerts auf 0,2 bedeutet, dass nur Gene angezeigt werden, die in höchstens 20 % der Zellen exprimiert werden, die NICHT im angegebenen Cluster vorhanden sind.
    6. Filtern Sie die p_val_adj Spalte nach Zahlen kleiner oder gleich 0,01. Diese Spalte bezieht sich auf den angepassten P-Wert oder die False Discovery Rate (FDR), die die statistische Stärke der angegebenen DEG angibt, und die Festlegung des Schwellenwerts auf 0,01 bedeutet, dass nur Gene mit einer FDR < 0,01 angezeigt werden.
       HINWEIS: Die ergänzende Tabelle 1 (JoVE_DEGs_cellMarkers.xlsx) enthält die vollständige Ausgabe der in Schritt 3.10 des Protokolls verwendeten differentiell exprimierten Gene in der Rangfolge. Ergänzende Tabelle 2 zeigt die Top 5 Gene für jeden Cluster in dieser Analyse, wobei die fett gedruckten Gene für nachfolgende Visualisierungen verwendet werden.
11. Verwenden Sie für eine unverzerrte Annotation von Clustern nach Zelltypen das webbasierte Anreicherungsanalysewerkzeug EnrichR.
    Verwenden Sie den Link: https://maayanlab.cloud/Enrichr/
12. Kopieren Sie die Listen der DEGs für jeden Cluster in ein separates EnrichR-Fenster und klicken Sie dann auf Analysieren. Das EnrichR-Tool lässt die Genliste durch Hunderte von kuratierten Datenbanken laufen und ordnet jeden angereicherten Begriff in jeder Kategorie ein.
13. Klicken Sie für die Annotation von Zelltypen oben auf die Registerkarte Zelltypen und konzentrieren Sie sich auf die Top 5 Anreicherungen in den drei von Zellmarkern kuratierten Datenbanken auf der linken Seite (Abbildung 3):
    CellMarker 2024 (http://bio-bigdata.hrbmu.edu.cn/CellMarker/)
    Tabula sapiens (https://tabula-sapiens-portal.ds.czbiohub.org/)
    PanglaoDB Erweitert (https://panglaodb.se/)
14. Bestätigen Sie auf der Grundlage der Anreicherungen der DEGs in diesen Datenbanken die wahrscheinliche Identität der 8 Cluster. Beachten Sie, dass es zwei Cluster (2, 6) gibt, die sich als Fibroblasten anreichern; Kombinieren Sie daher diese Cluster in Anmerkungen mit einem einzelnen Zellentyp. Weisen Sie die Zellentypidentitäten als Beschriftungen einer neuen Metadatenvariablen mit dem Namen cell_types zu.
    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. Visualisieren Sie die umbenannten Zellcluster als Anmerkungen in einem UMAP-Diagramm (Abbildung 4).
    DimPlot(dataset, group.by = "cell_types", raster = FALSE, label = FALSE)
16. Visualisieren Sie die Lokalisierung der oberen (fett gedruckten) Cluster-Markergene aus Tabelle 1 in einer Reihe von UMAP-Diagrammen (Abbildung 5).
    DefaultAssay(dataset) <- "RNA"
FeaturePlot(dataset, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4)
17. Visualisieren Sie die DEGs der Top-Cluster-Marker in einem Punktdiagramm, gruppiert nach den ursprünglichen Cluster-Nummern (Ergänzende Abbildung 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. Visualisieren Sie die DEGs der Top-Cluster-Marker in einem Punktdiagramm, gruppiert nach den mit Anmerkungen versehenen Zelltypen (Abbildung 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. Um Zeitreihenanalysen durchzuführen, vereinfachen Sie zunächst das Dataset, um die räumliche Komponente zu entfernen. Gruppieren Sie für Zeitverlaufsanalysen die Wundzeit-/Raum-Annotationen in Gesamttage nach der Verwundung (DPW) mit einer neuen Metadatenvariablen namens "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. Visualisieren Sie die neuen gewundenen Zeitverlaufsgruppierungen in einem UMAP-Diagramm (Ergänzende Abbildung 7).
    DimPlot(dataset, group.by = "DPW", raster = FALSE, label = FALSE)
21. Generieren Sie Tabellen, die anzeigen, wie viele Zellen jedes Typs in jedem DPW vorkommen.
    ​table(dataset$DPW, dataset$cell_types)
    1. Optionaler Schritt: Um auch DEG-Listen für verwundete Zeitverlaufsgruppen abzurufen, weisen Sie sie einer Variablen zu, und speichern Sie die Ausgabe als Textdatei mit Trennzeichen.
       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. Rechnen Sie die Zellzahlen in Anteile pro Kategorie um, um die relativen Veränderungen in der Zelltypzusammensetzung im Laufe des Heilungsverlaufs besser zu verstehen. Visualisieren Sie den Anteil von DPW in jedem Zelltyp (Ergänzende Abbildung 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. Visualisieren Sie den Anteil der Zelltypen in jedem DPW ( Abbildung 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. Speichern Sie das Seurat-Objekt des Datensatzes als RDS-Datei in das Arbeitsverzeichnis.
    saveRDS(dataset, "dataset_post_Method3.rds")

4. Analyse von Zellsubtypen mit Seurat

ANMERKUNG: Die Leistungsfähigkeit der Einzelzellanalyse ermöglicht die Entdeckung und Analyse seltener Subtypen innerhalb der wichtigsten Zelltypen, die oben analysiert wurden. Dieses Beispiel konzentriert sich auf Fibroblasten, die sich zunächst in zwei Seurat-Clustern gruppierten, bevor sie zu einer einzigen Kategorie zusammengefasst wurden. Dieser Teil des Protokolls konzentriert sich speziell auf Fibroblasten, schließt alle anderen Zelltypen aus und untersucht ihre Identitäten und zeitlichen Eigenschaften während der Wundheilung. Laden Sie als optionalen Schritt bei Bedarf die gespeicherte RDS-Datei als Seurat-Objekt.

dataset <- readRDS("dataset_post_Method3.rds")

1. Unterteilen Sie den ursprünglichen Datensatz entsprechend der Identität der Fibroblastenzelle.
   Idents (Datensatz) <- "cell_types"
   dataset_fibroblast <- subset(dataset, idents = "Fibroblast")
2. Führen Sie PCA für diesen kleineren Datensatz durch und visualisieren Sie das Ausmaß der Datensatzvariation in Bezug auf die PCA-Dimensionen (ergänzende Abbildung 9).
   dataset_fibroblast <- RunPCA(dataset_fibroblast, verbose = TRUE)
ElbowPlot(dataset_fibroblast, reduction = "pca", ndims = 50)
   ANMERKUNG: Ein Großteil der großen Variation findet innerhalb der ersten 9 Dimensionen statt.
3. Führen Sie das Zell-Clustering des Datensatzes mit einem festgelegten PCA-Dimensionsbereich von 1 bis 9 und einer festgelegten Auflösung von 0,1 durch.
   dataset_fibroblast <- FindNeighbors(dataset_fibroblast, verbose = TRUE, dims = 1:9)
dataset_fibroblast <- FindClusters(dataset_fibroblast, verbose = TRUE, resolution = 0.1)
   ANMERKUNG: Mit diesen Einstellungen unterscheidet der Algorithmus 3 einzigartige Fibroblastencluster. Diese Cluster werden automatisch einer Metadatenvariablen mit dem Namen seurat_clusters zugewiesen.
4. Führen Sie die UMAP-Bemaßungsreduzierung und die Suche nach Nachbarn anhand der ersten 9 PCA-Dimensionen durch. Fügen Sie die Startnummer 123 hinzu, um die Reproduzierbarkeit der resultierenden Datenprojektion sicherzustellen.
   dataset_fibroblast <- RunUMAP(dataset_fibroblast, verbose = TRUE, dims = 1:9, seed.use = 123)
5. Visualisieren Sie die Cluster-Bildung der Zellen in einem UMAP-Diagramm (Abbildung 8).
   DimPlot(dataset_fibroblast, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
6. Visualisieren Sie die Wund-Zeitverlauf-Annotation der Zellen in einem UMAP-Diagramm (Ergänzende Abbildung 10).
   DimPlot(dataset_fibroblast, group.by = "DPW", raster = FALSE, label = FALSE)
7. Holen Sie sich DEG-Listen für die drei Fibroblasten-Subtypen und speichern Sie sie in einer Textdatei im Arbeitsverzeichnis.
   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. Befolgen Sie ähnliche Schritte wie oben (Schritte 3.9-3.10) in einer Tabelle, um die DEGs nach den Top-Fibroblasten-Subtyp-Markern zu filtern.
9. Definieren Sie eine benutzerdefinierte Genliste als Variable, indem Sie die Top 5 der Fibroblasten-Subtyp-Marker für jeden der drei Cluster aus der DEG-Textdatei kopieren.
   FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
10. Visualisieren Sie die Gene in der Liste im Nur-Fibroblasten-Datensatz, indem Sie die Variable im Merkmalsparameter des Punktdiagramms aufrufen (Abbildung 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. Visualisieren Sie die Gene in der Liste im ursprünglichen Einzelzelldatensatz, indem Sie die Variable im Merkmalsparameter des Punktdiagramms aufrufen (Ergänzende Abbildung 11).
    DotPlot(dataset, group.by="cell_types", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
12. Visualisieren Sie den Anteil der Fibroblasten-Subtypen in jeder DPW-Kategorie (Ergänzende Abbildung 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. Visualisieren Sie den Anteil der DPW-Fibroblastenzellen in jeder Fibroblasten-Subtypkategorie (Ergänzende Abbildung 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. Speichern Sie das Seurat-Objekt des Datensatzes als RDS-Datei in das Arbeitsverzeichnis.
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

5. Beispiel für eine Nachverfolgungsanalyse über das Modul-Scoring

ANMERKUNG: Eine nützliche Methode zur Analyse von Einzelzelldatensätzen ist das sogenannte Modul-Scoring. In diesem Workflow kann man eine Genliste nach Vorwissen definieren und dann Modul-Scores berechnen, die potenzielle Anreicherungen der Genliste innerhalb jeder Zelle identifizieren können. Diese Bewertungen können über Zellanmerkungen hinweg gemittelt werden, um potenzielle Muster von Anreicherungen aufzudecken.

Verwenden Sie hier Genlisten aus einer zuvor veröffentlichten Studie², in der Wundheilungsphasen-spezifische Gene anhand von RNA-sequenzierten Massenproben aus dem gesamten Heilungskontinuum identifiziert wurden. Die Genlisten wurden in einer tabulatorgetrennten Textdatei (Ergänzungsdatei 2: JoVE_PhaseSpecificGenes.txt) gespeichert, die nun in das Arbeitsverzeichnis heruntergeladen werden kann und zur Erstellung von Genlisten verwendet werden kann, die die drei wichtigsten Heilungsphasen identifizieren.

1. Laden Sie Genlisten in eine Variable, indem Sie die TEXT-Datei lesen.
   PhaseSpecificGenes <- read_delim("JoVE_PhaseSpecificGenes.txt", delim="\t", col_names = T)
2. Trennen Sie die Spalten in einzelne Genlistenvariablen, und ändern Sie die Gene in Mausnamen, bei denen der erste Buchstabe groß geschrieben wird.
   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)
   HINWEIS: (Optional) Laden Sie bei Bedarf die gespeicherte RDS-Datei als Seurat Objekt.
   dataset <- readRDS("dataset_post_Method3.rds")
3. Verwenden Sie die Genlisten als Module, um jede Zelle im Datensatz nach den drei Phasen der Heilung zu bewerten.
   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. Visualisieren Sie die aggregierten Modulbewertungen pro Zellenkategorie, einschließlich DPW und Hauptzellentypen (Abbildung 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. Beispiel für eine Nachanalyse über CellChat

ANMERKUNG: Eine weitere nützliche und viel zitierte Methode zur Analyse von Einzelzelldatensätzen ist die Ableitung von Zell-Zell-Interaktionen. Verwenden Sie in diesem Arbeitsablauf das Paket CellChat, das die Zell-Zell-Kommunikation durch Analyse differentieller Liganden-Rezeptor-Wechselwirkungen zwischen Zellgruppen²² ableitet. Vor kurzem haben die Entwickler von CellChat ein detailliertes Schritt-für-Schritt-Protokoll für die allgemeine Verwendung^{veröffentlicht 24}, und dies ist eine hervorragende Ressource für Benutzer, wenn sie den folgenden Workflow durchlaufen und auf ihre Datensätze anwenden. Als Beispiel vergleicht der folgende Workflow die Interaktionen aller Hauptzellen in Wunden 1 und 14 Tage nach der Verwundung (DPW). Jeder Schritt wird nicht sehr detailliert beschrieben, da alle Schritte bereits in der offiziellen CellChat-Publikation²⁴ sowie in den hier verlinkten Tutorials beschrieben wurden:

Inferenz und Analyse der Zell-Zell-Kommunikation mittels
CellChat: https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html

Vergleichsanalyse mehrerer Datensätze mit CellChat:
https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html

Optionaler Schritt: Laden Sie bei Bedarf die gespeicherte RDS-Datei als Seurat Objekt:

dataset <- readRDS("dataset_post_Method3.rds")

1. Unterteilen Sie das ursprüngliche Dataset gemäß der DPW-Annotation in zwei Datasets.
   Idents(dataset) <- "DPW"
dataset_D1 <- subset(dataset, ident = "D1")
dataset_D14 <- subset(dataset, ident = "D14")
2. Definieren Sie die Anmerkung, mit der CellChat ausgeführt wird, --- verwenden Sie in diesem Fall die Hauptzellentypen.
   Idents(dataset_D1) <- "cell_types"
Idents(dataset_D14) <- "cell_types"
3. Erstellen Sie die CellChat-Objekte, und befolgen Sie den typischen CellChat-Arbeitsablauf. In den oben verlinkten Tutorials finden Sie detaillierte Referenzen für jeden Schritt.
   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. Visualisieren Sie die eingehenden und ausgehenden Interaktionsstärken in allen wichtigen Zelltypen zu jedem Zeitpunkt der Wundheilung (Ergänzende Abbildung 14).
   netAnalysis_signalingRole_scatter(cellchat_D1)
netAnalysis_signalingRole_scatter(cellchat_D14)
   Fibroblasten erhöhen ihre Wechselwirkungen zwischen D1 und D14 DPW dramatisch.
5. Zeigen Sie die Listen aller signifikanten abgeleiteten Zell-Zell-Kommunikationswege an.
   cellchat_D1@netP$pathways
cellchat_D14@netP$pathways
   Der Kollagenweg ist einer der signifikanten Signalwege sowohl auf D1 als auch auf D14 DPW.
6. Konzentrieren Sie sich auf den Kollagen-Signalweg und seine Interaktion mit Fibroblasten.
   pathways.show <- c("COLLAGEN")
7. Visualisieren Sie die Wechselwirkungen des Kollagen-Signalwegs zwischen Zelltypen mithilfe von Kreisdiagrammen (Ergänzende Abbildung 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. Visualisieren Sie die Wechselwirkungen des Kollagen-Signalwegs zwischen Zelltypen mithilfe von Sehnendiagrammen (Ergänzende Abbildung 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. Visualisieren Sie die Wechselwirkungen des COLLAGEN-Signalwegs mit Fibroblasten als Quellzellen (Ergänzende Abbildung 17).
   ANMERKUNG: Die Zelltypen in den cellchat-Objekten werden als IDs in der Reihenfolge aufgelistet, in der sie im ursprünglichen Seurat-Objekt zugewiesen wurden: 1 = Makrophage, 2 = Neutrophil, 3 = Fibroblast, 4 = Epithelzelle, 5 = Endothelzelle, 6 = T-Zelle, 7 = Glatte Muskelzelle.
   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. Visualisieren Sie die Beiträge jedes Liganden-Rezeptor-Paares im COLLAGEN-Signalweg mit Fibroblasten als Ausgangszellen.
    1. Verwendung von Blasendiagrammen (Ergänzende Abbildung 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. Verwendung von Sehnendiagrammen (Ergänzende Abbildung 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. Fokus auf die Col1a1-Cd44-Liganden-Rezeptor-Interaktion innerhalb des COLLAGEN-Signalwegs.
    ​LR.show <- "COL1A1_CD44"
12. Visualisieren Sie die Col1a1-Cd44-Liganden-Rezeptor-Wechselwirkungen zwischen Zelltypen mit Hilfe von Sehnendiagrammen (Ergänzende Abbildung 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. Führen Sie eine differenzielle CellChat-Analyse durch, indem Sie ein kombiniertes CellChat-Objekt generieren.
    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. Visualisieren Sie die Gesamtzahl und die relativen Stärken der Zell-Zell-Wechselwirkungen zu den Zeitpunkten der Wundheilung (Ergänzende Abbildung 21).
    gg1 <- compareInteractions(cellchat_D14_v_D1, show.legend = F)
gg2 <- compareInteractions(cellchat_D14_v_D1, show.legend = F, measure = "weight")
gg1 + gg2
15. Visualisieren Sie mit Hilfe eines Kreisdiagramms die unterschiedlichen Zell-Zell-Interaktionsstärken zwischen den einzelnen Zelltypen beim Übergang der Wunde von Tag 1 zu Tag 14 (Ergänzende Abbildung 22).
    netVisual_diffInteraction(cellchat_D14_v_D1, weight.scale = T, measure = "weight")
16. Visualisieren Sie anhand einer Heatmap die unterschiedlichen Zell-Zell-Interaktionsstärken zwischen den einzelnen Zelltypen beim Übergang der Wunde von Tag 1 zu Tag 14 (Ergänzende Abbildung 23).
    netVisual_heatmap(cellchat_D14_v_D1, measure = "weight")
17. Visualisieren Sie mit Hilfe eines Rangdiagramms die relativen Beiträge einzelner Signalwege zu Zell-Zell-Interaktionen mit Fibroblasten als Quellzellen an Tag 14 vs. Tag 1 (Ergänzende Abbildung 24).
    rankNet(cellchat_D14_v_D1, mode = "comparison", measure = "weight", sources.use = 3, targets.use = NULL, stacked = T, do.stat = TRUE)
18. Visualisieren Sie mit Hilfe von Blasendiagrammen die relativen Beiträge einzelner Liganden-Rezeptor-Paare im Kollagen-Signalweg mit Fibroblasten als Ausgangszellen an Tag 14 im Vergleich zu Tag 1 (Ergänzende Abbildung 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. Genau wie Seurat-Objekte können CellChat-Objekte als RDS-Dateien gespeichert und geöffnet werden.
    saveRDS(cellchat_D1, "cellchat_D1.rds")
saveRDS(cellchat_D14, "cellchat_D14.rds")
saveRDS(cellchat_D14_v_D1, "cellchat_D14_v_D1.rds")
20. Optionaler Schritt: CellChat-Objekte können auch aus RDS-Dateien geöffnet werden.
    cellchat_D1 <- readRDS("cellchat_D1.rds")
cellchat_D14 <- readRDS("cellchat_D14.rds")
cellchat_D14_v_D1 <- readRDS("cellchat_D14_v_D1.rds")

7. Beispiel für eine integrative Analyse durch Kombination mehrerer Einzelzelldatensätze

ANMERKUNG: Einzelzellige Datensätze werden oft in mehrere Dateien aufgeteilt, da sie in Stapeln oder Gruppen sequenziert wurden. Dieser Arbeitsablauf zeigt, wie zwei der fünf Chargen des Wundheilungsdatensatzes ²⁰ integriert werden. Die aktuellen Methoden zur Integration von Datensätzen werden durch die folgenden Seurat Vignetten beschrieben:

Einführung in die scRNA-seq-Integration:
https://satijalab.org/seurat/articles/integration_introduction

Integrative Analyse in Seurat v5:

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

Anmerkung: Es gibt zahlreiche Methoden zur Integration von Einzelzelldatensätzen, jede mit ihren eigenen Stärken und Schwächen. Weitere Informationen finden Sie im umfassenden Benchmark der Integrationsmethoden²⁵. Es ist wichtig, dass der Benutzer alle relevanten Dokumentationen liest, bevor er sich auf eine Integrationsmethode verlässt.

1. Wiederholen Sie alle Schritte in Methode 2 für einen weiteren Stapel des Datensatzes²⁰ von Hu et al. Im folgenden Protokoll wird Batch #3 verwendet. Die ergänzende R-Skriptdatei ist im Lieferumfang enthalten und kann zur Verarbeitung von Batch #3 verwendet werden (Ergänzende Datei 3: JoVE_Rscript_b3. R). Denken Sie daran, eine neue Variable für den Datensatz zu erstellen und zu verwenden --- im folgenden Code verwenden Sie "dataset_b3" für den Datensatz-Batch #3.
   1. Optionaler Schritt: Öffnen Sie bei Bedarf die beiden Datensätze als Seurat-Objekte aus ihren RDS-Dateien, die im Arbeitsverzeichnis gespeichert sind:
      dataset <- readRDS("dataset_post_Method2.rds")
dataset_b3 <- readRDS("dataset_b3_post_Method2.rds")
2. Weisen Sie jedem Datensatz eine neue Variable mit dem Namen "batch" zu, um den Ursprungsdatensatz in nachfolgenden Analysen zu beschriften.
   dataset@meta.data$batch <- "b1"
dataset_b3@meta.data$batch <- "b3"
3. Führen Sie eine Seurat-Zusammenführung der beiden Datasets durch, fügen Sie batchbasierte Zellen-ID-Anmerkungen hinzu, und führen Sie dann den standardmäßigen Seurat-Workflow für das zusammengeführte Dataset aus, wie in Methode 3 beschrieben.
   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. Führen Sie vor der Datenintegration Clustering und UMAP-Analysen für das kombinierte Dataset durch.
   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. Visualisieren Sie das UMAP-Diagramm nach Cluster- und Chargennummern (Ergänzende Abbildung 26).
   DimPlot(dataset_merged, reduction = "umap.unintegrated", group.by = c("seurat_clusters", "batch"))
6. Zeigen Sie die Verteilung der Zellennummern in jedem Cluster entsprechend der Batchnummer an.
   table(dataset_merged$batch, dataset_merged$seurat_clusters)
   ANMERKUNG: Aus dem UMAP-Diagramm und der Tabelle geht hervor, dass es für diese beiden Datensätze keine signifikanten Batch-Effekte gibt. Hinweise auf Batch-Effekte würden sich in unerwarteten Diskrepanzen in der Clusterverteilung zwischen den beiden Datensätzen manifestieren, was bedeuten könnte, dass es potenzielle technische Unterschiede zwischen den Datensätzen gibt, die tatsächliche biologische Ähnlichkeiten außer Kraft setzen.
7. Führen Sie die Integration von Seurat-Daten mit der RPCA-Methode durch. Weitere Informationen zu dieser und anderen Methoden der Datenintegration finden Sie in der oben verlinkten Seurat Vignette.
   dataset_merged <- IntegrateLayers(
object = dataset_merged, method = RPCAIntegration,
orig.reduction = "pca", new.reduction = "integrated.rpca",
verbose = TRUE
)
8. Führen Sie nach der Datenintegration Clustering und UMAP-Analysen für das kombinierte Dataset durch.
   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. Visualisieren Sie den UMAP-Plot nach der Integration nach Cluster- und Chargennummern (Ergänzende Abbildung 27).
   DimPlot(dataset_merged, reduction = "umap.rpca", group.by = c("seurat_clusters","batch"))
10. Zeigt die Verteilung der Zellennummern in jedem Cluster entsprechend der Batchnummer nach der Integration an.
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    Aus dem UMAP-Diagramm und der Tabelle der integrierten Daten ergibt sich nun eine hervorragende Überlappung zwischen den beiden Chargen über verschiedene Cluster hinweg. Interessanterweise führte die Integration der Daten zur Identifizierung eines zusätzlichen Clusters mit den gleichen Clustering-Parametern.
11. Nach der Integration des Datensatzes und vor nachgelagerten Analysen müssen die Layer des zusammengeführten Datensatzes verbunden werden.
    dataset_merged <- JoinLayers(dataset_merged)
12. Speichern Sie das Seurat-Objekt des Datensatzes als RDS-Datei in das Arbeitsverzeichnis.
    saveRDS(dataset_merged, "merged_dataset_post_Method7.rds")

Results

Beginnend mit Methode #2 führt das Protokoll die Schritte zum Laden und Durchführen von Qualitätskontrollschritten für einen Einzelzell-Wundheilungsdatensatz durch. Nach der Erstellung des Seurat-Objekts (Schritt 2.6.2) werden in einer Reihe von Schritten die beiden Assays innerhalb des Datensatzes (RNA und Protein; Schritte 2.6.3-2.6.7) zusammengeführt und der Proteinassay anhand von räumlich-zeitlichen Barcodes (Schritte 2.6.8-2.6.9) dekomplexiert. Die Decomplexing-Funktion weist jeder Zelle im Datensatz mehrere Metadatenbeschriftungen zu, darunter "barcodes_maxID", das den wahrscheinlichsten räumlich-zeitlichen Barcode jeder Zelle identifiziert (Schritt 2.6.10). In Schritt 2.6.11 wird die Violin-Plot-Funktion ausgeführt, um die Verteilung der detektierten Gene in Zellen anhand ihrer gemultiplexten Barcodes zu visualisieren. Das repräsentative Ergebnis dieses Schritts (Ergänzende Abbildung 1) zeigt, dass es eine ziemlich gleichmäßige Verteilung der detektierten Gene für jeden Barcode gibt, was für die Integrität des Datensatzes und die nachgelagerte Analyse von Wundheilungszeitpunkten wichtig ist. Nach der Zuweisung des entsprechenden Etiketts zu den Protein-Barcodes (Schritt 2.6.12) zeigt das Protokoll dann, wie Qualitätskontrollschritte am RNA-Assay des Datensatzes durchgeführt werden, beginnend mit der Berechnung des Prozentsatzes der mitochondrialen Gene in jeder Zelle (Schritt 2.9). In Schritt 2.10 wird die Funktion Merkmalsstreudiagramm ausgeführt, um die Verteilung der detektierten Gene, die Anzahl der RNA und den mitochondrialen Prozentsatz in allen Zellen zu visualisieren. Die repräsentativen Ergebnisse dieses Schritts (Ergänzende Abbildung 2) zeigen, dass es eine Reihe von Zellen mit hohem mitochondrialen Gehalt gibt, was mit niedrigen RNA-Zahlen korreliert und tote oder sterbende Zellen identifiziert. Nach der Entnahme von Zellen mit niedrigen RNA-Zahlen und großen mitochondrialen Gehalten (Schritt 2.11) wird in Schritt 2.12 eine weitere Merkmals-Scatterplot-Funktion für den Teilmengendatensatz durchgeführt, und das repräsentative Ergebnis dieses Schritts (Ergänzende Abbildung 3) zeigt, dass die Verteilung der nachgewiesenen Gene und des prozentualen Anteils an mitochondrialer RNA pro Zelle nun normaler ist, was den Weg für robuste nachgelagerte Analysen ebnet. Als nächstes beschreibt das Protokoll die Verwendung der scDblFinder-Funktion, um wahrscheinliche Dubletten im Datensatz zu identifizieren, und weist jeder Zelle eine neue Metadaten namens "scDblFinder.score" zu (Schritte 2.13-2.14). In Schritt 2.15 wird die Violin-Plot-Funktion ausgeführt, um die Verteilung der Dublett-Werte im Datensatz zu visualisieren, und das repräsentative Ergebnis dieses Schritts (Ergänzende Abbildung 4) zeigt, dass es eine Reihe von Zellen mit relativ hohen Dublett-Werten gibt und dass 0,25 ein natürlicher Grenzwert zu sein scheint, oberhalb dessen es eine Population wahrscheinlicher Dubletten gibt. Daher wird in den folgenden Schritten dieser Parameter verwendet, um das Dataset in Zellen unterhalb des Grenzwerts (Schritt 2.16) zu unterteilen, wodurch die Qualitätskontrollschritte für dieses einzellige Dataset abgeschlossen sind.

Beginnend mit Methode #3 führt das Protokoll durch die Schritte zur Analyse des qualitätskontrollierten Einzelzell-Wundheilungsdatensatzes unter Verwendung des Seurat-Pakets und -Workflows. Nach der Normalisierung und Skalierung der RNA-Daten wird die PCA-Analyse durchgeführt (Schritt 3.1). In Schritt 3.2 wird die Funktion elbow plot verwendet, um das Ausmaß der Variation des Datensatzes in Bezug auf die ersten 50 PCA-Dimensionen und das repräsentative Ergebnis dieses Schritts (Ergänzende Abbildung 5) zeigt, dass ein Großteil der größten Variation innerhalb der ersten 13 Dimensionen auftritt, die durch die Biegung im Diagramm gekennzeichnet sind. Das Protokoll zeigt dann, wie man Nachbarn findet und Zellclustering (Schritt 3.3) und UMAP-Dimensionsreduktion (Schritt 3.4) des Datensatzes unter Verwendung der ersten 13 PCA-Dimensionen und eines relativ niedrigen Clustering-Auflösungsparameters von 0,1 durchführt, die beide ausgewählt wurden, um die verallgemeinerbarsten Hauptzelltypen in Wunden zu identifizieren. In Schritt 3.5 wird die Funktion "Dimensionale Darstellung" ausgeführt, um die Clusterbildung der Zellen in einem UMAP-Diagramm zu visualisieren, und das repräsentative Ergebnis dieses Schritts (Abbildung 1) zeigt, dass alle Zellen im Datensatz um 8 wichtige farbcodierte Seurat-Clustergruppen gruppiert sind, mit leicht unterschiedlichen UMAP-Diagrammen, die von einem Computer mit Windows (links) und MacOS (rechts) erstellt wurden. In Schritt 3.6 wird eine weitere dimensionale Plot-Funktion ausgeführt, um die Wundzeit-/Raum-Annotation der Zellen und das repräsentative Ergebnis dieses Schritts (Abbildung 2) zeigt, dass alle Zellen im Datensatz entsprechend ihrem Zeit-/Raum-Ursprung verteilt sind, ohne offensichtliche Cluster-Bildung gemäß der Zeit-/Raum-Annotation. Das Protokoll beschreibt als nächstes, wie man Listen von differentiell exprimierten Genen erhält und sie in einer Textdatei speichert (Schritt 3.8), die Datentabelle in einer Tabelle öffnet und verschiedene Filterschritte durchführt, um die am besten bewerteten Clustermarker für jeden Zellcluster zu erhalten (Schritte 3.9-3.10.6). Das repräsentative Ergebnis dieser Schritte (Ergänzende Tabelle 1) ist die endgültige Tabellenkalkulationsdatei, die die vollständige Ausgabe der differentiell exprimierten Gene enthält, während ein weiteres repräsentatives Ergebnis (Ergänzende Tabelle 2) ist eine vereinfachte Tabelle, die die 5 wichtigsten hochregulierten und exprimierten Gene für jeden Seurat-Cluster zeigt. Das Protokoll beschreibt dann, wie ein webbasiertes Tool zur Analyse der funktionellen Anreicherung namens EnrichR verwendet wird, um mutmaßliche Zelltypen anhand der Top-Cluster-Markergene (Schritte 3.11-3.12) und der repräsentativen Ergebnisse dieser Schritte (Abbildung 3) sind beschnittene Screenshots der EnrichR-Ausgaben, die die am stärksten angereicherten Zelltypen für jeden der acht Zellcluster zeigen. Das Protokoll weist dann allen Zellen in den jeweiligen Seurat-Clustern ein neues Metadaten-Label namens "cell_types" zu, entsprechend ihren am stärksten angereicherten Zelltyp-Annotationen (Schritt 3.14). In Schritt 3.15 wird die Funktion "Dimensionales Zeichnen" ausgeführt, um die umbenannten Zellcluster als Zelltyp-Annotationen in einem UMAP-Diagramm zu visualisieren und die repräsentativen Ergebnisse dieses Schritts (Abbildung 4) zeigte, dass sich alle Zellen im Datensatz um die wichtigsten farbcodierten Zelltypen gruppierten. In Schritt 3.16 wurde die Feature-Plot-Funktion verwendet, um die Lokalisierung der Top-Cluster-Markergene (aus der ergänzenden Tabelle 2) auf einer Reihe von UMAP-Diagrammen zu visualisieren, und die repräsentativen Ergebnisse (Abbildung 5) sind ein Raster von UMAP-Diagrammen, die die hohe Expression der Top-Zell-Markergene innerhalb ihrer jeweiligen Hauptzelltyp-Cluster-Standorte zeigen. In den Schritten 3.17 und 3.18 wurde die Dotplot-Funktion durchgeführt, um die relativen Expressionsniveaus von Top-Cluster-Marker-Genen in Zellen zu visualisieren, zunächst gruppiert nach ihren ursprünglichen Seurat-Cluster-Nummern (Schritt 3.17) und dann gruppiert nach annotierten Zelltyp-Beschriftungen (Schritt 3.18). Die repräsentativen Ergebnisse dieser Schritte bestätigten die hohe Expression der Top-Zell-Markergene nur in ihren jeweiligen Seurat-Clustern (Ergänzende Abbildung 6) und nur in ihren jeweiligen Hauptzelltypen (Abbildung 6). Der nächste Schritt im Protokoll vereinfacht die ursprünglichen räumlich-temporalen Protein-basierten Markierungen in streng zeitliche Annotationen, die Zellen basierend auf den Tagen nach der Verletzung (DPW) identifizieren, aus denen sie stammen. In Schritt 3.20 wird die Funktion "Dimensionales Zeichnen" ausgeführt, um die Zellen als DPW-Annotationen in einem UMAP-Diagramm zu visualisieren, und die repräsentativen Ergebnisse dieses Schritts (Ergänzende Abbildung 7) zeigte die Lokalisierung von Wund-Zeitverlauf-Annotationen über den Einzelzell-Wundheilungsdatensatz. Erwartungsgemäß dominierten die Annotationen von Tag 1 (D1) in den Neutrophilen- und Makrophagenclustern, während die späteren Wundheilungszeitpunkte bei anderen Zelltypen stärker vertreten waren. In den folgenden Schritten des Protokolls wurden gestapelte Balkendiagramme verwendet, um zunächst die Anteile von DPW über verschiedene Zelltypen (Schritt 3.22) und dann die Anteile von Zelltypen über verschiedene Zeitpunkte hinweg zu visualisieren (3.23). Die repräsentativen Ergebnisse dieser Schritte sind Proportionsdiagramme, die die relative Anzahl von DPW-Zellen in jeder Hauptzellentypkategorie (Ergänzende Abbildung 8) und die relative Anzahl der Hauptzelltypen in jeder DPW-Kategorie (Abbildung 7). Diese Ergebnisse bestätigten die bekannte zelluläre Kaskade der Hautwundheilung, bei der Immunzellen (Neutrophile und Makrophagen) die frühen Zeitpunkte während der Entzündungsphase dominieren und die anderen Zelltypen (Epithelzellen und Endothelzellen) während der proliferativen Phase zu erscheinen beginnen, wobei Fibroblasten zu den späteren Zeitpunkten während der Wundauflösung besonders dominant sind.

Beginnend mit Methode #4 beschreibt das Protokoll die Schritte für die Verwendung von Seurat zur Fokussierung auf einen einzelnen Hauptzelltyp im Einzelzelldatensatz, um potenzielle zelluläre Subtypen während der Wundheilung zu identifizieren. Das Protokoll konzentriert sich auf Fibroblasten, die sich zunächst in zwei Seurat-Clustern gruppierten, bevor sie zu einer einzigen Kategorie zusammengefasst wurden, und beschreibt, wie ein neues Seurat-Objekt erstellt wird, das nur die Fibroblasten aus dem ursprünglichen Datensatz enthält (Schritt 4.1). Der Seurat-Workflow wird an diesem Fibroblasten-spezifischen Datensatz durchgeführt (Schritte 4.2-4.4), wobei Schritt 4.2 zu einem Ellbogendiagramm führt (Ergänzende Abbildung 9), das zeigt, dass ein Großteil der Hauptvariation im Fibroblasten-Datensatz innerhalb der ersten 9 PCA-Dimensionen auftritt. In Schritt 4.5 wird die dimensionale Plot-Funktion durchgeführt, um die Clusterbildung der Zellen in einem UMAP-Diagramm zu visualisieren, und die repräsentativen Ergebnisse dieses Schritts (Abbildung 8) zeigten die Fibroblasten im Datensatz, die um die 3 farbcodierten Zellsubtypen gruppiert waren. Die Visualisierung des Fibroblasten-Datensatzes gemäß ihrer DPW-Annotation (Schritt 4.6) führte zu einem UMAP-Diagramm (Ergänzende Abbildung 10), das die Fibroblasten im Datensatz gemäß ihrer DPW-Annotation verteilt. Das Protokoll beschreibt als nächstes, wie man Listen von differentiell exprimierten Genen erhält und sie in einer Textdatei speichert (Schritt 4.7), die Datentabelle in Excel öffnet und verschiedene Filterschritte durchführt, um die bestplatzierten Clustermarker für jeden Zellcluster zu erhalten (Schritt 4.8), und eine neue Variable zuweist, die die Top-Fibroblastenmarkergene mit dem Namen "FB_type_marker" auflistet (Schritt 4.9). In Schritt 4.10 wird die Dot-Plot-Funktion verwendet, um die Gene in der Liste im Nur-Fibroblasten-Datensatz zu visualisieren, indem die Variable "FB_type_marker" im Parameter features aufgerufen wird, und die repräsentativen Ergebnisse dieses Schritts (Abbildung 9) sind Dotplots, die eine hohe Expression von Fibroblasten-Subtyp-Markern nur in ihren jeweiligen Cluster-Kategorien (oben), aber fair verteilt über die DPW-Kategorien (unten) bestätigen. In Schritt 4.11 wird die Variable "Gleiche Merkmale" aufgerufen, um die Fibroblasten-Markergene im gesamten Wundheilungsdatensatz zu visualisieren, und das repräsentative Ergebnis (Ergänzende Abbildung 11) ist ein Punktdiagramm, das die hohe Expression von Fibroblasten-Subtyp-Markern hauptsächlich im ursprünglichen Fibroblasten bestätigt. In den folgenden Schritten des Protokolls wurden schließlich gestapelte Balkendiagramme verwendet, um zunächst die Anteile von DPW über die drei Fibroblasten-Subtypen hinweg zu visualisieren (Schritt 4.12) und dann die Anteile der Fibroblasten-Subtypen über verschiedene Zeitpunkte hinweg zu visualisieren (Schritt 4.13). Die repräsentativen Ergebnisse dieser Schritte sind Proportionsdiagramme, die die relative Anzahl der DPW-Zellen in jeder Fibroblasten-Subtyp-Kategorie (Ergänzende Abbildung 12) und die relative Anzahl der Fibroblasten-Subtypen in jeder DPW-Kategorie (Ergänzende Abbildung 13) zeigen. Diese Ergebnisse deuten auf eine signifikante Veränderung der Fibroblasten-Subtyp-Anteile im Zeitverlauf der Heilung hin, wobei der erste Fibroblasten-Subtyp (Cluster 0) in frühen Stadien der Wunde stark dominant war (D1 und D3), der zweite Subtyp (Cluster 1) während der Wundauflösung (D14) und der dritte Subtyp (Cluster 2) während der proliferativen Phase der Wundheilung am höchsten war (D7).

Beginnend mit Methode #5 führt das Protokoll die Schritte zur Analyse eines einzelligen Wundheilungsdatensatzes mit der Modulbewertungsfunktion in Seurat durch. Das Protokoll beschreibt zunächst die Schritte der Verwendung einer Tabulator-Textdatei zum Hochladen von Gensätzen in Variablen in R (Schritte 5.1-5.2), gefolgt von der Anwendung der Modul-Scoring-Funktion auf drei Gensätze, die sich auf die drei Hauptphasen der Wundheilung beziehen (Schritt 5.3). In Schritt 5.4 wird die Punktdiagrammfunktion verwendet, um die aggregierten Modulwerte über zwei verschiedene Metadatenkategorien hinweg zu visualisieren, und die repräsentativen Ergebnisse für diesen Schritt (Abbildung 10) sind Punktdiagramme, die die durchschnittliche Expression der wichtigsten Heilungsphasenmodule in Zellen in der Kategorie Tage nach der Verletzung (DPW, rechts) und in der Kategorie der Hauptzelltypen (links) zeigen. Diese Ergebnisse zeigen, dass die Anwendung von auf Massensequenzierung basierenden Genexpressionsprofilen auf Einzelzellexpressionsdatensätze in einer Pseudo-Bulk-Manier eine leistungsfähige Methode für vergleichende bioinformatische Ansätze ist, indem bereits veröffentlichte Datensätze im Bereich der Wundheilung verwendet werden.

Beginnend mit Methode #6 führt das Protokoll die Schritte zur Analyse eines von Seurat abgeleiteten Einzelzell-Wundheilungsdatensatzes unter Verwendung des CellChat-Pakets und -Workflows gemäß einer spezifischen wissenschaftlichen Fragestellung des Vergleichs von Zellen aus frühen im Vergleich zu Wunden in der späten Phase durch. Das Protokoll unterteilt zunächst den gesamten Seurat-Datensatz in zwei Zeitpunkte nach der Verletzung, einen während der Entzündungsphase (Tag 1 (D1)) und den anderen während der Wundauflösung (Tag 14 (D14)) (Schritt 6.1). Es werden zwei CellChat-Objekte erstellt, und das Protokoll durchläuft alle typischen Funktionen des CellChat-Protokolls, um alle mutmaßlichen Wechselwirkungen zwischen den Zelltypen zu berechnen, die in Methode #3 des Protokolls (Schritte 6.2-6.3) identifiziert wurden. In Schritt 6.4 wird die Signalsignal-Scatterplot-Funktion durchgeführt, um die ein- und ausgehenden Interaktionsstärken in allen wichtigen Zelltypen zu jedem Zeitpunkt der Wundheilung zu visualisieren. Die repräsentativen Ergebnisse dieses Schritts (Ergänzende Abbildung 14) sind Streudiagramme, die die Stärken der eingehenden (y-Achse) und ausgehenden (x-Achse) Wechselwirkungen für die wichtigsten Zelltypen zu den Zeitpunkten D1 (links) und D14 (rechts) zeigen. Diese Ergebnisse zeigten, dass Immunzellen wie Neutrophile und Makrophagen während der Entzündungsphase die höchsten Zell-Zell-Interaktionsstärken aufwiesen, während Fibroblasten die Zell-Zell-Interaktionen während der Wundauflösung dominierten, was jahrzehntelange Wundheilungsforschung bestätigt. Die folgenden Schritte konzentrieren sich bei der Analyse auf einen der signifikant angereicherten Signalwege, den Kollagenweg (Schritte 6.5-6.6). In Schritt 6.7 wird die Kreisdiagrammfunktion durchgeführt, um die Wechselwirkungen des Kollagen-Signalwegs zwischen den Zelltypen zu den beiden Zeitpunkten zu visualisieren. Die repräsentativen Ergebnisse dieses Schritts (Ergänzende Abbildung 15) sind Kreisdiagramme, die den abgeleiteten Kollagenweg zeigen, der Wechselwirkungen zwischen allen Zelltypen an D1 (links) und D14 (rechts) signalisiert. In Schritt 6.8 werden die gleichen Wechselwirkungen mit der Sehnendiagrammfunktion visualisiert, wobei die repräsentativen Ergebnisse (Ergänzende Abbildung 16) sind Sehnendiagramme, die den abgeleiteten Kollagenweg zeigen, der Wechselwirkungen zwischen allen Zelltypen zu jedem Zeitpunkt signalisiert. Erwartungsgemäß zeigten diese Ergebnisse, dass Fibroblasten die primären Quellzellen für den Kollagen-Signalweg waren, obwohl der Informationsfluss bei D1 im Vergleich zu D14 stärker auf Immunzellen beschränkt war. Um sich auf den Fibroblasten als Quellzelle in den Zell-Zell-Interaktionen zu konzentrieren, wiederholt Schritt 6.9 die Sehnendiagrammfunktion, indem er einen Quellzellparameter hinzufügt, und die repräsentativen Ergebnisse (Ergänzende Abbildung 17) sind Sehnendiagramme, die den abgeleiteten Kollagenweg zeigen, der Wechselwirkungen mit Fibroblasten als Quellzellen zu jedem Zeitpunkt signalisiert. In Schritt 6.10 werden zwei Funktionen ausgeführt, um die Beiträge jedes Liganden-Rezeptor-Paares im Kollagen-Signalweg mit Fibroblasten als Quellzellen zu visualisieren, eine mit Hilfe von Blasendiagrammen (Schritt 6.10.1) und die andere mit Hilfe von Sehnendiagrammen (Schritt 6.10.2). Die repräsentativen Ergebnisse zeigen die abgeleiteten Beiträge jedes Liganden-Rezeptor-Paares zur Signalgebung des Kollagensignalwegs mit Fibroblasten als Quellzellen zu den Zeitpunkten D1 (links) und D14 (rechts) unter Verwendung beider Blasendiagramme (Ergänzende Abbildung 18) und Sehnendiagramme (Ergänzende Abbildung 19). Diese Ergebnisse zeigten, dass bei D1 der Kollagenweg, der von Fibroblasten ausgeht, auf Neutrophile und Makrophagen mit einer Dominanz von Cd44- und Sdc4-Rezeptoren beschränkt war, während in D14 andere Zellen über eine Vielzahl von Rezeptoren, einschließlich Integrinen, als Empfänger fungierten. Um sich auf die Col1a1-Cd44-Liganden-Rezeptor-Interaktion zu konzentrieren, die starke Stärken bei Fibroblasten-Interaktionen zeigte, wird ein Parameter eingestellt (Schritt 6.11) und dann in Schritt 6.12 in einer Sehnendiagrammfunktion verwendet, um diese spezielle Ligand-Rezeptor-Interaktion zwischen allen Zelltypen zu visualisieren, mit den repräsentativen Ergebnissen (Ergänzende Abbildung 20) sind Sehnendiagramme, die die abgeleiteten Col1a1-Cd44-Liganden-Rezeptor-Wechselwirkungen zwischen allen Zelltypen zu den Zeitpunkten D1 (links) und D14 (rechts) zeigen. Diese Ergebnisse zeigten, dass diese Interaktion in D1 auf Fibroblasten als Quellzellen beschränkt ist, während in D14 Makrophagen und glatte Muskelzellen ebenfalls als Ausgangszellen fungieren. Als nächstes beschreibt das Protokoll, wie eine differentielle CellChat-Analyse durchgeführt wird, indem zunächst die CellChat-Objekte D1 und D14 zusammengeführt werden (Schritt 6.13). In Schritt 6.14 wird die Funktion "Wechselwirkungen vergleichen" durchgeführt, um die Gesamtzahl und die relativen Stärken der Zell-Zell-Wechselwirkungen zwischen den beiden Wundheilungszeitpunkten und die repräsentativen Ergebnisse (Ergänzende Abbildung 21) sind die resultierenden Balkendiagramme, die die Gesamtzahl (links) und Stärken (rechts) der abgeleiteten Wechselwirkungen in Zellen zeigen, die D1- und D14-Wunden umfassen, mit einer höheren Anzahl von Wechselwirkungen in D14 im Gegensatz zu höheren relativen Stärken von Wechselwirkungen in D1. In den Schritten 6.15 und 6.16 werden zwei Funktionen verwendet, um die differentiellen Zell-Zell-Interaktionsstärken zwischen den einzelnen Zelltypen beim Übergang der Wunde von Tag 1 zu Tag 14 mit ihren jeweiligen repräsentativen Ergebnissen zu visualisieren, wobei die erste ein Kreisdiagramm ist (Schritt 6.15, Ergänzende Abbildung 22) und zweitens eine Heatmap (Schritt 6.16, Ergänzende Abbildung 23), wobei erhöhte Wechselwirkungen in D14 im Vergleich zu D1 in rot und verminderte Wechselwirkungen in blau dargestellt wurden. Erwartungsgemäß sind die neutrophilen- und makrophagenvermittelten Wechselwirkungen in D1 und die Fibroblasten-vermittelten Wechselwirkungen in D14 erhöht. In Schritt 6.17 wird die Rangfolgefunktion verwendet, um ein Diagramm zu erstellen, das die relativen Beiträge einzelner Signalwege zu Zell-Zell-Interaktionen mit Fibroblasten als Quellzellen bei D14 im Vergleich zu D1 einstuft, und die repräsentativen Ergebnisse (Ergänzende Abbildung 24) zeigen das resultierende Rangdiagramm, wobei D1 oben rot und D14 unten blau dargestellt wird, wobei mehrere Pfade ausschließlich in D1 oder D14 dargestellt sind und viele andere einen Aktivierungsgradienten aufweisen. Schließlich werden in Schritt 6.18 zwei Blasendiagrammfunktionen verwendet, um die relativen Beiträge einzelner Liganden-Rezeptor-Paare im Kollagen-Signalweg mit Fibroblasten als Quellzellen bei D14 im Vergleich zu D1 zu zeigen, mit den entsprechenden repräsentativen Ergebnissen (Ergänzende Abbildung 25), die erhöhte (links) und verringerte (rechts) Signalpaarungen in D14 im Vergleich zu D1 über die vielen Zell-Zell-Interaktionen auf der x-Achse zeigen. Wie erwartet, wiesen Fibroblasten in D14-Wunden viel mehr erhöhte Liganden-Rezeptor-Paar-Interaktionen über mehrere Empfängerzellen hinweg auf, verglichen mit D1-Wunden, bei denen die Kommunikation während der Entzündungsphase stärker auf Neutrophile und Makrophagen beschränkt war.

Beginnend mit Methode #7 führt das Protokoll die Schritte zur Integration von zwei Einzelzell-Wundheilungsdatensätzen mit Seurat durch. Das Protokoll beschreibt zunächst die Schritte zum Zusammenführen von zwei Stapeln der veröffentlichten Einzelzelldatensätze und zum Anwenden des standardmäßigen Seurat-Workflows auf den zusammengeführten Datensatz (Schritte 7.1-7.4). In Schritt 7.5 wird die Funktion Dimensional Plot verwendet, um den UMAP-Plot anhand von Cluster- und Chargennummern des zusammengeführten, aber noch nicht integrierten Wundheilungsdatensatzes zu visualisieren. Die repräsentativen Ergebnisse aus diesem Schritt (Ergänzende Abbildung 26) sind UMAP-Diagramme, die die Verteilung von Seurat-Clustern (links) und Chargennummern (rechts) visualisieren und zeigen, dass es für diese beiden Datensätze vor der Datenintegration keine signifikanten Batch-Effekte zu geben scheint. Das Protokoll führt dann die Datenintegration mit der RPCA-Methode und dem anschließenden Seurat-Workflow des integrierten Datensatzes durch (Schritte 7.7-7.8). In Schritt 7.9 wird die Funktion Dimensional Plot verwendet, um den UMAP-Plot nach Cluster- und Chargennummern des integrierten Wundheilungsdatensatzes zu visualisieren. Die repräsentativen Ergebnisse aus diesem Schritt (Ergänzende Abbildung 27) sind UMAP-Diagramme, die die Verteilung von Seurat-Clustern (links) und Chargennummern (rechts) visualisieren und zeigen, dass es jetzt eine noch größere Überlappung zwischen den beiden Chargen über verschiedene Cluster hinweg gab. Die Ergebnisse zeigen auch die Entstehung eines zusätzlichen Clusters nach der Integration der Daten, was auf die erhöhte Fähigkeit hindeuten könnte, potenziell signifikante Zellsubtypen zu identifizieren, nachdem die technischen Effekte der Datenchargen kontrolliert wurden.

Abbildung 1: UMAP-Diagramm, das alle Zellen im Datensatz zeigt, die um 8 wichtige farbcodierte Clustergruppen gruppiert sind. Die Ergebnisse wurden von einem Computer mit Windows (links) und MacOS (rechts) erzielt. Diese Zahl entspricht Schritt 3.5. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 2: UMAP-Diagramm, das alle Zellen im Datensatz nach ihrem Zeit-/Raum-Ursprung ausbreitet, ohne offensichtliche Cluster-Bildung gemäß der Zeit-/Raum-Annotation. Diese Zahl entspricht Schritt 3.6. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 3: Beschnittene Screenshots der EnrichR-Ausgaben, die die am besten angereicherten Zelltypen für jeden Zellcluster zeigen. Diese Zahl entspricht Schritt 3.13. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 4: UMAP-Diagramm, das alle Zellen im Datensatz zeigt, die um die wichtigsten farbcodierten Zelltypen gruppiert sind. Diese Zahl entspricht Schritt 3.15. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 5: Raster von UMAP-Diagrammen, die die hohe Expression der Top-Zell-Markergene innerhalb der wichtigsten Zelltyp-Cluster zeigen. Diese Zahl entspricht Schritt 3.16. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 6: Punktdiagramme, die die hohe Expression der Top-Zell-Markergene nur in ihren jeweiligen Hauptzelltypen bestätigen. Diese Zahl entspricht Schritt 3.18. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 7: Proportionsdiagramm, das die relative Anzahl der wichtigsten Zellentypen in jeder DPW-Kategorie zeigt. Diese Zahl entspricht Schritt 3.23. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 8: UMAP-Diagramm mit Fibroblasten im Datensatz, die um die 3 farbcodierten Zellsubtypen gruppiert sind. Diese Zahl entspricht Schritt 4.5. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 9: Punktdiagramme, die eine hohe Expression von Fibroblasten-Subtyp-Markern nur in ihren jeweiligen Cluster-Kategorien, aber fair verteilt über die DPW-Kategorien bestätigen. Diese Zahl entspricht Schritt 4.10. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Abbildung 10: Punktdiagramme, die die durchschnittliche Expression der wichtigsten Heilungsphasenmodule über Zellen pro DPW und pro Hauptzelltyp zeigen. Diese Zahl entspricht Schritt 5.4. Bitte klicken Sie hier, um eine größere Version dieser Abbildung anzuzeigen.

Ergänzende Abbildung 1: Ergebnisse, die zeigen, dass es eine ziemlich gleichmäßige Verteilung der detektierten Gene für jeden Barcode gibt, was für die Integrität des Datensatzes und die nachgelagerte Analyse von Wundheilungszeitpunkten wichtig ist. Diese Zahl entspricht Schritt 2.6.11. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 2: Streudiagramme, die zeigen, dass es eine Reihe von Zellen mit hohem mitochondrialen Gehalt gibt, die mit niedrigen RNA-Zahlen korrelieren, --- es sich um tote oder sterbende Zellen handelt. Diese Zahl entspricht Schritt 2.10. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 3: Streudiagramme, die zeigen, dass die Verteilung der detektierten Gene und der prozentuale Anteil an mitochondrialer RNA pro Zelle nun normaler ist, was den Weg für robuste nachgelagerte Analysen ebnet. Diese Zahl entspricht Schritt 2.12. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 4: Violindiagramm, das zeigt, dass es eine Reihe von Zellen mit relativ hohem Dublettenwert gibt und dass 0,25 ein natürlicher Grenzwert zu sein scheint, oberhalb dessen es eine Population von wahrscheinlichen Dubletten gibt. Diese Zahl entspricht Schritt 2.15. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 5: Ellenbogendiagramm, das zeigt, dass ein Großteil der Hauptvariation innerhalb der ersten 13 Dimensionen auftritt. Diese Zahl entspricht Schritt 3.2. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 6: Punktdiagramm, das die hohe Expression der Top-Zell-Markergene nur in ihren jeweiligen Seurat-Clustern bestätigt. Diese Zahl entspricht Schritt 3.17. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 7: UMAP-Diagramm, das die Lokalisierung von Wundzeit-Verlauf-Annotationen im gesamten Wundheilungsdatensatz zeigt. Diese Zahl entspricht Schritt 3.20. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 8: Proportionsdiagramm, das die relative Anzahl der DPW-Zellen in jeder Hauptzellentypkategorie zeigt. Diese Zahl entspricht Schritt 3.22. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 9: Ellenbogendiagramm, das zeigt, dass ein Großteil der Hauptvariation im Fibroblastendatensatz innerhalb der ersten 9 Dimensionen auftritt. Diese Zahl entspricht Schritt 4.2. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 10: UMAP-Diagramm mit Fibroblasten im Datensatz, die entsprechend ihrer DPW-Annotation verteilt sind. Diese Zahl entspricht Schritt 4.6. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 11: Punktdiagramm, das eine hohe Expression von Fibroblasten-Subtyp-Markern hauptsächlich im ursprünglichen Fibroblasten-Cluster bestätigt. Diese Zahl entspricht Schritt 4.11. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 12: Proportionsdiagramm mit der relativen Anzahl von Fibroblasten-Subtypen in jeder DPW-Kategorie. Diese Zahl entspricht Schritt 4.12. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 13: Proportionsdiagramm, das die relative Anzahl von Fibroblasten in DPW in jeder Fibroblasten-Subtypkategorie zeigt. Diese Zahl entspricht Schritt 4.13. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 14: Streudiagramme, die die Stärken der eingehenden (y-Achse) und ausgehenden (x-Achse) Wechselwirkungen für die wichtigsten Zelltypen zu den Zeitpunkten Tag 1 (D1, links) und Tag 14 (D14, rechts) zeigen. Diese Zahl entspricht Schritt 6.4. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 15: Kreisdiagramme, die die abgeleiteten Wechselwirkungen des Kollagenwegs zwischen allen Zelltypen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.7. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 16: Sehnendiagramme, die die abgeleiteten Wechselwirkungen des Kollagenwegs zwischen allen Zelltypen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.8. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 17: Sehnendiagramme, die die abgeleiteten Wechselwirkungen des Kollagenwegs mit Fibroblasten als Quellzellen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.9. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 18: Blasendiagramme, die die abgeleiteten Beiträge jedes Liganden-Rezeptor-Paares zur Signalgebung des Kollagensignalwegs mit Fibroblasten als Quellzellen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.10.1. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 19: Sehnendiagramme, die die abgeleiteten Beiträge jedes Liganden-Rezeptor-Paares zur Signalübertragung des Kollagenwegs mit Fibroblasten als Quellzellen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.10.2. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 20: Sehnendiagramme, die die abgeleiteten Col1a1-Cd44-Liganden-Rezeptor-Wechselwirkungen zwischen allen Zelltypen in jeder DPW-Kategorie zeigen. Diese Zahl entspricht Schritt 6.12. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 21: Balkendiagramme, die die Anzahl (links) und Stärke (rechts) der abgeleiteten Wechselwirkungen bei Wunden an Tag 1 und Tag 14 zeigen. Diese Zahl entspricht Schritt 6.14. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 22: Kreisdiagramm, das die unterschiedlichen Zell-Zell-Interaktionsstärken zwischen den einzelnen Zelltypen beim Übergang der Wunde von Tag 1 (blau) zu Tag 14 (rot) DPW zeigt. Diese Zahl entspricht Schritt 6.15. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 23: Heatmap, die die unterschiedlichen Zell-Zell-Interaktionsstärken zwischen den einzelnen Zelltypen beim Übergang der Wunde von Tag 1 (blau) zu Tag 14 (rot) DPW zeigt. Diese Zahl entspricht Schritt 6.16. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 24: Rangdiagramm, das die relativen Beiträge einzelner Signalwege zu Zell-Zell-Interaktionen zwischen Fibroblasten und anderen Zelltypen an Tag 1 vs. Tag 14 DPW zeigt. Diese Zahl entspricht Schritt 6.17. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 25: Blasendiagramme, die die relativen Beiträge einzelner Liganden-Rezeptor-Paare im Kollagen-Signalweg mit Fibroblasten als Quellzellen an Tag 1 vs. Tag 14 DPW zeigen. Diese Zahl entspricht Schritt 6.18. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 26: UMAP-Diagramme, die die Verteilung der Seurat-Cluster (links) und der Chargennummern (rechts) vor der Datenintegration zeigen. Diese Zahl entspricht Schritt 7.5. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Abbildung 27: UMAP-Diagramme, die die Verteilung der Seurat-Cluster (links) und der Chargennummern (rechts) nach der Datenintegration zeigen. Diese Zahl entspricht Schritt 7.9. Bitte klicken Sie hier, um diese Abbildung herunterzuladen.

Ergänzende Datei 1: JoVE_Rscript.R: Hauptskriptdatei des R-Codes, die alle Schritte und Erklärungen enthält, die für alle Teile des Protokolls beschrieben sind. Bitte klicken Sie hier, um diese Datei herunterzuladen.

Ergänzende Akte 2: JoVE_PhaseSpecificGenes.txt. Tabulatorgetrennte Textdatei, die die Listen der Gene enthält, die in Schritt 5.1 des Protokolls geladen wurden. Bitte klicken Sie hier, um diese Datei herunterzuladen.

Zusatzdatei 3: JoVE_Rscript_b3.R. Ergänzende R-Codeskriptdatei, die alle Schritte und Erklärungen enthält, die zum Analysieren von Batch #3 des Datasets für die Verwendung in Schritt 7.1 des Protokolls erforderlich sind. Bitte klicken Sie hier, um diese Datei herunterzuladen.

Ergänzende Tabelle 1: JoVE_DEGs_cellMarkers.xlsx. Excel-Datei, die die vollständige Ausgabe der in Schritt 3.10 des Protokolls verwendeten differentiell exprimierten Gene enthält. Bitte klicken Sie hier, um diese Tabelle herunterzuladen.

Ergänzende Tabelle 2: Top 5 der hochregulierten und exprimierten Gene für jeden Seurat-Cluster. Bitte klicken Sie hier, um diese Tabelle herunterzuladen.

Discussion

In diesem Protokoll wird RStudio verwendet, um vorgeschriebene Codezeilen auszuführen, die eine grundlegende Analyse eines komplexen Einzelzelldatensatzes mit Seurat ermöglichen. Es werden mehrere Methoden vorgestellt, die für die Wundheilungsforschung relevant sind, einschließlich der Installation der R-Kodierungsumgebung, des Herunterladens eines zuvor veröffentlichten Einzelzell-Wundheilungsdatensatzes, der Durchführung kritischer Qualitätskontrollschritte und standardmäßiger Einzelzellanalyse-Workflows, einschließlich Visualisierungen, Annotationen der wichtigsten Zelltypen, Zellsubtyp-Analysen und integrativen Analysen mit Seurat sowie der Durchführung von Zell-Zell-Interaktionsanalysen mit CellChat.

Die hier vorgestellten Methoden sind vereinfachte Vignetten typischer Arbeitsabläufe für die Einzelzellanalyse mit R und seinen beliebten Open-Source-Wissenschaftspaketen Seurat²¹ und CellChat²². Tatsächlich ist der Workflow nur ein Beispiel für die Art von Analyse, die man mit einem komplexen Einzelzell-Wundheilungsdatensatz durchführen kann. Die möglichen Modifikationen dieser Methode sind nahezu unendlich, wobei die einzige Einschränkung die spezielle wissenschaftliche Untersuchung des Benutzers ist. Zum Beispiel kann der Benutzer einige der Schlüsselparameter, wie z. B. Zelltypen und Zeitpunkte, entsprechend den Forschungsfragen ändern, die er diesem Datensatz stellen möchte. Die Autoren hoffen auch, dass sich der Benutzer wohl genug fühlt, um diesen Workflow an seinen eigenen Einzelzelldatensatz anzupassen. Bei der Verwendung dieses Workflows zur Analyse anderer Datensätze ist jedoch Vorsicht geboten, da jedes Experiment technische Probleme und Probleme bei der Probenvorbereitung auf die Daten selbst übertragen kann. Daher ist es unerlässlich, dass der Benutzer alle experimentellen Details liest und versteht, bevor er die Ergebnisse zuvor veröffentlichter und erneut analysierter Einzelzelldatensätze interpretiert. Es ist wichtig, sich daran zu erinnern, dass bioinformatische Werkzeuge eine leistungsfähige Methode zur Erforschung biologischer Prozesse und zur Hypothesengenerierung sind und dass alle kritischen biologischen Interpretationen der Ergebnisse in Folgeexperimenten validiert werden müssen.

Beachten Sie im gesamten Protokoll, dass größere Änderungen an bestimmten Bereichen des Workflows vorgenommen werden können, um andere Aufgaben auszuführen. Die Einzelheiten aller möglichen Kombinationen von Änderungen am Arbeitsablauf gehen jedoch über den Rahmen dieses Manuskripts hinaus. Zum Beispiel sind die Auflösung, die für die Zellclusterbildung verwendet wird, und die Dimensionen, die für die UMAP-Analyse verwendet werden, notwendigerweise subjektiv, und die hier vorgestellten Werkzeuge ermöglichen sowohl groß angelegte Analysen (wie hier für weit gefasste Hauptzelltypen gezeigt wurde) als auch sehr spezifische Analysen, die die Unterclusterung von Zellen in seltenere Subpopulationen innerhalb des größeren Datensatzes beinhalten könnten. Für weitere Informationen zu diesem Aspekt der Einzelzellanalysemethode und für Details zu allen anderen Parametern, die in der Einzelzellanalyse-Pipeline geändert werden können, verweisen die Autoren den Benutzer auf die Seurat Publikationen^21,26 und die Website (https://satijalab.org/seurat/), wo die Autoren dieses sich weiterentwickelnden Tools ausführliche Erklärungen, Vignetten und Tutorials bereitstellen.

In diesem Manuskript wurden einige der in der Einzelzell-Transkriptomik-Literatur am häufigsten zitierten und verwendeten Werkzeuge vorgestellt, nämlich Seurat²¹ und CellChat²² für Einzelzell- bzw. Zell-Zell-Interaktionsanalysen. Es gibt jedoch auch andere Tools, die ähnliche Funktionen auf leicht unterschiedliche Weise ausführen. Für die Analyse von Einzelzelldatensätzen gibt es Scran²⁷, Scater ²⁸ und das Python-basierte ScanPy²⁹, die verschiedene Methoden für die Datensatzintegration^{verwenden 25}. In diesem Protokoll wurde die manuelle Annotation von Zelltypen demonstriert, die bei der Interpretation von Cluster-Zellmarker-Anreicherungen auf dem Urteilsvermögen des Benutzers beruht, aber es gibt jetzt verschiedene Werkzeuge, die eine automatisierte Klassifizierung von Zelltypen ermöglichen, wie z.B. SingleR³⁰ und scGate³¹. Für die Analyse der Zell-Zell-Kommunikation wurde CellChat in diesem Protokoll demonstriert, aber es gibt auch andere Werkzeuge zur Schätzung der Zell-Zell-Kommunikation, einschließlich CellPhoneDB³², Cytotalk³³ und andere Liganden-Rezeptor-Datenbanken, die innerhalb des LIANA-Konsensus-Frameworks³⁴ (LIgand-receptor ANalysis framework) implementiert sind. Alle Bioinformatik-Tools sind einzigartig und haben ihre eigenen Besonderheiten und modifizierbaren Parameter. Daher ist es wichtig, dass der Benutzer die zugehörige Dokumentation jedes Tools sorgfältig liest, um seine Nuancen zu verstehen, bevor er die aus ihrer Verwendung generierten Ergebnisse interpretiert. Unabhängig davon, welche Bioinformatik-Tools man verwendet, ist es wichtig, sich daran zu erinnern, dass sich solche Tools ständig weiterentwickeln und dass verschiedene Versionen der Pakete unterschiedliche Ergebnisse liefern können.

In R ist die Syntax von entscheidender Bedeutung, und ein falsch platziertes Satzzeichen, Anführungszeichen, eine eckige Klammer oder sogar ein falsch geschriebener Buchstabe führen zu einem Fehler. Daher ist es für den Nutzer unerlässlich, beim Tippen von Code auf Details zu achten und besonders achtsam beim Kopieren von Codezeilen zu sein, um ihn an neue wissenschaftliche Fragestellungen und Datensätze anzupassen. Zur Fehlerbehebung der spezifischen Fehler, auf die man stoßen kann, empfehlen die Autoren, die Fehlermeldung einfach in die bevorzugte Websuchmaschine des Benutzers zu kopieren und einzufügen und die Ergebnisse aus Bioinformatikforen wie GitHub und Stack Overflow zu durchsuchen, da die am häufigsten auftretenden Fehler wahrscheinlich bereits von einem sachkundigen Power-User beantwortet wurden. In einigen Foren werden die erfolgreichsten Antworten von anderen Benutzern "upvotet", die die Lösung für das Problem am besten gefunden haben. Der Benutzer muss darauf achten, dass er Codezeilen, die er im Internet gefunden hat, nicht einfach kopiert und auf seinen eigenen Computer einfügt (insbesondere, wenn eine Lösung die Änderung der Systemeinstellungen außerhalb der Programmiersprache R erfordert), da die Möglichkeit besteht, dass solche Programme bösartig sein können. Eine aufregende aufstrebende Methode zur Fehlerbehebung von Codierungsfehlern ist die Verwendung der leistungsstarken generativen KI-Modelle für große Sprachen wie ChatGPT von OpenAI, Copilot von Microsoft oder Gemini von Google. Diese Modelle haben sich als besonders nützlich für das Software-Engineering im Allgemeinen und die Fehlerbehebung im Besonderen erwiesen. Zu diesem Zweck kann der Benutzer ganze Zeilen seines Codes kopieren und einfügen, nachdem er dem Chatbot eine einfache Aufforderung über die Absicht des Benutzers für den Code gegeben hat. Der übliche Vorbehalt besteht darin, dass diese Modelle nicht narrensicher sind und der Benutzer möglicherweise mehrere Eingabeaufforderungen ausprobieren muss, um eine Antwort zu erhalten, die zur Lösung des Problems geeignet ist.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Laptop oder Desktop-Computer	N/A	N/A	Unter Windows oder MacOS
R	N/A	Version 4.4.1	Kostenloser Download von https://cran.rstudio.com/
Rstudio	Posit Software, PBC	Version 2024.09.0	Kostenloser Download von https://posit.co/download/rstudio-desktop/
Büro Excel	Microsoft	Jede Version	Für die Analyse von Tabellendaten
Internetbrowser	N/A	N/A	Zum Navigieren zu Websites
R Pakete	Repository	Version
devtools	CRAN	2.4.5
readxl	CRAN	1.4.3
OpenXLSX	CRAN	4.2.7.1
aufgeräumt	CRAN	2.0.0
scAnpassen	CRAN	2.1.2
BiocManager	Bioleiter	1.30.25
NMF	Bioleiter	0.28
KomplexeHeatmap	Bioleiter	2.20.0
BiocNeighbors (BiocNeighbors)	Bioleiter	1.22.0
SingleCellExperiment	Bioleiter	1.26.0
Zirkulieren	Bioleiter	0.4.16
edgeR	Bioleiter	4.2.1
scDblFinder	Bioleiter	1.18.0
Seurat	CRAN	5.1.0
CellChat (Englisch)	Github (Englisch)	2.1.2