Un flux de travail pratique pour l'analyse des données de transcriptomique spatiale : de l'acquisition des données aux analyses avancées

Résumé

Ce protocole présente un flux de travail reproductible pour l'analyse de données de transcriptomique spatiale, guidant les utilisateurs depuis l'acquisition de données publiques et le contrôle qualité basé sur Seurat jusqu'à l'intégration, la détection de caractéristiques spatiales, la déconvolution des types cellulaires, l'annotation de régions d'intérêt et l'analyse des communications cellule-cellule, avec des points de contrôle pratiques qui soutiennent une exécution transparente.

Résumé

La transcriptomique spatiale (ST) caractérise l'expression génique à l'échelle du génome tout en préservant le contexte spatial bidimensionnel des molécules d'ARNm au sein des sections de tissus, permettant ainsi l'étude de l'architecture tissulaire et de la biologie associée au microenvironnement. Toutefois, l'analyse ST demeure difficile, car l'importation des données, le contrôle de qualité, l'intégration, la déconvolution, les statistiques spatiales et la visualisation nécessitent souvent plusieurs environnements logiciels ainsi que des choix de paramètres reproductibles. Ce protocole présente un flux de travail computationnel pratique pour les jeux de données ST publics en R, en commençant par l'acquisition des données et la configuration des logiciels, puis en poursuivant avec le chargement des données, le contrôle de qualité, la normalisation, l'intégration multi-échantillons, le regroupement (clustering) et l'analyse des gènes à variabilité spatiale basés sur Seurat. Le flux de travail applique ensuite des stratégies complémentaires de déconvolution, notamment une analyse SPOTlight guidée par référence et une modélisation thématique STdeconvolve sans supervision, suivies d'une analyse de communication spatiale entre cellules basée sur Giotto et d'une sélection interactive de régions d'intérêt (ROI) à l'aide d'une application Python Dash personnalisée. En mettant l'accent sur l'exécution par script, les justifications explicites des paramètres, les résultats attendus et des points de contrôle pour le dépannage, ce protocole fournit un cadre adaptable aux jeux de données ST standards basés sur puces d'éléments ainsi qu'à des plateformes apparentées, après évaluation des paramètres spécifiques à chaque jeu de données et à chaque plateforme.

Introduction

La transcriptomique spatiale (ST) est une famille révolutionnaire de technologies qui mesure l'expression génique à l'échelle du génome tout en conservant les coordonnées spatiales des molécules d'ARN messager (ARNm) au sein des sections de tissus. Les méthodes de ST comprennent des approches fondées sur le séquençage utilisant des réseaux à codes-barres positionnels et des approches d'imagerie in situ qui cartographient les signaux transcriptionnels dans des microenvironnements tissulaires intacts^1,2. En préservant le contexte spatial, la ST permet l'analyse de l'architecture tissulaire, de l'organisation des voisinages cellulaires, de la communication entre cellules et des processus biologiques associés au microenvironnement, qui ne peuvent pas être entièrement résolus après la dissociation du tissu³.

La croissance rapide des dépôts publics de données de transcriptomique spatiale (ST) a créé des opportunités sans précédent pour l'analyse secondaire et le développement de méthodes³. Des ressources telles que la base de données CROST regroupent des centaines de jeux de données transcriptomiques à résolution spatiale provenant de multiples espèces et de diverses plates-formes technologiques, tandis que des collections spécialisées comme STOmicsDB se concentrent sur des méthodologies spécifiques telles que le Stereo-seq^4,5. Malgré l'abondance de données, l'analyse computationnelle demeure difficile en raison de la complexité des structures des données spatiales, de la diversité des outils d'analyse et des obstacles techniques liés à la mise en œuvre de flux de travail reproductibles^{6,7,8,9,10,11}.

Pour pallier les limites liées à l'utilisation d'un seul environnement logiciel, un flux de travail computationnel intégré utilisant des outils d'analyse complémentaires est présenté ici. Les écosystèmes d'analyse de données de séquençage spatial (ST) existants comprennent principalement Seurat, Giotto et des cadres basés sur Python tels que Squidpy^6,7,12. Bien que les outils basés sur Python comme Squidpy offrent des fonctionnalités étendues pour l'analyse de graphes spatiaux, regrouper le pipeline computationnel principal au sein d'un seul environnement de langage de programmation permet de minimiser les obstacles techniques liés à l'interopérabilité entre langages. En conséquence, le pipeline central est principalement implémenté en R afin de réduire ces difficultés techniques inter-langages. Dans ce flux de travail basé sur R, Seurat est utilisé pour le chargement des données, le contrôle de qualité, la normalisation, la réduction de dimensionnalité, la visualisation et l'intégration de plusieurs échantillons, reflétant ainsi son utilisation courante dans les pipelines de transcriptomique unicellulaire et spatiale. Giotto est ensuite utilisé pour la construction du réseau spatial et l'analyse des communications cellule-cellule basée sur les interactions ligand-récepteur. Ce pipeline établit donc un lien entre le prétraitement et l'intégration fondés sur Seurat et l'analyse spatiale fondée sur Giotto, tout en rendant explicite et reproductible le transfert des données entre ces deux ensembles d'outils.

Dans ce cadre, deux stratégies complémentaires de déconvolution sont mises en œuvre : SPOTlight, une méthode guidée par référence qui utilise les données de scRNA-séquençage pour estimer les proportions de types cellulaires, et STdeconvolve, une méthode de modélisation thématique non supervisée qui identifie des profils transcriptionnels latents^8,11. Les résultats fournissent des perspectives complémentaires sur l'hétérogénéité cellulaire spatiale, mais ne sont pas considérés comme une validation croisée quantitative, sauf si l'utilisateur effectue l'analyse de concordance facultative décrite dans le protocole. Une application Python Dash personnalisée, Select Spatial Spots, est intégrée pour l'annotation interactive de régions d'intérêt (ROI) et exporte des fichiers d'annotation standard basés sur des coordonnées, utilisables dans des analyses ultérieures en aval.

