microDwelleR.Rmd

---
title: "Microbiome convergence and deterministic community assembly along successional biocrust gradients on potash salt heaps"
author: "Roberto Siani"
date: "14.2.23"
---

DISCLAIMER:
This script is highly opinionated (choice of tools, normalization methods, decontamination methods, palettes, themes ...). It's up to you to use it as it is or modify it as it best suits you (implying, you will have to look into the code and make those changes).

# PRE-PROCESS

Installing and loading the required packages and processing amplicon sequences into a convenient phyloseq object

```{r, cache = TRUE, echo = FALSE, include = FALSE}

# pacman to install and load libraries

pacman::p_load(
        tidyverse,
        patchwork,
        pacman
)

# a minimal theme

my_theme <- theme_minimal() +
  theme(
    line = element_line(linewidth = 1.5),
    plot.margin = margin(1, 1, 1, 1, unit = "cm"),
    text = element_text(size = 16,
                        family = "Arial",
                        color = "#212121"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    legend.justification = c(1, 0),
    legend.margin = margin(0, 5, 3, 5),
    legend.title = element_blank(),
    panel.spacing.x = unit(3, "lines"),
    panel.spacing.y = unit(3, "lines"),
    axis.line = element_blank(),
    panel.border = element_rect(color = "#212121", linewidth = 1.5 , fill = NA),
    axis.text.x = element_text(size = 16,
                        family = "Arial",
                        color = "#212121"),
    axis.text.y = element_text(size = 16,
                        family = "Arial",
                        color = "#212121"),
    axis.title.x = element_text(size = 16,
                        family = "Arial",
                        color = "#212121", 
                        hjust = 1),
    axis.title.y = element_text(size = 16,
                        family = "Arial",
                        color = "#212121", 
                        hjust = 1),
    strip.text = element_text(size = 16,
                        family = "Arial",
                        color = "#212121")
  )

theme_set(my_theme)


# install/load the remainder of libraries

pacman::p_load(
  phyloseq,
  microbiome)

# specify the file path for QIIME2 results and the negative control samples

metadata = "halden_mapping_run_no.txt"
taxonomy = "taxonomy-dada2_240200e4.qza"
abundance = "table_all.qza"
ASVs = "rep-seqs_all.qza"
control_samples = c("Juliette14",
                    "Juliette20",
                    "Juliette21",
                    "Juliette24",
                    "Juliette25",
                    "Juliette27",
                    "Juliette28",
                    "Juliette29",
                    "Juliette30")

# import QIIME2 results, clean them and integrate in a phyloseq object

# metadata (mapping file). 
# I filter out 3 samples that were excluded due to sequencing issue

metaTab = read.delim2(metadata,
                      row.names = 1,
                      colClasses = "factor") %>% 
  select(!c(BarcodeSequence,
            LinkerPrimerSequence,
            ReversePrimer,
            Description)) %>% 
  filter(!Sample_num %in% "NTC1" &
           !Sample_num %in% "D3" &
           !Sample_num %in% "D7")

# taxonomy ("qiime feature-classifier")
# the taxonomy is parsed to a human readable format
# missing values are replaced by "U." (unknown) + last known taxonomy assignment
# e.g. U. Proteobacteriaceae

p_taxTab =  qiime2R::read_qza(taxonomy)$data %>%
  separate(Taxon,
           sep = ";",
           c("Domain",
             "Phylum",
             "Class",
             "Order",
             "Family",
             "Genus",
             "Species"),
           fill = "right",
           extra = "drop") %>%
  as.matrix() %>%
  gsub(x = .,
       pattern = "[A-Z]_[0-9]__",
       replacement = "") %>% 
  as.data.frame() %>%
  map_df(~ gsub(
    pattern = "metagenome|uncultured|unidentified|Unknown",
    replacement = NA,
    .x)) %>%
  mutate_if(is_character, str_trim) %>%
  column_to_rownames("Feature.ID") %>%
  mutate(Domain = ifelse(is.na(Domain),
                         "U. Domain",
                         Domain),
         Phylum = coalesce(Phylum,
                           ifelse(grepl("^U.", Domain),
                                  Domain,
                                  paste("U.", Domain))),
         Class = coalesce(Class,
                          ifelse(grepl("^U.", Phylum),
                                 Phylum,
                                 paste("U.", Phylum))),
         Order = coalesce(Order,
                          ifelse(grepl("^U.", Class),
                                 Class,
                                 paste("U.", Class))),
         Family = coalesce(Family,
                           ifelse(grepl("^U.", Order),
                                  Order,
                                  paste("U.", Order))),
         Genus = coalesce(Genus,
                          ifelse(grepl("^U.", Family),
                                 Family,
                                 paste("U.", Family))),
         Species = coalesce(Species,
                            ifelse(grepl("^U.", Genus),
                                   Genus,
                                   paste("U.", Genus))))


# ASVs counts ("qiime dada2 denoise-paired")

p_abuTab = qiime2R::read_qza(abundance)$data %>%
  as.data.frame()

# de-contamination and remove controls DOI: 10.1186/s40168-018-0605-2 
# filtering out chloroplast, mitochondria and non bacteria/archaea

c_abuTab = p_abuTab[!decontam::isContaminant(t(p_abuTab),
                                   neg = c(colnames(p_abuTab)
                                           %in% control_samples) ,
                                   method = "prevalence",
                                   threshold = 0.05,
                                   batch = metaTab$Run,
                                   normalize = T,
                                   detailed = F), ] %>%
  filter(rownames(.) %in%
           rownames(p_taxTab %>%
                      filter(!Order %in% "Chloroplast" &
                               !Order %in% "Mitochondria" &
                               !Domain %in% "Unassigned"))) %>%
  select(!all_of(control_samples))  %>%
  filter(rowSums(.) > 0)

# normalization by scaling with ranked subsampling
# DOI: 10.7717/peerj.9593

abuTab = SRS::SRS(c_abuTab,
             Cmin = min(colSums(c_abuTab))) %>%
  `rownames<-`(rownames(c_abuTab))

# filter taxonomy file from ASVs removed by the cleaning

taxTab = p_taxTab %>%
  filter(rownames(.) %in%
           rownames(abuTab)) %>%
  as.matrix()

# representative sequences

seqsTab = qiime2R::read_qza(ASVs)$data[rownames(abuTab)]

# compute phylogenetic tree

treeObj = 
  phangorn::phyDat(
    as(
      DECIPHER::AlignSeqs(seqsTab,
                      anchor = NA,
                      verbose = F),
            "matrix"),
         type = "DNA") %>%
  phangorn::dist.ml(.) %>% 
  phangorn::NJ(.)

# import into phyloseq-object

phytted =  phyloseq(otu_table(abuTab,
                              taxa_are_rows = T),
                    phyloseq::tax_table(taxTab),
                    sample_data(metaTab),
                    refseq(seqsTab),
                    phy_tree(treeObj))

# your metataxonomy at a glance

summarize_phyloseq(phytted)

# save phyloseq-obj as .RData to avoid repeating the lengthy computation steps

save(phytted,
     file = paste("phytted",
                  Sys.time(),
                  ".RData",
                  sep = ""))

```


