# [Minicurso 2025] 05. Processamento com DADA2

###### tags: `Minicurso25`


DADA2 Website: https://benjjneb.github.io/dada2/

---

[TOC]

## Processamento:

**Objetivos:**

- Sair das **sequências** para os **microrganismos**
- Obter contagens de sequências (ASVs) em cada amostra
- Remover contamitantes ambientais e/ou provenientes de um hospedeiro

**Entrada:**

- Dados pré-processados com qualidade assegurada;
- Podem estar ou não fundidos (Aqui assume que estão!);
- Bancos de dados de taxonômia especializados no fragmento.

**Saída:**

- Tabelas de **contagens** e **taxonômias**.


## 1. Procedimentos iniciais

### Acesso ao RStudio server:
1. Acessar: http://200.145.102.175:8787/
2. Logar com suas credenciais de `login` e `senha`.

### Criar novo script:
- "File" > "New File" > "R Script"

### Limpeza do ambiente R:
- "Resete" o ambiente pressionando: `ctrl` + `shift` + `f10` (Ou vá em: "Session" > "Restart R");
- Limpe o ambiente digitando `rm(list = ls())` no **console** do R.

### Definindo local de trabalho:
Nesse primeiro trecho de código, vamos dizer ao R de onde gostariamos que ele trabalhasse:

# Script DADA2 ------------------------------------------------------------

# Procedimentos iniciais --------------------------------------------------

# Definindo ambiente de trabalho:
getwd()
setwd("~/projeto/analises/")
getwd()

# Criando um atalho para o caminho que iremos salvar as saídas:
caminho <- paste0(getwd(), "/dada2/")
caminho

### Carregando as bibliotecas que serão utilizadas:

# Carregando bibliotecas:
library("tidyverse")    # Manipulação de tabelas
library("dada2")        # Processamento de ASVs
library("phyloseq")     # Análises de microbioma

### Criando funções customizadas:
Nem sempre os pacotes irão fornecer funções que fazem ***exatamente*** o que você quer. Nesse caso, é comum a utilização de funções customizadas para propósitos específicos.

# Definindo funções customizadas:

# Função: Resumo das taxonomias
summary.tax <- function(ordered_names, ...){
  
  tmp_names <- ordered_names
  tmp_list <- list(...)
  
  tmp_summary1 <- data.frame(Database = NA, 
                             Kingdom = NA, 
                             Phylum = NA, 
                             Class = NA, 
                             Order = NA,
                             Family = NA, 
                             Genus = NA, 
                             Species = NA)
  
  for (i in seq(length(tmp_list))) {
    
    tmp_db <- tmp_list[[i]] %>% select(Kingdom:Species)
    tmp_name <- tmp_names[[i]]
    
    tmp_ksum <- sum(table(tmp_db["Kingdom"]))
    tmp_psum <- sum(table(tmp_db["Phylum"]))
    tmp_csum <- sum(table(tmp_db["Class"]))
    tmp_osum <- sum(table(tmp_db["Order"]))
    tmp_fsum <- sum(table(tmp_db["Family"]))
    tmp_gsum <- sum(table(tmp_db["Genus"]))
    tmp_ssum <- sum(table(tmp_db["Species"]))
    
    tmp_summary1 <- tmp_summary1 %>% add_row(Database = tmp_name, 
                                             Kingdom = tmp_ksum, 
                                             Phylum = tmp_psum, 
                                             Class = tmp_csum, 
                                             Order = tmp_osum,
                                             Family = tmp_fsum, 
                                             Genus = tmp_gsum, 
                                             Species = tmp_ssum)
    
  }
  
  tmp_summary1 <- tmp_summary1[-1,]
  
  tmp_summary2 <- data.frame(Database = NA,
                             UnqKingdom = NA, 
                             UnqPhylum = NA, 
                             UnqClass = NA, 
                             UnqOrder = NA,
                             UnqFamily = NA, 
                             UnqGenus = NA, 
                             UnqSpecies = NA)
  
  for (i in seq(length(tmp_list))) {
    
    tmp_db <- tmp_list[[i]] %>% select(Kingdom:Species)
    tmp_name <- tmp_names[[i]]
    
    tmp_ksum_u <- sum(table(unique(tmp_db["Kingdom"])))
    tmp_psum_u <- sum(table(unique(tmp_db["Phylum"])))
    tmp_csum_u <- sum(table(unique(tmp_db["Class"])))
    tmp_osum_u <- sum(table(unique(tmp_db["Order"])))
    tmp_fsum_u <- sum(table(unique(tmp_db["Family"])))
    tmp_gsum_u <- sum(table(unique(tmp_db["Genus"])))
    tmp_ssum_u <- sum(table(unique(tmp_db["Species"])))
    
    tmp_summary2 <- tmp_summary2 %>% add_row(Database = tmp_name, 
                                             UnqKingdom = tmp_ksum_u, 
                                             UnqPhylum = tmp_psum_u, 
                                             UnqClass = tmp_csum_u, 
                                             UnqOrder = tmp_osum_u, 
                                             UnqFamily = tmp_fsum_u, 
                                             UnqGenus = tmp_gsum_u, 
                                             UnqSpecies = tmp_ssum_u)
    
  }
  
  tmp_summary2 <- tmp_summary2[-1,]
  
  return(list(tmp_summary1, tmp_summary2))
  
}