En ce qui concerne l'applicabilité pratique, ce flux de travail est conçu principalement pour les données spatiales standard basées sur des puces (par exemple, Visium avec une résolution de 55 µm) et pourrait être adapté à d'autres types de tissus après évaluation des paramètres. Certaines limitations importantes doivent être prises en compte avant l'analyse. Premièrement, le module de désconvolution guidée par référence dépend d'une référence de séquençage d'ARN monocellulaire (scRNA-seq) de haute qualité, appariée au tissu étudié. Deuxièmement, les plateformes subcellulaires ou quasi monocellulaires peuvent nécessiter un prétraitement modifié, un regroupement spatial des bins ou une segmentation cellulaire basée sur l'image avant l'intégration². Le jeu de données représentatif du côlon de souris est utilisé comme cas démonstratif pour illustrer comment le flux de travail permet d'évaluer les domaines spatiaux et l'organisation tissulaire définie par des marqueurs, et non comme preuve d'une compatibilité universelle avec la plateforme.

Protocole

Toutes les banques de données biologiques analysées dans ce protocole sont accessibles au public et utilisées uniquement à des fins de démonstration. Les numéros d'accès spécifiques aux données ainsi que les référentiels d'origine sont indiqués dans les étapes correspondantes. Les jeux de données originaux ont été produits par les chercheurs initiaux conformément aux lignes directrices éthiques institutionnelles applicables à chaque étude source. Consultez le Tableau des matériaux pour vérifier toutes les versions requises des logiciels et des packages R.

Exigences matérielles : La mémoire informatique requise pour ce flux de travail augmente en fonction du nombre d'échantillons et de points analysés. Pour un jeu de données typique de transcriptomique spatiale (par exemple, environ 3 000 points par échantillon sur un maximum de trois échantillons), une station de travail standard disposant d'au moins 16 Go de RAM suffit pour exécuter le pipeline. Toutefois, 32 Go de RAM ou plus sont fortement recommandés afin d'assurer des performances et une stabilité optimales, en particulier lors des étapes de calcul intensives en mémoire telles que la normalisation par SCTransform et la factorisation matricielle durant la déconvolution.

1. Acquisition des données et préparation de la structure des répertoires

1. Obtenir des jeux de données publics de transcriptomique spatiale (par exemple, l'accession GEO GSE169749, avec une date de soumission du 26 mars 2021 et une dernière mise à jour du 6 mars 2022) et des données de référence d'ARN-séquençage en cellule unique (par exemple, l'accession GEO GSE264408, avec une date de soumission du 19 avril 2024 et une dernière mise à jour du 10 décembre 2024).
2. Téléchargez l'archive des données brutes (par exemple, GSE169749_RAW.tar) depuis le référentiel et décompressez l'archive.
   REMARQUE : Pour l'ensemble de données représentatif utilisé dans ce protocole, le répertoire extrait contient plusieurs fichiers pour chaque échantillon, les noms de fichiers incluant généralement l'identifiant d'échantillon GEO (par exemple, GSM5213483).
3. Organisez les fichiers selon une structure de répertoires normalisée compatible avec Seurat Load10X_Espace fonction.
   1. Créez un répertoire principal pour les données (par exemple, ./data/).
   2. Créez un sous-répertoire dédié pour chaque échantillon (par exemple, ./data/sample_A1/).
   3. Transférez (copiez ou déplacez) les fichiers essentiels suivants pour chaque échantillon dans le sous-répertoire correspondant.
      1. Placer le filtered_feature_bc_matrix.h5 transférez dans le sous-répertoire principal de l'échantillon.
      2. Créer un spatial/ sous-dossier dans chaque répertoire d'échantillon.
      3. Placez les fichiers suivants dans le spatial/ sous-dossier : tissue_positions_list.csv.gz, scalefactors_json.json.gz et tissue_hires_image.png.gz.
   4. Décompressez le .gz fichiers dans le spatial/ dossier. Veillez à ce que les noms de fichiers d'origine restent exactement tels que requis par le Load10X_Espace fonction (par ex., tissue_positions_list.csv).
      REMARQUE : La structure finale et organisée pour un chargement sans problème doit suivre ce modèle, en prenant l'exemple de l'échantillon A1 :
      ./data/échantillon_A1/matrice_bc_filtre_avec_caractéristiques.h5
      ./data/échantillon_A1/spatial/positions_tissu_liste.csv
      ./data/sample_A1/spatial/facteurs_echelle_json.json
      ./data/échantillon_A1/spatial/image_tissu_haute_résolution.png

2. Configuration de l'environnement logiciel

1. Installez R (version 4.4.3 ou supérieure).
2. Obtenez les scripts d'analyse (1_ReadSpatialData.R à 8_Giotto_Communication.R) depuis le dépôt GitHub (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow, au commit d'empreinte 2d85e18 pour assurer la reproductibilité).
3. Installez les packages R requis depuis CRAN et Bioconductor en exécutant le script setup.R.
   NOTE : Une liste complète des versions spécifiques des packages utilisés dans ce flux de travail est fournie dans le Tableau des matériaux et est également disponible dans le fichier session_info.txt du dépôt GitHub associé (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow).
4. Installez la Giotto Suite en exécutant les commandes d'installation fournies dans la documentation officielle.
   1. Accédez à l'URL d'installation officielle pour récupérer les scripts de configuration : https://giottosuite.com/articles/installation.html
   2. Initialisez l'environnement Python requis et les dépendances système conformément aux instructions étape par étape.
5. Installez des packages supplémentaires pour une visualisation améliorée.
   1. Exécutez : install.packages(c("ggprism", "colorBlindness"))
6. Installez l'outil personnalisé Select Spatial Spots. Assurez-vous que le système d'exploitation est Windows, macOS ou Linux, et qu'une version de Python 3.8 ou supérieure est installée. Cette étape est facultative.
   1. Obtenez l'outil en accédant au dépôt GitHub (https://github.com/LeafLight/SelectSpatialSpots, au commit d'empreinte d20946e pour assurer la reproductibilité) et en téléchargeant le code source.
   2. Accédez au répertoire de l'outil et installez les dépendances Python : pip install -r requirements.txt
