RNAseq_DESeq2_tutorial_curso_omicas.rmd

---
title: "TRANSCRIPTÔMICA: ANÁLISE DE EXPRESSÃO GÊNICA DIFERENCIAL"
author: (Computational Systems Biology Laboratory)
output: 
  html_document:
    theme: simplex
    css: style.css
    number_sections: true
    df_print: paged
---

Tutorial de análise de expressão gênica diferencial utilizando o pacote **DESeq2** no R/RStudio preparado para a disciplina **ICB5747 - Ciências Ômicas em Doenças Infecciosas**.

# Chamada de pacotes e configurações básicas do R

------------------------------------------------------------------------

Para começar, precisamos carregar os pacotes necessários para a nossa análise, que além do próprio DESeq2, incluem pacotes para manipulação de dados e acesso à banco de dados com informações biológicas.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(webshot)
#webshot::install_phantomjs() # Ana, precisamos passar essa linha de código pro Carlos incluir no tutorial prévio, pois é necessário instalar esse pacote "webshot" e rodar essa linha de código aí pra visualizar páginas web no markdown
```

```{r, message=FALSE, warning=FALSE}
library(DESeq2)
library(data.table)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(biomaRt)
options(stringsAsFactors = F)
```

# Um pouco sobre os dados que vamos analisar

------------------------------------------------------------------------

O número de acesso para o estudo que usaremos como exemplo é o GSE52778 (Gene Expression Omnibus database - GEO), *Human Airway Smooth Muscle Transcriptome Changes in Response to Asthma Medications*. O experimento inclui 4 linhagens celulares diferentes que foram tratadas com o medicamento dexametasona durante 18 horas. Esse medicamento é usado por pacientes com asma para reduzir a inflamação nas vias aéreas. Para cada linhagem celular, o experimento contém uma amostra tratada e uma amostra não tratada (controle). O objetivo é entender as mudanças transcricionais que ocorrem durante o tratamento com dexametasona.

# Carregando as tabelas necessárias

------------------------------------------------------------------------

Agora vamos carregar as tabelas de contagem de genes e a tabela que contém a informação sobre as amostras (phenodata). Vamos converter também a coluna que indica o tratamento por dexametasona em fatores para indicar a ordem correta dos grupos (controle primeiro, tratado depois).

```{r}
count_data <- read.csv(file = "https://github.com/csbl-inovausp/RNAseq_DESeq2_R_tutorial/raw/main/data/counts_data.csv", header = T)
head(count_data)
```

```{r}
pheno_data <- read.csv(file = "https://github.com/csbl-inovausp/RNAseq_DESeq2_R_tutorial/raw/main/data/sample_info.csv", header = T)
pheno_data$dexametasona <- factor(pheno_data$dexametasona,levels = c("controle","tratado"))
head(pheno_data)
```

# Criando o objeto dds para DESeq2

------------------------------------------------------------------------

O próximo passo será criar o objeto dds, específico para as análises utilizando o DESeq2. Esse objeto irá conter a tabela de expressão, o phenodata e a indicação do design do experimento (tratado versus controle, por exemplo).

```{r, warning=FALSE}
dds <- DESeq2::DESeqDataSetFromMatrix(countData = count_data,
                                      colData = pheno_data,
                                      tidy = F,
                                      design = ~linhagem+dexametasona)
