nb-templates/step-diversity-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "15bcfa06",
   "metadata": {},
   "source": [
    " \n",
    "# Alpha and Beta Diversity\n",
    "\n",
    "Different higher-level measures are often used to describe the microbiome in a sample. These do not provide information on changes in the abundance of specific taxa but allow us to access a broader change or difference in the composition of microorganisms. Alpha and beta diversity are examples of such measures.\n",
    "\n",
    "Different measures exist to estimate diversity within a single sample, jointly called alpha diversity. The different measures reflect the richness (number) or distribution (evenness) of a microbial sample or aim to reflect a combination of both properties.\n",
    "\n",
    "Rarefaction curves are often used when calculating alpha diversity indices because increasing numbers of sequenced taxa allow increasingly accurate estimates of total population diversity. Rarefaction curves can therefore be used to estimate the full sample richness, as compared to the observed sample richness.\n",
    "\n",
    "While alpha diversity is a measure of microbiome diversity applicable to a single sample, beta diversity is a measure of the similarity or dissimilarity of two communities. As for alpha diversity, many indices exist, each reflecting different aspects of community heterogeneity. Key differences relate to how the indices value variation in rare species if they consider presence/absence only or incorporate abundance, and how they interpret shared absence. Bray-Curtis dissimilarity is a popular measure that considers both size (overall abundance per sample) and shape (abundance of each taxon) of the communities (Bray, 1957). Beta diversity is an essential measure for many popular statistical methods in ecology, such as ordination-based methods, and is widely used for studying the association between environmental variables and microbial composition.\n",
    "\n",
    "In summary, alpha diversity measures can be seen as a summary statistic of a single population (within-sample diversity), while beta diversity measures are estimates of similarity or dissimilarity between populations (between samples).\n",
    "\n",
    "**Source**: (https://biomcare.com/info/key-terms-in-microbiome-projects/)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e7f2330",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-22T19:36:33.828188Z",
     "start_time": "2022-09-22T19:36:33.817323Z"
    }
   },
   "source": [
    "### STEP : Diversity Analysis\n",
    "\n",
    "Using QIIME2 to create diversity analisys graphs and calculations.\n",
    "\n",
    "- [QIIME2 Workflow Overview](https://docs.qiime2.org/2022.8/tutorials/overview/)\n",
    "\n",
    "\n",
    "#### Methods\n",
    "- [diversity](https://docs.qiime2.org/2022.8/plugins/available/diversity/)\n",
    "- [diversity alpha](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha/)\n",
    "- [diversity alpha_phylogenetic](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-phylogenetic/)\n",
    "- [diversity beta](https://docs.qiime2.org/2022.8/plugins/available/diversity/beta/)\n",
    "- [diversity core_metrics](https://docs.qiime2.org/2022.8/plugins/available/diversity/core-metrics/)\n",
    "- [diversity alpha_group_significance](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-group-significance/)\n",
    "- [diversity beta_group_significance](https://docs.qiime2.org/2022.8/plugins/available/diversity/beta-group-significance/)\n",
    "- [feature_table core_features](https://docs.qiime2.org/2022.8/plugins/available/feature-table/core-features/)\n",
    "- [feature_table summarize](https://docs.qiime2.org/2022.8/plugins/available/feature-table/summarize/)\n",
    "- [taxa filter-table](https://docs.qiime2.org/2022.8/plugins/available/taxa/filter-table/)\n",
    "- [taxa collapse](https://docs.qiime2.org/2022.8/plugins/available/taxa/collapse/)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "843debbf",
   "metadata": {},
   "source": [
    "## Setup and settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "965e28e5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-11-03T20:35:26.198671Z",
     "start_time": "2022-11-03T20:35:06.473173Z"
    }
   },
   "outputs": [],
   "source": [
    "# Importing packages\n",
    "import os\n",
    "import pandas as pd\n",
    "from qiime2 import Artifact\n",
    "from qiime2 import Visualization\n",
    "from qiime2 import Metadata\n",
    "\n",
    "from qiime2.plugins.phylogeny.pipelines import align_to_tree_mafft_fasttree\n",
    "\n",
    "from qiime2.plugins.diversity.pipelines import alpha\n",
    "from qiime2.plugins.diversity.pipelines import beta\n",
    "from qiime2.plugins.diversity.pipelines import core_metrics\n",
    "from qiime2.plugins.diversity.pipelines import alpha_phylogenetic\n",
    "\n",
    "from qiime2.plugins.diversity.visualizers import alpha_group_significance\n",
    "from qiime2.plugins.diversity.visualizers import beta_group_significance\n",
    "from qiime2.plugins.diversity.visualizers import alpha_correlation\n",
    "from qiime2.plugins.diversity.visualizers import beta_rarefaction\n",
    "\n",
    "from qiime2.plugins.taxa.methods import filter_table\n",
    "from qiime2.plugins.taxa.methods import collapse\n",
    "\n",
    "from qiime2.plugins.feature_table.visualizers import tabulate_seqs\n",
    "from qiime2.plugins.feature_table.visualizers import summarize\n",
    "from qiime2.plugins.feature_table.visualizers import core_features\n",
    "from qiime2.plugins.diversity.pipelines import core_metrics_phylogenetic\n",
    "\n",
    "from qiime2.plugins.feature_table.methods import filter_samples\n",
    "from qiime2.plugins.feature_table.methods import filter_seqs\n",
    "\n",
    "from qiime2.plugins.alignment.methods import mafft\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "530ceb60",
   "metadata": {},
   "source": [
    "### Receiving the parameters\n",
    "\n",
    "The following cell can receive parameters using the [papermill](https://papermill.readthedocs.io/en/latest/) tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18d874ea",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:14.793270Z",
     "start_time": "2022-09-29T20:54:14.789367Z"
    },
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "base_dir = os.path.join('/', 'home')\n",
    "metadata_file = os.path.abspath(os.path.join(base_dir, 'data', 'metadatada.tsv'))\n",
    "experiment_name = ''\n",
    "class_col = ''\n",
    "replace_files = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38617d9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:14.804575Z",
     "start_time": "2022-09-29T20:54:14.794927Z"
    }
   },
   "outputs": [],
   "source": [
    "experiment_folder = os.path.abspath(os.path.join(base_dir, 'experiments', experiment_name))\n",
    "img_folder = os.path.abspath(os.path.join(experiment_folder, 'imgs'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea9b90fb",
   "metadata": {},
   "source": [
    "### Defining names, paths and flags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45727eef",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:14.816686Z",
     "start_time": "2022-09-29T20:54:14.806238Z"
    }
   },
   "outputs": [],
   "source": [
    "# QIIME2 Artifacts folder\n",
    "qiime_folder = os.path.join(experiment_folder, 'qiime-artifacts')\n",
    "\n",
    "# Input - DADA2 Artifacts\n",
    "dada2_tabs_path = os.path.join(qiime_folder, 'dada2-tabs.qza')\n",
    "dada2_reps_path = os.path.join(qiime_folder, 'dada2-reps.qza')\n",
    "dada2_stat_path = os.path.join(qiime_folder, 'dada2-stat.qza')\n",
    "\n",
    "# Input - Taxonaomic Artifacts\n",
    "taxonomy_path = os.path.join(qiime_folder, 'metatax.qza')\n",
    "\n",
    "# Create folder to store Alpha files\n",
    "alpha_path = os.path.join(qiime_folder, 'alpha-analysis')\n",
    "if not os.path.exists(alpha_path):\n",
    "    os.makedirs(alpha_path)\n",
    "    print(f'The new directory is created in {alpha_path}')\n",
    "    \n",
    "# Create folder to store Beta files\n",
    "beta_path = os.path.join(qiime_folder, 'beta-analysis')\n",
    "if not os.path.exists(beta_path):\n",
    "    os.makedirs(beta_path)\n",
    "    print(f'The new directory is created in {beta_path}')\n",
    "\n",
    "# Output -Diversity Artifacts\n",
    "alpha_diversity_path = os.path.join(alpha_path, 'alpha-diversity.qza')\n",
    "alpha_diversity_view_path = os.path.join(alpha_path, 'alpha-diversity.qzv')\n",
    "beta_diversity_path = os.path.join(beta_path, 'beta-diversity.qza')\n",
    "beta_diversity_view_path = os.path.join(beta_path, 'beta-diversity.qzv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "601ae263",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:14.829117Z",
     "start_time": "2022-09-29T20:54:14.818272Z"
    }
   },
   "outputs": [],
   "source": [
    "def filter_and_collapse(tab, seqs, tax, meta, lvl, exclude=True, exclude_list='uncultured,unidentified,metagenome'):\n",
    "    from qiime2.plugins.taxa.methods import collapse\n",
    "    from qiime2.plugins.taxa.methods import filter_table\n",
    "    from qiime2.plugins.feature_table.methods import filter_seqs\n",
    "    from qiime2.plugins.feature_table.visualizers import summarize\n",
    "    \n",
    "    to_include = ('d', 'p', 'c', 'o', 'f', 'g', 's')[lvl-1]\n",
    "    to_include += '__'\n",
    "    to_exclude = exclude_list if exclude else None\n",
    "    \n",
    "    filtered_tabs = filter_table(\n",
    "        table=tab, \n",
    "        taxonomy=tax,\n",
    "        include=to_include,\n",
    "        exclude=to_exclude,\n",
    "        mode='contains').filtered_table\n",
    "    \n",
    "    filtered_seqs = filter_seqs(\n",
    "        data = seqs,\n",
    "        table = filtered_tabs,\n",
    "    ).filtered_data\n",
    "    \n",
    "    collapsed_table = collapse(table=filtered_tabs, taxonomy=tax, level=lvl).collapsed_table\n",
    "    collapsed_table_view = summarize(table=collapsed_table, sample_metadata=meta).visualization\n",
    "    \n",
    "    return collapsed_table, collapsed_table_view, filtered_seqs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0999a1b2",
   "metadata": {},
   "source": [
    "## Step execution\n",
    "\n",
    "### Load input files\n",
    "\n",
    "This Step import the QIIME2 `FeatureTable[Frequency]` Artifact and the `Metadata` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f299d1a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:15.559837Z",
     "start_time": "2022-09-29T20:54:14.830405Z"
    }
   },
   "outputs": [],
   "source": [
    "#Load Metadata\n",
    "metadata_qa = Metadata.load(metadata_file)\n",
    "\n",
    "#Load FeatureTable[Frequency]\n",
    "tabs = Artifact.load(dada2_tabs_path)\n",
    "tabs_df = tabs.view(Metadata).to_dataframe().T\n",
    "\n",
    "# FeatureData[Sequence]\n",
    "reps = Artifact.load(dada2_reps_path)\n",
    "\n",
    "# FeatureData[Taxonomy]\n",
    "tax = Artifact.load(taxonomy_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8dfe530",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:15.564968Z",
     "start_time": "2022-09-29T20:54:15.561875Z"
    }
   },
   "outputs": [],
   "source": [
    "# Filter FeatureTable[Frequency | RelativeFrequency | PresenceAbsence | Composition] based on Metadata sample ID values\n",
    "tabs = filter_samples(\n",
    "    table=tabs,\n",
    "    metadata=metadata_qa,\n",
    ").filtered_table\n",
    "# Filter SampleData[SequencesWithQuality | PairedEndSequencesWithQuality | JoinedSequencesWithQuality] based on Metadata sample ID values; returns FeatureData[Sequence | AlignedSequence]\n",
    "reps = filter_seqs(\n",
    "    data=reps,\n",
    "    table=tabs,\n",
    ").filtered_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f43843a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-22T22:48:07.523952Z",
     "start_time": "2022-09-22T22:48:07.517565Z"
    }
   },
   "source": [
    "## Alpha diversity analysis\n",
    "\n",
    "#### Reference\n",
    "- [The Use and Types of Alpha-Diversity Metrics in Microbial NGS](https://www.cd-genomics.com/microbioseq/the-use-and-types-of-alpha-diversity-metrics-in-microbial-ngs.html)\n",
    "- [Alpha diversity metrics](http://scikit-bio.org/docs/0.2.0/generated/skbio.diversity.alpha.html)\n",
    "\n",
    "#### Methods\n",
    "- [diversity alpha](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha/): Computes a user-specified alpha diversity metric for all samples in a\n",
    "feature table.\n",
    "- [diversity alpha_phylogenetic](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-phylogenetic/): Computes a user-specified phylogenetic alpha diversity metric for all\n",
    "samples in a feature table.\n",
    "- [diversity alpha_correlation](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-correlation/): Determine whether numeric sample metadata columns are correlated with alpha diversity.\n",
    "- [diversity alpha_group_significance](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-group-significance/): Visually and statistically compare groups of alpha diversity values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "606b9c75",
   "metadata": {},
   "source": [
    "### Compute Alpha Diversity vectors\n",
    "- [diversity alpha](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha/): Computes a user-specified alpha diversity metric for all samples in a feature table.\n",
    "- [Alpha diversity metrics](http://scikit-bio.org/docs/0.2.0/generated/skbio.diversity.alpha.html)\n",
    " - Choices: ('ace', 'berger_parker_d', 'brillouin_d', 'chao1', 'chao1_ci', 'dominance', 'doubles', 'enspie', 'esty_ci', 'fisher_alpha', 'gini_index', 'goods_coverage', 'heip_e', 'kempton_taylor_q', 'lladser_pe', 'margalef', 'mcintosh_d', 'mcintosh_e', 'menhinick', 'michaelis_menten_fit', 'observed_features', 'osd', 'pielou_e', 'robbins', 'shannon', 'simpson', 'simpson_e', 'singles', 'strong')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "556e5264",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:20.790166Z",
     "start_time": "2022-09-29T20:54:20.451134Z"
    }
   },
   "outputs": [],
   "source": [
    "metrics = ('ace', 'berger_parker_d', 'brillouin_d', 'chao1', 'chao1_ci', 'dominance', 'doubles', 'enspie', 'esty_ci', 'fisher_alpha', 'gini_index', 'goods_coverage', 'heip_e', 'kempton_taylor_q', 'lladser_pe', 'margalef', 'mcintosh_d', 'mcintosh_e', 'menhinick', 'michaelis_menten_fit', 'observed_features', 'osd', 'pielou_e', 'robbins', 'shannon', 'simpson', 'simpson_e', 'singles', 'strong')\n",
    "\n",
    "# Sugestão de valores para Alpha diversity\n",
    "# chao1 e observed_features (riqueza); shannon e simpson (diversidade - que levam em consideração riqueza e equitabilidade).\n",
    "metrics = ('chao1', 'chao1_ci', 'observed_features', 'shannon', 'simpson', 'simpson_e')\n",
    "alpha_diversities = dict()\n",
    "for metric in metrics:\n",
    "    print(f\"Calculating alpha diversity: {metric}\")\n",
    "    try:\n",
    "        alpha_diversity = alpha(table=tabs, metric=metric).alpha_diversity\n",
    "        alpha_diversities[metric] = alpha_diversity\n",
    "        # Save SampleData[AlphaDiversity] Artifact\n",
    "        file_path = os.path.join(alpha_path, f'alpha-values-{metric}.qza')\n",
    "        alpha_diversity.save(file_path)\n",
    "        print(f\"DONE: Calculating alpha diversity: {metric}\")\n",
    "    except Exception as e:\n",
    "        print(f\"ERROR: Calculating alpha diversity: {metric}\")\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "184bcb7b",
   "metadata": {},
   "source": [
    "### Create Phylogenetic inference\n",
    "\n",
    "- [alignment align_to_tree_mafft_fasttree](https://docs.qiime2.org/2022.8/plugins/available/phylogeny/align-to-tree-mafft-fasttree/): Build a phylogenetic tree using fasttree and mafft alignment\n",
    "\n",
    "This pipeline will start by creating a sequence alignment using MAFFT,\n",
    "after which any alignment columns that are phylogenetically uninformative\n",
    "or ambiguously aligned will be removed (masked). The resulting masked\n",
    "alignment will be used to infer a phylogenetic tree and then subsequently\n",
    "rooted at its midpoint. Output files from each step of the pipeline will be\n",
    "saved. This includes both the unmasked and masked MAFFT alignment from\n",
    "q2-alignment methods, and both the rooted and unrooted phylogenies from\n",
    "q2-phylogeny methods.\n",
    "\n",
    "\n",
    "Returns\n",
    "- alignment : FeatureData[AlignedSequence] : The aligned sequences.\n",
    "- masked_alignment : FeatureData[AlignedSequence] : The masked alignment.\n",
    "- tree : Phylogeny[Unrooted] : The unrooted phylogenetic tree.\n",
    "- rooted_tree : Phylogeny[Rooted] : The rooted phylogenetic tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c8cfe47",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:35.139531Z",
     "start_time": "2022-09-29T20:54:20.792130Z"
    }
   },
   "outputs": [],
   "source": [
    "mafft_alignment, mafft_masked_alignment, mafft_tree, mafft_rooted_tree = align_to_tree_mafft_fasttree(\n",
    "    sequences=reps, n_threads=6, )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0768d598",
   "metadata": {},
   "source": [
    "### Compute Alpha Diversity (Phylogeny)\n",
    "- [diversity alpha_phylogenetic](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-phylogenetic/): Computes a user-specified phylogenetic alpha diversity metric for all samples in a feature table.\n",
    "- Metrics: Choices ('faith_pd')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0383aa19",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:35.310009Z",
     "start_time": "2022-09-29T20:54:35.141333Z"
    }
   },
   "outputs": [],
   "source": [
    "metrics = ('faith_pd', )\n",
    "alpha_diversities_phylogenetic = dict()\n",
    "for metric in metrics:\n",
    "    print(f\"Calculating alpha diversity: {metric}\")\n",
    "    try:\n",
    "        alpha_diversity = alpha_phylogenetic(table=tabs, phylogeny=mafft_rooted_tree, metric=metric).alpha_diversity\n",
    "        alpha_diversities_phylogenetic[metric] = alpha_diversity\n",
    "        # Save Artifact\n",
    "        file_path = os.path.join(alpha_path, f'alpha-phylogeny-{metric}.qza')\n",
    "        alpha_diversity.save(file_path)\n",
    "        print(f\"DONE: Calculating alpha phylogeny: {metric}\")\n",
    "    except Exception as e:\n",
    "        print(f\"ERROR: Calculating alpha phylogeny: {metric}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42db1ead",
   "metadata": {},
   "source": [
    "- [core-metrics-phylogenetic](https://docs.qiime2.org/2023.7/plugins/available/diversity/core-metrics-phylogenetic/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef726ab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_depth = int(tabs.view(pd.DataFrame).sum(axis=1).min())\n",
    "results = core_metrics_phylogenetic(\n",
    "    table = tabs,\n",
    "    phylogeny = mafft_rooted_tree,\n",
    "    sampling_depth = s_depth,\n",
    "    metadata = metadata_qa,\n",
    "    n_jobs_or_threads = 6,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19311218",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_info = [(\"rarefied_table\", \"FeatureTable[Frequency]\", \"The resulting rarefied feature table.\"),\n",
    "(\"faith_pd_vector\", \"SampleData[AlphaDiversity]\", \"Vector of Faith PD values by sample.\"),\n",
    "(\"observed_features_vector\", \"SampleData[AlphaDiversity]\", \"Vector of Observed Features values by sample.\"),\n",
    "(\"shannon_vector\", \"SampleData[AlphaDiversity]\", \"Vector of Shannon diversity values by sample.\"),\n",
    "(\"evenness_vector\", \"SampleData[AlphaDiversity]\", \"Vector of Pielou's evenness values by sample.\"),\n",
    "(\"unweighted_unifrac_distance_matrix\", \"DistanceMatrix\", \"Matrix of unweighted UniFrac distances between pairs of samples.\"),\n",
    "(\"weighted_unifrac_distance_matrix\", \"DistanceMatrix\", \"Matrix of weighted UniFrac distances between pairs of samples.\"),\n",
    "(\"jaccard_distance_matrix\", \"DistanceMatrix\", \"Matrix of Jaccard distances between pairs of samples.\"),\n",
    "(\"bray_curtis_distance_matrix\", \"DistanceMatrix\", \"Matrix of Bray-Curtis distances between pairs of samples.\"),\n",
    "(\"unweighted_unifrac_pcoa_results\", \"PCoAResults\", \"PCoA matrix computed from unweighted UniFrac distances between samples.\"),\n",
    "(\"weighted_unifrac_pcoa_results\", \"PCoAResults\", \"PCoA matrix computed from weighted UniFrac distances between samples.\"),\n",
    "(\"jaccard_pcoa_results\", \"PCoAResults\", \"PCoA matrix computed from Jaccard distances between samples.\"),\n",
    "(\"bray_curtis_pcoa_results\", \"PCoAResults\", \"PCoA matrix computed from Bray-Curtis distances between samples.\"),\n",
    "(\"unweighted_unifrac_emperor\", \"Visualization\", \"Emperor plot of the PCoA matrix computed from unweighted UniFrac.\"),\n",
    "(\"weighted_unifrac_emperor\", \"Visualization\", \"Emperor plot of the PCoA matrix computed from weighted UniFrac.\"),\n",
    "(\"jaccard_emperor\", \"Visualization\", \"Emperor plot of the PCoA matrix computed from Jaccard.\"),\n",
    "(\"bray_curtis_emperor\", \"Visualization\", \"Emperor plot of the PCoA matrix computed from Bray-Curtis.\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2264925a",
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_matrix = dict()\n",
    "for i, info in enumerate(results_info):\n",
    "    r_id, r_type, r_desc = info\n",
    "    #print(i, r_id, r_type)\n",
    "    file_name = f\"{r_id}.qzv\"\n",
    "    if r_type == \"FeatureTable[Frequency]\":\n",
    "        pass\n",
    "    elif r_type == \"DistanceMatrix\":\n",
    "        distance_matrix[r_id] = results[i]\n",
    "    elif r_id.endswith('emperor'):\n",
    "        print(i, r_id, r_type)\n",
    "        print(f\"--- {r_desc} ---\")\n",
    "        file_name = os.path.join(beta_path, file_name)\n",
    "        print(f'Saving emperor file at: {file_name}\\n')\n",
    "        results[i].save(filepath=file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7752cfd",
   "metadata": {},
   "source": [
    "### Alpha diversity correlation\n",
    "\n",
    "This method only process `numeric` columns.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2814f613",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:35.317618Z",
     "start_time": "2022-09-29T20:54:35.312094Z"
    }
   },
   "outputs": [],
   "source": [
    "methods = ('spearman', 'pearson')\n",
    "numerics_cols = metadata_qa.filter_columns(column_type='numeric')\n",
    "if numerics_cols.column_count > 0:\n",
    "    for metric, alpha_values in alpha_diversities.items():\n",
    "        for method in methods:\n",
    "            try:\n",
    "                corr_view = alpha_correlation(alpha_diversity=alpha_values, metadata=numerics_cols, \n",
    "                                          method=method, intersect_ids=True).visualization\n",
    "                view_path = os.path.join(alpha_path, f'alpha-correlation-{metric}-{method}.qzv')\n",
    "                corr_view.save(view_path)\n",
    "                corr_view\n",
    "                print(f\"DONE: Calculating alpha correlation: {metric} {method}\")\n",
    "            except Exception as e:\n",
    "                print(f\"ERROR: Calculating alpha correlation: {metric} {method}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ea83c8a",
   "metadata": {},
   "source": [
    "## Alpha diversity comparisons\n",
    "\n",
    "Visually and statistically compare groups of alpha diversity values.\n",
    "\n",
    "[diversity alpha_group_significance](https://docs.qiime2.org/2022.8/plugins/available/diversity/alpha-group-significance/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "276d1c18",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:36.089290Z",
     "start_time": "2022-09-29T20:54:35.318963Z"
    }
   },
   "outputs": [],
   "source": [
    "for metric, alpha_values in alpha_diversities.items():\n",
    "    print(f\"Processing alpha_group_significance: {metric}\")\n",
    "    try:\n",
    "        significance_view = alpha_group_significance(alpha_diversity=alpha_values, metadata=metadata_qa).visualization\n",
    "        view_path = os.path.join(alpha_path, f'alpha-group-significance-{metric}.qzv')\n",
    "        significance_view.save(view_path)\n",
    "        significance_view\n",
    "        print(f\"DONE: Calculating alpha group significance: {metric}\")\n",
    "    except Exception as e:\n",
    "        print(f\"ERROR: Calculating alpha group significance: {metric}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4349d357",
   "metadata": {},
   "source": [
    "## Beta diversity analysis\n",
    "\n",
    "#### Reference\n",
    "- [diversity beta](https://docs.qiime2.org/2022.8/plugins/available/diversity/beta/): Computes a user-specified beta diversity metric for all pairs of samples in a feature table.\n",
    "- [Beta diversity metrics](http://scikit-bio.org/docs/0.2.0/generated/skbio.diversity.beta.html)\n",
    "\n",
    "- Metric Choices('aitchison', 'braycurtis', 'canberra', 'canberra_adkins', 'chebyshev', 'cityblock', 'correlation', 'cosine', 'dice', 'euclidean', 'hamming', 'jaccard', 'jensenshannon', 'kulsinski', 'matching', 'minkowski', 'rogerstanimoto', 'russellrao', 'seuclidean', 'sokalmichener', 'sokalsneath', 'sqeuclidean', 'yule')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f358dd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:54:36.335038Z",
     "start_time": "2022-09-29T20:54:36.091288Z"
    }
   },
   "outputs": [],
   "source": [
    "metrics = ('aitchison', 'braycurtis', 'canberra', 'canberra_adkins', 'chebyshev', 'cityblock', 'correlation', 'cosine', 'dice', 'euclidean', 'hamming', 'jaccard', 'jensenshannon', 'kulsinski', 'matching', 'minkowski', 'rogerstanimoto', 'russellrao', 'seuclidean', 'sokalmichener', 'sokalsneath', 'sqeuclidean', 'yule')\n",
    "metrics = ('euclidean', 'dice', 'braycurtis', 'correlation', 'cosine', 'matching', 'jaccard')\n",
    "beta_diversities = dict()\n",
    "for metric in metrics:\n",
    "    print(f\"Calculating beta diversity: {metric}\")\n",
    "    try:\n",
    "        beta_diversity = beta(table=tabs, metric=metric, n_jobs=6, pseudocount=1).distance_matrix\n",
    "        beta_diversities[metric] = beta_diversity\n",
    "        # Save SampleData[BetaDiversity] Artifact\n",
    "        file_path = os.path.join(beta_path, f'beta-values-{metric}.qza')\n",
    "        beta_diversity.save(file_path)\n",
    "        print(f\"DONE: Calculating beta diversity: {metric}\")\n",
    "    except Exception as e:\n",
    "        print(f\"ERROR: Calculating beta diversity: {metric}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "200a2fa7",
   "metadata": {},
   "source": [
    "### Beta group significance\n",
    "\n",
    "- [diversity beta_group_significance](https://docs.qiime2.org/2022.8/plugins/available/diversity/beta-group-significance/): Determine whether groups of samples are significantly different from one another using a permutation-based statistical test.\n",
    "- Marti J Anderson. A new method for non-parametric multivariate analysis of variance. Austral ecology, 26(1):32–46, 2001. doi:https://doi.org/10.1111/j.1442-9993.2001.01070.pp.x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8dd42b6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-09-29T20:55:08.313860Z",
     "start_time": "2022-09-29T20:54:36.336993Z"
    }
   },
   "outputs": [],
   "source": [
    "methods = ('permanova', 'anosim', 'permdisp')\n",
    "for method in methods:\n",
    "    for metric, beta_diversity in beta_diversities.items():\n",
    "        print(f'Calculating beta group significance with method {method} and metric {metric}')\n",
    "        try:\n",
    "            beta_view = beta_group_significance(distance_matrix=beta_diversity, \n",
    "                                                metadata=metadata_qa.get_column(class_col), \n",
    "                                                pairwise=True, method=method).visualization\n",
    "            view_name = os.path.join(beta_path, f'beta-group-significance-{metric}-{method}.qzv')\n",
    "            beta_view.save(view_name)\n",
    "            print(f\"DONE: Calculating beta group significance: {method} {metric}\")\n",
    "        except Exception as e:\n",
    "            print(f\"ERROR: Calculating beta group significance: {method} {metric}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6366725e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Expand tests using UNIFRAC metrics\n",
    "methods = ('permanova', 'anosim', 'permdisp')\n",
    "for method in methods:\n",
    "    for metric, beta_diversity in distance_matrix.items():\n",
    "        print(f'Calculating beta group significance with method {method} and metric {metric}')\n",
    "        try:\n",
    "            beta_view = beta_group_significance(distance_matrix=beta_diversity, \n",
    "                                                metadata=metadata_qa.get_column(class_col), \n",
    "                                                pairwise=True, method=method).visualization\n",
    "            view_name = os.path.join(beta_path, f'beta-group-significance-{metric}-{method}.qzv')\n",
    "            beta_view.save(view_name)\n",
    "            print(f\"DONE: Calculating beta group significance: {method} {metric}\")\n",
    "        except Exception as e:\n",
    "            print(f\"ERROR: Calculating beta group significance: {method} {metric}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0de85cf",
   "metadata": {},
   "source": [
    "### Beta group Rarefaction\n",
    "\n",
    "- [diversity beta_rarefaction](https://docs.qiime2.org/2022.8/plugins/available/diversity/beta-rarefaction/): Repeatedly rarefy a feature table to compare beta diversity results within a given rarefaction depth.  For a given beta diversity metric, this visualizer will provide: an Emperor jackknifed PCoA plot, samples clustered by UPGMA or neighbor joining with support calculation, and a heatmap showing the correlation between rarefaction trials of that beta diversity metric."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}