Utilisation de R, Seurat et CellChat pour analyser un ensemble de données transcriptomiques unicellulaires de la cicatrisation des plaies cutanées de souris

La cicatrisation des plaies est l’un des processus les plus complexes de la biologie des mammifères, et elle implique un spectre de trois phases de cicatrisation : inflammatoire, proliférative et résolution ^1,2. Ces phases de cicatrisation classifient largement les actions coordonnées de dizaines de types de cellules et de centaines de leurs produits moléculaires à travers l’espace et le temps de la réparation des plaies³. Plusieurs décennies d’études histologiques et moléculaires basées sur l’échantillonnage de tissus de la plaie tout au long de la cicatrisation ont permis d’élucider les modèles cellulaires globaux de la réparation des tissus³, en particulier dans des modèles murins reproductibles de cicatrisation excisionnelle des plaies cutanées ^4,5,6. Ce n’est qu’au cours des deux dernières décennies qu’il est devenu possible d’apprécier plus pleinement la complexité de la cicatrisation des plaies, à commencer par l’avènement des analyses transcriptomiques à haut débit des plaies à l’échelle des tissus en vrac ^7,8,9 et des cellules 10,11,12,13,14. Plus récemment, plusieurs études ont profilé transcriptionnellement les plaies cutanées au niveau de la cellule unique, identifiant de nouveaux sous-types de cellules de plaie et montrant comment elles peuvent interagir les unes avec les autres pendant la cicatrisation 15,16,17,18,19,20. Hu et al. ont utilisé une approche innovante de séquençage spatial de l’ARN unicellulaire pour profiler les plaies cutanées tout au long du cours temporel de la cicatrisation à plusieurs distances radiales du centre de la plaie, ce qui a révélé de nouveaux « mouvements » intercellulaires et moléculaires à travers l’espace et le temps²⁰. De telles études dévoilent la complexité de la cicatrisation des plaies avec des détails sans précédent, et elles commencent à brosser un tableau d’une énorme hétérogénéité cellulaire et moléculaire.

Les progrès majeurs et récents dans les méthodes d’analyse bio-informatique permettent de donner un sens biologique aux ensembles de données multi-omiques complexes générés dans le domaine de la recherche sur la cicatrisation des plaies. Les progiciels d’analyse de cellules uniques comme Seurat fournissent des outils pour l’analyse robuste et l’intégration d’ensembles de données, y compris la classification des types de cellules dans des tissus complexes tels que les plaies²¹. Pour l’interprétation en aval des données unicellulaires, des outils tels que CellChat sont utilisés pour identifier les programmes d’interaction cellule-cellule présumés qui peuvent expliquer comment les cellules se coordonnent pour réparer les plaies²². Bien que ces outils soient bien documentés et bien cités, ils doivent être exécutés dans un environnement de codage informatique tel que R, un langage de programmation statistique et graphique qui est le plus couramment utilisé dans les domaines de la bioinformatique de la génomique et de la transcriptomique. Alors que les biologistes et les cliniciens dans le domaine de la cicatrisation des plaies utilisent de plus en plus d’approches unicellulaires pour étudier la réparation des tissus, peu d’entre eux ont la formation en bioinformatique nécessaire pour utiliser des outils tels que Seurat et CellChat directement dans leurs propres laboratoires. Un tel obstacle à l’utilisation de ces outils bioinformatiques empêche non seulement les scientifiques d’analyser plus en profondeur leurs propres ensembles de données sans l’aide de bioinformaticiens, mais il empêche également les scientifiques de réanalyser de manière fiable la richesse des données unicellulaires déjà publiées par d’autres groupes.

Par conséquent, un flux de travail étape par étape est présenté ici pour permettre aux scientifiques sans formation en bioinformatique d’analyser un ensemble de données de cicatrisation de plaies unicellulaires²⁰ précédemment publié et accessible au public. Le protocole utilise l’environnement de codage graphique R répandu et gratuit appelé RStudio, et il montre comment naviguer dans cet environnement pour exécuter des lignes de code prescrites qui permettent une analyse de base d’un ensemble de données complexe à cellule unique à l’aide de Seurat et CellChat. Dans ce protocole, les sept principales méthodes pertinentes pour la recherche sur la cicatrisation des plaies sont présentées, notamment : 1) l’installation de l’environnement de codage, 2) le téléchargement de l’ensemble de données et les étapes critiques de contrôle de la qualité, 3) les flux de travail d’analyse de cellules uniques, y compris les visualisations et les annotations de type de cellule, 4) les analyses de sous-types de cellules, 5) les analyses de notation de module, 6) les analyses d’interaction cellule-cellule et 7) les analyses intégratives de plusieurs ensembles de données. Dans chaque méthode, le code réel est fourni à l’utilisateur pour qu’il l’exécute côte à côte avec le protocole, et les résultats graphiques réels de chaque ligne de code sont affichés pour guider l’utilisateur tout au long du flux de travail. L’objectif principal de cette introduction guidée et visuelle à RStudio et à un flux de travail fondamental d’analyse de cellules uniques est de permettre à davantage de scientifiques spécialisés dans la cicatrisation des plaies d’utiliser directement ces outils puissants afin de permettre des progrès plus rapides dans le domaine de la recherche.

NOTE: Dans les flux de travail suivants, détaillant sept méthodes bioinformatiques, toutes les étapes des protocoles sont accompagnées de leurs blocs de code respectifs qui doivent être exécutés directement sur l’interface RStudio de l’utilisateur dans l’ordre dans lequel ils sont répertoriés. Pour rendre ce protocole aussi convivial que possible, un fichier de script R est inclus (Supplementary File 1 : JoVE_Rscript.R), qui peut être chargé directement dans la session RStudio de l’utilisateur, de sorte que chaque ligne de code peut être simplement exécutée. Cela évite à l’utilisateur d’avoir à taper ou copier-coller le code du document de protocole, ce qui pourrait introduire des erreurs. Toutes les instructions du protocole sont également incluses dans le fichier de script R sous forme de commentaires, indiqués par un symbole hashtag '#' au début de chaque ligne de commentaire.

1. Installation de R, RStudio et des packages R requis pour le flux de travail d’analyse de cellules uniques

1. Téléchargez et installez R (version 4.4.1) sur l’ordinateur. Utilisez le lien qui correspond au système d’exploitation de l’ordinateur.
   1. Si vous utilisez un ordinateur fonctionnant sous Microsoft Windows, utilisez ce lien : https://cran.rstudio.com/bin/windows/base/
   2. Si vous utilisez un ordinateur fonctionnant sous MacOS, utilisez ce lien : https://cran.rstudio.com/bin/macosx/
2. Installez la dernière version de RStudio sur l’ordinateur. Cliquez sur le lien suivant et suivez les instructions :
   https://www.rstudio.com/products/rstudio/download/#download
3. Installez les Rtools (version 4.4), qui permettront à R de compiler certains paquets. Cliquez sur le lien suivant et suivez les instructions :
   1. Si vous utilisez Windows, utilisez ce lien : https://cran.rstudio.com/bin/windows/Rtools/
   2. Si vous utilisez MacOS, utilisez ce lien :
      https://mac.r-project.org/tools/
4. Définir le répertoire de travail local ; Il s’agit du dossier de l’ordinateur à partir duquel tous les fichiers seront chargés et enregistrés. Définissez le répertoire de travail en sélectionnant Session dans la barre de menu RStudio et en cliquant sur Définir le répertoire de travail > Choisir le répertoire et en sélectionnant le dossier souhaité.
   1. Si vous utilisez un ordinateur Windows, utilisez la commande suivante pour définir le répertoire de travail. Remplacez [Directory] dans la ligne de code suivante par la structure de répertoire réelle. N’oubliez pas que le délimiteur de répertoire dans R est le caractère « / »
      setwd("C:/[Directory]")
   2. Si vous utilisez un ordinateur MacOS, la commande suivante définira également le répertoire de travail. Remplacez [Directory] dans la ligne de code suivante par la structure de répertoire réelle. N’oubliez pas que le délimiteur de répertoire dans R est le caractère « / »
      setwd("~/[Directory]")
   3. À tout moment au cours d’une session R, vérifiez le répertoire de travail à l’aide de la ligne de code suivante :
      getwd()
   4. Dans RStudio, explorez visuellement la structure des répertoires de travail, y compris tous les fichiers et dossiers qu’il contient, dans la fenêtre de droite de l’onglet Fichiers . Pour naviguer dans l’explorateur de fichiers de RStudio jusqu’au répertoire de travail, cliquez sur l’icône d’engrenage -> Aller au répertoire de travail.
5. Installez les packages suivants à partir du référentiel de packages R CRAN, qui sont des dépendances nécessaires pour le protocole. Pour installer ces packages, exécutez les commandes suivantes.
   install.packages("devtools")
install.packages("readxl")
install.packages("openxlsx")
install.packages("tidyverse")
install.packages("scCustomize")
   NOTE: Lors de l’installation de packages R, il est normal que différentes fenêtres apparaissent et disparaissent. Si une fenêtre s’affiche vous demandant de compiler un paquet, cliquez sur OUI. Si une fenêtre s’affiche vous demandant de redémarrer R avant d’installer le package, cliquez sur NON.
6. Installez les packages suivants à partir du référentiel de packages R organisé Bioconductor
   (https://bioconductor.org/), qui sont des dépendances nécessaires pour le protocole. Pour installer ces packages, exécutez les commandes suivantes.
   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. Installez les packages suivants, qui sont requis pour le flux de travail décrit dans ce manuscrit.
   install.packages("Seurat")
devtools::install_github("jinworks/CellChat")
8. Chargez dans chaque paquet pour confirmer que les installations ont réussi. Si l’un des paquets génère l’erreur « paquet introuvable », réinstallez-le en utilisant le code approprié ci-dessus.
   library(readxl)
library(openxlsx)
library(tidyverse)
library(scCustomize)
library(edgeR)
library(scDblFinder)
library(Seurat)
library(CellChat)

2. Chargement d’un ensemble de données de cicatrisation de plaies unicellulaire et exécution des étapes de contrôle qualité

NOTE: Pour ce flux de travail bioinformatique, une réanalyse d’une expérience spatio-temporelle de cicatrisation de plaies cutanées unicellulaires précédemment publiée est menée²⁰. Les fichiers de l’ensemble de données sont stockés dans le référentiel NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/).