# LOAD

Direct loading of the object after pre-processing is done

```{r, echo = FALSE, include = FALSE}

# load your phytted dataset

load("phytted2022-11-23 19:35:39.RData")

# again, take a good look at your results. Look at the sparsity and the percentage of singletons and filter if necessary

summarize_phyloseq(phytted)

# we filter anything belove 0.1 prevalence (10% samples)

phytted =
  phytted %>% 
  prune_taxa(microbiome::prevalence(., detection  = 0) >= 0.1, .)

summarize_phyloseq(phytted)

# relevel

sample_data(phytted)$type = factor(c(
  rep("biocrust", 3), rep("initial", 2), "biocrust", "initial", "biocrust", rep("initial", 3), rep("heap", 5), "biocrust"),
  levels = c("heap", "initial", "biocrust"))

# let's look at combinations

table(meta(phytted)$site,
      meta(phytted)$type)

# color palette

pal_type = c("#5CC8DF", "#FFCF04", "#008A61")


```

# OUTPUT

### Relative Abundance

```{r}

# select only top 10% taxa 

melted = phytted %>%
  metagMisc::phyloseq_filter_top_taxa(perc = 10) %>% 
  speedyseq::psmelt()

# plot relative abundance of most abundant phyla

plot_3 =
  ggplot(melted) +
  geom_bar(stat = "identity",
           position = "fill",
           color = "white",
           linewidth = 0.05,
           aes(x = type,
               y = Abundance,
               fill = Phylum),
           alpha = 0.9) +
  scale_fill_hue(l = 60, c = 60) +
  facet_wrap(~site,
             ncol = 4)

cairo_pdf("OUT/relativeAbundance.pdf",
    width = 16,
    height = 16,
    onefile = T)
plot_3
dev.off()

```

