Integração de múltiplas amostras de Single cell RNA-Seq

À medida que a complexidade dos dados de célula única aumenta, a integração de múltiplos conjuntos de dados tornou-se padrão. No entanto, os efeitos de lote — decorrentes de variações técnicas e biológicas — precisam ser corrigidos para uma análise precisa. Esses efeitos surgem de diferenças no manuseio das amostras, protocolos, plataformas de sequenciamento e fatores biológicos, como o histórico do doador ou a origem do tecido. Métodos computacionais ajudam a eliminar variações indesejadas, garantindo sinais biologicamente significativos. A correção de lote exige duas decisões principais: selecionar o método apropriado e seus parâmetros, e definir a covariável de lote com base no objetivo da integração. Neste notebook, exploramos conceitos centrais e métodos para integração de dados e correção de lote, com atividades práticas utilizando Seurat e Harmony. Além disso, realizamos benchmarking para comparar estratégias de integração, auxiliando na escolha do método mais eficaz, preservando ao mesmo tempo a relevância biológica.

Instalar as bibliotecas necessárias


download.file("https://github.com/eddelbuettel/r2u/raw/master/inst/scripts/add_cranapt_jammy.sh",
              "add_cranapt_jammy.sh")
Sys.chmod("add_cranapt_jammy.sh", "0755")
system("./add_cranapt_jammy.sh")
bspm::enable()
options(bspm.version.check=FALSE)

Criaremos uma função R para realizar todos os comandos necessários para seguir o notebook


shell_call <- function(command, ...) {
  result <- system(command, intern = TRUE, ...)
  cat(paste0(result, collapse = "\n"))
}

Instalar bibliotecas necessárias


# Install the R.utils package
install.packages("R.utils")

# Install a specific commit/version of the Seurat Wrappers package from GitHub
remotes::install_github('satijalab/seurat-wrappers@d28512f804d5fe05e6d68900ca9221020d52cf1d', upgrade = FALSE)

# Install the Seurat Data package from GitHub
remotes::install_github('satijalab/seurat-data')

# Check if the BiocManager package is installed; if not, install it quietly
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager", quiet = TRUE)

# Install the Harmony package for integrating single-cell datasets
install.packages("harmony")

Introdução

Integração é uma técnica usada para combinar dados de células individuais de diferentes experimentos ou condições em um único espaço de análise. Isso é especialmente útil quando dados de diferentes lotes (batches) ou condições podem ser integrados para identificar padrões biológicos comuns. Neste exemplo de integração de múltiplos conjuntos de dados representando, por exemplo, diferentes condições biológicas, manuseio de amostras ou protocolos experimentais, exploraremos a integração realizada por Seurat e Harmony.


# Carregar bibliotecas
library(Seurat)
library(SeuratData)
library(SeuratWrappers)
library(cowplot)
library(harmony)
# Se o tempo de execução for maior que 2000 segundos, a execução irá parar
options(timeout=2000)

Carregando dados

Here, we load the datasets necessary to perform an integration analysis.

Aqui, carregamos os conjuntos de dados necessários para realizar uma análise de integração. Este tutorial demonstra o processo de integração utilizando dados de células mononucleares do sangue periférico (peripheral blood mononuclear cells - PBMCs) Kang et al, 2017. O experimento são dois grupos distintos de PBMCs: um grupo controle e um grupo estimulado tratado com interferon beta. O tratamento induziu alterações na expressão gênica específicas nos tipos celulares, criando um desafio para a análise conjunta. Sem integração, as células podem se agrupar com base em sua condição de estimulação em vez de sua identidade intrínseca de tipo celular, complicando a interpretação dos resultados.

Para esclarecer, a integração visa alinhar os dados de ambos os grupos em um espaço comum, onde as células se agrupam por tipo, em vez de por condição de tratamento. Esse processo é essencial para identificar com precisão características compartilhadas e únicas entre os conjuntos de dados, minimizando a influência dos efeitos de lote introduzidos pelo delineamento experimental.

Etapas a seguir:

Carregar os conjuntos de dados de PBMCs: Utilize dados de expressão gênica pré-processados do estudo de Kang.

Inspecionar a composição do conjunto de dados: Verifique se os grupos controle e estimulado estão incluídos e verifique se há representação adequada do tipo de célula.

