analiseBioinfo.Rmd

---
title: "Notebook do projeto do morcegão"
author:
- Emerson
- Marcel
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
---

## Pacotes

```{r pacotes}
library(tidyverse)
library(Biobase)
library(limma)
library(AnnotationDbi)
library(clusterProfiler)
library(org.Rn.eg.db)
library(multiMiR)
```

## Objetivo

Encontrar miRNAs para validar com dados da Georgia.

## Meta


$$nrow(expMat) == nrow(adf)\\ncol(expMat) == nrow(sdrf)$$
- `expMat` é uma matriz!  
- `adf` é um data.frame anotado!  
- `sdrf` é um data.frame anotado!

## Importar dados

### Matriz de expressão (expMat)
```{r}
sample_files <- list.files(path = "data/",pattern = "sample",full.names = T)

sample_list <- list() #criar uma lista vazia

for (i in 1:length(sample_files)) { #importar dados
  
  sample_list[[i]] <- read_tsv(file = sample_files[i])
    
  print(length(sample_list))
}

```

```{r}
#Nomear os elementos

sample_files#ver nomes

#remove os caracteres indesejados
tmp1 <- str_replace(string = sample_files,pattern = "data/",replacement = "")
tmp2 <- str_replace(string = tmp1,pattern = "_sample_table.txt",replacement = "")
#gsub()

#Passar o tmp2 para os nomes da lista

names(sample_list) <- tmp2

```

### Descrição de amostras

```{r importante}
sdrf <- read_tsv(file = "data/E-GEOD-68913.sdrf.txt")

sdrf <- sdrf %>% 
  mutate(NomesDados = str_replace(string = `Source Name`,pattern = " 1",replacement = "")) %>% 
  dplyr::select(NomesDados,Source = "Comment [Sample_source_name]",Disease = "Characteristics [disease status]") %>% 
  arrange(NomesDados)
```

### Descrição de amostras

```{r}
adf <- read_tsv(file = "data/A-GEOD-16384.adf.txt",skip = 15)


adf <- adf %>% 
  filter(str_detect(string = `Reporter Database Entry [mirbase]`,pattern = "rno")) %>% 
  select(Rep = "Reporter Name", mirbase = `Reporter Database Entry [mirbase]`)
```


## Arrumar os dados

### Criar um unica tabela



```{r}
expTab <- sample_list %>% 
  purrr::reduce(left_join,by = "ID_REF") # junta todas as tabela em uma só baseada na colunas ID_REF


colnames(expTab)[2:13] <- names(sample_list) # Arruma os nomes de colunas

class(expTab) #retorna a classe do objeto

expMat <- expTab %>% 
  dplyr::filter(ID_REF %in% adf$Rep) %>% 
  column_to_rownames(var = "ID_REF") %>% 
  as.matrix() # Converte para matriz

ncol(expMat) == nrow(sdrf)
nrow(adf) == nrow(expMat)
(colnames(expMat) == sdrf$NomesDados) 
```

## Criar eset

```{r}
rownames(adf) <- adf$Rep
row.names(sdrf) <- sdrf$NomesDados


adf_final <- Biobase::AnnotatedDataFrame(adf)
sdrf_final <- Biobase::AnnotatedDataFrame(sdrf)

#Criar o conjunto de expressão
eset <- ExpressionSet(assayData = expMat,
                      phenoData = sdrf_final,
                      featureData = adf_final)

#eset é um objeto do S4 - programação orientada a objetos
# O R padrão é S3 - programação funcional
```

## LIMMA  

Modelo linear para os fatores que influenciam na expressão gênica/transcritos/mirna/etc.

$$f(x) = Par_1 + Par_2 + ...$$


$$Exp_{gene} = Source + Disease$$

$$Exp_{gene} = \beta0 + \beta1 Source + \beta2 Disease + \epsilon$$


### Densidades

```{r}
plotDensities(eset,legend = F)
```

aplicar log10

```{r}
# exprs(eset) #extrai a matriz de expressão
# pData(eset) #extrai os dados das amostras
# fData(eset) #extrai os parametros dos mirna

eset_log <- eset

min(exprs(eset_log))

#USAR QUANDO TIVER NUMERO NEGATIVO
exprs(eset_log) <- log10(exprs(eset) + abs(min(exprs(eset_log)))) #Abs = módulo

#USAR log(x+1) quando tiver inteiros 0,1,2,3,..


#replotar as densidades
plotDensities(eset_log,legend = F)
```

Escala normal
$$x = 0$$
Transformar para escala log
$$y = log(x+1)\\y= log(1)\\y=0$$

Difença entre escalas
-1000 e 100
$$1000-100 = 900$$
$$log_{10}(1000) - log_{10}(100)\\
3-2\\
1$$

### Normalizar

```{r}
eset_log_norm <- eset_log

exprs(eset_log_norm) <- normalizeBetweenArrays(exprs(eset_log),
                                               method = "quantile")

exprs(eset) <- normalizeBetweenArrays(exprs(eset),
                                      method = "quantile")

plotDensities(eset,legend = F)
```

### Expressão diferencial

