Praktyczny schemat analizy danych z transkryptomiki przestrzennej: od pozyskiwania danych po analizy zaawansowane

Podsumowanie

Niniejszy protokół przedstawia powtarzalny schemat analizy danych transkryptomiki przestrzennej, prowadząc użytkownika od pozyskania danych publicznych i kontroli jakości w oparciu o pakiet Seurat, poprzez integrację, detekcję cech przestrzennych, dekonwolucję typów komórek, adnotację obszarów zainteresowania, aż po analizę komunikacji międzykomórkowej, z uwzględnieniem praktycznych punktów kontrolnych wspierających przejrzystość wykonania.

Streszczenie

Transkryptomika przestrzenna (ST) profiluje ekspresję genów w całym genomie, zachowując dwuwymiarowy kontekst przestrzenny cząsteczek mRNA w przekrojach tkanek, co umożliwia badania architektury tkanek oraz biologii związanej z mikrośrodowiskiem. Analiza ST pozostaje jednak wyzwaniem, ponieważ import danych, kontrola jakości, integracja, dekonwolucja, statystyka przestrzenna i wizualizacja często wymagają wielu środowisk programistycznych i powtarzalnego doboru parametrów. Niniejszy protokół przedstawia praktyczny przepływ pracy obliczeniowej dla publicznych zbiorów danych ST w języku R, rozpoczynający się od pozyskania danych i konfiguracji oprogramowania, a następnie obejmujący ładowanie danych w pakiecie Seurat, kontrolę jakości, normalizację, integrację wielu próbek, klastrowanie oraz analizę genów zmiennych przestrzennie. Następnie przepływ pracy stosuje uzupełniające strategie dekonwolucji, w tym analizę SPOTlight z prowadzeniem referencyjnym oraz nienadzorowane modelowanie tematów STdeconvolve, a po nich analizę przestrzennej komunikacji międzykomórkowej w pakiecie Giotto oraz interaktywny wybór regionów zainteresowania (ROI) przy użyciu niestandardowej aplikacji Python Dash. Kładąc nacisk na wykonywanie za pomocą skryptów, jawne uzasadnienia parametrów, oczekiwane wyniki i punkty kontrolne rozwiązywania problemów, protokół zapewnia adaptacyjną strukturę dla standardowych zbiorów danych ST opartych na macierzach oraz powiązanych platform po ocenie parametrów specyficznych dla zbioru danych i platformy.

Wprowadzenie

Transkryptomika przestrzenna (ST) to przełomowa rodzina technologii umożliwiająca pomiar ekspresji genów w całym genomie przy jednoczesnym zachowaniu współrzędnych przestrzennych cząsteczek mRNA (messenger RNA) w przekrojach tkanek. Metody ST obejmują podejścia oparte na sekwencjonowaniu, które wykorzystują macierze z kodami kreskowymi określającymi pozycję, oraz podejścia obrazowe in situ, które mapują sygnały transkrypcyjne w nienaruszonych mikrośrodowiskach tkankowych^1,2. Dzięki zachowaniu kontekstu przestrzennego ST umożliwia analizę architektury tkanek, organizacji sąsiedztwa komórkowego, komunikacji międzykomórkowej oraz procesów biologicznych związanych z mikrośrodowiskiem, których nie można w pełni rozstrzygnąć po dysocjacji tkanek³.

Szybki wzrost publicznych repozytoriów danych ST stworzył bezprecedensowe możliwości dla wtórnych analiz i rozwoju nowych metod³. Zasoby takie jak baza danych CROST katalogują setki zestawów danych transkryptomicznych o rozdzielczości przestrzennej u wielu gatunków i na różnych platformach technologicznych, podczas gdy specjalistyczne kolekcje, takie jak STOmicsDB, koncentrują się na konkretnych metodologiach, takich jak Stereo-seq^4,5. Pomimo tej obfitości danych, analiza obliczeniowa pozostaje wyzwaniem ze względu na złożoność struktur danych przestrzennych, różnorodność narzędzi analitycznych oraz techniczne przeszkody w implementacji powtarzalnych przepływów pracy^{6,7,8,9,10,11}.

Aby rozwiązać ograniczenia wynikające z polegania na jednym środowisku programistycznym, w niniejszym opracowaniu przedstawiono zintegrowany przepływ obliczeniowy wykorzystujący komplementarne narzędzia analityczne. Istniejące kompleksowe ekosystemy analizy ST obejmują przede wszystkim Seurat, Giotto oraz frameworki oparte na języku Python, takie jak Squidpy^6,7,12. Chociaż narzędzia oparte na Pythonie, takie jak Squidpy, oferują szerokie funkcjonalności w zakresie analizy grafów przestrzennych, skonsolidowanie głównej ścieżki obliczeniowej w obrębie jednego języka programowania minimalizuje techniczne trudności związane z wymianą danych między językami. W związku z tym główny proces został zaimplementowany głównie w języku R. W ramach tego przepływu pracy opartego na R, Seurat jest wykorzystywany do ładowania danych, kontroli jakości, normalizacji, redukcji wymiarowości, wizualizacji oraz integracji wielu próbek, co odzwierciedla jego powszechne zastosowanie w analizach transkryptomiki jednokomórkowej i przestrzennej. Następnie Giotto jest wykorzystywane do konstrukcji sieci przestrzennych oraz analizy komunikacji międzykomórkowej w oparciu o interakcje ligand-receptor. Tym samym potok ten łączy wstępne przetwarzanie i integrację w Seurat z analizą przestrzenną w Giotto, zapewniając przy tym jawny i powtarzalny transfer danych między oboma zestawami narzędzi.

W ramach tego podejścia zaimplementowano dwie komplementarne strategie dekonwolucji: SPOTlight, metodę z wykorzystaniem referencji, która używa danych scRNA-seq do szacowania proporcji typów komórek, oraz STdeconvolve, nienadzorowaną metodę modelowania tematów, która identyfikuje utajone wzorce transkrypcyjne^8,11. Wyniki dostarczają komplementarnych spojrzeń na przestrzenną heterogeniczność komórkową, jednak nie są traktowane jako ilościowa walidacja krzyżowa, chyba że użytkownicy przeprowadzą opcjonalną analizę zgodności opisaną w protokole. Zintegrowano dedykowaną aplikację Python Dash, Select Spatial Spots, służącą do interaktywnej adnotacji obszarów zainteresowania (ROI), która eksportuje standardowe pliki adnotacji oparte na współrzędnych, mogące być wykorzystane w późniejszych analizach downstream.

