Nesta vinheta, usaremos o objeto Seurat processado na etapa anterior do curso para inferir e visualizar comunicação célula-célula em dados de transcriptômica espacial. O fluxo reúne duas abordagens complementares:
O objeto de entrada é obrigatoriamente data/brain.RDS,
gerado pela vinheta de controle de qualidade, pré-processamento e
clusterização. Embora ele contenha as seções anterior e posterior do
conjunto stxBrain, manteremos somente a região
anterior logo após o carregamento. Assim, todos os
cálculos, resultados e gráficos desta vinheta serão produzidos apenas
com essa região. O objeto também contém a expressão normalizada pelo
SCTransform() e os clusters do Seurat (Hao et al.
2024).
Ao final, a tabela de comunicações projetadas entre spots será salva como:
resultados/spot-spot_commu_anterior.rdsComo interpretar o resultado: o CellChat estima a comunicação no nível dos grupos, que nesta análise são os clusters do Seurat. A tabela spot–spot é uma projeção dessas estimativas sobre pares de spots espacialmente próximos e compatíveis com os clusters de origem e destino. Portanto, ela não representa uma inferência independente do CellChat para cada par individual de spots.
Antes de executar o código, defina o diretório de trabalho como a
pasta vinhetas/. No RStudio, use Session > Set
Working Directory > To Source File Location.
Curso_ST/
|-- data/
| `-- brain.RDS
|-- imagens/
|-- resultados/
`-- vinhetas/
|-- CellChat_Comunicacao_Celular_PathwaySpace.Rmd
`-- CellChat_Comunicacao_Celular_PathwaySpace.htmldiretorio_vinhetas <- getwd()
diretorio_curso <- dirname(diretorio_vinhetas)
diretorio_dados <- file.path(diretorio_curso, "data")
diretorio_imagens <- file.path(diretorio_curso, "imagens")
diretorio_resultados <- file.path(diretorio_curso, "resultados")
dir.create(diretorio_imagens, showWarnings = FALSE, recursive = TRUE)
dir.create(diretorio_resultados, showWarnings = FALSE, recursive = TRUE)A instalação de todos os pacotes e de suas dependências está
centralizada em Configuracao_do_ambiente_RStudio.Rmd.
Execute a primeira vinheta, reinicie a sessão do R e então use o bloco
abaixo apenas para confirmar que o ambiente está pronto.
pacotes_cellchat <- c(
"CellChat", "RGraphSpace", "PathwaySpace", "Seurat", "SeuratObject",
"dplyr", "tibble", "tidyr", "ggplot2", "RANN", "future"
)
pacotes_ausentes <- pacotes_cellchat[
!vapply(pacotes_cellchat, requireNamespace, logical(1), quietly = TRUE)
]
if (length(pacotes_ausentes) > 0) {
stop(
"Execute primeiro Configuracao_do_ambiente_RStudio.Rmd. ",
"Pacotes ausentes: ",
paste(pacotes_ausentes, collapse = ", ")
)
}Carregamos os pacotes usados na inferência, na manipulação dos dados e nas visualizações. Esta vinheta adota o Seurat 5.5.1 como referência.
library("CellChat")
library("RGraphSpace")
library("PathwaySpace")
library("Seurat")
library("SeuratObject")
library("dplyr")
library("tibble")
library("tidyr")
library("ggplot2")
options(stringsAsFactors = FALSE)
packageVersion("Seurat")Após carregar o objeto completo, o subset() descarta
todos os spots da região posterior. A imagem espacial correspondente
também é removida do objeto, mantendo apenas a seção
anterior1 para as etapas seguintes.
brain <- readRDS(file.path(diretorio_dados, "brain.RDS"))
brain <- subset(brain, subset = region == "anterior")
DefaultAssay(brain) <- "SCT"
Idents(brain) <- "seurat_clusters"
if (length(Images(brain)) != 1L) {
stop("O subset da região anterior deve manter exatamente uma imagem espacial.")
}
nome_imagem <- Images(brain)[1]
brain
table(brain$region)O resultado esperado é um objeto com 2.696 spots e uma única imagem,
denominada anterior1.
O CellChat requer quatro componentes principais para dados espaciais:
data.input: matriz normalizada de genes por spots;meta: rótulos dos clusters e das amostras;spatial.locs: coordenadas espaciais em resolução
completa; espatial.factors: fatores usados para converter pixels
em micrômetros.A camada normalizada do ensaio SCT fornece a expressão
usada pelo CellChat. Os metadados associam cada spot ao seu cluster do
Seurat e à única amostra espacial mantida.
data.input <- LayerData(brain, assay = "SCT", layer = "data")
spots_anterior <- Cells(brain[[nome_imagem]])
data.input <- data.input[, spots_anterior, drop = FALSE]
metadados_brain <- brain[[]]
meta <- data.frame(
labels = paste0("cluster_",
metadados_brain[colnames(data.input), "seurat_clusters"]),
samples = nome_imagem,
row.names = colnames(data.input)
)
meta$labels <- factor(meta$labels)
meta$samples <- factor(meta$samples, levels = nome_imagem)
table(meta$labels, meta$samples)No Seurat 5.5.1, GetTissueCoordinates() retorna as
posições não escaladas do VisiumV2 nas colunas
x e y. Depois de alinhar os spots às colunas
da matriz de expressão, estimamos o fator que converte pixels em
micrômetros. Para o Visium, adotamos 100 µm entre centros de spots
adjacentes e 55 µm como diâmetro teórico do spot.
spatial.locs <- GetTissueCoordinates(
brain, image = nome_imagem, scale = NULL)
spatial.locs <- as.data.frame(spatial.locs[, c("x", "y"), drop = FALSE])
spatial.locs <- spatial.locs[colnames(data.input), , drop = FALSE]
spot.size <- 55
distancia_centros <- 100
distancias <- computeCellDistance(coordinates = spatial.locs)
menor_distancia <- min(distancias[distancias > 0])
spatial.factors <- data.frame(
ratio = distancia_centros / menor_distancia,
tol = spot.size / 2,
row.names = nome_imagem
)
spatial.factors
head(spatial.locs)Com expressão, rótulos, coordenadas e escala na mesma ordem, construímos o objeto espacial que será usado na inferência.
cellchat <- createCellChat(
object = data.input, meta = meta, group.by = "labels",
datatype = "spatial", coordinates = spatial.locs,
spatial.factors = spatial.factors
)
cellchatComo o conjunto stxBrain é de camundongo, usaremos
CellChatDB.mouse. Para manter o tempo de execução adequado
à atividade didática, a análise será restrita à categoria
Secreted Signaling.
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)
CellChatDB.use <- subsetDB(CellChatDB,
search = "Secreted Signaling", key = "annotation")
cellchat@DB <- CellChatDB.useO CellChat mantém os genes presentes na base selecionada e identifica
genes e interações com expressão aumentada. A função plan()
com opção multisession permite distribuir essa etapa em até
quatro processos do R (cores).
cellchat <- subsetData(cellchat)
future::plan("multisession", workers = 4)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat, variable.both = FALSE)Nesta etapa, as interações e comunicações dependentes de contato foram restringidas a 100 µm, aproximadamente a distância entre os centros de spots adjacentes no Visium. Esse corte ajuda a acelerar a análise.
Para demonstração de aula, nboot foi limitado a 10 para
calcular mais rápido. nboot é número de testes de
permutação utilizado para calcular os p-valores.
Cuidado:* Etapa mais custosa da vinheta, a execução pode levar vários minutos.
cellchat <- computeCommunProb(
cellchat, type = "triMean", distance.use = TRUE,
interaction.range = 100, scale.distance = 0.01,
contact.range = 100, contact.dependent = TRUE,
contact.knn.k = 6, nboot = 10
)
cellchat <- filterCommunication(cellchat, min.cells = 10)Extraímos somente as interações com valor de p inferior a 0,05. A conversão dos grupos para texto facilitará as junções com as tabelas de spots.
df.net <- subsetCommunication(cellchat, thresh = 0.05)
df.net <- df.net %>%
mutate(source = as.character(source), target = as.character(target))
head(df.net)Agora identificaremos os pares de spots da região anterior separados por até 100 µm. Depois, cada par receberá os clusters de origem e destino e será combinado às interações inferidas pelo CellChat.
interaction.range <- 100
coordenadas_um <- spatial.locs * spatial.factors[nome_imagem, "ratio"]
matriz_distancias <- as.matrix(dist(coordenadas_um))
diag(matriz_distancias) <- Inf
indices <- which(matriz_distancias <= interaction.range, arr.ind = TRUE)
spot.pairs <- tibble(
source_spot = rownames(matriz_distancias)[indices[, 1]],
target_spot = colnames(matriz_distancias)[indices[, 2]],
distance_um = matriz_distancias[indices],
sample = nome_imagem
)
spot.groups <- tibble(spot = rownames(cellchat@meta),
group = as.character(cellchat@idents))
spot.pairs <- spot.pairs %>%
left_join(spot.groups, by = c("source_spot" = "spot")) %>%
rename(source = group) %>%
left_join(spot.groups, by = c("target_spot" = "spot")) %>%
rename(target = group)
head(spot.pairs)As interações inferidas para cada par de clusters são então associadas aos pares de spots compatíveis. A tabela resultante é ordenada e salva para reutilização.
spot.to.spot.communication <- spot.pairs %>%
inner_join(df.net, by = c("source", "target")) %>%
arrange(sample, source_spot, target_spot,
pathway_name, interaction_name)
head(spot.to.spot.communication)
saveRDS(
spot.to.spot.communication,
file = file.path(diretorio_resultados,
"spot-spot_commu_anterior.rds")
)O GraphSpace reúne as coordenadas, os metadados, a
expressão normalizada e a imagem da seção anterior. As coordenadas em
baixa resolução permitem sobrepor corretamente a rede à imagem do
tecido.
spots_visualizacao <- spots_anterior
coordenadas_gs <- GetTissueCoordinates(
brain, image = nome_imagem, scale = "lowres"
)
coordenadas_gs <- coordenadas_gs[
spots_visualizacao, c("x", "y"), drop = FALSE
]
metadados_gs <- brain[[]][spots_visualizacao, , drop = FALSE]
gs <- GraphSpace(cbind(coordenadas_gs, metadados_gs))
gs_fdata(gs) <- LayerData(
brain, assay = "SCT", layer = "data"
)[, spots_visualizacao, drop = FALSE]
gs_image(gs) <- GetImage(
brain[[nome_imagem]], mode = "raster"
)
gs <- normalizeGraphSpace(gs)Primeiro, conferimos a posição dos clusters sobre a imagem antes de acrescentar as relações de proximidade.
spatial_theme <- theme_gspace_coords(
theme = "th3", is_norm = TRUE,
xlab = "Coordenada espacial 1", ylab = "Coordenada espacial 2"
)
cpal1 <- DiscretePalette(nlevels(gs$seurat_clusters), "polychrome")
ggplot(gs) +
annotation_gspace_image(gs) +
geom_nodespace(mapping = aes(colour = seurat_clusters),
size = 1, pch = 19) +
scale_colour_manual(values = cpal1) +
theme_gspace_legend(discrete_colour = TRUE) +
spatial_themeA função abaixo conecta cada spot aos seis vizinhos mais próximos, respeitando o arranjo hexagonal esperado da lâmina Visium. Spots das bordas podem ter menos vizinhos.
proximityNetwork <- function(gs, nnei = 6) {
nodes <- gs@nodes
vizinhos <- RANN::nn2(nodes[, c("x", "y")],
nodes[, c("x", "y")], k = nnei + 1)
distancias <- as.data.frame(vizinhos$nn.dists[, 2:(nnei + 1)])
indices <- as.data.frame(vizinhos$nn.idx[, 2:(nnei + 1)])
distancia_tipica <- 1 / mean(1 / distancias$V1)
manter <- distancias
manter[, 1] <- TRUE
manter[, 2:ncol(manter)] <-
distancias[, 2:ncol(distancias)] >= 0.95 * distancia_tipica &
distancias[, 2:ncol(distancias)] <= 1.05 * distancia_tipica
matriz_arestas <- manter * indices
matriz_arestas$from <- row.names(matriz_arestas)
arestas <- matriz_arestas %>%
pivot_longer(cols = -from, names_to = "vizinho", values_to = "to") %>%
filter(to != 0) %>%
transmute(from = as.numeric(from), to = as.numeric(to),
v1 = pmin(from, to), v2 = pmax(from, to)) %>%
distinct(v1, v2) %>%
transmute(from = nodes$name[v1], to = nodes$name[v2])
gs_add_edges(gs) <- arestas
gs
}Antes de criar as arestas, ampliamos uma área central para observar apenas os spots. Essa visualização servirá de comparação para a rede construída no bloco seguinte.
ggplot(gs) +
annotation_gspace_image(gs) +
geom_edgespace() +
geom_nodespace(mapping = aes(colour = seurat_clusters),
size = 3, pch = 19) +
scale_colour_manual(values = cpal1) +
theme_gspace_legend(discrete_colour = TRUE) +
spatial_theme +
coord_cartesian(xlim = c(0.4, 0.6), ylim = c(0.4, 0.6))Agora aplicamos proximityNetwork() e repetimos o
recorte. As linhas mostram quais spots serão considerados vizinhos na
projeção das comunicações.
gs <- proximityNetwork(gs)
ggplot(gs) +
annotation_gspace_image(gs) +
geom_edgespace() +
geom_nodespace(mapping = aes(colour = seurat_clusters),
size = 3, pch = 19) +
scale_colour_manual(values = cpal1) +
theme_gspace_legend(discrete_colour = TRUE) +
spatial_theme +
coord_cartesian(xlim = c(0.4, 0.6), ylim = c(0.4, 0.6))As arestas do GraphSpace são não direcionadas, mas a
comunicação ligante–receptor possui uma origem e um destino. Por isso,
criaremos as duas orientações de cada aresta antes de associá-las às
comunicações do CellChat.
arestas <- gs@edges %>% select(name1, name2)
arestas_direcionadas <- bind_rows(
arestas,
arestas %>% transmute(name1 = name2, name2 = name1)
) %>%
distinct()
spotCom_dt <- arestas_direcionadas %>%
inner_join(
spot.to.spot.communication %>% filter(sample == nome_imagem),
by = c("name1" = "source_spot", "name2" = "target_spot")
)
spotCom_dt$arrowType <- 1
interacoes_disponiveis <- spotCom_dt %>%
select(ligand, receptor, interaction_name) %>%
filter(!is.na(ligand), !is.na(receptor)) %>%
distinct()
head(interacoes_disponiveis)A função communication_pathwayspace() mantém os spots
envolvidos em uma interação, constrói o objeto
PathwaySpace, soma a expressão do ligante e do receptor e
projeta esse sinal na rede. Nomes separados por sublinhado são tratados
como complexos formados por mais de um gene.
communication_pathwayspace <- function(gs, lig, recep, spotCom_dt) {
comunicacoes <- spotCom_dt %>% filter(ligand == lig, receptor == recep)
spots_selecionados <- unique(c(comunicacoes$name1, comunicacoes$name2))
gs_selecionado <- gs_subset_nodes(gs, spots_selecionados)
pspace_obj <- buildPathwaySpace(gs_selecionado)
distancia_spots <- getNearestNode(pspace_obj)
pdist <- mean(distancia_spots$dist)
pspace_obj <- silhouetteMapping(pspace_obj, baseline = 0.1)
resolver_features <- function(rotulo, features_disponiveis) {
if (rotulo %in% features_disponiveis) {
return(rotulo)
}
componentes <- strsplit(rotulo, "_", fixed = TRUE)[[1]]
features_disponiveis[match(tolower(componentes),
tolower(features_disponiveis))]
}
features_disponiveis <- colnames(pspace_obj@fdata)
features_ligante <- resolver_features(lig, features_disponiveis)
features_receptor <- resolver_features(recep, features_disponiveis)
if (anyNA(c(features_ligante, features_receptor))) {
stop("O ligante ou o receptor não foi encontrado nos dados de expressão.")
}
sinal_ligante <- rowSums(as.matrix(
pspace_obj@fdata[, features_ligante, drop = FALSE]))
sinal_receptor <- rowSums(as.matrix(
pspace_obj@fdata[, features_receptor, drop = FALSE]))
vertexSignal(pspace_obj) <- sinal_ligante + sinal_receptor
pspace_obj <- circularProjection(
pspace_obj, decay.fun = weibullDecay(decay = 0.5, pdist = pdist),
aggregate.fun = signalAggregation("wmean"),
k = gs_vcount(pspace_obj))
pspace_obj
}Usaremos a primeira interação disponível na região anterior. Se
preferir investigar um par específico, substitua
ligante_exemplo e receptor_exemplo pelos
símbolos desejados.
if (nrow(interacoes_disponiveis) == 0) {
stop("Nenhuma interação foi associada às arestas da região anterior.")
}
ligante_exemplo <- interacoes_disponiveis$ligand[1]
receptor_exemplo <- interacoes_disponiveis$receptor[1]
pspace_obj <- communication_pathwayspace(
gs = gs, lig = ligante_exemplo, recep = receptor_exemplo,
spotCom_dt = spotCom_dt)
grafico_exemplo <- ggplot(pspace_obj) +
annotation_gspace_image(pspace_obj) +
annotation_pspace_signal(pspace_obj, si.alpha = 0.25) +
spatial_theme + labs(
title = paste(ligante_exemplo, "-", receptor_exemplo))
grafico_exemploPara a atividade em sala, serão exportadas as três primeiras
interações. Para gerar todas as figuras, substitua a definição de
numero_interacoes por
nrow(interacoes_disponiveis).
numero_interacoes <- min(3, nrow(interacoes_disponiveis))
for (i in seq_len(numero_interacoes)) {
ligante <- interacoes_disponiveis$ligand[i]
receptor <- interacoes_disponiveis$receptor[i]
pspace_i <- communication_pathwayspace(
gs = gs,lig = ligante,recep = receptor,spotCom_dt = spotCom_dt)
grafico_i <- ggplot(pspace_i) +
annotation_gspace_image(pspace_i) +
annotation_pspace_signal(pspace_i, si.alpha = 0.25) +
spatial_theme + labs(title = paste(ligante, "-", receptor))
nome_arquivo <- paste0("CellChat_PathwaySpace_",nome_imagem,
"_",ligante,"-",receptor,".png")
nome_arquivo <- gsub("[^[:alnum:]_.-]","_",nome_arquivo)
ggsave(filename = file.path(diretorio_imagens, nome_arquivo),
plot = grafico_i,width = 9,height = 6,units = "in",dpi = 120,bg = "white")
rm(pspace_i, grafico_i)
gc()
}Ao concluir a vinheta, estarão disponíveis:
resultados/spot-spot_commu_anterior.rds: tabela de
comunicações projetadas entre spots da região anterior;imagens/CellChat_PathwaySpace_anterior1_<ligante>-<receptor>.png:
projeções das interações selecionadas; ecellchat, gs,
spot.to.spot.communication e spotCom_dt na
sessão do R.Repita essa análise de comunicação celular utilizando o dado produzido pelo RCTD
Em vez de utilizar a informação de clusters como nessa vinheta, utilize tipos celulares inferidos pela vinheta anterior de deconvolução.
Você pode utilizar o objeto Seurat com as anotações do RCTD salvo no
repositório na pasta resultados/. Você pode baixar em: stxbrain_Seurat_com_RCTD.rds
Por fim, registre as versões dos pacotes e do R usadas na execução.
sessionInfo()