Realizar a integração: Utilize ferramentas de integração como Seurat ou Harmony para combinar os conjuntos de dados, considerando os efeitos de lote.

Visualizar os resultados: Examine o agrupamento para garantir que as células sejam agrupadas principalmente por tipo, e não por condição de tratamento.


# Instalar o dataset
InstallData("ifnb")


# Carregar o dataset
LoadData("ifnb")


# Transferir os dados de ifnb para o dado sc_dataset
sc_datasets <- ifnb


# %>%: Este é o operador do pacote dplyr. Ele permite encadear operações de forma legível, passando o resultado de uma operação como entrada para a próxima.
# glimpse é uma função do pacote dplyr que fornece uma visualização rápida e concisa do conjunto de dados, mostrando os primeiros valores de cada coluna e as respectivas classes.

sc_datasets %>% dplyr::glimpse()

Processamento dos dados

Nós pré-processamos os dados até um ponto que é comum para a integração de dados do Seurat e do Harmony


# Atualiza o objeto Seurat para o formato mais recente, corrigindo problemas de compatibilidade com objetos Seurat mais antigos
sc_datasets <- UpdateSeuratObject(object = sc_datasets)
# Se esta etapa for ignorada, ocorrerá o erro: 'nenhum slot de nome "images" para este objeto da classe "Seurat"

# Normaliza os dados de expressão genética no objeto Seurat
sc_datasets <- Seurat::NormalizeData(sc_datasets, verbose = FALSE)

A função NormalizeData do Seurat simplifica os dados brutos de scRNA-seq em um formato padronizado, facilitando a comparação de células e a detecção de padrões biológicos significativos. Esse processo envolve três etapas principais:

Escalonamente de Contagens: Para cada gene, a contagem em uma determinada célula é dividida pelo número total de contagens de todos os genes naquela célula. Isso ajusta as diferenças na profundidade do sequenciamento entre as células, garantindo comparações justas.
Ajuste de Escalonamento: O resultado da primeira etapa é multiplicado por um fator de escala (comumente definido como 10.000). Isso torna os valores normalizados mais interpretáveis, assemelhando-se aproximadamente a contagens por 10.000 leituras.
Transformação Logarítmica: Por fim, as contagens escalonadas são transformadas usando uma função logarítmica. Essa etapa reduz o efeito de valores extremamente altos ou baixos, facilitando a detecção de padrões nos dados.

Este método é essencial porque garante que quaisquer diferenças observadas entre as células reflitam variações biológicas verdadeiras, e não diferenças técnicas do processo de sequenciamento.


# Visualize a estrutura do objeto Seurat, incluindo metadados e dados de recursos
dplyr::glimpse(sc_datasets)

Execute a integração de dados do Seurat

A integração aqui realizada com Seurat consiste na utilização de uma Análise de Correlação Canônica para identificar âncoras entre conjuntos de dados, conforme indicado no original paper:

Dividir conjunto de dados em uma lista de conjuntos de dados


# Dividir (split) o objeto Seurat em uma lista de objetos Seurat menores com base na coluna de metadados "stim"
sc_datasets.list <- SplitObject(sc_datasets, split.by = "stim")

Encontre características altamente variáveis para cada condição separadamente


# Identifique os 2.000 principais recursos variáveis para cada conjunto de dados na lista
sc_datasets.list <- lapply(X = sc_datasets.list, FUN = function(x) {
                    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000) })

Encontre os recursos que são repetidamente variáveis em conjuntos de dados para integração (anchor features).


# Identificar recursos de integração a serem usados para alinhar conjuntos de dados
features <- SelectIntegrationFeatures(object.list = sc_datasets.list)

Encontre as âncoras de integração ou o conjunto final de genes altamente variáveis mais frequentemente selecionados nos lotes (batches)


# Encontre âncoras de integração entre os conjuntos de dados usando os recursos selecionados

sc_datasets.anchors <- FindIntegrationAnchors(object.list = sc_datasets.list, anchor.features = features)

# Abaixo está o resumo do processo em execução que aparecerá na sua tela

Integrar conjuntos de dados - Cria um ensaio de dados integrado


