Практический рабочий процесс анализа данных пространственной транскриптомики: от получения данных до углубленного анализа

Краткое содержание

В данном протоколе представлен воспроизводимый рабочий процесс анализа данных пространственной транскриптомики, который ведет пользователя от получения общедоступных данных и контроля качества на базе Seurat до интеграции, обнаружения пространственных признаков, деконволюции типов клеток, аннотирования областей интереса и анализа межклеточных взаимодействий с практическими контрольными точками, обеспечивающими прозрачность выполнения.

Аннотация

Пространственная транскриптомика (ST) позволяет профилировать экспрессию генов по всему геному, сохраняя при этом двухмерный пространственный контекст молекул мРНК в срезах тканей, что делает возможным изучение архитектуры тканей и биологии, связанной с микроокружением. Однако анализ ST остается сложной задачей, поскольку импорт данных, контроль качества, интеграция, деконволюция, пространственная статистика и визуализация часто требуют использования нескольких программных сред и воспроизводимого выбора параметров. В данном протоколе представлен практический вычислительный рабочий процесс для общедоступных наборов данных ST в среде R, который начинается с получения данных и настройки программного обеспечения и продолжается загрузкой данных с помощью Seurat, контролем качества, нормализацией, интеграцией нескольких образцов, кластеризацией и анализом пространственно вариабельных генов. Затем в рабочем процессе применяются взаимодополняющие стратегии деконволюции, включая анализ SPOTlight с использованием референса и ненаправляемое тематическое моделирование STdeconvolve, после чего следует анализ пространственного взаимодействия клеток с помощью Giotto и интерактивный выбор областей интереса (ROI) с использованием специализированного приложения Python Dash. Делая акцент на выполнении с помощью скриптов, четком обосновании параметров, ожидаемых результатах и контрольных точках для устранения неполадок, данный протокол предоставляет адаптируемую структуру для стандартных наборов данных ST на основе массивов и соответствующих платформ после оценки параметров, специфичных для конкретного набора данных и платформы.

Введение

Пространственная транскриптомика (ST) представляет собой группу инновационных технологий, позволяющих измерять экспрессию генов во всем геноме, сохраняя при этом пространственные координаты молекул матричной РНК (mRNA) в срезах тканей. Методы ST включают подходы на основе секвенирования с использованием массивов с позиционным баркодингом и методы визуализации in situ, которые картируют транскрипционные сигналы внутри неповрежденных микроокружений ткани^1,2. Благодаря сохранению пространственного контекста, ST позволяет анализировать архитектуру ткани, организацию клеточного соседства, межклеточную коммуникацию и связанные с микроокружением биологические процессы, которые невозможно полностью изучить после диссоциации ткани³.

Стремительный рост публичных репозиториев данных ST создал беспрецедентные возможности для вторичного анализа и разработки новых методов³. Такие ресурсы, как база данных CROST, аккумулируют сотни наборов данных пространственной транскриптомики для различных видов и технологических платформ, в то время как специализированные коллекции, такие как STOmicsDB, сосредоточены на конкретных методологиях, например Stereo-seq^4,5. Несмотря на изобилие данных, вычислительный анализ остается сложной задачей из-за сложности структур пространственных данных, разнообразия аналитических инструментов и технических препятствий при внедрении воспроизводимых рабочих процессов^{6,7,8,9,10,11}.

Для преодоления ограничений, связанных с использованием одной программной среды, в данной работе представлен интегрированный вычислительный рабочий процесс, использующий взаимодополняющие аналитические инструменты. Существующие комплексные экосистемы для анализа ST включают преимущественно Seurat, Giotto и фреймворки на базе Python, такие как Squidpy^6,7,12. Хотя инструменты на базе Python, такие как Squidpy, предоставляют широкие возможности для анализа пространственных графов, объединение основного вычислительного конвейера в рамках одного языка программирования минимизирует технические сложности при взаимодействии между разными языками. Следовательно, основной конвейер реализован преимущественно на языке R для снижения этих технических трудностей. В данном рабочем процессе на базе R пакет Seurat используется для загрузки данных, контроля качества, нормализации, снижения размерности, визуализации и интеграции нескольких образцов, что отражает его широкое применение в рабочих процессах секвенирования единичных клеток и пространственной транскриптомики. Затем Giotto применяется для построения пространственных сетей и анализа межклеточных взаимодействий на основе лиганд-рецепторных пар. Таким образом, данный конвейер объединяет предварительную обработку и интеграцию на базе Seurat с пространственным анализом на базе Giotto, обеспечивая при этом явный и воспроизводимый перенос данных между двумя наборами инструментов.

В рамках данного подхода реализованы две взаимодополняющие стратегии деконволюции: SPOTlight — метод с использованием референса, который применяет данные scRNA-seq для оценки пропорций типов клеток, и STdeconvolve — метод обучения без учителя на основе тематического моделирования, который выявляет латентные транскрипционные паттерны^8,11. Результаты этих методов предоставляют взаимодополняющие представления о пространственной клеточной гетерогенности, однако они не рассматриваются как количественная перекрестная проверка, если пользователь не выполняет дополнительный анализ конкордантности, описанный в протоколе. Для интерактивной аннотации областей интереса (ROI) интегрировано специализированное приложение на Python Dash под названием Select Spatial Spots, которое экспортирует стандартные файлы аннотаций на основе координат, пригодные для последующего анализа.

