Nesta vinheta, duas seções sagitais de cérebro de camundongo
produzidas com a tecnologia Visium são carregadas do conjunto
stxBrain, reunidas em um objeto Seurat e submetidas a um
fluxo direto de:
SCTransform();Ao final, o objeto processado será salvo como
data/brain.RDS. Todas as figuras serão gravadas em formato
PNG na pasta imagens/.
A versão HTML desta vinheta apresenta o roteiro sem executar a
análise. Durante o curso, abra o arquivo .Rmd no RStudio e
execute os blocos na ordem em que aparecem.
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. Em seguida,
execute:
diretorio_vinhetas <- getwd()
diretorio_curso <- dirname(diretorio_vinhetas)
diretorio_dados <- file.path(diretorio_curso, "data")
diretorio_imagens <- file.path(diretorio_curso, "imagens")
dir.create(diretorio_dados, showWarnings = FALSE, recursive = TRUE)
dir.create(diretorio_imagens, showWarnings = FALSE, recursive = TRUE)library("Seurat")
library("SeuratData")
library("sctransform")
library("ggplot2")
library("patchwork")O conjunto stxBrain e as demais dependências são
instalados exclusivamente em
Configuracao_do_ambiente_RStudio.Rmd. Depois de executar a
primeira vinheta e reiniciar a sessão do R, use o bloco abaixo apenas
para confirmar que o conjunto está disponível.
if (!requireNamespace("stxBrain.SeuratData", quietly = TRUE)) {
stop(
"O conjunto stxBrain não está instalado. Execute primeiro ",
"Configuracao_do_ambiente_RStudio.Rmd."
)
}
message("O conjunto stxBrain está instalado.")Serão utilizadas uma seção anterior e uma seção posterior do cérebro de camundongo.
brain1 <- SeuratData::LoadData("stxBrain", type = "anterior1")
brain2 <- SeuratData::LoadData("stxBrain", type = "posterior1")
brain <- merge(brain1, brain2)
brainA coluna region, incluída no conjunto de dados,
identifica a região anatômica associada a cada spot. Com o código abaixo
é possível identificar o número de spots em cada região (2696 na
anterior e 3353 na posterior)
table(brain$region)Número de features (genes) identificados no objeto brain (anterior+posterior):
dim(brain)[1]
#31053As métricas principais são:
nCount_Spatial: número total de transcritos detectados
no spot;nFeature_Spatial: número de genes detectados; epercent.mt: porcentagem de transcritos provenientes de
genes mitocondriais.brain[["percent.mt"]] <- Seurat::PercentageFeatureSet(
brain,
pattern = "^mt-"
)vln_contagens <- VlnPlot(
brain,
features = "nCount_Spatial",
pt.size = 0
) + NoLegend()
mapa_contagens_anterior <- SpatialFeaturePlot(
brain,
features = "nCount_Spatial",
images = "anterior1"
)
mapa_contagens_posterior <- SpatialFeaturePlot(
brain,
features = "nCount_Spatial",
images = "posterior1"
)
painel_contagens <- vln_contagens /
(mapa_contagens_anterior | mapa_contagens_posterior) +
plot_annotation(title = "Transcritos por spot")
painel_contagens
ggsave(
filename = file.path(diretorio_imagens, "Brain_SpotQC_contagens.png"),
plot = painel_contagens,
width = 11,
height = 7,
units = "in",
dpi = 120,
bg = "white"
)vln_genes <- VlnPlot(
brain,
features = "nFeature_Spatial",
pt.size = 0
) + NoLegend()
mapa_genes_anterior <- SpatialFeaturePlot(
brain,
features = "nFeature_Spatial",
images = "anterior1"
)
mapa_genes_posterior <- SpatialFeaturePlot(
brain,
features = "nFeature_Spatial",
images = "posterior1"
)
painel_genes <- vln_genes /
(mapa_genes_anterior | mapa_genes_posterior) +
plot_annotation(title = "Genes detectados por spot")
painel_genes
ggsave(
filename = file.path(diretorio_imagens, "Brain_SpotQC_genes.png"),
plot = painel_genes,
width = 11,
height = 7,
units = "in",
dpi = 120,
bg = "white"
)vln_mitocondrial <- VlnPlot(
brain,
features = "percent.mt",
pt.size = 0
) + NoLegend()
mapa_mt_anterior <- SpatialFeaturePlot(
brain,
features = "percent.mt",
images = "anterior1"
)
mapa_mt_posterior <- SpatialFeaturePlot(
brain,
features = "percent.mt",
images = "posterior1"
)
painel_mitocondrial <- vln_mitocondrial /
(mapa_mt_anterior | mapa_mt_posterior) +
plot_annotation(title = "Porcentagem mitocondrial por spot")
painel_mitocondrial
ggsave(
filename = file.path(diretorio_imagens, "Brain_SpotQC_mitocondrial.png"),
plot = painel_mitocondrial,
width = 11,
height = 7,
units = "in",
dpi = 120,
bg = "white"
)Neste conjunto didático, o perfil das métricas é utilizado para inspeção e nenhum limiar adicional de filtragem será aplicado.
Caso você queira filtrar de acordo com alguma parâmetro dos metadados, corra o código abaixo. Para essa vinheta, ele é OPCIONAL, então é possível pulá-lo. Nessa exemplo é aplicado um filtro que passam apenas spots que tenham mais de 200 genes identificados; que possuam mais de 500 counts; e que possuam mais de 5% de mitocôndria.
brain_filt <- subset(brain, subset = nFeature_Spatial > 200 &
nCount_Spatial > 500 &
percent.mt > 5)
print("ANTES DE FILTRAR")
print(paste("Número de spots:",dim(brain)[2]))
print(paste("Número de genes:",dim(brain)[1]))
table(brain$region)
print("DEPOIS DE FILTRAR")
print(paste("Número de spots:",dim(brain_filt)[2]))
print(paste("Número de genes:",dim(brain_filt)[1]))
table(brain_filt$region)Com esse exemplo, foram excluídos 11 spots da região anterior e 298
da região posterior. Para comparar os novos parâmetros, rode os gráficos
acima com o dado filtrado brain_filt.
SCTransform() normaliza os dados e estabiliza a
variância. Em seguida, os componentes principais são utilizados para
construir UMAP, t-SNE, grafo de vizinhança e clusters.
Essas etapas são as mais demoradas da vinheta, principalmente SCTransform(). Curiosidade: ao rodar, será mostrado o tempo que levou para rodar as análises no seu computador.
set.seed(20260906)
inicio <- Sys.time()
brain <- SCTransform(
brain,
assay = "Spatial",
do.scale = TRUE,
return.only.var.genes = FALSE,
verbose = FALSE
)
fim <- Sys.time()
tempo <- fim - inicio
tempo
inicio <- Sys.time()
brain <- RunPCA(brain, verbose = FALSE)
brain <- RunUMAP(brain, dims = 1:10)
brain <- RunTSNE(brain, dims = 1:10)
brain <- FindNeighbors(brain, dims = 1:10)
brain <- FindClusters(brain, verbose = FALSE)
fim <- Sys.time()
tempo <- fim - inicio
tempoO objeto processado é salvo na pasta data/, onde poderá
ser reutilizado nas análises seguintes do curso.
saveRDS(
brain,
file = file.path(diretorio_dados, "brain.RDS")
)umap_brain <- DimPlot(brain, reduction = "umap") +
labs(title = "UMAP")
tsne_brain <- DimPlot(brain, reduction = "tsne") +
labs(title = "t-SNE")
painel_reducoes <- umap_brain | tsne_brain
painel_reducoes
ggsave(
filename = file.path(diretorio_imagens, "Brain_clustering_reducoes.png"),
plot = painel_reducoes,
width = 11,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)mapa_clusters <- SpatialDimPlot(brain) +
plot_annotation(title = "Distribuição espacial dos clusters")
mapa_clusters
ggsave(
filename = file.path(diretorio_imagens, "Brain_clustering_espacial.png"),
plot = mapa_clusters,
width = 11,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)distribuicao_clusters <- as.data.frame(
table(brain$seurat_clusters, brain$orig.ident)
)
colnames(distribuicao_clusters) <- c("Cluster", "Secao", "Spots")
p_numero_spots <- ggplot(
distribuicao_clusters,
aes(x = Secao, y = Spots, fill = Cluster)
) +
geom_col() +
coord_flip() +
labs(x = "Seção", y = "Número de spots")
p_proporcao_spots <- ggplot(
distribuicao_clusters,
aes(x = Secao, y = Spots, fill = Cluster)
) +
geom_col(position = "fill") +
coord_flip() +
labs(x = "Seção", y = "Proporção de spots")
painel_composicao <- p_numero_spots | p_proporcao_spots
painel_composicao
ggsave(
filename = file.path(diretorio_imagens, "Brain_clustering_composicao.png"),
plot = painel_composicao,
width = 10,
height = 5,
units = "in",
dpi = 120,
bg = "white"
)Os genes abaixo destacam populações neuronais e não neuronais presentes no tecido. Sua distribuição ajuda a relacionar os clusters aos compartimentos anatômicos, mas não substitui uma anotação celular completa.
genes <- c(
"Hpca", "Nrgn", # neurônios
"Plp1", # oligodendrócitos
"Aif1", # micróglia
"Gfap", # astrócitos
"Ttr", # células do plexo coroide
"Mgp" # fibroblastos
)
painel_genes_marcadores <- SpatialFeaturePlot(
brain,
features = genes,
ncol = 2
)
painel_genes_marcadores
ggsave(
filename = file.path(diretorio_imagens, "Brain_Genes.png"),
plot = painel_genes_marcadores,
width = 12,
height = 18,
units = "in",
dpi = 120,
bg = "white",
limitsize = FALSE
)A função FindMarkers() do Seurat que
realiza um cálculo de expressão diferencial entre dois fatores
pré-anotados. Para esse exemplo, serão utilizados dois clusteres que
foram calculados anteriormente, a informação de clusteres fica
armazenada nos metadados do objeto Seurat. A função pode
ser utilizada para qualquer fator pré-anotado, isso inclui clusteres,
PCs, tipos celulares e afins.
Para prosseguir com o exemplo dessas análises, iremos trabalhar
apenas com região anterior. Para isso, será feito um subset
do objeto brain.
# subset da região
brain_ant <- subset(brain, subset = region == "anterior")
# Foram selecionados os clusteres 1 e 10 para a análise do FindMarkers
SpatialDimPlot(brain_ant,
cells.highlight = CellsByIdentities(object = brain, idents = c(1,10)),
facet.highlight = TRUE,ncol = 2)# Prepara os diferentes modelos SCT para expressão diferencial
brain_ant <- PrepSCTFindMarkers(
brain_ant,
assay = "SCT"
)
# Procura "marcadores" das regiões selecionadas
de_markers <- FindMarkers(
brain_ant,
ident.1 = "1",
ident.2 = "10",
assay = "SCT"
)
de_markersO data frame de_markers obtido elenca os genes
diferencialmente expressos entre os fatores fornecidos, nesse caso dois
clusteres. As colunas indicam:
p_val : p-valor
avg_log2FC : log fold-change da expressão média entre os dois “grupos”. Valores positivos indicam que o feature (gene) está mais expresso no primeiro grupo. Valores negativos indicam que o feature (gene) está mais expresso no segundo grupo.
pct.1 : O percentual de células em que o feature (gene) é detectado no primeiro grupo
pct.2 : O percentual de células em que o feature (gene) é detectado no segundo grupo
p_val_adj : p-valor ajustado, o ajuste é realizado com a correção de bonferroni usando todos os features do dataset.
Encontre mais informações sobre a função FindMarkers()
na vinheta do Seurat: Differential
expression testing.
g <- SpatialFeaturePlot(object = brain_ant,
features = rownames(de_markers)[1:3],
alpha = c(0.1, 1), ncol = 3)
g
# Salvar figura
ggsave(
filename = file.path(diretorio_imagens, "FindMarkers.png"),
plot = g,
width = 6,
height = 3,
units = "in",
dpi = 120,
bg = "white",
limitsize = FALSE
)A imagem gerada apresenta três genes diferencialmente expressos do cluster 10 em relação ao cluster 1, ou seja, genes que estão enriquecidos no cluster 10.
SVGs são genes que a distribuição da sua expressão gênica é dependente da localização. Em outras palavras, a distribuição da expressão de um SVG não é aleatória ao longo das coordenadas XY.
A função FindSpatiallyVariableFeatures() do
Seurat utiliza métodos estatísticas que calculam quanto
quea expressão de determinado gene é dependente da localização
(coordenadas XY). Abaixo, vamos utilizar o I de Moran, uma
estatística de pode ser entendida como a correlação da expressão em
relação à localizações próximas, também chamada de autocorrelação
espacial.
Para fins de demonstração, vamos calcular a estatística para apenas 100 genes. Você pode alterar o valor de acordo com a sua vontade, mas quanto maior for, mais tempo vai levar.
brain_ant_I <- FindSpatiallyVariableFeatures(brain_ant,
features = VariableFeatures(brain)[1:100],
selection.method = "moransi")# Seis maiores SVGs identificados
top.features <- head(SpatiallyVariableFeatures(brain_ant_I, method = "moransi"), 6)
# Apesar de análise ter sido feita apenas com a região anterior,
# podemos plotar o gráfico com as duas regiões.
g <- SpatialFeaturePlot(brain,
features = top.features, ncol = 2, alpha = c(0.1, 1)) &
theme(plot.margin = unit(c(1, 1, 1, 1), "mm"))
g
# Salvar figura
ggsave(
filename = file.path(diretorio_imagens, "SVGs.png"),
plot = g,
width = 10,
height = 12,
units = "in",
dpi = 120,
bg = "white",
limitsize = FALSE
)Note que nesse caso não se aplica o conceito de Gene Diferencialmente Expresso (DEG), pois não é uma comparação entre caso vs. controle. Ainda, não se aplica o conceito de Fold Change, por não ter uma referência ou por não ser uma comparação entre dois genes. Métodos que identificam SVGs descrevem a natureza da distribuição da expressão de um gene em relação ao espaço geográfico. Existem diversos métodos na literatura que identificam SVGs, como: I de Moran,SpatialDE,nnSVG, SPARK-X, MERINGUE.
Ao concluir a vinheta, os seguintes arquivos estarão disponíveis:
data/brain.RDS: objeto Seurat normalizado, reduzido e
clusterizado;imagens/Brain_SpotQC_contagens.png;imagens/Brain_SpotQC_genes.png;imagens/Brain_SpotQC_mitocondrial.png;imagens/Brain_clustering_reducoes.png;imagens/Brain_clustering_espacial.png;imagens/Brain_clustering_composicao.png; eimagens/Brain_Genes.png.sessionInfo()