author: ""
date: '`r format(Sys.time(), "%d %m %Y")`'
header-includes:
     - \usepackage{color, fancyvrb}
output:
    rmdformats::readthedown:
        highlight: kate
        number_sections : yes
    pdf_document:
        toc: yes
        toc_depth: 2
        number_sections : yes
---

```{r load-packages, include=FALSE}

#install.packages(c("picante", "rmdformats"))
#mamba install -c conda-forge freetype libpng harfbuzz fribidi
#mamba install -c conda-forge r-systemfonts r-svglite r-kableExtra freetype fontconfig harfbuzz fribidi libpng
library(knitr)
library(rmdformats)
library(readxl)
library(dplyr)
library(kableExtra)
library(openxlsx)
library(DESeq2)
library(writexl)

options(max.print="75")
knitr::opts_chunk$set(fig.width=8,
                                            fig.height=6,
                                            eval=TRUE,
                                            cache=TRUE,
                                            echo=TRUE,
                                            prompt=FALSE,
                                            tidy=FALSE,
                                            comment=NA,
                                            message=FALSE,
                                            warning=FALSE)
opts_knit$set(width=85)
# Phyloseq R library
#* Phyloseq web site : https://joey711.github.io/phyloseq/index.html
#* See in particular tutorials for
#    - importing data: https://joey711.github.io/phyloseq/import-data.html
#    - heat maps: https://joey711.github.io/phyloseq/plot_heatmap-examples.html
#rmarkdown::render('Phyloseq.Rmd',output_file='Phyloseq.html')
#options(max.print = 1e6)
```

# Data

Import raw data and assign sample key:

```{r, echo=FALSE, warning=FALSE}
#extend qiime2_metadata_for_qza_to_phyloseq.tsv with Diet and Flora
#setwd("~/DATA/Data_Laura_16S_2/core_diversity_e4753")
#map_corrected <- read.csv("qiime2_metadata_for_qza_to_phyloseq.tsv", sep="\t", row.names=1)
#knitr::kable(map_corrected) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
```

# Prerequisites to be installed

* R : https://pbil.univ-lyon1.fr/CRAN/
* R studio : https://www.rstudio.com/products/rstudio/download/#download

```R
install.packages("dplyr")     # To manipulate dataframes
install.packages("readxl")    # To read Excel files into R
install.packages("ggplot2")   # for high quality graphics
install.packages("heatmaply")
source("https://bioconductor.org/biocLite.R")
biocLite("phyloseq")
```

```{r libraries, echo=TRUE, message=FALSE}
#mamba install -c conda-forge r-ggplot2 r-vegan r-data.table
#BiocManager::install("microbiome")
#install.packages("ggpubr")
#install.packages("heatmaply")
library("readxl") # necessary to import the data from Excel file
library("ggplot2") # graphics
library("picante")
library("microbiome") # data analysis and visualisation
library("phyloseq") # also the basis of data object. Data analysis and visualisation
library("ggpubr") # publication quality figures, based on ggplot2
library("dplyr") # data handling, filter and reformat data frames
library("RColorBrewer") # nice color options
library("heatmaply")
library(vegan)
library(gplots)
#install.packages("openxlsx")
library(openxlsx)
```

# Read the data and create phyloseq objects

Three tables are needed

* OTU
* Taxonomy
* Samples

## Create analysis-specific phyloseq objects

We maintain **one filtered “base” phyloseq object** and then derive **analysis-specific** objects from it.
This avoids accidental overwriting and ensures each analysis uses the appropriate data scale (counts vs. relative abundance vs. rarefied counts).

- **`ps_raw`** → Raw imported phyloseq object (**integer counts**; import stage only)
- **`ps_base`** → `ps_raw` with **taxonomy + sample metadata** properly aligned (the clean master object before any filtering)
- **`ps_pruned`** → Optional **sample subset** of `ps_base` (e.g., drop unwanted samples by ID/pattern); still **integer counts**
- **`ps_filt`** → The shared filtered backbone: **low-depth samples removed** + taxa with **zero total counts dropped**; remains **integer counts**

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

        library(tidyr)

        # For QIIME1
        #ps.ng.tax <- import_biom("./exported_table/feature-table.biom", "./exported-tree/tree.nwk")

        # For QIIME2
        #install.packages("remotes")
        #remotes::install_github("jbisanz/qiime2R")
        #"core_metrics_results/rarefied_table.qza", rarefying performed in the code, therefore import the raw table.
        library(qiime2R)
        ps_raw <- qza_to_phyloseq(
            features =  "dada2_tests/test_59_f235_r245/table.qza",
            tree = "rooted-tree.qza",
            metadata = "qiime2_metadata_for_qza_to_phyloseq.tsv"
        )

        # Refresh/ensure sample_data (optional but keeps things explicit)
        sample_df <- read.csv("./qiime2_metadata_for_qza_to_phyloseq.tsv", sep="\t", row.names=1, check.names=FALSE)
        SAM <- sample_data(sample_df, errorIfNULL = TRUE)

        # Add taxonomy table (exported from QIIME2)
        taxonomy <- read.delim("./exported-taxonomy/taxonomy.tsv", sep="\t", header=TRUE)

        # Separate taxonomy string into separate ranks
        taxonomy_df <- taxonomy %>% separate(Taxon, into = c("Domain","Phylum","Class","Order","Family","Genus","Species"), sep = ";", fill = "right", extra = "drop")
        # Use Feature.ID as rownames
        rownames(taxonomy_df) <- taxonomy_df$Feature.ID
        taxonomy_df <- taxonomy_df[, -c(1, ncol(taxonomy_df))]  # Drop Feature.ID and Confidence
        # Create tax_table
        tax_table_final <- phyloseq::tax_table(as.matrix(taxonomy_df))

        # Base object: raw integer counts + metadata + taxonomy
        ps_base <- merge_phyloseq(ps_raw, SAM, tax_table_final)
        print(ps_base)

        #colnames(phyloseq::tax_table(ps_base)) <- c("Domain","Phylum","Class","Order","Family","Genus","Species")
        saveRDS(ps_base, "./ps_base.rds")
```

Visualize data
```{r, echo=TRUE, warning=FALSE}
# Inspect the base object (raw integer counts)
sample_names(ps_base)
rank_names(ps_base)
sample_variables(ps_base)