Что касается практической применимости, данный рабочий процесс предназначен прежде всего для стандартных данных ST на основе массивов (например, Visium с разрешением 55 µm) и может быть адаптирован для других типов тканей после оценки параметров. Перед началом анализа следует учитывать основные ограничения. Во-первых, модуль деконволюции с использованием референса зависит от наличия высококачественного scRNA-seq референса, соответствующего типу ткани. Во-вторых, платформы с субклеточным разрешением или разрешением, близким к одноклеточному, могут потребовать модифицированной предварительной обработки, агрегации пространственных бинов или сегментации клеток на основе изображений перед интеграцией². Репрезентативный набор данных толстой кишки мыши используется в качестве демонстрационного примера, чтобы показать, как данный рабочий процесс может оценивать пространственные домены и организацию тканей, определяемую маркерами, а не как доказательство универсальной совместимости платформы.

Протокол

Все биологические наборы данных, проанализированные в данном протоколе, находятся в открытом доступе и используются исключительно в демонстрационных целях. Ссылки на конкретные записи данных и репозитории источников приведены в соответствующих этапах. Исходные наборы данных были созданы исследователями в соответствии с институциональными этическими нормами, применимыми к каждому исходному исследованию. Обратитесь к Таблице материалов, чтобы проверить все необходимые версии программного обеспечения и пакетов R.

Требования к оборудованию: объем оперативной памяти, необходимый для данного рабочего процесса, масштабируется в зависимости от количества анализируемых образцов и спотов. Для типичного набора данных пространственной транскриптомики (например, примерно 3 000 спотов на образец для трех образцов) для выполнения конвейера достаточно стандартной рабочей станции с минимум 16 GB оперативной памяти. Однако для обеспечения оптимальной производительности и стабильности настоятельно рекомендуется использовать 32 GB оперативной памяти или более, особенно на этапах ресурсоемких вычислений, таких как нормализация SCTransform и факторизация матриц при деконволюции.

1. Сбор данных и подготовка структуры каталогов

1. Получите общедоступные наборы данных пространственной транскриптомики (например, GEO accession GSE169749, дата подачи 26 марта 2021 г., дата последнего обновления 6 марта 2022 г.) и референсные данные секвенирования РНК единичных клеток (например, GEO accession GSE264408, дата подачи 19 апреля 2024 г., дата последнего обновления 10 декабря 2024 г.).
2. Загрузите архив с необработанными данными (например, GSE169749_RAW.tar) из репозитория и извлеките содержимое архива.
   ПРИМЕЧАНИЕ: Для репрезентативного набора данных, используемого в данном протоколе, извлеченный каталог содержит несколько файлов для каждого образца, причем имена файлов обычно включают номер доступа образца в GEO (например, GSM5213483).
3. Организуйте файлы в стандартной структуре каталогов, совместимой с функцией Seurat Load10X_Spatial.
   1. Создайте основной каталог данных (например, ./data/).
   2. Создайте отдельный подкаталог для каждого образца (например, ./data/sample_A1/).
   3. Перенесите (скопируйте или переместите) следующие основные файлы для каждого образца в соответствующий подкаталог.
      1. Поместите файл filtered_feature_bc_matrix.h5 в основной подкаталог образца.
      2. Создайте подпапку spatial/ внутри каталога каждого образца.
      3. Поместите следующие файлы в подпапку spatial/: tissue_positions_list.csv.gz, scalefactors_json.json.gz и tissue_hires_image.png.gz.
   4. Распакуйте файлы .gz в папке spatial/. Убедитесь, что исходные имена файлов остались точно такими, как того требует функция Load10X_Spatial (например, tissue_positions_list.csv).
      ПРИМЕЧАНИЕ: Итоговая организованная структура для корректной загрузки должна соответствовать следующей схеме (на примере образца A1):
      ./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. Настройка программной среды