1. Accédez aux fichiers de l’ensemble de données à partir de GEO en utilisant le numéro d’acquisition GSE204777. Utilisez le lien suivant et examinez le plan expérimental de l’étude : https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204777
2. Sur la page du référentiel GEO, il y a cinq lots individuels de données obtenues à partir de cinq voies de séquençage. Cliquez sur le premier jeu de données, intitulé GSM6190913. Ce qui suit est un lien direct vers l’exemple : https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM6190913
3. Faites défiler vers le bas de la page et téléchargez les trois fichiers suivants, en utilisant les liens ftp ou html . Dans l’explorateur de fichiers de l’ordinateur, déplacez ces trois fichiers dans un répertoire appelé b1. Assurez-vous que le dossier b1 se trouve dans le répertoire de travail défini à l’étape 1.4.
   Nom du fichier : GSM6190913_b1_barcodes.tsv.gz / Taille du fichier : 18.5 Mb
   Nom du fichier : GSM6190913_b1_features.tsv.gz / Taille du fichier : 254.1 Kb
   Nom du fichier : GSM6190913_b1_matrix.mtx.gz / Taille du fichier : 151.2 Mb
4. Obtenez les informations de répertoire pour les fichiers de séquençage à cellule unique téléchargés à l’étape 2.3.
   b1_data_dir <- file.path(getwd(), "b1")
5. Chargez les fichiers de séquençage à cellule unique. Le paramètre gene.column spécifie le nom de gène/caractéristique utilisé. Dans ce cas, utilisez gene.column = 2 pour les symboles de gènes (gene.column = 1 pour les noms de gènes Ensembl).
   b1 <- Read10X_GEO(data_dir = b1_data_dir, gene.column = 2)
   NOTE: La plupart des jeux de données à cellule unique n’ont pas de données multiplexées supplémentaires, de sorte que les fichiers 10x générés à l’aide de cette étape n’auraient pas plusieurs couches. Pour le jeu de données multiplexé actuel, passez à l’étape 2.6. Pour les jeux de données sans multiplexage des données, passez à l’étape 2.7.
6. Démultiplexez l’ensemble de données à cellule unique à l’aide de codes-barres spatio-temporels.
   1. Pour l’ensemble de données de travail, séparez les données d’expression génique et de HTO (multiplexage).
      dataset_rna_counts <- b1$GSM6190913_b1_$`Gene Expression`
​dataset_barcodes <- b1$GSM6190913_b1_$`Antibody Capture`
   2. Créez un objet Seurat à l’aide des données d’expression génique, tout en filtrant immédiatement les gènes exprimés dans moins de 5 cellules et les cellules avec moins de 200 gènes détectés.
      dataset <- CreateSeuratObject(counts = dataset_rna_counts, min.cells = 5, min.features = 200)
   3. Créez un ensemble de données d’expression génique sous forme de couche et générez une liste de cellules et de codes-barres communs aux deux tests.
      dataset_rna_counts2 <- LayerData(dataset, data = "RNA")
​dataset_joint.barcodes <- intersect(colnames(dataset_rna_counts2), colnames(dataset_barcodes))
   4. Sous-ensemble de l’expression génique et du nombre de HTO par des codes-barres de cellules conjointes.
      dataset_rna_counts2 <- dataset_rna_counts2[, dataset_joint.barcodes]
​dataset_barcodes2 <- as.matrix(dataset_barcodes[, dataset_joint.barcodes])
   5. Vérifiez que l’OCT a les noms de code-barres attendus.
      ​rownames(dataset_barcodes2)
   6. Créez un nouveau test pour stocker les informations de code-barres et ajoutez-le à l’objet Seurat précédemment créé.
      dataset_barcode_assay <- CreateAssayObject(counts = dataset_barcodes2)
​dataset[["barcodes"]] <- dataset_barcode_assay
   7. Vérifiez que l’objet contient désormais plusieurs analyses.
      DefaultAssay(dataset)
   8. Normalisez les données des codes-barres et effectuez le démultiplexage via la fonction HTODemux. Cette méthode est décrite en détail dans la vignette Seurat suivante : https://satijalab.org/seurat/articles/hashing_vignette
      ​dataset <- HTODemux(dataset, assay = "barcodes", positive.quantile = 0.99)
   9. Regroupez les cellules en fonction des résultats de la classification globale et supprimez les cellules sans classification par code-barres.
      Idents(dataset) <- "barcodes_classification.global"
​dataset <- subset(dataset, idents = "Negative", invert = TRUE)
   10. Regroupez les cellules en fonction du signal HTO maximal.
       Idents(dataset) <- "barcodes_maxID"
   11. Visualisez la distribution des gènes détectés dans les cellules à l’aide de leurs codes-barres multiplexés (figure supplémentaire 1).
       VlnPlot(dataset, features = "nFeature_RNA", pt.size = 0.1, log = TRUE)
   12. Renommez les codes-barres en fonction de leur temps de blessure (jours après la blessure) et de leur espace (2-8 mm) (tirés du manuscrit original) et attribuez-les à une nouvelle variable de métadonnées appelée time_space.
       Idents(dataset) <- "barcodes_maxID"
levels(dataset)
dataset <- RenameIdents(dataset,
"Barcode-1" ="D01_2mm",
"Barcode-2" ="D01_4mm",
"Barcode-3" ="D01_6mm",
"Barcode-4" ="D01_8mm",
"Barcode-5" ="D03_2mm",
"Barcode-6" ="D03_4mm",
"Barcode-7" ="D03_6mm",
"Barcode-8" ="D03_8mm",
"Barcode-9" ="D07_2mm",
"Barcode-10" ="D07_4mm",
"Barcode-11" ="D07_6mm",
"Barcode-12" ="D07_8mm",
"Barcode-13" ="D14_2mm",
"Barcode-14" ="D14_4mm",
"Barcode-15" ="D14_6mm",
"Barcode-16" ="D14_8mm",
"Barcode-17" ="UW")
levels(dataset)
dataset[["time_space"]] <- Idents(dataset)
7. Pour les jeux de données sans multiplexage des données : créez un objet Seurat, tout en filtrant immédiatement les gènes exprimés dans moins de 5 cellules et les cellules avec moins de 200 gènes détectés.
   dataset <- CreateSeuratObject(counts = b1, min.cells = 5, min.features = 200)
8. Passez à l’analyse du test d’expression génique de l’ensemble de données.
   DefaultAssay(dataset) <- "RNA"
Idents(dataset) <- "data"
9. Comme étape importante du contrôle de la qualité, calculez le pourcentage de gènes mitochondriaux dans chaque cellule et attribuez-le en tant que variable de métadonnées. Cette méthode est décrite en détail dans la vignette Seurat suivante : https://satijalab.org/seurat/articles/pbmc3k_tutorial#qc-and-selecting-cells-for-further-analysis
   dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^mt-")
10. Visualisez la distribution des gènes détectés, le nombre d’ARN et le pourcentage mitochondrial dans toutes les cellules (Figure supplémentaire 2).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
11. Il existe un certain nombre de cellules avec un contenu mitochondrial important, ce qui est corrélé à un faible nombre d’ARN ; Ce sont des cellules mortes ou mourantes. Supprimez ces cellules de mauvaise qualité de l’ensemble de données en utilisant une limite équitable. Dans ce cas, utilisez les valeurs de l’ensemble de données d’origine décrites dans l’étude²⁰ précédemment publiée, dans laquelle des cellules contenant plus de 25 % de gènes mitochondriaux sont retirées.
    dataset <- subset(dataset, subset = nFeature_RNA > 200 & percent.mt < 25)
12. Visualisez la distribution des gènes détectés, le nombre d’ARN et le pourcentage mitochondrial dans toutes les cellules après l’élimination des cellules de faible qualité (Figure supplémentaire 3).
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
13. Comme étape supplémentaire importante du contrôle de la qualité, détectez les doublets probables dans l’ensemble de données. Ce sont des cellules qui ont été jointes lors du séquençage des gouttelettes et qui entraîneront ainsi des expressions géniques qui ne sont pas au niveau de la cellule unique. Pour contrer ce problème, plusieurs outils ont été développés. Il est important de noter que chaque outil fait certaines hypothèses sur les données d’une seule cellule, il est donc important que l’utilisateur lise toute la documentation pertinente avant d’utiliser un seul outil. Dans ce flux de travail, implémentez une méthode appelée scDblFinder²³. Veuillez noter que cette méthode utilise un algorithme qui impose un taux de doublet attendu fixe. Utilisez les commandes suivantes pour exécuter le pipeline scDblFinder.
    DefaultAssay(dataset) <- "RNA"
sce <- scDblFinder(GetAssayData(dataset, assay="RNA", slot="counts"), samples=Idents(dataset))
    NOTE: Les ensembles de données multiplexes à cellule unique tels que celui-ci peuvent également être filtrés pour les doublets en supprimant les cellules qui expriment plusieurs codes à barres. Ce flux de travail n’a pas démontré cette méthode, car la plupart des jeux de données à cellule unique ne disposent pas de cette caractéristique unique. Au lieu de cela, un pipeline plus généralisable pour le dépistage des doublets est montré à l’aide de la méthode scDblFinder.
14. Attribuez le score du doublet à une nouvelle variable de métadonnées.
    dataset$scDblFinder.score <- sce$scDblFinder.score
15. Visualisez la distribution du score du doublet dans toutes les cellules (Figure supplémentaire 4).
    VlnPlot(dataset, features = "scDblFinder.score", raster=FALSE, pt.size=0.5)
16. Supprimez les cellules au-dessus du seuil de score double de 0,25. Ce seuil a été choisi sur la base du graphique de violon généré ci-dessus, qui a montré que la plupart des cellules de l’ensemble de données pouvaient être attribuées à des doubles scores très élevés ou très faibles, et 0,25 est un seuil raisonnable pour cet ensemble de données qui supprimerait la grande majorité des doublets probables sans supprimer de nombreux doublets improbables.
    dataset <- subset(dataset, scDblFinder.score < 0.25)
17. Enregistrez l’objet Seurat du jeu de données en tant que fichier RDS dans le répertoire de travail.
    saveRDS(dataset, "dataset_post_Method2.rds")

3. Analyse d’un ensemble de données de cicatrisation de plaies unicellulaires à l’aide de Seurat

REMARQUE : (Étape facultative) Si vous démarrez le flux de travail ici, chargez le fichier RDS enregistré en tant qu’objet Seurat.

dataset <- readRDS("dataset_post_Method2.rds")

1. Exécutez le flux de travail Seurat standard pour la normalisation et la mise à l’échelle des ensembles de données à cellule unique, la mise à l’échelle et l’analyse en composantes principales (ACP). Ce flux de travail standard est décrit dans les vignettes Seurat suivantes :
   Seurat - Tutoriel sur le clustering guidé : https://satijalab.org/seurat/articles/pbmc3k_tutorial
   Liste de commandes Seurat : https://satijalab.org/seurat/articles/essential_commands
   DefaultAssay(dataset) <- « ARN »
   dataset <- NormalizeData(object = dataset)