# Integre os conjuntos de dados em um único objeto Seurat usando as âncoras de integração
sc_datasets.combined <- IntegrateData(anchorset = sc_datasets.anchors)

# Abaixo está o resumo do processo em execução que aparecerá na sua tela


# Veja a estrutura do objeto Seurat integrado
dplyr::glimpse(sc_datasets.combined)

Execute o fluxo de trabalho padrão para visualização e agrupamento.


# Escalone os dados integrados, execute PCA, execute UMAP para redução de dimensionalidade, encontre vizinhos e agrupe células
sc_datasets.combined <- ScaleData(sc_datasets.combined, verbose = FALSE) %>%
                        RunPCA(npcs = 30, verbose = FALSE) %>%
                        RunUMAP(reduction = "pca", dims = 1:30) %>%
                        FindNeighbors(reduction = "pca", dims = 1:30) %>%
                        FindClusters(resolution = 0.5)


# Defina o tamanho dos gráficos para visualização
options(repr.plot.height = 5, repr.plot.width = 16)

# Visualização: Crie gráficos UMAP com diferentes opções de agrupamento ou rotulagem
p1 <- DimPlot(sc_datasets.combined, reduction = "umap", group.by = "stim")  # Agrupe células por metadados "stim"
p2 <- DimPlot(sc_datasets.combined, reduction = "umap", label = TRUE, repel = TRUE) # Rotule clusters no UMAP
p3 <- DimPlot(sc_datasets.combined, reduction = "umap", group.by = "seurat_annotations")  # Agrupe por anotações

p1 + p2 + p3 # Combine os gráficos

Execute Harmony

A partir do artigo original do Harmony: "O PCA incorpora células em um espaço com dimensionalidade reduzida. O Harmony aceita as coordenadas da célula neste espaço reduzido e executa um algoritmo iterativo para ajustar os efeitos específicos do conjunto de dados

O Harmony usa agrupamentos difusos para atribuir cada célula a múltiplos agrupamentos, enquanto um termo de penalidade garante que a diversidade dos conjuntos de dados dentro de cada agrupamento seja maximizada.

O Harmony calcula uma centralidade global para cada agrupamento, bem como centralidades específicos do conjunto de dados para cada agrupamento.

Dentro de cada agrupamento, o Harmony calcula um fator de correção para cada conjunto de dados com base nas centralidades.

Por fim, o Harmony corrige cada célula com um fator específico da célula: uma combinação linear de fatores de correção do conjunto de dados ponderados pelas atribuições suaves de agrupamento da célula feitas na etapa a. O Harmony repete as etapas de a a d até a convergência. A dependência entre a atribuição de clusters e o conjunto de dados diminui a cada rodada. Os conjuntos de dados são representados com cores e os tipos de células com formatos diferentes."


# Defina o tamanho dos gráficos para visualização
options(repr.plot.height = 4, repr.plot.width = 6)

# Encontre recursos variáveis (Variable Features), escalone dados, execute PCA e integre usando Harmony
sc_datasets.harmony <- FindVariableFeatures(sc_datasets, selection.method = "vst", nfeatures = 2000) %>%
                       ScaleData(verbose = FALSE) %>%
                       RunPCA(npcs = 30, verbose = FALSE) %>%
                       RunHarmony("stim", plot_convergence = TRUE)  # Integrar usando "stim" como uma variável de lote

Para acessar diretamente os embeddings harmonizados:


# Extraia embeddings Harmony (representações de baixa dimensão) e exiba as primeiras 5 linhas e colunas
harmony_embeddings <- Embeddings(sc_datasets.harmony, 'harmony')
harmony_embeddings[1:5, 1:5]


# Execute UMAP usando embeddings Harmony, encontre vizinhos e agrupe células
sc_datasets.harmony <- RunUMAP(sc_datasets.harmony, reduction = "harmony", dims = 1:30) %>%
                       FindNeighbors(reduction = "harmony", dims = 1:30) %>%
                       FindClusters(resolution = 0.5)


#  Defina o tamanho dos gráficos para visualização
options(repr.plot.height = 5, repr.plot.width = 16)