# Optional: prune to a naming pattern (avoids hard-coding long sample lists)
#samples_keep <- sample_names(ps_base)
#samples_keep <- samples_keep[grepl("^sample-[A-H]", samples_keep)]
samples_keep <- c("O23092004","O24010402","O23083101","A23080101","PC_8","O23092702","O24010401","O23092001","O23082402","U25011701","O23120702","O23121301","O23112304","U24111801","PC01","A24080201","A23060602","A23051102","UR009768","O24031305","O23090801","UR009909","U24121801","O23120103","U23090801","O23121501","O23110901","A23051103","O24013102","O24011801","O23091403","O23102704","O24020903","A24030402","O23100401","U24101801","O24011105","O24010302","O23121502","O23092005","O24021402","O23102602","NTC_3","O23111502","PC_7","A24071701","O23111501","O24011005","A23112001","U24080201","O24010404","O23112205","O23092803","NTC_6","NTC_8","O24011001","O24010301","O24013104","O24020103","O24011201","O23102601","NTC_2","A23112002","O23102502","O24020901","O24021502","O24021503","O23112901","O23112303","O23112201","O23112902","O24020703","O24010403","O23092003","U23052401","O24021501","O23120701","NTC_11","O24022101","O24011202","O23110902","O23092701","O24011002","A24031201","U23052201","O24021401","O23091303","O24021601","NTC_13","O23092104","O24011102","O23092102","NTC_10","O24022302","A24072901","O24022301","O24022202","O24022901","A24062801","O24011004","O23090602","O23112204","O23100701","O24012403","O24011003","A24071901","O24020802","O24021403","O24021602","NTC_14","NTC_9","PC_1","O24011901","O24013103","O24021404","O24030601","O24020701","O24011204","O24011103","A24030401","NTC_4","O24020902","O23092201","NTC_7","U23071901","O24020904","O23110903","U25011702","U23071701","O24031304","O24030602","O23091501","O24012401","O24030603","O24020704","PC_4","O23110904","O24020804","O23100801","O24030604","O23100802","NTC_15","O24022201","O24012501","O23100706","O23090601","O23100803","O24031301","O23082303","O24010501","O24022902","O23102703","O23110203","O23092202","NTC_12","O23102701","O23082302","O23121304","NTC_16","O24020803","O23090802","O24020101","U25020701","O23121306","O23082301","O23092203","O23110302","O23120101","O24020905","PC_5","O23092802","U25022101","O23100702","PC_6","O23110301","O23081801","O24020705","O24011207","O23091302","O23083001","O23102503","O23121303","O24011206","O24011802","O23102504","O23110204","O24022102","O23120104","O23092002","O23102603","O23092101","O23100503","O23092103","O23111601","O24012402","O24013105","PC_3","O23092801","O23091301","O23100402","O23121503","O24011104","O23112206","O23100703","O23102702","O24022801","O24011203","O23091401","O23090803","O23100705","O23102604","O23091402","PC_2","O23121302","U24070401","O23092006","O23100704","O23112203","O24020702","O24020102","O23112301","O23100804","O24020805","O23121305","O23112202","O24011205","O24013101","O24020801","O23110202","O24031302","O23111602","O23102501","O23110101","O23090701","O24011101","O23091305","O23083102","O24021603","O23112302")
ps_pruned <- prune_samples(samples_keep, ps_base)

# Drop taxa absent after pruning
ps_pruned <- prune_taxa(taxa_sums(ps_pruned) > 0, ps_pruned)

# Quick sanity checks
nsamples(ps_base); ntaxa(ps_base)
nsamples(ps_pruned); ntaxa(ps_pruned)
```

Preprocessing statistics for each sample
```{r, echo=TRUE, warning=FALSE}
 denoising_stats <- read.csv("dada2_tests/test_59_f235_r245/data_stats.tsv", sep="\t")
 # Display the table
 kable(denoising_stats, caption = "Preprocessing statistics for each sample") %>%
     kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
```

In the QC filtering step, we removed 14 samples that either fell below the minimum sequencing depth (library size < 12,201 reads) or were pre-specified for exclusion via sample_ids (as defined in the previous step). The filtered object (ps_filt) therefore contains only samples meeting the depth cutoff, and taxa were re-pruned to retain only those with nonzero total abundance across the retained samples.

```{r, echo=TRUE, warning=FALSE}
# ------------------------------------------------------------
#   Filter low-depth samples (recommended for all analyses)
# ------------------------------------------------------------
min_depth <- 12201  # <-- adjust to your data / study design, keeps all!
ps_filt <- prune_samples(sample_sums(ps_pruned) >= min_depth, ps_pruned)
ps_filt <- prune_taxa(taxa_sums(ps_filt) > 0, ps_filt)
ps_filt

# Keep a depth summary for reporting / QC
depth_summary <- summary(sample_sums(ps_filt))
depth_summary
```

**Differential abundance (DESeq2)** → **`ps_deseq`**: **non-rarefied integer counts** derived from `ps_filt`, with optional **count-based** taxon prefilter
    *(default: taxa total counts ≥ 10 across all samples)*

From `ps_filt` (e.g. 5669 taxa and 239 samples), we branch into analysis-ready objects in two directions:

* Direction 1 for diversity analyses

    - **Alpha diversity**: `ps_rarefied` ✅ (common)
    - **Beta diversity**:
        - **Unweighted UniFrac / Jaccard**: `ps_rarefied` ✅ (often recommended)
        - **Bray–Curtis / ordination on abundances**: `ps_rel` or Hellinger ✅ (rarefaction optional)
        - **Aitchison (CLR)**: CLR-transformed (non-rarefied) ✅ (no rarefaction)

Rarefaction

```{r, echo=TRUE, warning=FALSE}

# RAREFACTION
set.seed(9242)  # This will help in reproducing the filtering and nomalisation.
ps_rarefied <- rarefy_even_depth(ps_filt, sample.size = 12201)

# # NORMALIZE number of reads in each sample using median sequencing depth.
# total = median(sample_sums(ps.ng.tax))
# #> total
# #[1] 42369
# standf = function(x, t=total) round(t * (x / sum(x)))
# ps.ng.tax = transform_sample_counts(ps.ng.tax, standf)
# ps_rel <- microbiome::transform(ps.ng.tax, "compositional")
#
# saveRDS(ps.ng.tax, "./ps.ng.tax.rds")
```

* Direction 2 for taxonomic composition plots

    - **Taxonomic composition** → **`ps_rel`**: **relative abundance** (compositional) computed **after sample filtering** (e.g. 5669 taxa and 239 samples)
    - **Optional cleaner composition plots** → **`ps_abund` / `ps_abund_rel`**: taxa filtered for plotting (e.g., keep taxa with **mean relative abundance > 0.1%**); (e.g. 95 taxa and 239 samples)
    `ps_abund` = **counts**, `ps_abund_rel` = **relative abundance** *(use for visualization, not DESeq2)*

For the heatmaps, we focus on the most abundant OTUs by first converting counts to relative abundances within each sample. We then filter to retain only OTUs whose mean relative abundance across all samples exceeds 0.1% (0.001). We are left with 199 OTUs which makes the reading much more easy.

```{r, echo=FALSE, warning=FALSE}
# 1) Convert to relative abundances
ps_rel <- transform_sample_counts(ps_filt, function(x) x / sum(x))

# 2) Get the logical vector of which OTUs to keep (based on relative abundance)
keep_vector <- phyloseq::filter_taxa(
    ps_rel,
    function(x) mean(x) > 0.001,
    prune = FALSE
)

# 3) Use the TRUE/FALSE vector to subset absolute abundance data
ps_abund <- prune_taxa(names(keep_vector)[keep_vector], ps_filt)