```{r}
pData(eset)

(design <- model.matrix(~ 0 + factor(pData(eset)$Source) : factor(pData(eset)$Disease)))

(colnames(design) <- c("mv_CKD","vsm_CKD",
                       "mv_normal","vsm_normal"))

design
# model.matrix(~0 + Source * Disease, data = pData(eset))
fit <- lmFit(exprs(eset),
            design)

cm <- makeContrasts(b1 = "mv_CKD + vsm_CKD - mv_normal - vsm_normal", #Doença
                    b2 = "mv_CKD + mv_normal - vsm_CKD - vsm_normal", #Fonte
                    b3 = "mv_normal - vsm_normal", #Fonte - comparando normal
                    b4 = "vsm_normal - vsm_CKD", #Doença - comparando nas VSM
                    levels = design)


fit2 <- contrasts.fit(fit, cm)
fit2 <- eBayes(fit2)

#10 primeiras linhas dos coeficientes
fit2$coefficients[1:10,]

#Para o coeficiente 1, ou seja, efeito da doença na expressão de miRNA
topTable(fit2,
         coef = 1,#o primeiro contraste avaliado
         adjust.method = "BH",#Benjamin hocheback - Falsos positivos corrigidos
         number = 10 #numero de linhas para serem exibidas
         )
```

### Volcano plot

$$|log_2{FC}| > 1 \\|FC|>2$$

$$-log_{10}{adjPval} > 2\\
adjPval < 10^{-2}\\adjPval <0.01$$

**Organizar o res (renomear) e rodar para o coef 2**

```{r}
res <- topTable(fit2,
         coef = 1,#o primeiro contraste avaliado
         adjust.method = "BH",#Benjamin hocheback - Falsos positivos corrigidos
         number = Inf #numero de linhas para serem exibidas
         )

#Volcano plot - log2FC
res %>% 
  rownames_to_column(var = "miRNA") %>% 
  mutate(Cond = case_when(logFC > 1 & adj.P.Val < 0.01 ~ "Up",
                          logFC < -1 & adj.P.Val < 0.01 ~ "Down",
                          TRUE ~ "Not Sig")) %>%  #cria uma coluna
  ggplot(aes(x = logFC, y = -log10(adj.P.Val),color = Cond,label = miRNA)) +
  geom_point() +
  # geom_vline(xintercept = c(-1,1)) +
  # geom_hline(yintercept = 2) +
  labs(title = "Expressao diferencial - Volcano plot B1") +
  scale_color_manual(values = c("blue","black","red"))

#plotly::ggplotly()
```

Extrair os miRNA diferencialmente regulados

```{r}
res %>% 
  rownames_to_column(var = "miRNA") %>% 
  mutate(Cond = case_when(logFC > 1 & adj.P.Val < 0.01 ~ "Up",
                          logFC < -1 & adj.P.Val < 0.01 ~ "Down",
                          TRUE ~ "Not Sig")) %>% 
  filter(Cond != "Not Sig") %>% 
  select(miRNA, logFC,adj.P.Val,Cond)#Opcional
```

## Enriquecimento funcional(EF)

Para analise de EF de miRNAs usamos seus alvos.

**Gerar os targets dos miRNAs com multiMir e então rodar o enriquecimento**

```{r}
mirs <- res %>% 
  rownames_to_column(var = "miRNA") %>% 
  mutate(Cond = case_when(logFC > 1 & adj.P.Val < 0.01 ~ "Up",
                          logFC < -1 & adj.P.Val < 0.01 ~ "Down",
                          TRUE ~ "Not Sig")) %>% 
  filter(Cond != "Not Sig") %>% 
  pull(miRNA)

str_replace(string = mirs,
            pattern = "-star_st|_st",
            replacement = "")

res_target <- get_multimir(org = "rno",
             mirna = str_replace(string = mirs,
            pattern = "-star_st|_st",
            replacement = ""),
             table = "predicted",
             limit = 50)
res_target@data

library(miRNAtap)
library(miRNAtap.db)
g2Predicts <- function(x, db = c('pictar','diana','targetscan','miranda','mirdb'), spc = "rno" ){
  tab <- NULL
  miRNA_targets <- NULL
  for(i in 1:length(x)){
    target <- getPredictedTargets(x[i], sources = db, both_strands = TRUE,species = spc)
    if(is.null(target[i])){
      dat <- data.frame(miRNA = x[i], targetID = NA, Database = "Targetscan, Pictar")
    } else dat<- data.frame(miRNA = x[i], targetID = rownames(target), Database = "Targetscan, Pictar")
    tab <- rbind(tab, dat)
  }
  return(tab)
}


tmpmirs <- g2Predicts(x = str_replace(string = mirs,
            pattern = "-star_st|_st",
            replacement = ""))

tmpmirs %>% 
  group_by(miRNA) %>% 
  count()
```

#### Trocar anotações

```{r}
rno <- org.Rn.eg.db
targets_id <- AnnotationDbi::select(x= rno, 
                     keys = unique(tmpmirs$targetID),
                     columns = c("SYMBOL","ENTREZID"),
                     keytype = "ENTREZID")

targets_id$SYMBOL

```

#### KEGG pathways 

```{r}
targets_id_NA <- targets_id %>% 
  filter(!is.na(ENTREZID)) %>% 
  distinct()

kegg <- enrichKEGG(gene = (targets_id_NA$ENTREZID),
                   organism     = "rno",
                   #nPerm        = 1000,
                   minGSSize    = 3,
                   maxGSSize    = 800,
                   pvalueCutoff = 0.05,
                   pAdjustMethod = "BH", ##no p-value correction or BH
                   keyType       = "ncbi-geneid")
```