Master_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f1de1b45",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "\n",
    "# Eunicea spp. Microbiome Analysis\n",
    "\n",
    "This is a synopsis of the bioinformatic processing of microbiome data collected in the doctoral work of Adriana Sarmiento. Microbiome samples of the coral genus *Eunicea* spp. were collected 2019 in the Caribbean Sea and sequenced at Iowa State University.\n",
    "<br>This Jupyter-Notebook is available at\n",
    "https://github.com/julianselke/Eunicea_spp_Microbiome_Analysis .<br><br>\n",
    "\n",
    "Maintainer: Julian Selke<br>\n",
    "Contact: julianselke@web.de<br>\n",
    "Last updated: 28.05.2022<br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddeacfa6",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "# Setup\n",
    "\n",
    "This notebook comprises the entire analysis of the *Eunicea* spp. microbiome dataset. It somewhat still presents a medley of various scripts produced as analysis progressed. One inconsistency remains: Some functions in ```R``` require a transposed _feature table_. If any error occurs, transposing the respective _feature table_ might present a quick fix. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f74706a4",
   "metadata": {},
   "source": [
    "In order to run ```R``` commands in this notebook, either the ```R```-kernel has to be running or ```rpy2``` has to be loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6b3dc06",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47c51bf8",
   "metadata": {},
   "source": [
    "Now, ```%%R``` in a cell enables *R* code to be translated to ```python``` and to be run, and ```%%bash``` enables *bash* code to be run. <br>\n",
    "However, some functions cannot be translated and therefore require a running ```R```-kernel. The kernel has to be changed manually and you will be advised when to do so."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1ea9d51",
   "metadata": {},
   "source": [
    "Check current workind directory. Please note that ```!``` at the beginning of a line initiates the ```sh``` shell, not ```bash```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f7beb95",
   "metadata": {},
   "outputs": [],
   "source": [
    "! realpath ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8399b08",
   "metadata": {},
   "source": [
    "Check if conda environment is working correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37c1a937",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "qiime --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b030bc3",
   "metadata": {},
   "source": [
    "Line numbers can be displayed by pressing ```SHIFT + L``` when no cell is activated."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6c5fa91",
   "metadata": {},
   "source": [
    "To adjust plotting options of ```R``` cells add and run the follwing cell: \n",
    "\n",
    "```options(repr.plot.width = X, repr.plot.height = Y, repr.plot.res = Z)```\n",
    "\n",
    "Replace ```X```, ```Y```, and ```Z``` with desired values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1157dbf6",
   "metadata": {},
   "source": [
    "Press ```Shift```+```O``` to see all output of the selected cell. Useful, e.g., when tables are cut off. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1c009f1",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "# Content\n",
    "\n",
    "[fastp](#fastp) <br>\n",
    "[Qiime2](#Qiime2) <br>\n",
    "[Preparation](#Preparation) <br>\n",
    "[Exploration](#Exploration) <br>\n",
    "> [Euclidian Distance](#EuclidianDistance) <br>\n",
    "[Bray-Curtis Dissimilarity](#Bray-CurtisDissimilarity) <br>\n",
    "[Jaccard Distance](#JaccardDistance) <br>\n",
    "[Unweighted UniFrac Distance](#UnweightedUniFracDistance) <br>\n",
    "[Weighted UniFrac Distance](#WeightedUniFracDistance) <br>\n",
    "[Aitchison Distance](#AitchisonDistance) <br>\n",
    "[PhILR Distance](#PhILRDistance) <br>\n",
    "\n",
    "[Statistics](#Statistics) <br>\n",
    "[t-SNE](#t-SNE) <br>\n",
    "[Composition](#Composition) <br>\n",
    "[iNEXT](#iNEXT) <br>\n",
    "[Core Microbiome](#CoreMicrobiome) <br>\n",
    "[Host Phylogeny](#HostPhylogeny) <br>\n",
    "[ALDEx2](#ALDEx2) <br>\n",
    "\n",
    "<br><hr><br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "868f689a",
   "metadata": {},
   "source": [
    "# fastp\n",
    "\n",
    "Preprocessing with fastp enables base correction in overlapping regions based on PHRED scores. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fd6e463",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "echo '\n",
    "\n",
    "#!/usr/bin/env bash\n",
    "\n",
    "# Default values of arguments\n",
    "OUTPUT_DIRECTORY=./fastp_output\n",
    "\n",
    "# Loop through arguments and process them\n",
    "for arg in \"$@\"\n",
    "do\n",
    "    case $arg in\n",
    "        -i|--input_directory)\n",
    "\t# get list of all files in directory\n",
    "        INPUT_DIRECTORY=$2\n",
    "        shift # Remove --input_directory from processing\n",
    "        ;;\n",
    "        -o|--output_directory)\n",
    "\tOUTPUT_DIRECTORY=$3\n",
    "        shift # Remove --output_directory from processing\n",
    "\t;;\n",
    "        -h|--help)\n",
    "\techo\n",
    "\tprintf \"\\033[1;32m  Process files of a directory with fastp. \\033[0m\\n\"\n",
    "        echo\n",
    "        echo \"  Syntax: bash $0 [-i -o|-h]\"\n",
    "        echo\n",
    "        echo \"  options:\"\n",
    "        echo\n",
    "        echo \"   \\033[1;30m-i     The input directory containing fastq files.\"\n",
    "        echo \"          \\033[0mdefault: \\e[47m.\\e[49m\"\n",
    "        echo\n",
    "        echo \"   \\033[1;30m-o     The output directory to save processed files to.\"\n",
    "        echo \"          \\033[0mdefault: \\e[47m./fastp_output\\e[49m\" \n",
    "        echo\n",
    "        echo \"   \\033[1;30m-h     Print this Help.\"\n",
    "        echo\n",
    "\texit 0\n",
    "    esac\n",
    "done\n",
    "\n",
    "INPUT_FILES_PATH=$(find $INPUT_DIRECTORY -type f | sed \"s/R._001\\.fastq\\.gz//g\" | uniq)\n",
    "\n",
    "mkdir $OUTPUT_DIRECTORY\n",
    "\n",
    "for file in $INPUT_FILES_PATH\n",
    "do\n",
    "  echo \"Processing $file\"\n",
    "  fastp -i \"$file\"R1_001.fastq.gz \\\n",
    "        -I \"$file\"R2_001.fastq.gz \\\n",
    "        -o \"$OUTPUT_DIRECTORY\"/$(basename $file)R1_001.fastq.gz \\\n",
    "        -O \"$OUTPUT_DIRECTORY\"/$(basename $file)R2_001.fastq.gz -c  \n",
    "done\n",
    "\n",
    "' > preprocess_fastp.sh "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "934187c7",
   "metadata": {},
   "source": [
    "Give permission to execute script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "622ad777",
   "metadata": {},
   "outputs": [],
   "source": [
    "! chmod +x preprocess_fastp.sh "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d153b29",
   "metadata": {},
   "source": [
    "See help options of script. See also: https://github.com/OpenGene/fastp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8d4832e",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ./preprocess_fastp.sh -h"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469cf01f",
   "metadata": {},
   "source": [
    "JSON and HTML reports are overwritten in each run. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d00c6ba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ./preprocess_fastp.sh -i ./raw_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1e992db",
   "metadata": {},
   "source": [
    "<br><hr><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bed595f",
   "metadata": {},
   "source": [
    "# Qiime2\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "485dc7e0",
   "metadata": {},
   "source": [
    "### manifest\n",
    "\n",
    "Write the manifest file to import forward and reverse sequences into Qiime2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6f8328b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "ls -p /root/biommar/fastp_output/ | grep -v / | \\\n",
    "awk -F'_' 'BEGIN {OFS=\"\"; print \"sample-id,absolute-filepath,direction\"} \\\n",
    "{ if ($0 ~ /.+_R1_.+/) {print $1,\",\",\"/root/biommar/fastp_output/\",$0,\",forward\" \\\n",
    "} else {print $1,\",\",\"/root/biommar/fastp_output/\",$0,\",reverse\"}}' \\\n",
    "> /root/biommar/data/eunicea_PE_manifest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba08d111",
   "metadata": {},
   "source": [
    "### Import raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc8b390c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime tools import \\\n",
    "  --type 'SampleData[PairedEndSequencesWithQuality]' \\\n",
    "  --input-path /root/biommar/data/eunicea_PE_manifest \\\n",
    "  --output-path /root/biommar/data/eunicea_paired_end_demux.qza \\\n",
    "  --input-format PairedEndFastqManifestPhred33"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cfcd351",
   "metadata": {},
   "source": [
    "### Trimming\n",
    "\n",
    "Doc.: https://docs.qiime2.org/2022.2/plugins/available/dada2/denoise-paired/\n",
    "\n",
    "Denoise and dereplicate the sequences, filters chimeras and merge the paired-end sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "849e2d43",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime dada2 denoise-paired \\\n",
    "  --i-demultiplexed-seqs /root/biommar/data/eunicea_paired_end_demux.qza \\\n",
    "  --p-trunc-len-f 230 \\\n",
    "  --p-trunc-len-r 165 \\\n",
    "  --p-n-threads 4 \\\n",
    "  --o-table /root/biommar/data/dada_table.qza \\\n",
    "  --o-representative-sequences /root/biommar/data/dada_rep_seqs.qza \\\n",
    "  --o-denoising-stats /root/biommar/data/dada_stats.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7308d6e6",
   "metadata": {},
   "source": [
    "### Filtering\n",
    "\n",
    "Filter features to at leat 10 occurrences across all samples and sequnces based on featurs subsequently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac8316fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime feature-table filter-features \\\n",
    "  --i-table /root/biommar/data/dada_table.qza \\\n",
    "  --p-min-frequency 10 \\\n",
    "  --o-filtered-table /root/biommar/data/dada_filtered_table.qza\n",
    "\n",
    "qiime feature-table filter-seqs \\\n",
    "  --i-data /root/biommar/data/dada_rep_seqs.qza\\\n",
    "  --i-table /root/biommar/data/dada_filtered_table.qza \\\n",
    "  --o-filtered-data /root/biommar/data/dada_filtered_rep_seqs.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfa29a40",
   "metadata": {},
   "source": [
    "### Phylogeny\n",
    "\n",
    "Fragment Insertion\n",
    "\n",
    "Doc.: https://library.qiime2.org/plugins/q2-fragment-insertion/16/\n",
    "\n",
    "Paper: SEPP: SATé-Enabled Phylogenetic Placement\n",
    "\n",
    "Phylogenetic placement of retained sequences.\n",
    "In order to explore the phylogenetic composition of our samples, we place the probes into a prebuild phylogenetic tree. The utilized tool, SEPP, has proven to be advantagious in working with sOTUs, since it is using a controlled reference tree, but also allows non-perfect hits to its tips.\n",
    "The phylogeny of the sample may offer biologically relevant information and allows are more sophisticated comparision of microbial communities.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ac93779",
   "metadata": {},
   "source": [
    "<hr>\n",
    "Download the reference database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b52918a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd /root/biommar/data/\n",
    "! wget https://data.qiime2.org/2022.2/common/sepp-refs-silva-128.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a3027d3",
   "metadata": {},
   "source": [
    "<hr> \n",
    "\n",
    "Check if any errors occurred while downloading. The warning ``` no properly formatted MD5 checksum lines found ``` can be ignored. You should see an **OK** in the second line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83e06dbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash \n",
    "\n",
    "md5sum -c <<<\"7879792a6f42c5325531de9866f5c4de  /root/biommar/data/sepp-refs-silva-128.qza\" \\\n",
    "- /root/biommar/data/sepp-refs-silva-128.qza | tail -n -1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77210660",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "Run the fragment insertion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "817778e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime fragment-insertion sepp \\\n",
    "  --i-representative-sequences /root/biommar/data/dada_filtered_rep_seqs.qza  \\\n",
    "  --i-reference-database /root/biommar/data/sepp-refs-silva-128.qza \\\n",
    "  --p-threads 4 \\\n",
    "  --o-tree /root/biommar/data/sepp_phylogeny.qza \\\n",
    "  --o-placements /root/biommar/data/sepp_placements.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e182986",
   "metadata": {},
   "source": [
    "### --> not enough RAM (16GB insufficient)\n",
    "--> run on cluster instead"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8efdf9c3",
   "metadata": {},
   "source": [
    "### Read Classification\n",
    "\n",
    "Doc.: https://docs.qiime2.org/2022.2/plugins/available/feature-classifier/classify-sklearn/\n",
    "\n",
    "Classification of reads by taxon using a pre-fitted classifier.\n",
    "Without the knowledge of the species, genera, etc. present in the samples, the biological interpretation of the communities is very limited. As an example, no mutualist/pathogen relations can be estimated, probale metabolic processes and therefore biological function of the species remain veiled. Furthermore, this can inform us about contamination still present in the samples, e.g. species present in the human microbiome but not in the soil may indicate contamination of respective samples while collecting or processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ea3ef5",
   "metadata": {},
   "source": [
    "<hr>\n",
    "Download the prefitted classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c4497f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd /root/biommar/data/\n",
    "! wget https://data.qiime2.org/2022.2/common/silva-138-99-515-806-nb-classifier.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5631618a",
   "metadata": {},
   "source": [
    "<hr> \n",
    "\n",
    "Check if any errors occurred while downloading. The warning ``` no properly formatted MD5 checksum lines found ``` can be ignored. You should see an **OK** in the second line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8923c4c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash \n",
    "\n",
    "md5sum -c <<<\"e05afad0fe87542704be96ff483824d4  /root/biommar/data/silva-138-99-515-806-nb-classifier.qza\" \\\n",
    "- /root/biommar/data/silva-138-99-515-806-nb-classifier.qza | tail -n -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37779ec1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime feature-classifier classify-sklearn \\\n",
    "  --i-classifier /root/biommar/data/silva-138-99-515-806-nb-classifier.qza \\\n",
    "  --i-reads /root/biommar/data/dada_filtered_rep_seqs.qza \\\n",
    "  --o-classification /root/biommar/data/filtered_taxonomy_silva.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17506ea1",
   "metadata": {},
   "source": [
    "### --> not enough RAM (16GB insufficient)\n",
    "--> run on cluster instead"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05368dbe",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "Filter features to match metadata.\n",
    "\n",
    "Doc.: https://docs.qiime2.org/2021.11/plugins/available/feature-table/filter-samples/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28d0f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime feature-table filter-samples \\\n",
    "  --i-table /root/biommar/data/dada_filtered_table.qza \\\n",
    "  --m-metadata-file /root/biommar/data/Metadata_Eunicea_Bacteria.tsv \\\n",
    "  --o-filtered-table /root/biommar/data/dada_filtered_table_metadata_filtered.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4820037f",
   "metadata": {},
   "source": [
    "### Filter out Chloroplast and Mitochindria\n",
    "\n",
    "Doc.: https://docs.qiime2.org/2022.2/plugins/available/feature-table/filter-features/\n",
    "\n",
    "Contamination belonging to these two taxa/organelles is present in our data and blurring the distance metrics. A SQLite command is used to specify the filter applied on taxa.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8820ab9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime feature-table filter-features \\\n",
    "  --i-table /root/biommar/data/dada_filtered_table_metadata_filtered.qza  \\\n",
    "  --m-metadata-file /root/biommar/data/filtered_taxonomy_silva.qza \\\n",
    "  --p-where  \"(Taxon NOT LIKE '%c__Chloroplast%') AND (Taxon NOT LIKE '%f__mitochondria%')\" \\\n",
    "  --o-filtered-table /root/biommar/data/dada_table_no_mitochloro_silva.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89f45535",
   "metadata": {},
   "source": [
    "### Extract Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e51a037",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "cd  /root/biommar/data\n",
    "unzip -qq dada_table_no_mitochloro_silva.qza\n",
    "cp $(echo $(ls -td -- */ | head -n 1) $(echo \"data/feature-table.biom\") | sed 's/\\ //g') .\n",
    "biom convert -i feature-table.biom -o feature-table.tsv --to-tsv\n",
    "rm -r $(ls -td -- */ | head -n 1)\n",
    "unzip -qq filtered_taxonomy_silva.qza\n",
    "cp $(echo $(ls -td -- */ | head -n 1) $(echo \"data/taxonomy.tsv\") | sed 's/\\ //g') . \n",
    "rm -r $(ls -td -- */ | head -n 1)\n",
    "unzip -qq sepp_phylogeny.qza \n",
    "cp $(echo $(ls -td -- */ | head -n 1) $(echo \"data/tree.nwk\") | sed 's/\\ //g') .\n",
    "rm -r $(ls -td -- */ | head -n 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efdd1f6a",
   "metadata": {},
   "source": [
    "<br><hr><br><br><hr><br>\n",
    "\n",
    "# Analysis in R\n",
    "\n",
    "<span style=\"color:red\">**PLEASE CHANGE KERNEL ACCORDINGLY!**</span>\n",
    "\n",
    "Select ```Kernel``` in the menu bar at the top, click on ```Change kernel```, and select ```R```.\n",
    "\n",
    "<br><hr><br><br><hr><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeffbabc",
   "metadata": {},
   "outputs": [],
   "source": [
    "library( ALDEx2,        quietly = TRUE)\n",
    "library( ape,           quietly = TRUE)\n",
    "library( aplot,         quietly = TRUE)\n",
    "library( cluster,       quietly = TRUE)\n",
    "library( ggdendro,      quietly = TRUE)\n",
    "library( grid,          quietly = TRUE)\n",
    "library( glmnet,        quietly = TRUE)\n",
    "library( ggplot2,       quietly = TRUE)\n",
    "library( ggfortify,     quietly = TRUE)\n",
    "library( ggthemes,      quietly = TRUE)\n",
    "library( gridExtra,     quietly = TRUE)\n",
    "library( ggpmisc,       quietly = TRUE)\n",
    "library( ggsci,         quietly = TRUE)\n",
    "library( ggtree,        quietly = TRUE)\n",
    "library( tidyr,         quietly = TRUE)\n",
    "library( metagenomeSeq, quietly = TRUE)\n",
    "library( microbiome,    quietly = TRUE)\n",
    "library( phyloseq,      quietly = TRUE)\n",
    "library( philr,         quietly = TRUE)\n",
    "library( mia,           quietly = TRUE)\n",
    "library( vegan,         quietly = TRUE)\n",
    "library( multcompView,  quietly = TRUE)\n",
    "library( plyr,          quietly = TRUE)\n",
    "library( dplyr,         quietly = TRUE)\n",
    "library( Rtsne,         quietly = TRUE)\n",
    "library( plotly,        quietly = TRUE)\n",
    "library( zgtools,       quietly = TRUE)\n",
    "library( reshape,       quietly = TRUE)\n",
    "library( onewaytests,   quietly = TRUE)\n",
    "library( iNEXT,         quietly = TRUE)\n",
    "library( phytools,      quietly = TRUE)\n",
    "library( reshape2,      quietly = TRUE)\n",
    "library( repr,          quietly = TRUE)\n",
    "library( kableExtra,    quietly = TRUE)\n",
    "library( IRdisplay,     quietly = TRUE)\n",
    "library( jsonlite,      quietly = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa06e369",
   "metadata": {},
   "source": [
    "<br><hr><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2621985e",
   "metadata": {},
   "source": [
    "# Preparation\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addc96cd",
   "metadata": {},
   "source": [
    "### Import Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5958e741",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_feature_table <- data.frame(t(read.delim(\"/root/biommar/data/feature-table.tsv\",\n",
    "                                             skip        = 1, \n",
    "                                             row.names   = 1, \n",
    "                                             check.names = FALSE)),\n",
    "                                check.names = FALSE)\n",
    "\n",
    "bac_metadata <- read.table(\"/root/biommar/data/Metadata_Eunicea_Bacteria.tsv\",\n",
    "                           header      = T, \n",
    "                           sep         = \"\\t\", \n",
    "                           check.names = FALSE)[, 1:12]\n",
    "\n",
    "bac_taxonomy <- read.table(\"/root/biommar/data/taxonomy.tsv\",\n",
    "                           header    = T,\n",
    "                           row.names = 1,\n",
    "                           sep       = \"\\t\")\n",
    "\n",
    "# filter out outgroups\n",
    "bac_metadata <- bac_metadata[!bac_metadata$genus %in% c(\"Water sample\", \"Muricea\"),]\n",
    "\n",
    "# assign row names as sample ID to facilitate merging on row names\n",
    "rownames(bac_metadata) <- bac_metadata$sampleid\n",
    "\n",
    "# make sure to exclude features from other objects as well\n",
    "bac_feature_table <- bac_feature_table[rownames(bac_feature_table) %in% rownames(bac_metadata),]\n",
    "\n",
    "bac_taxonomy <- bac_taxonomy[rownames(bac_taxonomy) %in% colnames(bac_feature_table),]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba55d989",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Prepare Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "872ff4fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter taxonomy to remove discarded ASVs (Chloroplasts and Mitochondria)\n",
    "# !!! asymmetric operation !!! setdiff(b, a) != setdiff(a, b)\n",
    "out <- setdiff(rownames(bac_taxonomy), colnames(bac_feature_table))\n",
    "# remove missing samples\n",
    "# bac_taxonomy <- bac_taxonomy[-c(which(rownames(bac_taxonomy) %in% out)), ]\n",
    "\n",
    "# same for metadata\n",
    "out <- setdiff(rownames(bac_metadata), rownames(bac_feature_table)) \n",
    "# remove missing samples\n",
    "bac_metadata <- bac_metadata[-c(which(rownames(bac_metadata) %in% out)), ]\n",
    "\n",
    "# make all entries numeric\n",
    "bac_feature_table <- mutate_all(bac_feature_table, function(x) as.numeric(as.character(x)))\n",
    "\n",
    "# split taxonomy\n",
    "taxa <- separate(\n",
    "  data = bac_taxonomy,\n",
    "  col  = Taxon,\n",
    "  into = c(\n",
    "    \"kingdom\",\n",
    "    \"phylum\",\n",
    "    \"class\",\n",
    "    \"order\",\n",
    "    \"family\",\n",
    "    \"genus\",\n",
    "    \"species\"\n",
    "  ),\n",
    "  sep = \";\"\n",
    ")\n",
    "\n",
    "# store row names\n",
    "taxa_names <- rownames(taxa)\n",
    "# denote all cells without information as \"unidentified\" \n",
    "taxa[is.na(taxa)] <- \"unidentified\"\n",
    "taxa <- as.data.frame(lapply(taxa, function(x) gsub(\"[a-z]__$\", \"unidentified\", x)))\n",
    "# crop leading pattern\n",
    "taxa <- as.data.frame(lapply(taxa, function(x) gsub(\"[a-z]__\", \"\", x)))\n",
    "# remove leading whitespaces\n",
    "taxa <- as.data.frame(sapply(taxa, function(x) {gsub(\"^ \", \"\", x)}))\n",
    "# re-assign row names\n",
    "rownames(taxa) <- taxa_names\n",
    "# cleanup\n",
    "rm(taxa_names)\n",
    "\n",
    "write.table(taxa[,1:7],\n",
    "            file  = \"/root/biommar/data/taxa.txt\",\n",
    "            sep   = \"\\t\",\n",
    "            quote = F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ea5eb00",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "407e3454",
   "metadata": {},
   "outputs": [],
   "source": [
    "# redundant, since this filtering has been done in Qiime2 before\n",
    "# subset the feature table to contain only features that have at least ten reads\n",
    "bac_feature_table <- bac_feature_table[, (colSums(bac_feature_table) > 9)]\n",
    "\n",
    "bac_taxonomy <- bac_taxonomy[rownames(bac_taxonomy) %in% colnames(bac_feature_table), ]\n",
    "\n",
    "# subset taxa as well\n",
    "taxa <- taxa[rownames(taxa) %in% colnames(bac_feature_table), ]\n",
    "\n",
    "# filter data: remove samples with <1000 reads\n",
    "NumReads <- rowSums(bac_feature_table)\n",
    "# remove respective samples\n",
    "bac_feature_table <- bac_feature_table[NumReads > 1000, ]\n",
    "# cleanup\n",
    "rm(NumReads)\n",
    "\n",
    "# make sure filtered features are excluded\n",
    "bac_metadata <- bac_metadata[rownames(bac_metadata) %in% rownames(bac_feature_table),]\n",
    "\n",
    "bac_taxonomy <- bac_taxonomy[rownames(bac_taxonomy) %in% colnames(bac_feature_table),]\n",
    "\n",
    "taxa <- \n",
    "  taxa[rownames(taxa) %in% colnames(bac_feature_table),]\n",
    "\n",
    "# calculate sparcity, i.e. percentage of cells with zero count\n",
    "sparcity <- length(bac_feature_table[bac_feature_table == 0]) / \n",
    "  (nrow(bac_feature_table) * ncol(bac_feature_table)) * 100\n",
    "\n",
    "print(\"Sparcity, i.e. percentage of cells in feature table with zero-count is:\")\n",
    "print(sparcity)\n",
    "\n",
    "# filter for taxa without occurrences in remaining samples \n",
    "taxa_bool <- ifelse(bac_feature_table > 0, 1, 0)\n",
    "# aggregate by feature\n",
    "taxa_occ <- as.vector(colSums(taxa_bool))\n",
    "taxa_occ_abs <- as.vector(colSums(bac_feature_table))\n",
    "\n",
    "# remove respective taxa\n",
    "bac_feature_table <- bac_feature_table[, taxa_occ > 0]\n",
    "\n",
    "# aggregate by sample\n",
    "sample_occ <- as.vector(rowSums(taxa_bool))\n",
    "sample_occ_abs <- as.vector(rowSums(bac_feature_table))\n",
    "\n",
    "# keep only samples with more than one feature\n",
    "bac_feature_table <- bac_feature_table[sample_occ > 1, ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b90d993",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d0d76c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set option number display\n",
    "options(scipen=999)\n",
    "options(repr.plot.width = 4, repr.plot.height = 4, repr.plot.res = 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e940a3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_taxa_occ <- ggplot(data.frame(taxa_occ), aes(x = taxa_occ)) + \n",
    "  geom_histogram(color = \"gray10\", fill = \"#00468BFF\", bins = 30) + \n",
    "  scale_y_log10() +\n",
    "  theme_bw() +\n",
    "  ylab(\"Number of Features\") +\n",
    "  xlab(\"Occurrences [Number of Samples]\") +\n",
    "  ggtitle(\"Feature Occurrences\") +\n",
    "  theme(plot.title = element_text(hjust=0.5))\n",
    "\n",
    "plot_taxa_occ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682f6e32",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_taxa_occ_abs <- ggplot(data.frame(taxa_occ_abs), aes(x = taxa_occ_abs)) + \n",
    "  geom_histogram(color = \"gray10\", fill = \"#00468BFF\", bins = 30) + \n",
    "  scale_x_log10(breaks = c(1,10,100,1000,10000,100000,1000000)) +\n",
    "  theme_bw() +\n",
    "  ylab(\"Number of Features\") +\n",
    "  xlab(\"Count across all Samples\") +\n",
    "  ggtitle(\"Feature Counts\") +\n",
    "  theme(plot.title = element_text(hjust=0.5))\n",
    "\n",
    "plot_taxa_occ_abs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79c268b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_taxa_occ_abs <- ggplot(data.frame(taxa_occ_abs), aes(x = taxa_occ_abs)) + \n",
    "  geom_histogram(color = \"gray10\", fill = \"#00468BFF\", bins = 30) + \n",
    "  scale_x_log10(breaks = c(1,10,100,1000,10000,100000,1000000)) +\n",
    "  theme_bw() +\n",
    "  ylab(\"Number of Features\") +\n",
    "  xlab(\"Count across all Samples\") +\n",
    "  ggtitle(\"Feature Counts\") +\n",
    "  theme(plot.title = element_text(hjust=0.5))\n",
    "\n",
    "plot_taxa_occ_abs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e822579f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sample_occ_abs <- ggplot(data.frame(sample_occ_abs), aes(x = sample_occ_abs)) + \n",
    "  geom_histogram(color = \"gray10\", fill = \"#00468BFF\", bins = 30) + \n",
    "  scale_x_log10() +\n",
    "  theme_bw() +\n",
    "  xlab(\"Count of Features\") +\n",
    "  ylab(\"Number of Samples\") +\n",
    "  ggtitle(\"Feature Counts per Sample\") +\n",
    "  theme(plot.title = element_text(hjust=0.5))\n",
    "\n",
    "\n",
    "plot_sample_occ_abs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b19754",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Write Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "645ede47",
   "metadata": {},
   "outputs": [],
   "source": [
    "write.table(bac_feature_table,\n",
    "            file      = \"/root/biommar/data/B_asv.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)\n",
    "\n",
    "write.table(bac_metadata,\n",
    "            file      = \"/root/biommar/data/B_asv_metadata.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)\n",
    "\n",
    "write.table(bac_taxonomy,\n",
    "            file      = \"/root/biommar/data/B_asv_taxonomy.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99eb8bc3",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling (CSS) normalization to remove biases in count data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c58e061",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a MRexperiment object from count data.\n",
    "table_metagenomeseq = newMRexperiment(t(bac_feature_table))\n",
    "# calculate the percentile for which to sum counts up to and scale by\n",
    "pthq = cumNormStatFast(table_metagenomeseq)\n",
    "# calculate each column's quantile and calculate the sum up to and \n",
    "# including that quantile\n",
    "table_metagenomeseq_norm = cumNorm(table_metagenomeseq, p = pthq)\n",
    "\n",
    "bac_feature_table_css <- data.frame(MRcounts(table_metagenomeseq_norm, norm=TRUE))\n",
    "bac_feature_table_css <- as.data.frame(t(bac_feature_table_css))\n",
    "\n",
    "# cleanup\n",
    "rm(pthq, table_metagenomeseq_norm, table_metagenomeseq)\n",
    "\n",
    "# make sure to exclude features from other objects as well\n",
    "bac_taxonomy_css <- bac_taxonomy[rownames(bac_taxonomy) %in% rownames(bac_feature_table_css),]\n",
    "\n",
    "bac_metadata_css <- bac_metadata[rownames(bac_metadata) %in% colnames(bac_feature_table_css),]\n",
    "\n",
    "bac_taxa_css <- taxa[rownames(taxa) %in% rownames(bac_feature_table_css),]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b49e6a",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Write CSS Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3241ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "write.table(bac_feature_table_css,\n",
    "            file      = \"/root/biommar/data/CSS_B_asv.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)\n",
    "\n",
    "write.table(bac_metadata_css,\n",
    "            file      = \"/root/biommar/data/CSS_B_asv_metadata.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)\n",
    "\n",
    "write.table(bac_taxonomy_css,\n",
    "            file      = \"/root/biommar/data/CSS_B_asv_taxonomy.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a7e5e12",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### PhILR "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e321b058",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.file  <- file.path(\"/root/biommar/data/B_asv.txt\")\n",
    "tax.file  <- file.path(\"/root/biommar/data/taxa.txt\")\n",
    "meta.file <- file.path(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "bac_ps <- read_phyloseq(\n",
    "  otu.file      = otu.file,\n",
    "  taxonomy.file = tax.file,\n",
    "  metadata.file = meta.file,\n",
    "  type          = \"simple\",\n",
    "  sep           = \"\\t\"\n",
    ")\n",
    "\n",
    "meta <- read.table(meta.file, row.names = 1, header = T, sep=\"\\t\")\n",
    "\n",
    "# import tree\n",
    "tree <- read_tree(\"/root/biommar/data/tree.nwk\")\n",
    "\n",
    "bac_ps <- merge_phyloseq(bac_ps, tree)\n",
    "\n",
    "tse <- mia::makeTreeSummarizedExperimentFromPhyloseq(bac_ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a483f189",
   "metadata": {},
   "outputs": [],
   "source": [
    "tse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2941aec0",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Filter Extremely Low-Abundance OTUs\n",
    "\n",
    "Prevalence is a measurement, which describes in how many samples certain\n",
    "microbes were detected.\n",
    "\n",
    "A 0.01% relative abundance threshold is set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60319c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "tse <-  tse %>% subsetByPrevalentTaxa(\n",
    "  detection   = 0.01/100,\n",
    "  prevalence  = 0,\n",
    "  as_relative = TRUE)\n",
    "\n",
    "# Collapse the tree\n",
    "# Otherwise the original tree with all nodes is kept\n",
    "# (including those that were filtered out from rowData)\n",
    "tree <- ape::keep.tip(phy = rowTree(tse), \n",
    "                      tip = rowLinks(tse)$nodeNum)\n",
    "rowTree(tse) <- tree\n",
    "\n",
    "## Add a new assay with a pseudocount \n",
    "assays(tse)$counts.shifted <- assays(tse)$counts + 1e-12\n",
    "\n",
    "# enumerate nodes\n",
    "#tree <- makeNodeLabel(tree, method=\"number\", prefix='n')\n",
    "\n",
    "# check if tree is rooted and if tree is binary\n",
    "if(!is.rooted(tree)){ \n",
    "  system(\"echo '\\033[1;31mError: Phylogenetic tree is not rooted!\\033[0m'\")}\n",
    "if(!is.binary(tree)){ \n",
    "  system(\"echo '\\033[1;31mError: Phylogenetic tree is not binary!\\033[0m'\")}\n",
    "# Add the modified tree back to the (`TreeSE`) data object \n",
    "rowTree(tse) <- tree\n",
    "\n",
    "bac_feature_table_philr <- as.data.frame(philr(tse, \n",
    "                                               part.weights = 'uniform', \n",
    "                                               ilr.weights  = 'uniform', \n",
    "                                               return.all   = FALSE, \n",
    "                                               abund_values = \"counts.shifted\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd5ca574",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Write PhILR Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a349439e",
   "metadata": {},
   "outputs": [],
   "source": [
    "write.table(bac_feature_table_philr,\n",
    "            file      = \"/root/biommar/data/PhILR_B_asv.txt\",\n",
    "            sep       = \"\\t\",\n",
    "            row.names = T,\n",
    "            quote     = F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f09481d",
   "metadata": {},
   "source": [
    "<br><hr><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f629470",
   "metadata": {},
   "source": [
    "# Exploration\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc93d1c3",
   "metadata": {},
   "source": [
    "\n",
    "## Remark:\n",
    "<br>\n",
    "The following distance and dissimilarity metrics are applied to count data of samples.\n",
    "\n",
    "Two filters were applied a priori: ASVs with less than 10 counts across all samples and samples with less than 1000 reads where excluded.\n",
    "\n",
    "Where indicated, counts where transformed by **C**umulative **S**um **S**caling.<br>\n",
    "\n",
    "https://doi.org/10.1038/nmeth.2658\n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c612deb",
   "metadata": {},
   "source": [
    "<a id='EuclidianDistance'></a>\n",
    "\n",
    "<br><hr><br>\n",
    "# Euclidian Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "\n",
    "The euclidian distance between two points is defined as\n",
    "<br><br>\n",
    "$$d_e(p,q) = \\sqrt{\\sum_{i=1}^{n}(p_i-q_i)^2}$$\n",
    "<br>\n",
    "where *n* is the number of dimensions.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eff9eb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_feature_table <- read.delim(\"/root/biommar/data/B_asv.txt\", check.names = F)\n",
    "bac_metadata <- read.delim(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "# transpose and set class of data for computation\n",
    "bac_feature_table <- as.data.frame(bac_feature_table, check.names = F)\n",
    "bac_feature_table <- mutate_all(bac_feature_table, \n",
    "                                    function(x) as.numeric(as.character(x)))\n",
    "\n",
    "bac_feature_table_css <- read.delim(\"/root/biommar/data/CSS_B_asv.txt\", check.names = F)\n",
    "bac_metadata_css <- read.delim(\"/root/biommar/data/CSS_B_asv_metadata.txt\")\n",
    "\n",
    "# transpose and set class of data for computation\n",
    "bac_feature_table_css <- as.data.frame(t(bac_feature_table_css), check.names = F)\n",
    "bac_feature_table_css <- mutate_all(bac_feature_table_css, \n",
    "                                    function(x) as.numeric(as.character(x)))\n",
    "\n",
    "# calculate distance matrix\n",
    "bac_euclidian_dist <- vegan::vegdist(bac_feature_table, method = \"euclidian\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_euclidian_pcoa <- ape::pcoa(bac_euclidian_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_euclidian_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_euclidian_pcoa$vectors)),\n",
    "    X = bac_euclidian_pcoa$vectors[, 1],\n",
    "    Y = bac_euclidian_pcoa$vectors[, 2],\n",
    "    Y = bac_euclidian_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_euclidian_pcoa_data$Sample <- rownames(bac_euclidian_pcoa_data)\n",
    "\n",
    "colnames(bac_euclidian_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_euclidian_pcoa_data <-\n",
    "  left_join(bac_euclidian_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "                                    \n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[1]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[2]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[3]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdf68733",
   "metadata": {},
   "source": [
    "##### Remark: \n",
    "If necessary, please adjust ```width``` and ```height``` based on your screens dimensions.<br>\n",
    "Uncomment ```%>% zgtools::plotly_elegant()``` to view in dark theme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5aa112",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_euclidian_pcoa_data,\n",
    "    x      = ~X,\n",
    "    y      = ~Y,\n",
    "    z      = ~Z,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Euclidian Distances\") #%>% zgtools::plotly_elegant()\n",
    "                                    \n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d80d41",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cefbf0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_euclidian_nmds_2D <- metaMDS(bac_feature_table,k=2,trymax=20,distance=\"euclidian\")\n",
    "bac_euclidian_nmds_3D <- metaMDS(bac_feature_table,k=3,trymax=20,distance=\"euclidian\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f48ed059",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 7, repr.plot.height = 7, repr.plot.res = 200)\n",
    "\n",
    "cat(\"\\n\", \"Stress 3D: \", bac_euclidian_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_euclidian_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_euclidian_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a3b6429",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_euclidian_nmds_3D_plot_data <- as.data.frame(bac_euclidian_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_euclidian_nmds_3D_plot_data$sampleid <- rownames(bac_euclidian_nmds_3D_plot_data)\n",
    "\n",
    "bac_euclidian_nmds_3D_plot_data <-\n",
    "  left_join(bac_euclidian_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_euclidian_nmds_3D_plot_data$values$Relative_eig[[1]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_euclidian_nmds_3D_plot_data$values$Relative_eig[[2]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_euclidian_nmds_3D_plot_data$values$Relative_eig[[3]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d13c178e",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_euclidian_nmds_3D_plot_data,\n",
    "    x      = ~MDS1,\n",
    "    y      = ~MDS2,\n",
    "    z      = ~MDS3,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nNMDS on Euclidian Distances\") #%>% zgtools::plotly_elegant()\n",
    "                                    \n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6be1df5",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19b41ef0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate distance matrix\n",
    "bac_euclidian_dist <- vegan::vegdist(bac_feature_table_css, method = \"euclidian\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_euclidian_pcoa <- ape::pcoa(bac_euclidian_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_euclidian_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_euclidian_pcoa$vectors)),\n",
    "    X = bac_euclidian_pcoa$vectors[, 1],\n",
    "    Y = bac_euclidian_pcoa$vectors[, 2],\n",
    "    Z = bac_euclidian_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_euclidian_pcoa_data$Sample <- rownames(bac_euclidian_pcoa_data)\n",
    "\n",
    "colnames(bac_euclidian_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_euclidian_pcoa_data <-\n",
    "  left_join(bac_euclidian_pcoa_data, bac_metadata_css, by = \"sampleid\")\n",
    "  \n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[1]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[2]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_euclidian_pcoa$values$Relative_eig[[3]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38d9b527",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_euclidian_pcoa_data,\n",
    "    x      = ~X,\n",
    "    y      = ~Y,\n",
    "    z      = ~Z,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    text   = bac_euclidian_pcoa_data$sampleid,\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Euclidian Distances [CSS]\") #%>% zgtools::plotly_elegant()\n",
    "                                    \n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "793b07cf",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f48e36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_euclidian_nmds_2D <- metaMDS(bac_feature_table_css,k=2,trymax=20,distance=\"euclidian\")\n",
    "bac_euclidian_nmds_3D <- metaMDS(bac_feature_table_css,k=3,trymax=20,distance=\"euclidian\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "285a8cb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_euclidian_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_euclidian_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_euclidian_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aa8c357",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_euclidian_nmds_3D_plot_data <- as.data.frame(bac_euclidian_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_euclidian_nmds_3D_plot_data$sampleid <- rownames(bac_euclidian_nmds_3D_plot_data)\n",
    "\n",
    "bac_euclidian_nmds_3D_plot_data <-\n",
    "  left_join(bac_euclidian_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "theme = list(tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_euclidian_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Euclidian Distances [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fafa2b62",
   "metadata": {},
   "source": [
    "\n",
    "<br><hr><br>\n",
    "\n",
    "<a id='Bray-CurtisDissimilarity'></a>\n",
    "\n",
    "# Bray-Curtis Dissimilarity\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The Bray-Curtis dissimilarity between two sites is defined as\n",
    "<br><br>\n",
    "$$BC(p,q) = 1-\\frac{2 C_{pq}}{S_p+S_q}$$\n",
    "<br>\n",
    "where $C_{pq}$ is the sum of the lesser values for only those species in common between both sites. ${S_p}$ and ${S_q}$ are the total number of specimens counted at both sites. \n",
    "<br><br>\n",
    "https://doi.org/10.2307/1942268\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf2f6a1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate distance matrix\n",
    "bac_braycurtis_dist <- vegan::vegdist(bac_feature_table, method = \"bray\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_braycurtis_pcoa <- ape::pcoa(bac_braycurtis_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_braycurtis_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_braycurtis_pcoa$vectors)),\n",
    "    X = bac_braycurtis_pcoa$vectors[, 1],\n",
    "    Y = bac_braycurtis_pcoa$vectors[, 2],\n",
    "    Z = bac_braycurtis_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_braycurtis_pcoa_data$Sample <- rownames(bac_braycurtis_pcoa_data)\n",
    "\n",
    "colnames(bac_braycurtis_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_braycurtis_pcoa_data <-\n",
    "  left_join(bac_braycurtis_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_braycurtis_pcoa$values$Relative_eig[[1]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_braycurtis_pcoa$values$Relative_eig[[2]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "                                    \n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_braycurtis_pcoa$values$Relative_eig[[3]] * 100, 2), \n",
    "                             \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9990ad90",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_braycurtis_pcoa_data,\n",
    "    x      = ~X,\n",
    "    y      = ~Y,\n",
    "    z      = ~Z,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    text   = bac_braycurtis_pcoa_data$sampleid,\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Bray-Curtis Dissimilarities\") #%>% zgtools::plotly_elegant()\n",
    "                                    \n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f40ae57",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a806a748",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_braycurtis_nmds_2D <- metaMDS(bac_feature_table,k=2,trymax=20,distance=\"bray\")\n",
    "bac_braycurtis_nmds_3D <- metaMDS(bac_feature_table,k=3,trymax=20,distance=\"bray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76ac107b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_braycurtis_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_braycurtis_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_braycurtis_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f78dd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_braycurtis_nmds_3D_plot_data <- as.data.frame(bac_braycurtis_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_braycurtis_nmds_3D_plot_data$sampleid <- rownames(bac_braycurtis_nmds_3D_plot_data)\n",
    "\n",
    "bac_braycurtis_nmds_3D_plot_data <-\n",
    "  left_join(bac_braycurtis_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_braycurtis_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Bray-Curtis Dissimilarities\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f75adcb",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d309a33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate distance matrix\n",
    "bac_braycurtis_dist <- vegan::vegdist(bac_feature_table_css, method = \"bray\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_braycurtis_pcoa <- ape::pcoa(bac_braycurtis_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_braycurtis_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_braycurtis_pcoa$vectors)),\n",
    "    X = bac_braycurtis_pcoa$vectors[, 1],\n",
    "    Y = bac_braycurtis_pcoa$vectors[, 2],\n",
    "    Z = bac_braycurtis_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_braycurtis_pcoa_data$Sample <- rownames(bac_braycurtis_pcoa_data)\n",
    "\n",
    "colnames(bac_braycurtis_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_braycurtis_pcoa_data <-\n",
    "  left_join(bac_braycurtis_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "# plotting\n",
    "axx <- list(title = paste(\"PCo 1 - \", round(\n",
    "  bac_braycurtis_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axy <- list(title = paste(\"PCo 2 - \", round(\n",
    "  bac_braycurtis_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axz <- list(title = paste(\"PCo 3 - \", round(\n",
    "  bac_braycurtis_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_braycurtis_pcoa_data,\n",
    "    x = ~X,\n",
    "    y = ~Y,\n",
    "    z = ~Z,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz),\n",
    "         title = \"\\nPCoA on Bray-Curtis Dissimilarities [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beb98d61",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cebeb280",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_braycurtis_nmds_2D <- metaMDS(bac_feature_table_css,k=2,trymax=20,distance=\"bray\")\n",
    "bac_braycurtis_nmds_3D <- metaMDS(bac_feature_table_css,k=3,trymax=20,distance=\"bray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8535c3a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_braycurtis_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_braycurtis_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_braycurtis_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce527a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_braycurtis_nmds_3D_plot_data <- as.data.frame(bac_braycurtis_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_braycurtis_nmds_3D_plot_data$sampleid <- rownames(bac_braycurtis_nmds_3D_plot_data)\n",
    "\n",
    "bac_braycurtis_nmds_3D_plot_data <-\n",
    "  left_join(bac_braycurtis_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_braycurtis_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Bray-Curtis Dissimilarities [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19147a06",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='JaccardDistance'></a>\n",
    "\n",
    "# Jaccard Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The Jaccard distance between two sites is defined as\n",
    "<br><br>\n",
    "$$d_j(P,Q) = \\frac{|P \\cup Q| - |P \\cap Q|}{P \\cap Q}$$\n",
    "<br>\n",
    "i.e. by dividing the difference of the sizes of the union and the intersection of two sets by the size of the union.\n",
    "<br><br>\n",
    "https://doi.org/10.1111/j.1469-8137.1912.tb05611.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e2384ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate distance matrix\n",
    "bac_jaccard_dist <- vegan::vegdist(bac_feature_table, method = \"jaccard\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_jaccard_pcoa <- ape::pcoa(bac_jaccard_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_jaccard_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_jaccard_pcoa$vectors)),\n",
    "    X = bac_jaccard_pcoa$vectors[, 1],\n",
    "    Y = bac_jaccard_pcoa$vectors[, 2],\n",
    "    Z = bac_jaccard_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_jaccard_pcoa_data$Sample <- rownames(bac_jaccard_pcoa_data)\n",
    "\n",
    "colnames(bac_jaccard_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_jaccard_pcoa_data <-\n",
    "  left_join(bac_jaccard_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title = paste(\"PCo 1 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axy <- list(title = paste(\"PCo 2 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axz <- list(title = paste(\"PCo 3 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c47a0a37",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_jaccard_pcoa_data,\n",
    "    x = ~X,\n",
    "    y = ~Y,\n",
    "    z = ~Z,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz),\n",
    "         title = \"\\nPCoA on Jaccard Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcc13d1d",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de54d23a",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_jaccard_nmds_2D <- metaMDS(bac_feature_table,k=2,trymax=20,distance=\"jaccard\")\n",
    "bac_jaccard_nmds_3D <- metaMDS(bac_feature_table,k=3,trymax=20,distance=\"jaccard\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad0bd02c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_jaccard_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_jaccard_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_jaccard_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "421d90c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_jaccard_nmds_3D_plot_data <- as.data.frame(bac_jaccard_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_jaccard_nmds_3D_plot_data$sampleid <- rownames(bac_jaccard_nmds_3D_plot_data)\n",
    "\n",
    "bac_jaccard_nmds_3D_plot_data <-\n",
    "  left_join(bac_jaccard_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_jaccard_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Jaccard Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46585a8a",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa63470f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate distance matrix\n",
    "bac_jaccard_dist <- vegan::vegdist(bac_feature_table_css, method = \"jaccard\")\n",
    "\n",
    "# compute the PCoA\n",
    "bac_jaccard_pcoa <- ape::pcoa(bac_jaccard_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_jaccard_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_jaccard_pcoa$vectors)),\n",
    "    X = bac_jaccard_pcoa$vectors[, 1],\n",
    "    Y = bac_jaccard_pcoa$vectors[, 2],\n",
    "    Z = bac_jaccard_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_jaccard_pcoa_data$Sample <- rownames(bac_jaccard_pcoa_data)\n",
    "\n",
    "colnames(bac_jaccard_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\",\"Z\")\n",
    "\n",
    "bac_jaccard_pcoa_data <-\n",
    "  left_join(bac_jaccard_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title = paste(\"PCo 1 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')),color = \"#aaaaaa\")\n",
    "axy <- list(title = paste(\"PCo 2 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axz <- list(title = paste(\"PCo 3 - \", round(\n",
    "  bac_jaccard_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43dc7fa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_jaccard_pcoa_data,\n",
    "    x = ~X,\n",
    "    y = ~Y,\n",
    "    z = ~Z,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz),\n",
    "         title = \"\\nPCoA on Jaccard Distances [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e629c537",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3281afba",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_jaccard_nmds_2D <- metaMDS(bac_feature_table_css,k=2,trymax=20,distance=\"jaccard\")\n",
    "bac_jaccard_nmds_3D <- metaMDS(bac_feature_table_css,k=3,trymax=20,distance=\"jaccard\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd9fd2c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_jaccard_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_jaccard_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_jaccard_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92ad5ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_jaccard_nmds_3D_plot_data <- as.data.frame(bac_jaccard_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_jaccard_nmds_3D_plot_data$sampleid <- rownames(bac_jaccard_nmds_3D_plot_data)\n",
    "\n",
    "bac_jaccard_nmds_3D_plot_data <-\n",
    "  left_join(bac_jaccard_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_jaccard_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Jaccard Distances [CSS]\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae996c2",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='UnweightedUniFracDistance'></a>\n",
    "\n",
    "# Unweighted UniFrac Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The unweighted UniFrac distance between two samples is defined as\n",
    "<br><br>\n",
    "$$ d^U=\\sum_{i=1}^{n} {\\frac{b_i|I(p_i^A>0)-I(p_i^B>0)|}{\\sum_{i=1}^{n}b_i}} $$\n",
    "<br>\n",
    "where $n$ is the number of branches in a rooted phylogenetic tree, $b_i$ is the length of branch $i$, $p_i^A$ and $p_i^B$ are the taxa proportions descending from the branch $i$ for community $A$ and $B$, respectively.  $I(.)$ is the indicator function and only presence/absence of species of branch $i$, $I(p_i^A > 0)$ and $I(p_i^B > 0)$.\n",
    "\n",
    "More generalized it can be described as\n",
    "<br><br>\n",
    "$$ {\\displaystyle \\left({\\frac {sum~of~unshared~branch~lengths}{sum~of~all~tree~branch~lengths}}\\right){=}fraction~of~total~unshared~branch~lengths} $$\n",
    "<br><br>\n",
    "https://doi.org/10.1093/bioinformatics/bts342"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eeec5ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.file <- file.path(\"/root/biommar/data/B_asv.txt\")\n",
    "tax.file <- file.path(\"/root/biommar/data/B_asv_taxonomy.txt\")\n",
    "meta.file <- file.path(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "bac_ps <- read_phyloseq(\n",
    "  otu.file      = otu.file,\n",
    "  taxonomy.file = tax.file,\n",
    "  metadata.file = meta.file,\n",
    "  type          = \"simple\",\n",
    "  sep           = \"\\t\"\n",
    ")\n",
    "\n",
    "# import tree\n",
    "tree <- read_tree(\"/root/biommar/data/tree.nwk\")\n",
    "\n",
    "bac_ps <- merge_phyloseq(bac_ps, tree)\n",
    "\n",
    "bac_u_unifrac_dist <- UniFrac(bac_ps, weighted = FALSE, normalized = TRUE,\n",
    "  parallel = FALSE, fast = TRUE)\n",
    "  \n",
    "# compute the PCoA\n",
    "bac_u_unifrac_pcoa <- ape::pcoa(bac_u_unifrac_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_u_unifrac_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_u_unifrac_pcoa$vectors)),\n",
    "    X = bac_u_unifrac_pcoa$vectors[, 1],\n",
    "    Y = bac_u_unifrac_pcoa$vectors[, 2],\n",
    "    Z = bac_u_unifrac_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_u_unifrac_pcoa_data$Sample <- rownames(bac_u_unifrac_pcoa_data)\n",
    "\n",
    "colnames(bac_u_unifrac_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_u_unifrac_pcoa_data <-\n",
    "  left_join(bac_u_unifrac_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title = paste(\"PCo 1 - \", round(\n",
    "  bac_u_unifrac_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')),color = \"#aaaaaa\")\n",
    "axy <- list(title = paste(\"PCo 2 - \", round(\n",
    "  bac_u_unifrac_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axz <- list(title = paste(\"PCo 3 - \", round(\n",
    "  bac_u_unifrac_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11f45cde",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_u_unifrac_pcoa_data,\n",
    "    x = ~X,\n",
    "    y = ~Y,\n",
    "    z = ~Z,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz),\n",
    "         title = \"\\nPCoA on Unweighted UniFrac Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d491a03",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd638a8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_2D <- metaMDS(bac_u_unifrac_dist,k=2,trymax=20)\n",
    "bac_u_unifrac_nmds_3D <- metaMDS(bac_u_unifrac_dist,k=3,trymax=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d801a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_u_unifrac_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_u_unifrac_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_u_unifrac_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b472e29d",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_3D_plot_data <- as.data.frame(bac_u_unifrac_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data$sampleid <- rownames(bac_u_unifrac_nmds_3D_plot_data)\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data <-\n",
    "  left_join(bac_u_unifrac_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_nmds_3D_plot_data,\n",
    "    x          = ~MDS1,\n",
    "    y          = ~MDS2,\n",
    "    z          = ~MDS3,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = theme,\n",
    "                      yaxis = theme,\n",
    "                      zaxis = theme),\n",
    "         title = \"\\nNMDS on Unweighted UniFrac Distances\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9ea8db",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49bcb483",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.file <- file.path(\"/root/biommar/data/CSS_B_asv.txt\")\n",
    "tax.file <- file.path(\"/root/biommar/data/CSS_B_asv_taxonomy.txt\")\n",
    "meta.file <- file.path(\"/root/biommar/data/CSS_B_asv_metadata.txt\")\n",
    "\n",
    "bac_ps_css <- read_phyloseq(\n",
    "  otu.file      = otu.file,\n",
    "  taxonomy.file = tax.file,\n",
    "  metadata.file = meta.file,\n",
    "  type          = \"simple\",\n",
    "  sep           = \"\\t\"\n",
    ")\n",
    "\n",
    "# import tree\n",
    "tree <- read_tree(\"/root/biommar/data/tree.nwk\")\n",
    "\n",
    "bac_ps_css <- merge_phyloseq(bac_ps_css, tree)\n",
    "\n",
    "bac_u_unifrac_dist <- UniFrac(bac_ps_css, weighted = FALSE, normalized = TRUE,\n",
    "  parallel = FALSE, fast = TRUE)\n",
    "\n",
    "# compute the PCoA\n",
    "bac_u_unifrac_pcoa <- ape::pcoa(bac_u_unifrac_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_u_unifrac_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_u_unifrac_pcoa$vectors)),\n",
    "    X = bac_u_unifrac_pcoa$vectors[, 1],\n",
    "    Y = bac_u_unifrac_pcoa$vectors[, 2],\n",
    "    Z = bac_u_unifrac_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_u_unifrac_pcoa_data$Sample <- rownames(bac_u_unifrac_pcoa_data)\n",
    "\n",
    "colnames(bac_u_unifrac_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_u_unifrac_pcoa_data <- left_join(bac_u_unifrac_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9738ce60",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_pcoa_data,\n",
    "    x          = ~X,\n",
    "    y          = ~Y,\n",
    "    z          = ~Z,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Unweighted UniFrac Distances [CSS]\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74488d0f",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ab9f95",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_2D <- metaMDS(bac_u_unifrac_dist, k = 2, trymax = 20)\n",
    "bac_u_unifrac_nmds_3D <- metaMDS(bac_u_unifrac_dist, k = 3, trymax = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4008cf1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_u_unifrac_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_u_unifrac_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_u_unifrac_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74cb9fd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_3D_plot_data <- as.data.frame(bac_u_unifrac_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data$sampleid <- rownames(bac_u_unifrac_nmds_3D_plot_data)\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data <-left_join(bac_u_unifrac_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_nmds_3D_plot_data,\n",
    "    x          = ~MDS1,\n",
    "    y          = ~MDS2,\n",
    "    z          = ~MDS3,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = theme,\n",
    "                      yaxis = theme,\n",
    "                      zaxis = theme),\n",
    "         title = \"\\nNMDS on Unweighted UniFrac Distances [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0b8a6f6",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='WeightedUniFracDistance'></a>\n",
    "\n",
    "# Weighted UniFrac Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The weighted UniFrac distance between two samples $(A,B)$ is defined as\n",
    "<br><br>\n",
    "$$ d^W={\\frac{\\sum_{i=1}^{n}b_i|p_i^A-p_i^B|}{\\sum_{i=1}^{n}b_i(p_i^A+p_i^B)}} $$\n",
    "<br><br>\n",
    "https://doi.org/10.1093/bioinformatics/bts342"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37418611",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.file  <- file.path(\"/root/biommar/data/B_asv.txt\")\n",
    "tax.file  <- file.path(\"/root/biommar/data/B_asv_taxonomy.txt\")\n",
    "meta.file <- file.path(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "bac_ps <- read_phyloseq(\n",
    "  otu.file      = otu.file,\n",
    "  taxonomy.file = tax.file,\n",
    "  metadata.file = meta.file,\n",
    "  type          = \"simple\",\n",
    "  sep           = \"\\t\"\n",
    ")\n",
    "\n",
    "# import tree\n",
    "tree <- read_tree(\"/root/biommar/data/tree.nwk\")\n",
    "\n",
    "bac_ps <- merge_phyloseq(bac_ps, tree)\n",
    "\n",
    "bac_u_unifrac_dist <- UniFrac(bac_ps, weighted = TRUE, normalized = TRUE, parallel = FALSE, fast = TRUE)\n",
    "  \n",
    "# compute the PCoA\n",
    "bac_u_unifrac_pcoa <- ape::pcoa(bac_u_unifrac_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_u_unifrac_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_u_unifrac_pcoa$vectors)),\n",
    "    X      = bac_u_unifrac_pcoa$vectors[, 1],\n",
    "    Y      = bac_u_unifrac_pcoa$vectors[, 2],\n",
    "    Z      = bac_u_unifrac_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_u_unifrac_pcoa_data$Sample <- rownames(bac_u_unifrac_pcoa_data)\n",
    "\n",
    "colnames(bac_u_unifrac_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_u_unifrac_pcoa_data <- left_join(bac_u_unifrac_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_u_unifrac_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960f1d71",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_pcoa_data,\n",
    "    x          = ~X,\n",
    "    y          = ~Y,\n",
    "    z          = ~Z,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Weighted UniFrac Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa2b2aec",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15938e68",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_2D <- metaMDS(bac_u_unifrac_dist, k = 2, trymax = 20)\n",
    "bac_u_unifrac_nmds_3D <- metaMDS(bac_u_unifrac_dist, k = 3, trymax = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "874562f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_u_unifrac_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_u_unifrac_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_u_unifrac_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4774b5c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_3D_plot_data <- as.data.frame(bac_u_unifrac_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data$sampleid <- rownames(bac_u_unifrac_nmds_3D_plot_data)\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data <- left_join(bac_u_unifrac_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_nmds_3D_plot_data,\n",
    "    x          = ~MDS1,\n",
    "    y          = ~MDS2,\n",
    "    z          = ~MDS3,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = theme,\n",
    "                      yaxis = theme,\n",
    "                      zaxis = theme),\n",
    "         title = \"NMDS on Weighted UniFrac Distances\") %>% \n",
    "  zgtools::plotly_elegant()\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85f48029",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Cumulative Sum Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53645bcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_w_unifrac_dist <- UniFrac(bac_ps_css, weighted = TRUE, normalized = TRUE, parallel = FALSE, fast = TRUE)\n",
    "\n",
    "# compute the PCoA\n",
    "bac_w_unifrac_pcoa <- ape::pcoa(bac_w_unifrac_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_w_unifrac_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_w_unifrac_pcoa$vectors)),\n",
    "    X      = bac_w_unifrac_pcoa$vectors[, 1],\n",
    "    Y      = bac_w_unifrac_pcoa$vectors[, 2],\n",
    "    Z      = bac_w_unifrac_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_w_unifrac_pcoa_data$Sample <- rownames(bac_w_unifrac_pcoa_data)\n",
    "\n",
    "colnames(bac_w_unifrac_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_w_unifrac_pcoa_data <- left_join(bac_w_unifrac_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "# plotting\n",
    "axx <- list(title    = paste(\"PCo 1 - \", round( bac_w_unifrac_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axy <- list(title    = paste(\"PCo 2 - \", round(bac_w_unifrac_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axz <- list(title    = paste(\"PCo 3 - \", round(bac_w_unifrac_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_w_unifrac_pcoa_data,\n",
    "    x          = ~X,\n",
    "    y          = ~Y,\n",
    "    z          = ~Z,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"PCoA on Weighted UniFrac Distances [CSS]\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aeb43c8",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c67dd4f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_2D <- metaMDS(bac_w_unifrac_dist,k=2,trymax=20)\n",
    "bac_u_unifrac_nmds_3D <- metaMDS(bac_w_unifrac_dist,k=3,trymax=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1425d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_u_unifrac_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_u_unifrac_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_u_unifrac_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5036d889",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_u_unifrac_nmds_3D_plot_data <- as.data.frame(bac_u_unifrac_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data$sampleid <- rownames(bac_u_unifrac_nmds_3D_plot_data)\n",
    "\n",
    "bac_u_unifrac_nmds_3D_plot_data <-\n",
    "  left_join(bac_u_unifrac_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_u_unifrac_nmds_3D_plot_data,\n",
    "    x          = ~MDS1,\n",
    "    y          = ~MDS2,\n",
    "    z          = ~MDS3,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = theme,\n",
    "                      yaxis = theme,\n",
    "                      zaxis = theme),\n",
    "         title = \"NMDS on Weighted UniFrac Distances [CSS]\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1c0c0a",
   "metadata": {},
   "source": [
    "### Remark:<br>\n",
    "\n",
    "Weighted UniFrac distance incorporates the branch lengths as absolute positive values. The phylogenetic tree on which it operates is produced in *Qiime2* with the </span><span class=\"code\">fragment-insertion</span><span class=\"x\"> plugin. Internally, SEPP places the reads on a phylogenetic reference tree (</span><span class=\"code\">silva128</span><span class=\"x\">). Unknown clades may be placed with considerable branch lengths and therefore dominate feature segregation in weighted UniFrac plots.\n",
    "\n",
    "https://doi.org/10.1128/mSystems.00021-18"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2499815",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='AitchisonDistance'></a>\n",
    "\n",
    "# Aitchison Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The Aitchison distance between two $n$-dimensional compositions, $x_i$ and $x_j$ is defined as\n",
    "<br><br>\n",
    "$$ d^A(x_i,x_j)=\\sqrt{\\sum_{g=1}^n\\displaystyle \\left[ln\\frac{x_{gi}}{g(x_i)}-ln\\frac{x_{gj}}{g(x_j)}\\right]^2} $$\n",
    "<br>\n",
    "where $g(.)$ is the geometric mean. The Aitchison distance is equal to the Euclidean distance between clr-transformed compositions. It has scale invariance, perturbation invariance, permutation invariance and sub-compositional dominance.\n",
    "\n",
    "https://doi.org/10.1093/bioinformatics/bty175"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7a51da",
   "metadata": {},
   "outputs": [],
   "source": [
    "aitchison_data <- bac_feature_table\n",
    "\n",
    "# add pseudo count \n",
    "aitchison_data[aitchison_data == 0] <- 1e-12\n",
    "\n",
    "aitchison_data <- t(apply(aitchison_data, 1, \n",
    "                          function(x)log(x/exp(mean(log(x))))))\n",
    "bac_aitchison_dist <- as.matrix(dist(aitchison_data))\n",
    "\n",
    "# compute the PCoA\n",
    "bac_aitchison_pcoa <- ape::pcoa(bac_aitchison_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_aitchison_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_aitchison_pcoa$vectors)),\n",
    "    X = bac_aitchison_pcoa$vectors[, 1],\n",
    "    Y = bac_aitchison_pcoa$vectors[, 2],\n",
    "    Z = bac_aitchison_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_aitchison_pcoa_data$Sample <- rownames(bac_aitchison_pcoa_data)\n",
    "\n",
    "colnames(bac_aitchison_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_aitchison_pcoa_data <-\n",
    "  left_join(bac_aitchison_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "# plotting\n",
    "axx <- list(title = paste(\"PCo 1 - \", round(\n",
    "  bac_aitchison_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axy <- list(title = paste(\"PCo 2 - \", round(\n",
    "  bac_aitchison_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "axz <- list(title = paste(\"PCo 3 - \", round(\n",
    "  bac_aitchison_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "  tickfont=list(color=c('#323232')), color = \"#aaaaaa\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data       = bac_aitchison_pcoa_data,\n",
    "    x          = ~X,\n",
    "    y          = ~Y,\n",
    "    z          = ~Z,\n",
    "    color      = ~department,\n",
    "    colors     = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker     = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"PCoA on Aitchison Distances\")\n",
    "                          \n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e808dd78",
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_aitchison_pcoa_data,\n",
    "    x      = ~X,\n",
    "    y      = ~Y,\n",
    "    z      = ~Z,\n",
    "    color  = ~Clade,\n",
    "    colors = c(\"dodgerblue2\", \"#E31A1C\", \"green4\", \"#6A3D9A\", \n",
    "               \"#FF7F00\", \"gold1\", \"maroon\", \"orchid1\", \"deeppink1\", \n",
    "               \"blue1\", \"steelblue4\", \"darkturquoise\", \"green1\", \n",
    "               \"yellow4\",\"darkgreen\",\"black\", \"thistle4\"),\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"PCoA on Aitchison Distances\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae14d2a3",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Non-Metric Multidimensional Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "508411db",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_aitchison_nmds_2D <- metaMDS(bac_aitchison_dist,k=2,trymax=20)\n",
    "bac_aitchison_nmds_3D <- metaMDS(bac_aitchison_dist,k=3,trymax=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec5f993",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_aitchison_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_aitchison_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_aitchison_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73adb726",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_aitchison_nmds_3D_plot_data <- as.data.frame(bac_aitchison_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_aitchison_nmds_3D_plot_data$sampleid <- rownames(bac_aitchison_nmds_3D_plot_data)\n",
    "\n",
    "bac_aitchison_nmds_3D_plot_data <-\n",
    "  left_join(bac_aitchison_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_aitchison_nmds_3D_plot_data,\n",
    "    x      = ~MDS1,\n",
    "    y      = ~MDS2,\n",
    "    z      = ~MDS3,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = theme,\n",
    "                      yaxis = theme,\n",
    "                      zaxis = theme),\n",
    "         title = \"NMDS on Aitchison Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdd55ca2",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='PhILRDistance'></a>\n",
    "\n",
    "# PhILR Distance\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "\n",
    "The **Ph**ylogenetic **I**sometric **L**og-**R**atio transformation combines evolutionary models with the isometric log-ratio transform and projects the data into a cartesian space.\n",
    "\n",
    "https://doi.org/10.7554/eLife.21887.001\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dac958dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.file <- file.path(\"/root/biommar/data/B_asv.txt\")\n",
    "tax.file <- file.path(\"/root/biommar/data/B_asv_taxonomy.txt\")\n",
    "meta.file <- file.path(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "bac_ps <- read_phyloseq(\n",
    "  otu.file = otu.file,\n",
    "  taxonomy.file = \"/root/biommar/data/taxa.txt\",\n",
    "  metadata.file = meta.file,\n",
    "  type = \"simple\",\n",
    "  sep = \"\\t\"\n",
    ")\n",
    "\n",
    "meta <- read.table(meta.file, row.names = 1, header = T, sep=\"\\t\")\n",
    "\n",
    "## import tree\n",
    "tree <- read_tree(\"/root/biommar/data/tree.nwk\")\n",
    "\n",
    "bac_ps <- merge_phyloseq(bac_ps, tree)\n",
    "\n",
    "tse <- mia::makeTreeSummarizedExperimentFromPhyloseq(bac_ps)\n",
    "\n",
    "tse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cb82f57",
   "metadata": {},
   "outputs": [],
   "source": [
    "#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "#\n",
    "# Filter Extremely Low-Abundance OTUs\n",
    "#\n",
    "# Taxa that were not seen with more than 3 counts in at least 20% of samples are\n",
    "# filtered. Subsequently, those with a coefficient of variation ≤ 3 are\n",
    "# filtered. Finally we add a pseudocount of 1 to the remaining OTUs to avoid\n",
    "# calculating log-ratios involving zeros.\n",
    "\n",
    "## Select prevalent taxa \n",
    "tse <-  tse %>% subsetByPrevalentTaxa(\n",
    "  detection = 0.1,\n",
    "  prevalence = 0.001,\n",
    "  as_relative = FALSE)\n",
    "\n",
    "## Pick taxa that have notable abundance variation across sammples\n",
    "variable.taxa <- apply(assays(tse)$counts, 1, function(x) sd(x)/mean(x) > 2.0)\n",
    "tse <- tse[variable.taxa,]\n",
    "\n",
    "tse\n",
    "\n",
    "# Collapse the tree!\n",
    "# Otherwise the original tree with all nodes is kept\n",
    "# (including those that were filtered out from rowData)\n",
    "tree <- ape::keep.tip(phy = rowTree(tse), tip = rowLinks(tse)$nodeNum)\n",
    "rowTree(tse) <- tree\n",
    "\n",
    "## Add a new assay with a pseudocount \n",
    "assays(tse)$counts.shifted <- assays(tse)$counts + 1e-12\n",
    "\n",
    "tree <- makeNodeLabel(tree, method=\"number\", prefix='n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a05fe5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "is.rooted(tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0232ff7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "is.binary(tree)\n",
    "# if FALSE\n",
    "# tree <- ape::multi2di(tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fb0c920",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the modified tree back to the (`TreeSE`) data object \n",
    "rowTree(tse) <- tree\n",
    "\n",
    "# Extract taxonomy table from the TreeSE object\n",
    "tax <- rowData(tse)[,taxonomyRanks(tse)]\n",
    "\n",
    "# Get name balances\n",
    "name.balance(tree, tax, 'n1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2252fe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "otu.table <- t(as(assays(tse)$counts.shifted, \"matrix\"))\n",
    "tree <- rowTree(tse)\n",
    "metadata <- colData(tse)\n",
    "tax <- rowData(tse)[,taxonomyRanks(tse)]\n",
    "\n",
    "data_philr <- philr(tse, \n",
    "                    part.weights = 'uniform',\n",
    "                    ilr.weights  = 'uniform', \n",
    "                    return.all   = FALSE,\n",
    "                    abund_values = \"counts.shifted\")\n",
    "\n",
    "# Distances between samples based on philr transformed data\n",
    "bac_philr_dist <- dist(data_philr, method = \"euclidean\") \n",
    "\n",
    "# compute the PCoA\n",
    "bac_philr_pcoa <- ape::pcoa(data_dist)\n",
    "\n",
    "# save PCo vector data for plot\n",
    "bac_philr_pcoa_data <-\n",
    "  data.frame(\n",
    "    Sample = as.character(rownames(bac_philr_pcoa$vectors)),\n",
    "    X = bac_philr_pcoa$vectors[, 1],\n",
    "    Y = bac_philr_pcoa$vectors[, 2],\n",
    "    Z = bac_philr_pcoa$vectors[, 3]\n",
    "  )\n",
    "\n",
    "# add metadata\n",
    "bac_philr_pcoa_data$Sample <- rownames(bac_philr_pcoa_data)\n",
    "\n",
    "colnames(bac_philr_pcoa_data) <- c(\"sampleid\", \"X\", \"Y\", \"Z\")\n",
    "\n",
    "bac_philr_pcoa_data <- left_join(bac_philr_pcoa_data, bac_metadata, by = \"sampleid\")\n",
    "  \n",
    "# plotting\n",
    "axx <- list(title    = paste(\"PCo 1 - \", \n",
    "                             round(bac_philr_pcoa$values$Relative_eig[[1]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axy <- list(title    = paste(\"PCo 2 - \", \n",
    "                             round(bac_philr_pcoa$values$Relative_eig[[2]] * 100, 2),\"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "axz <- list(title    = paste(\"PCo 3 - \", \n",
    "                             round(bac_philr_pcoa$values$Relative_eig[[3]] * 100, 2), \"%\", sep = \"\"),\n",
    "            tickfont = list(color = c('#323232')), \n",
    "            color    = \"#aaaaaa\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data   = bac_philr_pcoa_data,\n",
    "    x      = ~X,\n",
    "    y      = ~Y,\n",
    "    z      = ~Z,\n",
    "    color  = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(size = 5, opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis = axx,\n",
    "                      yaxis = axy,\n",
    "                      zaxis = axz),\n",
    "         title = \"\\nPCoA on Euclidian Distances\") \n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55283382",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_philr_nmds_2D <- metaMDS(bac_philr_dist, k = 2, trymax = 20)\n",
    "bac_philr_nmds_3D <- metaMDS(bac_philr_dist, k = 3, trymax = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5d262a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat(\"\\n\", \"Stress 3D: \", bac_philr_nmds_3D[[\"stress\"]], \"\\n\",\n",
    "    \"Stress 2D: \", bac_philr_nmds_2D[[\"stress\"]])\n",
    "\n",
    "stressplot(bac_philr_nmds_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f5d2666",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_philr_nmds_3D_plot_data <- as.data.frame(bac_philr_nmds_3D[[\"points\"]])\n",
    "\n",
    "bac_philr_nmds_3D_plot_data$sampleid <- rownames(bac_philr_nmds_3D_plot_data)\n",
    "\n",
    "bac_philr_nmds_3D_plot_data <-\n",
    "  left_join(bac_philr_nmds_3D_plot_data, bac_metadata, by = \"sampleid\")\n",
    "\n",
    "p <- plot_ly(width = 900, height = 700) %>%\n",
    "  add_trace(\n",
    "    type='scatter3d',\n",
    "    mode='markers',\n",
    "    data = bac_philr_nmds_3D_plot_data,\n",
    "    x = ~MDS1,\n",
    "    y = ~MDS2,\n",
    "    z = ~MDS3,\n",
    "    color = ~department,\n",
    "    colors = c(\"#00468B\", \"#ED0000\", \"#42B540\", \"#0099B4\", \"#925E9F\"),\n",
    "    marker = list(\n",
    "      size = 5,\n",
    "      opacity = 1.0),\n",
    "    showlegend = T) %>% \n",
    "  layout(scene = list(xaxis=theme,yaxis=theme,zaxis=theme),\n",
    "         title = \"\\nNMDS on Euclidian Distances\")\n",
    "\n",
    "embed_notebook(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a60986b",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='Statistics'></a>\n",
    "\n",
    "# Statistics\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d52241d1",
   "metadata": {},
   "source": [
    "### Pairwise ADONIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69d9d663",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://github.com/pmartinezarbizu/pairwiseAdonis/blob/master/pairwiseAdonis/R/pairwise.adonis.R\n",
    "\n",
    "pairwise.adonis <- function(x, factors, sim.function = 'vegdist', sim.method = 'bray', \n",
    "                            p.adjust.m ='bonferroni', reduce=NULL, perm=999){\n",
    "  require(\"vegan\")\n",
    "  co <- combn(unique(as.character(factors)),2)\n",
    "  pairs <- c()\n",
    "  Df <- c()\n",
    "  SumsOfSqs <- c()\n",
    "  F.Model <- c()\n",
    "  R2 <- c()\n",
    "  p.value <- c()\n",
    "  for(elem in 1:ncol(co)){\n",
    "    if(inherits(x, 'dist')){\n",
    "      x1=as.matrix(x)[factors %in% c(as.character(co[1,elem]), as.character(co[2,elem])),\n",
    "                      factors %in% c(as.character(co[1,elem]), as.character(co[2,elem]))]\n",
    "    }\n",
    "    else  (\n",
    "      if (sim.function == 'daisy'){\n",
    "        x1 = daisy(x[factors %in% c(co[1,elem],co[2,elem]),],metric=sim.method)\n",
    "      } \n",
    "      else{x1 = vegdist(x[factors %in% c(co[1,elem],co[2,elem]),], method=sim.method)}\n",
    "    )\n",
    "    x2 = data.frame(Fac = factors[factors %in% c(co[1,elem],co[2,elem])])\n",
    "    ad <- adonis2(x1 ~ Fac, data = x2, permutations = perm);\n",
    "    pairs <- c(pairs,paste(co[1,elem],'-',co[2,elem]));\n",
    "    Df <- c(Df,ad$Df[1])\n",
    "    SumsOfSqs <- c(SumsOfSqs,ad$SumOfSqs[1])\n",
    "    F.Model <- c(F.Model,ad$F[1]);\n",
    "    R2 <- c(R2,ad$R2[1]);\n",
    "    p.value <- c(p.value,ad$`Pr(>F)`[1])\n",
    "  }\n",
    "  p.adjusted <- p.adjust(p.value,method=p.adjust.m)\n",
    "  sig = c(rep('',length(p.adjusted)))\n",
    "  sig[p.adjusted <= 0.1] <-'.'\n",
    "  sig[p.adjusted <= 0.05] <-'* '\n",
    "  sig[p.adjusted <= 0.01] <-'**'\n",
    "  sig[p.adjusted <= 0.001] <-'***'\n",
    "  pairw.res <- data.frame(pairs,Df,SumsOfSqs,F.Model,R2,p.value,p.adjusted,sig)\n",
    "  if(!is.null(reduce)){\n",
    "    pairw.res <- subset (pairw.res, grepl(reduce,pairs))\n",
    "    pairw.res$p.adjusted <- p.adjust(pairw.res$p.value,method=p.adjust.m)\n",
    "    sig = c(rep('',length(pairw.res$p.adjusted)))\n",
    "    sig[pairw.res$p.adjusted <= 0.1] <-'.'\n",
    "    sig[pairw.res$p.adjusted <= 0.05] <-'* '\n",
    "    sig[pairw.res$p.adjusted <= 0.01] <-'**'\n",
    "    sig[pairw.res$p.adjusted <= 0.001] <-'***'\n",
    "    pairw.res <- data.frame(pairw.res[,1:7],sig)\n",
    "  }\n",
    "  class(pairw.res) <- c(\"pwadonis\", \"data.frame\")\n",
    "  return(pairw.res)\n",
    "} \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83e37371",
   "metadata": {},
   "outputs": [],
   "source": [
    "philr.df <- as.data.frame(data_philr)\n",
    "meta <- read.table(meta.file, row.names = 1, header = T, sep=\"\\t\")\n",
    "#meta$country[is.na(meta$country)] <- \"Curazao\"\n",
    "\n",
    "philr.dist <- vegdist(philr.df, \"euclidean\")\n",
    "\n",
    "disp.country <- permutest(betadisper(philr.dist, meta$country), \n",
    "                          pairwise = TRUE, permutations = how(nperm = 999))\n",
    "\n",
    "disp.Clade <- permutest(betadisper(philr.dist, meta$Clade), \n",
    "                        pairwise = TRUE, permutations = how(nperm = 999))\n",
    "\n",
    "disp.locality <- permutest(betadisper(philr.dist, meta$locality), \n",
    "                           pairwise = TRUE, permutations = how(nperm = 999))\n",
    "\n",
    "disp.department <- permutest(betadisper(philr.dist, meta$department), \n",
    "                           pairwise = TRUE, permutations = how(nperm = 999))\n",
    "\n",
    "philr.dist <- as.data.frame(as.matrix(vegdist(philr.df, \"euclidean\")))\n",
    "\n",
    "x1 <- capture.output(print(disp.country), file = NULL)\n",
    "x2 <- capture.output(print(disp.Clade), file = NULL)\n",
    "x3 <- capture.output(print(disp.locality), file = NULL)\n",
    "x4 <- capture.output(print(disp.department), file = NULL)\n",
    "\n",
    "gsub(\" \\\\(t\\\\)\", \"\" , names(disp.country[[\"statistic\"]][2:length(disp.country[[\"statistic\"]])]))\n",
    "\n",
    "#paste(x, collapse = '\\n') %>% cat()\n",
    "\n",
    "pwa.country <- pairwise.adonis(data_philr , meta$country, \n",
    "                               p.adjust.m = \"bonferroni\",\n",
    "                               sim.method = \"euclidean\", perm = 999)\n",
    "\n",
    "pwa.Clade <- pairwise.adonis(data_philr , meta$Clade, \n",
    "                             p.adjust.m = \"bonferroni\",\n",
    "                             sim.method = \"euclidean\", perm = 999)\n",
    "\n",
    "pwa.locality <- pairwise.adonis(data_philr , meta$locality, \n",
    "                                p.adjust.m = \"bonferroni\",\n",
    "                                sim.method = \"euclidean\", perm = 999)\n",
    "\n",
    "pwa.department <- pairwise.adonis(data_philr , meta$department, \n",
    "                                  p.adjust.m = \"bonferroni\",\n",
    "                                  sim.method = \"euclidean\", perm = 999)\n",
    "\n",
    "rownames(pwa.country) <- names(disp.country[[\"pairwise\"]][[\"permuted\"]])\n",
    "\n",
    "\n",
    "pwa.country$disp <- disp.country[[\"pairwise\"]][[\"permuted\"]]\n",
    "pwa.country$disp.sig <- \n",
    "  ifelse (pwa.country$disp <= 0.001, \"***\",\n",
    "          ifelse (pwa.country$disp <= 0.01, \"**\",\n",
    "                  ifelse (pwa.country$disp <= 0.05, \"*\",\n",
    "                          ifelse(pwa.country$disp <= 0.1, \".\", \"\"))))\n",
    "\n",
    "pwa.Clade$disp <- disp.Clade[[\"pairwise\"]][[\"permuted\"]]\n",
    "pwa.Clade$disp.sig <- \n",
    "  ifelse (pwa.Clade$disp <= 0.001, \"***\",\n",
    "          ifelse (pwa.Clade$disp <= 0.01, \"**\",\n",
    "                  ifelse (pwa.Clade$disp <= 0.05, \"*\",\n",
    "                          ifelse(pwa.Clade$disp <= 0.1, \".\", \"\"))))\n",
    "\n",
    "pwa.locality$disp <- disp.locality[[\"pairwise\"]][[\"permuted\"]]\n",
    "pwa.locality$disp.sig <- \n",
    "  ifelse (pwa.locality$disp <= 0.001, \"***\",\n",
    "          ifelse (pwa.locality$disp <= 0.01, \"**\",\n",
    "                  ifelse (pwa.locality$disp <= 0.05, \"*\",\n",
    "                          ifelse(pwa.locality$disp <= 0.1, \".\", \"\"))))\n",
    "\n",
    "\n",
    "pwa.department$disp <- disp.department[[\"pairwise\"]][[\"permuted\"]]\n",
    "pwa.department$disp.sig <- \n",
    "  ifelse (pwa.department$disp <= 0.001, \"***\",\n",
    "          ifelse (pwa.department$disp <= 0.01, \"**\",\n",
    "                  ifelse (pwa.department$disp <= 0.05, \"*\",\n",
    "                          ifelse(pwa.department$disp <= 0.1, \".\", \"\"))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60256737",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwa.country %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c77f76d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwa.Clade %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac24f9f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwa.locality %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf91165",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwa.department %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb3782c2",
   "metadata": {},
   "source": [
    "### Pairwise ANOSIM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91583b15",
   "metadata": {},
   "outputs": [],
   "source": [
    "pairwise_anosim <- function(x, meta_data, group, \n",
    "                            permutations = 999, distance = \"bray\")\n",
    "{\n",
    "  .group <- unique(meta_data[, group])\n",
    "  .combs <- combn(.group, 2)\n",
    "  .res <- data.frame(matrix(ncol = 4, nrow = length(.combs)/2))\n",
    "  require(\"vegan\")\n",
    "  cat(paste(\"Pairwise ANOSIM of\", length(.combs)/2, \"combinations.\\n   \n",
    "              Permutations: \", permutations))\n",
    "  for(i in seq(length(.combs)/2)) {\n",
    "    .this_data <- x[meta_data[, group] %in% .combs[,i], ]\n",
    "    .this_meta <- meta_data[meta_data[, group] %in% .combs[,i], ]\n",
    "    .tmp <- anosim(.this_data, .this_meta[, group], \n",
    "                   permutations = permutations, distance = distance)\n",
    "    .var <- as.character(unique(.tmp[[\"class.vec\"]]))\n",
    "    .sig <- as.character(.tmp[[\"signif\"]])\n",
    "    .R <- as.character(.tmp[[\"statistic\"]])\n",
    "    .res[i,] <- c(.var[-which(.var==\"Between\")], .sig, .R) \n",
    "    cat(paste(\"\\r    \",i, \"of\", length(.combs)/2, \"done\"), fill = TRUE)\n",
    "  }\n",
    "  colnames(.res) <- c(\"variable.1\", \"variable.2\", \"p.value\", \"R.value\")\n",
    "  .res$p.value <- as.numeric(.res$p.value)\n",
    "  .res <- .res[order(.res$p.value), ]\n",
    "  .res$adj.p <- .res$p.value * rev(seq(1:nrow(.res)))\n",
    "  .dis <- vegdist(x, method = distance)\n",
    "  .bet <- permutest(betadisper(.dis, meta_data[,group]), pairwise = TRUE)\n",
    "  .p <- .bet[[\"tab\"]][[\"Pr(>F)\"]][1]\n",
    "  .res$dispersion.p <- .p\n",
    "  return(.res)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6bb641",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.country <- pairwise_anosim(data_philr , meta, \"country\", \n",
    "                               distance = \"euclidean\",\n",
    "                               permutations = 999)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "765aa0a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.country %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "027e5df6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.department <- pairwise_anosim(data_philr , meta, \"department\", \n",
    "                                  distance = \"euclidean\", \n",
    "                                  permutations = 999)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddacf6af",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.department %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa520410",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.Clade <- pairwise_anosim(data_philr , meta, \"Clade\", \n",
    "                                  distance = \"euclidean\", \n",
    "                                  permutations = 999)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dffd754",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.Clade %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88687d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.morpho_specie <- pairwise_anosim(data_philr , meta, \"morpho_specie\", \n",
    "                                  distance = \"euclidean\", \n",
    "                                  permutations = 999)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9ab9187",
   "metadata": {},
   "outputs": [],
   "source": [
    "pwn.morpho_specie %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6802fa00",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a46644d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pairwise_anosim_plot <- function(x, meta_data, group, \n",
    "                                 permutations = 999, distance = \"bray\")\n",
    "{\n",
    "  require(\"ggplot2\")\n",
    "  require(\"reshape\")\n",
    "  require(\"vegan\")\n",
    "  require(\"multcompView\")\n",
    "  require(\"dplyr\")\n",
    "  .pdat <- pairwise_anosim(x, meta_data, group, \n",
    "                           permutations = permutations, distance = distance)\n",
    "  .beta <- unique(.pdat[,ncol(.pdat)])\n",
    "  .pdat <- .pdat[,-ncol(.pdat)]\n",
    "  .tmp <- anosim(x, meta_data[, group], permutations = permutations, distance = distance)\n",
    "  # extract ranks to plot on y axis and respective variable to plot on x axis\n",
    "  .xdat <- data.frame(t(as.data.frame(rbind(\n",
    "    as.character(.tmp[[\"class.vec\"]]), .tmp[[\"dis.rank\"]]))))\n",
    "  .xdat$X2 <- as.numeric(.xdat$X2)\n",
    "  # make symmetric matrix of P-values for compact letter display\n",
    "  .x1 <- reshape::cast(.pdat, variable.1 ~ variable.2, value = \"adj.p\", fill = 0)\n",
    "  .x2 <- reshape::cast(.pdat, variable.2 ~ variable.1, value = \"adj.p\", fill = 0)\n",
    "  rownames(.x1) <- .x1$variable.1 \n",
    "  rownames(.x2) <- .x2$variable.2 \n",
    "  # casting produces X*X-1 matrix; solve by merging X*X-1 and X-1*X matrices\n",
    "  .z1 <- base::merge(.x1, .x2, by = \"row.names\", all = T)\n",
    "  rownames(.z1) <- .z1$Row.names\n",
    "  .z1 <- .z1[, 2:ncol(.x1) + 1]\n",
    "  .z2 <- merge(.x2, .x1, by =\"row.names\", all = T)\n",
    "  rownames(.z2) <- .z2$Row.names\n",
    "  .z2 <- .z2[, 2:ncol(.x1) + 1]\n",
    "  .z <- merge(.z1, .z2, by = \"row.names\", all = T)\n",
    "  rownames(.z) <- .z$Row.names\n",
    "  .z <- .z[ , -1]\n",
    "  .z <- t(.z)\n",
    "  # get rid of duplicate names\n",
    "  .out <- .z[grep(\".x.y\", row.names(.z)), ]\n",
    "  .out <- rownames(.out)\n",
    "  .z <- subset(.z, !(rownames(.z) %in% .out))\n",
    "  .z <- as.data.frame(.z)\n",
    "  rownames(.z) <- gsub(\".x.x$\", \"\", rownames(.z))\n",
    "  .z <- .z[order(rownames(.z)), ]\n",
    "  .z <- .z[order(colnames(.z)), ]\n",
    "  diag(.z) <- 1\n",
    "  .z[.z == 0] <- NA\n",
    "  # fill empty cells with diagonal pair\n",
    "  for (i in 1:ncol(.z)){\n",
    "    for (j in 1:nrow(.z)){\n",
    "      if (!is.na(.z[j,i])){\n",
    "        .z[i,j] <- .z[j,i]\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "  .let <- multcompView::multcompLetters(.z)\n",
    "  .let <- as.data.frame(.let[[\"Letters\"]])\n",
    "  .let$X1 <- rownames(.let)\n",
    "  .plt <- merge(.xdat, .let, by =\"X1\", all = T)\n",
    "  .p <- ggplot(.plt, aes(x = X1, y = X2, fill = X1)) + \n",
    "    geom_jitter(aes(color = X1), alpha = 0.2, size = 2) +\n",
    "    geom_boxplot(position = position_dodge2(),\n",
    "                 alpha = 0.9, varwidth = TRUE) +  \n",
    "    geom_label(data = .plt %>% \n",
    "               group_by(X1) %>% \n",
    "               summarise(y = max(X2), l = `.let[[\"Letters\"]]`),\n",
    "               aes(y = y, label = l), position = position_dodge(width = 1), \n",
    "               size = 5, stat = \"unique\", parse = TRUE, vjust = -.8) +\n",
    "    theme_bw() +\n",
    "    theme(legend.position =\"none\",\n",
    "          plot.title = element_text(size = 25, hjust = 0.5, face = \"bold\", \n",
    "                                    color = \"black\"),\n",
    "          plot.subtitle = element_text(size = 18, hjust = 0.5, face = \"italic\", \n",
    "                                       color = \"gray15\"),\n",
    "          text = element_text(size = 18),\n",
    "          axis.text = element_text(size = 12, face =\"bold\"),\n",
    "          axis.title = element_text(size = 18, face = \"bold\"),\n",
    "          plot.caption = element_text(size = 10))+\n",
    "    xlab(paste(\"\\n\",group)) +\n",
    "    ylab(\"Rank\\n\") +\n",
    "    labs(title = \"Pairwise ANOSIM\",\n",
    "         subtitle = \"Holm-Bonferroni correction of P-values for \\\n",
    "         compact letter display\\n\",\n",
    "         caption = paste(\"Beta Dispersion\", .beta, \n",
    "                         \"\\tPermutations:\", permutations)) +\n",
    "    expand_limits(y = c(0, max(.plt$X2)*1.1))\n",
    "  return(.p)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "462fd975",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_specie <- pairwise_anosim_plot(data_philr, \n",
    "                                   bac_metadata, \n",
    "                                   \"morpho_specie\",\n",
    "                                   distance = \"euclidian\",\n",
    "                                   permutations = 999) + \n",
    "  scale_shape_manual(values=rep(c(21,22,23,24,25), 5)) \n",
    "#+\n",
    "  #scale_colour_manual(values = \n",
    "   #                     c(\"dodgerblue2\", \"#E31A1C\", \"green4\", \"#6A3D9A\", \n",
    "    #                      \"#FF7F00\", \"gold1\", \"maroon\", \"orchid1\", \"deeppink1\", \n",
    "     #                     \"blue1\", \"steelblue4\", \"darkturquoise\", \"green1\", \n",
    "      #                    \"yellow4\",\"green\",\"blue\")) +\n",
    "  #scale_fill_manual(values = \n",
    "   #                   c(\"dodgerblue2\", \"#E31A1C\", \"green4\", \"#6A3D9A\", \n",
    "    #                    \"#FF7F00\", \"gold1\", \"maroon\", \"orchid1\", \"deeppink1\", \n",
    "     #                   \"blue1\", \"steelblue4\", \"darkturquoise\", \"green1\", \n",
    "      #                  \"yellow4\",\"green\",\"blue\"))\n",
    "    \n",
    "plt_specie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e3ba2f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_department <- pairwise_anosim_plot(data_philr, \n",
    "                                       bac_metadata, \n",
    "                                       \"department\",\n",
    "                                       distance = \"euclidian\",\n",
    "                                       permutations = 999) + \n",
    "  scale_fill_lancet() + \n",
    "  scale_color_lancet()\n",
    "\n",
    "plt_department"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5877a145",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_clade <- pairwise_anosim_plot(data_philr, \n",
    "                                  bac_metadata, \n",
    "                                  \"Clade\",\n",
    "                                  distance = \"euclidian\",\n",
    "                                  permutations = 999) +\n",
    "  scale_shape_manual(values=rep(c(21,22,23,24,25), 5)) + \n",
    "  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))\n",
    "\n",
    "plt_clade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f81be910",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep <- pairwise.adonis(bac_feature_table_css, \n",
    "                           bac_metadata_css$department, \n",
    "                           perm = 999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "126a6124",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "306c451c",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep <- pairwise.adonis(bac_feature_table_css, \n",
    "                           bac_metadata_css$department, \n",
    "                           perm = 999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfda6d3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02071ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep <- pairwise.adonis(bac_feature_table_css, \n",
    "                           bac_metadata_css$department, \n",
    "                           perm = 999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb3bba2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.dep %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:100%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = F, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a390dc65",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "\n",
    "# t-SNE\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ebc9cf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 6, repr.plot.height = 4.5, repr.plot.res = 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "065e44d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsne_indata <- merge(bac_feature_table, bac_metadata, by = 'row.names')\n",
    "rownames(tsne_indata) <- tsne_indata$Row.names\n",
    "tsne_indata <- tsne_indata[,-1]\n",
    "\n",
    "bac_tsne_dat <- Rtsne::Rtsne(tsne_indata[, 1:(ncol(tsne_indata)-12)])\n",
    "\n",
    "# extract data\n",
    "tsne_plot <- data.frame(x = bac_tsne_dat$Y[, 1],\n",
    "                        y = bac_tsne_dat$Y[, 2],\n",
    "                        tsne_indata[, (ncol(tsne_indata)-11):ncol(tsne_indata)])\n",
    "\n",
    "  ggplot(tsne_plot) + \n",
    "  geom_point(aes(x = x, y = y, color = department, size = depth)) +\n",
    "  scale_color_lancet() +\n",
    "  theme_bw() +\n",
    "  theme(\n",
    "    # text = element_text(size = 18),\n",
    "    # axis.text = element_text(size = 12),\n",
    "    # axis.title = element_text(size = 14, face = \"bold\"),\n",
    "    plot.title = element_text(hjust = 0.5)\n",
    "  ) +\n",
    "  ggtitle(\"T-Distributed Stochastic Neighbor Embedding \\n\") +\n",
    "  xlab(\"\") +\n",
    "  ylab(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bbc99b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsne_indata <- merge(bac_feature_table_css, bac_metadata_css, by = 'row.names')\n",
    "rownames(tsne_indata) <- tsne_indata$Row.names\n",
    "tsne_indata <- tsne_indata[,-1]\n",
    "\n",
    "bac_tsne_dat <- Rtsne::Rtsne(tsne_indata[, 1:(ncol(tsne_indata)-12)])\n",
    "\n",
    "# extract data\n",
    "tsne_plot <- data.frame(x = bac_tsne_dat$Y[, 1],\n",
    "                        y = bac_tsne_dat$Y[, 2],\n",
    "                        tsne_indata[, (ncol(tsne_indata)-11):ncol(tsne_indata)])\n",
    "\n",
    "  ggplot(tsne_plot) + \n",
    "  geom_point(aes(x = x, y = y, color = department, size = depth)) +\n",
    "  scale_color_lancet() +\n",
    "  theme_bw() +\n",
    "  theme(\n",
    "    # text = element_text(size = 18),\n",
    "    # axis.text = element_text(size = 12),\n",
    "    # axis.title = element_text(size = 14, face = \"bold\"),\n",
    "    plot.title = element_text(hjust = 0.5)\n",
    "  ) +\n",
    "  ggtitle(\"T-Distributed Stochastic Neighbor Embedding (CSS)\\n\") +\n",
    "  xlab(\"\") +\n",
    "  ylab(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "870d2ef3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare data\n",
    "A.tsne_indata <- as.matrix(tsne_indata[, 1:(ncol(tsne_indata)-12)])\n",
    "# add pseudo count \n",
    "A.tsne_indata[A.tsne_indata == 0] <- 1e-12\n",
    "# transform to Aitchison distance\n",
    "A.vars <- t(apply(A.tsne_indata, 1, function(x)log(x/exp(mean(log(x))))))\n",
    "A.dist <- as.matrix(dist(A.vars))\n",
    "\n",
    "# compute t-SNE with arbitrary parameters\n",
    "bac_tsne_dat <- Rtsne(A.dist, \n",
    "                      initial_dims = 5, # manipulate\n",
    "                      perplexity = 25,   # \"\n",
    "                      max_iter = 50000,   # \"\n",
    "                      is_distance = T)\n",
    "\n",
    "# extract data for plotting\n",
    "tsne_plot <- data.frame(x = bac_tsne_dat$Y[, 1],\n",
    "                        y = bac_tsne_dat$Y[, 2],\n",
    "                        tsne_indata[, (ncol(tsne_indata)-11):ncol(tsne_indata)])\n",
    "\n",
    "  ggplot(tsne_plot) + \n",
    "  geom_point(aes(x = x, y = y, color = department, size = depth)) +\n",
    "  scale_color_lancet() +\n",
    "  theme_bw() +\n",
    "  theme(\n",
    "    #text = element_text(size = 18),\n",
    "    #axis.text = element_text(size = 12),\n",
    "    #axis.title = element_text(size = 14, face = \"bold\"),\n",
    "    plot.title = element_text(hjust = 0.5, size = 16)\n",
    "  ) +\n",
    "  ggtitle(\"T-Distributed Stochastic Neighbor Embedding (Aitchison Distance)\\n\") +\n",
    "  xlab(\"\") +\n",
    "  ylab(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043ea289",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(tsne_plot) + \n",
    "  geom_point(aes(x = x, y = y, color = Clade, \n",
    "                 fill = Clade, shape = Clade), size = 3) +\n",
    "  scale_shape_manual(values=rep(c(21,22,23,24,25), 5)) +\n",
    "  scale_colour_manual(values = \n",
    "                        c(\"dodgerblue2\", \"#E31A1C\", \"green4\", \"#6A3D9A\", \n",
    "                          \"#FF7F00\", \"gold1\", \"maroon\", \"orchid1\", \"deeppink1\", \n",
    "                          \"blue1\", \"steelblue4\", \"darkturquoise\", \"green1\", \n",
    "                          \"yellow4\",\"darkgreen\", \"black\", \"thistle4\")) +\n",
    "  scale_fill_manual(values = \n",
    "                      c(\"dodgerblue2\", \"#E31A1C\", \"green4\", \"#6A3D9A\", \n",
    "                        \"#FF7F00\", \"gold1\", \"maroon\", \"orchid1\", \"deeppink1\", \n",
    "                        \"blue1\", \"steelblue4\", \"darkturquoise\", \"green1\", \n",
    "                        \"yellow4\",\"darkgreen\",\"black\", \"thistle4\")) +\n",
    "  theme_bw() +\n",
    "  theme(\n",
    "    #text = element_text(size = 18),\n",
    "    #axis.text = element_text(size = 12),\n",
    "    #axis.title = element_text(size = 14, face = \"bold\"),\n",
    "    plot.title = element_text(hjust = 0.5, size = 18)\n",
    "  ) +\n",
    "  ggtitle(\"T-Distributed Stochastic Neighbor Embedding (Aitchison Distance)\\n\") +\n",
    "  xlab(\"\") +\n",
    "  ylab(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39f1c614",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "\n",
    "# Composition\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0886b036",
   "metadata": {},
   "outputs": [],
   "source": [
    "to_incidence_freq <- function(feature_table=bac_feature_table, \n",
    "                              meta_data=bac_metadata, var_group){\n",
    "  \n",
    "  group <- deparse(substitute(var_group))\n",
    "  a <- as.data.frame(ifelse(feature_table > 0, 1, 0))\n",
    "  a$sampleid <- rownames(a)\n",
    "  b <- dplyr::left_join(a, meta_data[,c(\"sampleid\", group)], by = \"sampleid\")\n",
    "  b <- b[,c(1:(ncol(b)-2),ncol(b))]\n",
    "  c <- b %>% dplyr::group_by({{ var_group }}) %>% dplyr::summarise(dplyr::across(dplyr::everything(), sum))\n",
    "  class(c) <- \"data.frame\"\n",
    "  d <- as.list(as.data.frame(t(c)))\n",
    "  names(d) <- base::unique(meta_data[,group])\n",
    "  d <- lapply(d, function(x) \n",
    "    as.numeric(c(length(b[,group][b[,group] == x[1]]), x[2:length(x)])))\n",
    "  d <- lapply(d, function(x) {x[x!=0]})\n",
    "  return(d)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163de0c0",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "\n",
    "<h2>Coverage - Clade</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ff7fa35",
   "metadata": {},
   "outputs": [],
   "source": [
    "next_data_Cla <- to_incidence_freq(var_group=Clade)\n",
    "\n",
    "next_res_Cla <- iNEXT(next_data_Cla, q=c(0,1,2), datatype=\"incidence_freq\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d72c97b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 10, repr.plot.height = 4, repr.plot.res = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9a40807",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_next_cov_Cla <- ggiNEXT(next_res_Cla, type=2, se=TRUE, facet.var=\"site\", \n",
    "                         color.var=\"site\", grey=TRUE) +\n",
    "  theme_bw() +\n",
    "  theme(legend.position = \"none\")\n",
    "  \n",
    "plot_next_cov_Cla"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd85e231",
   "metadata": {},
   "source": [
    "Sample coverage is in most cases high enough for reliable analysis, i.e. coverage of >= 50% (Chao and Jost 2012).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6373d733",
   "metadata": {},
   "source": [
    "\n",
    "<hr><br>\n",
    "\n",
    "<h2>Richness - Clade</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57e0ea56",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_next_rich_Cla <- ggiNEXT(next_res_Cla, type=1, se=TRUE, facet.var=\"site\", \n",
    "                           color.var=\"site\") +\n",
    "  scale_color_lancet() +\n",
    "  scale_fill_lancet() +\n",
    "  theme_bw() +\n",
    "  theme(legend.position = \"bottom\")\n",
    "\n",
    "plot_next_rich_Cla[[\"guides\"]][[\"colour\"]][[\"title\"]] <- \"Hill Number\"\n",
    "plot_next_rich_Cla[[\"guides\"]][[\"fill\"]][[\"title\"]] <- \"Hill Number\"\n",
    "plot_next_rich_Cla[[\"guides\"]][[\"shape\"]][[\"title\"]] <- \"Hill Number\"\n",
    "\n",
    "plot_next_rich_Cla"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfb5f1f4",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "\n",
    "<h2>Coverage - Department</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc6e4d04",
   "metadata": {},
   "outputs": [],
   "source": [
    "next_data_dep <- to_incidence_freq(var_group=department)\n",
    "\n",
    "next_res_dep <- iNEXT(next_data_dep, q=c(0,1,2), datatype=\"incidence_freq\")\n",
    "\n",
    "\n",
    "plot_next_cov_dep <- ggiNEXT(next_res_dep, type=2, se=TRUE, facet.var=\"none\", \n",
    "                             color.var=\"site\") +\n",
    "  scale_color_lancet() +\n",
    "  scale_fill_lancet()\n",
    "\n",
    "data1 <- rbind(t(plot_next_cov_dep[[\"layers\"]][[1]][[\"data\"]][[\"x\"]]),\n",
    "               t(plot_next_cov_dep[[\"layers\"]][[1]][[\"data\"]][[\"y\"]])) %>% \n",
    "  as.matrix() %>% \n",
    "  t() %>% \n",
    "  as.data.frame() %>% \n",
    "  \"colnames<-\"(c(\"x\",\"y\")) %>% \n",
    "  mutate(department = next_res_dep[[\"DataInfo\"]][[\"site\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a1b959",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_next_cov_dep <- \n",
    "  plot_next_cov_dep + \n",
    "  geom_point(data = data1, aes(color = department, fill = department,\n",
    "                              shape = department), size =6) +\n",
    "  scale_shape_manual(values=c(21:25)) +\n",
    "  theme_bw()\n",
    "\n",
    "plot_next_cov_dep"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3be717",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "\n",
    "<h2>Richness - Department</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c29e889b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_next_rich_dep <- ggiNEXT(next_res_dep, type=1, se=TRUE, facet.var=\"order\", \n",
    "                          color.var=\"site\") +\n",
    "  scale_color_lancet() +\n",
    "  scale_fill_lancet()\n",
    "\n",
    "\n",
    "data2 <- rbind(t(plot_next_rich_dep[[\"layers\"]][[1]][[\"data\"]][[\"x\"]]),\n",
    "               t(plot_next_rich_dep[[\"layers\"]][[1]][[\"data\"]][[\"y\"]])) %>% \n",
    "  as.matrix() %>% \n",
    "  t() %>% \n",
    "  as.data.frame() %>% \n",
    "  \"colnames<-\"(c(\"x\",\"y\")) %>% \n",
    "  mutate(department = rep(next_res_dep[[\"DataInfo\"]][[\"site\"]],each = 3)) %>% \n",
    "  mutate(order = rep(0:2,5))\n",
    "\n",
    "plot_next_rich_dep <- \n",
    "  plot_next_rich_dep + \n",
    "  geom_point(data = data2, aes(color = department, fill = department,\n",
    "                              shape = department), size = 6) +\n",
    "  scale_shape_manual(values=c(21:25)) +\n",
    "  theme_bw()\n",
    "\n",
    "plot_next_rich_dep "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98b679cb",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "\n",
    "<h2>Evenness</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7c26b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "dat_even <- as.data.frame(evenness(bac_ps))\n",
    "\n",
    "dat_even <- merge(dat_even, bac_metadata_css, by = \"row.names\")\n",
    "\n",
    "plot_evenness <-\n",
    "  ggplot(dat_even, aes(x=department, y=pielou, color=department, \n",
    "                       fill=department)) + \n",
    "  geom_violin(trim = FALSE, alpha = 0.6) +\n",
    "  geom_boxplot(width = 0.1, color = \"black\", outlier.color = \"black\",\n",
    "               outlier.alpha = 0.5, outlier.size = 2.5) +\n",
    "  theme_bw() +\n",
    "  scale_fill_lancet() +\n",
    "  scale_color_lancet() + \n",
    "  xlab(\"\\ndepartment\") +\n",
    "  ylab(\"Pielou's Evenness Index\\n\") +\n",
    "  theme(\n",
    "    text = element_text(size = 14),\n",
    "    axis.text = element_text(size = 10),\n",
    "    axis.text.x = element_text(angle = -20),\n",
    "    axis.title = element_text(size = 12, face = \"bold\"),\n",
    "    plot.subtitle = element_text(vjust = 0.5, hjust = 0.5, size = 10),\n",
    "    plot.title = element_text(hjust = 0.5, size = 18),\n",
    "    legend.position = \"none\")\n",
    "\n",
    "\n",
    "# bf.test(pielou ~ department, data = dat_even)\n",
    "\n",
    "# hist(dat_even$pielou)\n",
    "\n",
    "# shapiro.test(dat_even$pielou)\n",
    "\n",
    "# oneway.test(pielou ~ department, data = dat_even)\n",
    "\n",
    "# kruskal.test(pielou ~ department, data = dat_even)\n",
    "\n",
    "even_aov <- aov(pielou ~ department, data = dat_even)\n",
    "\n",
    "# plot(even_aov)\n",
    "\n",
    "aov_test <- summary(even_aov)\n",
    "\n",
    "even_aov_post <- TukeyHSD(even_aov)\n",
    "\n",
    "plot_evenness <- \n",
    "  plot_evenness + labs(subtitle = \n",
    "                         paste(\"AOV test:\",\n",
    "                               \"  F value = \", \n",
    "                               signif(aov_test[[1]]$`F value`[1], 4), \n",
    "                               \",  Df = \",\n",
    "                               aov_test[[1]]$Df[1], \n",
    "                               \",  Pr(>F) = \", \n",
    "                               signif(aov_test[[1]]$`Pr(>F)`[1], 2), \n",
    "                               \"\\n\", sep = \"\", collapse = \"\"),\n",
    "                       title = \"Evenness\\n\") +\n",
    "  geom_segment(aes(x=3,y=1.25,xend=5,yend=1.25), color=\"black\") +\n",
    "  annotate(\"text\", x = 4, y = 1.3, label = \"**\") +\n",
    "  ggtitle(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bf35cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 10, repr.plot.height = 6, repr.plot.res = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17ac7a41",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_evenness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a7fc5bf",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "\n",
    "<h2>Clade X Department</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d36bdc21",
   "metadata": {},
   "outputs": [],
   "source": [
    "x <- as.data.frame(cbind(bac_metadata$Clade, bac_metadata$department))\n",
    "x$V3 = 1\n",
    "y <- cast(x, V1~V2, value=\"V3\")\n",
    "y %>% \n",
    "  mutate(`Bocas del Toro` = \n",
    "           cell_spec(`Bocas del Toro`, background = \n",
    "                       ifelse(`Bocas del Toro` != 0, \"white\",\"#ff8178\"))) %>% \n",
    "  mutate(`Cartagena` = \n",
    "           cell_spec(`Cartagena`, background = \n",
    "                       ifelse(`Cartagena` != 0,\"white\",\"#ff8178\"))) %>% \n",
    "  mutate(`Curazao` = \n",
    "           cell_spec(`Curazao`, background = \n",
    "                       ifelse(`Curazao` != 0, \"white\",\"#ff8178\"))) %>% \n",
    "  mutate(`Long Keys Florida` = \n",
    "           cell_spec(`Long Keys Florida`, background = \n",
    "                       ifelse(`Long Keys Florida` != 0, \"white\",\"#ff8178\"))) %>% \n",
    "  mutate(`San Andres Isla` = \n",
    "           cell_spec(`San Andres Isla`, background = \n",
    "                       ifelse(`San Andres Isla` != 0, \"white\",\"#ff8178\"))) %>% \n",
    "  kable(align = \"c\", format = \"html\", \n",
    "        table.attr = \"style='width:60%;'\", escape = FALSE) %>%\n",
    "  column_spec(1, bold = T) %>% \n",
    "  column_spec(c(1:ncol(y)), width = \"18em\") %>% \n",
    "  kable_styling(c(\"striped\", \"bordered\"), full_width = T, font_size = 14, \n",
    "                position = \"center\") %>%\n",
    "as.character() %>%\n",
    "display_html()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "471b644d",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "\n",
    "# iNEXT\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7f3205f",
   "metadata": {},
   "source": [
    "\n",
    "<h2>Point Estimate (50% Coverage)</h2>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68e739ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "next_50_cov_Cla <- estimateD(next_data_Cla, datatype=\"incidence_freq\",\n",
    "          base=\"coverage\", level=0.5, conf=0.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d1b6a7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(next_50_cov_Cla, aes(x=as.factor(order), y=qD, fill=site, ymin=qD.LCL, \n",
    "                            ymax=qD.UCL)) +\n",
    "  geom_col(stat=\"identity\", position =\"dodge2\") +\n",
    "  facet_grid(~as.factor(order), scale=\"free_x\") +\n",
    "  geom_errorbar(stat=\"identity\", position =\"dodge2\") +\n",
    "  geom_text(aes(label=site, y=15), \n",
    "            position = position_dodge(width = 0.9), angle = 90, hjust=0.05) +\n",
    "  theme_bw() +\n",
    "  theme(legend.position = \"none\") +\n",
    "  scale_fill_grey() + xlab(\"\") + ylab(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f80093d9",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "\n",
    "<h2>Point Estimate (90% Coverage)</h2>\n",
    "<br>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26f2e223",
   "metadata": {},
   "outputs": [],
   "source": [
    "next_90_cov_dep <- estimateD(next_data_dep, datatype=\"incidence_freq\",\n",
    "                         base=\"coverage\", level=0.9, conf=0.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b8f8e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(next_90_cov_dep, aes(x=as.factor(order), y=qD, fill=site, ymin=qD.LCL, \n",
    "                            ymax=qD.UCL)) +\n",
    "  geom_col(stat=\"identity\", position =\"dodge2\") +\n",
    "  facet_wrap(~as.factor(order), scale=\"free_x\") +\n",
    "  geom_errorbar(stat=\"identity\", position =\"dodge2\") +\n",
    "  theme_bw() +\n",
    "  scale_fill_grey() + xlab(\"\") + ylab(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5db87a4d",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='CoreMicrobiome'></a>\n",
    "\n",
    "# Core Microbiome\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef46845",
   "metadata": {},
   "outputs": [],
   "source": [
    "taxa <- read.table(\"/root/biommar/data/taxa.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49cbaa6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make data relative\n",
    "bac_feature_table_relative <- as.data.frame(apply(\n",
    "  bac_feature_table, 1, function(x) { x / sum(x)}\n",
    "))\n",
    "\n",
    "# add metadata\n",
    "bac_feature_table_relative_taxa <- merge(bac_feature_table_relative, \n",
    "                                         taxa, by = \"row.names\") \n",
    "\n",
    "rownames(bac_feature_table_relative_taxa) <- \n",
    "  bac_feature_table_relative_taxa$Row.names\n",
    "\n",
    "bac_feature_table_relative_taxa  <- bac_feature_table_relative_taxa[,-1]\n",
    "\n",
    "\n",
    "# filter Mycoplasmatales and Oceanospirillales\n",
    "bac_sub_tax <-\n",
    "  bac_feature_table_relative_taxa[bac_feature_table_relative_taxa$order %in%\n",
    "                                    c(\"Mycoplasmatales\",\"Oceanospirillales\"), ]\n",
    "\n",
    "# merge data based on order (exclude last 8 columns which are taxa info)\n",
    "abund <- \n",
    "  cbind(\n",
    "  colSums(bac_sub_tax[bac_sub_tax$order %in% \n",
    "                        c(\"Mycoplasmatales\"), c(1:(ncol(bac_sub_tax)-8))]),\n",
    "  colSums(bac_sub_tax[bac_sub_tax$order %in% \n",
    "                        c(\"Oceanospirillales\"), c(1:(ncol(bac_sub_tax)-8))])\n",
    "  )\n",
    "\n",
    "colnames(abund) <- c(\"Mycoplasmatales\",\"Oceanospirillales\")\n",
    "\n",
    "abund <- as.data.frame(abund)\n",
    "\n",
    "\n",
    "ggplot(abund) +\n",
    "  geom_density(aes(x = Mycoplasmatales, fill = \"Mycoplasmatales\"),\n",
    "               color = \"gray10\", alpha = 0.5) +\n",
    "  geom_density(aes(x = Oceanospirillales, fill = \"Oceanospirillales\"),\n",
    "               color = \"gray10\", alpha = 0.5) +\n",
    "  labs(x = \"Relative abundance\", title = \"\") +\n",
    "  theme_classic() + # Changes the background\n",
    "  scale_x_continuous(label = scales::percent) +\n",
    "  scale_fill_lancet() +\n",
    "  guides(fill=guide_legend(title=\"Order\")) +\n",
    "  theme(\n",
    "    text = element_text(size = 18),\n",
    "    axis.text = element_text(size = 12),\n",
    "    axis.text.x = element_text(angle = -20),\n",
    "    axis.title = element_text(size = 14, face = \"bold\"),\n",
    "    plot.subtitle = element_text(vjust = 1.5, hjust = 0.85, size = 12),\n",
    "    plot.title = element_text(hjust = 0.5, size = 24))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2eb0d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove zero counts\n",
    "bac_ps <- prune_taxa(taxa_sums(bac_ps) > 0, bac_ps)\n",
    "\n",
    "# In transformation typ, the 'compositional' abundances are returned as\n",
    "# relative abundances in [0, 1]\n",
    "bac_ps_comp <- microbiome::transform(bac_ps, \"compositional\")\n",
    "\n",
    "\n",
    "plot_heatmap(bac_ps_comp, method = \"NMDS\", distance = \"jaccard\") +\n",
    "  geom_raster(aes(fill = log10(Abundance))) + \n",
    "  scale_fill_viridis_c(option=\"magma\", \n",
    "                       na.value = \"gray90\", \n",
    "                       name = \"log10(Abundance)\",\n",
    "                       direction = -1) +\n",
    "  theme(axis.text.y = element_text(size=5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a26c5206",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define core microbiome based on detection and prevalence\n",
    "bac_core_taxa <- core_members(bac_ps_comp,\n",
    "                              detection = 0.001,\n",
    "                              prevalence = 0.1)\n",
    "\n",
    "bac_core_taxa <- as.data.frame(bac_core_taxa)\n",
    "colnames(bac_core_taxa) <- c(\"sampleid\")\n",
    "\n",
    "write.table(bac_core_taxa,\n",
    "            file = \"/root/biommar/data/bac_core_taxa.tsv\",\n",
    "            sep = \"\\t\",\n",
    "            quote = F,\n",
    "            row.names = F)\n",
    "\n",
    "# add metadata\n",
    "colnames(bac_core_taxa) <- c(\"ID\")\n",
    "rownames(bac_core_taxa) <- bac_core_taxa$ID\n",
    "bac_core_taxa <- merge(bac_core_taxa, taxa, by = \"row.names\")\n",
    "\n",
    "prevalences <- seq(.1, .5, .05)\n",
    "detections <- round(10 ^ seq(log10(1e-3), log10(.2), length = 10), 5)\n",
    "\n",
    "plot_bac_ps_core <- \n",
    "  plot_core(\n",
    "    bac_ps_comp,\n",
    "    plot.type = \"heatmap\",\n",
    "    colours = gray,\n",
    "    prevalences = prevalences,\n",
    "    detections = detections,\n",
    "    min.prevalence = .1) +\n",
    "  xlab(\"\\nDetection Threshold (Relative Abundance (%))\") +\n",
    "  theme_bw() +\n",
    "  #ggtitle(paste(\"Identification Level :\", filter_tax)) +\n",
    "  ylab(\"ASVs\") +\n",
    "  scale_fill_viridis_c()\n",
    "\n",
    "# replace ID by taxonomic name\n",
    "names <- \n",
    "  plyr:::mapvalues(plot_bac_ps_core[[\"data\"]][[\"Taxa\"]], \n",
    "                   from = as.character(bac_core_taxa$ID), \n",
    "                   to = as.character(bac_core_taxa$genus))\n",
    "\n",
    "plot_bac_ps_core + scale_y_discrete(labels = names) +\n",
    "  theme(axis.text.y = element_text(size = 6),\n",
    "        panel.border = element_blank(),\n",
    "        panel.grid = element_blank())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e38fb18e",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "### Save R Environment\n",
    "\n",
    "Saving the objects in the R environment to a file lets you access them faster after restarting the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40132f02",
   "metadata": {},
   "outputs": [],
   "source": [
    "save.image(\"/root/biommar/data/environment.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "493a393b",
   "metadata": {},
   "source": [
    "<br><hr><br><br><hr><br>\n",
    "\n",
    "# Interlude in bash\n",
    "\n",
    "<span style=\"color:red\">**PLEASE CHANGE KERNEL ACCORDINGLY!**</span>\n",
    "\n",
    "Select ```Kernel``` in the menu bar at the top, click on ```Change kernel```, and select ```Python 3 (ipykernel)```.\n",
    "\n",
    "<br><hr><br><br><hr><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "817c6c38",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime feature-table filter-seqs \\\n",
    "    --i-data /root/biommar/data/dada_filtered_rep_seqs.qza \\\n",
    "    --m-metadata-file /root/biommar/data/bac_core_taxa.tsv \\\n",
    "    --o-filtered-data /root/biommar/data/core_sequences.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "085a6ad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "source activate qiime2-2022.2\n",
    "\n",
    "qiime phylogeny align-to-tree-mafft-fasttree \\\n",
    "    --i-sequences /root/biommar/data/core_sequences.qza \\\n",
    "    --p-n-threads 4 \\\n",
    "    --o-alignment /root/biommar/data/core_alignemnt.qza \\\n",
    "    --o-masked-alignment /root/biommar/data/core_masked_alignment.qza \\\n",
    "    --o-tree /root/biommar/data/core_unrooted_tree.qza \\\n",
    "    --o-rooted-tree /root/biommar/data/core_rooted_tree.qza\n",
    "\n",
    "cd  /root/biommar/data/\n",
    "unzip -qq core_rooted_tree.qza\n",
    "cp $(echo $(ls -td -- /root/biommar/data/*/ | head -n1) $(echo \"data/tree.nwk\") | sed 's/\\ //g') \\\n",
    "   /root/biommar/data/core_tree.nwk\n",
    "rm -r $(ls -td -- /root/biommar/data/* | head -n 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d95fa532",
   "metadata": {},
   "source": [
    "<br><hr><br><br><hr><br>\n",
    "\n",
    "# Continue in R\n",
    "\n",
    "<span style=\"color:red\">**PLEASE CHANGE KERNEL ACCORDINGLY!**</span>\n",
    "\n",
    "Select ```Kernel``` in the menu bar at the top, click on ```Change kernel```, and select ```R```.\n",
    "\n",
    "<br><hr><br><br><hr><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "062fbd17",
   "metadata": {},
   "outputs": [],
   "source": [
    "library( ALDEx2,        quietly = TRUE)\n",
    "library( ape,           quietly = TRUE)\n",
    "library( aplot,         quietly = TRUE)\n",
    "library( cluster,       quietly = TRUE)\n",
    "library( ggdendro,      quietly = TRUE)\n",
    "library( grid,          quietly = TRUE)\n",
    "library( glmnet,        quietly = TRUE)\n",
    "library( ggplot2,       quietly = TRUE)\n",
    "library( ggfortify,     quietly = TRUE)\n",
    "library( ggthemes,      quietly = TRUE)\n",
    "library( gridExtra,     quietly = TRUE)\n",
    "library( ggpmisc,       quietly = TRUE)\n",
    "library( ggsci,         quietly = TRUE)\n",
    "library( ggtree,        quietly = TRUE)\n",
    "library( tidyr,         quietly = TRUE)\n",
    "library( metagenomeSeq, quietly = TRUE)\n",
    "library( microbiome,    quietly = TRUE)\n",
    "library( phyloseq,      quietly = TRUE)\n",
    "library( philr,         quietly = TRUE)\n",
    "library( mia,           quietly = TRUE)\n",
    "library( vegan,         quietly = TRUE)\n",
    "library( multcompView,  quietly = TRUE)\n",
    "library( plyr,          quietly = TRUE)\n",
    "library( dplyr,         quietly = TRUE)\n",
    "library( Rtsne,         quietly = TRUE)\n",
    "library( plotly,        quietly = TRUE)\n",
    "library( zgtools,       quietly = TRUE)\n",
    "library( reshape,       quietly = TRUE)\n",
    "library( onewaytests,   quietly = TRUE)\n",
    "library( iNEXT,         quietly = TRUE)\n",
    "library( phytools,      quietly = TRUE)\n",
    "library( reshape2,      quietly = TRUE)\n",
    "library( repr,          quietly = TRUE)\n",
    "library( kableExtra,    quietly = TRUE)\n",
    "library( IRdisplay,     quietly = TRUE)\n",
    "library( jsonlite,     quietly = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3fdb7d4",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "### Load R Environment\n",
    "\n",
    "Load the environment to restore objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cc9eaa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "load(\"/root/biommar/data/environment.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f20fb08",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='HostPhylogeny'></a>\n",
    "\n",
    "# Host Phylogeny\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8d9b3ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 10, repr.plot.height = 12, repr.plot.res = 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f53851d",
   "metadata": {},
   "source": [
    "### Eunicea phylogeny plot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6cfbd33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the tree produced in Qiime2\n",
    "host_tree <- read_tree(\"/root/biommar/data/subd_params-data4.phy.treefile\")\n",
    "\n",
    "host_tree[[\"tip.label\"]] <- as.character(host_tree[[\"tip.label\"]])\n",
    "\n",
    "# get an overview of the phylogeny\n",
    "\n",
    "ggtree(host_tree) + \n",
    "  geom_tiplab(size=2, align=TRUE, linesize=.5) +\n",
    "  ggplot2::xlim(0, 0.05)\n",
    "\n",
    "# crop off prefix and suffix added to ID during construction\n",
    "host_tree[[\"tip.label\"]] <- gsub(\"eunicea_\", \"\", host_tree[[\"tip.label\"]]) \n",
    "host_tree[[\"tip.label\"]] <- gsub(\"_ctadpt1\", \"\", host_tree[[\"tip.label\"]])\n",
    "\n",
    "# match metadata\n",
    "host_tree <- keep.tip(host_tree, \n",
    "                      intersect(bac_metadata_css$sampleid, \n",
    "                                host_tree[[\"tip.label\"]]))\n",
    "\n",
    "host_plot <- ggtree(host_tree) + \n",
    "  geom_tiplab(size=2, align=TRUE, linesize=.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5845b7e",
   "metadata": {},
   "source": [
    " combine Eunicea and core microbiome\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0242a16",
   "metadata": {},
   "outputs": [],
   "source": [
    "core_tree <- read.tree(\"/root/biommar/data/core_tree.nwk\")\n",
    "\n",
    "bac_core_taxa <- as.matrix(read.csv(\"/root/biommar/data/bac_core_taxa.tsv\"))\n",
    "\n",
    "# filter data to core microbiome members\n",
    "core_tax <- taxa[rownames(taxa) %in% bac_core_taxa,]\n",
    "core_feature <- bac_feature_table[,colnames(bac_feature_table) %in% \n",
    "                                        bac_core_taxa]\n",
    "core_meta <- bac_metadata[rownames(bac_metadata) %in% \n",
    "                                rownames(core_feature),]\n",
    "\n",
    "# reorder and subsitute tip labels to bear taxonomic labels\n",
    "core_tree[[\"tip.label\"]] <- core_tax[match(core_tree[[\"tip.label\"]], \n",
    "                                           rownames(core_tax)),]$family\n",
    "\n",
    "core_tree[[\"tip.label\"]] <- as.character(core_tree[[\"tip.label\"]])\n",
    "core_feature <- as.data.frame(t(core_feature))\n",
    "\n",
    "\n",
    "#core_feature <- core_feature[,match(get_taxa_name(host_plot), \n",
    "#                                   colnames(core_feature))]\n",
    "\n",
    "core_feature$id <- rownames(core_feature)\n",
    "molten <- melt(core_feature, id.vars=c(\"id\"))\n",
    "molten$sampleid <- molten$variable\n",
    "data_merged <- merge(molten, core_meta[,c(1:12)], by=\"sampleid\")\n",
    "\n",
    "data_merged <- data_merged[data_merged$sampleid %in% host_tree[[\"tip.label\"]],]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83b8ee09",
   "metadata": {},
   "source": [
    "department data preparation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c73d6d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_merged_department <- data_merged[order(data_merged[,\"department\"]),]\n",
    "data_merged_department$sampleid <- \n",
    "  factor(data_merged_department$sampleid, \n",
    "         levels=unique((data_merged_department$sampleid)[order(\n",
    "           data_merged_department[, \"department\"])]))\n",
    "\n",
    "core_meta[, \"department\"] <- as.factor(core_meta[, \"department\"])\n",
    "core_meta_department <- base::transform(core_meta, \n",
    "                                     count = table(department)[department])\n",
    "\n",
    "core_meta_department <- core_meta_department[order(\n",
    "  core_meta_department[, \"department\"]),]\n",
    "core_feature_department <- core_feature[order(\n",
    "  core_meta_department[, \"department\"]),]\n",
    "\n",
    "names <- get_taxa_name(host_plot)\n",
    "\n",
    "tmp_df <- data_merged_department$sampleid[match(rep(names, \n",
    "                                                    each=nrow(bac_core_taxa)), \n",
    "                                    data_merged_department$sampleid)]\n",
    "\n",
    "tmp_df <- as.data.frame(tmp_df)\n",
    "tmp_df$num <- rownames(tmp_df)\n",
    "colnames(tmp_df) <- c(\"sampleid\",\"num\")\n",
    "\n",
    "data_merged_department <- merge(data_merged_department, tmp_df, by=\"sampleid\")\n",
    "data_merged_department$num <- as.numeric(data_merged_department$num)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18168afe",
   "metadata": {},
   "source": [
    "### department plotting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "743ffcba",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph.tree <- ggtree(core_tree) + \n",
    "  geom_tiplab(size=2, align=TRUE, linesize=.5) +\n",
    "  ggplot2::xlim(0, 1.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd0c0186",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph.heat_department <- ggplot(data_merged_department) + \n",
    "  geom_tile(aes(x=id, y=as.character(reorder(sampleid,num)), fill=log10(value+1)), \n",
    "            color = \"#000000\") +\n",
    "  scale_fill_viridis_c() + theme_void() + scale_y_discrete(limit=(names)) +\n",
    "  coord_flip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca694ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph.meta_department <- ggplot(data_merged_department) + \n",
    "  geom_raster(aes(y=1, x=sampleid, fill=department), alpha = 1) +\n",
    "  theme_void() +\n",
    "  scale_fill_manual(values=c(\"#FFC300\",\"#FF5733\",\"#C70039\",\"#900C3F\",\"#581845\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca96f53c",
   "metadata": {},
   "outputs": [],
   "source": [
    "host_plot <- ggtree(host_tree) + \n",
    "ggplot2::xlim(-0.05, 0.05)\n",
    "  \n",
    "host_plot <- host_plot  %<+% core_meta + \n",
    "  layout_dendrogram() +\n",
    "    geom_tippoint(aes(color=Clade)) +\n",
    "    scale_color_manual(values=c(\"#00468BFF\", \"#ED0000FF\", \"#42B540FF\", \n",
    "                                \"#0099B4FF\", \"#925E9FFF\", \"#FDAF91FF\", \n",
    "                                \"#AD002AFF\", \"#ADB6B6FF\", \"#1B1919FF\", \n",
    "                                \"#374E55FF\", \"#DF8F44FF\", \"#900C3F\", \n",
    "                                \"#FFC300\", \"#79AF97FF\", \"green\", \"blue\")) #+\n",
    "  #geom_label(aes(label=node), size = 2) +\n",
    "  #geom_cladelabel(node=194, label=\"mammosa complex\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#AD002AFF\", hjust = 0.2,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  #\n",
    "  #geom_cladelabel(node=215, label=\"pinta\", align=T, \n",
    "  #                angle=90, fontsize=3 , color = \"#1B1919FF\", hjust = 0.2,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=222, label=\"flexuosa\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#0099B4FF\", hjust = 0.5,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  #\n",
    "  #geom_cladelabel(node=174, label=\"asperula - fusca\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#00468BFF\", hjust = 0.7,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=225, label=\"tayrona\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#FFC300\", hjust = 0.6,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=232, label=\"laciniata\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#FDAF91FF\", hjust = 0.2,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=256, label=\"knighti\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=267, label=\"calyclata\", align=T, \n",
    "  #                angle=90, fontsize=3, color = \"#ED0000FF\", hjust = 0.5,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  #\n",
    "  #geom_cladelabel(node=268, label=\"tourneforti\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  #\n",
    "  #geom_cladelabel(node=280, label=\"sp.\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=302, label=\"clavigera\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  #\n",
    "  #geom_cladelabel(node=312, label=\"sp. varias\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5) +\n",
    "  # \n",
    "  #geom_cladelabel(node=315, label=\"sp. curazao\", align=T, \n",
    "  #                angle=90, fontsize=3, color =  \"#925E9FFF\", hjust = 0.3,\n",
    "  #                offset.text=.001, barsize=1.5)\n",
    " \n",
    "#161 - 157\n",
    "#host_plot <- host_plot + geom_tiplab(angle=90, hjust=1, offset=-10, show.legend=FALSE)\n",
    "\n",
    "plot_core_host <-\n",
    "  ph.heat_department %>% insert_top(host_plot, \n",
    "                                   height = 0.5) %>% \n",
    "    insert_left(ph.tree, width=0.5) %>%  insert_bottom(ph.meta_department, \n",
    "                                                      height = 0.1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70924d18",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 10, repr.plot.height = 8, repr.plot.res = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9cc3d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_core_host"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eae5a63d",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "<a id='ALDEx2'></a>\n",
    "# ALDEx2\n",
    "\n",
    "\n",
    "[&#8593; back to content &#8593;](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87082a1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "bac_feature_table <- read.delim(\"/root/biommar/data/B_asv.txt\", check.names = F)\n",
    "bac_metadata <- read.delim(\"/root/biommar/data/B_asv_metadata.txt\")\n",
    "\n",
    "deps <- unique(bac_metadata$department)\n",
    "deps_comb <- combn(deps, 2)\n",
    "\n",
    "for(i in 1:ncol(deps_comb)) {\n",
    "  comb <- deps_comb[,i]\n",
    "  cond <- bac_metadata$department[bac_metadata$department %in% comb]\n",
    "  dat <- t(bac_feature_table[\n",
    "      rownames(bac_feature_table) %in% \n",
    "        rownames(bac_metadata[bac_metadata$department %in% comb,])\n",
    "      ,])\n",
    "  z <- aldex(dat, cond, mc.samples=100, test=\"t\", effect=TRUE,\n",
    "        include.sample.summary=FALSE, denom=\"all\", verbose=FALSE)\n",
    "  assign(paste(\"x\",i, collapse=\"\", sep=\"\"), z) \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "491c4a3b",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "### Function: _aldex.ggplot_\n",
    "\n",
    "ALDEx2 provides a plotting function that uses R base plots. The function below produces similar plots using the ggplot package. Some scaling has to be adjusted manually. In general, ggplot should provide more options for manipulation of plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95d77baa",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot <- \n",
    "  function(x, ..., type = c(\"MW\", \"MA\"), all.col = rgb(0, 0, 0, 0.2), \n",
    "           all.pch = 20,all.cex = 2, called.col = \"red\", called.pch = 20, \n",
    "           called.cex = 2.5, thres.line.col = \"darkgrey\", thres.lwd = 0.5, \n",
    "           test = \"welch\", cutoff.pval = 0.05, cutoff.effect = 1, \n",
    "           rare.col = \"black\", rare = 0, rare.pch = 20, rare.cex = 2){\n",
    "  require(grid)\n",
    "  require(ggplot2)\n",
    "  called <- x$we.eBH <= cutoff.pval\n",
    "  # names must not contain \".\"\n",
    "  cond.1 <- strsplit(colnames(x)[2], \".\", fixed=TRUE)[[1]][3] \n",
    "  cond.2 <- strsplit(colnames(x)[3], \".\", fixed=TRUE)[[1]][3]\n",
    "  \n",
    "  if (type == \"MW\") {\n",
    "  plot <- \n",
    "    ggplot() +\n",
    "    theme_bw() +\n",
    "    geom_abline(aes(intercept=0, slope=1), \n",
    "                color=thres.line.col, linetype=2, size=thres.lwd) +\n",
    "    geom_abline(aes(intercept=0, slope=-1), \n",
    "                color=thres.line.col, linetype=2, size=thres.lwd) +\n",
    "    geom_point(data=x, aes(x=diff.win, y=diff.btw),\n",
    "               color=all.col, shape=all.pch, size=all.cex) +\n",
    "    xlab(expression(\"Median\" ~ ~Log[2] ~ ~\"Dispersion\")) +\n",
    "    ylab(expression(atop(\"Median\" ~ ~Log[2] ~ ~\"Difference\"))) +\n",
    "    scale_color_manual(values=c(rare.col, called.col))\n",
    "  if(nrow(x[x$rab.all < rare,]) > 0) {\n",
    "    plot <- plot + \n",
    "      geom_point(data=x[x$rab.all < rare,], \n",
    "                 aes(x=diff.win, y=diff.btw), color=rare.col, \n",
    "                 shape=rare.pch, size=rare.cex)\n",
    "  }\n",
    "  if(nrow(x[called,]) > 0) {\n",
    "    plot <- plot + \n",
    "      geom_point(data=x[called,], aes(x=diff.win, y=diff.btw), \n",
    "                 color=called.col, shape=called.pch, size=called.cex)\n",
    "  }\n",
    "  plot <- plot + ggtitle(\"\")\n",
    "  print(plot)\n",
    "  grid.text(\n",
    "    cond.1,\n",
    "    rot = 90,\n",
    "    gp=gpar(fontsize=8, col=\"gray\", fontface=\"bold\"),\n",
    "    x=unit(0, \"npc\")+unit(7,'mm'),\n",
    "    y=unit(0, \"npc\")+unit(12,'mm'),\n",
    "    just=c(\"left\", \"top\"))\n",
    "  grid.text(\n",
    "    cond.2,\n",
    "    rot = 90,\n",
    "    gp=gpar(fontsize=8, col=\"gray\", fontface=\"bold\"),\n",
    "    x=unit(0, \"npc\")+unit(7,'mm'),\n",
    "    y=unit(1, \"npc\")-unit(10,'mm'),\n",
    "    just=c(\"right\", \"top\"))\n",
    "  } \n",
    "  if (type == \"MA\") {\n",
    "    plot <- \n",
    "      ggplot() +\n",
    "      theme_bw() +\n",
    "      geom_point(data=x, aes(x=rab.all, y=diff.btw),\n",
    "                 color=all.col, shape=all.pch, size=all.cex) +\n",
    "      xlab(expression(\"Median\" ~ ~Log[2] ~ ~\"Relative Abundance\")) +\n",
    "      ylab(expression(atop(\"Median\" ~ ~Log[2] ~ ~\"Difference\"))) +\n",
    "      scale_color_manual(values=c(rare.col, called.col))\n",
    "    if(nrow(x[x$rab.all < rare,]) > 0) {\n",
    "      plot <- plot + \n",
    "        geom_point(data=x[x$rab.all < rare,], \n",
    "                   aes(x=rab.all, y=diff.btw), color=rare.col, \n",
    "                   shape=rare.pch, size=rare.cex)\n",
    "    }\n",
    "    if(nrow(x[called,]) > 0) {\n",
    "      plot <- plot + \n",
    "        geom_point(data=x[called,], aes(x=rab.all, y=diff.btw), \n",
    "                   color=called.col, shape=called.pch, size=called.cex)\n",
    "    }\n",
    "    plot <- plot + ggtitle(\"\")\n",
    "    print(plot)\n",
    "    grid.text(\n",
    "      cond.1,\n",
    "      rot = 90,\n",
    "      gp=gpar(fontsize=8, col=\"gray\", fontface=\"bold\"),\n",
    "      x=unit(0, \"npc\")+unit(7,'mm'),\n",
    "      y=unit(0, \"npc\")+unit(12,'mm'),\n",
    "      just=c(\"left\", \"top\"))\n",
    "    grid.text(\n",
    "      cond.2,\n",
    "      rot = 90,\n",
    "      gp=gpar(fontsize=8, col=\"gray\", fontface=\"bold\"),\n",
    "      x=unit(0, \"npc\")+unit(7,'mm'),\n",
    "      y=unit(1, \"npc\")-unit(10,'mm'),\n",
    "      just=c(\"right\", \"top\"))\n",
    "  } \n",
    "  }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8cfa89f",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 4, repr.plot.height = 4, repr.plot.res = 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d68dc02b",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Cartagena vs. Curazao</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aac7ce0",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x1, type=\"MW\") \n",
    "aldex.ggplot(x1, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fce3fae1",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Cartagena vs. Long Keys Florida</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfeb9d52",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x2, type=\"MW\")\n",
    "aldex.ggplot(x2, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d7e1232",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Bocas del Toro vs. Cartagena</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "838f458f",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x3, type=\"MW\")\n",
    "aldex.ggplot(x3, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cf18822",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Cartagena vs. San Andres Isla</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78562e02",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x4, type=\"MW\")\n",
    "aldex.ggplot(x4, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b3b97fc",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Curazao vs. Long Keys Florida</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d842ec15",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x5, type=\"MW\")\n",
    "aldex.ggplot(x5, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a6e3db7",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Bocas el Toro vs. Curazao</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f99185e",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x6, type=\"MW\")\n",
    "aldex.ggplot(x6, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d176828a",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Curazao vs. San Andres Isla</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd201811",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x7, type=\"MW\")\n",
    "aldex.ggplot(x7, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "388bd0fb",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Bocas del Toro vs. Long Keys Florida</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7116f57",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x8, type=\"MW\")\n",
    "aldex.ggplot(x8, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "224f9223",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Long Keys Florida vs. San Andres Isla</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fbb174a",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x9, type=\"MW\")\n",
    "aldex.ggplot(x9, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa698d4",
   "metadata": {},
   "source": [
    "<hr><br>\n",
    "<h3>Bocas del Toro vs. San Andres Isla</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af12ee3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "aldex.ggplot(x10, type=\"MW\")\n",
    "aldex.ggplot(x10, type=\"MA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e43de5",
   "metadata": {},
   "source": [
    "<br><hr><br>\n",
    "### Save R Environment\n",
    "\n",
    "Saving the objects in the R environment to a file lets you access them faster after restarting the notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "684c3ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "save.image(\"/root/biommar/data/environment.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09516ca2",
   "metadata": {},
   "source": [
    "When closing the notebook and shutting down the kernel, objects in the environment are removed. To restore objects, run the following cells upon restart. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a66ec356",
   "metadata": {},
   "outputs": [],
   "source": [
    "load(\"/root/biommar/data/environment.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec9b62e3",
   "metadata": {},
   "source": [
    "To see all objects in the current environment, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b43554e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [],
      "text/latex": [],
      "text/markdown": [],
      "text/plain": [
       "character(0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ls()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0466498",
   "metadata": {},
   "source": [
    "<br><hr>\n",
    "[&#8593; back to content &#8593;](#Content)\n",
    "<hr><hr><hr><br>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}