dataset <- FindVariableFeatures(object = dataset)
dataset <- ScaleData(object = dataset)
dataset <- RunPCA(object = dataset)
2. Visualisez l’ampleur de la variation de l’ensemble de données par rapport aux 50 premières dimensions de l’ACP (figure supplémentaire 5).
   ElbowPlot(dataset, reduction = "pca", ndims = 50)
   Une grande partie de la variation majeure se produit dans les 13 premières dimensions.
3. Effectuez le clustering de cellules du jeu de données à l’aide d’une plage de dimensions PCA définie de 1 à 13 et d’une résolution définie de 0,1.
   dataset <- FindNeighbors(dataset, verbose = TRUE, dims = 1:13)
dataset <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
   NOTE: Ce flux de travail se concentre sur les différences à grande échelle entre les types de cellules. Par conséquent, il utilise un paramètre assez conservateur pour la plage de dimensions PCA, dans lequel les 13 premières dimensions représentent la grande majorité de la variation au sein de l’ensemble de données. Pour la discrimination des cellules en sous-types plus petits et plus rares, l’utilisateur peut utiliser un nombre plus élevé de dimensions pour l’analyse en aval, car ces sous-types de cellules rares sont probablement à l’origine de niveaux plus faibles de variation de l’ensemble de données. Le paramètre de résolution est compris entre 0 et 1 et détermine l’ampleur de la séparation catégorielle imposée à l’ensemble de données. Le réglage de ce paramètre dépend de la question de recherche de l’utilisateur. Pour regrouper les cellules en de nombreux sous-types petits et rares, utilisez des résolutions plus élevées pour l’analyse en aval. Étant donné que ce flux de travail vise à explorer des différences plus larges entre les principaux types de cellules, il utilise une valeur de résolution assez faible de 0,1, qui devrait regrouper les cellules en groupes moins nombreux et plus grands. L’utilisation des paramètres mentionnés à l’étape 3.3 permet de discriminer 8 groupes de cellules uniques. Ces clusters sont automatiquement attribués à une variable de métadonnées appelée « seurat_clusters ».
