Expressão diferencial com Sleuth - lmigueel/Bioinformatica GitHub Wiki
1. Sobre
Sleuth é uma ferramenta rápida e leve que usa a saída de estimativas de abundância transcrita de algoritmos de pseudoalinhamento que usam amostragem bootstrap, como Sailfish, Salmon e Kallisto, para realizar análises de expressão diferencial de isoformas de genes.
sleuth fornece ferramentas para análise exploratória de dados utilizando Shiny no RStudio e implementa algoritmos estatísticos para análise diferencial que potencializam as estimativas boostrap de kallisto. Ele também executa análises temporais com modelos LRTs, mas deixemos para complementaçoes.
Acesse o manual do sleuth aqui.
2. Instalação
O sleuth é um pacote em R, logo precisa de alguns pacotes complementares. Observe a versão do seu R para a instalação do pacote rhdf5.
#versoes antigas usando o BiocLite
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")
#Versoes novas R > 3.8
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install(version = "3.13") #ou 3.12
#BiocManager::install("rhdf5")
install.packages("devtools")
devtools::install_github("pachterlab/sleuth")
library("sleuth")
3. Explicação rápida
Sleuth modela diferentes fontes de variação para prever transcritos e genes diferencialmente expressos. A variação biológica (biol. Var.) proveniente das diferenças no conteúdo de RNA entre as réplicas e da bioquímica estocástica durante a preparação da biblioteca, enquanto a variação inferencial surge do sequenciamento aleatório e da análise computacional dos reads.
Sleuth depende da decomposição de variância para identificar diferenças biológicas na transcrição ou expressão gênica, enquanto usa uma estratégia padrão de shrinkage para estabilizar as estimativas de variância de poucas amostras, o sleuth é capaz de alavancar avanços recentes na quantificação para obter estimativas de erro que podem ser usadas para desacoplar a variação biológica da variação inferencial antes do shrinkage. A decomposição de variância é importante por causa da diversidade de estimativas de variância entre os genes que surgem ao quantificar abundâncias.
Em um exemplo no artigo, DESeq2 e voom foram executados nos mesmos dados com resumos de featureCounts, e essas ferramentas identificaram um gene como diferencialmente expresso a uma taxa de descoberta falsa (FDR) limite de 0,10 (FDR relatado 8,81e−21 e 5,56e−10, respectivamente); enquanto o sleuth não encontrou diferenças significativas entre as condições (FDR relatado 0,156) por causa da alta variância inferencial. Alguns métodos tentaram utilizar estimativas de erros de quantificação6, mas esses métodos são limitados por longos tempos de execução e falta de robustez para mapear reads de forma ambígua. Aproveitando a quantificação rápida de kallisto e a estimativa de variância, o sleuth supera esses problemas e fornece uma ferramenta estatisticamente rigorosa, flexível e eficiente para análise de RNA-Seq.
4. Pipeline básico (par-a-par)
Aqui vamos realizar uma comparação entre tratado (HS) e controle (CT). VocÊ deve ter quantificado as bibliotecas com o kallisto (VEJA AQUI) e o nome das pastas igual ao SRA ID. Assim estão descritos:
| condição | SRA ID |
|---|---|
| HS1 | SRR9696660 |
| HS2 | SRR9696664 |
| HS3 | SRR9696668 |
| CT1 | SRR9696658 |
| CT2 | SRR9696662 |
| CT3 | SRR9696666 |
Primeiro, carregue todas as bibliotecas R que serão usadas para as análises de hoje:
library(limma)
library(edgeR)
library(sleuth)
Criamos o dataset contendo todas as informações
base_dir <- "/opt/lucas/"
sample_id <- list('SRR9696658','SRR9696662','SRR9696666',
'SRR9696660','SRR9696664','SRR9696668')
paths <- list(paste(base_dir,"/sleuth/SRR9696658",sep=""),
paste(base_dir,"/sleuth/SRR9696662",sep=""),
paste(base_dir,"/sleuth/SRR9696666",sep=""),
paste(base_dir,"/sleuth/SRR9696660",sep=""),
paste(base_dir,"/sleuth/SRR9696664",sep=""),
paste(base_dir,"/sleuth/SRR9696668",sep=""))
names(paths) <- sample_id
Agora precisamos criar o arquivo de amostras, que criará dois modelos e os comparará: (1) Modelo base, que não depende do tratamento; (2) Modelo base que depende do tratamento. Assim, o sleuth verificará quais genes se alteram devido a inserção do tratamento. Abaixo está o arquivo amostras.txt que está dentro do base_dir, caso contrário altere a primeira linha do código abaixo:
| sample | condition | reps |
|---|---|---|
| SRR9696658 | Controle | 1 |
| SRR9696662 | Controle | 2 |
| SRR9696666 | Controle | 3 |
| SRR9696660 | Tratado | 1 |
| SRR9696664 | Tratado | 2 |
| SRR9696668 | Tratado | 3 |
s2c <- read.table(file.path(base_dir, "amostras.txt"), header = TRUE, stringsAsFactors=FALSE)
#s2c <- read.table(file.path(paste(base_dir,"/sleuth/",sep=""), "amostras.txt"), header = TRUE, stringsAsFactors=FALSE)
s2c <- dplyr::select(s2c, sample = sample, condition, reps)
s2c <- dplyr::mutate(s2c, path = paths)
s2c
s2c[] <- lapply(s2c, as.character)
#s2c <- data.frame(lapply(s2c, as.character), stringsAsFactors=FALSE)
Caso você tenha um arquivo transcrito-gene (aqui seria o arquivo t2g.txt), você pode passar para o Sleuth para que ele agregre por gene. O cabeçalho deve ser target_id e ens_gene. Vou descrever essas duas opçoes, mas rode um ou o outro.
t2g <- read.table("t2g.txt", header = TRUE, stringsAsFactors=FALSE)
======== Transcritos
so <- sleuth_prep(s2c, ~condition, extra_bootstrap_summary = TRUE)
so <- sleuth_fit(so)
so <- sleuth_wt(so, "conditionTratado")
======= Agregação em gene
so <- sleuth_prep(s2c, ~condition,aggregation_column="ens_gene", target_mapping = t2g,extra_bootstrap_summary = TRUE, gene_mode = TRUE)
so <- sleuth_fit(so)
so <- sleuth_wt(so, "conditionTratado")
Acesse o Shiny dos resultados (precisa de tela gráfica):
sleuth_live(so)
Analisaremos os resultados gerados
1. Diferenciais
results_table <- sleuth_results(so, test='conditionTratado', test_type = 'wald')
sleuth_significant <- dplyr::filter(results_table, qval <= 0.05)
head(sleuth_significant, 20)
sleuth_list <- sleuth_significant[,1]
write.table(sleuth_significant,file="diferenciais_sleuth.txt")
2. Volcano plot
pdf("Sleuth_Volcano.pdf")
plot(results_table$b, -1*log10(results_table$qval), col=ifelse(results_table$qval, "red", "black"),xlab="log(qval)", ylab="beta", title="Volcano plot", pch=20)
dev.off()
3. PCA
plot_pca(so, color_by = 'condition')
#ou
#caso queira deixar somente pontos e não labels: text_labels = FALSE
png("pca-transcritos.png",height = 8, width = 10,
units = 'in',res=600)
plot_pca(so,pc_x = 1L, pc_y = 2L, use_filtered = TRUE,
units = "est_counts", text_labels = TRUE, color_by = "condition",
point_size = 3)
dev.off()
4. Densidade
plot_group_density(so, use_filtered = TRUE, units = "est_counts",trans = "log", grouping = "condition", offset = 1)
5. Boxplots
Caso você queira saber coo estão dispostos os boxplots de um gene/transcrito específico, por exemplo, do geneX, basta fazer:
plot_bootstrap(so, "geneX", units = "est_counts", color_by = "condition")
plot_bootstrap(so, "geneX", units = "tpm", color_by = "condition")
6. Salvando resultados
Sempre é bom salvar o objeto 'so' do sleuth, assim podemos carregar ele a qualquer momento e usá-lo. Dica muito importante para aquele que utiliza o RStudio, por exemplo. Não esqueça de carregar a library sleuth após carregar o objeto salvo para utilizar as funções novamente.
# Salva o objeto 'so'
saveRDS(so, "so_sleuth.rds")
# Restaura ele em qualquer momento (nome qualquer)
so <- readRDS("so_sleuth.rds")
5. Análise temporal
O primeiro problema é definir o que queremos dizer com expresso diferencialmente no contexto de uma série temporal. Para um experimento do tipo tratamento versus controle, é simples: quais genes são mais ou menos expressos na condição de tratamento em comparação com o controle? A primeira generalização que vem à mente seria encontrar genes cujo nível aumenta ou diminui com os pontos de tempo. Isso significa encontrar correlação com pontos de tempo e se resume a regressão linear de expressão com pontos de tempo. No entanto, isso irá mascarar genes potencialmente interessantes que podem, por exemplo, atingir o pico no meio da série temporal. O que queremos é descobrir se a expressão de um gene segue um padrão geral em um grau mais elevado do que apenas o ruído.
Em uma estrutura de modelagem linear, como Sleuth, uma solução comum para isso é usar splines naturais. Dado um número de graus de liberdade para um modelo spline natural, nós serão colocados ao longo dos quantis das observações do eixo do tempo, que definirão os polinômios de base com suporte local.
A ideia agora é comparar o modelo de expressão de um gene que inclui os termos polinomiais que dependem do tempo, com um modelo que inclui apenas o termo ruído (base).
O script da análise temporal pode ser acessado aqui.
Citação
Pimentel, H., Bray, N., Puente, S. et al. Differential analysis of RNA-seq incorporating quantification uncertainty. Nat Methods 14, 687–690 (2017). https://doi.org/10.1038/nmeth.4324