# 4) Normalize the final subset to relative abundances per sample
ps_abund_rel <- transform_sample_counts(
    ps_abund,
    function(x) x / sum(x)
)
```

# Heatmaps

```{r, echo=FALSE, warning=FALSE}
datamat_ = as.data.frame(otu_table(ps_abund))
#datamat <- datamat_[c("1","2","5","6","7",  "8","9","10","12","13","14",    "15","16","17","18","19","20",  "21","22","23","24","25","26","27","28",    "29","30","31","32",  "33","34","35","36","37","38","39","51",    "40","41","42","43","44","46",  "47","48","49","50","52","53","55")]
datamat <- datamat_[c("O23092004","O24010402","O23083101","A23080101","PC_8","O23092702","O24010401","O23092001","O23082402","U25011701","O23120702","O23121301","O23112304","U24111801","PC01","A24080201","A23060602","A23051102","UR009768","O24031305","O23090801","UR009909","U24121801","O23120103","U23090801","O23121501","O23110901","A23051103","O24013102","O24011801","O23091403","O23102704","O24020903","A24030402","O23100401","U24101801","O24011105","O24010302","O23121502","O23092005","O24021402","O23102602","NTC_3","O23111502","PC_7","A24071701","O23111501","O24011005","A23112001","U24080201","O24010404","O23112205","O23092803","NTC_6","NTC_8","O24011001","O24010301","O24013104","O24020103","O24011201","O23102601","NTC_2","A23112002","O23102502","O24020901","O24021502","O24021503","O23112901","O23112303","O23112201","O23112902","O24020703","O24010403","O23092003","U23052401","O24021501","O23120701","NTC_11","O24022101","O24011202","O23110902","O23092701","O24011002","A24031201","U23052201","O24021401","O23091303","O24021601","NTC_13","O23092104","O24011102","O23092102","NTC_10","O24022302","A24072901","O24022301","O24022202","O24022901","A24062801","O24011004","O23090602","O23112204","O23100701","O24012403","O24011003","A24071901","O24020802","O24021403","O24021602","NTC_14","NTC_9","PC_1","O24011901","O24013103","O24021404","O24030601","O24020701","O24011204","O24011103","A24030401","NTC_4","O24020902","O23092201","NTC_7","U23071901","O24020904","O23110903","U25011702","U23071701","O24031304","O24030602","O23091501","O24012401","O24030603","O24020704","PC_4","O23110904","O24020804","O23100801","O24030604","O23100802","NTC_15","O24022201","O24012501","O23100706","O23090601","O23100803","O24031301","O23082303","O24010501","O24022902","O23102703","O23110203","O23092202","NTC_12","O23102701","O23082302","O23121304","NTC_16","O24020803","O23090802","O24020101","U25020701","O23121306","O23082301","O23092203","O23110302","O23120101","O24020905","PC_5","O23092802","U25022101","O23100702","PC_6","O23110301","O23081801","O24020705","O24011207","O23091302","O23083001","O23102503","O23121303","O24011206","O24011802","O23102504","O23110204","O24022102","O23120104","O23092002","O23102603","O23092101","O23100503","O23092103","O23111601","O24012402","O24013105","PC_3","O23092801","O23091301","O23100402","O23121503","O24011104","O23112206","O23100703","O23102702","O24022801","O24011203","O23091401","O23090803","O23100705","O23102604","O23091402","PC_2","O23121302","U24070401","O23092006","O23100704","O23112203","O24020702","O24020102","O23112301","O23100804","O24020805","O23121305","O23112202","O24011205","O24013101","O24020801","O23110202","O24031302","O23111602","O23102501","O23110101","O23090701","O24011101","O23091305","O23083102","O24021603","O23112302")]
# Remove rows with zero variance
datamat <- datamat[apply(datamat, 1, var) > 0, ]
# Remove cols with zero variance
#datamat <- datamat[, apply(datamat, 2, var) > 0]

# (optional) replace with your actual column names
#colnames(datamat) <- c(
#    "A24040201",...
#  )

# ---------- 0) Sample names from datamat ----------
samplenames <- sub("^sample-", "", colnames(datamat))

# ---------- 1) Read metadata ----------
meta_path <- "qiime2_metadata.tsv"
meta <- read.delim(
    meta_path,
    sep = "\t",
    header = TRUE,
    stringsAsFactors = FALSE,
    check.names = FALSE,
    comment.char = ""
)

# ---------- 2) Identify SampleID + Group columns ----------
sample_id_col <- c("#SampleID","SampleID","sample-id","sampleid")
sample_id_col <- sample_id_col[sample_id_col %in% colnames(meta)][1]
if (is.na(sample_id_col)) sample_id_col <- colnames(meta)[1]

group_col <- c("Group","group","GROUP")
group_col <- group_col[group_col %in% colnames(meta)][1]
if (is.na(group_col)) stop("No 'Group' column found in metadata.")

# ---------- 3) Build lookup: sample -> group ----------
meta_ids <- sub("^sample-", "", meta[[sample_id_col]])
meta_groups <- trimws(as.character(meta[[group_col]]))
group_map <- setNames(meta_groups, meta_ids)

# Map datamat columns to group labels
groups <- unname(group_map[samplenames])
groups[is.na(groups)] <- "unknown"

# ---------- 4) Color mapping for YOUR labels ----------
# (Adjust colors if you prefer different ones)
color_map <- c(
    "1" = "#a6cee3",
    "2" = "#1f78b4",
    "3" = "#b2df8a",
    "4" = "#33a02c",
    "5" = "#fb9a99",
    "negative control" = "#6a3d9a",
    "positive control" = "#ff7f00",
    "unknown" = "GREY"
)

# Assign colors safely
sampleCols <- unname(color_map[groups])
sampleCols[is.na(sampleCols)] <- "GREY"
names(sampleCols) <- samplenames

# ---------- 5) Checks ----------
cat("Unique groups found in datamat:\n")
print(sort(unique(groups)))
cat("\nCounts per group:\n")
print(table(groups, useNA = "ifany"))
cat("\nFirst 10 sample colors:\n")
print(head(sampleCols, 10))

# Optional: list any samples that didn't match metadata
unmatched <- samplenames[groups == "unknown"]
if (length(unmatched) > 0) {
    cat("\nWARNING: Unmatched samples (showing up to 20):\n")
    print(head(unmatched, 20))
}
```

```{r, echo=TRUE, warning=FALSE, fig.cap="Heatmap", out.width = '100%', fig.align= "center"}
## --- 1) order columns by group (and then by sample name) ---
group_order <- c("1","2","3","4","5","negative control","positive control","unknown")
groups_fac  <- factor(groups, levels = group_order, ordered = TRUE)

col_order <- order(groups_fac, samplenames)

datamat_ord     <- datamat[, col_order, drop = FALSE]
groups_ord      <- groups[col_order]
samplenames_ord <- samplenames[col_order]
sampleCols_ord  <- sampleCols[col_order]

stopifnot(identical(colnames(datamat_ord), samplenames_ord))

## group separators
grp_counts <- table(factor(groups_ord, levels = group_order))
grp_breaks <- cumsum(as.integer(grp_counts[grp_counts > 0]))

## --- 2) cluster ROWS using the *ordered* matrix (columns don't matter, but be consistent) ---
hr   <- hclust(as.dist(1 - cor(t(datamat_ord), method = "pearson")), method = "complete")
mycl <- cutree(hr, h = max(hr$height) / 1.08)

mycol_palette <- c("YELLOW","DARKBLUE","DARKORANGE","DARKMAGENTA","DARKCYAN","DARKRED",
                                     "MAROON","DARKGREEN","LIGHTBLUE","PINK","MAGENTA","LIGHTCYAN",
                                     "LIGHTGREEN","BLUE","ORANGE","CYAN","RED","GREEN")
mycol <- mycol_palette[as.vector(mycl)]

## --- 3) plot using datamat_ord and sampleCols_ord; keep column order fixed ---
library(RColorBrewer)
heatmap_colors <- colorRampPalette(brewer.pal(9, "Blues"))(100)

png("figures/heatmap.png", width = 1800, height = 2400)
heatmap.2(
    as.matrix(datamat_ord),
    Rowv = as.dendrogram(hr),
    Colv = NA,                        # IMPORTANT: do NOT cluster columns
    dendrogram = "row",
    scale = "row",
    trace = "none",
    col = heatmap_colors,
    cexRow = 1.2,
    cexCol = 0.8,
    RowSideColors = mycol,
    ColSideColors = sampleCols_ord,   # IMPORTANT: use ordered colors
    srtCol = 85,
    labRow = row.names(datamat_ord),
    labCol = samplenames_ord,         # optional but explicit
    key = TRUE,
    margins = c(10, 15),
    lhei = c(0.7, 15),
    lwid = c(1, 8),
    colsep = grp_breaks,              # optional separators
    sepcolor = "black",
    sepwidth = c(0.02, 0.02)
)
dev.off()

knitr::include_graphics("./figures/heatmap.png")
```

\pagebreak

```{r, echo=FALSE, warning=FALSE}
    library(stringr)
#for id in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100  101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199; do
#for id in 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300; do
#for id in 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382; do
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Domain\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Domain\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Phylum\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Phylum\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Class\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Class\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Order\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Order\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Family\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Family\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Genus\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Genus\"], \"__\")[[1]][2]"
#  echo "phyloseq::tax_table(ps_abund_rel)[${id},\"Species\"] <- str_split(phyloseq::tax_table(ps_abund_rel)[${id},\"Species\"], \"__\")[[1]][2]"
#done