4. Effectuez une analyse de réduction dimensionnelle UMAP et de recherche de voisins à l’aide des 13 premières dimensions PCA. Ajoutez le numéro de départ 123 pour assurer la reproductibilité de la projection de données résultante.
   dataset <- RunUMAP(dataset, verbose = TRUE, dims = 1:13, seed.use = 123)
   NOTE: L’algorithme UMAP est stochastique et introduit du caractère aléatoire dans la réduction dimensionnelle (voir « Stabilité et reproductibilité » dans https://cran.r-project.org/web/packages/umap/vignettes/umap.html). Bien que l’utilisation d’une graine cohérente fournisse un « niveau minimal de reproductibilité » à l’algorithme, le graphique résultant peut tout de même être légèrement différent de ce qui est montré dans les figures représentatives et les résultats en aval. Les tests ont révélé que les résultats seront particulièrement différents entre les ordinateurs fonctionnant sous Windows et ceux fonctionnant sous MacOS, probablement en raison de différentes implémentations de l’aléatoire dans ces systèmes d’exploitation.
5. Visualisez le regroupement des cellules sur un graphique UMAP (Figure 1).
   DimPlot(dataset, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
   NOTE: Le caractère aléatoire de l’algorithme UMAP peut générer des tracés légèrement différents, comme le montre la figure montrant des tracés UMAP alternatifs générés à l’aide du même code que ci-dessus sur des ordinateurs fonctionnant sous Windows et MacOS ; Notez les légères différences dans les formes des grappes. Il est donc impératif que l’utilisateur enregistre et horodate toutes les données et tous les graphiques au fur et à mesure qu’ils sont générés, et que toutes les analyses en aval sur les grappes soient effectuées avec soin et en tenant compte de la compréhension biologique, comme il est décrit ci-dessous pour l’annotation du type de cellule.
6. Étant donné que les étiquettes de l’expérience originale qui font référence à l’endroit et au moment de l’origine des cellules pendant la cicatrisation des plaies sont incluses, visualisez l’annotation temps/espace des cellules sur un graphique UMAP (Figure 2).
   DimPlot(dataset, group.by = "time_space", raster = FALSE, label = FALSE)
7. Générez une table d’association de groupes de cellules avec l’annotation temps/espace de l’enroulé.
   table(dataset$time_space, dataset$seurat_clusters)
8. Déterminez l’identité des principaux types de cellules dans le jeu de données. Pour ce faire, calculez les gènes exprimés différentiellement (DEG) entre tous les clusters. Récupérez les listes DEG pour les clusters, attribuez-les à une variable et enregistrez la sortie sous forme de fichier texte délimité dans le répertoire de travail.
   NOTE: Cette étape est gourmande en ressources CPU et peut prendre beaucoup de temps, en fonction du matériel de l’utilisateur.
   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. Téléchargez et ouvrez le fichier dataset_cluster_markers.txt inclus dans une feuille de calcul (par exemple, Excel) en copiant le contenu du fichier texte et en utilisant l’assistant d’importation de texte pour spécifier le délimiteur de virgule et l’identité des colonnes de nom de gène en tant que texte. Il est important d’indiquer que les noms de gènes sont 'Texte', sinon Excel convertira automatiquement certains noms de gènes en dates, par exemple Sept7 changeant en September 7..
10. Dans une feuille de calcul, filtrez les résultats en fonction des paramètres recommandés suivants :
    1. Classez la colonne avg_log2FC de la plus grande à la plus petite afin d’organiser toutes les lignes en fonction des variations décroissantes de log₂ fois (log2FC).
    2. Classez la colonne du cluster de la plus petite à la plus grande pour organiser toutes les lignes en fonction du nombre croissant de clusters Serat.
    3. Filtrez la colonne avg_log2FC pour les nombres supérieurs ou égaux à 2,5 afin d’afficher uniquement les gènes les plus différentiellement exprimés (DEG) dans le cluster indiqué par rapport aux autres clusters.
    4. Filtrez la colonne pct.1 pour les nombres supérieurs ou égaux à 0,4. Cette colonne fait référence au pourcentage de cellules de la grappe indiquée exprimant le gène indiqué (Grappe %) et la définition du seuil à 0,4 signifie que seuls les gènes exprimés dans au moins 40 % des cellules de la grappe indiquée sont représentés.
    5. Filtrez la colonne pct.2 pour les nombres inférieurs ou égaux à 0,2. Cette colonne fait référence au pourcentage de cellules NON dans le groupe indiqué exprimant le gène indiqué (Non-cluster %), et la définition du seuil à 0,2 signifie que seuls les gènes exprimés dans au plus 20 % des cellules NON dans le cluster indiqué sont indiqués.
    6. Filtrez la colonne p_val_adj pour les nombres inférieurs ou égaux à 0,01. Cette colonne fait référence à la valeur P ajustée ou au taux de fausses découvertes (FDR), indiquant l’intensité statistique de la DEG indiquée, et la définition du seuil à 0,01 signifie que seuls les gènes avec un FDR < 0,01 sont indiqués.
       REMARQUE : Le tableau supplémentaire 1 (JoVE_DEGs_cellMarkers.xlsx) contient la sortie complète des gènes différentiellement exprimés classés utilisés à l’étape 3.10 du protocole. Le tableau supplémentaire 2 montre les 5 principaux gènes de chaque groupe dans cette analyse, les gènes en gras étant utilisés pour les visualisations ultérieures.
11. Pour une annotation non biaisée de type de cellule des clusters, utilisez l’outil d’analyse d’enrichissement basé sur le Web EnrichR.
    Utilisez le lien : https://maayanlab.cloud/Enrichr/
12. Copiez les listes de DEG pour chaque cluster dans une fenêtre EnrichR distincte, puis cliquez sur Analyser. L’outil EnrichR parcourt la liste de gènes à travers des centaines de bases de données organisées et classe chaque terme enrichi dans chaque catégorie.
13. Pour les besoins de l’annotation de type de cellule, cliquez sur l’onglet Types de cellules ci-dessus et concentrez-vous sur les 5 principaux enrichissements dans les trois bases de données organisées par marqueurs de cellule sur le côté gauche (Figure 3) :
    CellMarker 2024 (http://bio-bigdata.hrbmu.edu.cn/CellMarker/)
    Tabula Sapiens (https://tabula-sapiens-portal.ds.czbiohub.org/)
    PanglaoDB Augmented (https://panglaodb.se/)
14. Sur la base des enrichissements des DEG dans ces bases de données, confirmer l’identité probable des 8 clusters. Notons qu’il y a deux amas (2, 6) qui s’enrichissent sous forme de fibroblastes ; Par conséquent, combinez ces clusters en annotations de type de cellule unique. Attribuez les identités de type de cellule sous forme d’étiquettes à une nouvelle variable de métadonnées appelée cell_types.
    Idents(dataset) <- "seurat_clusters"
dataset[["cell_types"]] <- Idents(dataset)
Idents(dataset) <- "cell_types"
dataset <- RenameIdents(dataset,
"0" = "Macrophage",
"1" = "Neutrophil",
"2" = "Fibroblast",
"3" = "Epithelial cell",
"4" = "Endothelial cell",
"5" = "T cell",
"6" = "Fibroblast",
"7" = "Smooth muscle cell"
)
levels(dataset)
dataset[["cell_types"]] <- Idents(dataset)
15. Visualisez les clusters de cellules renommés sous forme d’annotations sur un graphique UMAP (Figure 4).
    DimPlot(dataset, group.by = "cell_types", raster = FALSE, label = FALSE)
16. Visualisez la localisation des gènes marqueurs de cluster supérieurs (en gras) du tableau 1 sur une série de graphiques UMAP (Figure 5).
    DefaultAssay(dataset) <- "RNA"
FeaturePlot(dataset, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4)
17. Visualisez les DEG des marqueurs de cluster supérieurs sur un graphique à points, regroupés par les numéros de cluster d’origine (Figure supplémentaire 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. Visualisez les DEG des marqueurs de cluster supérieurs sur un graphique à points, regroupés par les types de cellules annotés (Figure 6).
    DotPlot(dataset, group.by = "cell_types", features = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "blue", mid = "grey", high = "red")
19. Pour effectuer des analyses de séries chronologiques, simplifiez d’abord le jeu de données en supprimant la composante spatiale. Pour les analyses temporelles, regroupez les annotations temps/espace de la blessure en jours post-blessure globaux (DPW) avec une nouvelle variable de métadonnées appelée « 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. Visualisez les nouveaux regroupements temps-course de la blessure sur un graphique UMAP (Figure supplémentaire 7).
    DimPlot(dataset, group.by = "DPW", raster = FALSE, label = FALSE)
21. Générez des tableaux indiquant le nombre de cellules de chaque type dans chaque DPW.
    ​table(dataset$DPW, dataset$cell_types)
    1. Étape facultative : Pour obtenir également des listes DEG pour les groupes de temps de blessure, attribuez-les à une variable et enregistrez la sortie sous forme de fichier texte délimité.
       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. Convertissez les numéros de cellules en proportions par catégorie pour mieux comprendre les changements relatifs dans la composition des types de cellules tout au long de la guérison. Visualisez la proportion de DPW dans chaque type de cellule (figure supplémentaire 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. Visualisez la proportion de types de cellules dans chaque DPW ( figure 7).
    pt2 <- table(dataset$cell_types, dataset$DPW)
pt2 <- as.data.frame(pt2)
pt2$Var1 <- as.character(pt2$Var1)
ggplot(pt2, aes(x = Var2, y = Freq, fill = Var1)) +
theme_bw(base_size = 15) +
geom_col(position = "fill", width = 0.5) +
xlab("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
24. Enregistrez l’objet Seurat du jeu de données en tant que fichier RDS dans le répertoire de travail.
    saveRDS(dataset, "dataset_post_Method3.rds")

4. Analyse des sous-types cellulaires à l’aide de Seurat

NOTE: La puissance de l’analyse unicellulaire permet de découvrir et d’analyser des sous-types rares au sein des principaux types de cellules analysés ci-dessus. Cet exemple se concentre sur les fibroblastes, qui se sont initialement regroupés en deux amas de Seurat avant d’être combinés en une seule catégorie. Cette partie du protocole se concentre spécifiquement sur les fibroblastes, à l’exclusion de tous les autres types de cellules, et explore leurs identités et leurs propriétés temporelles pendant la cicatrisation des plaies. Si nécessaire, chargez le fichier RDS enregistré en tant qu’objet Seurat.

dataset <- readRDS("dataset_post_Method3.rds")

1. Sous-ensemble de l’ensemble de données d’origine en fonction de l’identité de la cellule du fibroblaste.
   Idents(dataset) <- « cell_types »
   dataset_fibroblast <- subset(dataset, idents = "Fibroblast")
2. Effectuez l’ACP sur cet ensemble de données plus petit et visualisez l’ampleur de la variation de l’ensemble de données par rapport aux dimensions de l’ACP (figure supplémentaire 9).
   dataset_fibroblast <- RunPCA(dataset_fibroblast, verbose = TRUE)
ElbowPlot(dataset_fibroblast, reduction = "pca", ndims = 50)
   NOTE: Une grande partie de la variation majeure se produit dans les 9 premières dimensions.
3. Effectuez le clustering de cellules du jeu de données à l’aide d’une plage de dimensions PCA définie de 1 à 9 et d’une résolution définie de 0,1.
   dataset_fibroblast <- FindNeighbors(dataset_fibroblast, verbose = TRUE, dims = 1:9)
dataset_fibroblast <- FindClusters(dataset_fibroblast, verbose = TRUE, resolution = 0.1)
   NOTE: À l’aide de ces paramètres, l’algorithme discrimine 3 groupes de fibroblastes uniques. Ces clusters sont automatiquement attribués à une variable de métadonnées appelée seurat_clusters.
4. Effectuez une analyse de réduction dimensionnelle UMAP et de recherche de voisins à l’aide des 9 premières dimensions PCA. Ajoutez le numéro de départ 123 pour assurer la reproductibilité de la projection de données résultante.
   dataset_fibroblast <- RunUMAP(dataset_fibroblast, verbose = TRUE, dims = 1:9, seed.use = 123)
5. Visualisez le regroupement des cellules sur un graphique UMAP (Figure 8).
   DimPlot(dataset_fibroblast, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
6. Visualisez l’annotation temps-course de l’enroulement des cellules sur un graphique UMAP (Figure supplémentaire 10).
   DimPlot(dataset_fibroblast, group.by = "DPW", raster = FALSE, label = FALSE)
7. Obtenez les listes DEG pour les trois sous-types de fibroblastes et enregistrez-les dans un fichier texte dans le répertoire de travail.
   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. Suivez les étapes similaires à celles ci-dessus (étapes 3.9-3.10) dans une feuille de calcul pour filtrer les DEG dans les marqueurs de sous-type de fibroblaste supérieurs.
9. Définissez une liste de gènes personnalisée en tant que variable en copiant les 5 principaux marqueurs de sous-type de fibroblastes pour chacun des trois groupes à partir du fichier texte DEG.
   FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
10. Visualisez les gènes de la liste dans l’ensemble de données uniquement sur les fibroblastes en appelant la variable dans le paramètre features du dotplot (Figure 9).
    DotPlot(dataset_fibroblast, group.by="seurat_clusters", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
DotPlot(dataset_fibroblast, group.by="DPW", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
11. Visualisez les gènes de la liste dans l’ensemble de données unicellulaire d’origine en appelant la variable dans le paramètre features du dotplot (Figure supplémentaire 11).
    DotPlot(dataset, group.by="cell_types", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
12. Visualisez la proportion de sous-types de fibroblastes dans chaque catégorie de DPW (figure supplémentaire 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. Visualisez la proportion de cellules de fibroblastes DPW dans chaque catégorie de sous-type de fibroblaste (Figure supplémentaire 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. Enregistrez l’objet Seurat du jeu de données en tant que fichier RDS dans le répertoire de travail.
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

5. Exemple d’analyse de suivi via le scoring des modules

NOTE: Une méthode utile pour analyser des ensembles de données à cellule unique est appelée notation de module. Dans ce flux de travail, il est possible de définir une liste de gènes en fonction des connaissances antérieures, puis de calculer les scores des modules, ce qui permet d’identifier les enrichissements potentiels de la liste de gènes au sein de chaque cellule. Ces scores peuvent être moyennés à travers les annotations de cellule pour révéler des modèles potentiels d’enrichissement.

Ici, utilisez des listes de gènes d’une étude précédemment publiée², où des gènes spécifiques à la phase de cicatrisation des plaies ont été identifiés à l’aide d’échantillons séquencés par l’ARN en vrac de l’ensemble du continuum de guérison. Les listes de gènes ont été enregistrées dans un fichier texte délimité par des tabulations (Fichier supplémentaire 2 : JoVE_PhaseSpecificGenes.txt) qui peut maintenant être téléchargé dans le répertoire de travail et utilisé pour générer des listes de gènes qui identifient les trois principales phases de guérison.

1. Chargez des listes de gènes dans une variable en lisant le fichier TEXT.
   PhaseSpecificGenes <- read_delim("JoVE_PhaseSpecificGenes.txt", delim="\t", col_names = T)
2. Séparez les colonnes en variables de liste de gènes individuelles et remplacez les gènes par des noms de souris dont la première lettre est en majuscule.
   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)
   REMARQUE : (Facultatif) Si nécessaire, chargez le fichier RDS enregistré en tant qu’objet Serat.
   dataset <- readRDS("dataset_post_Method3.rds")
3. Utilisez les listes de gènes comme modules pour noter chaque cellule de l’ensemble de données en fonction des trois phases de guérison.
   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. Visualisez les scores agrégés des modules par catégorie de cellules, y compris le DPW et les principaux types de cellules (Figure 10).
   DotPlot(dataset, group.by="DPW", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
DotPlot(dataset, group.by="cell_types", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")

6. Exemple d’une analyse de suivi via CellChat

NOTE: Une autre méthode utile et bien citée pour l’analyse d’ensembles de données unicellulaires consiste à déduire les interactions cellule-cellule. Dans ce flux de travail, utilisez le package CellChat, qui déduit les communications cellule-cellule en analysant les interactions différentielles ligand-récepteur entre les groupes de cellules²². Récemment, les développeurs de CellChat ont publié un protocole détaillé étape par étape pour son utilisation généralisée²⁴, et il s’agit d’une excellente ressource pour les utilisateurs lorsqu’ils travaillent sur le flux de travail suivant et l’appliquent à leurs ensembles de données. À titre d’exemple, le flux de travail suivant compare les interactions de toutes les cellules principales des plaies à 1 ou 14 jours après la blessure. Chaque étape n’est pas décrite en détail, car toutes les étapes ont déjà été décrites dans la publication officielle de CellChat²⁴ ainsi que dans les tutoriels, liés ici :

Inférence et analyse de la communication cellule-cellule à l’aide de
CellChat : https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html

Analyse comparative de plusieurs ensembles de données à l’aide de CellChat :
https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html

Étape facultative : Si nécessaire, chargez le fichier RDS enregistré en tant qu’objet Seurat :

dataset <- readRDS("dataset_post_Method3.rds")

1. Sous-ensemble du jeu de données d’origine en deux jeux de données en fonction de l’annotation DPW.
   Idents(dataset) <- "DPW"
dataset_D1 <- subset(dataset, ident = "D1")
dataset_D14 <- subset(dataset, ident = "D14")
2. Définissez l’annotation par laquelle CellChat sera exécuté --- dans ce cas, utilisez les principaux types de cellules.
   Idents(dataset_D1) <- "cell_types"
Idents(dataset_D14) <- "cell_types"
3. Créez les objets CellChat et suivez le flux de travail CellChat typique. Voir les tutoriels liés ci-dessus comme références détaillées pour chaque étape.
   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. Visualisez les forces d’interaction entrantes et sortantes dans tous les principaux types de cellules à chaque point de temps de cicatrisation des plaies (Figure supplémentaire 14).
   netAnalysis_signalingRole_scatter(cellchat_D1)
netAnalysis_signalingRole_scatter(cellchat_D14)
   Les fibroblastes augmentent considérablement leurs interactions entre D1 et D14 DPW.
5. Affichez les listes de toutes les voies de communication significatives inférées entre les cellules.
   cellchat_D1@netP$pathways
cellchat_D14@netP$pathways
   La voie du collagène est l’une des voies importantes sur les DPW D1 et D14.
6. Concentrez-vous sur la voie de signalisation du collagène et son interaction avec les fibroblastes.
   pathways.show <- c("COLLAGEN")
7. Visualisez les interactions de la voie de signalisation du collagène entre les types de cellules à l’aide de diagrammes circulaires (figure supplémentaire 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. Visualisez les interactions de la voie de signalisation du collagène entre les types de cellules à l’aide de diagrammes de cordes (Figure supplémentaire 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. Visualisez les interactions de la voie de signalisation COLLAGEN avec les fibroblastes en tant que cellules sources (Figure supplémentaire 17).
   NOTE: Les types de cellules dans les objets cellchat sont répertoriés en tant qu’ID dans l’ordre dans lequel ils ont été attribués dans l’objet Seurat d’origine : 1 = Macrophage, 2 = Neutrophile, 3 = Fibroblaste, 4 = Cellule épithéliale, 5 = Cellule endothéliale, 6 = Cellule T, 7 = Cellule musculaire lisse.
   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. Visualisez les contributions de chaque paire ligand-récepteur dans la voie de signalisation COLLAGEN avec des fibroblastes comme cellules sources.
    1. À l’aide de graphiques à bulles (figure supplémentaire 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. À l’aide de diagrammes d’accords (figure supplémentaire 19) :
       par(mfrow=c(1,2))
strwidth <- function(x) {0.4}
netVisual_chord_gene(cellchat_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, show.legend= F)
netVisual_chord_gene(cellchat_D14, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, legend.pos.x = 60)
​par(mfrow=c(1,1))
11. Focus sur l’interaction ligand-récepteur Col1a1-Cd44 dans la voie de signalisation COLLAGEN
    ​LR.show <- "COL1A1_CD44"
12. Visualisez les interactions ligand-récepteur Col1a1-Cd44 entre les types de cellules à l’aide de diagrammes de cordes (Figure supplémentaire 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. Effectuez une analyse différentielle CellChat en générant un objet CellChat combiné.
    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. Visualisez le nombre total et les forces relatives des interactions cellule-cellule dans les points temporels de cicatrisation des plaies (Figure supplémentaire 21).
    gg1 <- compareInteractions(cellchat_D14_v_D1, show.legend = F)
gg2 <- compareInteractions(cellchat_D14_v_D1, show.legend = F, measure = "weight")
gg1 + gg2
15. Visualisez à l’aide d’un graphique circulaire les forces d’interaction différentielles entre chaque type de cellule lorsque la plaie passe du jour 1 au jour 14 (figure supplémentaire 22).
    netVisual_diffInteraction(cellchat_D14_v_D1, weight.scale = T, measure = "weight")
16. Visualisez à l’aide d’une carte thermique les forces d’interaction différentielles entre chaque type de cellule lorsque la plaie passe du jour 1 au jour 14 (figure supplémentaire 23).
    netVisual_heatmap(cellchat_D14_v_D1, measure = "weight")
17. Visualisez à l’aide d’un graphique de rang les contributions relatives des voies individuelles aux interactions cellule-cellule avec les fibroblastes comme cellules sources au jour 14 par rapport au jour 1 ( figure supplémentaire 24).
    rankNet(cellchat_D14_v_D1, mode = "comparison", measure = "weight", sources.use = 3, targets.use = NULL, stacked = T, do.stat = TRUE)
18. Visualisez à l’aide de graphiques à bulles les contributions relatives des paires ligand-récepteur individuelles dans la voie de signalisation du collagène avec les fibroblastes comme cellules sources au jour 14 par rapport au jour 1 (figure supplémentaire 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. Tout comme les objets Seurat, les objets CellChat peuvent être enregistrés et ouverts en tant que fichiers RDS.
    saveRDS(cellchat_D1, "cellchat_D1.rds")
saveRDS(cellchat_D14, "cellchat_D14.rds")
saveRDS(cellchat_D14_v_D1, "cellchat_D14_v_D1.rds")
20. Étape facultative : Les objets CellChat peuvent également être ouverts à partir de fichiers RDS.
    cellchat_D1 <- readRDS("cellchat_D1.rds")
cellchat_D14 <- readRDS("cellchat_D14.rds")
cellchat_D14_v_D1 <- readRDS("cellchat_D14_v_D1.rds")

7. Exemple d’analyse intégrative par combinaison de plusieurs ensembles de données unicellulaires

NOTE: Les jeux de données à cellule unique sont souvent séparés en plusieurs fichiers, car ils ont été séquencés en lots ou en groupes. Ce flux de travail montre comment intégrer deux des cinq lots de l’ensemble de données de cicatrisation ²⁰. Les méthodes actuelles d’intégration des jeux de données sont décrites par les vignettes Seurat suivantes :

Introduction à l’intégration du scRNA-seq :
https://satijalab.org/seurat/articles/integration_introduction

L’analyse intégrative dans Seurat v5 :

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

Note: Il existe de nombreuses méthodes d’intégration d’ensembles de données à cellule unique, chacune ayant ses propres forces et faiblesses. Pour plus de détails, voir le benchmark complet des méthodes d’intégration²⁵. Il est important que l’utilisateur lise toute la documentation pertinente avant de s’appuyer sur une méthode d’intégration.

1. Répétez toutes les étapes de la méthode 2 pour un autre lot de l’ensemble de données Hu et al.²⁰. Dans le protocole suivant, le lot #3 est utilisé. Le fichier de script R supplémentaire est inclus et peut être utilisé pour traiter le lot #3 (Fichier supplémentaire 3 : JoVE_Rscript_b3. N’oubliez pas de créer et d’utiliser une nouvelle variable pour le jeu de données --- dans le code ci-dessous, utilisez « dataset_b3 » pour le lot de données #3.
   1. Étape facultative : Si nécessaire, ouvrez les deux jeux de données en tant qu’objets Seurat à partir de leurs fichiers RDS enregistrés dans le répertoire de travail :
      dataset <- readRDS("dataset_post_Method2.rds")
dataset_b3 <- readRDS("dataset_b3_post_Method2.rds")
2. Attribuez une nouvelle variable à chaque jeu de données nommé « batch » pour étiqueter le jeu de données d’origine dans les analyses ultérieures.
   dataset@meta.data$batch <- "b1"
dataset_b3@meta.data$batch <- "b3"
3. Effectuez la fusion Seurat des deux jeux de données, en ajoutant des annotations d’ID de cellule basées sur des lots, puis exécutez le flux de travail Seurat standard pour le jeu de données fusionné, comme décrit dans la méthode 3.
   dataset_merged <- merge(x = dataset, y = c(dataset_b3), add.cell.ids = c("b1", "b3"), merge.data = TRUE)
DefaultAssay(dataset_merged) <- "RNA"
dataset_merged <- NormalizeData(dataset_merged)
dataset_merged <- FindVariableFeatures(dataset_merged)
dataset_merged <- ScaleData(dataset_merged)
dataset_merged <- RunPCA(dataset_merged)
ElbowPlot(dataset_merged, reduction = "pca", ndims = 50)
4. Effectuez une analyse de clustering et UMAP sur l’ensemble de données combiné avant l’intégration des données.
   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. Visualisez le tracé UMAP en fonction des numéros de cluster et de lot (Figure supplémentaire 26).
   DimPlot(dataset_merged, reduction = "umap.unintegrated", group.by = c("seurat_clusters", "batch"))
6. Affichez la répartition des numéros de cellules dans chaque cluster en fonction du numéro de lot.
   table(dataset_merged$batch, dataset_merged$seurat_clusters)
   NOTE: D’après le graphique UMAP et le tableau, il ne semble pas y avoir d’effets de lot significatifs pour ces deux ensembles de données. Les preuves d’effets de lot se manifesteraient par des écarts inattendus dans la distribution des grappes entre les deux ensembles de données, ce qui pourrait signifier qu’il existe des différences techniques potentielles entre les ensembles de données qui l’emportent sur les similitudes biologiques réelles.
7. Effectuez l’intégration des données Seurat à l’aide de la méthode RPCA. Pour plus d’informations à ce sujet et sur d’autres méthodes d’intégration de données, veuillez lire la vignette Seurat liée ci-dessus.
   dataset_merged <- IntegrateLayers(
object = dataset_merged, method = RPCAIntegration,
orig.reduction = "pca", new.reduction = "integrated.rpca",
verbose = TRUE
)
8. Effectuez une analyse de clustering et UMAP sur le jeu de données combiné après l’intégration des données.
   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. Visualisez le tracé UMAP en fonction des numéros de cluster et de lot après l’intégration (Figure supplémentaire 27).
   DimPlot(dataset_merged, reduction = "umap.rpca", group.by = c("seurat_clusters","batch"))
10. Afficher la distribution des numéros de cellules dans chaque cluster en fonction du numéro de lot après l’intégration.
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    À partir du graphique UMAP et de la table des données intégrées, il y a maintenant un excellent chevauchement entre les deux lots sur différents clusters. Il est intéressant de noter que l’intégration des données a permis d’identifier un cluster supplémentaire en utilisant les mêmes paramètres de clustering.
11. Après l’intégration du jeu de données et avant les analyses en aval, les couches du jeu de données fusionné doivent être jointes.
    dataset_merged <- JoinLayers(dataset_merged)
12. Enregistrez l’objet Seurat du jeu de données en tant que fichier RDS dans le répertoire de travail.
    saveRDS(dataset_merged, "merged_dataset_post_Method7.rds")

À partir de la méthode #2, le protocole parcourt les étapes de chargement et d’exécution des étapes de contrôle de la qualité sur un ensemble de données de cicatrisation de plaies unicellulaires. Après avoir créé l’objet Seurat (étape 2.6.2), une série d’étapes fusionne les deux tests au sein de l’ensemble de données (ARN et protéine ; étapes 2.6.3-2.6.7) et effectue le décomplexage du test de protéine selon des codes-barres spatio-temporels (étapes 2.6.8-2.6.9). La fonction de décomplexification attribue plusieurs étiquettes de métadonnées à chaque cellule du jeu de données, dont « barcodes_maxID », qui identifie le code-barres spatio-temporel le plus probable de chaque cellule (étape 2.6.10). À l’étape 2.6.11, la fonction de graphique en violon est exécutée pour visualiser la distribution des gènes détectés dans les cellules en fonction de leurs codes-barres multiplexés. Le résultat représentatif de cette étape (figure supplémentaire 1) montre qu’il existe une distribution assez uniforme des gènes détectés pour chaque code-barres, ce qui est important pour l’intégrité de l’ensemble de données et l’analyse en aval des points temporels de cicatrisation des plaies. Après avoir attribué l’étiquette appropriée aux codes-barres des protéines (étape 2.6.12), le protocole montre ensuite comment effectuer des étapes de contrôle de la qualité sur le test d’ARN de l’ensemble de données, en commençant par calculer le pourcentage de gènes mitochondriaux dans chaque cellule (étape 2.9). À l’étape 2.10, la fonction de nuage de points des caractéristiques est exécutée pour visualiser la distribution des gènes détectés, le nombre d’ARN et le pourcentage mitochondrial dans toutes les cellules. Les résultats représentatifs de cette étape (figure supplémentaire 2) montrent qu’il existe un certain nombre de cellules avec un contenu mitochondrial important, ce qui est corrélé à un faible nombre d’ARN et identifie les cellules mortes ou mourantes. Après avoir retiré les cellules avec un faible nombre d’ARN et un contenu mitochondrial important (étape 2.11), à l’étape 2.12, une autre fonction de nuage de points de caractéristiques est exécutée sur l’ensemble de données du sous-ensemble, et le résultat représentatif de cette étape (figure supplémentaire 3) montre que la distribution des gènes détectés et le pourcentage d’ARN mitochondrial par cellule sont maintenant plus normaux, ouvrant la voie à des analyses robustes en aval. Ensuite, le protocole décrit l’utilisation de la fonction scDblFinder pour identifier les doublets probables dans l’ensemble de données et attribue une nouvelle métadonnée appelée « scDblFinder.score » à chaque cellule (étapes 2.13-2.14). À l’étape 2.15, la fonction de graphique du violon est exécutée pour visualiser la distribution des scores des doublets dans l’ensemble de données, et le résultat représentatif de cette étape (figure supplémentaire 4) montre qu’il y a un certain nombre de cellules avec des scores de doublets relativement élevés, et que 0,25 semble être un seuil naturel au-dessus duquel il y a une population de doublets probables. Par conséquent, les étapes suivantes utilisent ce paramètre pour sous-définir le jeu de données en cellules inférieures à la limite (étape 2.16), complétant ainsi les étapes de contrôle qualité pour ce jeu de données à cellule unique.

À partir de la méthode #3, le protocole décrit les étapes d’analyse de l’ensemble de données de cicatrisation des plaies unicellulaires contrôlées par la qualité à l’aide du package et du flux de travail Seurat. Après la normalisation et la mise à l’échelle des données d’ARN, une analyse PCA est effectuée (étape 3.1). À l’étape 3.2, la fonction de tracé coudé est utilisée pour visualiser la quantité de variation de l’ensemble de données par rapport aux 50 premières dimensions de l’ACP, et le résultat représentatif de cette étape (Figure supplémentaire 5) montre qu’une grande partie de la variation majeure se produit dans les 13 premières dimensions identifiées par la courbure dans le graphique. Le protocole montre ensuite comment trouver des voisins et effectuer le clustering de cellules (étape 3.3) et la réduction dimensionnelle UMAP (étape 3.4) de l’ensemble de données en utilisant les 13 premières dimensions PCA et un paramètre de résolution de clustering relativement faible de 0,1, qui ont tous deux été choisis afin d’identifier les principaux types de cellules les plus généralisables dans les plaies. À l’étape 3.5, la fonction de tracé dimensionnel est exécutée pour visualiser le regroupement des cellules sur un graphique UMAP, et le résultat représentatif de cette étape (Graphique 1) montre que toutes les cellules de l’ensemble de données sont regroupées autour de 8 grands groupes de clusters Seurat codés par couleur, avec des tracés UMAP légèrement différents obtenus à partir d’un ordinateur fonctionnant sous Windows (à gauche) et MacOS (à droite). À l’étape 3.6, une autre fonction de tracé dimensionnel est exécutée pour visualiser l’annotation temps/espace de l’enroulement des cellules, et le résultat représentatif de cette étape (Graphique 2) montre que toutes les cellules de l’ensemble de données sont réparties en fonction de leur origine spatio-temporelle, sans regroupement apparent en fonction de l’annotation temps/espace. Le protocole décrit ensuite comment obtenir des listes de gènes exprimés différentiellement et les enregistrer dans un fichier texte (étape 3.8), ouvrir la table de données dans une feuille de calcul et effectuer diverses étapes de filtrage afin d’obtenir les marqueurs de cluster les mieux classés pour chaque cluster de cellules (étapes 3.9-3.10.6). Le résultat représentatif de ces mesures (Tableau supplémentaire 1) est le fichier final de la feuille de calcul contenant la sortie complète des gènes exprimés différentiellement classés, tandis qu’un autre résultat représentatif (Tableau supplémentaire 2) est un tableau simplifié montrant les 5 principaux gènes régulés et exprimés à la hausse pour chaque cluster de Seurat. Le protocole décrit ensuite comment utiliser un outil d’analyse d’enrichissement fonctionnel basé sur le Web appelé EnrichR pour identifier les types de cellules présumés en fonction des gènes marqueurs de cluster supérieurs (étapes 3.11-3.12), et les résultats représentatifs de ces étapes (Graphique 3) sont des captures d’écran recadrées des sorties EnrichR montrant les types de cellules enrichis les plus élevés pour chacun des huit groupes de cellules. Le protocole attribue ensuite une nouvelle étiquette de métadonnées appelée « cell_types » à toutes les cellules des clusters Seurat respectifs en fonction de leurs annotations de type de cellule les plus enrichies (étape 3.14). À l’étape 3.15, la fonction de tracé dimensionnel est exécutée pour visualiser les grappes de cellules renommées en tant qu’annotations de type de cellule sur un tracé UMAP, et les résultats représentatifs de cette étape (Graphique 4) a montré que toutes les cellules de l’ensemble de données étaient regroupées autour des principaux types de cellules codés par couleur. À l’étape 3.16, la fonction de diagramme de caractéristiques a été utilisée pour visualiser la localisation des gènes marqueurs de grappe supérieurs (à partir du tableau supplémentaire 2) sur une série de graphiques UMAP, et les résultats représentatifs (Graphique 5) sont une grille de graphiques UMAP montrant la forte expression des gènes marqueurs cellulaires supérieurs dans leurs principaux groupes de types cellulaires respectifs. Dans les étapes 3.17 et 3.18, la fonction de diagramme à points a été exécutée pour visualiser les niveaux d’expression relatifs des gènes marqueurs de cluster supérieurs dans les cellules, d’abord regroupés par leur numéro de cluster Seurat d’origine (étape 3.17) et ensuite regroupés par des étiquettes de type de cellule annotées (étape 3.18). Les résultats représentatifs de ces étapes ont confirmé le niveau élevé d’expression des gènes marqueurs cellulaires supérieurs uniquement dans leurs grappes Seurat respectives (Figure supplémentaire 6) et seulement dans leurs principaux types cellulaires respectifs (Graphique 6). L’étape suivante du protocole simplifie les étiquettes spatio-temporelles originales basées sur les protéines en annotations strictement temporelles, qui identifient les cellules en fonction des jours post-blessure (DPW) d’où elles proviennent. À l’étape 3.20, la fonction de tracé dimensionnel est exécutée pour visualiser les cellules sous forme d’annotations DPW sur un graphique UMAP, et les résultats représentatifs de cette étape (Figure supplémentaire 7) a montré la localisation des annotations temporelles de l’évolution de la plaie dans l’ensemble de données de cicatrisation des plaies unicellulaires. Comme prévu, les annotations du jour 1 (D1) dominaient les amas de neutrophiles et de macrophages, tandis que les points temporels de cicatrisation des plaies ultérieurs étaient plus représentés dans d’autres types de cellules. Les étapes suivantes du protocole ont utilisé des diagrammes à barres empilées pour visualiser d’abord les proportions de DPW dans différents types de cellules (étape 3.22), puis pour visualiser les proportions de types de cellules à différents points temporels (3.23). Les résultats représentatifs de ces étapes sont des diagrammes de proportion montrant le nombre relatif de cellules DPW dans chaque grande catégorie de type de cellule (Figure supplémentaire 8) et le nombre relatif de principaux types de cellules dans chaque catégorie de MTP (Graphique 7). Ces résultats ont confirmé la cascade cellulaire connue de la cicatrisation des plaies cutanées, dans laquelle les cellules immunitaires (neutrophiles et macrophages) dominent les premiers points temporels de la phase inflammatoire, et les autres types de cellules (cellules épithéliales et cellules endothéliales) commencent à apparaître pendant la phase proliférative, les fibroblastes étant particulièrement dominants dans les derniers points temporels de la résolution des plaies.

À partir de la méthode #4, le protocole décrit les étapes d’utilisation de Seurat pour se concentrer sur un type de cellule majeur individuel dans l’ensemble de données unicellulaire afin d’identifier les sous-types cellulaires potentiels pendant la cicatrisation des plaies. Le protocole se concentre sur les fibroblastes, qui se sont initialement regroupés en deux groupes Seurat avant d’être combinés en une seule catégorie, et décrit comment créer un nouvel objet Seurat qui ne contient que les fibroblastes de l’ensemble de données d’origine (étape 4.1). Le flux de travail Seurat est effectué sur cet ensemble de données spécifique aux fibroblastes (étapes 4.2-4.4), l’étape 4.2 aboutissant à un graphique coudé (figure supplémentaire 9) montrant qu’une grande partie de la variation majeure de l’ensemble de données sur les fibroblastes se produit dans les 9 premières dimensions de l’ACP. À l’étape 4.5, la fonction de tracé dimensionnel est exécutée pour visualiser le regroupement des cellules sur un graphique UMAP, et les résultats représentatifs de cette étape (Figure 8) ont montré que les fibroblastes de l’ensemble de données étaient regroupés autour des 3 sous-types de cellules codés par couleur. La visualisation de l’ensemble de données sur les fibroblastes en fonction de leur annotation DPW (étape 4.6) a donné lieu à un graphique UMAP (figure supplémentaire 10) montrant les fibroblastes de l’ensemble de données répartis dans l’ensemble en fonction de leur annotation DPW. Le protocole décrit ensuite comment obtenir des listes de gènes exprimés différentiellement et les enregistrer dans un fichier texte (étape 4.7), ouvrir la table de données dans Excel et effectuer diverses étapes de filtrage afin d’obtenir les marqueurs de cluster les mieux classés pour chaque cluster de cellules (étape 4.8), et attribuer une nouvelle variable répertoriant les gènes marqueurs de fibroblastes supérieurs nommés « FB_type_marker » (étape 4.9). À l’étape 4.10, la fonction de diagramme à points est utilisée pour visualiser les gènes de la liste dans l’ensemble de données uniquement sur les fibroblastes en appelant la variable « FB_type_marker » dans le paramètre features, et les résultats représentatifs de cette étape (Figure 9) sont des diagrammes à points confirmant l’expression élevée des marqueurs de sous-type de fibroblastes uniquement dans leurs catégories de grappes respectives (en haut), mais équitablement répartis dans les catégories de DPW (en bas). À l’étape 4.11, la même variable de caractéristiques est appelée pour visualiser les gènes marqueurs des fibroblastes dans l’ensemble des données de cicatrisation des plaies, et le résultat représentatif (figure supplémentaire 11) est un graphique à points qui a confirmé l’expression élevée des marqueurs de sous-type de fibroblaste principalement dans le fibroblaste d’origine. Enfin, les étapes suivantes du protocole ont utilisé des diagrammes à barres empilées pour visualiser d’abord les proportions de DPW dans les trois sous-types de fibroblastes (étape 4.12), puis pour visualiser les proportions de sous-types de fibroblastes à différents points temporels (étape 4.13). Les résultats représentatifs de ces étapes sont des diagrammes de proportion montrant le nombre relatif de cellules DPW dans chaque catégorie de sous-types de fibroblastes (figure supplémentaire 12) et le nombre relatif de sous-types de fibroblastes dans chaque catégorie de DPW (figure supplémentaire 13). Ces résultats indiquent un changement significatif dans les proportions de sous-types de fibroblastes au cours de la cicatrisation, le premier sous-type de fibroblastes (groupe 0) étant fortement dominant dans les plaies à un stade précoce (D1 et D3), le deuxième sous-type (groupe 1) dominant pendant la résolution de la plaie (D14) et le troisième sous-type (groupe 2) étant le plus élevé pendant la phase proliférative de la cicatrisation (D7).

À partir de la méthode #5, le protocole décrit les étapes d’analyse d’un ensemble de données de cicatrisation de plaies unicellulaire à l’aide de la fonction de notation du module dans Seurat. Le protocole décrit d’abord les étapes d’utilisation d’un fichier texte séparé par des tabulations pour télécharger des ensembles de gènes dans des variables dans R (étapes 5.1-5.2), puis l’application de la fonction de notation du module à trois ensembles de gènes relatifs aux trois phases principales de la cicatrisation des plaies (étape 5.3). À l’étape 5.4, la fonction de diagramme à points est utilisée pour visualiser les scores agrégés des modules dans deux catégories de métadonnées différentes, et les résultats représentatifs de cette étape (Figure 10) sont des diagrammes à points montrant l’expression moyenne des principaux modules de la phase de cicatrisation dans les cellules de la catégorie des jours suivant la blessure (DPW, à droite) et dans la catégorie des principaux types de cellules (à gauche). Ces résultats montrent que l’application de profils d’expression génique basés sur le séquençage en vrac à des ensembles de données d’expression unicellulaire de manière pseudo-vrac est une méthode puissante pour les approches bioinformatiques comparatives en utilisant des ensembles de données précédemment publiés dans le domaine de la cicatrisation des plaies.

À partir de la méthode #6, le protocole décrit les étapes d’analyse d’un ensemble de données de cicatrisation de plaies unicellulaires dérivées de Seurat à l’aide du package et du flux de travail CellChat selon une question scientifique spécifique consistant à comparer les cellules dérivées de plaies précoces et avancées. Le protocole subdivise d’abord l’ensemble des données Seurat en deux points temporels après la blessure, l’un pendant la phase inflammatoire (jour 1 (J1)) et l’autre pendant la résolution de la plaie (jour 14 (J14)) (étape 6.1). Deux objets CellChat sont créés, et le protocole passe en revue toutes les fonctions typiques du protocole CellChat pour calculer toutes les interactions présumées entre les types de cellules identifiés dans la méthode #3 du protocole (étapes 6.2-6.3). À l’étape 6.4, la fonction de nuage de points de signalisation est exécutée pour visualiser les forces d’interaction entrantes et sortantes dans tous les principaux types de cellules à chaque point temporel de cicatrisation des plaies. Les résultats représentatifs de cette étape (Figure supplémentaire 14) sont des nuages de points qui montrent l’intensité des interactions entrantes (axe des y) et sortantes (axe des x) pour les principaux types de cellules aux points temporels D1 (à gauche) et J14 (à droite). Ces résultats ont montré que les cellules immunitaires telles que les neutrophiles et les macrophages avaient les plus fortes d’interaction cellule-cellule pendant la phase inflammatoire, mais que les fibroblastes dominaient les interactions cellule-cellule pendant la résolution des plaies, ce qui confirme des décennies de recherche sur la cicatrisation des plaies. Les étapes suivantes concentrent l’analyse sur l’une des voies significativement enrichies, la voie du collagène (étapes 6.5-6.6). À l’étape 6.7, la fonction de diagramme circulaire est exécutée pour visualiser les interactions de la voie de signalisation du collagène entre les types de cellules aux deux points temporels. Les résultats représentatifs de cette étape (Figure supplémentaire 15) sont des diagrammes circulaires qui montrent les interactions de signalisation de la voie du collagène déduites entre tous les types de cellules à D1 (à gauche) et J14 (à droite). À l’étape 6.8, les mêmes interactions sont visualisées à l’aide de la fonction de diagramme d’accords, avec les résultats représentatifs (Figure supplémentaire 16) étant des diagrammes d’accords montrant les interactions de signalisation de la voie du collagène déduites entre tous les types de cellules à chaque point temporel. Comme prévu, ces résultats ont montré que les fibroblastes étaient les cellules sources primaires de la voie de signalisation du collagène, bien que le flux d’information soit plus restreint aux cellules immunitaires à J1 qu’à J14. Pour se concentrer sur le fibroblaste en tant que cellule source dans les interactions cellule-cellule, l’étape 6.9 répète la fonction de diagramme de corde en ajoutant un paramètre de cellule source, et les résultats représentatifs (Figure supplémentaire 17) sont des diagrammes de cordes montrant les interactions de signalisation de la voie du collagène inférées avec les fibroblastes en tant que cellules sources à chaque point temporel. À l’étape 6.10, deux fonctions sont exécutées pour visualiser les contributions de chaque paire ligand-récepteur dans la voie de signalisation du collagène avec des fibroblastes comme cellules sources, l’une à l’aide de diagrammes à bulles (étape 6.10.1) et l’autre à l’aide de diagrammes d’accords (étape 6.10.2). Les résultats représentatifs montrent les contributions déduites de chaque paire ligand-récepteur dans la signalisation de la voie du collagène avec les fibroblastes comme cellules sources aux points temporels D1 (à gauche) et D14 (à droite) en utilisant les deux graphiques à bulles (Figure supplémentaire 18) et les diagrammes d’accords (Figure supplémentaire 19). Ces résultats ont montré qu’à D1, la voie du collagène provenant des fibroblastes était limitée aux neutrophiles et aux macrophages avec une dominance des récepteurs Cd44 et Sdc4, mais que dans D14, d’autres cellules agissaient comme récepteurs via une variété de récepteurs, y compris les intégrines. Pour se concentrer sur l’interaction ligand-récepteur Col1a1-Cd44, qui a montré de fortes forces dans les interactions fibroblastes, un paramètre est défini (étape 6.11) puis utilisé à l’étape 6.12 dans une fonction de diagramme d’accord pour visualiser cette interaction ligand-récepteur particulière entre tous les types de cellules, avec les résultats représentatifs (Figure supplémentaire 20) étant des diagrammes d’accords montrant les interactions ligand-récepteur Col1a1-Cd44 inférées entre tous les types de cellules aux points temporels D1 (à gauche) et D14 (à droite). Ces résultats ont montré que si dans D1, cette interaction est limitée aux fibroblastes en tant que cellules sources, dans D14, les macrophages et les cellules musculaires lisses agissent également comme cellules sources. Ensuite, le protocole décrit comment effectuer une analyse différentielle CellChat en fusionnant d’abord les objets CellChat D1 et D14 (étape 6.13). À l’étape 6.14, la fonction de comparaison des interactions est exécutée pour visualiser le nombre total et les forces relatives des interactions cellule-cellule entre les deux points de temps de cicatrisation, et les résultats représentatifs (Figure supplémentaire 21) sont les diagrammes à barres résultants montrant le nombre total (à gauche) et l’intensité (à droite) des interactions inférées dans les cellules comprenant des plaies D1 et D14, avec un nombre plus élevé d’interactions dans D14 par opposition à des forces relatives plus élevées dans D1. Dans les étapes 6.15 et 6.16, deux fonctions sont utilisées pour visualiser les forces d’interaction différentielles cellule-cellule entre chaque type de cellule lorsque la plaie passe du jour 1 au jour 14 avec leurs résultats représentatifs respectifs, la première étant un graphique circulaire (étape 6.15, Figure supplémentaire 22) et la seconde étant une carte thermique (étape 6.16, Figure supplémentaire 23), où l’augmentation des interactions en D14 par rapport à D1 indiquée en rouge et celles qui sont diminuées en bleu. Comme prévu, les interactions médiées par les neutrophiles et les macrophages sont augmentées dans D1, et les interactions médiées par les fibroblastes sont augmentées dans D14. À l’étape 6.17, la fonction de classement est utilisée pour créer un graphique qui classe les contributions relatives des voies individuelles aux interactions cellule-cellule avec les fibroblastes en tant que cellules sources à D14 par rapport à D1, et les résultats représentatifs (Figure supplémentaire 24) montrent le graphique de classement résultant avec D1 représenté en haut en rouge et D14 en bas en bleu, avec plusieurs voies étant exclusivement représentées en D1 ou D14 et beaucoup d’autres montrant un gradient d’activation. Enfin, à l’étape 6.18, deux fonctions de diagramme à bulles sont utilisées pour montrer les contributions relatives des paires ligand-récepteur individuelles dans la voie de signalisation du collagène avec des fibroblastes comme cellules sources à D14 par rapport à D1, avec les résultats représentatifs correspondants (Figure supplémentaire 25) montrant une augmentation (à gauche) et une diminution (à droite) des paires de signalisation dans D14 par rapport à D1 pour les nombreuses interactions cellule-cellule sur l’axe des x. Comme prévu, les fibroblastes présentaient beaucoup plus d’interactions ligand-récepteur sortantes à travers plusieurs cellules réceptrices dans les plaies D14 par rapport aux plaies D1, où la communication était plus limitée vers les neutrophiles et les macrophages pendant la phase inflammatoire.

À partir de la méthode #7, le protocole décrit les étapes d’intégration de deux ensembles de données de cicatrisation de plaies unicellulaires à l’aide de Seurat. Le protocole décrit d’abord les étapes de fusion de deux lots d’ensembles de données monocellulaires publiés et l’application du flux de travail Seurat standard à l’ensemble de données fusionné (étapes 7.1 à 7.4). À l’étape 7.5, la fonction de tracé dimensionnel est utilisée pour visualiser le tracé UMAP en fonction des numéros de cluster et de lot de l’ensemble de données de cicatrisation des plaies fusionné mais pas encore intégré. Les résultats représentatifs de cette étape (figure supplémentaire 26) sont des graphiques UMAP qui visualisent la distribution des clusters Seurat (à gauche) et des numéros de lot (à droite), montrant qu’il ne semble pas y avoir d’effets de lot significatifs pour ces deux ensembles de données avant l’intégration des données. Le protocole effectue ensuite l’intégration des données à l’aide de la méthode RPCA et du flux de travail Seurat de suivi de l’ensemble de données intégré (étapes 7.7-7.8). À l’étape 7.9, la fonction de tracé dimensionnel est utilisée pour visualiser le tracé UMAP en fonction des numéros de cluster et de lot de l’ensemble de données intégré de cicatrisation. Les résultats représentatifs de cette étape (figure supplémentaire 27) sont des graphiques UMAP qui visualisent la distribution des grappes Seurat (à gauche) et des numéros de lot (à droite), montrant qu’il y avait maintenant un chevauchement encore plus grand entre les deux lots à travers différentes grappes. Les résultats montrent également l’émergence d’un cluster supplémentaire à la suite de l’intégration des données, ce qui peut indiquer une capacité accrue à identifier des sous-types de cellules potentiellement significatifs après que les effets techniques des lots de données sont contrôlés.

Figure 1 : Graphique UMAP montrant que toutes les cellules de l’ensemble de données sont regroupées autour de 8 grands groupes de clusters codés par couleur. Résultats obtenus à partir d’un ordinateur fonctionnant sous Windows (à gauche) et MacOS (à droite). Ce chiffre correspond à l’étape 3.5. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 2 : Graphique UMAP montrant toutes les cellules de l’ensemble de données réparties en fonction de leur origine temporelle/spatiale, sans regroupement apparent selon l’annotation temps/espace. Ce chiffre correspond à l’étape 3.6. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 3 : Captures d’écran recadrées des sorties EnrichR, montrant les types de cellules enrichis les plus élevés pour chaque groupe de cellules. Ce chiffre correspond à l’étape 3.13. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 4 : Graphique UMAP montrant toutes les cellules de l’ensemble de données regroupées autour des principaux types de cellules codés par couleur. Ce chiffre correspond à l’étape 3.15. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 5 : Grille de graphiques UMAP montrant la forte expression des gènes marqueurs cellulaires supérieurs dans les principaux groupes de types cellulaires. Ce chiffre correspond à l’étape 3.16. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 6 : Graphiques à points confirmant le niveau élevé d’expression des gènes marqueurs cellulaires supérieurs uniquement dans leurs principaux types cellulaires respectifs. Ce chiffre correspond à l’étape 3.18. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 7 : Graphique des proportions montrant le nombre relatif des principaux types de cellules dans chaque catégorie de DPW. Ce chiffre correspond à l’étape 3.23. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 8 : Graphique UMAP montrant les fibroblastes de l’ensemble de données regroupés autour des 3 sous-types de cellules codés par couleur. Ce chiffre correspond à l’étape 4.5. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 9 : Graphiques à points confirmant l’expression élevée des marqueurs de sous-type de fibroblastes uniquement dans leurs catégories de grappes respectives, mais équitablement répartis dans les catégories DPW. Ce chiffre correspond à l’étape 4.10. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure 10 : Graphiques à points montrant l’expression moyenne des principaux modules de la phase de cicatrisation dans les cellules par DPW et par type de cellule majeur. Ce chiffre correspond à l’étape 5.4. Veuillez cliquer ici pour voir une version agrandie de cette figure.

Figure supplémentaire 1 : Résultats montrant qu’il existe une distribution assez uniforme des gènes détectés pour chaque code-barres, ce qui est important pour l’intégrité de l’ensemble de données et l’analyse en aval des points temporels de cicatrisation des plaies. Cette figure correspond à l’étape 2.6.11. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 2 : Nuages de points montrant qu’il existe un certain nombre de cellules avec un contenu mitochondrial important, ce qui est corrélé à un faible nombre d’ARN --- il s’agit de cellules mortes ou mourantes. Ce chiffre correspond à l’étape 2.10. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 3 : Nuages de points montrant que la distribution des gènes détectés et le pourcentage d’ARN mitochondrial par cellule sont maintenant plus normaux, ouvrant la voie à des analyses robustes en aval. Cette figure correspond à l’étape 2.12. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 4 : Graphique en violon montrant qu’il y a un certain nombre de cellules avec un score de doublet relativement élevé, et que 0,25 semble être une limite naturelle, au-dessus de laquelle il y a une population de doublets probables. Ce chiffre correspond à l’étape 2.15. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 5 : Graphique en coude montrant qu’une grande partie de la variation majeure se produit dans les 13 premières dimensions. Ce chiffre correspond à l’étape 3.2. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 6 : Graphique à points confirmant le niveau élevé d’expression des gènes marqueurs cellulaires supérieurs uniquement dans leurs groupes Seurat respectifs. Ce chiffre correspond à l’étape 3.17. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 7 : Graphique UMAP montrant la localisation des annotations de l’évolution temporelle de la plaie dans l’ensemble de données de cicatrisation des plaies. Ce chiffre correspond à l’étape 3.20. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 8 : Graphique des proportions montrant le nombre relatif de cellules DPW dans chaque grande catégorie de type de cellule. Ce chiffre correspond à l’étape 3.22. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 9 : Graphique du coude montrant qu’une grande partie de la variation majeure dans l’ensemble de données sur les fibroblastes se produit dans les 9 premières dimensions. Ce chiffre correspond à l’étape 4.2. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 10 : Graphique UMAP montrant les fibroblastes de l’ensemble de données répartis dans l’ensemble en fonction de leur annotation DPW. Ce chiffre correspond à l’étape 4.6. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 11 : Graphique à points confirmant l’expression élevée des marqueurs de sous-type de fibroblastes, principalement dans le groupe de fibroblastes d’origine. Ce chiffre correspond à l’étape 4.11. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 12 : Graphique des proportions montrant le nombre relatif de sous-types de fibroblastes dans chaque catégorie de DPW. Ce chiffre correspond à l’étape 4.12. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 13 : Graphique des proportions montrant le nombre relatif de fibroblastes dans l’ensemble du DPW dans chaque catégorie de sous-type de fibroblastes. Ce chiffre correspond à l’étape 4.13. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 14 : Nuages de points montrant l’intensité des interactions entrantes (axe des y) et sortantes (axe des x) pour les principaux types de cellules aux points temporels du jour 1 (D1, à gauche) et du jour 14 (D14, à droite). Ce chiffre correspond à l’étape 6.4. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 15 : Graphiques circulaires montrant les interactions de signalisation de la voie du collagène déduites entre tous les types de cellules dans chaque catégorie DPW. Ce chiffre correspond à l’étape 6.7. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 16 : Diagrammes d’accords montrant les interactions de signalisation de la voie du collagène déduites entre tous les types de cellules dans chaque catégorie de DPW. Ce chiffre correspond à l’étape 6.8. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 17 : Diagrammes d’accords montrant les interactions de signalisation de la voie du collagène déduites avec les fibroblastes en tant que cellules sources dans chaque catégorie DPW. Ce chiffre correspond à l’étape 6.9. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 18 : Graphiques à bulles montrant les contributions inférées de chaque paire ligand-récepteur dans la signalisation de la voie du collagène avec des fibroblastes comme cellules sources dans chaque catégorie DPW. Cette figure correspond à l’étape 6.10.1. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 19 : Diagrammes d’accords montrant les contributions inférées de chaque paire ligand-récepteur dans la signalisation de la voie du collagène avec les fibroblastes comme cellules sources dans chaque catégorie DPW. Cette figure correspond à l’étape 6.10.2. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 20 : Diagrammes d’accords montrant les interactions ligand-récepteur Col1a1-Cd44 déduites entre tous les types de cellules dans chaque catégorie de DPW. Ce chiffre correspond à l’étape 6.12. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 21 : Diagrammes à barres montrant le nombre (à gauche) et la force (à droite) des interactions inférées dans les blessures du jour 1 et du jour 14. Ce chiffre correspond à l’étape 6.14. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 22 : Graphique circulaire montrant les forces d’interaction différentielles entre les cellules entre chaque type de cellule lorsque la plaie passe du jour 1 (bleu) au jour 14 (rouge) DPW. Ce chiffre correspond à l’étape 6.15. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 23 : Carte thermique montrant les forces d’interaction différentielles entre les cellules entre chaque type de cellule lorsque la plaie passe du jour 1 (bleu) au jour 14 (rouge) DPW. Ce chiffre correspond à l’étape 6.16. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 24 : Graphique des rangs montrant les contributions relatives des voies individuelles aux interactions cellule-cellule entre les fibroblastes et d’autres types de cellules au jour 1 par rapport au jour 14 DPW. Ce chiffre correspond à l’étape 6.17. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 25 : Graphiques à bulles montrant les contributions relatives des paires ligand-récepteur individuelles dans la voie de signalisation du collagène avec les fibroblastes comme cellules sources au jour 1 par rapport au jour 14 DPW. Ce chiffre correspond à l’étape 6.18. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 26 : Graphiques UMAP montrant la distribution des clusters Seurat (à gauche) et les numéros de lot (à droite) avant l’intégration des données. Ce chiffre correspond à l’étape 7.5. Veuillez cliquer ici pour télécharger cette figure.

Figure supplémentaire 27 : Graphiques UMAP montrant la distribution des clusters Seurat (à gauche) et les numéros de lot (à droite) après intégration des données. Ce chiffre correspond à l’étape 7.9. Veuillez cliquer ici pour télécharger cette figure.

Dossier supplémentaire 1 : JoVE_Rscript.R : Fichier de script de code R principal, qui comprend toutes les étapes et explications décrites pour toutes les parties du protocole. Veuillez cliquer ici pour télécharger ce fichier.

Dossier supplémentaire 2 : JoVE_PhaseSpecificGenes.txt. Fichier texte délimité par des tabulations, qui contient les listes des gènes chargés à l’étape 5.1 du protocole. Veuillez cliquer ici pour télécharger ce fichier.

Dossier supplémentaire 3 : JoVE_Rscript_b3.R. Fichier de script de code R supplémentaire, qui comprend toutes les étapes et explications requises pour analyser le lot #3 de l’ensemble de données à utiliser à l’étape 7.1 du protocole. Veuillez cliquer ici pour télécharger ce fichier.

Tableau supplémentaire 1 : JoVE_DEGs_cellMarkers.xlsx. Fichier Excel, qui contient la sortie complète des gènes exprimés différentiellement classés utilisés à l’étape 3.10 du protocole. Veuillez cliquer ici pour télécharger ce tableau.

Tableau supplémentaire 2 : Les 5 principaux gènes régulés et exprimés à la hausse pour chaque groupe de seurat. Veuillez cliquer ici pour télécharger ce tableau.

Les auteurs n’ont aucun conflit d’intérêts à divulguer.