7. Installez le package glmGamPoi en exécutant BiocManager::install("glmGamPoi") afin d'accélérer la normalisation SCTransform.
   NOTE : Le flux de travail principal nécessite uniquement les packages listés aux étapes 2.1–2.3 et 2.6. Le package Rfast2 est utilisé pour accélérer le calcul de la statistique de Moran's I. Les packages facultatifs de l'étape 2.5 permettent de générer des graphiques prêts pour publication avec un thème de style prism (ggprism) et d'accéder à une palette de couleurs adaptée aux personnes malvoyantes (colorBlindness). L'outil de l'étape 2.6 active la fonctionnalité de sélection interactive des spots. Exécutez la commande sessionInfo() et enregistrez la sortie complète de la console dans un fichier texte. Cela documente précisément les versions des logiciels et les dépendances des packages, éléments essentiels pour garantir une reproductibilité à long terme.

3. Chargement des données spatiales et contrôle de qualité (1_ReadSpatialData.R, 2_SpatialDataQC.R)

1. Lire les données spatiales dans des objets Seurat.
   1. Utiliser Read10X_Image pour charger manuellement l'image tissulaire haute résolution, en spécifiant le image.dir et image.name.
   2. Utiliser Load10X_Spatial avec le paramètre image défini sur l'objet image créé à l'étape 3.1.1 afin de créer l'objet Seurat.
2. Calculer les métriques de contrôle qualité.
   1. Calculer le pourcentage de lectures mitochondriales à l'aide de PercentageFeatureSet avec le motif ^mt-.
3. Visualiser et interpréter les données selon les métriques de contrôle qualité.
   1. Générer des diagrammes en violon de nCount_Spatial, nFeature_Spatial et percent.mt à l'aide de VlnPlot.
   2. Créer des diagrammes spatiaux des caractéristiques pour ces métriques à l'aide de SpatialFeaturePlot afin d'identifier les spots situés en dehors de la zone tissulaire.
      NOTE : Point de contrôle visuel : les spots situés en dehors de la zone tissulaire affichent généralement un faible nombre de lectures UMI (nCount_Spatial < 500) et une faible détection de gènes (nFeature_Spatial < 200).
   3. À des fins de démonstration, appliquer des filtres pour supprimer les spots de faible qualité (par exemple, subset(seurat_obj, subset = nFeature_Spatial > 200 & nCount_Spatial > 500)). Cette étape est facultative.
      NOTE : L'objectif principal du contrôle qualité en transcriptomique spatiale est d'identifier et d'annoter les artefacts techniques, tels que les spots situés en dehors du tissu. Pour les données Visium du côlon utilisées en démonstration, les seuils nFeature_Spatial > 200 et nCount_Spatial > 500 permettent d'éliminer les spots à faible complexité ou correspondant au bruit de fond. Pour les tissus à faible teneur en ARN ou dégradés, examiner les diagrammes en violon et les diagrammes spatiaux des caractéristiques avant d'augmenter les seuils ; pour les tissus denses ou à nombre élevé d'UMI, des seuils plus stricts peuvent être appropriés. Éviter de filtrer uniquement sur des valeurs d'expression élevées, sauf si des doublets, des plis tissulaires ou des artefacts évidents sont visibles. Le filtrage n'est souvent pas recommandé car il supprime physiquement des emplacements spatiaux, ce qui peut perturber la continuité de la structure tissulaire pour les analyses spatiales ultérieures.

4. Prétraitement, intégration et regroupement des données (3_IntegrationAndClustering.R)

1. Normaliser et prétraiter les échantillons individuels.
   1. Appliquer la normalisation SCTransform à chaque échantillon séparément avec assay = Spatial.
2. Intégrer plusieurs échantillons.
   1. Préparer la liste des objets normalisés par SCT pour l'intégration. S'assurer que chaque objet possède un assay « RNA » en copiant l'assay « Spatial » : spatial_list[[1]][[« RNA »]] <- spatial_list[[1]][[« Spatial »]]
      ​REMARQUE : Copier l'assay « Spatial » dans un emplacement d'assay « RNA » standard constitue une solution nécessaire pour assurer la compatibilité avec les fonctions d'intégration en aval initialement conçues pour les objets Seurat unicellulaires.
   2. Utiliser SelectIntegrationFeatures et PrepSCTIntegration pour identifier les caractéristiques variables partagées.
   3. Trouver les ancres d'intégration à l'aide de FindIntegrationAnchors avec normalization.method = « SCT ».
   4. Intégrer les données à l'aide de IntegrateData avec normalization.method = « SCT ».