dds
```

## Filtrando os genes pouco expressos

A nossa tabela de contagem contem muitas linhas (genes) em que não há contagem de reads em nenhuma amostra e também genes em que a média de contagens nas amostras é muito baixa. É desejável filtrar esses genes com contagem zero ou muito baixa para aumentar a velocidade do processamento e também para evitar resultados que não apresentam validade. Primeiro vamos filtrar as linhas com contagem zero em todas as colunas.

```{r}
nrow(dds)
```

```{r}
keep <- rowSums(DESeq2::counts(dds)) > 1
dds_filtrado <- dds[keep,]
nrow(dds_filtrado)
```

Em seguida, filtramos as linhas em que pelo menos três amostras tenham contagem de 10 ou superior.

```{r}
keep_10 <- rowSums(DESeq2::counts(dds_filtrado) >= 10) >= 3
dds_filtrado_final <- dds_filtrado[keep_10,]
nrow(dds_filtrado_final)
```

## Visualizando as amostras usando Análise de Componentes Principais (PCA)

A próxima etapa é a visualizar as nossas amostras utilizando a Análise de Componentes Principais (PCA). Essa análise irá nos mostrar se as amostras controles e amostras tratadas são diferentes entre os grupos e similares dentro do mesmo grupo, levando em conta a contagem de todos os genes. Essa é uma análise exploratória inicial muito importante para verificar se os grupos experimentais foram bem separados.

```{r, fig.align = 'center'}
vsd <- DESeq2::vst(dds_filtrado_final, blind = FALSE)
DESeq2::plotPCA(vsd, intgroup = "dexametasona")
```

# Executando a análise de expressão diferencial (DESeq2)

------------------------------------------------------------------------

Vimos que as amostras controle e tratadas estão bem separadas umas das outras e mais próximas entre si. Podemos seguir adiante para a análise de expressão gênica diferencial entre os grupos. A função usada para essa análise é bem simples:

```{r, message=FALSE}
dds_de <- DESeq2::DESeq(dds_filtrado_final)
dds_de
```

## Obtendo os resultados de análise de expressão diferencial

Pronto, a análise já foi feita. Agora precisamos criar a nossa tabela de genes diferencialmente expressos contendo os valores de log2FoldChange e valores de p (p-value) e valores de p ajustados para múltiplas comparações (p-adjusted).

```{r}
res <- DESeq2::results(dds_de)
res[1:50,]
```

## Redefinindo os grupo experimental e controle (caso necessário)

Caso você não tenha definido previamente o grupo controle e o grupo tratado (ou tenha por engano errado a ordem), nada está perdido! Ainda é possível definir manualmente quais são os grupos que você deseja comparar. Isso também pode ser interessante quando temos mais de um grupo tratado ou diferentes tipos de controle, os quais desejamos comparar entre si de forma combinada. Para isso, basta adicionar o argumento `contrast` na função `result` e criar um vetor contendo o nome da coluna do phenodata que tem a informação dos grupos, o nome do grupo tratado (primeiro) e o nome do grupo controle (depois).

```{r}
res <- DESeq2::results(dds_de, contrast = c("dexametasona","tratado","controle"))
res[1:50,]
```

## Contando o número de genes diferencialmente expressos

Vamos agora verificar o número de genes que tiveram expressão diferencial significativa de acordo com o teste estatístico. O DESeq2 tem por padrão o valor de corte de p-adjusted \< 0.1 como significativo.

```{r}
DESeq2::summary(res)
```

## Ajustando o rigor estatístico

O número de genes diferencialmente expresso com p-adjusted \< 0.1 foi bastante alto. Vamos aumentar a restrição e considerar significativos apenas aqueles genes com p-adjusted \< 0.05, que é o padrão visto em muitos artigos na literatura científica.

```{r}
res.05 <- results(dds_de, alpha = 0.05)
DESeq2::summary(res.05)
```

# Visualizando os resultados

------------------------------------------------------------------------

Vamos utilizar o pacote ggplot2 para criar o volcano plot (gráfico de vulcão). Nesse gráfico, o eixo x representa o valor de log2FoldChange, que como visto durante o curso, representa o tamanho do efeito da expressão diferencial. Valores positivos de log2FoldChange representam genes cuja expressão está aumentada no grupo tratado em relação ao grupo controle, enquanto valores de log2FoldChange negativo representam o contrário, genes cuja expressão no grupo tratado é menor que a expressão no grupo controle. Além do filtro de p-adjusted \< 0.05, vamos também criar um filtro de log2FoldChange. Apenas genes com p-adjusted \< 0.05 e log2FoldChange maior que 1 ou menor que -1. Esses valores de log2FoldChange representam uma diferença entre tratado e controle que significa que o tratado tem o dobro da expressão do controle (log2FC = 1) ou que o tratado tem a metade da expressão do control (log2FC = -1).

## Criando a tabela para plotar o volcano plot

Primeiro, vamos criar a tabela de resultados que vamos usar para criar o volcano plot.

```{r}
res_df <- as.data.frame(res)
res_df[1:50,]
```

## Plotando o volcano plot

Agora, vamos criar uma coluna no dataframe indicando se cada gene é significativo (padj \< 0.05 & \|logFC\| \> 1) - positivo ou negativo. Já vamos também definir a cor de cada gene pro volcano plot: azul, significativo e negativo, vermelho, significativo e positivo e cinza, não significativo, independente da direção.

```{r}
plot_df <- res_df %>%
  mutate(significativo=ifelse(padj<0.05 & abs(log2FoldChange) > 1,
                              yes = "sim",
                              no = "não"),
         dir=ifelse(log2FoldChange>0,yes = "positivo",no = "negativo"),
         DE=ifelse(significativo=="sim" & dir=="positivo",
                   yes = "positivo",
                   no = ifelse(significativo=="sim" & dir=="negativo",
                               yes = "negativo",
                               no = "não significativo")),
         DE=factor(DE, levels = c("negativo","não significativo","positivo")))