### Alpha Diversity

```{r}

# compute alpha diversity metric indexes

alpha_frame = 
  phytted %>%
  microbiome::alpha(
    index = c("chao1",
              "evenness_pielou", 
              "diversity_shannon",
              "dominance_simpson")) %>%
  merge(meta(phytted),
        by = 0) %>% 
  `colnames<-`(gsub("_",
                    " ",
                    colnames(.))) %>% 
  pivot_longer(cols = 2:5)

# write down test results

ggpubr::compare_means(
  formula = value ~ type,
  method = "wilcox.test",
  data = alpha_frame,
  group.by = c("name"),
  p.adjust.method = "fdr") %>% write_csv("compare_type.csv")

ggpubr::compare_means(
  formula = value ~ type,
  method = "wilcox.test",
  data = alpha_frame,
  group.by = c("name", "site"),
  p.adjust.method = "fdr") %>% write_csv("compare_type_site.csv")

# boxplots

div_list =
  alpha_frame %>% 
  group_by(name) %>% 
  nest() %>% 
  mutate(plotted = 
           map(data, 
               ~ ggplot(.x,
                        aes(x = type,
                            y = value,
                            fill = type)) +
                 scale_shape_manual(values = c(21, 23)) +
                 geom_boxplot(aes(fill = type),
                              outlier.shape = NA, alpha = 1,
                              linewidth = 1.5,
                              color = "#212121") +
                 scale_fill_manual(values = pal_type) +
                 theme(legend.position = "none",
                       axis.title.y = element_blank(),
                       axis.title.x = element_blank()) +
                 facet_wrap(~ site) +
                 ggtitle(name)))

div_list$plotted[[1]] + div_list$plotted[[2]] +
  div_list$plotted[[3]] + div_list$plotted[[4]] +
  plot_layout(ncol = 2, nrow = 2)

a = 
  ggplot(alpha_frame %>% 
         filter(
           name %in%
             c("chao1",
               "evenness pielou", 
               "diversity shannon",
               "dominance simpson")),
       aes(x = type,
           y = value,
           fill = type,
           color = type,
           title = type)) +
  geom_path(aes(linetype = site,
                group = site),
            stat = "summary", color = "#757575", linewidth = 1.5,
            position = position_dodge(0.5)) +
  scale_shape_manual(values = c(21, 23)) +
  geom_errorbar(stat = "summary",
                aes(group = site),
            position = position_dodge(0.5),
            color = "#757575",
            linewidth = 1.5,
            width = .1) +
  geom_point(aes(shape = site), stat = "summary", size = 9,
            position = position_dodge(0.5)) +
  scale_color_manual(values = pal_type, aesthetics = c("color", "fill")) +
  ggpubr::stat_compare_means(
    aes(label = after_stat(p.adj)),
    ref.group = "heap", 
    hide.ns = T,
    method.args = list(p.adjust.methods = "fdr"),
    method = "wilcox.test",
    size = 5,
    vjust = 1,
    family = "Arial"
  ) +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  facet_wrap(
    ~ name,
    scales = "free", ncol = 2) 

# save all in a pdf

tiff("OUT/fig3.tiff", width = 360, height = 360, units = "mm", res = 1000)
div_list$plotted[[1]] + div_list$plotted[[2]] +
  div_list$plotted[[3]] + div_list$plotted[[4]] +
  plot_layout(ncol = 2, nrow = 2)
dev.off()

```

