DNA de vírus: alinhamento e distância

Conhecimentos básicos para lidar com sequências de DNA, RNA e proteínas utilizando ferramentas de R

Dica

Nota: Este material foi escrito para práticas em R usando o RStudio.

Conceitos fundamentais

Polímeros biológicos como sequências

DNA, RNA e proteínas são polímeros formados por monômeros ordenados:

• DNA e RNA: cadeias de nucleotídeos

• Proteínas: cadeias de aminoácidos

A ordem dos monômeros codifica informação biológica e determina função e estrutura: a sequência de DNA codifica RNAs e proteínas; a sequência de aminoácidos determina estrutura tridimensional das proteínas e sua atividade; a sequência de nucleotídeos de um RNA, codificador ou não de proteínas, pode determinar regiões de auto pareamentos, que resultam em estruturas secundárias e terciárias com capacidade catalítica, semelhantemente às proteínas.

O formato fasta

O fasta armazena sequências (de DNA, RNA ou proteínas).

Como os nucleotídeos e aminoácidos são moléculas definidas e fácilmente representadas por letras, a International Union of Pure and Applied Chemistry IUPAC padronizou o código de uma letra, que toda a área da Bioinformática de sequências utiliza em suas análises de genética, genômica, transcriptômica, proteômica e outras áras afins.

O arquivo fasta é muito simples, composto por múltiplas linhas que podem ser apenas de 2 tipos:

• Cabeçalho: linha iniciada por > com um identificador único da sequência e, opcionalmente, metadados.
• Sequência: linhas subsequentes com as letras do alfabeto correspondente (aminoácidos ou nucleotídeos).

Exemplo (fasta de DNA):

>seq1 organismo=Homo_sapiens gene=MT-ND1
ATGCCCTCGGCTTCTGGTAATACACCGCGGACATCGGCGTAAAGAGTGGTTAGGAAGTCTTTCAAACTACGGTTACGGAC
ATCGGCGTAAAGAGTGGTTAGGAAGTCTTTCAAACTATACGAATACGGACATCGGCGTAAAGAGTGGTTAGGAAGTCTTT
CAAACTACGGACATCGGCGTAAAGAGTGGTTAGGAAGTCTTTCAAACTAGGCCCAAGTTGATGAACATCGGCGTAAAAAG
TGGTTAGGAACAGAGCTAGATACTAAA

Repositórios de sequências

Muitos Bancos de Dados são dedicados para o armazenamento de sequências de DNA, RNA e proteínas. O principal destes é o NCBI nucleotide. Este banco guarda virtualmente todas sequências biológicas já registradas pela ciência, com livre acesso e possibilidade de download. Todo cientista que realiza algum estudo genético deve depositar os dados genéticos lá, para contribuir para o conhecimento coletivo.

Prática

Vamos conhecer um arquivo fasta? Siga as instruções a seguir no seu Rstudio. Copie, cole e execute os comandos abaixo

Preparativos

Carregando pacotes

Todos pacotes públicos criados para o R são de livre acesso. Qualquer usuário poide baixar e instalar com poucos comandos. Existem vários repositórios de pacotes (packages) e na bioinforática utilizamos práticamente dois:

• CRAN (pacotes “gerais” do R)
  Repositório amplo e genérico (estatística, web, gráficos, etc.). Publicação via install.packages(). Ciclo de versões independente entre pacotes; revisão é mais básica e não força integração entre eles.
• Bioconductor (instalado/gerenciado via BiocManager)
  Ecossistema focado em bioinformática/ômicas (RNA-seq, ChIP-seq, variantes, anotação, análise de dados de sequenciamento, …). Apresenta larga padronização em formatos de resultados. Possui annotation/data packages (ex.: genomas, bancos de genes) e serviços como AnnotationHub/ExperimentHub. Tem releases semestrais sincronizadas com a versão do R e checagem de compatibilidade entre pacotes, visando reprodutibilidade.

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")
}

::: Dica: caso receba alguma mensagem de que o pacote não foi encontrado, você pode instalá-lo com:

install.packages("pacote")

::: Dica: Alguns pacotes de bioinformática são fornecidos pelo Bioconductor, e podem ser instalados com o gerenciador de pacotes BiocManager

install.packages("BiocManager")
BiocManager::install("pacote")

Definindo seu diretório de trabalho

Para cada projeto, sempre crie uma pasta onde você vai salvar todos arquivos e scripts, assim você mantém seu trabalho sempre organizado. No comando abaixo, coloque dentro das aspas o caminho para o seu diretório de trabalho.

::: Sempre utilize a tecla TAB para autocompletar os caminhos a partir do ~ ou / (ou ****, caso no windows). Para isso, coloque o cursor (|, a barra onde sua digitação aparece) entre as ” “.

# 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

Baixe o arquivo .fasta que utilizaremos.

# 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

Agora você vai carregar para o Ambiente (Environment) do R o arquivo que você baixou, o salvando em um objeto. Vamos fazer isso com o pacote Biostrings. O pacote Biostrings é um dos pilares do ecossistema Bioconductor, amplamente usado em bioinformática, genômica e análise de sequências biológicas. Ele fornece estruturas de dados e funções eficientes para representar, manipular e analisar sequências de DNA, RNA e proteínas em R.

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

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

fasta_seqs

alinhamento de sequências

Vamos explorar esse objeto. As sequências parecem muito diferentes, mas na verdade só estão desorientadas. Como o DNA é dupla fita, você pode usar a informação em qualquer uma das fitas (sempre no sentido 5’ -> 3’) mas, para analisar vários DNAs juntos, é necessário colocar todos na mesma orientação.

O pacote DECIPHER oferece ferramentas para realizar a orientação, alinhamentos e visulização 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 = "JC",
                                            processors = 20)

# esta é a matriz de distâncias genéticas entre essas sequências
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"))


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

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

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

fasta_sub_char <- Biostrings::subseq(fasta_subseqs_alg,1,99) %>% 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

E por hoje ficamos por aqui!

---

Voltar para a página inicial