# DNA de vírus: alinhamento e distância
#
# Script da prática: https://heronoh.github.io/bioinfo/aulas/r/Pratica_3b-Alinhamento_e_distancia-Virus.html
# Gerado automaticamente a partir da página da prática; para mudar o
# código, edite o .qmd correspondente.

# Prática ----

## Preparativos ----

### Carregando pacotes ----

#### CRAN - Instalando pacotes ----

# 00a - instalar pacotes CRAN----

#instale apenas caso você ainda não tenha o pacote ----
#    DICA: só é preciso instalar cada pacote uma única vez

# lista de pacotes CRAN que você quer garantir instalados
cran_pkgs <- c("pheatmap",
               "stringr",
               "BiocManager",
               "ggseqlogo",
               "tidyverse",
               "RColorBrewer")

# instala apenas os que não estão presentes
for (pkg in cran_pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
}

#### Bioconductor - Instalando pacotes ----

# 00b - instalar pacotes Bioconductor----

#instale apenas caso você ainda não tenha o pacote ----
#    DICA: só é preciso instalar cada pacote uma única vez

# lista de pacotes Bioconductor
bio_pkgs <- c("Biostrings", "DECIPHER")

# instala apenas os que faltam
for (pkg in bio_pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    BiocManager::install(pkg, ask = FALSE, update = FALSE)
  }
}

#### Carregando pacotes ----

# 01 - carregando pacotes ----
{
library("tidyverse")
library("Biostrings")
library("stringr")
library("dplyr")
library("DECIPHER")
library("tibble")
library("RColorBrewer")
library("pheatmap")
library("ggseqlogo")
}

### Definindo seu diretório de trabalho ----

# 02 - diretório de trabalho ----

# verifique onde stamos trabalhando e onde os arquivos estão sendo salvos
getwd()

#caso queira mudar de diretório, você pode usar 
#  setwd()
#utilize o tab para autocompletar com o cursor entre as aspas

### Obter arquivos ----

# 03 - baixar arquivo fasta com sequencias ----

download.file(
  "https://raw.githubusercontent.com/heronoh/bioinfo/main/aulas/r/arquivos/seqs_biomed_DNA.fasta",
  destfile = "seqs_biomed_DNA.fasta",
  mode = "wb",            # seguro no Windows/macOS/Linux
  method = "libcurl"      # segue redirecionamentos HTTPS direitinho
)

### carregar o .fasta como objeto no R ----

# 04 - ler arquivo fasta como objeto no R ----

fasta_seqs <- Biostrings::readDNAStringSet(filepath = "seqs_biomed_DNA.fasta")

fasta_seqs

### alinhamento de sequências ----

# 05 - alinhamentos ----

### a - orientar seqs ----
# orientar todas sequências na mesma direção (sobrescrevendo o objeto fasta_seqs)
#              (se alguma estiver no complemento reverso?)

fasta_seqs <- DECIPHER::OrientNucleotides(fasta_seqs)

DECIPHER::BrowseSeqs(fasta_seqs)

### b - alinhar seqs ----

fasta_seqs_alg <- DECIPHER::AlignSeqs(myXStringSet = fasta_seqs, 
                                      refinements = 100,
                                      iterations = 100,
                                      verbose = TRUE)

fasta_seqs_alg

### c - ver seqs ----

#veja seu alinhamento 
DECIPHER::BrowseSeqs(fasta_seqs_alg)

# onde começam e terminam as posições compartilhadas entre todos?

# 06 - Core alignment ----

# devemos utilizar apenas a parte comum a todos alinhamentos, jogar o resto fora!

#veja o objeto biostrings
fasta_subseqs_alg <-  Biostrings::subseq(fasta_seqs_alg,
                                         start = 1,   #substitua pela primeira posição compartilhada
                                         end = 2500)  #substitua pela ultima posição compartilhada

DECIPHER::BrowseSeqs(fasta_subseqs_alg)

### Matriz de distâncias ----

# 07 - Distâncias ----

### a - calcular a matriz distâncias ----

fasta_seqs_dist <- DECIPHER::DistanceMatrix(myXStringSet = fasta_subseqs_alg,
                                            includeTerminalGaps = FALSE,
                                            correction = "JC69",   # Jukes-Cantor (nas versões antigas do DECIPHER: "Jukes-Cantor")
                                            processors = 20)

# esta é a matriz de distâncias genéticas entre essas sequências
round(fasta_seqs_dist, 3)

View(fasta_seqs_dist)

### Visualizando a matriz de distâncias ----

# 08 - gráfico de distância ----

# vamos ver como fica a matriz

pheatmap::pheatmap(fasta_seqs_dist,
                   fontsize = 10,
                   color = RColorBrewer::brewer.pal(n = 9,name = "Greens"))

### Logo de sequências ----

# 09 - Gráfico de proporção ----

## plotar proporção de bases para uma subsequência (ggseq logo)

### a - selecionar apenas uma pequena região parcialmente conservada

fasta_sub_char <- Biostrings::subseq(fasta_subseqs_alg,886,926) %>% as.character()

seqs_logo_plot <- ggplot2::ggplot() +
  ggseqlogo::geom_logo(data = fasta_sub_char,
                       method = "probability") +
  ggplot2::annotate(geom = "rect",
                    xmin = 12.5,xmax = 33.5,
                    ymin = -0.1,ymax = 1.1,
                    fill = "#3cff00", alpha = 0.25)

seqs_logo_plot

ggplot2::ggsave(filename = "seqs_logo.png",
                device = "png",
                width = 30, height = 10,
                units = "cm")

# Referências ----
# https://www.bioconductor.org/help/course-materials/2016/BioC2016/ConcurrentWorkshops2/Wright/BigBioSeqData.pdf
# https://www.bioconductor.org/packages/release/bioc/vignettes/DECIPHER/inst/doc/DECIPHERing.pdf
# http://www2.decipher.codes/RLessons/RLesson14.html
# https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12687
# http://www2.decipher.codes/OligoDesign.html
# http://www2.decipher.codes/Bioinformatics/BigBioSeqData/BigBioSeqData2.html
# https://cran.r-project.org/web/packages/ggdendro/vignettes/ggdendro.html