plot_df[1:50,]
```

Então, plotamos o gráfico:

```{r, fig.align='center'}
plot <- plot_df %>%
  ggplot(aes(x=log2FoldChange,y = -log(padj))) +
  geom_point(aes(color=DE)) +
  scale_color_manual(values = c("blue","grey","red3")) +
  geom_hline(yintercept = -log(0.05),linetype="dashed") +
  geom_vline(xintercept = c(-1,1),linetype="dashed") +
  theme_bw()
plot
```

## Salvando a figura em PDF

Podemos agora salvar a figura em um arquivo png ou pdf:

```{r}
# pdf
ggsave(path = "figures",filename = "volcano_padj0.05_logFC1.pdf",
       plot = plot,device = "pdf",dpi = 300)
# png
ggsave(path = "figures",filename = "volcano_padj0.05_logFC1.png",
       plot = plot,device = "png",dpi = 300)
```

## Convertendo os nomes dos genes

Repare que o nome dos genes está no formato "ENSG..." que é o formato Ensemble gene id. Podemos buscar os nomes desses genes no Google como estão ou converter os nomes dos genes em official gene symbol, que é são os nomes de genes que estamos acostumados a ler, como IL6 e GAPDH. Para isso, vamos usar o pacote biomaRt.

```{r}
ensembl_genes <- rownames(res_df)
mart <- biomaRt::useMart("ensembl",dataset="hsapiens_gene_ensembl")

converted_genes <- biomaRt::getBM(attributes = c("ensembl_gene_id","external_gene_name"),
                                  filters = "ensembl_gene_id",
                                  values = ensembl_genes,
                                  mart = mart)

converted_genes[1:100,]
```

## Criando a tabela de DEGs com nomes de genes oficiais

Agora com a tabela contendo a relação direta entre os ids ensembl e os nomes oficiais dos genes, vamos criar uma nova tabela de reusltados substituindo os ids ensembl pelos nomes. Repare ao final que muitos genes possuem id ensembl, mas não possuem nome de gene oficial (exemplo: ENSG00000269772). Isso acontece, pois novos genes são primeiro catalogados pelo Ensembl, que não é responsável pela criação do nome oficial do gene. Portanto, existem genes que são descobertos e não nomeados, mas ainda assim recebem ids no catálogo Ensembl.

```{r}
res_convert <- res_df %>%
        tibble::rownames_to_column("ensembl_gene_id") %>%
        dplyr::left_join(converted_genes,by = "ensembl_gene_id") %>%
        dplyr::select(ensembl_gene_id,external_gene_name,everything())