# Função: Parseamento das taxonomias
fix_tax <- function(tax) {
  
  tmp <- as.data.frame(tax) %>% dplyr::mutate(across(Kingdom:Species, ~ str_replace_na(., "Unclassified")))
  
  tmp[2] <- str_remove_all(tmp[[2]],"k__")
  tmp[3] <- str_remove_all(tmp[[3]],"p__")
  tmp[4] <- str_remove_all(tmp[[4]],"c__")
  tmp[5] <- str_remove_all(tmp[[5]],"o__")
  tmp[6] <- str_remove_all(tmp[[6]],"f__")
  tmp[7] <- str_remove_all(tmp[[7]],"g__")
  tmp[8] <- str_remove_all(tmp[[8]],"s__")
  
  tmp[2] <- paste0("k__",tmp[[2]])
  tmp[3] <- paste0("p__",tmp[[3]])
  tmp[4] <- paste0("c__",tmp[[4]])
  tmp[5] <- paste0("o__",tmp[[5]])
  tmp[6] <- paste0("f__",tmp[[6]])
  tmp[7] <- paste0("g__",tmp[[7]])
  tmp[8] <- paste0("s__",tmp[[8]]) 
  
  names(tmp) <- c("ID", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
  
  tmp[grep("_Unclassified$", tmp$Species), 8] <- "s__Unclassified"
  
  tmp$Species <- sub("s__\\(.*", NA, tmp[,8])
  tmp$Species <- gsub("\\(.*", "", tmp[,8])
  
  return(tmp)
  
}

## 2. Processamento com DADA2

### Carregamento das Bibliotecas:
Bibliotecas de reads pré-processadas e fundidas que se encontram no diretório de saída do programa `flash`.

# Carregando bibliotecas:
reads <- list.files("./flash/", pattern = ".fastq", full.names=TRUE)
reads


### Separando nomes das amostras:

# Separando nome das amostras:
nomes <- sapply(strsplit(basename(reads), "_joined"), `[`, 1)
nomes

### Criando caminhos para bibliotecas filtradas:

# Criando caminhos para bibliotecas filtradas:
reads_filt <- file.path(caminho, "filt", paste0(nomes, "_filt.fastq"))
names(reads_filt) <- nomes
reads_filt

### (Mais um) Controle de qualidade:

Vamos dar uma breve visualizada nas qualidades dos dados:

# QC

plotQualityProfile(reads[1:4]) # Visualização pré-QC
plotQualityProfile(reads[21:24])


Agora, o CQ propriamente dito:

qc <- filterAndTrim(fwd = reads, 
                    filt = reads_filt, 
                    truncLen = 250, 
                    maxEE = 2, 
                    maxN = 0, 
                    rm.phix = T, 
                    compress = F)
qc

Note o valor de truncagem baseado no valor que obtivemos com o relatório do `USEARCH`.

:::danger
Com sequências de **ITS** (Fungos), é desaconselhável realizar a **truncagem**, pois há grande variação de tamanho nos diferentes grupos de fungos.
:::

Olhando, novamente, as qualidades:

plotQualityProfile(reads_filt[1:4]) # Visualização pós-QC
plotQualityProfile(reads_filt[21:24])


### Correção dos Erros:

- Diferencial do DADA2: Algoritmo de correção de erros;
- Processo de *machine-learning* que estima as substituições de bases mais recorrentes em cada qualidade;
- É baseado nos próprios dados. Iteração até achar melhor resultado.

**Relembrando: ASVs vs. OTUs**



:::warning
**Copie o código, mas não rode! O processo de correção de erros é uma etapa lenta e não será possível executar em tempo real.**
:::

# Aprendizado:
#erros <- learnErrors(reads_filt, multithread = T) # Demorado: All = ~6 min; 2T = ~17 min; 1T = ~30 min

# Visualização:
plotErrors(erros, nominalQ = TRUE)

# Dereplicação:
#reads_derep <- derepFastq(reads_filt) 

# Denoising:
#reads_dada <- dada(reads_derep, err = erros, multithread = T) # Razoavelmente demorado: All = ~1.5 min; 2T = ~4 min; 1T = ~7 min

:::success
- **Pegando um atalho**
Copie e cole a linha a seguir no console do R:

    load("/srv/cibioinfo/checkpoints/checkpoint_1.RData")

:::

### Finalizando o processamento:

- Quantificando ASVS

# Obtendo contagens de ASVs:
asvs <- makeSequenceTable(reads_dada)
dim(asvs)

- Removendo quimeras :lion_face::snake:
O método `dada` corrige os erros de substituição e indel, mas sequências quiméricas não. Portanto, há um método específico para a remoção dessas sequências, cujas extremidades

# Removendo quimeras:
asvs_nochim <- removeBimeraDenovo(asvs, method = "consensus", multithread = T) # Rapido: 1T = ~3.5 min
dim(asvs_nochim)

- Resumo das contagens

# Resumo dos processos:

# Função para contar ASVs por amostra
getN <- function(x) sum(getUniques(x))

# Resumo dos processamentos no DADA
cont_dada <- data.frame(cbind(nomes, qc[,2], sapply(reads_dada, getN), rowSums(asvs_nochim)))
cont_dada

colnames(cont_dada) <- c("ID", "Filt", "Denoised", "ASVs")
rownames(cont_dada) <- NULL

# Carregando resumo do pre-processamento:
cont_preproc <- read.table("../info/resumo_preproc.txt", header = T, sep = "\t")
cont_preproc

# Fundindo os resumos preproc + DADA:
cont_proc <- merge(cont_preproc, cont_dada, by = "ID") %>% mutate_at(-1, as.numeric)
cont_proc

# Salvando tabela
write.table(cont_proc, file = "../info/resumo_proc.txt", sep = "\t", col.names = T, row.names = F)

## 3. Atribuição taxonômica

- Designar microrganismo a qual a sequênica pertence;
- Contraste contra bancos de dados de **sequências-referência**;
- Capacidade da atribuição é variável - Região, tamanho, qualidade, grupo, etc...
- Bancos já treinados ([**Link**](https://benjjneb.github.io/dada2/training.html)) ou customizados.

  

:::warning
**Copie o código, mas não rode! O processo de anotação taxonômica é uma etapa lenta e não será possível executar em tempo real.**
:::

# Anotação taxonômica -----------------------------------------------------

# Processo demorado: ~30 min (All); ~5h (1T) [CADA]

# GTDB (Genome Taxonomy Database): 
#taxGTDB <- assignTaxonomy(asvs_nochim, "/srv/cibioinfo/tax_dbs/GTDB_bac120_arc122_ssu_r202_fullTaxo.fa.gz", multithread = T, tryRC = T)
#taxGTDB <- addSpecies(taxGTDB, "/srv/cibioinfo/tax_dbs/GTDB_bac120_arc122_ssu_r202_Species.fa.gz", tryRC = T)

# RDP (Ribosomal Database Project):
#taxRDP <- assignTaxonomy(asvs_nochim, "/srv/cibioinfo/tax_dbs/rdp_train_set_18.fa.gz", multithread = T, tryRC = T)
#taxRDP <- addSpecies(taxRDP, "/srv/cibioinfo/tax_dbs/rdp_species_assignment_18.fa.gz", tryRC = T)

# RefRDP (NCBI RefSeq 16S rRNA database supplemented by RDP):
#taxRefRDP <- assignTaxonomy(asvs_nochim, "/srv/cibioinfo/tax_dbs/RefSeq_16S_6-11-20_RDPv16_fullTaxo.fa.gz", multithread = T, tryRC = T)
#taxRefRDP <- addSpecies(taxRefRDP, "/srv/cibioinfo/tax_dbs/RefSeq_16S_6-11-20_RDPv16_Species.fa.gz", tryRC = T)

# SILVA: 
#taxSilva <- assignTaxonomy(asvs_nochim, "/srv/cibioinfo/tax_dbs/silva_nr99_v138.1_train_set.fa.gz", multithread = T, tryRC = T)
#taxSilva <- addSpecies(taxSilva, "/srv/cibioinfo/tax_dbs/silva_species_assignment_v138.1.fa.gz", tryRC = T)

:::success
- **Pegando um atalho**
Copie e cole a linha a seguir no console do R:

    load("/srv/cibioinfo/checkpoints/checkpoint_2.RData")

:::

- Vamos visualizar a estrutura dos arquivos de taxonômias:

# Visualização das taxonomias:
View(taxGTDB)
View(taxRDP)
View(taxRefRDP)
View(taxSilva)

## 4. Filtragens e Transformações das tabelas

- Agora temos em mão todos os arquivos essenciais para a análise do microbioma. Contudo, ainda são necessários alguns ajustes:

- Filtragem de contaminantes;
- Seleção do banco de taxonomia;
- Formatação das tabelas.

- Vamos começar por padronizar as tabelas de taxonomia:

# Transformações das tabelas de taxonomia:
taxGTDB_df <- taxGTDB %>% 
  data.frame() %>% # Transformando de matriz para data frame
  mutate(seq = rownames(.)) %>% # Criando uma coluna para a ASV
  select(seq, Kingdom:Species) %>% # Alterando a ordem das colunas
  mutate(Species = str_remove_all(Species, "\\(.*$")) %>%  # Removendo alguns apêndices das espécies. Ids, cepas...
  mutate(Species = str_replace_all(Species, " ", "_")) %>% # Trocando espaços por "_" 
  mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) # Trocando células vazias por "NA"

taxRDP_df <- taxRDP %>% 
  data.frame() %>% 
  mutate(seq = rownames(.)) %>%
  mutate(Species2 = paste0(Genus, "_", Species)) %>% 
  mutate(Species2 = str_remove_all(Species2, ".*NA$")) %>%
  mutate(Species = Species2) %>% 
  select(seq, Kingdom:Species) %>% 
  mutate(across(Kingdom:Species, ~ gsub("^$", NA, .)))

taxRefRDP_df <- taxRefRDP %>% 
  data.frame() %>% 
  mutate(seq = rownames(.)) %>% 
  select(seq, Kingdom:Species) %>% 
  mutate(across(Kingdom:Species, ~ gsub("^$", NA, .))) %>% 
  mutate(Species = str_remove_all(Species, "\\(.*$")) %>% 
  mutate(Species = str_replace_all(Species, " ", "_"))

taxSilva_df <- taxSilva %>% 
  data.frame() %>% 
  mutate(seq = rownames(.)) %>%
  mutate(Species2 = paste0(Genus, "_", Species)) %>% 
  mutate(Species2 = str_remove_all(Species2, ".*NA$")) %>%
  mutate(Species = Species2) %>% 
  select(seq, Kingdom:Species) %>% 
  mutate(across(Kingdom:Species, ~ gsub("^$", NA, .)))

- Deixando a tabela de contagens no mesmo padrão:

# Obtendo tabela de contagem:
contagens_df <- data.frame(t(asvs_nochim)) %>% 
  mutate(seq = rownames(.)) %>% 
  select(seq, everything())

- Fundindo contagens e taxonomias:

# Fundindo tudo:
taxAll_cont_df <- 
  merge(taxGTDB_df %>% mutate(GTDB = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, GTDB), 
        merge(taxRDP_df %>% mutate(RDP = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, RDP),
              merge(taxRefRDP_df %>% mutate(RefRDP = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, RefRDP),
                    merge(taxSilva_df %>% mutate(Silva = paste(Kingdom, Phylum, Class, Order, Family, Genus, Species, sep = "; ")) %>% select(seq, Silva),
                          contagens_df, 
                          by = "seq"), 
                    by = "seq"), 
              by = "seq"), 
        by = "seq")

- Filtrando contaminantes:

# Estabelecendo contaminantes:
contaminantes <- c("Chloroplast;", "Mitochondria;", "^Unclassified;", "^Eukaryota;")

# Extraindo seqs. não contaminantes:
seqs_boas <- taxAll_cont_df %>% 
  filter_at(vars(-seq), all_vars(str_detect(., paste(contaminantes, collapse = "|"), negate = T))) %>% 
  pull(seq)

# Filtrando seqs. contaminantes:
taxGTDB_filt_df <- taxGTDB_df %>% filter(seq %in% seqs_boas)
taxRDP_filt_df <- taxRDP_df %>% filter(seq %in% seqs_boas)
taxRefRDP_filt_df <- taxRefRDP_df %>% filter(seq %in% seqs_boas)
taxSilva_filt_df <- taxSilva_df %>% filter(seq %in% seqs_boas)
contagens_filt_df <- contagens_df %>% filter(seq %in% seqs_boas)

- Selecionando banco de taxonomia:

# Comparação das classificações - Função customizada:
ClassTax <- summary.tax(ordered_names = c("GTDB", "RDP", "RefRDP", "Silva"), 
                           taxGTDB_filt_df, taxRDP_filt_df, taxRefRDP_filt_df, taxSilva_filt_df)

# Seqs. classificadas até determinado nível:
ClassTax_cont <- ClassTax[[1]]
ClassTax_cont

# Táxons únicos até determinado nível:
ClassTax_unic <- ClassTax[[2]]
ClassTax_unic

- Visualizando os melhores:

# Visualizando os melhores:
View(taxGTDB_filt_df)
View(taxSilva_filt_df)

- **Pontos a serem levados em conta:**
*- Quantidade;*
*- Qualidade;*
*- Precisão;*
*- Atualização;*
*- Aceitação*

- Construção da tabela "mestra":
*Tabela que contém todas as informações centralizadas.*

# Construção da "master table":
master <- merge(taxSilva_filt_df, contagens_filt_df) %>%  # Fusão das taxonomias e contagens
  mutate("CountSum" = rowSums(.[9:ncol(.)])) %>% # Somatória de abundancia de cada ASV
  mutate("CountPrev" = rowSums(.[9:ncol(.)] != 0) - 1) %>%  # Somatória de prev. nas amostras (descontando a somatória)
  arrange(-CountSum) %>% # Ordenando tabela pela abundancia
  mutate("ID" = paste0("ASV_", seq(1, nrow(.)))) %>% # Adicionando identificador único
  select(ID, everything()) # Reordenando...

- Obtendo algumas informações finais...

# Contagens finais:
cont_proc$ASVs_clean <- master %>% select(MP_BS_H_1:VD_RZ_L_3) %>% colSums()
cont_proc$Usable_perc <- round((cont_proc$ASVs_clean / cont_proc$Raw) * 100, 2)
View(cont_proc)

write.table(cont_proc, file = "../info/resumo_proc.txt", sep = "\t", col.names = T, row.names = F)

# Porcentagens de classificações:
sum(master %>% filter(!is.na(Family)) %>% pull(CountSum)) / sum(master$CountSum) # % fam
sum(master %>% filter(!is.na(Genus)) %>% pull(CountSum)) / sum(master$CountSum)  # % gen
sum(master %>% filter(!is.na(Species)) %>% pull(CountSum)) / sum(master$CountSum)  # % sp

- Tabelas formatadas para o [MicrobiomeAnalyst](https://dev.microbiomeanalyst.ca/MicrobiomeAnalyst/upload/OtuUploadView.xhtml):
*O site pede uma formatação específica, que iremos atender a seguir:*

# Salvando tabelas para o MicrobiomeAnalyst:
tb_meta <- read.delim(file = "../info/metadados.txt", header = T, sep = "\t") %>%
  arrange(Sample_name) %>% 
  mutate(Group = str_remove(Sample_name, "_[0-9]$")) %>% 
  select(Sample_name, Group, everything()) %>% 
  rename("#NAME" = "Sample_name")

tb_cont <- master %>% 
  select(ID, starts_with("MP_"), starts_with("VD_")) %>% 
  rename("#NAME" = "ID")

tb_taxa <- master %>% 
  select(ID, Kingdom:Species) %>% 
  fix_tax() %>% # Função customizada
  rename("#TAXONOMY" = "ID")

View(tb_meta)
View(tb_cont)
View(tb_taxa)

write.table(tb_cont, file = paste0(caminho, "MA01_contagens.txt"), quote = F, sep = "\t", col.names = T, row.names = F)
write.table(tb_meta, file = paste0(caminho, "MA02_metadados.txt"), quote = F, sep = "\t", col.names = T, row.names = F)
write.table(tb_taxa, file = paste0(caminho, "MA03_taxonomias.txt"), quote = F, sep = "\t", col.names = T, row.names = F)

# Salvando tabelas para o MicrobiomeAnalyst (Funcional):
tb_meta_fun <- read.delim(file = "../info/metadados.txt", header = T, sep = "\t") %>%
  arrange(Sample_name) %>% 
  mutate(Group = str_remove(Sample_name, "_[0-9]$")) %>% 
  select(Sample_name, Group, everything()) %>% 
  rename("#NAME" = "Sample_name")

tb_cont_fun <- master %>% 
  select(ID = seq, starts_with("MP_"), starts_with("VD_")) %>% 
  rename("#NAME" = "ID")

tb_taxa_fun <- master %>% 
  select(ID = seq, Kingdom:Species) %>% 
  fix_tax() %>% # Função customizada
  rename("#TAXONOMY" = "ID")

View(tb_meta_fun)
View(tb_cont_fun)
View(tb_taxa_fun)

write.table(tb_cont_fun, file = paste0(caminho, "funMA01_contagens.txt"), quote = F, sep = "\t", col.names = T, row.names = F)
write.table(tb_meta_fun, file = paste0(caminho, "funMA02_metadados.txt"), quote = F, sep = "\t", col.names = T, row.names = F)
write.table(tb_taxa_fun, file = paste0(caminho, "funMA03_taxonomias.txt"), quote = F, sep = "\t", col.names = T, row.names = F)

## pacote 'phyloseq'
O pacote 'phyloseq' é uma ferramenta para importar, armazenar, analisar e exibir graficamente dados de sequenciamento de microbiomas.

**Tutorial:** https://vaulot.github.io/tutorials/Phyloseq_tutorial.html#step-by-step

# Perdeu algo? Carregue as tabelas finais:
# Perdeu algo? Carregue as tabelas finais:
load("/srv/cibioinfo/checkpoints/checkpoint_3.RData")

library("RColorBrewer") 

# Criando objeto 'phyloseq'
row.names(tb_cont) <- tb_cont$`#NAME`
row.names(tb_taxa) <- tb_taxa$`#TAXONOMY`
row.names(tb_meta) <- tb_meta$`#NAME`

tb_cont <- tb_cont %>% select(-`#NAME`)
tb_taxa <- tb_taxa %>% select(-`#TAXONOMY`)
tb_meta <- tb_meta %>% mutate(Sample = factor(`#NAME`, levels = unique(pull(., `#NAME`)))) %>% select(-`#NAME`)

tb_cont <- as.matrix(tb_cont)
tb_taxa <- as.matrix(tb_taxa)

ps_cont <- otu_table(tb_cont, taxa_are_rows = TRUE)
ps_taxa <- tax_table(tb_taxa)
ps_meta <- sample_data(tb_meta)

ps_table <- phyloseq(ps_cont, ps_taxa, ps_meta)
ps_table

# Gráficos:

# Div. alfa:
plot_richness(ps_table, measures=c("Chao1", "Shannon", "Simpson"), x = "Group", color = "Trial")

# Div. beta:
bray <- ordinate(ps_table, method = "PCoA", distance="bray")

plot_ordination(ps_table, bray, color="Trial", shape = "Env") + 
  geom_point(size = 3)

# Perfil taxonômico:
plot_bar(ps_table, fill = "Phylum")

# Perfil taxonômico relativo:
ps_table_rel <- transform_sample_counts(ps_table, function(x) x/sum(x)*100) # Transforma abund. absoluta em relativa
ps_table_rel

plot_bar(ps_table_rel, fill = "Phylum")

ntaxa(ps_table_rel) # n. de taxa

ps_table_rel_glom_phylum <- tax_glom(ps_table_rel, "Phylum") # Soma contagens em um determinado nível taxonomico
ps_table_rel_glom_phylum

ntaxa(ps_table_rel_glom_phylum) # n. de taxa aglomerado

# Default
plot_bar(ps_table_rel_glom_phylum, fill = "Phylum")

# Paletas de cores pré definidas
plot_bar(ps_table_rel_glom_phylum, fill = "Phylum") + 
  scale_fill_manual(values = colorRampPalette(brewer.pal(9, "Set1"))(ntaxa(ps_table_rel_glom_phylum)))

# ...ou customizar ainda mais
plot_bar(ps_table_rel_glom_phylum, fill = "Phylum") + 
  scale_fill_manual(values = 
                      colorRampPalette(
                      c("#DA4453FF", "#E95546FF", "#F6BA59FF", "#8BC163FF", 
                        "#34BC9DFF", "#3BB0D6FF", "#4B8AD6FF", "#977BD5FF", 
                        "#D870A9FF", "#E6E9EDFF", "#AAB2BCFF", "#434A53FF"))(ntaxa(ps_table_rel_glom_phylum))
                    ) +
  theme_bw() + 
  facet_grid(. ~ Trial + Env, scales = "free_x") + 
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(
    legend.title = element_text(face = "bold"),
    strip.text.x = element_text(face = "bold.italic"),
    axis.text.x = element_text(angle = 90, vjust = .5, hjust = 0, size = 6),
    axis.title.y = element_text(face = "bold", margin = margin(r = 10)),
    axis.title.x = element_text(face = "bold", margin = margin(t = 10))
  )

# Apenas os Tops
top_taxa <- names(sort(taxa_sums(ps_table_rel_glom_phylum), TRUE)[1:10])
top_taxa

ps_top_taxa <- prune_taxa(top_taxa, ps_table_rel_glom_phylum) # seleciona de acordo com um vetor
ps_top_taxa

ntaxa(ps_top_taxa) # n. de taxa

plot_bar(ps_top_taxa, fill = "Phylum") + 
  scale_fill_manual(values = 
                      colorRampPalette(
                        c("#DA4453FF", "#E95546FF", "#F6BA59FF", "#8BC163FF", 
                          "#34BC9DFF", "#3BB0D6FF", "#4B8AD6FF", "#977BD5FF", 
                          "#D870A9FF", "#E6E9EDFF", "#AAB2BCFF", "#434A53FF"))(ntaxa(ps_top_taxa))
  ) +
  theme_bw() + 
  facet_grid(. ~ Trial + Env, scales = "free_x") + 
  scale_x_discrete(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100), breaks = seq(0, 100, 10)) +
  theme(
    legend.title = element_text(face = "bold"),
    strip.text.x = element_text(face = "bold.italic"),
    axis.text.x = element_text(angle = 90, vjust = .5, hjust = 0, size = 6),
    axis.title.y = element_text(face = "bold", margin = margin(r = 10)),
    axis.title.x = element_text(face = "bold", margin = margin(t = 10))
  )