phyloseq::tax_table(ps_abund_rel)[1,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[1,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[1,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[2,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[2,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[3,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[3,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[4,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[4,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[5,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[5,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[6,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[6,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[7,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[7,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[8,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[8,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[9,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[9,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[10,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[10,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[11,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[11,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[12,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[12,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[13,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[13,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[14,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[14,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[15,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[15,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[16,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[16,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[17,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[17,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[18,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[18,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[19,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[19,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[20,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[20,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[21,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[21,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[22,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[22,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[23,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[23,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[24,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[24,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[25,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[25,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[26,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[26,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[27,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[27,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[28,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[28,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[29,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[29,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[30,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[30,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[31,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[31,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[32,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[32,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[33,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[33,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[34,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[34,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[35,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[35,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[36,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[36,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[37,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[37,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[38,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[38,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[39,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[39,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[40,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[40,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[41,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[41,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[42,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[42,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[43,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[43,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[44,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[44,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[45,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[45,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[46,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[46,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[47,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[47,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[48,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[48,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[49,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[49,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[50,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[50,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[51,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[51,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[52,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[52,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[53,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[53,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[54,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[54,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[55,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[55,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[56,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[56,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[57,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[57,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[58,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[58,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[59,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[59,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[60,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[60,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[61,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[61,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[62,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[62,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[63,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[63,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[64,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[64,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[65,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[65,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[66,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[66,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[67,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[67,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[68,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[68,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[69,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[69,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[70,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[70,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[71,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[71,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[72,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[72,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[73,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[73,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[74,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[74,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[75,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[75,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[76,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[76,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[77,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[77,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[78,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[78,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[79,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[79,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[80,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[80,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[81,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[81,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[82,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[82,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[83,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[83,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[84,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[84,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[85,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[85,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[86,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[86,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[87,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[87,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[88,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[88,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[89,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[89,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[90,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[90,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[91,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[91,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[92,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[92,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[93,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[93,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[94,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[94,"Species"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Domain"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Domain"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Phylum"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Phylum"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Class"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Class"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Order"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Order"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Family"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Family"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Genus"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Genus"], "__")[[1]][2]
phyloseq::tax_table(ps_abund_rel)[95,"Species"] <- str_split(phyloseq::tax_table(ps_abund_rel)[95,"Species"], "__")[[1]][2]
```

# Taxonomic summary

## Bar plots in phylum level

```{r, fig.width=16, fig.height=8, echo=TRUE, warning=FALSE}
    sample_order <- c(
            "U24080201","U25020701","O23092004","U24101801","U25022101","O23102703","A24062801","O23112205","U23071901","A23112002","U24111801","O23110101","U24121801","O23120101","O24011202","O23090803","A23060602","A24030402","U25011701","O24011901","U23090801","O24011201","O24011003","O23092202","O23082301","O23091403","O23112901","O23092201","O24013103","O24021403","O24010402","O23092005","O23092203","O24010302","O23090701","O23091501","O23092701","O24022202","O23092802","O23090601","O23100401","O24022102","O23081801","O23092006","O23100503","O23090602","O24013104","O24020803","O24010301","O24010404","O23090802","O23092801","O24022801","O23100706","O23102602","O24021601","O24012401","O24021603","O24022901","O24021501","O23110902","O23102601","O23102704","O23100803","O23102701","O24021401","O24022101","O24030603","O23110901","O23110903","O23110301","O24022301","O23102502","O23111501","O23111602","O24020705","O24021502","O24022201","O23110202","O23090801","U23071701","U23052401","U23052201","U24070401","O24011801","O23092003","A24071901","A24072901","O24011102","O23121501","O23092104","O23092001","O23121301","O24020701","O23112201","O23100701","O23100801","O24020903","O24020901","O24020703","O23112204","O23100802","O24011205","O23092002","O24011207","O23092103","O23102501","O24011005","A24030401","O24011004","A23051102","U25011702","O24011204","O23121502","O23120702","O24011206","O24021404","O23092101","O24010403","O23112303","O23083001","O23082302","O24010401","O24022302","O24010501","O23112902","O23082303","O23083102","O24013101","O23100402","O24020801","O23120701","O23121304","O24021602","O24011802","O23121306","O23120103","O24020905","O24012403","O24013102","O24021503","O24020904","O23102504","O24013105","O24030601","O23100705","O24030604","O23111601","O24020103","O24030602","O23110302","O23102603","O24031304","O24021402","O24020101","O24012501","O24020804","O23100804","O23102503","O24022902","O24020704","O23110904","O24020102","O24012402","O23102702","O23102604","O23110204","O23110203","O23083101","O23092702","O23112304","A23051103","A24071701","A23080101","A24031201","A24080201","O24011105","O23091305","O23121302","O23092803","O23091303","O23112301","O24011203","A23112001","O24011001","O24011002","O23091302","O24020902","O23092102","O23091401","O23121503","O23091402","O24020702","O23091301","O23112206","O24011103","O23121305","O24011101","O23121303","O24011104","O23120104","O23100704","O23112302","O23112203","O23100703","O24020805","O24020802","O23112202","O24031302","O23111502","O23100702","O24031301","O24031305","O23082402","NTC_2","NTC_3","NTC_16","NTC_4","NTC_6","NTC_7","NTC_8","NTC_9","NTC_11","NTC_12","NTC_13","NTC_14","NTC_15","NTC_10","UR009768","UR009909","PC01","PC_1","PC_2","PC_8","PC_4","PC_3","PC_5","PC_6","PC_7"
    )

    # create a sample ID column in sample_data (or overwrite an existing one)
    sample_data(ps_abund_rel)$SampleID <- sample_names(ps_abund_rel)

    # set the order as a factor with the desired levels
    sample_data(ps_abund_rel)$SampleID <- factor(
        sample_data(ps_abund_rel)$SampleID,
        levels = sample_order
    )

    #aes(color="Phylum", fill="Phylum") --> aes()
    #ggplot(data=data, aes(x=Sample, y=Abundance, fill=Phylum))
    my_colors <- c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")

    plot_bar(ps_abund_rel, x = "SampleID", fill = "Phylum") +
        geom_bar(stat = "identity", position = "stack") +
        scale_fill_manual(values = my_colors) +
        theme(
            axis.text.x = element_text(
                angle = 85,           # 85° rotation
                hjust = 1,            # right-justified so labels are neat
                vjust = 1,
                size  = 5,
                colour = "black"
            ),
            axis.text.y = element_text(size = 7, colour = "black"),
            legend.position = "bottom"
        ) +
        guides(fill = guide_legend(nrow = 2))
```

### Aggregate samples by group and normalize read counts within each group to correct for differences in sequencing depth.

```{r, echo=TRUE, warning=FALSE}
    # merge + normalize
    ps_abund_rel_group  <- merge_samples(ps_abund_rel, "Group")
    ps_abund_rel_group_ <- transform_sample_counts(
        ps_abund_rel_group,
        function(x) x / sum(x)
    )

    # desired order on x-axis
    group_order <- c("1", "2", "3", "4", "5", "negative control", "positive control")

    plot_bar(ps_abund_rel_group_, fill = "Phylum") +
        geom_bar(stat = "identity", position = "stack") +
        scale_x_discrete(limits = group_order) +          # <-- set order here
        scale_fill_manual(values = my_colors) +
        labs(x = "Group") +   # <- change x-axis label from "Sample" to "Group"
        theme(axis.text = element_text(angle = 0, size = 10, colour="black"), axis.text.x = element_text(angle = -15),hjust = 10,vjust = 2)
```

## Bar plots in class level

```{r, fig.width=16, fig.height=8, echo=TRUE, warning=FALSE}
    my_colors <- c("darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink", "khaki2", "firebrick", "brown1", "darkorange1", "cyan1", "royalblue4", "darksalmon", "darkblue","royalblue4", "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen", "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1", "brown1", "darkorange1", "cyan1", "darkgrey")
    plot_bar(ps_abund_rel, x = "SampleID", fill = "Class") +
        geom_bar(stat = "identity", position = "stack") +
        scale_fill_manual(values = my_colors) +
        theme(
            axis.text.x = element_text(
                angle = 85,
                hjust = 1,
                vjust = 1,
                size  = 5,
                colour = "black"
            ),
            axis.text.y = element_text(size = 7, colour = "black"),
            legend.position = "bottom"
        ) +
        guides(fill = guide_legend(nrow = 3))
```

### Aggregate samples by group and normalize read counts within each group to correct for differences in sequencing depth.

```{r, echo=TRUE, warning=FALSE}
    plot_bar(ps_abund_rel_group_, fill="Class") + geom_bar(aes(), stat="identity", position="stack") + scale_x_discrete(limits = group_order) +
    scale_fill_manual(values = my_colors) + labs(x = "Group") + theme(axis.text = element_text(angle = 0, size = 10, colour="black"), axis.text.x = element_text(angle = -15),hjust = 10,vjust = 2)
```
\pagebreak

## Bar plots in order level

```{r, fig.width=16, fig.height=8, echo=TRUE, warning=FALSE}
    my_colors <- c(
        "darkblue", "darkgoldenrod1", "darkseagreen", "darkorchid",
        "darkolivegreen1", "lightskyblue", "darkgreen", "deeppink",
        "khaki2", "firebrick", "brown1", "darkorange1", "cyan1",
        "royalblue4", "darksalmon", "darkblue","royalblue4",
        "dodgerblue3", "steelblue1", "lightskyblue", "darkseagreen",
        "darkgoldenrod1", "darkseagreen", "darkorchid", "darkolivegreen1",
        "brown1", "darkorange1", "cyan1", "darkgrey"
    )

    plot_bar(ps_abund_rel, x = "SampleID", fill = "Order") +
        geom_bar(stat = "identity", position = "stack") +
        scale_fill_manual(values = my_colors) +
        theme(
            axis.text.x = element_text(
                angle = 85,
                hjust = 1,
                vjust = 1,
                size  = 5,
                colour = "black"
            ),
            axis.text.y = element_text(size = 7, colour = "black"),
            legend.position = "bottom"
        ) +
        guides(fill = guide_legend(nrow = 4))
```

### Aggregate samples by group and normalize read counts within each group to correct for differences in sequencing depth.

```{r, echo=TRUE, warning=FALSE}
    plot_bar(ps_abund_rel_group_, fill="Order") + geom_bar(aes(), stat="identity", position="stack") + scale_x_discrete(limits = group_order) +
    scale_fill_manual(values = my_colors) + labs(x = "Group") + theme(axis.text = element_text(angle = 0, size = 10, colour="black"), axis.text.x = element_text(angle = -15),hjust = 10,vjust = 2)
```
\pagebreak

## Bar plots in family level

```{r, fig.width=16, fig.height=8, echo=TRUE, warning=FALSE}
    my_colors <- c(
                    "#FF0000", "#000000", "#0000FF", "#C0C0C0", "#FFFFFF", "#FFFF00", "#00FFFF", "#FFA500", "#00FF00", "#808080", "#FF00FF", "#800080", "#FDD017", "#0000A0", "#3BB9FF", "#008000", "#800000", "#ADD8E6", "#F778A1", "#800517", "#736F6E", "#F52887", "#C11B17", "#5CB3FF", "#A52A2A", "#FF8040", "#2B60DE", "#736AFF", "#1589FF", "#98AFC7", "#8D38C9", "#307D7E", "#F6358A", "#151B54", "#6D7B8D", "#FDEEF4", "#FF0080", "#F88017", "#2554C7", "#FFF8C6", "#D4A017", "#306EFF", "#151B8D", "#9E7BFF", "#EAC117", "#E0FFFF", "#15317E", "#6C2DC7", "#FBB917", "#FCDFFF", "#15317E", "#254117", "#FAAFBE", "#357EC7"
                )
    plot_bar(ps_abund_rel, x = "SampleID", fill = "Family") +
        geom_bar(stat = "identity", position = "stack") +
        scale_fill_manual(values = my_colors) +
        theme(
            axis.text.x = element_text(
                angle = 85,
                hjust = 1,
                vjust = 1,
                size  = 5,
                colour = "black"
            ),
            axis.text.y = element_text(size = 7, colour = "black"),
            legend.position = "bottom"
        ) +
        guides(fill = guide_legend(nrow = 8))
```

### Aggregate samples by group and normalize read counts within each group to correct for differences in sequencing depth.

```{r, echo=TRUE, warning=FALSE}
    plot_bar(ps_abund_rel_group_, fill="Family") + geom_bar(aes(), stat="identity", position="stack") + scale_x_discrete(limits = group_order) +
    scale_fill_manual(values = my_colors) + labs(x = "Group") + theme(axis.text = element_text(angle = 0, size = 10, colour="black"), axis.text.x = element_text(angle = -15),hjust = 10,vjust = 2)
```
\pagebreak

```{r, fig.width=16, fig.height=8, echo=FALSE, warning=FALSE}
# !!!!NOT_USED!!!!: #Export Relative abundances of Phylum, Class, Order, and Family levels across all samples in Excel files!
library(phyloseq)
library(writexl)
library(dplyr)

# Function to check for NA or empty values in a taxonomic rank
check_taxa_names <- function(tax_table, rank) {
    tax_values <- tax_table[[rank]]
    na_count <- sum(is.na(tax_values) | tax_values == "")
    cat("Number of NA or empty values in", rank, ":", na_count, "\n")
    if (na_count > 0) {
        cat("Taxa with NA or empty", rank, ":\n")
        print(tax_values[is.na(tax_values) | tax_values == ""])
    }
}

# Function to create and save relative abundance table for a given taxonomic rank with normalization
save_taxa_abundance <- function(ps, rank, output_file) {
    # Check for NA or empty values in the taxonomy table
    tax_table_df <- as.data.frame(tax_table(ps))
    check_taxa_names(tax_table_df, rank)

    # Aggregate OTUs by taxonomic rank, removing taxa with NA at the specified rank
    ps_glom <- tax_glom(ps, taxrank = rank, NArm = TRUE)

    # Extract OTU table (relative abundances)
    otu_table <- as.data.frame(otu_table(ps_glom))

    # Normalize each column to sum to 1
    otu_table_normalized <- apply(otu_table, 2, function(x) x / sum(x))

    # Convert matrix to data frame
    otu_table_normalized <- as.data.frame(otu_table_normalized)

    # Verify column sums are approximately 1.0
    col_sums <- colSums(otu_table_normalized)
    if (any(abs(col_sums - 1) > 1e-6)) {
        warning("Column sums in ", rank, " table do not equal 1.0: ", paste(col_sums, collapse = ", "))
    } else {
        cat("Column sums for ", rank, " table are all approximately 1.0\n")
    }

    # Extract taxonomy table and get the specified rank for taxa names
    tax_table_glom <- as.data.frame(tax_table(ps_glom))
    taxa_names <- tax_table_glom[[rank]]

    # Replace NA or empty strings with "Unclassified"
    taxa_names <- ifelse(is.na(taxa_names) | taxa_names == "", paste0("Unclassified_", rank), taxa_names)

    # Ensure unique row names
    taxa_names <- make.unique(taxa_names)

    # Set row names to taxa names (for internal reference)
    rownames(otu_table_normalized) <- taxa_names

    # Add taxa names as a column
    otu_table_normalized[[rank]] <- taxa_names

    # Reorder to move rank column to the first position
    otu_table_normalized <- otu_table_normalized[, c(rank, setdiff(names(otu_table_normalized), rank))]

    # Rename sample columns by removing "sample-" prefix
    names(otu_table_normalized)[-1] <- sub("sample-", "", names(otu_table_normalized)[-1])

    # Write the data frame to Excel, including the rank column
    write_xlsx(otu_table_normalized, path = output_file)
    cat("Saved", output_file, "\n")
}

# Verify column sums of ps_abund_rel
col_sums <- colSums(otu_table(ps_abund_rel))
cat("Column sums of ps_abund_rel:\n")
summary(col_sums)

## Generate Excel files for Phylum, Class, Order, and Family levels with normalization and renamed sample names
#save_taxa_abundance(ps_abund_rel, "Phylum", "relative_abundance_phylum_old.xlsx")
#save_taxa_abundance(ps_abund_rel, "Class", "relative_abundance_class_old.xlsx")
#save_taxa_abundance(ps_abund_rel, "Order", "relative_abundance_order_old.xlsx")
#save_taxa_abundance(ps_abund_rel, "Family", "relative_abundance_family_old.xlsx")
```

```{r, fig.width=16, fig.height=8, echo=FALSE, warning=FALSE}
library(phyloseq)
library(writexl)
library(dplyr)

# Function to check for NA or empty values in a taxonomic rank
check_taxa_names <- function(tax_table, rank) {
    tax_values <- tax_table[[rank]]
    na_count <- sum(is.na(tax_values) | tax_values == "")
    cat("Number of NA or empty values in", rank, ":", na_count, "\n")
    if (na_count > 0) {
        cat("Taxa with NA or empty", rank, ":\n")
        print(tax_values[is.na(tax_values) | tax_values == ""])
    }
}

# Function to create and save relative abundance table for a given taxonomic rank with normalization
save_taxa_abundance <- function(ps, rank, output_file) {
    # Clean the taxonomy table by removing D_[level]__ prefixes
    tax_table_df <- as.data.frame(tax_table(ps))
    tax_table_df[[rank]] <- ifelse(is.na(tax_table_df[[rank]]) | tax_table_df[[rank]] == "",
                                                                 paste0("Unclassified_", rank),
                                                                 sub("^D_[0-9]+__(.+)", "\\1", tax_table_df[[rank]]))
    tax_table(ps) <- as.matrix(tax_table_df)  # Update taxonomy table with cleaned names

    # Check for NA or empty values in the taxonomy table
    check_taxa_names(tax_table_df, rank)

    # Aggregate OTUs by taxonomic rank, removing taxa with NA at the specified rank
    ps_glom <- tax_glom(ps, taxrank = rank, NArm = TRUE)

    # Extract OTU table (relative abundances)
    otu_table <- as.data.frame(otu_table(ps_glom))

    # Normalize each column to sum to 1
    otu_table_normalized <- apply(otu_table, 2, function(x) x / sum(x))

    # Convert matrix to data frame
    otu_table_normalized <- as.data.frame(otu_table_normalized)

    # Verify column sums are approximately 1.0
    col_sums <- colSums(otu_table_normalized)
    if (any(abs(col_sums - 1) > 1e-6)) {
        warning("Column sums in ", rank, " table do not equal 1.0: ", paste(col_sums, collapse = ", "))
    } else {
        cat("Column sums for ", rank, " table are all approximately 1.0\n")
    }

    # Extract taxonomy table and get the specified rank for taxa names
    tax_table_glom <- as.data.frame(tax_table(ps_glom))
    taxa_names <- tax_table_glom[[rank]]

    # Ensure unique row names
    taxa_names <- make.unique(taxa_names)

    # Set row names to taxa names (for internal reference)
    rownames(otu_table_normalized) <- taxa_names

    # Add taxa names as a column
    otu_table_normalized[[rank]] <- taxa_names

    # Reorder to move rank column to the first position
    otu_table_normalized <- otu_table_normalized[, c(rank, setdiff(names(otu_table_normalized), rank))]

    # Rename sample columns by removing "sample-" prefix
    names(otu_table_normalized)[-1] <- sub("sample-", "", names(otu_table_normalized)[-1])

    # Write the data frame to Excel, including the rank column
    write_xlsx(otu_table_normalized, path = output_file)
    cat("Saved", output_file, "\n")
}

# Verify column sums of ps_abund_rel
col_sums <- colSums(otu_table(ps_abund_rel))
cat("Column sums of ps_abund_rel:\n")
summary(col_sums)

## Generate Excel files for Phylum, Class, Order, and Family levels with normalization and renamed sample names
#save_taxa_abundance(ps_abund_rel, "Phylum", "relative_abundance_phylum.xlsx")
#save_taxa_abundance(ps_abund_rel, "Class", "relative_abundance_class.xlsx")
#save_taxa_abundance(ps_abund_rel, "Order", "relative_abundance_order.xlsx")
#save_taxa_abundance(ps_abund_rel, "Family", "relative_abundance_family.xlsx")

#Sum up the last two colums with the same row.names to a new column, export the file as csv, then delete the two rows before last, then merge them with csv2xls to a Excel-file, adapt the sheet-names.
#~/Tools/csv2xls-0.4/csv_to_xls.py relative_abundance_phylum.csv relative_abundance_order.csv relative_abundance_family.csv -d$'\t' -o relative_abundance_phylum_order_family.xls;
```

\pagebreak

# Alpha diversity
Plot Chao1 richness estimator, Observed OTUs, Shannon index, and Phylogenetic diversity.
Regroup together samples from the same group.
```{r, echo=FALSE, warning=FALSE}
# using rarefied data
#gunzip table_even4753.biom.gz
#alpha_diversity.py -i table_even42369.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering/rep_set.tre
#gunzip table_even4753.biom.gz
#alpha_diversity.py -i table_even4753.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering_stool/rep_set.tre
#gunzip table_even4753.biom.gz
#alpha_diversity.py -i table_even4753.biom --metrics chao1,observed_otus,shannon,PD_whole_tree -o adiv_even.txt -t ../clustering_swab/rep_set.tre
```

```{r, echo=TRUE, warning=FALSE}
hmp.meta <- meta(ps_rarefied)
hmp.meta$sam_name <- rownames(hmp.meta)

# ---- enforce Group order (edit if you have different labels) ----
group_order <- c("1","2","3","4","5","negative control", "positive control")
hmp.meta$Group <- factor(as.character(hmp.meta$Group), levels = group_order)

# for QIIME2: Lesen der Metriken
shannon <- read.table("exported_alpha/shannon/alpha-diversity.tsv", header=TRUE, sep="\t")  #cp -r ../Data_Karoline_16S_2025/exported_alpha/ .
faith_pd <- read.table("exported_alpha/faith_pd/alpha-diversity.tsv", header=TRUE, sep="\t")
observed <- read.table("exported_alpha/observed_features/alpha-diversity.tsv", header=TRUE, sep="\t")
#chao1 <- read.table("exported_alpha/chao1/alpha-diversity.tsv", header=TRUE, sep="\t")    #TODO: Check the correctness of chao1-calculation.

# Umbenennen für Klarheit
colnames(shannon) <- c("sam_name", "shannon")
colnames(faith_pd) <- c("sam_name", "PD_whole_tree")
colnames(observed) <- c("sam_name", "observed_otus")
#colnames(chao1) <- c("sam_name", "chao1")

# Merge alles in ein DataFrame
div.df <- Reduce(function(x, y) merge(x, y, by="sam_name"),
                                    list(shannon, faith_pd, observed))

# Meta-Daten einfügen
div.df <- merge(div.df, hmp.meta, by="sam_name")

# Reformat
div.df2 <- div.df[, c("sam_name", "Group", "shannon", "observed_otus", "PD_whole_tree")]
colnames(div.df2) <- c("Sample name", "Group", "Shannon", "OTU", "Phylogenetic Diversity")
write.csv(div.df2, file="alpha_diversities.txt")
knitr::kable(div.df2) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

#https://uc-r.github.io/t_test
#We can perform the test with t.test and transform our data and we can also perform the nonparametric test with the wilcox.test function.
stat.test.Shannon <- compare_means(
 Shannon ~ Group, data = div.df2,
 method = "t.test"
)
knitr::kable(stat.test.Shannon) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

div_df_melt <- reshape2::melt(div.df2)
#head(div_df_melt)

#https://plot.ly/r/box-plots/#horizontal-boxplot
#http://www.sthda.com/english/wiki/print.php?id=177
#https://rpkgs.datanovia.com/ggpubr/reference/as_ggplot.html
#http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/82-ggplot2-easy-way-to-change-graphical-parameters/
#https://plot.ly/r/box-plots/#horizontal-boxplot
#library("gridExtra")
#par(mfrow=c(4,1))
p <- ggboxplot(div_df_melt, x = "Group", y = "value",
                            facet.by = "variable",
                            scales = "free",
                            width = 0.5,
                            fill = "gray", legend= "right")
#ggpar(p, xlab = FALSE, ylab = FALSE)
lev <- levels(factor(div_df_melt$Group)) # get the variables
#FITTING4: delete H47(1) in lev
#lev <- lev[-c(3)]
# make a pairwise list that we want to compare.
#my_stat_compare_means
#https://stackoverflow.com/questions/47839988/indicating-significance-with-ggplot2-in-a-boxplot-with-multiple-groups
L.pairs <- combn(seq_along(lev), 2, simplify = FALSE, FUN = function(i) lev[i]) #%>% filter(p.signif != "ns")
my_stat_compare_means  <- function (mapping = NULL, data = NULL, method = NULL, paired = FALSE,
        method.args = list(), ref.group = NULL, comparisons = NULL,
        hide.ns = FALSE, label.sep = ", ", label = NULL, label.x.npc = "left",
        label.y.npc = "top", label.x = NULL, label.y = NULL, tip.length = 0.03,
        symnum.args = list(), geom = "text", position = "identity",
        na.rm = FALSE, show.legend = NA, inherit.aes = TRUE, ...)
{
        if (!is.null(comparisons)) {
                method.info <- ggpubr:::.method_info(method)
                method <- method.info$method
                method.args <- ggpubr:::.add_item(method.args, paired = paired)
                if (method == "wilcox.test")
                        method.args$exact <- FALSE
                pms <- list(...)
                size <- ifelse(is.null(pms$size), 0.3, pms$size)
                color <- ifelse(is.null(pms$color), "black", pms$color)
                map_signif_level <- FALSE
                if (is.null(label))
                        label <- "p.format"
                if (ggpubr:::.is_p.signif_in_mapping(mapping) | (label %in% "p.signif")) {
                        if (ggpubr:::.is_empty(symnum.args)) {
                                map_signif_level <- c(`****` = 1e-04, `***` = 0.001,
                                    `**` = 0.01, `*` = 0.05, ns = 1)
                        } else {
                             map_signif_level <- symnum.args
                        }
                        if (hide.ns)
                                names(map_signif_level)[5] <- " "
                }
                step_increase <- ifelse(is.null(label.y), 0.12, 0)
                ggsignif::geom_signif(comparisons = comparisons, y_position = label.y,
                        test = method, test.args = method.args, step_increase = step_increase,
                        size = size, color = color, map_signif_level = map_signif_level,
                        tip_length = tip.length, data = data)
        } else {
                mapping <- ggpubr:::.update_mapping(mapping, label)
                layer(stat = StatCompareMeans, data = data, mapping = mapping,
                        geom = geom, position = position, show.legend = show.legend,
                        inherit.aes = inherit.aes, params = list(label.x.npc = label.x.npc,
                                label.y.npc = label.y.npc, label.x = label.x,
                                label.y = label.y, label.sep = label.sep, method = method,
                                method.args = method.args, paired = paired, ref.group = ref.group,
                                symnum.args = symnum.args, hide.ns = hide.ns,
                                na.rm = na.rm, ...))
        }
}

# Rotate the x-axis labels to 45 degrees and adjust their position
p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust=1, size=8))
p2 <- p +
stat_compare_means(
    method="t.test",
    comparisons = list(c("1", "2"), c("1", "3"), c("1", "4"), c("1", "5"), c("2", "3"), c("2", "4"), c("2", "5"), c("3", "4"), c("3", "5"), c("4", "5")),
    label = "p.signif",
    symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("****", "***", "**", "*", "ns"))
)
#comparisons = L.pairs,
#symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05), symbols = c("****", "***", "**", "*")),
#stat_pvalue_manual
print(p2)
#https://stackoverflow.com/questions/20500706/saving-multiple-ggplots-from-ls-into-one-and-separate-files-in-r
ggsave("./figures/alpha_diversity_Group.png", device="png", height = 10, width = 15)
ggsave("./figures/alpha_diversity_Group.svg", device="svg", height = 10, width = 15)
```

```{r, echo=FALSE, warning=FALSE, fig.cap="Alpha diversity", out.width = '100%', fig.align= "center"}
## MANUALLY selected alpha diversities unter host-env after 'cp alpha_diversities.txt selected_alpha_diversities.txt'
#knitr::include_graphics("./figures/alpha_diversity_Group.png")
#selected_alpha_diversities<-read.csv("selected_alpha_diversities.txt",sep="\t")
#knitr::kable(selected_alpha_diversities) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
```

```{r, echo=FALSE, warning=FALSE, out.width = '100%', fig.align= "center"}
#!!# Beta diversity (Bray-Curtis distance)
#!!## Group1 vs Group2

#fig.cap="Beta diversity",

#for QIIME1: file:///home/jhuang/DATA/Data_Marius_16S/core_diversity_e42369/bdiv_even42369_Group/weighted_unifrac_boxplots/Group_Stats.txt

# -- for QIIME2: MANUALLY filter permanova-pairwise.csv and save as permanova-pairwise_.csv
# #grep "Permutations" exported_beta_group/permanova-pairwise.csv > permanova-pairwise_.csv
# #grep "Group1,Group2" exported_beta_group/permanova-pairwise.csv >> permanova-pairwise_.csv
# #grep "Group3,Group4" exported_beta_group/permanova-pairwise.csv >> permanova-pairwise_.csv
# beta_diversity_group_stats<-read.csv("permanova-pairwise_.csv",sep=",")
# #beta_diversity_group_stats <- beta_diversity_group_stats[beta_diversity_group_stats$Group.1 == "Group1" & beta_diversity_group_stats$Group.2 == "Group2", ]
# #beta_diversity_group_stats <- beta_diversity_group_stats[beta_diversity_group_stats$Group.1 == "Group3" & beta_diversity_group_stats$Group.2 == "Group4", ]
# knitr::kable(beta_diversity_group_stats) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

#NOTE: Run this Phyloseq.Rmd, then run the code of MicrobiotaProcess.R to manually generate Comparison_of_Bray_Distances_Group1_vs_Group2.png and Comparison_of_Bray_Distances_Group3_vs_Group4.png, then run this Phyloseq.Rmd!

#knitr::include_graphics("./figures/Comparison_of_Bray_Distances_Group1_vs_Group2.png")

```

# Principal coordinates analysis (PCoA) based on Bray–Curtis dissimilarity

Global PERMANOVA on the weighted UniFrac distance matrix indicated a significant effect of Group on overall community composition (adonis2: R² = 0.1606, F = 7.397, p = 1×10⁻⁴; 9,999 permutations).

```{r, echo=FALSE, results='asis'}
# --- Global beta-diversity (PERMANOVA) ---
cat("```text\n")
cat(
"[PERMANOVA result]\n",
"The object contained internal attribute: PCoA ADONIS\n",
"Permutation test for adonis under reduced model\n",
"Permutation: free\n",
"Number of permutations: 9999\n\n",
"vegan::adonis2(formula = .formula, data = sampleda, permutations = permutations, method = distmethod)\n",
"          Df SumOfSqs      R2      F Pr(>F)\n",
"Model      6   11.446 0.16058 7.3971  1e-04 ***\n",
"Residual 232   59.829 0.83942\n",
"Total    238   71.274 1.00000\n",
"---\n",
"Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
sep = ""
)
cat("```\n")
```

Pairwise PERMANOVA tests were performed on Bray–Curtis distance matrices to compare bacterial community composition between all pairs of sample groups (metadata column Group). For each pairwise comparison, the distance matrix was subset to samples from the two groups only, and significance was assessed using vegan::adonis2 with 9,999 permutations. Resulting p-values were adjusted for multiple testing using both Benjamini–Hochberg (BH/FDR) and Bonferroni corrections.

```{r, echo=FALSE, warning=FALSE, out.width = '100%', fig.align= "center"}
#, and the full results were exported to figures/Bray_pairwise_PERMANOVA.csv
# --- Pairwise PERMANOVA results ---

Bray_pairwise_PERMANOVA <- read.csv("figures/Bray_pairwise_PERMANOVA.csv", sep = ",")
knitr::kable(Bray_pairwise_PERMANOVA, caption = "Pairwise PERMANOVA results (distance-based community differences among Group levels).") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# --- Ordination figures ---

#knitr::include_graphics("./PCoA.png")
knitr::include_graphics("./PCoA2.png")
knitr::include_graphics("./PCoA3.png")
```

# Differential abundance analysis

Differential abundance analysis aims to find the differences in the abundance of each taxa between two groups of samples, assigning a significance value to each comparison.

```{r, echo=FALSE, warning=FALSE}
# ------------------------------------------------------------
#  DESeq2: non-rarefied integer counts + optional taxon prefilter
# ------------------------------------------------------------
ps_deseq <- ps_filt

Group1<-c("U24080201","U25020701","O23092004","U24101801","U25022101","O23102703","A24062801","O23112205","U23071901","A23112002","U24111801","O23110101","U24121801","O23120101","O24011202","O23090803","A23060602","A24030402","U25011701","O24011901","U23090801","O24011201","O24011003","O23092202","O23082301","O23091403","O23112901","O23092201","O24013103","O24021403","O24010402","O23092005","O23092203","O24010302","O23090701","O23091501","O23092701","O24022202","O23092802","O23090601","O23100401","O24022102","O23081801","O23092006","O23100503","O23090602","O24013104","O24020803","O24010301","O24010404","O23090802","O23092801","O24022801","O23100706","O23102602","O24021601","O24012401","O24021603","O24022901","O24021501","O23110902","O23102601","O23102704","O23100803","O23102701","O24021401","O24022101","O24030603","O23110901","O23110903","O23110301","O24022301","O23102502","O23111501","O23111602","O24020705","O24021502","O24022201","O23110202","O23090801")
Group2<-c("U23071701","U23052401","U23052201","U24070401","O24011801","O23092003","A24071901","A24072901","O24011102","O23121501","O23092104","O23092001","O23121301","O24020701","O23112201","O23100701","O23100801","O24020903","O24020901","O24020703","O23112204")
Group3<-c("O23100802","O24011205","O23092002","O24011207","O23092103","O23102501","O24011005","A24030401","O24011004","A23051102","U25011702","O24011204","O23121502","O23120702","O24011206","O24021404","O23092101","O24010403","O23112303","O23083001","O23082302","O24010401","O24022302","O24010501","O23112902","O23082303","O23083102","O24013101","O23100402","O24020801","O23120701","O23121304","O24021602","O24011802","O23121306","O23120103","O24020905","O24012403","O24013102","O24021503","O24020904","O23102504","O24013105","O24030601","O23100705","O24030604","O23111601","O24020103","O24030602","O23110302","O23102603","O24031304","O24021402","O24020101","O24012501","O24020804","O23100804","O23102503","O24022902","O24020704","O23110904","O24020102","O24012402","O23102702","O23102604","O23110204","O23110203","O23083101","O23092702")
Group4<-c("O23112304","A23051103","A24071701","A23080101","A24031201","A24080201","O24011105","O23091305","O23121302","O23092803")
Group5<-c("O23091303","O23112301","O24011203","A23112001","O24011001","O24011002","O23091302","O24020902","O23092102","O23091401","O23121503","O23091402","O24020702","O23091301","O23112206","O24011103","O23121305","O24011101","O23121303","O24011104","O23120104","O23100704","O23112302","O23112203","O23100703","O24020805","O24020802","O23112202","O24031302","O23111502","O23100702","O24031301","O24031305","O23082402")
NTC<-c("NTC_2","NTC_3","NTC_16","NTC_4","NTC_6","NTC_7","NTC_8","NTC_9","NTC_11","NTC_12","NTC_13","NTC_14","NTC_15","NTC_10")
PC<-c("PC_1","PC_2","PC_8","PC_4","PC_3","PC_5","PC_6","PC_7","UR009768","UR009909","PC01")
```

## Group 1 vs Group 2

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group1,Group2)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "2")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group1_vs_Group2"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 1 vs Group 3

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group1,Group3)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "3")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group1_vs_Group3"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 1 vs Group 4

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group1,Group4)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "4")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group1_vs_Group4"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 1 vs Group 5

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group1,Group5)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "5")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group1_vs_Group5"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 2 vs Group 3

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group2,Group3)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "3")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group2_vs_Group3"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 2 vs Group 4

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group2,Group4)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "4")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group2_vs_Group4"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 2 vs Group 5

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group2,Group5)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "5")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group2_vs_Group5"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 3 vs Group 4

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group3,Group4)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "4")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group3_vs_Group4"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 3 vs Group 5

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group3,Group5)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "5")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group3_vs_Group5"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```

## Group 4 vs Group 5

```{r, echo=TRUE, warning=FALSE}
ps_deseq_sel <- data.table::copy(ps_deseq)
otu_table(ps_deseq_sel) <- otu_table(ps_deseq)[,c(Group4,Group5)]
diagdds = phyloseq_to_deseq2(ps_deseq_sel, ~Group)
diagdds$Group <- relevel(diagdds$Group, "5")
diagdds <- DESeq(
    diagdds,
    test   = "Wald",
    fitType = "parametric",
    sfType  = "poscounts"  # <- important
)
resultsNames(diagdds)

res = results(diagdds, cooksCutoff = FALSE)
alpha = 0.05
sigtab = res[which(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(phyloseq::tax_table(ps_deseq_sel)[rownames(sigtab), ], "matrix"))
# file base name
fname <- "DEGs_Group4_vs_Group5"
write.xlsx(sigtab, file = paste0(fname, ".xlsx"), rowNames = TRUE)

kable(sigtab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

library("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
        scale_fill_brewer(palette = palname, ...)
}
x = tapply(sigtab$log2FoldChange, sigtab$Order, function(x) max(x))
x = sort(x)
sigtab$Order = factor(as.character(sigtab$Order), levels=names(x))
x = tapply(sigtab$log2FoldChange, sigtab$Family, function(x) max(x))
x = sort(x)
sigtab$Family = factor(as.character(sigtab$Family), levels=names(x))

#ggplot(sigtab, aes(x=log2FoldChange, y=Family, color=Order)) + geom_point(aes(size=padj)) + scale_size_continuous(name="padj",range=c(8,4))+
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust=0.5))
# build the plot
p <- ggplot(sigtab, aes(x = log2FoldChange, y = Family, color = Order)) +
    geom_point(aes(size = padj)) +
    scale_size_continuous(name = "padj", range = c(8, 4)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = -25, hjust = 0, vjust = 0.5))
# SVG (svglite gives crisp text)
if (!requireNamespace("svglite", quietly = TRUE)) install.packages("svglite")
ggplot2::ggsave(paste0(fname, ".svg"), plot = p, device = svglite::svglite,
                                width = 8, height = 6, units = "in", dpi = 300)
# PNG
ggplot2::ggsave(paste0(fname, ".png"), plot = p, device = "png",
                                width = 8, height = 6, units = "in", dpi = 300)
knitr::include_graphics(paste0(fname, ".png"))
```