3. Réaliser une réduction de dimensionnalité et un regroupement sur l'assay intégré.
   1. Exécuter l'ACP sur les données intégrées à l'aide de RunPCA.
   2. Déterminer le nombre optimal de composantes principales (CP) pour l'analyse en aval en calculant la variance cumulée expliquée. Identifier le point du coude de façon programmée (par exemple, le point où la variance cumulée dépasse 90 % et où le gain marginal tombe en dessous de 0,1 %), en utilisant le code dans 3_IntegrationAndClustering.R, lignes 36-38. Le nombre résultant de CP sera désormais désigné par pc.use.
   3. Exécuter RunUMAP avec dims = 1:pc.use.
   4. Regrouper les cellules à l'aide de FindNeighbors avec dims = 1:pc.use et FindClusters avec resolution = 0,5. Ajuster la résolution uniquement après avoir examiné la stabilité des groupes et la cohérence des gènes marqueurs.
   5. Effectuer une analyse d'expression différentielle entre les groupes cibles (par exemple, « B1_colon_d14 » versus « A1_colon_d0 ») à l'aide de la fonction FindMarkers.
      ​REMARQUE : Point de contrôle visuel : une intégration réussie affichera un mélange approprié des échantillons sur le graphique UMAP tout en conservant des groupes biologiquement distincts.
4. Identifier les gènes variables spatialement.
   1. Pour chaque échantillon d'origine, exécuter FindSpatiallyVariableFeatures en utilisant la méthode « moransi » sur l'assay « SCT » afin de calculer l'autocorrélation spatiale.

5. Prétraitement des données de référence à cellule unique (4_scDataPreProcessing.R)

1. Lire la matrice de comptage d'ARN-séquençage en cellule unique à l'aide de Read10X et créer un objet Seurat.
2. Effectuer le contrôle qualité standard, la normalisation et le regroupement en clusters.
   1. Calculer le pourcentage de lectures mitochondriales (percent.mt) et filtrer les cellules (par exemple, nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 25).
   2. Normaliser les données à l'aide de SCTransform, en définissant vars.to.regress = "percent.mt".
   3. Exécuter l'ACP, l'UMAP et regrouper les cellules en utilisant la méthode dynamique de sélection des composantes principales décrite à l'étape 4.3.2.
3. Annoter les types cellulaires.
   1. Calculer les scores de module pour les gènes marqueurs canoniques des types cellulaires à l'aide de AddModuleScore.
   2. Annoter les clusters en fonction des scores de module et des connaissances biologiques établies. En alternative, importer des annotations pré-calculées à partir des métadonnées.

6. Décodage guidé par référence avec SPOTlight (5_SPOTlight_Deconv.R)

1. Préparer les données pour SPOTlight.
   1. Convertir l'objet Seurat monocellulaire annoté et l'objet Seurat spatial en objets SingleCellExperiment.
   2. Normaliser les données monocellulaires par logarithme à l’aide de logNormCounts.
2. Exécuter la déconvolution SPOTlight.
   1. Identifier les gènes hypervariables (HVG) dans les données monocellulaires à l’aide de modelGeneVar et getTopHVGs.
   2. Calculer les gènes marqueurs spécifiques aux types cellulaires à l’aide de scoreMarkers et filtrer les marqueurs de haute qualité (par exemple, mean.AUC > 0,8).
   3. Sous-échantillonner la référence monocellulaire pour chaque type cellulaire à un nombre gérable (par exemple, 50 cellules) afin de réduire le temps de calcul.
   4. Exécuter la déconvolution à l’aide de la fonction SPOTlight avec weight_id = "mean.AUC", group_id = "cluster", et gene_id = "gene".
3. Visualiser et exporter les résultats.
   1. Extraire la matrice de résultats de déconvolution (proportions des types cellulaires par spot).
   2. Utiliser plotSpatialScatterpie pour visualiser la composition cellulaire sur les coordonnées spatiales.
   3. Ajouter les résultats de déconvolution aux métadonnées de l’objet Seurat spatial à l’aide de AddMetaData. Résultat attendu : une matrice de déconvolution comportant des colonnes de proportions de types cellulaires pour chaque spot spatial, des graphiques en camembert spatiaux montrant la composition cellulaire locale, et un objet Seurat contenant les proportions de déconvolution comme métadonnées.

7. Dédconvolution sans référence avec STdeconvolve (7_STdeconvolve.R)

1. Préparer les données spatiales.
   1. Extraire la matrice de comptage brute à partir de l'objet spatial Seurat en utilisant GetAssayData avec slot = "comptages".
   2. Supprimer les spots et les gènes de faible qualité à l'aide de cleanCounts de STdeconvolve (par ex., min.lib.size = 100).
2. Identifier les types cellulaires latents.
   1. Filtrer davantage le corpus pour les gènes exprimés dans une fraction minimale de spots en utilisant restreindreCorpus (par exemple, removeAbove=1.0, removeBelow = 0,05).
   2. Ajuster un modèle d’allocation latente de Dirichlet (LDA) pour une gamme de nombres potentiels de sujets (K) (par ex. Ks = seq(2, 9, par = 1)) en utilisant fitLDA.
   3. Sélectionner le modèle optimal en fonction de la perplexité minimale en utilisant modèle optimal avec opt = "min".