### CCA
 
```{r}

## read additional metadata

sampling_environment = 
  read_tsv("20200823_abiotic.txt") %>% 
  filter(seq.id != "Juliette7") %>% 
  select(-c(`DON`, `TNb_g`, `Nitrit_g`, `Nitrat_g`, `Ammonium_g`, `water_perc`))


sampling_environment$type = factor(c("biocrust", "biocrust", "biocrust", "initial", "initial", "initial", "heap", "heap", "biocrust", "biocrust", "biocrust", "initial", "initial", "initial", "heap", "heap", "heap"),
                                   levels = c("heap", "initial", "biocrust"))

# scale chosen env variables

sampling_environment2 =
  sampling_environment %>% 
  select(6:10) %>% 
  mutate_all(scale) %>% 
  `colnames<-`(c("pH", "EC", "Chla", "TDN", "DOC"))

count_table =
  microbiome::abundances(phytted, 
             "compositional") %>% 
  as.data.frame()

# either all together or subsets

Y = 
  sampling_environment2 %>%
  `rownames<-`(sampling_environment$seq.id)

Y_OD = 
  sampling_environment2[9:17,] %>% 
  `rownames<-`(sampling_environment[9:17,]$seq.id)

Y_WT = 
  sampling_environment2[1:8,] %>%
  `rownames<-`(sampling_environment[1:8,]$seq.id)

X =
  count_table[, match(sampling_environment$seq.id, colnames(count_table))] %>% 
  t() %>% 
  as.data.frame() %>% 
  `colnames<-`(rownames(count_table))

X_OD =
  count_table[, match(sampling_environment[9:17,]$seq.id, colnames(count_table))] %>% 
  t() %>% 
  as.data.frame() %>% 
  `colnames<-`(rownames(count_table))

X_WT =
  count_table[, match(sampling_environment[1:8,]$seq.id, colnames(count_table))] %>% 
  t() %>% 
  as.data.frame() %>% 
  `colnames<-`(rownames(count_table))

## CCA

res_cca = 
  vegan::cca(X ~ ., data = Y, scale = T)

coordinates = 
  res_cca$CCA$wa %>% 
  as.data.frame() %>% 
  rownames_to_column("seq.id") %>% 
  left_join(sampling_environment,
            by = "seq.id")

coordinates$type = 
  factor(coordinates$type, 
         levels = c("heap", "initial", "biocrust"))

arrows = 
  res_cca$CCA$biplot %>% 
  as.data.frame() %>%
  rownames_to_column("variable")

plot_1 =
  ggplot(data = coordinates) +
  geom_point(aes(x = CCA1,
                 y = CCA2,
                 color = type,
                 shape = site),
             size = 9) +
  scale_color_manual(values = pal_type) +
  geom_segment(data = arrows,
               aes(x = 0,
                   y = 0,
                   xend = CCA1,
                   yend = CCA2),
               arrow = arrow(length = unit(3, "pt")),
               linejoin = "mitre",
               lineend = "butt",
               linewidth = 1.5,
               color = "#757575") +
  ggrepel::geom_text_repel(data = arrows,
            aes(x = CCA1,
                y = CCA2,
                label = variable),
            force = 2,
            size = 5,
            color = "#333333") +
  scale_shape_manual(values = c(18, 20))

res_cca1 = 
  vegan::cca(X_OD ~ ., data = Y_OD, scale = T)

coordinates1 = 
  res_cca1$CCA$wa %>% 
  as.data.frame() %>% 
  rownames_to_column("seq.id") %>% 
  left_join(sampling_environment,
            by = "seq.id")

coordinates1$type = 
  factor(coordinates1$type, 
         levels = c("heap", "initial", "biocrust"))

arrows1 = 
  res_cca1$CCA$biplot %>% 
  as.data.frame() %>%
  rownames_to_column("variable")

plot_1b =
  ggplot(data = coordinates1)  +
  geom_point(aes(x = CCA1,
                 y = CCA2,
                 color = type),
             size = 9) +
  scale_color_manual(values = pal_type) +
  geom_segment(data = arrows,
               aes(x = 0,
                   y = 0,
                   xend = CCA1,
                   yend = CCA2),
               arrow = arrow(length = unit(3, "pt")),
               linejoin = "mitre",
               lineend = "butt",
               linewidth = 1.5,
               color = "#757575") +
  ggrepel::geom_text_repel(data = arrows,
            aes(x = CCA1,
                y = CCA2,
                label = variable),
            force = 2,
            size = 5,
            color = "#333333") +
  scale_shape_manual(values = c(18))

res_cca2 = 
  vegan::cca(X_WT ~ ., data = Y_WT, scale = T)

coordinates2 = 
  res_cca2$CCA$wa %>% 
  as.data.frame() %>% 
  rownames_to_column("seq.id") %>% 
  left_join(sampling_environment,
            by = "seq.id")

coordinates2$type = 
  factor(coordinates2$type, 
         levels = c("heap", "initial", "biocrust"))

arrows2 = 
  res_cca2$CCA$biplot %>% 
  as.data.frame() %>%
  rownames_to_column("variable")

plot_1c =
  ggplot(data = coordinates2)  +
  geom_point(aes(x = CCA1,
                 y = CCA2,
                 color = type),
             size = 9) +
  scale_color_manual(values = pal_type) +
  geom_segment(data = arrows,
               aes(x = 0,
                   y = 0,
                   xend = CCA1,
                   yend = CCA2),
               arrow = arrow(length = unit(3, "pt")),
               linejoin = "mitre",
               lineend = "butt",
               linewidth = 1.5,
               color = "#757575") +
  ggrepel::geom_text_repel(data = arrows,
            aes(x = CCA1,
                y = CCA2,
                label = variable),
            force = 2,
            size = 5,
            color = "#333333") +
  scale_shape_manual(values = c(20))

list_plots  = list(plot_1, plot_1b, plot_1c)


tiff("OUT/fig2.tiff", 360, 360, units = "mm", res = 1000,)
plot_1
dev.off()

## checking for phylogenetic signal

summary(res_cca)

#read in otu table from phyloseq object 
otu = vegan::scores(res_cca, choices = c(1:5))$species

## read in the phylogeny
phylo = phy_tree(phytted);
phylo; # a summary of the phylogeny
#plot.phylo(phylo,typ="fan"); # a quick plot
rooted.phy = ape::root(phy_tree(phytted), 582, resolve.root = T)


res.phylosignal_cca = 
  picante::multiPhylosignal(otu[rooted.phy$tip.label,], rooted.phy, reps = 999, checkdata=TRUE)


write.csv(res.phylosignal_cca, "multiPhylosignal.csv")

```