res_convert[1:10,]
```

# Interpretando os resultados

------------------------------------------------------------------------

## Analisando os top 20 genes up e down

Com os nomes de genes convertidos, vamos agora dar uma olhada nos top genes diferencialmente expressos para começarmos a interpretar os resultados. Para isso vamos criar uma tabela com os top 20 genes mais upregulados e outra com os top 20 genes mais downregulados.

```{r}
top_up <- res_convert %>%
  dplyr::filter(external_gene_name != "",
                !is.na(external_gene_name)) %>%
  dplyr::top_n(n = 20,wt = log2FoldChange)
top_up
```

```{r}
top_down <- res_convert %>%
  dplyr::filter(external_gene_name != "",
                !is.na(external_gene_name)) %>%
  dplyr::top_n(n = 20,wt = -log2FoldChange)
top_down
```

Vamos também salvar tabelas fora do R com os top DEGs up e down.

```{r}
write.csv(x = top_up,file = "data/top_up.csv",row.names = F,quote = F)

write.csv(x = top_down,file = "data/top_down.csv",row.names = F,quote = F)
```


## Buscando informações sobre top DEGs no Google

O primeiro passo para começar a interpretar os resultados de expressão diferencial é buscar o nomes dos genes mais diferencialmente expressos no Google e começar a ler e estudar sobre esses genes. Vamos ainda realizar uma análise de enriquecimento funcional que irá nos dizer quais vias ou processos bioquímicos estão mais representados entre os genes diferencialmente expressos, mas esse primeiro passo de fazer uma busca inicial é importante para começar a entender os resultados.

Um dos sites que costuma aparecer primeiro na busca é o [GeneCards](https://www.genecards.org/). Esse site, assim como a Wikipedia, é muito confiável para começarmos a aprender mais sobre os genes que temos em mãos. Veja abaixo o site GeneCards para o gene FKBP5, que foi o que teve o maior valor positivo de log2FC na nossa análise.

## Top gene up (FKBP5)

```{r, echo=FALSE}
knitr::include_url("https://www.genecards.org/cgi-bin/carddisp.pl?gene=FKBP5")
```

Veja que a entrada para esse gene diz que: "The protein encoded by this gene is a member of the immunophilin protein family, which play a role in immunoregulation and basic cellular processes involving protein folding and trafficking". Lembre-se que estamos analisando dados de células da via aérea tratadas com um medicamento antiinflamatório. É interessante notar que o gene cuja expressão mais aumentou com esse tratamento foi um gene envolvido com imunoregulação. Vejamos agora o que o GeneCards tem a dizer sobre o gene cuja expressão foi mais reduzida nas células tratadas.

## Top gene down (CYP24A1)

```{r, echo=FALSE}
knitr::include_url("https://www.genecards.org/cgi-bin/carddisp.pl?gene=CYP24A1")
```

Veja que esse é um gene que codifica uma enzima envolvida no metabolismo de drogas, mas também na síntese de esteróides. A droga utilizada no experimento que estamos analisando, a dexametasona, é um esteróide sintético. Portanto, faz sentido que células tratadas com uma substância esteróide passem a expressar menor quantidade dos genes que estão envolvidos na produção dessas moléculas. Provavelmente trata-se de um mecanismo de feedback negativo devido ao tratamento. Podemos investigar essa hipótese mais a fundo olhando para outros genes relacionados à via de síntese de esteróides nos nosso resultados.

## Enriquecimento funcional utilizando o EnrichR

Agora que já obtivemos a lista de DEGs completa, podemos fazer análise de enriquecimento funcional usando a plataforma [enrichR](https://maayanlab.cloud/Enrichr/). Vamos acessar a tabela de top DEGs up e down, utilizando o Excel ou outro editor de planilhas e seguir os mesmos passos que fizemos ontem na aula teórica.