index.qmd

---
title: "MaRe - Resultados preliminares"
author:
  - name: João V. F. Cavalcante
    orcid: 0000-0001-7513-7376
  - name: Rodrigo Dalmolin
    orcid: 0000-0002-1688-6155
format:
  html:
    embed-resources: true
    toc: true
bibliography: references.bib
lang: pt
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE,
  message = FALSE,
  warning = FALSE,
  results = TRUE
)

library(tidyverse)
library(reactable)
```

# Descrição das amostras


```{r}
metadata <- read.csv("data/metadata.csv")

reactable(
  metadata |> distinct() |> arrange(location) |> setNames(c("Amostra", "Local")),
  bordered = TRUE,
  highlight = TRUE,
  striped = TRUE,
  compact = TRUE
)
```

# Metodologia

A análise foi realizada com o pipeline nf-core/mag v4.0.0 [@nf-core-mag], orquestrado pelo gerenciador de fluxos de trabalho Nextflow e utilizando o perfil Singularity. O script utilizado para a análise pode ser encontrado [aqui](https://github.com/dalmolingroup/mare/blob/main/scripts/run_mag.sh).

# Controle de Qualidade

## Contagens de reads

Cada amostra possui por volta de 60M a 100M de *reads* não-duplicadas.

Abaixo, valores obtidos com o software FastQC [@FastQC].

![Contagens do número absoluto de *reads* duplicadas e não duplicadas em cada amostra](./results/figures/fastqc_sequence_counts_plot.png){fig-align="center"}

![Porcentagem de *reads* duplicadas e não duplicadas em cada amostra](./results/figures/fastqc_sequence_counts_plot_percentage.png){fig-align="center"}

## Filtragem de *reads*

Utilizando fastp [@Chen2018], por volta de 99% das *reads* foram mantidas em todas as amostras - uma retenção alta e adequada pra fazer uma caracterização profunda do microbioma.

![Porcetagem de *reads* retidas em cada amostra](./results/figures/fastp_filtered_reads_plot.png){fig-align="center"}

## Classificação taxonômica

Utilizando o banco de dados *standardpluspf* (Bacteria, Archaea, Protozoa e Fungi) do Kraken2 [@Wood2019] para observar a classificação taxonômica das leituras, observamos que por volta de 50-60% das *reads* resultaram em nenhuma classificação, o que é esperado dado a profundidade do sequenciamento.

![Resultado da classificação taxonômica com Kraken2, a nível de reino](./results/figures/Kraken2-top-n-plot_kingdom.png){fig-align="center"}

![Resultado da classificação taxonômica com Kraken2, a nível de filo](./results/figures/Kraken2-top-n-plot.png){fig-align="center"}

![Resultado da classificação taxonômica com Kraken2, a nível de família](./results/figures/Kraken2-top-n-plot_familia.png){fig-align="center"}

# Montagem

Após montagem com o MEGAHIT [@Li2015], observamos que a maioria dos *contigs* possuem um comprimento de até 1000 pares de base.

![Porcentagem de *contigs* de cada comprimento nas montagens de cada amostra](results/figures/quast_num_contigs.png){fig-align="center"}

Valores de N50 e N75 mostram resultados bons.

![Distribuições de estatísticas relativas às montagens de cada amostra.](results/figures/quast_table.png){fig-align="center"}

# Binning

Após uma abordagem híbrida de binning utilizando os softwares MaxBin [@Wu2016], CONCOCT [@Alneberg2014] e MetaBAT2 [@Kang2019], seguida de refinamento com o DASTool [@Sieber2018], pudemos obter inicialmente **1,549 *bins*, que foram refinados para 137 *metagenome-assembled genomes* (MAGs)**.

Após o refinamento, notamos melhorias quanto ao comprimento dos *contigs* quando comparados às montagens brutas.

![Porcentagem de *contigs* de cada comprimento nas MAGs obtidas](results/figures/quast_num_contigs-1.png){fig-align="center"}

## Predição de domínio

Realizamos um passo de predição de sequências eucarióticas, procarióticas e de archaea com o software Tiara [@Karlicki2022], mostrando que quase todas as MAGs obtidas são procarióticas.

```{r}
tiara_res <- read_table("results/nfcore_mag_results/GenomeBinning/Tiara/tiara_summary.tsv") |> 
  select(BinID, classification)