## PCoA

```{r}

## make unique names for the genera. Will come handy for the plotting

new_tax_table = 
  phytted@tax_table@.Data %>% 
  as.data.frame() %>% 
  mutate(Genus = vctrs::vec_as_names(Genus, repair = "unique") %>% 
           sub("...", "_", fixed = T, .))

phytted_id = 
  phyloseq(
    otu_table(phytted@otu_table, taxa_are_rows = T),
    sample_data(phytted@sam_data),
    tax_table(new_tax_table %>% 
                as.matrix),
    phy_tree(phytted@phy_tree))

saveRDS(phytted_id, "phytted_id.RDS")

count_table =
  abundances(phytted) %>% 
  as.data.frame()

# PCoA using unifrac distance

set.seed(1)
res_pcoa =
  ordinate(phytted,
           "PCoA",
           "unifrac",
           weighted = T)

covar = as.factor(paste(meta(phytted)$site, meta(phytted)$type))

# we first select a subset of taxa using a univariate filter

set.seed(1)

kruskal_filter = 
  map_dbl(t(count_table) %>% as.data.frame,
          ~ kruskal.test(.x ~ covar, data = count_table)$p.value) %>% 
  p.adjust("BH")


# then we deploy a random forest classification to extract the variable importance

p_load(caret)
seedlings = vector(mode = "list", length = 18)
set.seed(3)
for (i in 1:17) seedlings[[i]] <- sample.int(1000, 15)  
seedlings[[18]] <- sample.int(1000, 1)

res_rf =
  train(
    x = t(count_table[kruskal_filter <= 0.05,]),
    y = covar,
    preProcess = "center",
    tuneLength = 15,
    trControl = trainControl(method = "cv",
                             number = 999, 
                             seeds=  seedlings))

# we extract all the results for the plots

varImp_rf =
  varImp(res_rf) %>%
  pluck("importance") %>%
  rownames_to_column("ASVs") 

samples =
  plot_ordination(transform(phytted, "compositional"),
                  res_pcoa,
                  type = "samples",
                  justDF = T)

species =
  plot_ordination(transform(phytted, "compositional"),
                  res_pcoa,
                  type = "species",
                  justDF = T) %>% 
  rownames_to_column("ASVs") %>% 
  left_join(varImp_rf) %>% 
  replace_na(list(Overall = 0)) %>% 
  mutate(qval_KW = kruskal_filter,
         zscore = (Overall - mean(Overall))/sd(Overall))

# the first plot is an ordination

a = ggplot(species,
       aes(x = Axis.1,
           y = Axis.2,
           color = ifelse(kruskal_filter <= 0.05, "DA", "NS"),
           alpha = Overall)) + 
  geom_point(shape = 20) +
  scale_color_manual(values = c("#FF82AB", "#666666")) +
  ggrepel::geom_text_repel(
    aes(Axis.1, Axis.2,
        label = ifelse(zscore >= 5, Genus, "")),
    show.legend = F,
    max.overlaps = Inf,
    color = "#333333",
    force = 2,
    size = 7.5) +
  ggnewscale::new_scale_colour() +
  scale_color_manual(values = pal_type) +
  geom_point(data = samples,
             mapping = aes(x = Axis.1,
                           y = Axis.2,
                           shape = site,
                           color = type),
             size = 9,
             alpha = 2/3,
             inherit.aes = F) +
  scale_shape_manual(values = c(18, 19)) +
  labs(x = paste("Axis.1:", round(res_pcoa$values[1, 2], 2)),
       y = paste("Axis.2:", round(res_pcoa$values[2, 2], 2)))


# setting the order of the labels right

Genus_order = paste(tax_tibble(phytted) %>% 
                      arrange(Phylum) %>% pull(Genus),
                    tax_tibble(phytted) %>% 
                      arrange(Phylum) %>% pull(FeatureID))

# extract the relative abundances for the ASVs

da_abundances =
  species %>% 
  filter(zscore >= 5) %>% 
  left_join(phytted %>% 
              transform("compositional") %>% 
              speedyseq::psmelt() %>% 
              group_by(site, type, OTU) %>% 
              summarise(relAbu = mean(Abundance)) %>% 
              rename(ASVs = OTU)) %>%
  arrange(Phylum) %>% 
  mutate(Genus_2 = paste(Genus, ASVs) %>% factor(levels = Genus_order, 
                                                 ordered = T))


write_tsv(da_abundances %>% 
  pivot_wider(
    values_from = relAbu,
    names_from = c("site", "type"),
    names_sep = "_"),
  "OUT/DA_rab.tsv"
)

# heatmap

b = 
  ggplot(da_abundances,
           aes(y = Genus_2 %>% gsub(".{31}$", "", .),
               x = site)) +
  geom_point(
    aes(size = relAbu,
        fill = type),
    alpha = 2/3,
    show.legend = T,
    shape = 21,
    stroke = 0.67,
    color = "#333333") +
  scale_fill_manual(values = pal_type) +
  scale_size(range = c(0, 20), labels = scales::percent_format()) +
  facet_grid(cols = vars(type),
             rows = vars(Phylum),
             scales = "free",
             space = "free") +
  theme(legend.position = "bottom",
        strip.text = element_blank(),
        axis.title.y = element_blank(),
        panel.border = element_blank(),
        legend.text = element_text(family = "Arial", size = 16),
        panel.spacing.y = unit(0.5, "lines"),
        panel.spacing.x = unit(0, "lines")) +
  guides(fill = "none")


c =
  ggplot(da_abundances,
         aes(y = Genus_2 %>% gsub(".{31}$", "", .),
             x = 1)) +
  geom_tile(aes(fill = Phylum),
            width = 0.5) +
  hrbrthemes::scale_fill_ipsum() +
  theme_void() +
  facet_grid(rows = vars(Phylum), scales = "free", space = "free") +
  theme(legend.position = "right",
        legend.box = "horizontal",
        legend.text = element_text(family = "Arial", size = 16),
        legend.title = element_blank(),
        strip.text = element_blank(),
        panel.spacing.y = unit(0.5, "lines")) +
  guides(fill = guide_legend(ncol = 1))

bc = b + c +
  plot_layout(widths = c(4, 1))


tiff("OUT/fig4a.tiff", 360, 360, "mm", res = 1000)
a
dev.off()

tiff("OUT/fig5c.tiff", 360, 360, "mm", res = 1000)
bc
dev.off()

cairo_pdf("OUT/fig5c.pdf", 1600, 1600)
bc
dev.off()


```