W odniesieniu do praktycznego zastosowania, niniejszy schemat postępowania jest przeznaczony przede wszystkim dla standardowych danych ST opartych na matrycach (np. Visium o rozdzielczości 55 µm) i może być dostosowany do innych typów tkanek po ewaluacji parametrów. Przed przystąpieniem do analizy należy rozważyć kluczowe ograniczenia. Po pierwsze, moduł dekonwolucji kierowanej referencyjnie zależy od wysokiej jakości referencji scRNA-seq, dopasowanej do danej tkanki. Po drugie, platformy o rozdzielczości subkomórkowej lub zbliżonej do pojedynczej komórki mogą wymagać zmodyfikowanego wstępnego przetwarzania, agregacji binów przestrzennych lub segmentacji komórek na podstawie obrazów przed integracją². Reprezentatywny zbiór danych z jelita grubego myszy służy jako przykład demonstracyjny, aby pokazać, w jaki sposób schemat postępowania może oceniać domeny przestrzenne i organizację tkanki zdefiniowaną przez markery, a nie jako dowód na uniwersalną kompatybilność platformy.

Protokół

Wszystkie zbiory danych biologicznych analizowane w niniejszym protokole są ogólnodostępne i wykorzystane wyłącznie w celach demonstracyjnych. Konkretne numery dostępowe do danych oraz repozytoria źródłowe podano w odpowiednich krokach. Oryginalne zbiory danych zostały wygenerowane przez pierwotnych badaczy zgodnie z instytucjonalnymi wytycznymi etycznymi obowiązującymi dla każdego badania źródłowego. Należy zapoznać się z Tabelą Materiałów, aby zweryfikować wszystkie wymagane wersje oprogramowania i pakietów R.

Wymagania sprzętowe: Pamięć obliczeniowa wymagana dla tego przepływu pracy skaluje się wraz z liczbą analizowanych próbek i punktów. W przypadku typowego zbioru danych transkryptomiki przestrzennej (np. około 3 000 punktów na próbkę dla maksymalnie trzech próbek), standardna stacja robocza z minimum 16 GB pamięci RAM jest wystarczająca do wykonania potoku przetwarzania. Jednakże, w celu zapewnienia optymalnej wydajności i stabilności, szczególnie podczas pamięciożernych etapów obliczeniowych, takich jak normalizacja SCTransform oraz faktoryzacja macierzy podczas dekonwolucji, wysoce zalecane jest 32 GB pamięci RAM lub więcej.

1. Pozyskiwanie danych i przygotowanie struktury katalogów

1. Pobrać publiczne zbiory danych transkryptomiki przestrzennej (np. numer dostępu GEO GSE169749, data przesłania 26 marca 2021 r. i data ostatniej aktualizacji 6 marca 2022 r.) oraz referencyjne dane scRNA-seq (np. numer dostępu GEO GSE264408, data przesłania 19 kwietnia 2024 r. i data ostatniej aktualizacji 10 grudnia 2024 r.).
2. Pobrać archiwum z surowymi danymi (np. GSE169749_RAW.tar) z repozytorium i rozpakować archiwum.
   UWAGA: W przypadku reprezentatywnego zbioru danych użytego w niniejszym protokole, rozpakowany katalog zawiera wiele plików dla każdej próbki, a nazwy plików zazwyczaj zawierają numer dostępu do próbki GEO (np. GSM5213483).
3. Zorganizować pliki w standaryzowanej strukturze katalogów kompatybilnej z funkcją Seurat Load10X_Spatial.
   1. Utworzyć główny katalog danych (np. ./data/).
   2. Utworzyć dedykowany podkatalog dla każdej próbki (np. ./data/sample_A1/).
   3. Przenieść (skopiować lub przenieść) następujące niezbędne pliki dla każdej próbki do odpowiedniego podkatalogu.
      1. Umieścić plik filtered_feature_bc_matrix.h5 w głównym podkatalogu próbki.
      2. Utworzyć podfolder spatial/ w każdym katalogu próbki.
      3. Umieścić następujące pliki w podfolderze spatial/: tissue_positions_list.csv.gz, scalefactors_json.json.gz oraz tissue_hires_image.png.gz.
   4. Rozpakować pliki .gz w folderze spatial/. Upewnić się, że oryginalne nazwy plików pozostały dokładnie takie, jakich wymaga funkcja Load10X_Spatial (np. tissue_positions_list.csv).
      UWAGA: Końcowa, zorganizowana struktura dla sprawnego ładowania musi odpowiadać temu wzorowi, przyjmując próbkę A1 jako przykład:
      ./data/sample_A1/filtered_feature_bc_matrix.h5
      ./data/sample_A1/spatial/tissue_positions_list.csv
      ./data/sample_A1/spatial/scalefactors_json.json
      ./data/sample_A1/spatial/tissue_hires_image.png

2. Konfiguracja środowiska programistycznego

1. Zainstaluj program R (wersja 4.4.3 lub wyższa).
2. Pobierz skrypty analityczne (od 1_ReadSpatialData.R do 8_Giotto_Communication.R) z repozytorium GitHub (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow, suma kontrolna commitu 2d85e18 w celu zapewnienia powtarzalności).
3. Zainstaluj wymagane pakiety R z CRAN i Bioconductor, uruchamiając skrypt setup.R.
   UWAGA: Pełna lista konkretnych wersji pakietów wykorzystanych w tym przepływie pracy znajduje się w Tabeli Materiałów oraz w pliku session_info.txt w powiązanym repozytorium GitHub (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow).
4. Zainstaluj Giotto Suite, wykonując polecenia instalacyjne zawarte w oficjalnej dokumentacji.
   1. Przejdź do oficjalnego adresu URL instalacji, aby pobrać skrypty konfiguracyjne: https://giottosuite.com/articles/installation.html
   2. Zainicjuj wymagane środowisko Python oraz zależności systemowe zgodnie z instrukcjami krok po kroku.