# Visualização: Crie gráficos UMAP com diferentes opções de agrupamento ou rotulagem
p1 <- DimPlot(sc_datasets.harmony, reduction = "umap", group.by = "stim")  # Agrupar por "stim"
p2 <- DimPlot(sc_datasets.harmony, reduction = "umap", label = TRUE, repel = TRUE)  # Rotular de agrupamentos (clusters)
p3 <- DimPlot(sc_datasets.harmony, reduction = "umap", group.by = "seurat_annotations")  # Agrupar por anotação
p1 + p2 + p3  # Combine os gráficos

Métricas de Integração

Inspecionar a mistura dentro do cluster

Usaremos uma função para criar um gráfico da composição do lote de cada cluster.

Esta função foi retirada de: https://github.com/cellgeni/scRNA.seq.course/blob/master/course_files/utils/custom_seurat_functions.R


# Defina uma função personalizada para visualizar a distribuição de clusters em conjuntos de dados
plot_integrated_clusters = function(srat, batchcolumn) {
# Pegue um objeto Seurat integrado e plote as distribuições na "batchcolumn"
  library(Seurat)
  library(patchwork)
  library(ggplot2)
  library(reshape2)
  library(RColorBrewer)

  # Cria uma tabela de contagens de células por cluster e conjunto de dados
  count_table <- table(srat@meta.data$seurat_clusters, srat@meta.data[[batchcolumn]])
  count_mtx   <- as.data.frame.matrix(count_table)  # Converte a tabela em um dataframe
  count_mtx$cluster <- rownames(count_mtx)  # Adiciona identificadores de cluster como uma coluna
  melt_mtx    <- melt(count_mtx)  # Remodela os dados para ggplot
  melt_mtx$cluster <- as.factor(melt_mtx$cluster)  # Converte a coluna do cluster em um fator

  # Calcule o tamanho dos agrupamentos (clusters)
  cluster_size   <- aggregate(value ~ cluster, data = melt_mtx, FUN = sum)

  # Organize os agrupamentos por tamanho
  sorted_labels <- paste(sort(as.integer(levels(cluster_size$cluster)), decreasing = TRUE))
  cluster_size$cluster <- factor(cluster_size$cluster, levels = sorted_labels)
  melt_mtx$cluster <- factor(melt_mtx$cluster, levels = sorted_labels)

  colnames(melt_mtx)[2] <- "dataset" # Renomear coluna para maior clareza

 # Trace a distribuição de células por cluster (escala logarítmica)
  p1 <- ggplot(cluster_size, aes(y = cluster, x = value)) +
        geom_bar(position = "dodge", stat = "identity", fill = "grey60") +
        theme_bw() + scale_x_log10() + xlab("Cells per cluster, log10 scale") + ylab("")

  # Trace uma fração de células em cada conjunto de dados para cada agrupamento
  p2 <- ggplot(melt_mtx, aes(x = cluster, y = value, fill = dataset)) +
        geom_bar(position = "fill", stat = "identity") +
        theme_bw() + coord_flip() +
        scale_fill_brewer(palette = "Set2") +
        ylab("Fraction of cells in each dataset") + xlab("Cluster number") + theme(legend.position = "top")

  p2 + p1 + plot_layout(widths = c(3, 1))  # Combine os gráficos
}


# Visualizar agrupamentos para conjuntos de dados integrados
plot_integrated_clusters(sc_datasets.combined, 'stim')


# Visualizar agrupamentos pelo Harmony
plot_integrated_clusters(sc_datasets.harmony, 'stim')

Isso mede o quão bem misturado é um conjunto de dados composto. De acordo comeste link, pontuações mais baixas representam uma melhor mistura. Veja também este código

Esta métrica avalia se a vizinhança de uma célula está bem misturada. Em outras palavras, se ela contém um pequeno número de células de cada conjunto de dados (por exemplo, k = 5).


# Calcular a métrica de mistura para dados integrados de Seurat
seurat_mixing <- MixingMetric(sc_datasets.combined,
                              'seurat_clusters', # Coluna para avaliar a mistura
                              reduction = "pca", # Usar embeddings de PCA
                              dims = 1:2,        # Avaliar as duas primeiras dimensões
                              k = 5,             # Número de vizinhos mais próximos
                              max.k = 300,       # Número máximo de vizinhos
                              eps = 0,           # Tolerância para busca de vizinhos
                              verbose = TRUE     # Exibir progresso
                             )