tiara_res$BinID <- gsub("_", ".", tiara_res$BinID)
bin_summary <- read_tsv("results/nfcore_mag_results/GenomeBinning/bin_summary.tsv") |> 
  select(bin, Name)
bin_summary$unrefined_name <- gsub("Refined|\\.fa|_sub", "", bin_summary$bin)

count_df <- tiara_res %>%
  count(classification, name = "number_of_bins")

ggplot(count_df,
       aes(
         x = reorder(classification, -number_of_bins),
         y = number_of_bins,
         fill = classification
       )) +
  geom_col(show.legend = FALSE, width = 0.6) +
  geom_text(
    aes(label = number_of_bins),
    vjust = -0.5,
    size = 4,
    color = "black"
  ) + # Adds count labels on top of bars
  scale_fill_brewer(palette = "Set2") + # A nice color palette
  labs(title = "Número de Bins não-refinados por reino (predito)", x = "Classificação", y = "# bins") +
  theme_minimal(base_size = 14) + # A clean theme for the plot
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid.major.x = element_blank(),
    # Removes vertical grid lines
    panel.grid.minor.y = element_blank()
  )
```

Considerando apenas bins pós-refinamento, a predominância de *prokarya* permanece.

```{r}
classifs <- inner_join(tiara_res, bin_summary, by = c("BinID" = "unrefined_name"))

count_df <- classifs %>%
  count(classification, name = "number_of_bins")

ggplot(count_df,
       aes(
         x = reorder(classification, -number_of_bins),
         y = number_of_bins,
         fill = classification
       )) +
  geom_col(show.legend = FALSE, width = 0.6) +
  geom_text(
    aes(label = number_of_bins),
    vjust = -0.5,
    size = 4,
    color = "black"
  ) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Número de Bins refinados por reino (predito)", x = "Classificação", y = "# bins") +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid.major.x = element_blank(),
    # Removes vertical grid lines
    panel.grid.minor.y = element_blank()
  )
```

## Classificação taxonômica dos bins

Realizamos a classificação taxonômica dos bins para **táxons de Bacteria e Archaea** de acordo com a referência GTDB-Tk [@Chaumeil2022].

```{r}
gtdb_summary <- read_tsv("results/nfcore_mag_results/Taxonomy/GTDB-Tk/gtdbtk_summary.tsv")
```

```{r}
#| fig-cap: "Distribuição dos principais filos identificados a partir de genomas montados de metagenomas (MAGs). O gráfico de barras exibe o número total de MAGs classificados em diferentes filos com base no Genome Taxonomy Database (GTDB). A análise revela que *Pseudomonadota* e *Bacteroidota* são os filos mais representados no conjunto de dados. Mais de 40 MAGs não puderam ser classificadas, 8 pertencendo a Eukarya de acordo com o Tiara."

phylum_counts <- gtdb_summary %>%
  # Split the classification string by ';' and get the second element (phylum)
  separate(classification, into = c("domain", "phylum"), sep = ";", extra = "drop") %>%
  # Clean up the phylum name by removing the 'p__' prefix
  mutate(phylum = str_replace(phylum, "p__", "")) %>%
  # Count the occurrences of each phylum
  count(phylum, sort = TRUE)

# Create the plot
ggplot(phylum_counts, aes(x = n, y = reorder(phylum, n))) +
  geom_bar(stat = "identity", fill = "#2c7fb8") +
  labs(
    title = "Distribuição taxonômica das MAGs",
    subtitle = "Contagens a nível de filo",
    x = "Número de MAGs",
    y = NULL
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    axis.text.y = element_text(face = "italic")
  )
```

# Passos futuros

- Classificação das MAGs eucarióticas.
- Anotação funcional das montagens
  - Busca de clusters gênicos biosintéticos
  - Identificação de genes de resistência a antibióticos
  - Identificação de peptídeos antimicrobianos