1. Установите R (версии 4.4.3 или выше).
2. Загрузите скрипты для анализа (от 1_ReadSpatialData.R до 8_Giotto_Communication.R) из репозитория GitHub (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow, хэш коммита 2d85e18 для обеспечения воспроизводимости).
3. Установите необходимые пакеты R из CRAN и Bioconductor, запустив скрипт setup.R.
   ПРИМЕЧАНИЕ: Полный список конкретных версий пакетов, использованных в данном рабочем процессе, приведен в Таблице материалов, а также доступен в файле session_info.txt в соответствующем репозитории GitHub (https://github.com/LeafLight/SpatialTranscriptomicsWorkFlow).
4. Установите Giotto Suite, выполнив команды установки, указанные в официальной документации.
   1. Перейдите по официальному URL-адресу установки для получения скриптов настройки: https://giottosuite.com/articles/installation.html
   2. Инициализируйте необходимую среду Python и системные зависимости согласно пошаговым инструкциям.
5. Установите дополнительные пакеты для улучшенной визуализации.
   1. Выполните: install.packages(c("ggprism", "colorBlindness"))
6. Установите специализированный инструмент Select Spatial Spots. Убедитесь, что используется операционная система Windows, macOS или Linux и установлен Python версии 3.8 или выше. Этот шаг является необязательным.
   1. Получите инструмент, перейдя в репозиторий GitHub (https://github.com/LeafLight/SelectSpatialSpots, хэш коммита d20946e для обеспечения воспроизводимости), и загрузите исходный код.
   2. Перейдите в директорию инструмента и установите зависимости Python: pip install -r requirements.txt
7. Установите пакет glmGamPoi , выполнив команду BiocManager::install("glmGamPoi") для ускорения нормализации SCTransform.
   ПРИМЕЧАНИЕ: Основной рабочий процесс требует только пакетов, перечисленных в пунктах 2.1–2.3 и 2.6. Пакет Rfast2 используется для ускорения вычисления статистики I Морана. Необязательные пакеты в пункте 2.5 предназначены для создания графиков для публикаций с темой в стиле Prism (ggprism) и для доступа к палитре, подходящей для людей с нарушением цветовосприятия (colorBlindness). Инструмент в пункте 2.6 обеспечивает функциональность интерактивного выбора спотов. Выполните команду sessionInfo() и сохраните полный вывод консоли в текстовый файл. Это позволит задокументировать точные версии программного обеспечения и зависимости пакетов, что критически важно для обеспечения долгосрочной воспроизводимости.

3. Загрузка пространственных данных и контроль качества (1_ReadSpatialData.R, 2_SpatialDataQC.R)

1. Загрузите пространственные данные в объекты Seurat.
   1. Применение Изображение_Read10X для ручной загрузки изображения ткани с высоким разрешением с указанием директория изображения и имя изображения.
   2. Использование Загрузка10X_Spatial с помощью изображение параметр, присваиваемый объекту изображения, созданному на этапе 3.1.1, для создания объекта Seurat.
2. Рассчитайте показатели контроля качества.
   1. Вычислите процент митохондриальных прочтений с помощью Процентный набор признаков с паттерном ^mt-.
3. Визуализируйте и интерпретируйте данные на основе показателей контроля качества.
   1. Постройте скрипичные диаграммы для nCount_Spatial, nFeature_Spatialи процент мт-ДНК использование VlnPlot.
   2. Создайте графики пространственных признаков для этих метрик с помощью SpatialFeaturePlot для идентификации пятен за пределами области ткани.
      ПРИМЕЧАНИЕ: Визуальный контроль: участки за пределами области ткани обычно характеризуются низким количеством UMI (nCount_Spatial < 500) и низкий уровень детекции гена (nFeature_Spatial < 200)
   3. Для демонстрационных целей примените фильтры для удаления пятен низкого качества (например, subset(seurat_obj, subset = nFeature_Spatial) > 200 & nCount_Spatial > 500)). Этот шаг является необязательным.
      ПРИМЕЧАНИЕ: Основной целью контроля качества (КК) в пространственной транскриптомике является выявление и аннотирование технических артефактов, таких как споты вне пределов ткани. Для демонстрационных данных Visium по толстой кишке параметр nFeature_Spatial > 200 и nCount_Spatial > 500 удалите пятна с низкой сложностью/фоновым уровнем. Для тканей с низким содержанием РНК или деградировавших срезов перед повышением пороговых значений изучите скрипичные диаграммы (violin plots) и графики пространственных признаков; для плотных тканей или образцов с высоким количеством UMI могут потребоваться более строгие пороги. Избегайте фильтрации исключительно по высоким значениям экспрессии, если только не наблюдаются дуплеты, складки ткани или явные артефакты. Фильтрация часто не рекомендуется, так как она физически удаляет пространственные локации, что может нарушить непрерывность структуры ткани для последующего пространственного анализа.

4. Предобработка, интеграция и кластеризация данных (3_IntegrationAndClustering.R)

1. Нормализуйте и проведите предварительную обработку отдельных образцов.
   1. Примените нормализацию SCTransform к каждому образцу отдельно с параметром assay = Spatial.
2. Интегрируйте несколько образцов.
   1. Подготовьте список SCT-нормализованных объектов для интеграции. Убедитесь, что каждый объект имеет анализ "RNA", скопировав анализ "Spatial": spatial_list[[1]][["RNA"]] <- spatial_list[[1]][["Spatial"]]
      ​ПРИМЕЧАНИЕ: Копирование анализа "Spatial" в стандартный слот анализа "RNA" является необходимым обходным решением для обеспечения совместимости с последующими функциями интеграции, изначально разработанными для одноклеточных объектов Seurat.
   2. Используйте функции SelectIntegrationFeatures и PrepSCTIntegration для определения общих вариабельных признаков.
   3. Найдите якоря интеграции с помощью FindIntegrationAnchors при normalization.method = "SCT".
   4. Интегрируйте данные с помощью IntegrateData при normalization.method = "SCT".
3. Выполните снижение размерности и кластеризацию интегрированного анализа.
   1. Запустите PCA для интегрированных данных с помощью RunPCA.
   2. Определите оптимальное количество главных компонентов (PCs) для последующего анализа путем вычисления кумулятивной объясненной дисперсии. Определите точку перегиба («локоть») программно (например, точку, где кумулятивная дисперсия превышает 90%, а предельный прирост падает ниже 0,1%), используя код в файле 3_IntegrationAndClustering.R, строки 36–38. Полученное количество PCs далее будет обозначаться как pc.use.
   3. Запустите RunUMAP с параметром dims = 1:pc.use.
   4. Проведите кластеризацию клеток с помощью FindNeighbors при dims = 1:pc.use и FindClusters при resolution = 0.5. Корректируйте разрешение только после проверки стабильности кластеров и согласованности маркерных генов.
   5. Проведите анализ дифференциальной экспрессии между целевыми группами (например, "B1_colon_d14" против "A1_colon_d0") с помощью функции FindMarkers.
      ​ПРИМЕЧАНИЕ: Визуальный контроль: при успешной интеграции на графике UMAP будет наблюдаться надлежащее смешивание образцов при сохранении биологически различных кластеров.
4. Идентифицируйте пространственно вариабельные гены.
   1. Для каждого исходного образца запустите FindSpatiallyVariableFeatures, используя метод "moransi" в анализе "SCT" для вычисления пространственной автокорреляции.

5. Предварительная обработка референсных данных одноклеточного секвенирования (4_scDataPreProcessing.R)

1. Считайте матрицу количеств RNA-seq единичных клеток с помощью Read10X и создайте объект Seurat.
2. Выполните стандартный контроль качества (КК), нормализацию и кластеризацию.
   1. Рассчитайте процент митохондриальных ридов (percent.mt) и отфильтруйте клетки (например, nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 25).
   2. Нормализуйте данные с помощью SCTransform, установив vars.to.regress = "percent.mt".
   3. Выполните PCA и UMAP, а затем проведите кластеризацию клеток, используя метод динамического выбора главных компонент (PC), описанный в шаге 4.3.2.
3. Аннотируйте типы клеток.
   1. Рассчитайте показатели модулей (module scores) для канонических ген-маркеров типов клеток с помощью AddModuleScore.
   2. Аннотируйте кластеры на основе показателей модулей и известных биологических данных. В качестве альтернативы импортируйте предварительно рассчитанные аннотации из метаданных.

6. Деконволюция с использованием референсов с помощью SPOTlight (5_SPOTlight_Deconv.R)

1. Подготовка данных для SPOTlight.
   1. Преобразовать аннотированный объект Seurat с данными одноклеточного секвенирования и пространственный объект Seurat в объекты SingleCellExperiment.
   2. Выполнить лог-нормализацию одноклеточных данных с помощью logNormCounts.
2. Запуск деконволюции SPOTlight.
   1. Определить гипервариабельные гены (HVGs) в одноклеточных данных с помощью modelGeneVar и getTopHVGs.
   2. Рассчитать маркерные гены типов клеток с помощью scoreMarkers и отфильтровать высококачественные маркеры (например, mean.AUC > 0.8).
   3. Выполнить даунсэмплинг одноклеточной референсной выборки для каждого типа клеток до приемлемого количества (например, 50 клеток) для сокращения времени вычислений.
   4. Выполнить деконволюцию с помощью функции SPOTlight с параметрами weight_id = "mean.AUC", group_id = "cluster" и gene_id = "gene".
3. Визуализация и экспорт результатов.
   1. Извлечь матрицу результатов деконволюции (пропорции типов клеток на один спот).
   2. Использовать plotSpatialScatterpie для визуализации клеточного состава в пространственных координатах.
   3. Добавить результаты деконволюции в метаданные пространственного объекта Seurat с помощью AddMetaData. Ожидаемый результат: матрица деконволюции со столбцами пропорций типов клеток для каждого пространственного спота, диаграммы пространственного распределения (scatterpie plots), отображающие локальный клеточный состав, и объект Seurat, содержащий пропорции деконволюции в виде метаданных.

7. Деконволюция без использования референсов с помощью STdeconvolve (7_STdeconvolve.R)

1. Подготовка пространственных данных.
   1. Извлеките матрицу необработанных количественных данных (raw count matrix) из пространственного объекта Seurat с помощью GetAssayData с параметром slot = "counts".
   2. Удалите низкокачественные споты и гены, используя функцию cleanCounts из пакета STdeconvolve (например, min.lib.size = 100).
2. Идентификация латентных типов клеток.
   1. Дополнительно отфильтруйте корпус генов, оставив те, которые экспрессируются в минимальной фракции спотов, используя restrictCorpus (например, removeAbove=1.0, removeBelow = 0.05).
   2. Постройте модель латентного размещения Дирихле (LDA) для диапазона потенциальных количеств тем (K) (например, Ks = seq(2, 9, by = 1)), используя функцию fitLDA.
   3. Выберите оптимальную модель на основе минимальной перплексии с помощью optimalModel с параметром opt = "min".
3. Анализ и визуализация результатов.
   1. Извлеките пропорции типов клеток (тем) (theta) и профили экспрессии генов (beta) из оптимальной модели с помощью getBetaTheta.
   2. Для облегчения биологической интерпретации деконволивированных тем импортируйте аннотации областей интереса (ROI), созданные с помощью инструмента "Select Spatial Spots" (см. шаг 9.3.3). Используйте эти аннотации в качестве параметра groups в функции vizAllTopics, чтобы спроецировать деконволивированные пропорции типов клеток обратно на пространственные координаты и раскрасить споты в соответствии с их ROI. Этот шаг является необязательным.
      ПРИМЕЧАНИЕ: Ожидаемый результат: STdeconvolve должен вернуть значения theta, представляющие пропорции тем на один спот, и значения beta, представляющие профили генов для каждой темы. Если пользователи хотят количественно сравнить результаты ненаправляемой деконволюции STdeconvolve с результатами SPOTlight, основанными на референсе (из шага 6), они могут экспортировать обе матрицы пропорций (theta из STdeconvolve и матрицу SPOTlight) и вычислить корреляцию или показатели согласованности для каждого спота (например, корреляцию Пирсона или Спирмена), используя стандартные функции R. Количественное сравнение SPOTlight и STdeconvolve является необязательным и не требуется для выполнения основного рабочего процесса; однако соответствующие выходные матрицы из шагов 6.3.1 и 7.3.1 доступны пользователям, желающим провести такой анализ.

8. Пространственное межклеточное взаимодействие с использованием Giotto(8_Giotto_Communication.R)

1. Преобразуйте объект Seurat в объект Giotto.
   1. Используйте функцию createGiottoObject, указав матрицу необработанных значений экспрессии (raw count matrix) и пространственные координаты.
2. Выполните предварительную обработку объекта Giotto и добавьте результаты деконволюции.
   1. Нормализуйте данные с помощью normalizeGiotto.
   2. Добавьте аннотации типов клеток (например, основной тип клеток, полученный в результате деконволюции SPOTlight) в метаданные клеток с помощью addCellMetadata.
3. Проведите анализ пространственно-обусловленного межклеточного взаимодействия.
   1. Создайте пространственную сеть с помощью createSpatialNetwork с параметрами method = "Delaunay" и name = "Delaunay_network".
   2. Загрузите базу данных лигандов и рецепторов в среду R. Убедитесь, что загруженный объект представлен в виде таблицы (data frame) с двумя отдельными столбцами, содержащими символы генов лиганда и рецептора (например, скомпилированную сеть для мыши из Zenodo: https://zenodo.org/api/records/15168114/files/lr_network_mouse.csv/content, версия v6, опубликована 7 апреля 2025 г.).
   3. Запустите exprCellCellcom с параметром cluster_column = "celltype_major" для выявления значимых взаимодействий лиганд-рецептор между типами клеток, находящимися в пространственной близости. Ожидаемый результат: таблица, содержащая пары лиганд-рецептор, комбинации типов клеток-источников и клеток-мишеней, значения log2 fold-change и скорректированные значения P; для визуализации оставьте только значимые взаимодействия (например, p.adj < 0.05).

9. Интерактивный выбор спотов с помощью функции select spatial spots (6_SelectSpatialSpots.R)

1. Подготовьте данные для интерактивного инструмента.
   1. Извлеките пространственные координаты из объекта Seurat с помощью функции GetTissueCoordinates.
2. Отформатируйте и экспортируйте данные.
   1. Приведите фрейм данных с координатами к виду, содержащему следующие точные заголовки столбцов: CELL_ID, X и Y.
   2. Экспортируйте отформатированный фрейм данных в CSV-файл.
3. Проведите анализ областей интереса (ROI).
   1. Запустите специализированное приложение Dash Select Spatial Spots и загрузите CSV-файл.
   2. Интерактивно выберите споты на основе их пространственного расположения.
   3. Экспортируйте список выбранных спотов и присвоенные им метки групп/ROI в новый CSV-файл. Ожидаемый результат: CSV-файл с идентификаторами спотов, координатами X/Y и присвоенными метками групп/ROI, которые могут быть сопоставлены с объектом Seurat по CELL_ID.

Результаты

Реализация рабочего процесса и интеграция данных иллюстрируют основные особенности ткани

Вычислительный рабочий процесс был применен к данным пространственной транскриптомики толстой кишки мыши для иллюстрации ожидаемых результатов на различных этапах анализа. Как показано на схеме рабочего процесса (Рисунок 1), конвейер начинался со сбора данных и контроля качества, где графики пространственных признаков очерчивали границы ткани (Рисунок 2A,B). Затем для уменьшения технических эффектов батча при сохранении интерпретируемой биологической вариабельности использовался рабочий процесс интеграции на основе «якорей» (anchors) Seurat. Визуализации UMAP продемонстрировали выравнивание образцов и паттерны пространственной кластеризации после интеграции (Рисунок 2C,D). Для управления снижением размерности и последующей кластеризацией был внедрен количественный динамический выбор главных компонентов (PCs) на основе кумулятивной дисперсии (см. Дополнительный рисунок 1). Анализ тепловых карт маркерных генов выявил отчетливые транскрипционные профили, лежащие в основе пространственных кластеров (Рисунок 2E).

Чтобы оценить, соответствуют ли вычислительные кластеры известной анатомической архитектуре гистологии толстой кишки, были проанализированы профили экспрессии канонических маркерных генов специфических слоев. В слое слизистого эпителия наблюдалась экспрессия маркеров эпителиальных клеток, включая Epcam и Krt8, наряду с маркером бокаловидных клеток Muc2. Мезенхимальные и стромальные маркеры, такие как Col1a1 и Vim, маркировали области собственной пластинки и подслизистого слоя, тогда как внешний слой мышечной оболочки определялся по структурным генам гладких мышц, таким как Acta2 и Tagln. Пространственная ограниченность этих ассоциированных с определенными линиями маркеров подтверждает интерпретацию, согласно которой рабочий процесс интеграции и кластеризации сохранил основные гистологические слои ткани толстой кишки вдоль оси от слизистой оболочки к мышечной (см. Дополнительный рисунок 2).

После валидации кластеров был проведен последующий анализ дифференциальной экспрессии для выявления дифференциально экспрессируемых генов (DEGs) между экспериментальными условиями (Рисунок 2F,G). Кроме того, с помощью статистики I Морана были определены пространственно вариабельные гены, что позволило выделить гены со значимым неслучайным пространственным распределением в ткани (Рисунок 2H).

Клеточная деконволюция и сети пространственных взаимодействий раскрывают микроорганизацию ткани

Обработка референсных данных scRNA-seq позволила получить аннотации, подтвержденные фильтрацией по контролю качества (Рисунок 3A), ненаправленной кластеризацией (Рисунок 3B), валидацией по marker-генам (Рисунок 3C) и соответствием независимым аннотациям (Рисунок 3D). Клеточный состав (Рисунок 3E) лег в основу стратегии даунсемплинга для деконволюции. С помощью SPOTlight были оценены пропорции типов клеток с использованием референса в пространственных спотах (Рисунок 4A,B), в то время как STdeconvolve обеспечил ненаправленный анализ пространственных клеточных паттернов с помощью тематического моделирования (Рисунок 5B). Специализированный инструмент Select Spatial Spots позволил сопоставить эти паттерны с гистологическим контекстом (Рисунок 5A). Наконец, на основе результатов деконволюции типов клеток был проведен анализ пространственной коммуникации, в ходе которого выявлены лиганд-рецепторные взаимодействия между пространственно близко расположенными группами типов клеток (Рисунок 6A,B).

Анализ проблем, выявленных в ходе оптимизации протокола

В ходе оптимизации протокола был выявлен ряд проблем, которые легли в основу практических контрольных точек. Неоптимальные результаты деконволюции наблюдались в случаях плохого соответствия одноклеточных референсов контексту ткани, что указывает на необходимость использования данных scRNA-seq, соответствующих типу ткани и виду организма, при их наличии. Первоначальные попытки кластеризации с параметрами по умолчанию не всегда позволяли выделить ожидаемые биологические структуры; анализ выбора главных компонентов (PC), разрешения кластеризации и когерентности маркерных генов помог определить пространственно интерпретируемые домены, соответствующие анатомии ткани. Эти наблюдения служат практическими примерами того, как пользователи могут диагностировать типичные аналитические проблемы при выполнении данного рабочего процесса.

Рисунок 1Рабочий процесс интегрированного анализа пространственной транскриптомики. Схематическое представление аналитического конвейера: от сбора и предварительной обработки данных до углубленного пространственного анализа. Основные этапы включают: (1) загрузку данных, контроль качества и интеграцию нескольких образцов с помощью Seurat; (2) пространственную кластеризацию и выявление пространственно-вариабельных генов; (3) деконволюцию типов клеток. посредством методы на основе референса (SPOTlight) и методы без учителя (STdeconvolve); (4) анализ пространственных межклеточных взаимодействий с помощью Giotto и интерактивный выбор областей интереса с использованием специализированного инструмента; Выбор пространственных точекРезультаты всех модулей синтезируются для получения биологических данных об архитектуре тканей и клеточном микроокружении. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этого рисунка.

Рисунок 2: Интеграция данных, кластеризация и анализ дифференциальной экспрессии. (A,B) Метрики контроля качества для пространственных образцов A1 и B1, демонстрирующие распределение количества генов, количества UMI и процента митохондриальных генов. (C) Визуализация UMAP интегрированных данных пространственной транскриптомики, расцвеченная по происхождению образца (слева) и идентификатору кластера (справа). (D) Пространственная проекция идентификаторов кластеров на срезы тканей. (E) Тепловая карта основных маркерных генов для каждого пространственного кластера. (F) График-вулкан (volcano plot), отображающий дифференциально экспрессируемые гены между условиями A1_colon_d0 и B1_colon_d14. (G) Пространственные паттерны экспрессии репрезентативных дифференциально экспрессируемых генов на срезах тканей. (H) Карты пространственной экспрессии основных пространственно-вариабельных генов, идентифицированных via I-статистику Морана, где две левые панели отображают гены из образца A1_colon_d0, а две правые панели — гены из образца B1_colon_d14. Пожалуйста, нажмите здесь, чтобы просмотреть этот рисунок в большем размере.

Рисунок 3: Обработка и аннотация эталонных данных одноклеточного секвенирования. (A) Показатели контроля качества эталонных данных scRNA-seq до и после фильтрации. (B) Визуализация данных scRNA-seq с помощью UMAP, где цвета соответствуют кластерам, определенным методом обучения без учителя. (C) Точечный график (dot plot), отображающий уровни экспрессии канонических маркерных генов типов клеток в различных кластерах. (D) Аннотированная визуализация данных scRNA-seq с помощью UMAP с указанием основных типов клеток. (E) Клеточный состав эталонного набора данных scRNA-seq. Красная пунктирная линия указывает порог даунсэмплинга (n = 50 клеток на тип), примененный при деконволюции SPOTlight для баланса между вычислительной эффективностью и репрезентативностью типов клеток. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этого рисунка.

Рисунок 4: Пространственная деконволюция клеточной гетерогенности. (A,B) Пространственные круговые диаграммы (scatterpie plots), полученные с помощью деконволюции SPOTlight, демонстрирующие пропорциональный состав основных типов клеток в каждом споте для образцов A1 (A) и B1 (B). (C) Репрезентативное пространственное распределение B-клеток в образцах A1 (слева) и B1 (справа), демонстрирующее паттерны локализации конкретной популяции иммунных клеток с пространственным разрешением, выявленные в результате деконволюции. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этого рисунка.

Рисунок 5: Интерактивный анализ выбранных областей интереса и сравнение результатов ненаправленной деконволюции. (A) Интерфейс специализированного инструмента «Select Spatial Spots», демонстрирующий интерактивный выбор областей, соответствующих проксимальному отделу ободочной кишки, дистальному отделу ободочной кишки и другим тканевым доменам. (B) Визуализация результатов ненаправленной деконволюции (STdeconvolve) для образца A1 в виде пространственной диаграммы scatterpie, где споты окрашены в соответствии с областями, размеченными вручную на панели (A), что иллюстрирует соответствие между аннотацией на основе гистологии и расчетными распределениями клеточных тем. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этого рисунка.

Рисунок 6: Сети межклеточного взаимодействия с учетом пространственной организации. (A,B) Сети взаимодействия лиганд-рецептор, реконструированные с помощью Giotto для образцов A1 (A) и B1 (B). Узлы представляют типы клеток, ребра — значимые пары лиганд-рецептор (FDR < 0.05), а толщина ребер соответствует силе взаимодействия. Для обеспечения сопоставимости и четкости визуализации для всех образцов был применен единый порог значимости (FDR < 0.05), и для каждого состояния отображены 20 наиболее сильных взаимодействий, ранжированных по log2FC. Сети демонстрируют специфические для типов клеток паттерны коммуникации в пространственном контексте ткани толстой кишки. Пожалуйста, нажмите здесь, чтобы просмотреть увеличенную версию этого рисунка.

Дополнительный рисунок 1: Количественная оценка оптимизации параметров для снижения размерности.График «локтя» демонстрирует программный подход рабочего процесса к динамическому выбору оптимального количества главных компонент (PCs). Выбор рассчитывается на основе кумулятивного стандартного отклонения и пороговых значений предельной дисперсии, обозначенных красной вертикальной линией, чтобы зафиксировать биологическую вариативность при одновременном снижении технического шума перед последующей кластеризацией.Пожалуйста, нажмите здесь, чтобы скачать этот файл.

Дополнительный рисунок 2: Валидация пространственной кластеризации с использованием канонических маркеров специфических слоев толстой кишки. (A) Точечная диаграмма, демонстрирующая повышенную экспрессию маркеров эпителия, стромы и гладких мышц в вычислительных кластерах. (B) Графики пространственных признаков, отображающие репрезентативные маркеры (Epcam, Col1a1, Acta2) на координатах ткани.Пожалуйста, нажмите здесь, чтобы скачать этот файл.

Обсуждение

В данном протоколе представлен комплексный вычислительный рабочий процесс для анализа данных пространственной транскриптомики, который сочетает в себе глубину анализа с практической доступностью. Пошаговое руководство проводит исследователя через весь аналитический конвейер — от первоначального получения данных до углубленного пространственного анализа, уделяя особое внимание критическим точкам принятия решений и потенциальным ошибкам.

Несколько этапов протокола требуют особого внимания ввиду их влияния на последующие результаты. Пороги контроля качества и фильтрации должны быть тщательно адаптированы к конкретным типам тканей и технологическим платформам, поскольку излишне строгая фильтрация может привести к удалению биологически значимых спотов, в то время как слишком мягкие пороги могут привести к появлению технического шума. Выбор метода нормализации существенно влияет на результаты последующей кластеризации и анализа дифференциальной экспрессии. Например, в данном конвейере используется SCTransform, а не стандартная логарифмическая нормализация, поскольку было показано, что SCTransform позволяет моделировать и снижать техническую вариабельность, связанную с глубиной секвенирования и другими техническими эффектами в одноклеточных и пространственных рабочих процессах^14,15. При интеграции выбор признаков интеграции и параметров разрешения требует тщательного рассмотрения для достижения баланса между силой интеграции и сохранением биологического сигнала. Для определения пространственно-вариабельных генов был выбран индекс Морана (Moran's I) из-за его вычислительной масштабируемости и пригодности для анализа пространственной автокорреляции в данных на основе микрочипов, хотя в зависимости от размера набора данных и целей исследования могут быть рассмотрены альтернативные инструменты, такие как SPARK¹⁶.

Раздел по поиску и устранению неисправностей в данном рабочем процессе сосредоточен на типичных проблемах совместимости, которые возникают при передаче выходных данных из одного пакета в другой. Вместо использования одного общего этапа конвертации, рабочий процесс предусматривает специализированное преобразование форматов на каждом программном интерфейсе: данные пространственного подсчета копируются в стандартный слот анализа РНК перед интеграцией в Seurat; объекты Seurat и референсные объекты одноклеточного анализа преобразуются в объекты SingleCellExperiment для SPOTlight; матрицы пространственных подсчетов переформатируются для STdeconvolve; полученные в Seurat данные о подсчетах, координаты и метаданные типов клеток конвертируются в объект Giotto для анализа межклеточных взаимодействий; а аннотации областей интереса (ROI) из Select Spatial Spots экспортируются в виде CSV-файлов с указанием CELL_ID, координат X/Y и меток групп для их последующего сопоставления с объектом Seurat. Эти этапы помогают пользователям выявлять и исправлять распространенные ошибки, такие как несовместимые слоты анализа, несоответствие идентификаторов спотов, отсутствие столбцов метаданных, неправильный формат координат и несоответствие символов генов лигандов и рецепторов.

Дополнительной особенностью данного рабочего процесса является использование двух взаимодополняющих стратегий деконволюции вместо опоры на один алгоритм. SPOTlight использует имеющиеся референсные данные scRNA-seq для оценки пропорций типов клеток с использованием эталона, в то время как STdeconvolve позволяет выявлять латентные транскрипционные темы без использования референса. В демонстрационном наборе данных темы STdeconvolve совпали с гистологически определенными областями (Рисунок 5B), что подтверждает биологическую интерпретируемость. Прямое количественное сравнение этих двух методов в данном протоколе не проводилось; пользователи, желающие провести их сравнительный анализ на собственных данных, могут воспользоваться алгоритмом, описанным в шаге протокола 7.3.2. Кроме того, протокол обеспечивает удобство использования и прослеживаемость за счет применения скриптов для основных функций при наличии графического интерфейса (Select Spatial Spots) для интуитивно понятной изоляции ROI.

Связывание Seurat и Giotto позволяет объединить в рабочем процессе предварительную обработку и интеграцию на базе Seurat с пространственной статистикой и сетевым анализом на базе Giotto. Seurat предоставляет привычную среду для исследователей с опытом работы с scRNA-seq и поддерживает интеграцию нескольких образцов, в то время как перенос гармонизированных данных в Giotto позволяет осуществлять построение пространственных сетей и анализ лиганд-рецепторных взаимодействий. Такая архитектура позволяет пользователям использовать задокументированные преимущества обеих платформ, при этом данные инструменты не сопоставлялись в данном исследовании со Squidpy или другими фреймворками. Для анализа межклеточных коммуникаций использовался метод вывода лиганд-рецепторных взаимодействий Giotto на основе сети Делоне, поскольку он учитывает пространственную смежность в той же аналитической среде; такие фреймворки, как CellChat, предоставляют обширные базы данных сигнальных путей, но не оценивались в данном протоколе¹⁷.

Необходимо учитывать ряд ограничений. Подход к деконволюции на основе референсов критически зависит от наличия высококачественных соответствующих референсов одноклеточного секвенирования. Кроме того, фундаментальные допущения, лежащие в основе анализа межклеточного взаимодействия, опираются на коэкспрессию транскриптов как косвенный признак физического белкового взаимодействия, что требует последующей экспериментальной валидации. Наконец, данный рабочий процесс был разработан преимущественно для стандартных наборов данных Visium. По мере развития области в сторону технологий с разрешением, близким к одноклеточному, таких как Visium HD, аналитические подходы изменятся; данные более высокого разрешения могут потребовать иных параметров предварительной обработки и снизить абсолютную необходимость в деконволюции спотов. Чтобы адаптировать этот рабочий процесс для наборов данных Visium HD, бины высокого разрешения можно вычислительно объединить в более крупные пространственные бины или пропустить модули деконволюции спотов в пользу сегментации клеток на основе изображений¹⁸.

Эти особенности указывают на потенциальную применимость метода в различных областях биологии, включая биологию развития, нейронауки, исследования рака и иммунологию^19,20,21,22. Модульная архитектура позволяет исследователям адаптировать отдельные компоненты под свои задачи, будь то идентификация пространственных доменов, изучение межклеточных взаимодействий или региональной специализации. По мере развития пространственных технологий и расширения наборов данных данный протокол служит основой, которую можно дополнять новыми аналитическими методами для решения возникающих биологических вопросов.

Материалы

Список материалов, использованных в этой статье

Имя	Компания	Каталожный номер	Комментарии
ggplot2	Posit Software, PBC	v4.0.0(CRAN)	Продвинутая визуализация данных
Giotto	Dries Lab	v4.2.2 (GitHub)	Анализ пространственных сетей и межклеточных взаимодействий
patchwork	Thomas Lin Pedersen	v1.3.2 (CRAN)	Композиция и компоновка графиков
Программное обеспечение R	R Foundation for Statistical Computing	v4.4.3	Основная среда выполнения (macOS aarch64)
scater	Davis McCarthy et al.	v1.34.1 (Bioconductor)	Контроль качества и визуализация данных единичных клеток
scran	Aaron Lun et al.	v1.34.0 (Bioconductor)	Моделирование дисперсии и поиск маркеров в данных единичных клеток
Select Spatial Spots (пользовательский инструмент Python)	LeafLight	v1.0.0 (GitHub)	Интерактивный выбор пространственных областей интереса (ROI) (https://github.com/LeafLight/SelectSpatialSpots)
Seurat	Satija Lab	v5.3.0 (CRAN)	Препроцессинг, интеграция и кластеризация пространственных данных
SeuratObject	Satija Lab	v5.2.0 (CRAN)	Структуры данных для анализа единичных клеток и пространственных данных
SingleCellExperiment	Bioconductor Core Team	v1.28.1 (Bioconductor)	Стандартизированный контейнер данных для scRNA-seq
SPOTlight	Marc Elosua-Bayes et al.	v1.10.0 (Bioconductor)	Пространственная деконволюция с использованием референса
Stdeconvolve	Jean Fan Lab	v1.3.2 (Bioconductor)	Неконтролируемое моделирование латентных тем
tidyverse	Posit Software, PBC	v2.0.0 (CRAN)	Основной набор инструментов для манипулирования и форматирования данных