# Calcular a média da métrica de mistura de Seurat
mean(seurat_mixing)


# Calcular o desvio padrão da métrica de mistura de Seurat
sd(seurat_mixing)


# Calcular a métrica de mistura para dados integrados do Harmony
harmony_mixing <- MixingMetric(sc_datasets.harmony,
                               'seurat_clusters',
                               reduction = "harmony",
                               dims = 1:2,
                               k = 5,
                               max.k = 300,
                               eps = 0,
                               verbose = TRUE
                             )


# Calcular a média da métrica de mistura Harmony
mean(harmony_mixing)


# Calcular o desvio padrão da métrica de mistura Harmony
sd(harmony_mixing)

Vamos comparar quando usamos nenhuma integração, Seurat e Harmony:


# Identifique as 2.000 características mais variáveis (genes) usando o método "vst".
# Essas características são usadas para análise a jusante (downstream analysis)
sc_datasets <- FindVariableFeatures(sc_datasets, selection.method = "vst", nfeatures = 2000) %>%
                
# Escalone e centralize os dados para as características variáveis identificadas.
# Isso padroniza os valores da expressão entre as células.
                ScaleData(verbose = FALSE) %>%
                
# Performe Principal Component Analysis (PCA) para reduzir a dimensionalidade.
# Retenha os primeiros 30 componentes principais para análise a jusante (downstream analysis).
                RunPCA(npcs = 30, verbose = FALSE) %>%

## Aplique a Aproximação e Projeção Uniforme de Variedades (UMAP) para redução da dimensionalidade não linear.
# Use os resultados do PCA como arquivo de entrada (input) e os primeiros 30 componentes principais.
                RunUMAP(reduction = "pca", dims = 1:30) %>%
                
# Calcula um grafo de vizinho mais próximo compartilhado (shared nearest neighbor - SNN) para agrupamento.
# Esta etapa identifica relacionamentos entre células com base em seus embeddings UMAP.
                FindNeighbors(reduction = "pca", dims = 1:30) %>%
                
# Executa o agrupamento de células usando o grafo SNN.
# A resolução controla a granularidade dos clusters (maior resolução = mais clusters).
                FindClusters(resolution = 0.5)


# Defina o tamanho da visualização de saída.
options(repr.plot.height = 5, repr.plot.width = 16)

# Crie um gráfico UMAP agrupando as células pela coluna de metadados "stim".
p1 <- DimPlot(sc_datasets, reduction = "umap", group.by = "stim")

# Crie um gráfico UMAP para o conjunto de dados integrado, agrupando as células pela coluna de metadados "stim".
p2 <- DimPlot(sc_datasets.combined, reduction = "umap", group.by = "stim")

# Crie um gráfico UMAP para o conjunto de dados corrigido por Harmonia, agrupando as células pela coluna de metadados "stim".
p3 <- DimPlot(sc_datasets.harmony, reduction = "umap", group.by = "stim")

# Combine os três gráficos UMAP em uma única visualização.
p1 + p2 + p3


# Defina a altura e a largura da visualização de saída.
options(repr.plot.height = 5, repr.plot.width = 16)

# Crie um gráfico UMAP para o objeto sc_datasets original.
# O parâmetro 'label' adiciona rótulos aos clusters no gráfico.
# O parâmetro 'repel' impede a sobreposição de rótulos de texto.
p1 <- DimPlot(sc_datasets, reduction = "umap", label = TRUE, repel = TRUE)

# Crie um gráfico UMAP para o objeto sc_datasets combinado.
p2 <- DimPlot(sc_datasets.combined, reduction = "umap", label = TRUE, repel = TRUE)

# Crie um gráfico UMAP para o objeto sc_datasets com a harmonia corrigida.
# Este gráfico também usa UMAP, com rótulos e repeling para evitar sobreposições.
p3 <- DimPlot(sc_datasets.harmony, reduction = "umap", label = TRUE, repel = TRUE)

p1 + p2 + p3 #Combine os gráficos

Perguntas

Qual a diferença entre os resultados da integração de Seurat e Harmony?

Por que isso pode acontecer?

Você consegue identificar quais marcadores mudaram nos clusters gerados por um método em comparação ao outro?