### CYANO ONLY

```{r}

# phylogenetic tree of the cyanobacteria in the metataxonomy

cyano = 
  subset_taxa(phytted, Phylum == "Cyanobacteria") %>% 
  transform("compositional") %>% 
  prune_taxa(taxa_sums(.) > 0, .) %>% 
  speedyseq::tree_glom(resolution = .03) %>% 
  metagMisc::phyloseq_add_max_tax_rank()

cyano@tax_table$Genus =
  coalesce(tax_table(cyano)$Genus, 
           paste(tax_table(cyano) %>% select(Rank), 
                 tax_table(cyano)$LowestTaxon,
                 sep = "_"))

p_load(ggtree)


p = ggtree(cyano)

nodes = 
  p$data %>% 
  select(parent, Genus) %>% 
  na.omit() %>% 
  distinct()

treplot = 
  ggtree(cyano,
         layout = "rectangular",
         size = .5) + 
  geom_tiplab(aes(label = coalesce(Genus, paste(RankName, LowestTaxon, sep = "_"))),
              size = 5,
              linetype = 3,
              align = T,
              offset = .02) +
  geom_tippoint(
    aes(color = type,
        shape = site,
        size = Abundance),
    position = position_dodge(width = 0.02),
    alpha = .9) +
  scale_color_manual(values = pal_type) +
  scale_shape_manual(values = c(15, 19)) +
  scale_size(range = c(3, 6)) +
  xlim(0, 0.55)

cairo_pdf("OUT/cyano_tree_silva_OTU.pdf",
          width = 20, height = 20)
treplot
dev.off()


```