stxBrain com RCTDEsta vinheta demonstra como estimar misturas de tipos celulares na
seção anterior do cérebro de camundongo stxBrain com o
método Robust Cell Type Decomposition (RCTD), usando o Seurat 5.5.1 e a
versão 2.2.1 do pacote spacexr disponível no GitHub (Cable et al.
2022).
O objeto espacial será obtido diretamente pelo pacote
SeuratData, usando a seção anterior1 (Hao et al.
2024). Nenhum objeto espacial previamente processado será
reutilizado. Essa separação é importante porque o RCTD modela contagens
brutas de expressão e, por padrão, exige valores inteiros. Matrizes
normalizadas, como as produzidas por NormalizeData() ou
SCTransform(), não devem ser fornecidas como entrada do
método.
Ao final do roteiro, você será capaz de:
SeuratData as contagens brutas e as
coordenadas do Visium;Um spot do Visium é maior que uma única célula e pode capturar RNA proveniente de vários tipos celulares. O RCTD aprende perfis de expressão a partir de dados anotados de scRNA-seq (single-cell) ou snRNA-seq (single-nucleus), corrige diferenças entre as plataformas da referência e da amostra espacial e estima a contribuição de cada tipo celular em cada spot (Cable et al. 2022).
A análise possui duas entradas e uma saída principal:
| Componente | Objeto utilizado | Informação necessária |
|---|---|---|
| Amostra espacial | stxBrain, seção anterior1, carregada por
SeuratData |
contagens brutas de genes por spots e coordenadas dos spots |
| Referência single-cell | córtex de camundongo do Allen Institute | contagens brutas de genes por células e um rótulo de tipo celular por célula |
| Resultado do RCTD | objeto da classe RCTD |
classificação e pesos de até dois tipos celulares por spot |
O spacexr oferece três modos:
doublet: ajusta no máximo dois tipos celulares por
spot;multi: ajusta até MAX_MULTI_TYPES tipos
celulares por spot; efull: não impõe um limite fixo ao número de tipos
celulares.Usaremos o modo doublet, por ser mais leve. No modo
doublet, o RCTD atribui um ou dois tipos celulares a cada
spot. Isso é uma simplificação do modelo e não significa que um spot
Visium contenha, biologicamente, no máximo duas células ou dois tipos
celulares.
Como referência de desempenho, a execução validada para a etapa de
cálculo do algortimo do RCTD levou 14,22 minutos em um servidor de alta
performance com 160 GB de RAM e 30 núcleos de CPU disponíveis para
processamento paralelo. O tempo pode variar bastante de acordo com o
equipamento, a carga do servidor e o número de tipos celulares da
referência. O parâmetro max_cores controla o paralelismo na
CPU; esta análise não utiliza aceleração por GPU.
Além disso, a paralelização empregada neste fluxo está disponível
apenas em Linux e macOS. No Windows, max_cores deve ser
igual a 1, o que pode aumentar consideravelmente o tempo de
processamento. Diferentemente das demais, esta vinheta inclui figuras
previamente geradas para que o fluxo possa ser acompanhado mesmo quando
não for possível executar a análise completa durante o curso.
Durante o curso, recomenda-se executar os blocos de código na ordem em que aparecem. Assim, cada objeto intermediário pode ser conferido antes de prosseguir.
O conjunto stxBrain não precisa ser copiado para o
repositório: ele é instalado pelo SeuratData na biblioteca
de pacotes do R. O único arquivo de entrada local é
allen_cortex.rds, que será baixado mais adiante e
armazenado em data/. Os resultados serão gravados em
resultados/ e as figuras em imagens/.
Curso_ST/
|-- data/
| `-- allen_cortex.rds
|-- imagens/
`-- resultados/
A instalação de todos os pacotes e de suas dependências está
centralizada em Configuracao_do_ambiente_RStudio.Rmd.
O roteiro pressupõe esse ambiente recém-instalado e adota o Seurat 5.5.1
como referência, sem caminhos alternativos para versões anteriores. Para
a deconvolução, utiliza a versão 2.2.1 do repositório original do
spacexr no GitHub, cuja interface oferece
create.RCTD() e run.RCTD() e não é
intercambiável com a interface do Bioconductor.
pacotes_rctd <- c(
"Seurat", "SeuratObject", "SeuratData", "Matrix", "ggplot2", "spacexr"
)
pacotes_rctd_ausentes <- pacotes_rctd[
!vapply(pacotes_rctd, requireNamespace, logical(1), quietly = TRUE)
]
if (length(pacotes_rctd_ausentes) > 0) {
stop(
"Execute primeiro Configuracao_do_ambiente_RStudio.Rmd. ",
"Pacotes ausentes: ",
paste(pacotes_rctd_ausentes, collapse = ", ")
)
}
Para que cada bloco de código possa ser executado separadamente, o
fluxo usa a pasta de trabalho atual como ponto de partida. Antes de
continuar, confirme que ela corresponde à pasta vinhetas/
do curso.
Siga estes passos:
getwd() no Console. O R deverá mostrar um
caminho terminado em Curso_ST/vinhetas.diretorio_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_dados, showWarnings = FALSE, recursive = TRUE)
dir.create(diretorio_imagens, showWarnings = FALSE, recursive = TRUE)
dir.create(diretorio_resultados, showWarnings = FALSE, recursive = TRUE)
diretorio_curso
diretorio_dados
diretorio_imagens
diretorio_resultados
library("Seurat")
library("SeuratData")
library("Matrix")
library("ggplot2")
library("spacexr")
options(stringsAsFactors = FALSE)
Confira as versões carregadas. Os resultados esperados são Seurat
5.5.1 e spacexr 2.2.1 ou
maior:
packageVersion("Seurat")
packageVersion("spacexr")
stxBrain diretamente do SeuratDataNesta vinheta, o objeto Seurat stxBrain será carregado
novamente porque o RCTD precisa das contagens brutas, e não do objeto
normalizado produzido nas etapas anteriores do curso. O bloco abaixo
apenas confirma que o conjunto de dados já está instalado. Caso ele
ainda não esteja disponível, retorne à vinheta Configuracao_do_ambiente_RStudio.Rmd.
if (!requireNamespace("stxBrain.SeuratData", quietly = TRUE)) {
stop(
"O conjunto stxBrain não está instalado.\nExecute primeiro a vinheta ",
"Configuracao_do_ambiente_RStudio.Rmd: \n",
"https://flaviogckessler.github.io/Curso_ST/vinhetas/Configuracao_do_ambiente_RStudio.html"
)
}
message("O conjunto stxBrain já foi instalado.")
Com o objeto instalado, basta realizar o seu carregamento:
stxbrain <- SeuratData::LoadData("stxBrain", type = "anterior1")
stxbrain
O RCTD receberá somente a camada counts do ensaio
Spatial, que conserva as contagens inteiras e não
normalizadas. As camadas data e scale.data,
assim como um ensaio SCT, não devem ser usadas como
entrada.
contagens_espaciais <- SeuratObject::LayerData(
object = stxbrain,
assay = "Spatial",
layer = "counts"
)
print(paste("Número de spots:",dim(contagens_espaciais)[2]))
print(paste("Número de genes:",dim(contagens_espaciais)[1]))
contagens_espaciais[
seq_len(min(5L, nrow(contagens_espaciais))),
seq_len(min(5L, ncol(contagens_espaciais))),
drop = FALSE]
O total de transcritos varia entre os spots. Este resumo permite examinar a profundidade de sequenciamento e a ordem de grandeza das bibliotecas antes da deconvolução.
summary(Matrix::colSums(contagens_espaciais))
No Seurat 5.5.1, a imagem espacial do objeto stxBrain
fornece as coordenadas não escaladas nas colunas x e
y. O identificador de cada spot permanece nos nomes das
linhas, permitindo alinhar diretamente a tabela de coordenadas às
colunas da matriz de contagens antes de criar o objeto do
spacexr.
nome_imagem <- SeuratObject::Images(stxbrain)[1]
coordenadas_espaciais <- SeuratObject::GetTissueCoordinates(
object = stxbrain,
image = nome_imagem,
scale = NULL
)
coordenadas_espaciais <- as.data.frame(
coordenadas_espaciais[, c("x", "y"), drop = FALSE]
)
coordenadas_espaciais <- coordenadas_espaciais[
colnames(contagens_espaciais), , drop = FALSE
]
head(coordenadas_espaciais)
O RCTD requer uma referência de RNA-seq biologicamente adequada e anotada. Usaremos a taxonomia celular do córtex adulto de camundongo do Allen Institute, também empregada no tutorial espacial do Seurat (Tasic et al. 2016). Como a referência e a amostra espacial são de camundongo, os genes podem ser associados por símbolo gênico.
Baixe a referência single-cell com o código abaixo. O arquivo será
salvo como data/allen_cortex.rds. Atenção:
o download possui cerca de 1 GB e precisa ser realizado
apenas uma vez.
arquivo_referencia <- file.path(diretorio_dados, "allen_cortex.rds")
download.file(
url = paste0(
"https://www.dropbox.com/scl/fi/sgmfzk72cdoc945ap19c9/",
"allen_cortex.rds?rlkey=qfou47stxczo14m37vmea54wl&dl=1"
),
destfile = arquivo_referencia,
mode = "wb"
)
A referência publicada armazena os rótulos de tipos celulares na
coluna subclass.
Cuidado: o carregamento do arquivo de referência, células corticais de camundongo disponibilizadas pelo Allen Institute, é realizado no código abaixo. O arquivo é grande e vai ocupar 3GB de memória quando for carregado para o ambiente de trabalho.
arquivo_referencia <- file.path(diretorio_dados, "allen_cortex.rds")
# Confere se arquivo existe
if (!file.exists(arquivo_referencia)) {
stop(
"A referência single-cell não foi encontrada em: ",
arquivo_referencia,
"\nBaixe o arquivo como descrito acima ou copie o arquivo fornecido pelo instrutor para data/."
)
}
# Carrega o arquivo de referência
referencia_allen <- readRDS(arquivo_referencia)
#Extrai metadados das células do arquivo de referência
metadados_referencia <- referencia_allen[[]]
# Extrai contagens do arquivo de referência
contagens_referencia <- SeuratObject::LayerData(
object = referencia_allen,
assay = "RNA",
layer = "counts"
)
head(contagens_referencia[,1:5])
# Organiza os metadados para seguir a mesma ordem da matriz de contagem
metadados_referencia <- metadados_referencia[
colnames(contagens_referencia), , drop = FALSE
]
head(metadados_referencia)
# Extrai coluna dos tipos celulares dos metadados
tipos_celulares_referencia <- as.character(
metadados_referencia$subclass
)
# Corrige strings dos nomes das células
tipos_celulares_referencia <- gsub(
"/", "-", tipos_celulares_referencia, fixed = TRUE
)
celulas_validas <- !is.na(tipos_celulares_referencia) &
nzchar(tipos_celulares_referencia)
tamanhos_tipos <- table(tipos_celulares_referencia[celulas_validas])
# Filtra tipos celulares que estão presentes em menos de 25 células
tipos_elegiveis <- names(tamanhos_tipos[tamanhos_tipos >= 25])
celulas_validas <- celulas_validas &
tipos_celulares_referencia %in% tipos_elegiveis
# Retira células com tipo celular raro da contagem de referência
contagens_referencia <- contagens_referencia[
, celulas_validas, drop = FALSE
]
tipos_celulares_referencia <- droplevels(
factor(tipos_celulares_referencia[celulas_validas])
)
print(paste("N° de células de referência: ",ncol(contagens_referencia)))
print("Células de referência:")
levels(tipos_celulares_referencia)
O filtro padrão do RCTD exige pelo menos 25 células por tipo celular. Os tipos raros são removidos explicitamente acima, evitando que desapareçam de forma silenciosa durante o pré-processamento.
Vamos analisar os dados gerais que iremos trabalhar
genes_comuns <- intersect(
rownames(contagens_espaciais),
rownames(contagens_referencia)
)
resumo_genes <- data.frame(
fonte = c("amostra espacial", "referência single-cell", "interseção"),
numero_de_genes = c(
nrow(contagens_espaciais),
nrow(contagens_referencia),
length(genes_comuns)
)
)
resumo_genes
if (length(genes_comuns) < 1000L) {
stop(
"Foram encontrados menos de 1.000 genes compartilhados. Verifique o organismo, ",
"os identificadores gênicos e se ambas as matrizes contêm contagens brutas."
)
}
contagens_espaciais <- contagens_espaciais[
genes_comuns, , drop = FALSE
]
contagens_referencia <- contagens_referencia[
genes_comuns, , drop = FALSE
]
# Retirar dados pesados do ambiente de trabalho (liberar RAM)
rm(resumo_genes, referencia_allen)
invisible(gc())
O dado espacial (stxBrain) possui 31053 genes
identificados. A referência single-cell com os tipos celulares anotados
possui 34617 genes identificados. O dado espacial e a referência possuem
uma intersecção de 19293 genes identificados.
spacexrO spacexr recebe a referência single-cell em um objeto
Reference e a amostra espacial em um objeto
SpatialRNA. Em ambos, as somas por coluna representam o
número total de UMIs de cada célula ou spot.
names(tipos_celulares_referencia) <- colnames(contagens_referencia)
# Construir objeto do pacote spacexr para o dado da referência
referencia_rctd <- spacexr::Reference(
contagens_referencia,
tipos_celulares_referencia,
Matrix::colSums(contagens_referencia)
)
# Construir objeto do pacote spacexr para o dado espacial
consulta_rctd <- spacexr::SpatialRNA(
as.data.frame(coordenadas_espaciais),
contagens_espaciais,
Matrix::colSums(contagens_espaciais)
)
data.frame(objeto = c("referencia_rctd", "consulta_rctd"),
genes = c(nrow(referencia_rctd@counts),nrow(consulta_rctd@counts)),
unidades = c(ncol(referencia_rctd@counts),ncol(consulta_rctd@counts)),
tipo_de_unidade = c("células", "spots"))
# Retirar dados pesados do ambiente de trabalho (liberar RAM)
rm(contagens_espaciais,contagens_referencia,coordenadas_espaciais,
metadados_referencia,tipos_celulares_referencia,genes_comuns)
invisible(gc())
create.RCTD() prepara os dados, seleciona genes
informativos e aprende os perfis dos tipos celulares da referência. No
código abaixo foi atribuído 1 núcleo de processamento para
evitar exaurir o seu computador. Contudo, isso irá deixar a tarefa muito
demorado. Consulte o número de cores que o CPU do seu computador possui
e assine um número que seja seguro. A parelelização está apenas
disponíveis em sistemas Linux e macOS. No Windows, use obrigatoriamente
1.
set.seed(20260904)
numero_nucleos <- 1
RCTD <- spacexr::create.RCTD(
consulta_rctd,
referencia_rctd,
max_cores = numero_nucleos
)
rm(consulta_rctd, referencia_rctd)
invisible(gc())
O modo doublet procura um ou dois tipos celulares por
spot. O tempo de execução varia de acordo com o servidor e o número de
células da referência.
Cuidado: Essa é a etapa mais demorada da análise!
inicio_rctd <- Sys.time()
RCTD <- spacexr::run.RCTD(
RCTD,
doublet_mode = "doublet"
)
fim_rctd <- Sys.time()
tempo_rctd <- fim_rctd - inicio_rctd
tempo_rctd
saveRDS(
RCTD,
file = file.path(
diretorio_resultados,
"stxbrain_RCTD_doublet_results.rds"
)
)
rm(inicio_rctd, fim_rctd, tempo_rctd, numero_nucleos)
invisible(gc())
No modo doublet, os resultados principais ficam em dois
componentes do objeto RCTD:
results_df: classificação do spot e tipos celulares
atribuídos;weights_doublet: pesos associados ao primeiro e ao
segundo tipo.resultado_spots <- RCTD@results$results_df
pesos_doublet <- as.data.frame(RCTD@results$weights_doublet)
colnames(pesos_doublet) <- c("peso_primeiro_tipo", "peso_segundo_tipo")
resumo_spots_rctd <- cbind(resultado_spots, pesos_doublet)
head(resumo_spots_rctd)
table(resumo_spots_rctd$spot_class)
As classes produzidas pelo RCTD são:
singlet: um tipo celular foi atribuído com
confiança;doublet_certain: dois tipos celulares foram atribuídos
com confiança;doublet_uncertain: somente o primeiro tipo é confiável;
ereject: o método não produziu uma atribuição
confiável.Vamos acrescentar a tabela do RCTD ao metadado do
stxbrain. O Seurat utiliza os códigos de barras presentes
nos nomes das linhas para posicionar cada resultado no spot correto.
stxbrain <- SeuratObject::AddMetaData(
object = stxbrain,
metadata = resumo_spots_rctd
)
rm(RCTD, resultado_spots, pesos_doublet)
invisible(gc())
Primeiro, observe onde o RCTD identificou singlet,
doublet_certain, doublet_uncertain ou
reject. Ao executar os blocos manualmente, todas as figuras
serão salvas como PNG em imagens/, com resolução moderada
de 120 dpi.
p_classes <- Seurat::SpatialDimPlot(
object = stxbrain,
images = nome_imagem,
group.by = "spot_class",
pt.size.factor = 1.6,
crop = TRUE
) +
labs(title = "Classificação dos spots pelo RCTD")
p_classes
ggplot2::ggsave(
filename = file.path(
diretorio_imagens,
"rctd_classificacao_spots.png"
),
plot = p_classes,
width = 7,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)
rm(p_classes)
invisible(gc())
A figura abaixo foi produzida durante a execução manual validada da vinheta e é apresentada aqui como exemplo. Ela é incorporada ao HTML a partir do arquivo PNG já salvo, sem executar novamente o código.
Em seguida, visualize o primeiro e o segundo tipo celular atribuídos.
O segundo tipo deve ser interpretado somente junto com
spot_class: ele é confiável quando o spot foi classificado
como doublet_certain.
p1 <- Seurat::SpatialDimPlot(
object = stxbrain,
images = nome_imagem,
group.by = "first_type",
pt.size.factor = 1.6,
crop = TRUE
) +
labs(title = "Primeiro tipo celular")
p2 <- Seurat::SpatialDimPlot(
object = stxbrain,
images = nome_imagem,
group.by = "second_type",
pt.size.factor = 1.6,
crop = TRUE
) +
labs(title = "Segundo tipo celular")
p_tipos <- p1 + p2
p_tipos
ggplot2::ggsave(
filename = file.path(
diretorio_imagens,
"rctd_primeiro_tipo_celular.png"
),
plot = p1,
width = 7,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)
ggplot2::ggsave(
filename = file.path(
diretorio_imagens,
"rctd_segundo_tipo_celular.png"
),
plot = p2,
width = 7,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)
ggplot2::ggsave(
filename = file.path(
diretorio_imagens,
"rctd_tipos_celulares_combinados.png"
),
plot = p_tipos,
width = 10,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)
rm(p1, p2, p_tipos)
invisible(gc())
Os mapas individuais permitem examinar com mais detalhe a distribuição espacial de cada atribuição.
O painel combinado facilita a comparação direta entre o primeiro e o segundo tipo celular.
write.csv(
resumo_spots_rctd,
file = file.path(
diretorio_resultados,
"stxbrain_RCTD_doublet_spot_summary.csv"
),
row.names = TRUE
)
saveRDS(
stxbrain,
file = file.path(
diretorio_resultados,
"stxbrain_Seurat_com_RCTD.rds"
)
)
Os arquivos possuem finalidades complementares:
stxbrain_RCTD_doublet_results.rds: objeto completo
produzido pelo RCTD;stxbrain_RCTD_doublet_spot_summary.csv: classificação,
primeiro e segundo tipos celulares e seus pesos por spot; estxbrain_Seurat_com_RCTD.rds: objeto Seurat com os
resultados do RCTD no metadado.As figuras rctd_classificacao_spots.png,
rctd_primeiro_tipo_celular.png,
rctd_segundo_tipo_celular.png e
rctd_tipos_celulares_combinados.png serão gravadas em
imagens/.
Antes de formular conclusões biológicas, verifique:
reject como atribuições
confiáveis.doublet_uncertain, somente o primeiro tipo deve ser tratado
como confiável.As atribuições e os pesos do RCTD são estimativas baseadas em um modelo. Eles não correspondem diretamente ao número de células presentes no spot.
sessionInfo()