5. Zainstaluj dodatkowe pakiety w celu ulepszonej wizualizacji.
   1. Wykonaj: install.packages(c("ggprism", "colorBlindness"))
6. Zainstaluj niestandardowe narzędzie Select Spatial Spots. Upewnij się, że system operacyjny to Windows, macOS lub Linux, a zainstalowana jest wersja języka Python 3.8 lub wyższa. Krok ten jest opcjonalny.
   1. Pobierz narzędzie, przechodząc do repozytorium GitHub (https://github.com/LeafLight/SelectSpatialSpots, suma kontrolna commitu d20946e w celu zapewnienia powtarzalności) i pobierając kod źródłowy.
   2. Przejdź do katalogu narzędzia i zainstaluj zależności Pythona: pip install -r requirements.txt
7. Zainstaluj pakiet glmGamPoi , wykonując polecenie BiocManager::install("glmGamPoi"), aby przyspieszyć normalizację SCTransform.
   UWAGA: Podstawowy przepływ pracy wymaga jedynie pakietów wymienionych w krokach 2.1–2.3 i 2.6. Pakiet Rfast2 służy do przyspieszenia obliczeń statystyki I Morana. Opcjonalne pakiety z kroku 2.5 służą do generowania wykresów gotowych do publikacji w stylu prism (ggprism) oraz do korzystania z palety kolorów dostosowanej do osób z zaburzeniami rozpoznawania barw (colorBlindness). Narzędzie z kroku 2.6 umożliwia interaktywny wybór punktów (spotów). Wykonaj polecenie sessionInfo() i zapisz pełny wynik konsoli w pliku tekstowym. Dokumentuje to dokładne wersje oprogramowania i zależności pakietów, co jest kluczowe dla zapewnienia długoterminowej powtarzalności.

3. Ładowanie danych przestrzennych i kontrola jakości (1_ReadSpatialData.R, 2_SpatialDataQC.R)

1. Wczytaj dane przestrzenne do obiektów Seurat.
   1. Użyj Read10X_Image, aby ręcznie wczytać obraz tkanki o wysokiej rozdzielczości, określając image.dir i image.name.
   2. Użyj funkcji Load10X_Spatial z parametrem image ustawionym na obiekt obrazu utworzony w kroku 3.1.1, aby stworzyć obiekt Seurat.
2. Oblicz wskaźniki kontroli jakości.
   1. Oblicz procent odczytów mitochondrialnych za pomocą funkcji PercentageFeatureSet z wzorcem ^mt-.
3. Wizualizuj i interpretuj dane na podstawie wskaźników QC.
   1. Wygeneruj wykresy skrzypcowe (violin plots) dla nCount_Spatial, nFeature_Spatial oraz percent.mt za pomocą funkcji VlnPlot.
   2. Utwórz przestrzenne wykresy cech (spatial feature plots) tych wskaźników za pomocą SpatialFeaturePlot, aby zidentyfikować punkty znajdujące się poza obszarem tkanki.
      UWAGA: Punkt kontrolny wizualizacji: Punkty poza obszarem tkanki zazwyczaj wykazują niską liczbę UMI (nCount_Spatial < 500) i niską wykrywalność genów (nFeature_Spatial < 200)
   3. W celach demonstracyjnych zastosuj filtry w celu usunięcia punktów o niskiej jakości (np. subset(seurat_obj, subset = nFeature_Spatial > 200 & nCount_Spatial > 500)). Jest to krok opcjonalny.
      UWAGA: Głównym celem QC w transkryptomice przestrzennej jest identyfikacja i adnotacja artefaktów technicznych, takich jak punkty poza tkanką. W przypadku demonstracyjnych danych Visium z jelita grubego, warunki nFeature_Spatial > 200 i nCount_Spatial > 500 usuwają punkty o niskiej złożoności/tła. W przypadku tkanek o niskiej zawartości RNA lub zdegradowanych przekrojów, sprawdź wykresy skrzypcowe i przestrzenne wykresy cech przed podniesieniem progów; dla tkanek gęstych lub o wysokim poziomie UMI mogą być odpowiednie bardziej rygorystyczne progi. Unikaj filtrowania wyłącznie na podstawie wysokich wartości ekspresji, chyba że widoczne są dublety, fałdy tkanki lub wyraźne artefakty. Filtrowanie jest często niewskazane, ponieważ fizycznie usuwa lokalizacje przestrzenne, co może zakłócić ciągłość struktury tkanki w dalszej analizie przestrzennej.

4. Wstępne przetwarzanie, integracja i klasteryzacja danych (3_IntegrationAndClustering.R)

1. Znormalizuj i przeprowadź wstępną obróbkę poszczególnych próbek.
   1. Zastosuj normalizację SCTransform do każdej próbki oddzielnie, ustawiając assay = Spatial.
2. Zintegruj wiele próbek.
   1. Przygotuj listę obiektów znormalizowanych metodą SCT do integracji. Upewnij się, że każdy obiekt posiada assay „RNA”, kopiując assay „Spatial”: spatial_list[[1]][["RNA"]] <- spatial_list[[1]][["Spatial"]]
      ​UWAGA: Kopiowanie assay „Spatial” do standardowego slotu assay „RNA” służy jako niezbędne obejście w celu zapewnienia kompatybilności z funkcjami integracji w dalszych etapach, które zostały pierwotnie zaprojektowane dla obiektów Seurat z danych jednokomórkowych.
   2. Użyj funkcji SelectIntegrationFeatures oraz PrepSCTIntegration, aby zidentyfikować wspólne cechy zmienne.
   3. Znajdź kotwice integracji (integration anchors) za pomocą funkcji FindIntegrationAnchors z parametrem normalization.method = "SCT".
   4. Zintegruj dane za pomocą funkcji IntegrateData z parametrem normalization.method = "SCT".
3. Przeprowadź redukcję wymiarowości i klasteryzację na zintegrowanym assay.
   1. Uruchom PCA na zintegrowanych danych, korzystając z funkcji RunPCA.
   2. Wyznacz optymalną liczbę głównych składowych (PCs) do dalszych analiz poprzez obliczenie skumulowanej wariancji. Zidentyfikuj programowo punkt przegięcia (np. punkt, w którym wariancja skumulowana przekracza 90%, a przyrost marginalny spada poniżej 0,1%), korzystając z kodu w pliku 3_IntegrationAndClustering.R, linie 36-38. Otrzymaną liczbę PCs określaj dalej jako pc.use.
   3. Uruchom funkcję RunUMAP z parametrem dims = 1:pc.use.
   4. Pogrupuj komórki w klastry, używając funkcji FindNeighbors z parametrem dims = 1:pc.use oraz funkcji FindClusters z parametrem resolution = 0.5. Dostosuj rozdzielczość (resolution) dopiero po sprawdzeniu stabilności klastrów i spójności genów markerowych.
   5. Przeprowadź analizę różnicowej ekspresji między grupami docelowymi (np. „B1_colon_d14” versus „A1_colon_d0”) za pomocą funkcji FindMarkers.
      ​UWAGA: Punkt kontrolny wizualny: Prawidłowa integracja powinna wykazać odpowiednie wymieszanie próbek na wykresie UMAP przy jednoczesnym zachowaniu biologicznie odrębnych klastrów.
4. Zidentyfikuj geny zmienne przestrzennie.
   1. Dla każdej oryginalnej próbki uruchom funkcję FindSpatiallyVariableFeatures, stosując metodę „moransi” na assay „SCT”, aby obliczyć autokorelację przestrzenną.

5. Przetwarzanie wstępne danych referencyjnych z pojedynczych komórek (4_scDataPreProcessing.R)

1. Wczytaj macierz zliczeń scRNA-seq za pomocą funkcji Read10X i utwórz obiekt Seurat.
2. Przeprowadź standardową kontrolę jakości (QC), normalizację i klastrowanie.
   1. Oblicz procent odczytów mitochondrialnych (percent.mt) i odfiltruj komórki (np. nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 25).
   2. Znormalizuj dane za pomocą funkcji SCTransform, ustawiając vars.to.regress = "percent.mt".
   3. Uruchom PCA, UMAP i przeprowadź klastrowanie komórek, stosując metodę dynamicznego wyboru PC opisaną w kroku 4.3.2.
3. Wykonaj adnotację typów komórek.
   1. Oblicz wyniki modułów (module scores) dla kanonicznych genów markerowych typów komórek za pomocą funkcji AddModuleScore.
   2. Zadnotuj klastry na podstawie wyników modułów i znanej biologii. Alternatywnie zaimportuj wcześniej obliczone adnotacje z metadanych.

6. Dekonwolucja z wykorzystaniem referencji w programie SPOTlight (5_SPOTlight_Deconv.R)

1. Przygotowanie danych dla SPOTlight.
   1. Konwersja adnotowanego obiektu Seurat z danych jednokomórkowych oraz obiektu Seurat z danych przestrzennych na obiekty SingleCellExperiment.
   2. Log-normalizacja danych jednokomórkowych przy użyciu logNormCounts.
2. Uruchomienie dekonwolucji SPOTlight.
   1. Identyfikacja genów hiperzmiennych (HVGs) w danych jednokomórkowych przy użyciu funkcji modelGeneVar i getTopHVGs.
   2. Obliczenie genów markerowych dla typów komórek przy użyciu scoreMarkers i odfiltrowanie markerów wysokiej jakości (np. mean.AUC > 0,8).
   3. Downsampling referencji jednokomórkowej dla każdego typu komórek do zarządzalnej liczby (np. 50 komórek) w celu skrócenia czasu obliczeń.
   4. Wykonanie dekonwolucji przy użyciu funkcji SPOTlight z parametrami weight_id = "mean.AUC", group_id = "cluster" oraz gene_id = "gene".
3. Wizualizacja i eksport wyników.
   1. Ekstrakcja macierzy wyników dekonwolucji (proporcje typów komórek na spot).
   2. Użycie funkcji plotSpatialScatterpie do wizualizacji składu komórkowego w współrzędnych przestrzennych.
   3. Dodanie wyników dekonwolucji do metadanych przestrzennego obiektu Seurat przy użyciu AddMetaData. Oczekiwany wynik: macierz dekonwolucji z kolumnami proporcji typów komórek dla każdego spotu przestrzennego, wykresy spatial scatterpie pokazujące lokalny skład komórkowy oraz obiekt Seurat zawierający proporcje dekonwolucji w metadanych.

7. Dekonwolucja bezreferencyjna za pomocą STdeconvolve (7_STdeconvolve.R)

1. Przygotowanie danych przestrzennych.
   1. Wyodrębnij macierz surowych zliczeń z przestrzennego obiektu Seurat, korzystając z funkcji GetAssayData z parametrem slot = "counts".
   2. Usuń niskiej jakości punkty (spots) i geny za pomocą funkcji cleanCounts z pakietu STdeconvolve (np. min.lib.size = 100).
2. Identyfikacja ukrytych typów komórek.
   1. Dalej filtruj korpus pod kątem genów ekspresjonowanych w minimalnej frakcji punktów, korzystając z funkcji restrictCorpus (np. removeAbove=1.0, removeBelow = 0.05).
   2. Dopasuj model Latent Dirichlet Allocation (LDA) dla zakresu potencjalnych liczb tematów (K) (np. Ks = seq(2, 9, by = 1)), używając funkcji fitLDA.
   3. Wybierz optymalny model na podstawie minimalnej perpleksji, korzystając z funkcji optimalModel z parametrem opt = "min".
3. Analiza i wizualizacja wyników.
   1. Wyodrębnij proporcje typów komórek (tematów) (theta) oraz profile genów (beta) z optymalnego modelu za pomocą funkcji getBetaTheta.
   2. Aby wspomóc biologiczną interpretację zdekonwolwowanych tematów, zaimportuj adnotacje obszarów zainteresowania (ROI) wygenerowane przez narzędzie „Select Spatial Spots” (patrz krok 9.3.3). Użyj tych adnotacji jako parametru groups w funkcji vizAllTopics, aby rzutować zdekonwolwowane proporcje typów komórek z powrotem na współrzędne przestrzenne i zakodować kolorami punkty według ich ROI. Jest to krok opcjonalny.
      UWAGA: Oczekiwany wynik: STdeconvolve powinno zwrócić wartości theta reprezentujące proporcje tematów na punkt oraz wartości beta reprezentujące profile genów dla każdego tematu. Jeśli użytkownicy chcą ilościowo porównać wyniki nienadzorowanego STdeconvolve z wynikami SPOTlight kierowanymi przez referencje (z Kroku 6), mogą wyeksportować obie macierze proporcji (theta z STdeconvolve oraz macierz SPOTlight) i obliczyć korelacje na poziomie punktów lub metryki zgodności (np. korelację Pearsona lub Spearmana) przy użyciu standardowych funkcji R. Ilościowe porównanie SPOTlight i STdeconvolve jest opcjonalne i nie jest wymagane do wykonania podstawowego przepływu pracy; jednak odpowiednie macierze wyników są dostępne w Krokach 6.3.1 i 7.3.1 dla użytkowników, którzy chcą przeprowadzić taką analizę.

8. Przestrzenna komunikacja międzykomórkowa z wykorzystaniem Giotto(8_Giotto_Communication.R)

1. Konwertuj obiekt Seurat na obiekt Giotto.
   1. Użyj funkcji createGiottoObject, podając macierz surowych zliczeń oraz współrzędne przestrzenne.
2. Przeprowadź wstępne przetwarzanie obiektu Giotto i dodaj wyniki dekonwolucji.
   1. Znormalizuj dane za pomocą funkcji normalizeGiotto.
   2. Dodaj adnotacje typów komórek (np. główny typ komórki z dekonwolucji SPOTlight) do metadanych komórkowych za pomocą funkcji addCellMetadata.
3. Wywnioskuj komunikację międzykomórkową z uwzględnieniem informacji przestrzennych.
   1. Utwórz sieć przestrzenną za pomocą funkcji createSpatialNetwork z parametrami method = "Delaunay" oraz name = "Delaunay_network".
   2. Załaduj bazę danych ligand-receptor do środowiska R. Upewnij się, że załadowany obiekt ma format ramki danych (data frame) z dwiema odrębnymi kolumnami reprezentującymi symbole genów ligandu i receptora (np. wstępnie skompilowana sieć dla myszy z Zenodo: https://zenodo.org/api/records/15168114/files/lr_network_mouse.csv/content, wersja v6, opublikowana 7 kwietnia 2025 r.).
   3. Uruchom funkcję exprCellCellcom z parametrem cluster_column = "celltype_major", aby zidentyfikować istotne interakcje ligand-receptor między typami komórek znajdującymi się w bliskim sąsiedztwie przestrzennym. Oczekiwany wynik: tabela zawierająca pary ligand-receptor, kombinacje typów komórek źródłowych i docelowych, wartości log2 fold-change oraz skorygowane wartości P; zachowaj istotne interakcje (np. p.adj < 0.05) do wizualizacji.

9. Interaktywny wybór obszarów za pomocą funkcji select spatial spots (6_SelectSpatialSpots.R)

1. Przygotowanie danych do narzędzia interaktywnego.
   1. Wyodrębnienie współrzędnych przestrzennych z obiektu Seurat przy użyciu funkcji GetTissueCoordinates.
2. Formatowanie i eksport danych.
   1. Sformatowanie ramki danych współrzędnych tak, aby zawierała następujące nagłówki kolumn: CELL_ID, X oraz Y.
   2. Eksport sformatowanej ramki danych do pliku CSV.
3. Przeprowadzenie analizy obszaru zainteresowania (ROI).
   1. Uruchomienie dedykowanej aplikacji Dash Select Spatial Spots i wczytanie pliku CSV.
   2. Interaktywny wybór punktów na podstawie ich lokalizacji przestrzennej.
   3. Eksport listy wybranych punktów wraz z przypisanymi etykietami grup/ROI do nowego pliku CSV. Oczekiwany wynik: plik CSV z identyfikatorami punktów, współrzędnymi X/Y oraz przypisanymi etykietami grup/ROI, które można powiązać z obiektem Seurat za pomocą CELL_ID.

Wyniki

Wdrożenie przepływu pracy i integracja danych ilustrują główne cechy tkanek

Przepływ obliczeniowy zastosowano do danych z transkryptomiki przestrzennej jelita grubego myszy, aby zilustrować oczekiwane wyniki na poszczególnych etapach analitycznych. Jak przedstawiono na schemacie przepływu (Rysunek 1), potok rozpoczął się od pozyskania danych i kontroli jakości, gdzie wykresy cech przestrzennych wyznaczyły granice tkanek (Rysunek 2A,B). Następnie wykorzystano proces integracji oparty na kotwicach (anchors) w programie Seurat w celu zredukowania technicznych efektów partii przy zachowaniu interpretowalnej zmienności biologicznej. Wizualizacje UMAP po integracji wykazały wyrównanie próbek oraz wzorce klastrowania przestrzennego (Rysunek 2C,D). W celu sterowania redukcją wymiarowości i późniejszym klastrowaniem zaimplementowano ilościowy, dynamiczny wybór głównych składowych (PCs) w oparciu o wariancję skumulowaną (patrz Rysunek uzupełniający 1). Analiza mapy ciepła genów markerowych wykazała odrębne profile transkrypcyjne leżące u podstaw klastrów przestrzennych (Rysunek 2E).

Aby ocenić, czy klastry obliczeniowe były zgodne ze znaną architekturą anatomiczną histologii jelita grubego, przeanalizowano profile ekspresji kanonicznych markerowych genów specyficznych dla poszczególnych warstw. Warstwa nabłonka śluzówki wykazała ekspresję markerów komórek nabłonkowych, w tym Epcam i Krt8, wraz z markerem komórek kubkowych Muc2. Markery mezenchymalne i stropowe, takie jak Col1a1 i Vim, oznakowały obszary błony właściwej i podśluzowej, natomiast zewnętrzna warstwa mięśniówki gładkiej była wskazana przez strukturalne geny mięśni gładkich, takie jak Acta2 i Tagln. Przestrzenne ograniczenie tych markerów związanych z liniami komórkowymi potwierdza interpretację, że proces integracji i klastrowania zachował główne warstwy histologiczne tkanki jelita grubego wzdłuż osi od śluzówki do mięśniówki (patrz Supplementary Figure 2).

Po walidacji klastrów przeprowadzono analizę różnicowej ekspresji w celu zidentyfikowania genów o różnej ekspresji (DEGs) pomiędzy warunkami eksperymentalnymi (Rysunek 2F,G). Ponadto, przy użyciu statystyki Moran's I, zidentyfikowano geny o zmienności przestrzennej, wskazując geny o istotnym, nierandomowym rozkładzie przestrzennym w obrębie tkanki (Rysunek 2H).

Dekonwolucja komórkowa i przestrzenne sieci interakcji ujawniają mikroorganizację tkanki

Przetwarzanie referencyjnych danych sekwencjonowania RNA pojedynczych komórek (single-cell RNA-seq) pozwoliło uzyskać adnotacje potwierdzone poprzez filtrowanie QC (Rycina 3A), klasteryzację nienadzorowaną (Rycina 3B), walidację genów markerowych (Rycina 3C) oraz zgodność z niezależnymi adnotacjami (Rycina 3D). Skład komórkowy (Rycina 3E) posłużył do opracowania strategii downsamplingu dla dekonwolucji. Program SPOTlight oszacował proporcje typów komórek prowadzone przez referencje w poszczególnych spotach przestrzennych (Rycina 4A,B), podczas gdy STdeconvolve dostarczył nienadzorowany wgląd w przestrzenne wzorce komórkowe oparte na modelowaniu tematów (Rycina 5B). Niestandardowe narzędzie Select Spatial Spots zapewniło kontekst histologiczny dla tych wzorców (Rycina 5A). Na koniec, wykorzystując przypisania typów komórek uzyskane z dekonwolucji, analiza komunikacji przestrzennej zidentyfikowała oddziaływania ligand-receptor pomiędzy grupami typów komórek znajdujących się blisko siebie w przestrzeni (Rycina 6A,B).

Uwagi dotyczące rozwiązywania problemów wynikające z optymalizacji protokołu

Podczas optymalizacji protokołu zidentyfikowano kilka problemów, które stały się podstawą do ustalenia praktycznych punktów kontrolnych. Suboptymalne wyniki dekonwolucji występowały w sytuacjach, gdy referencje dla pojedynczych komórek były słabo dopasowane do kontekstu tkankowego, co wskazało na konieczność stosowania danych scRNA-seq dopasowanych do tkanki i gatunku, jeśli są one dostępne. Początkowe próby klasteryzacji z parametrami domyślnymi nie zawsze pozwalały na rozdzielenie oczekiwanych struktur biologicznych; analiza wyboru PC, rozdzielczości klasteryzacji oraz spójności genów markerowych pomogła zidentyfikować interpretowalne przestrzennie domeny zgodne z anatomią tkanki. Obserwacje te stanowią praktyczne przykłady tego, jak użytkownicy mogą diagnozować typowe problemy analityczne podczas realizacji schematu postępowania.

Rycina 1: Schemat zintegrowanej analizy transkryptomiki przestrzennej. Schematyczna reprezentacja potoku analitycznego, od pozyskania i wstępnego przetwarzania danych po zaawansowane analizy przestrzenne. Kluczowe etapy obejmują: (1) ładowanie danych, kontrolę jakości oraz integrację wielu próbek przy użyciu Seurat; (2) klastrowanie przestrzenne i detekcję genów o zmienności przestrzennej; (3) dekonwolucję typów komórek via metody oparte na referencji (SPOTlight) oraz metody nienadzorowane (STdeconvolve); (4) analizę przestrzennej komunikacji międzykomórkowej z wykorzystaniem Giotto oraz interaktywny wybór regionów zainteresowania za pomocą autorskiego narzędzia Select Spatial Spots. Wyniki ze wszystkich modułów są syntetyzowane w celu uzyskania biologicznych wniosków dotyczących architektury tkanki i mikrośrodowiska komórkowego. Proszę kliknąć tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 2Integracja danych, klasteryzacja oraz analiza różnicowej ekspresji. (A,B) Metryki kontroli jakości dla próbek przestrzennych A1 i B1, przedstawiające rozkłady liczby genów, liczby UMI oraz procentowej zawartości genów mitochondrialnych. (C) Wizualizacja UMAP zintegrowanych danych z transkryptomiki przestrzennej, z kolorowaniem według pochodzenia próbki (lewo) i tożsamości klastra (prawo). (D) Projekcja przestrzenna tożsamości klastrów na przekroje tkankowe. (E) Mapa ciepła głównych genów markerowych dla każdego klastra przestrzennego. (F) Wykres wulkaniczny (volcano plot) przedstawiający geny różnicowo wyrażone pomiędzy warunkami A1_colon_d0 a B1_colon_d14. (G) Wzorce ekspresji przestrzennej reprezentatywnych genów o zróżnicowanej ekspresji w przekrojach tkanek. (H) Mapy ekspresji przestrzennej zidentyfikowanych genów o najwyższej zmienności przestrzennej poprzez Statystyka I Morana, gdzie dwa lewe panele przedstawiają geny z próbki A1_colon_d0, a dwa prawe panele przedstawiają geny z próbki B1_colon_d14. Kliknij tutaj, aby wyświetlić powiększoną wersję tej figury.

Rycina 3: Przetwarzanie i adnotacja referencyjnych danych z pojedynczych komórek. (A) Metryki kontroli jakości dla referencyjnych danych scRNA-seq przed i po filtrowaniu. (B) Wizualizacja UMAP danych scRNA-seq z kolorowaniem według klastrów nienadzorowanych. (C) Wykres punktowy (dot plot) przedstawiający wyniki ekspresji kanonicznych genów markerowych typów komórek w poszczególnych klastrach. (D) Adnotowana wizualizacja UMAP danych scRNA-seq z oznaczeniami głównych typów komórek. (E) Skład komórkowy referencyjnego zbioru danych scRNA-seq. Czerwona przerywana linia wskazuje próg downsamplingu (n = 50 komórek na typ) zastosowany podczas dekonwolucji SPOTlight w celu zrównoważenia wydajności obliczeniowej i reprezentacji typów komórek. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rysunek 4: Przestrzenna dekonwolucja heterogeniczności komórkowej.(A,B) Przestrzenne wykresy typu scatterpie z dekonwolucji SPOTlight, przedstawiające proporcjonalny skład głównych typów komórek w każdym punkcie (spot) dla próbek A1 (A) i B1 (B). (C) Reprezentatywny rozkład przestrzenny limfocytów B w próbkach A1 (lewo) i B1 (prawo), obrazujący wzorce lokalizacji z rozdzielczością przestrzenną dla konkretnej populacji komórek odpornościowych zidentyfikowanej za pomocą dekonwolucji. Aby wyświetlić powiększoną wersję tego rysunku, kliknij tutaj.

Rysunek 5: Interaktywna analiza obszarów zainteresowania i porównanie dekonwolucji nienadzorowanej. (A) Interfejs niestandardowego narzędzia „Select Spatial Spots”, pokazujący interaktywny wybór obszarów odpowiadających bliższej części okrężnicy, dalszej części okrężnicy oraz innym domenom tkankowym. (B) Wizualizacja przestrzenna typu scatterpie wyników dekonwolucji nienadzorowanej (STdeconvolve) dla próbki A1, w której punkty są pokolorowane zgodnie z ręcznie naniesionymi adnotacjami obszarów z (A), co ilustruje zgodność między adnotacją opartą na histologii a obliczeniowo wyprowadzonymi rozkładami tematów komórkowych. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 6: Przestrzenne sieci komunikacji międzykomórkowej.(A,B) Sieci interakcji ligand-receptor wyznaczone za pomocą Giotto dla próbek A1 (A) i B1 (B). Węzły reprezentują typy komórek, krawędzie reprezentują istotne pary ligand-receptor (FDR < 0,05), a grubość krawędzi odpowiada sile interakcji. Aby zapewnić porównywalność i przejrzystość wizualizacji, we wszystkich próbkach zastosowano jednolity próg istotności (FDR < 0,05) i wyświetlono 20 najważniejszych interakcji uszeregowanych według log2FC dla każdego stanu. Sieci uwydatniają specyficzne dla typów komórek wzorce komunikacji w kontekście przestrzennym tkanki jelita grubego. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rysunek uzupełniający 1: Ilościowa ocena optymalizacji parametrów w celu redukcji wymiarowości.Wykres łokciowy przedstawia programistyczne podejście procesu analizy do dynamicznego wyboru optymalnej liczby głównych składowych (PCs). Wybór jest obliczany na podstawie skumulowanego odchylenia standardowego oraz progów wariancji krańcowej, reprezentowanych przez czerwoną linię pionową, aby uchwycić wariancję biologiczną przy jednoczesnym ograniczeniu szumu technicznego przed późniejszą klastryzacją.Proszę kliknąć tutaj, aby pobrać ten plik.

Rycina uzupełniająca 2: Walidacja klastrowania przestrzennego z wykorzystaniem kanonicznych markerów specyficznych dla poszczególnych warstw okrężnicy. (A) Wykres punktowy (dot plot) przedstawiający wzbogaconą ekspresję markerów nabłonkowych, zrębowych i mięśni gładkich w klastrach obliczeniowych. (B) Przestrzenne wykresy cech mapujące reprezentatywne markery (Epcam, Col1a1, Acta2) z powrotem na współrzędne tkanki.Proszę kliknąć tutaj, aby pobrać ten plik.

Dyskusja

Niniejszy protokół przedstawia kompleksowy przepływ pracy obliczeniowej dla analizy danych z transkryptomiki przestrzennej, który równoważy głębię analityczną z praktyczną dostępnością. Podejście krok po kroku prowadzi badaczy przez cały proces analityczny, od początkowego pozyskania danych po zaawansowane analizy przestrzenne, kładąc nacisk na kluczowe punkty decyzyjne oraz potencjalne pułapki.

Kilka etapów protokołu wymaga szczególnej uwagi ze względu na ich wpływ na dalsze wyniki. Progi kontroli jakości i filtrowania powinny być starannie dostosowane do konkretnych typów tkanek i platform technologicznych, ponieważ zbyt rygorystyczne filtrowanie może usunąć biologicznie istotne plamki, natomiast zbyt łagodne progi mogą wprowadzić szum techniczny. Wybór metody normalizacji znacząco wpływa na późniejszą klasteryzację i wyniki analizy ekspresji różnicowej. Na przykład w potoku przetwarzania wykorzystano SCTransform zamiast standardowej normalizacji logarytmicznej, ponieważ wykazano, że SCTransform modeluje i redukuje zmienność techniczną związaną z głębokością sekwencjonowania oraz inne efekty techniczne w analizach pojedynczych komórek i analizach przestrzennych^14,15. Podczas integracji wybór cech integracyjnych oraz parametrów rozdzielczości wymaga uważnego rozważenia, aby zrównoważyć siłę integracji z zachowaniem sygnału biologicznego. Do identyfikacji genów zmiennych przestrzennie wybrano wskaźnik Moran's I ze względu na skalowalność obliczeniową i przydatność do analizy autokorelacji przestrzennej w danych opartych na macierzach, choć w zależności od wielkości zbioru danych i celów badania można rozważyć alternatywne narzędzia, takie jak SPARK¹⁶.

Komponent rozwiązywania problemów w ramach tego przepływu pracy koncentruje się na powszechnych problemach z interoperacyjnością, które pojawiają się podczas przekazywania wyników z jednego pakietu do kolejnego. Zamiast polegać na pojedynczym, generycznym kroku konwersji, przepływ pracy wykorzystuje dostosowane konwersje formatów na każdym interfejsie oprogramowania: przestrzenne dane zliczeń są kopiowane do standardowego slotu badania RNA przed integracją w Seurat; obiekty Seurat oraz obiekty referencyjne pojedynczych komórek są konwertowane na obiekty SingleCellExperiment dla SPOTlight; macierze zliczeń przestrzennych są reformatowane dla STdeconvolve; zliczenia, współrzędne i metadane typów komórek pochodzące z Seurat są konwertowane na obiekt Giotto w celu analizy komunikacji międzykomórkowej; a adnotacje ROI z Select Spatial Spots są eksportowane jako pliki CSV z CELL_ID, współrzędnymi X/Y i etykietami grup, aby można było je ponownie przypisać do obiektu Seurat. Kroki te pomagają użytkownikom zidentyfikować i skorygować typowe problemy, takie jak niekompatybilne sloty badań, niedopasowane identyfikatory spotów, brakujące kolumny metadanych, błędne formatowanie współrzędnych oraz niezgodności symboli genów liganda i receptora.

Dodatkową cechą tego przepływu pracy jest zastosowanie dwóch komplementarnych strategii dekonwolucji zamiast polegania na pojedynczym algorytmie. SPOTlight wykorzystuje wcześniejsze referencje scRNA-seq do referencyjnego szacowania proporcji typów komórek, podczas gdy STdeconvolve umożliwia bezreferencyjne wykrywanie ukrytych tematów transkrypcyjnych. W zbiorze danych demonstracyjnych tematy STdeconvolve były zgodne z regionami zdefiniowanymi histologicznie (Rycina 5B), co wspiera interpretowalność biologiczną. W niniejszym protokole nie przeprowadzono bezpośredniego porównania ilościowego obu metod; użytkownicy chcący przeprowadzić ich walidację na własnych danych mogą postępować zgodnie z ramami przedstawionymi w kroku 7.3.2 protokołu. Ponadto protokół wspiera użyteczność i identyfikowalność poprzez zastosowanie skryptów do funkcji podstawowych, oferując jednocześnie interfejs graficzny (Select Spatial Spots) do intuicyjnej izolacji ROI.

Powiązanie programów Seurat i Giotto pozwala na stworzenie przepływu pracy łączącego preprocessing oraz integrację opartą na Seurat z analizą sieci i statystyką przestrzenną w Giotto. Seurat zapewnia sprawdzone środowisko dla badaczy z doświadczeniem w scRNA-seq i wspiera integrację wielu próbek, natomiast przeniesienie zharmonizowanych danych do Giotto umożliwia konstrukcję sieci przestrzennych oraz analizę ligand-receptor. Taka konstrukcja pozwala użytkownikom wykorzystać udokumentowane zalety obu platform, bez sugerowania, że zostały one tutaj zestawione w ramach benchmarku z Squidpy lub innymi ramami analitycznymi. W przypadku komunikacji międzykomórkowej zastosowano analizę ligand-receptor w Giotto opartą na sieci Delaunaya, ponieważ uwzględnia ona sąsiedztwo przestrzenne w tym samym środowisku analitycznym; inne narzędzia, takie jak CellChat, oferują obszerne bazy danych o sygnalizacji, ale nie zostały one ocenione w ramach niniejszego protokołu¹⁷.

Należy wziąć pod uwagę kilka ograniczeń. Podejście do dekonwolucji oparte na referencji zależy krytycznie od dostępności wysokiej jakości, dopasowanych referencji jednokomórkowych. Ponadto podstawowe założenia analizy komunikacji międzykomórkowej opierają się na współekspresji transkrypcyjnej jako przybliżeniu fizycznej interakcji białek, co wymaga późniejszej walidacji eksperymentalnej. Wreszcie, niniejszy schemat analizy został opracowany przede wszystkim dla standardowych zbiorów danych Visium. Wraz z rozwojem dziedziny w kierunku technologii o rozdzielczości zbliżonej do pojedynczej komórki, takich jak Visium HD, kwestie analityczne ulegną zmianie; dane o wyższej rozdzielczości mogą wymagać innych parametrów wstępnego przetwarzania i zmniejszyć bezwzględną konieczność dekonwolucji spotów. Aby dostosować ten schemat do zbiorów danych Visium HD, binowe dane o wysokiej rozdzielczości można obliczeniowo zagregować w większe bin przestrzenne lub pominąć moduły dekonwolucji spotów na rzecz segmentacji komórek w oparciu o obraz¹⁸.

Cechy te sugerują potencjalną możliwość zastosowania w różnych dziedzinach biologii, w tym w biologii rozwoju, neuronauce, badaniach nad rakiem oraz immunologii^19,20,21,22. Modułowa konstrukcja pozwala badaczom na dostosowanie poszczególnych komponentów do ich potrzeb, niezależnie od tego, czy skupiają się oni na identyfikacji domen przestrzennych, komunikacji międzykomórkowej czy specjalizacji regionalnej. W miarę rozwoju technologii przestrzennych i rozszerzania się zbiorów danych, niniejszy protokół stanowi fundament, który może zostać rozbudowany o nowe metody analityczne w celu odpowiedzi na pojawiające się pytania biologiczne.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
ggplot2	Posit Software, PBC	v4.0.0(CRAN)	Zaawansowana wizualizacja danych
Giotto	Dries Lab	v4.2.2 (GitHub)	Analiza sieci przestrzennych i komunikacji międzykomórkowej
patchwork	Thomas Lin Pedersen	v1.3.2 (CRAN)	Kompozycja i układ wykresów
R software	R Foundation for Statistical Computing	v4.4.3	Podstawowe środowisko wykonawcze (macOS aarch64)
scater	Davis McCarthy et al.	v1.34.1 (Bioconductor)	Kontrola jakości i wizualizacja danych jednokomórkowych
scran	Aaron Lun et al.	v1.34.0 (Bioconductor)	Modelowanie wariancji i detekcja markerów w danych jednokomórkowych
Select Spatial Spots (Własne narzędzie Python)	LeafLight	v1.0.0 (GitHub)	Interaktywny wybór przestrzennych obszarów zainteresowania (ROI) (https://github.com/LeafLight/SelectSpatialSpots)
Seurat	Satija Lab	v5.3.0 (CRAN)	Preprocessing, integracja i klasteryzacja danych przestrzennych
SeuratObject	Satija Lab	v5.2.0 (CRAN)	Struktury danych dla danych jednokomórkowych i przestrzennych
SingleCellExperiment	Bioconductor Core Team	v1.28.1 (Bioconductor)	Standaryzowany kontener danych dla scRNA-seq
SPOTlight	Marc Elosua-Bayes et al.	v1.10.0 (Bioconductor)	Dekonwolucja przestrzenna z wykorzystaniem referencji
Stdeconvolve	Jean Fan Lab	v1.3.2 (Bioconductor)	Nienadzorowane modelowanie tematów latentnych
tidyverse	Posit Software, PBC	v2.0.0 (CRAN)	Podstawowy pakiet narzędzi do manipulacji i formatowania danych