3. Analyser et visualiser les résultats.
   1. Extraire les proportions de types cellulaires (thêta) et les profils géniques (bêta) à partir du modèle optimal en utilisant getBetaTheta.
   2. Pour faciliter l'interprétation biologique des thèmes déconvolués, importez des annotations de régions d'intérêt (ROI) générées par le "Sélectionner les spots spatiaux" outil (voir étape 9.3.3). Utilisez ces annotations comme paramètre de groupes dans le vizTousLesSujets fonction permettant de projeter les proportions de types cellulaires déconvoluées sur les coordonnées spatiales et de coder par couleur les points selon leur région d'intérêt (ROI). Cette étape est facultative.
      REMARQUE : Le résultat attendu est que STdeconvolve renvoie des valeurs theta représentant les proportions de thèmes par point et des valeurs beta représentant les profils génétiques pour chaque thème. Si les utilisateurs souhaitent comparer quantitativement les résultats non supervisés de STdeconvolve avec les sorties de SPOTlight guidées par référence (de l'étape 6), ils peuvent exporter les deux matrices de proportions (theta issue de STdeconvolve et la matrice de SPOTlight) et calculer des mesures de corrélation ou de concordance par point (par exemple, corrélation de Pearson ou de Spearman) à l'aide de fonctions R standards. La comparaison quantitative entre SPOTlight et STdeconvolve est facultative et n'est pas requise pour exécuter le flux de travail principal ; toutefois, les matrices de sortie pertinentes sont disponibles aux étapes 6.3.1 et 7.3.1 pour les utilisateurs souhaitant réaliser une telle analyse.

8. Communication cellulaire spatiale à l'aide de Giotto(8_Giotto_Communication.R)

1. Convertir l'objet Seurat en un objet Giotto.
   1. Utiliser la fonction createGiottoObject, en fournissant la matrice brute de comptages et les coordonnées spatiales.
2. Prétraiter l'objet Giotto et ajouter les résultats de la déconvolution.
   1. Normaliser les données à l'aide de normalizeGiotto.
   2. Ajouter les annotations des types cellulaires (par exemple, le type cellulaire majeur issu de la déconvolution SPOTlight) aux métadonnées cellulaires à l'aide de addCellMetadata.
3. Inférer la communication cellule-cellule informée spatialement.
   1. Créer un réseau spatial à l'aide de createSpatialNetwork avec method = "Delaunay" et name = "Delaunay_network".
   2. Charger une base de données de ligands-récepteurs dans l'environnement R. S'assurer que l'objet chargé est formaté comme une trame de données avec deux colonnes distinctes représentant les symboles génétiques du ligand et du récepteur (par exemple, un réseau souris pré-établi provenant de Zenodo : https://zenodo.org/api/records/15168114/files/lr_network_mouse.csv/content, version v6, publié le 7 avril 2025).
   3. Exécuter exprCellCellcom avec cluster_column = "celltype_major" pour identifier les interactions significatives ligand-récepteur entre les types cellulaires situés en proximité spatiale. Résultat attendu : un tableau contenant des paires ligand-récepteur, des combinaisons de types cellulaires sources et cibles, des valeurs de log2 fold-change et des valeurs de P ajustées ; conserver les interactions significatives (par exemple, p.adj < 0,05) pour la visualisation.

9. Sélection interactive de taches avec la fonction select spatial spots (6_SelectSpatialSpots.R)

1. Préparer les données pour l'outil interactif.
   1. Extraire les coordonnées spatiales de l'objet Seurat à l'aide de GetTissueCoordinates.
2. Formater et exporter les données.
   1. Formater le tableau de données des coordonnées afin qu'il contienne exactement les en-têtes de colonne suivants : CELL_ID, X et Y.
   2. Exporter le tableau de données formaté vers un fichier CSV.
3. Effectuer l'analyse de la région d'intérêt (ROI).
   1. Lancer l'application personnalisée Dash Select Spatial Spots et charger le fichier CSV.
   2. Sélectionner interactivement des points en fonction de leur emplacement spatial.
   3. Exporter la liste des points sélectionnés et leurs étiquettes de groupe/ROI attribuées sous forme d'un nouveau fichier CSV. Résultat attendu : un fichier CSV contenant les identifiants des points, les coordonnées X/Y et les étiquettes de groupe/ROI attribuées, pouvant être appariés à l'objet Seurat via CELL_ID.

Résultats

Mise en œuvre du flux de travail et intégration des données illustrant les caractéristiques majeures des tissus

Le flux de travail informatique a été appliqué à des données de transcriptomique spatiale du côlon de souris afin d'illustrer les résultats attendus à chaque étape analytique. Comme indiqué dans le schéma du flux de travail (Figure 1), le pipeline débutait par l'acquisition des données et le contrôle de qualité, où les graphiques spatiaux des caractéristiques délimitaient les frontières du tissu (Figure 2A,B). Le flux de travail d'intégration basé sur des ancres de Seurat a ensuite été utilisé pour réduire les effets techniques liés aux lots tout en préservant la variation biologique interprétable. Les visualisations UMAP montraient l'alignement des échantillons et les motifs de regroupement spatial après intégration (Figure 2C,D). Une sélection quantitative et dynamique des composantes principales (PC) fondée sur la variance cumulée a été mise en œuvre afin de guider la réduction dimensionnelle et le regroupement ultérieur (voir Figure supplémentaire 1). L'analyse du heatmap des gènes marqueurs a révélé des profils transcriptionnels distincts sous-tendant les groupes spatiaux (Figure 2E).

Afin d'évaluer si les groupes calculés étaient cohérents avec l'architecture anatomique connue de l'histologie du côlon, les profils d'expression de gènes marqueurs canoniques spécifiques aux couches ont été analysés. La couche de l'épithélium muqueux présentait une expression de marqueurs de cellules épithéliales, notamment Epcam et Krt8, ainsi que du marqueur des cellules caliciformes Muc2. Les marqueurs mésenchymateux et stromaux tels que Col1a1 et Vim identifiaient la lame propre et les régions sous-muqueuses, tandis que la couche externe de la musculeuse propre était indiquée par des gènes structuraux des muscles lisses tels que Acta2 et Tagln. La restriction spatiale de ces marqueurs associés aux lignées soutient l'interprétation selon laquelle le flux de travail d'intégration et de regroupement a préservé les principales laminations histologiques du tissu colique selon l'axe muqueuse-musculeuse (voir Supplementary Figure 2).

Après la validation des groupes, une analyse différentielle de l'expression a été réalisée afin d'identifier les gènes différentiellement exprimés (DEG) entre les conditions expérimentales (Figure 2F,G). De plus, les gènes présentant une variabilité spatiale ont été identifiés à l'aide de la statistique de Moran I, mettant en évidence les gènes dont la distribution spatiale dans le tissu est significativement non aléatoire (Figure 2H).

La désconvolution cellulaire et les réseaux d'interaction spatiale révèlent la microorganisation tissulaire

Le traitement des données de référence de séquençage de l'ARN unicellulaire a produit des annotations étayées par un filtrage de contrôle qualité (Figure 3A), un regroupement non supervisé (Figure 3B), une validation par gènes marqueurs (Figure 3C) et une concordance avec des annotations indépendantes (Figure 3D). La composition cellulaire (Figure 3E) a guidé la stratégie de sous-échantillonnage pour la déconvolution. SPOTlight a estimé les proportions de types cellulaires guidées par référence à travers les spots spatiaux (Figure 4A,B), tandis que STdeconvolve a fourni une vue, basée sur un modèle thématique non supervisé, des motifs cellulaires spatiaux (Figure 5B). L'outil personnalisé Select Spatial Spots a fourni un contexte histologique à ces motifs (Figure 5A). Enfin, à l'aide des attributions de types cellulaires déconvoluées, l'analyse de communication spatiale a identifié des interactions ligand-récepteur entre des groupes de types cellulaires proches spatialement (Figure 6A,B).

Observations de dépannage issues de l'optimisation du protocole

Lors de l'optimisation du protocole, plusieurs problèmes ont été identifiés et ont permis d'établir des points de contrôle pratiques. Des résultats de déconvolution sous-optimaux ont été observés lorsque les références unicellulaires étaient mal appariées au contexte tissulaire, ce qui souligne la nécessité d'utiliser des données scRNA-seq appariées au tissu et à l'espèce lorsque cela est possible. Les premières tentatives de regroupement avec les paramètres par défaut ne permettaient pas toujours de distinguer les structures biologiques attendues ; l'analyse du choix des composantes principales, de la résolution du regroupement et de la cohérence des gènes marqueurs a permis d'identifier des domaines interprétables spatialement, conformes à l'anatomie du tissu. Ces observations fournissent des exemples pratiques de la manière dont les utilisateurs peuvent diagnostiquer des problèmes analytiques courants lors de l'exécution du flux de travail.

Figure 1: Flux de travail pour l'analyse intégrée de la transcriptomique spatiale. Représentation schématique du pipeline analytique, allant de l'acquisition des données et du prétraitement jusqu'aux analyses spatiales avancées. Les étapes clés comprennent : (1) le chargement des données, le contrôle de qualité et l'intégration multi-échantillons à l'aide de Seurat ; (2) le regroupement spatial et la détection des gènes à variabilité spatiale ; (3) la déconvolution des types cellulaires via méthodes basées sur une référence (SPOTlight) et méthodes non supervisées (STdeconvolve) ; (4) Analyse de la communication spatiale entre cellules à l'aide de Giotto et sélection interactive de régions d'intérêt à l'aide d'un outil personnalisé, Sélectionner les spots spatiauxLes résultats de tous les modules sont synthétisés afin de dégager des informations biologiques sur l'architecture tissulaire et le microenvironnement cellulaire. Veuillez cliquer ici pour afficher une version agrandie de cette figure.

Figure 2 : Intégration des données, regroupement et analyse d'expression différentielle. (A,B) Métriques de contrôle qualité pour les échantillons spatiaux A1 et B1, montrant les distributions des nombres de gènes, des nombres d'UMI et des pourcentages de gènes mitochondriaux. (C) Visualisation UMAP des données d'expression transcriptionnelle spatiale intégrées, colorées selon l'origine des échantillons (à gauche) et l'identité des regroupements (à droite). (D) Projection spatiale des identités des groupes sur les sections tissulaires. (E) Carte thermique des gènes marqueurs les plus significatifs pour chaque groupe spatial. (F) Graphique en volcan illustrant les gènes différentiellement exprimés entre les conditions A1_colon_d0 et B1_colon_d14. (G) Profils d'expression spatiale de gènes différentiellement exprimés représentatifs à travers les sections tissulaires. (H) Cartes d'expression spatiale des gènes les plus variables spatialement identifiés via la statistique de Moran's I, les deux panneaux de gauche montrant des gènes de l'échantillon A1_colon_d0 et les deux panneaux de droite montrant des gènes de l'échantillon B1_colon_d14. Veuillez cliquer ici pour consulter une version agrandie de cette figure.

Figure 3 : Traitement et annotation des données de référence à cellule unique. (A) Métriques de contrôle de qualité des données de référence scRNA-seq avant et après filtration.  (B) Visualisation UMAP des données scRNA-seq colorées selon les clusters non supervisés.  (C) Graphique en points montrant les scores d'expression des gènes marqueurs canoniques des types cellulaires selon les clusters.  (D) Visualisation UMAP annotée des données scRNA-seq avec les principaux types cellulaires identifiés.  (E) Composition cellulaire du jeu de données de référence scRNA-seq. La ligne rouge pointillée indique le seuil de sous-échantillonnage (n = 50 cellules par type) appliqué lors de la déconvolution SPOTlight afin d'équilibrer efficacité computationnelle et représentation des types cellulaires. Veuillez cliquer ici pour visualiser une version agrandie de cette figure.

Figure 4 : Déconvolution spatiale de l'hétérogénéité cellulaire. (A,B) Graphiques en camembert spatiaux issus de la déconvolution SPOTlight, montrant la composition proportionnelle des principaux types cellulaires à chaque position pour les échantillons A1 (A) et B1 (B). (C) Répartition spatiale représentative des cellules B dans les échantillons A1 (à gauche) et B1 (à droite), illustrant les profils de localisation spatialement résolus d'une population spécifique de cellules immunitaires identifiée par déconvolution. Veuillez cliquer ici pour visualiser une version agrandie de cette figure.

Figure 5 : Analyse interactive de régions d’intérêt et comparaison de la déconvolution non supervisée. (A) Interface de l'outil personnalisé « Sélectionner les points spatiaux » montrant la sélection interactive de régions correspondant au côlon proximal, au côlon distal et à d'autres domaines tissulaires. (B) Visualisation spatiale en camembert des résultats de déconvolution non supervisée (STdeconvolve) pour l'échantillon A1, avec des points colorés selon les régions annotées manuellement sur (A), illustrant la correspondance entre l'annotation histologique et les distributions de thèmes cellulaires dérivées par calcul. Veuillez cliquer ici pour visualiser une version agrandie de cette figure.

Figure 6 : Réseaux de communication entre cellules informés par la localisation spatiale. (A,B) Réseaux d'interactions ligand-récepteur inférés par Giotto pour les échantillons A1 (A) et B1 (B). Les nœuds représentent les types cellulaires, les arêtes représentent des paires ligand-récepteur significatives (FDR < 0,05), et l'épaisseur des arêtes correspond à l'intensité de l'interaction. Afin d'assurer la comparabilité et la clarté de la visualisation, un seuil de significativité uniforme (FDR < 0,05) a été appliqué à tous les échantillons, et les 20 premières interactions classées par log2FC sont affichées pour chaque condition. Ces réseaux mettent en évidence des profils de communication spécifiques aux types cellulaires dans le contexte spatial du tissu colique. Veuillez cliquer ici pour visualiser une version agrandie de cette figure.

Figure supplémentaire 1 : Évaluation quantitative de l'optimisation des paramètres pour la réduction de dimensionnalité. Le graphique du coude illustre l'approche algorithmique du flux de travail pour sélectionner dynamiquement le nombre optimal de composantes principales (PCs). Ce choix est calculé à partir des seuils d'écart-type cumulé et de variance marginale, représentés par la ligne verticale rouge, afin de capturer la variance biologique tout en atténuant le bruit technique avant le regroupement en amont.Veuillez cliquer ici pour télécharger ce fichier.

Figure supplémentaire 2 : Validation du regroupement spatial à l'aide de marqueurs spécifiques aux couches coliques canoniques. (A) Graphique en points montrant l'expression enrichie de marqueurs épithéliaux, stromaux et musculaires lisses à travers les groupes obtenus par calcul informatique. (B) Cartes spatiales représentant la localisation de marqueurs caractéristiques (Epcam, Col1a1, Acta2) sur les coordonnées tissulaires.Veuillez cliquer ici pour télécharger ce fichier.

Discussion

Ce protocole fournit un flux de travail computationnel complet pour l'analyse de données de transcriptomique spatiale, équilibrant profondeur analytique et accessibilité pratique. L'approche étape par étape accompagne les chercheurs tout au long du pipeline analytique complet, de l'acquisition initiale des données jusqu'aux analyses spatiales avancées, en mettant l'accent sur les points de décision cruciaux et les pièges potentiels.

Plusieurs étapes du protocole méritent une attention particulière en raison de leur impact sur les résultats ultérieurs. Les seuils de contrôle qualité et de filtrage doivent être soigneusement adaptés aux types de tissus et aux plates-formes technologiques spécifiques, car un filtrage trop strict pourrait éliminer des spots biologiquement pertinents, tandis que des seuils trop permissifs peuvent introduire un bruit technique. Le choix de la méthode de normalisation influence considérablement les résultats de regroupement et d'expression différentielle en aval. Par exemple, le pipeline utilise SCTransform plutôt que la normalisation logarithmique standard, car SCTransform a été signalé pour modéliser et réduire la variation technique liée à la profondeur de séquençage et à d'autres effets techniques dans les flux de travail unicellulaires et spatiaux^14,15. Lors de l'intégration, la sélection des caractéristiques d'intégration et des paramètres de résolution nécessite une réflexion attentive afin d'équilibrer la force d'intégration et la préservation du signal biologique. Pour la détection des gènes spatialement variables, l'indice de Moran a été choisi pour sa scalabilité computationnelle et son adéquation à l'analyse d'autocorrélation spatiale dans les données basées sur des puces, bien que d'autres méthodes de détection alternatives telles que SPARK puissent être envisagées selon la taille du jeu de données et les objectifs de l'étude¹⁶.

Le composant de dépannage de ce flux de travail se concentre sur les problèmes courants d'interopérabilité qui surviennent lors du transfert des sorties d'un package à un autre. Plutôt que de s'appuyer sur une seule étape de conversion générique, le flux de travail utilise des conversions de format personnalisées à chaque interface logicielle : les données de comptage spatiales sont copiées dans un emplacement d'essai RNA standard avant l'intégration avec Seurat ; les objets de référence Seurat et monocellulaires sont convertis en objets SingleCellExperiment pour SPOTlight ; les matrices de comptage spatiales sont reformulées pour STdeconvolve ; les comptages, coordonnées et métadonnées de type cellulaire provenant de Seurat sont convertis en un objet Giotto pour l'analyse de la communication intercellulaire ; et les annotations de régions d'intérêt (ROI) provenant de Select Spatial Spots sont exportées sous forme de fichiers CSV contenant les identifiants CELL_ID, les coordonnées X/Y et les étiquettes de groupe, afin de pouvoir être replacées sur l'objet Seurat. Ces étapes aident les utilisateurs à identifier et corriger des problèmes fréquents tels que des emplacements d'essai incompatibles, des identifiants de spots non concordants, des colonnes de métadonnées manquantes, un formatage incorrect des coordonnées et des discordances entre les symboles génétiques de ligands et de récepteurs.

Une caractéristique supplémentaire de ce flux de travail est l'utilisation de deux stratégies complémentaires de déconvolution plutôt que de s'appuyer sur un seul algorithme. SPOTlight utilise des références préexistantes issues de séquençage de l'ARNm au niveau cellulaire unique (scRNA-seq) pour une estimation guidée par référence des proportions de types cellulaires, tandis que STdeconvolve permet la découverte sans référence de thèmes transcriptionnels latents. Dans le jeu de données de démonstration, les thèmes identifiés par STdeconvolve correspondaient à des régions définies histologiquement (Figure 5B), ce qui soutient leur interprétabilité biologique. Une comparaison quantitative directe entre les deux méthodes n'a pas été réalisée dans ce protocole ; les utilisateurs souhaitant les évaluer sur leurs propres données peuvent suivre le cadre fourni à l'étape 7.3.2 du protocole. En outre, le protocole favorise l'utilisabilité et la traçabilité en utilisant des scripts pour les fonctions principales tout en proposant une interface graphique (Sélectionner les Spots Spatiaux) pour une isolation intuitive des régions d'intérêt (ROI).

La liaison entre Seurat et Giotto permet au flux de travail de combiner un prétraitement et une intégration basés sur Seurat avec des analyses spatiales statistiques et des réseaux basés sur Giotto. Seurat offre un environnement éprouvé pour les chercheurs expérimentés en scRNA-seq et prend en charge l'intégration de plusieurs échantillons, tandis que le transfert des données harmonisées vers Giotto permet la construction de réseaux spatiaux et l'analyse des interactions ligand-récepteur. Cette conception permet aux utilisateurs de tirer parti des fonctionnalités documentées des deux plateformes, sans impliquer qu'elles aient été comparées ici à Squidpy ou à d'autres cadres. Pour la communication entre cellules, l'inférence de ligand-récepteur basée sur le réseau de Delaunay de Giotto a été utilisée, car elle intègre l'adjacence spatiale au sein du même environnement d'analyse ; des cadres tels que CellChat offrent des bases de données de signalisation étendues, mais ne sont pas évalués dans ce protocole¹⁷.

Plusieurs limites doivent être prises en compte. L'approche de déconvolution basée sur une référence dépend fortement de la disponibilité de références unicellulaires de haute qualité et appariées. De plus, les hypothèses fondamentales sous-jacentes à l'analyse des communications entre cellules reposent sur la co-expression des transcrits comme proxy d'une interaction physique entre protéines, ce qui nécessite une validation expérimentale ultérieure. Enfin, ce flux de travail a été principalement conçu pour des jeux de données Visium standards. À mesure que le domaine progresse vers des technologies approchant la résolution unicellulaire, telles que Visium HD, les considérations analytiques évolueront ; des données à plus haute résolution pourraient nécessiter des paramètres de prétraitement différents et réduire l'obligation absolue de déconvolution des spots. Pour adapter ce flux de travail aux jeux de données Visium HD, les bins à haute résolution peuvent être regroupés informatiquement en bins spatiaux plus larges, ou les modules de déconvolution des spots peuvent être contournés au profit d'une segmentation cellulaire basée sur l'image¹⁸.

Ces caractéristiques suggèrent une applicabilité potentielle dans divers domaines biologiques, notamment la biologie du développement, les neurosciences, la recherche sur le cancer et l'immunologie^19,20,21,22. La conception modulaire permet aux chercheurs d'adapter des composants spécifiques à leurs besoins, qu'ils se concentrent sur l'identification de domaines spatiaux, la communication cellulaire ou la spécialisation régionale. Alors que les technologies spatiales continuent d'évoluer et que les ensembles de données s'élargissent, ce protocole fournit une base pouvant être étendue pour intégrer de nouvelles méthodes d'analyse et répondre à des questions biologiques émergentes.

Matériaux

Liste des matériaux utilisés dans cet article

Nom	Entreprise	Numéro de catalogue	Commentaires
ggplot2	Posit Software, PBC	v4.0.0(CRAN)	Visualisation avancée des données
Giotto	Dries Lab	v4.2.2 (GitHub)	Analyse des réseaux spatiaux et des communications cellule-à-cellule
patchwork	Thomas Lin Pedersen	v1.3.2 (CRAN)	Composition et disposition des graphiques
R software	Fondation R pour l'informatique statistique	v4.4.3	Environnement d'exécution principal (macOS aarch64)
scater	Davis McCarthy et al.	v1.34.1 (Bioconductor)	Contrôle qualité et visualisation en cellule unique
scran	Aaron Lun et al.	v1.34.0 (Bioconductor)	Modélisation de la variance en cellule unique et détection de marqueurs
Select Spatial Spots (Outil Python personnalisé)	LeafLight	v1.0.0 (GitHub)	Sélection interactive de régions d'intérêt spatiales (ROI) (https://github.com/LeafLight/SelectSpatialSpots)
Seurat	Satija Lab	v5.3.0 (CRAN)	Prétraitement, intégration et regroupement des données spatiales
SeuratObject	Satija Lab	v5.2.0 (CRAN)	Structures de données pour les données en cellule unique et spatiales
SingleCellExperiment	Équipe centrale de Bioconductor	v1.28.1 (Bioconductor)	Conteneur de données standardisé pour la scRNA-séq
SPOTlight	Marc Elosua-Bayes et al.	v1.10.0 (Bioconductor)	Déconvolution spatiale guidée par référence
Stdeconvolve	Jean Fan Lab	v1.3.2 (Bioconductor)	Modélisation de sujets latents non supervisée
tidyverse	Posit Software, PBC	v2.0.0 (CRAN)	Suite principale de manipulation et de formatage des données