Submodule03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Submodule 3: Taxonomic Classification and Diversity Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Overview\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objectives\n",
    "+ Perform taxonomic classification of 16S rRNA sequences\n",
    "+ Evaluate and explore rarefaction curves to assess sampling depth and species richness\n",
    "+ Analyze alpha and beta diversity within microbial communities\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Import Data\n",
    "\n",
    "### Files\n",
    " \n",
    " 1. **Food Frequency Questionnaire (FFQ):**\n",
    "    - This CSV file contains data from a food frequency questionnaire completed by study participants. It includes columns for participant identifiers (e.g., SampleID), demographic information (e.g., SEX, AGE), and food consumption frequencies and quantities (e.g., BREAKFASTSANDWICHFREQ, EGGSFREQ). This data will allow us to analyze dietary habits and link them with microbiome profiles.\n",
    " 2. **Lifestyle Questionnaire:**\n",
    "    - This text file contains responses to lifestyle questions, such as physical activity levels, smoking status, and other lifestyle factors. It has a row for each question, an answer key, and responses for individual participants in columns.\n",
    "   \n",
    " 3.  **FASTQ Files:**\n",
    "    - These files contain raw sequencing data from 16S rRNA sequencing. Each read in a FASTQ file has a sequence identifier, nucleotide sequence, and quality scores. \n",
    "\n",
    "Let's start by loading the survey files into DataFrames and performing some initial quality checks to ensure data integrity.\n",
    "\n",
    "## Import Survey Data\n",
    "\n",
    "### Step 1: Load the Survey Data\n",
    "Let's load both survey files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Install and load required packages\n",
    "if (!requireNamespace(\"dplyr\", quietly = TRUE)) {\n",
    "    install.packages(\"dplyr\")\n",
    "}\n",
    "\n",
    "if (!requireNamespace(\"ggplot2\", quietly = TRUE)) {\n",
    "    install.packages(\"ggplot2\")\n",
    "}\n",
    "\n",
    "if (!requireNamespace(\"vegan\", quietly = TRUE)) {\n",
    "    install.packages(\"vegan\")\n",
    "}\n",
    "\n",
    "if (!requireNamespace(\"ggpubr\", quietly = TRUE)) {\n",
    "    install.packages(\"ggpubr\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(\"dplyr\")\n",
    "library(\"ggplot2\")\n",
    "library(\"vegan\")\n",
    "library(\"ggpubr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#Load in the data\n",
    "ffq <- read.csv('FFQ_Data.csv')\n",
    "lsq <- read.csv(\"MetaData.txt\", sep = '\\t')\n",
    "# Display the first few rows of the dataframe to understand its structure\n",
    "head(ffq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our Food Frequency Questionaire has columns representing the Sample ID, Respondent ID number, and all the questions asked. The heading of the columns represents abbreviations for the questions. Now let's look at the dimensions of our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Full dimensions\n",
    "dim(ffq)\n",
    "# Number of row\n",
    "nrow(ffq)\n",
    "# Number of columns\n",
    "ncol(ffq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we have 100 participants (every row is a sample) and 1074 columns which means over 1000 food questions! Now let's look at the lifestyle survey. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Display the first few rows of the dataframe to understand its structure\n",
    "head(lsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our Lifestyle Questionaire has columns representing the Sample ID and all the lifestyle questions asked. The heading of the columns represents abbreviations for the questions. Now let's look at the dimensions of our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Full dimensions\n",
    "dim(lsq)\n",
    "# Number of row\n",
    "nrow(lsq)\n",
    "# Number of columns\n",
    "ncol(lsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have less questions in this survey, only 33, and 96 participants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Quick Quality Control Checks\n",
    "Let's perform several initial checks to assess data quality:\n",
    "\n",
    " 1. **Collect Intersection of Samples:** We need to find the Sample IDs in the FFQ data rows and in the LSQ data columns. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Find sample intersection\n",
    "samples <- intersect(ffq$SampleID, lsq$SampleID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "length(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " 2. **Basic Summary Statistics:** We can get an overview of the data, which can help identify outliers or unexpected values. We will look at the demographic questions in the FFQ and Lifestyle Questionnaire separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Summary statistics for the first 10 columns\n",
    "summary(ffq[, 1:9])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Summary statistics for 3 columns\n",
    "summary(lsq[, 3:5])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. **Visualize Data:** Ensure the data looks accurate with exploratory visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Load necessary library\n",
    "library(ggplot2)\n",
    "\n",
    "# Set up the plotting theme\n",
    "theme_set(theme_minimal())\n",
    "\n",
    "# 1. Bar plot for Sex distribution\n",
    "# Bar plot for Sex distribution without modifying the original data\n",
    "ggplot(ffq, aes(x = factor(SEX, levels = c(1, 2), labels = c(\"Male\", \"Female\")))) +\n",
    "  geom_bar(fill = \"#0073C2FF\") +  # Blue color for bars\n",
    "  labs(title = \"Distribution of Sex\", x = \"Sex\", y = \"Count\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "    axis.text = element_text(size = 18))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(lsq, aes(x = factor(Sex))) +\n",
    "  geom_bar(fill = \"#0073C2FF\") +  # Blue color for bars\n",
    "  labs(title = \"Distribution of Sex\", x = \"Sex\", y = \"Count\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "        axis.text = element_text(size = 18))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# 2. Histogram for Age distribution\n",
    "# Dropping NA values for AGE\n",
    "ggplot(ffq, aes(x = AGE)) +\n",
    "  geom_histogram(bins = 10, fill = \"lightgreen\", color = \"black\") +\n",
    "  labs(title = \"Age Distribution\", x = \"Age\", y = \"Frequency\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "    axis.text = element_text(size = 18))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(lsq, aes(x = Age)) +\n",
    "  geom_histogram(bins = 10, fill = \"lightgreen\", color = \"black\") +\n",
    "  labs(title = \"Age Distribution\", x = \"Age\", y = \"Frequency\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "        axis.text = element_text(size = 18))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# 3. Histogram for Weight distribution\n",
    "# Dropping NA values for WEIGHT\n",
    "ggplot(ffq, aes(x = WEIGHT)) +\n",
    "  geom_histogram(bins = 15, fill = \"lightcoral\", color = \"black\") +\n",
    "  labs(title = \"Weight Distribution\", x = \"Weight\", y = \"Frequency\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "    axis.text = element_text(size = 18))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(lsq, aes(x = Weight)) +\n",
    "  geom_histogram(bins = 15, fill = \"lightcoral\", color = \"black\") +\n",
    "  labs(title = \"Weight Distribution\", x = \"Weight\", y = \"Frequency\") +\n",
    "  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size = 20),\n",
    "        axis.text = element_text(size = 18))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Combine Metadata\n",
    "\n",
    "We have two different surveys, let's combine them into one so we only have one metadata file to work with. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Load in data from Submodule 2\n",
    "sample.names <- readRDS('sampleNames.rds')\n",
    "taxa         <- readRDS(\"taxa.rds\")\n",
    "reads        <- readRDS(\"seqtabnochim.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Merge data\n",
    "metadata <- merge(ffq, lsq, by = \"SampleID\", all = TRUE)\n",
    "# Rename samples to match FastQ files\n",
    "metadata$SampleID <- gsub(\"WP_\", \"WP-\", metadata$SampleID)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a lot more samples in our metadata than we do sequenced, therefore we need to pull out only the samples who have sequenced reads from our metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Add in only if they have sequences\n",
    "dfseqIndo <- data.frame(sampleNames = sample.names,\n",
    "                       inSeq       = rep('yes', length(sample.names)))\n",
    "# Only keep those who do\n",
    "metadata <- merge(dfseqIndo, metadata, by.x = 'sampleNames', by.y = 'SampleID', all.x = TRUE)\n",
    "colnames(metadata)[1] <- 'SampleID'\n",
    "rownames(metadata)    <- metadata$SampleID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "dim(metadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have 42 samples with metadata and sequenced reads.\n",
    "\n",
    "When merging metadata from two different surveys, it is common to encounter discrepancies due to human error, such as participants providing slightly different responses across surveys for variables like age or weight. While these inconsistencies highlight the importance of careful data handling, for this analysis, we will rely on the FFQ data as it is more standardized and consistent. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Flag discrepancies and missing data\n",
    "discrepancies <- metadata %>%\n",
    "  mutate(\n",
    "    age_discrepancy = ifelse(Age != AGE, TRUE, FALSE),\n",
    "    weight_discrepancy = ifelse(Weight != WEIGHT, TRUE, FALSE)\n",
    "  )\n",
    "\n",
    "\n",
    "# Summarize discrepancies and missing data\n",
    "summary <- data.frame(\n",
    "  age_discrepancy = sum(discrepancies$age_discrepancy, na.rm = TRUE),\n",
    "  weight_discrepancy = sum(discrepancies$weight_discrepancy, na.rm = TRUE)\n",
    ")\n",
    "\n",
    "summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data shows some inconsistencies. We will use the FFQ data, as it is a standardized survey. Now, let's examine our data more by creating a table to get an overview of our data by sex and BMI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Order BMI\n",
    "metadata$BMI_ordinal = factor(metadata$BMI_ordinal, levels = c('Healthy weight', 'Overweight', 'Obese'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Group data by Sex\n",
    "table1 <- metadata %>%\n",
    "  group_by(Sex) %>%\n",
    "  summarize(\n",
    "    Count = n(),\n",
    "    Mean_Age = mean(AGE, na.rm = TRUE),\n",
    "    Mean_Weight = mean(WEIGHT, na.rm = TRUE),\n",
    "    Mean_Height = mean(Height, na.rm = TRUE),\n",
    "    Mean_BMI = mean(BMI, na.rm = TRUE),\n",
    "    Mean_Num_Animals = mean(num_animals, na.rm = TRUE),\n",
    "    Mean_Kcal = mean(DT_KCAL, na.rm = TRUE),\n",
    "    Mean_Protein = mean(DT_PROT, na.rm = TRUE),\n",
    "    Mean_Carbohydrate = mean(DT_CARB, na.rm = TRUE),\n",
    "    Mean_Fat = mean(DT_TFAT, na.rm = TRUE),\n",
    "    Mean_Alcohol = mean(DT_ALCO, na.rm = TRUE),\n",
    "    Mean_Sugars = mean(DT_SUG_T, na.rm = TRUE),\n",
    "    Mean_Fiber = mean(DT_FIBE, na.rm = TRUE)\n",
    "  )\n",
    "\n",
    "# View the table\n",
    "table1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Group data by BMI\n",
    "table2 <- metadata %>%\n",
    "  group_by(BMI_ordinal) %>%\n",
    "  summarize(\n",
    "    Count = n(),\n",
    "    Mean_Age = mean(AGE, na.rm = TRUE),\n",
    "    Mean_Weight = mean(WEIGHT, na.rm = TRUE),\n",
    "    Mean_Height = mean(Height, na.rm = TRUE),\n",
    "    Male_Percent = mean(SEX == 1, na.rm = TRUE) * 100,\n",
    "    Female_Percent = mean(SEX == 2, na.rm = TRUE) * 100,\n",
    "    Mean_Num_Animals = mean(num_animals, na.rm = TRUE),\n",
    "    Mean_Kcal = mean(DT_KCAL, na.rm = TRUE),\n",
    "    Mean_Protein = mean(DT_PROT, na.rm = TRUE),\n",
    "    Mean_Carbohydrate = mean(DT_CARB, na.rm = TRUE),\n",
    "    Mean_Fat = mean(DT_TFAT, na.rm = TRUE),\n",
    "    Mean_Alcohol = mean(DT_ALCO, na.rm = TRUE),\n",
    "    Mean_Sugars = mean(DT_SUG_T, na.rm = TRUE),\n",
    "    Mean_Fiber = mean(DT_FIBE, na.rm = TRUE)\n",
    "  )\n",
    "\n",
    "# View the table\n",
    "table2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now categorize dietary intake so we can eaisly view differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Create metadata variables with levels based on dietary standards\n",
    "metadata <- metadata %>%\n",
    "  mutate(\n",
    "    # Calories (kcal/day, general range: 1500-3000 for adults)\n",
    "    Kcal_Level = case_when(\n",
    "      DT_KCAL < 1500 ~ \"Low\",\n",
    "      DT_KCAL >= 1500 & DT_KCAL <= 3000 ~ \"Moderate\",\n",
    "      DT_KCAL > 3000 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    # Protein (g/day, general range: 46-56 g for adults)\n",
    "    Protein_Level = case_when(\n",
    "      DT_PROT < 46 ~ \"Low\",\n",
    "      DT_PROT >= 46 & DT_PROT <= 56 ~ \"Moderate\",\n",
    "      DT_PROT > 56 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    # Carbohydrates (g/day, general range: 225-325 g based on 2000 kcal diet)\n",
    "    Carb_Level = case_when(\n",
    "      DT_CARB < 225 ~ \"Low\",\n",
    "      DT_CARB >= 225 & DT_CARB <= 325 ~ \"Moderate\",\n",
    "      DT_CARB > 325 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    # Total Fat (g/day, general range: 44-77 g based on 2000 kcal diet)\n",
    "    Fat_Level = case_when(\n",
    "      DT_TFAT < 44 ~ \"Low\",\n",
    "      DT_TFAT >= 44 & DT_TFAT <= 77 ~ \"Moderate\",\n",
    "      DT_TFAT > 77 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    # Dietary Fiber (g/day, general recommendation: 25-38 g)\n",
    "    Fiber_Level = case_when(\n",
    "      DT_FIBE < 25 ~ \"Low\",\n",
    "      DT_FIBE >= 25 & DT_FIBE <= 38 ~ \"Moderate\",\n",
    "      DT_FIBE > 38 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    \n",
    "    # Categorize alcohol intake\n",
    "    Alcohol_Level = case_when(\n",
    "      DT_ALCO == 0 ~ \"No Alcohol\",\n",
    "      DT_ALCO > 0 & DT_ALCO <= 14 ~ \"Low\",\n",
    "      DT_ALCO > 14 & DT_ALCO <= 30 ~ \"Moderate\",\n",
    "      DT_ALCO > 30 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    \n",
    "    # Categorize sugar intake\n",
    "    Sugar_Level = case_when(\n",
    "      DT_SUG_T <= 25 ~ \"Low\",\n",
    "      DT_SUG_T > 25 & DT_SUG_T <= 50 ~ \"Moderate\",\n",
    "      DT_SUG_T > 50 ~ \"High\",\n",
    "      TRUE ~ NA_character_\n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "# Reorder Fiber_Level\n",
    "metadata$Fiber_Level <- factor(metadata$Fiber_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Reorder Alcohol_Level\n",
    "metadata$Alcohol_Level <- factor(metadata$Alcohol_Level, levels = c(\"No Alcohol\", \"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Reorder Sugar_Level\n",
    "metadata$Sugar_Level <- factor(metadata$Sugar_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "\n",
    "# Reorder Kcal_Level\n",
    "metadata$Kcal_Level <- factor(metadata$Kcal_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Reorder Protein_Level\n",
    "metadata$Protein_Level <- factor(metadata$Protein_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Reorder Carbohydrate_Level\n",
    "metadata$Carb_Level <- factor(metadata$Carb_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Reorder Fat_Level\n",
    "metadata$Fat_Level <- factor(metadata$Fat_Level, levels = c(\"Low\", \"Moderate\", \"High\"))\n",
    "\n",
    "# Treat the number of animals as discrete not continuous\n",
    "metadata$num_animals = factor(metadata$num_animals, levels = c('0', '1', '2', '3', '4'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view the distriibution of our `Protien_Level`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(metadata, aes(x = Protein_Level, fill = Protein_Level)) +\n",
    "  geom_bar() +\n",
    "  theme_classic() +\n",
    "  labs(\n",
    "    title = \"Distribution of Protein Levels\",\n",
    "    x = \"Protein Level\",\n",
    "    y = \"Count\"\n",
    "  ) +\n",
    "  scale_fill_manual(values = c('Low' = '#f7b801',  \n",
    "                                'Moderate' = '#f35b04', \n",
    "                                'High' = '#52006A')) +\n",
    "  theme(\n",
    "    axis.text.x = element_text(angle = 45, hjust = 1),\n",
    "    legend.position = \"none\"\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view protein intake by sex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(metadata, aes(x = Protein_Level, fill = Protein_Level)) +\n",
    "  geom_bar() +\n",
    "  theme_classic() +\n",
    "  labs(\n",
    "    title = \"Distribution of Protein Levels\",\n",
    "    x = \"Protein Level\",\n",
    "    y = \"Count\"\n",
    "  ) +\n",
    "  scale_fill_manual(values = c('Low' = '#f7b801',  \n",
    "                                'Moderate' = '#f35b04', \n",
    "                                'High' = '#52006A')) +\n",
    "  facet_wrap(~SEX, labeller = as_labeller(c(\"1\" = \"Male\", \"2\" = \"Female\"))) +\n",
    "  theme(\n",
    "    axis.text.x = element_text(angle = 45, hjust = 1),\n",
    "    legend.position = \"none\"\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now by BMI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(metadata, aes(x = Protein_Level, fill = Protein_Level)) +\n",
    "  geom_bar() +\n",
    "  theme_classic() +\n",
    "  labs(\n",
    "    title = \"Distribution of Protein Levels\",\n",
    "    x = \"Protein Level\",\n",
    "    y = \"Count\"\n",
    "  ) +\n",
    "  scale_fill_manual(values = c('Low' = '#f7b801',  \n",
    "                                'Moderate' = '#f35b04', \n",
    "                                'High' = '#52006A')) +\n",
    "  facet_wrap(~BMI_ordinal) +\n",
    "  theme(\n",
    "    axis.text.x = element_text(angle = 45, hjust = 1),\n",
    "    legend.position = \"none\"\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Assign Taxonomy to our Sequencing Data\n",
    "\n",
    "We will be using the package `phyloseq` to analyze our microbiomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Install and load libraries\n",
    "if (!requireNamespace(\"phyloseq\", quietly = TRUE)) {\n",
    "    install.packages(\"phyloseq\")\n",
    "}\n",
    "\n",
    "library(phyloseq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Assigning taxonomy\n",
    "ps <- phyloseq(otu_table(reads, taxa_are_rows = F),\n",
    "               tax_table(taxa),\n",
    "               sample_data(metadata))\n",
    "\n",
    "ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "write.table(otu_table(reads, taxa_are_rows = F), \n",
    "            \"ASV_table.txt\", row.names = TRUE, quote = FALSE, sep = '\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Data Analysis and Visualization \n",
    "\n",
    "#### Rarefaction Curves\n",
    "\n",
    "Rarefaction curves are essential tools in microbiome analysis used to assess the richness of microbial communities and the adequacy of sequencing depth. These curves plot the number of observed species against the number of sequencing reads, providing insight into species diversity within a sample.\n",
    "\n",
    "The shape of a rarefaction curve indicates whether sequencing captured the community’s full diversity: a plateau suggests that most species have been identified, while a steep incline indicates that additional sequencing may uncover more species. Rarefaction also enables comparisons across samples with varying sequencing depths, ensuring that observed differences in diversity are not artifacts of unequal sampling effort.\n",
    "\n",
    "Let's start by looking at the differences between sex and protein intake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rare <- rarecurve(reads,step = 1000, tidy = TRUE)\n",
    "rare <- inner_join(rare, metadata, by=c(\"Site\"=\"SampleID\"))\n",
    "\n",
    "ggplot(rare, aes(x=Sample, y=Species, col=Protein_Level, linetype=Sex, group=Site))+\n",
    "  geom_line(linewidth = 2, alpha = 0.8)+\n",
    "  scale_color_manual(values = c('Low' = '#f7b801',  \n",
    "                                'Moderate' = '#f35b04', \n",
    "                                'High' = '#52006A')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"Sample Size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot BMI and protein intake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(rare, aes(x=Sample, y=Species, col=Protein_Level, linetype=BMI_ordinal, group=Site))+\n",
    "  geom_line(linewidth = 2, alpha = 0.8) +\n",
    "  scale_color_manual(values = c('Low' = '#f7b801',  \n",
    "                                'Moderate' = '#f35b04', \n",
    "                                'High' = '#52006A')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"Sample Size\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Diversity in Microbiome\n",
    "\n",
    "Alpha and beta diversity are key concepts in microbiome studies used to understand the composition of microbial communities.\n",
    " - **Alpha diversity** refers to the diversity within a single sample, capturing the richness (number of species) and evenness (distribution of species). It provides insights into how diverse a community is at a specific site or condition. Common metrics include species richness, Simpson’s index, and the Shannon index.\n",
    " - **Beta diversity**, on the other hand, measures the variation in microbial communities between samples. It highlights differences in composition and can reveal patterns of microbial distribution across environments or conditions. Common beta diversity metrics include Bray-Curtis dissimilarity and UniFrac distances.\n",
    "\n",
    "##### Alpha Diversity\n",
    "\n",
    "Let's look at the alpha diversity in our samples with the **Shannon index**. The Shannon index measures both the richness and evenness of species within a sample. A higher Shannon index indicates greater diversity, while a value of 0 means only one species is present in the community."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "shannon <- diversity(reads, index = \"shannon\") %>%\n",
    "  as_tibble(rownames = \"SampleID\")%>%\n",
    "  inner_join(., metadata, by=\"SampleID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(shannon, aes(x=Sex, y=value, col=Sex))+\n",
    "  geom_boxplot(outlier.shape = NA)+\n",
    "  geom_jitter(width = 0.2) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Shannon Weaver Index\") +\n",
    "  stat_compare_means(aes(x=Sex, y=value, col=Sex), hide.ns = F, method=\"t.test\",label = \"p.signif\",\n",
    "                     label.x = 1.5, label.y= 4.8, show.legend = F) +\n",
    "  theme(legend.position=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at BMI within sex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Instal and load ggstance for better visualization\n",
    "if (!requireNamespace(\"ggstance\", quietly = TRUE)) {\n",
    "    install.packages(\"ggstance\")\n",
    "}\n",
    "\n",
    "\n",
    "library(ggstance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(shannon, aes(x=BMI_ordinal, y=value, col=BMI_ordinal))+\n",
    "  geom_boxplot()+  \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.5), alpha=0.5) +     \n",
    "  scale_color_manual(values = c('Healthy weight' = '#f7b801', 'Overweight' = '#f18701', \n",
    "                                'Obese' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Shannon Weaver Index\") +\n",
    "  theme(legend.position=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about alpha diversity by protein intake?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(shannon, aes(x=Protein_Level, y=value, col=Protein_Level))+\n",
    "  geom_boxplot()+  \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha=0.5) +     \n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Shannon Weaver Index\") +\n",
    "  theme(legend.position=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at the combination of sex or BMI by protein intake. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(shannon, aes(x=Sex, y=value, col=Protein_Level))+\n",
    "  geom_boxplot()+ \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha=0.5) +     \n",
    "  facet_grid(col=vars(Sex), scales = \"free_x\")+\n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Shannon Weaver Index\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(shannon, aes(x=BMI_ordinal, y=value, col=Protein_Level))+\n",
    "  geom_boxplot()+   \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha=0.5) + \n",
    "  facet_grid(col=vars(BMI_ordinal), scales = \"free_x\")+\n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Shannon Weaver Index\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's look at the alpha diversity with the **Inverse Simpson Index**. The Inverse Simpson Index reflects both species richness and evenness within a community. It accounts for the number of species present and their relative abundances, emphasizing the dominance of common species. A higher Inverse Simpson value indicates greater diversity, as it signifies more species evenly distributed in the community. This index is particularly useful for understanding community composition and detecting dominance patterns in microbiome studies. We will focus on our experimental questions, looking at the combination of variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "invsimp<-diversity(reads, index = \"invsimpson\")%>%\n",
    "  as_tibble(rownames = \"SampleID\")%>%\n",
    "  inner_join(., metadata, by=\"SampleID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(invsimp, aes(x=Sex, y=value, col=Protein_Level))+\n",
    "  geom_boxplot()+ \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha=0.5) +     \n",
    "  facet_grid(col=vars(Sex), scales = \"free_x\")+\n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Inverse Simpson Index\") \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggplot(invsimp, aes(x=BMI_ordinal, y=value, col=Protein_Level))+\n",
    "  geom_boxplot()+ \n",
    "  geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.5), alpha=0.5) +     \n",
    "  facet_grid(col=vars(BMI_ordinal), scales = \"free_x\")+\n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\") +\n",
    "  labs(x=\"\",y=\"Inverse Simpson Index\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Beta Diversity\n",
    "\n",
    "Non-metric Multidimensional Scaling (NMDS) is an ordination method used to visualize differences in community composition between samples, making it a key tool for analyzing beta diversity. It represents pairwise dissimilarities (e.g., Bray-Curtis distances) in a low-dimensional space, preserving the rank order of distances rather than their exact values.\n",
    "\n",
    "The result is a visual representation where samples with similar microbial communities appear closer together, while dissimilar communities are farther apart. NMDS is particularly useful for identifying patterns and grouping structures in complex ecological or microbiome datasets. The quality of the ordination is assessed using a stress value, with lower values indicating a better representation of the original data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Transform data to proportions as appropriate for Bray-Curtis distances\n",
    "ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))\n",
    "ord.nmds.bray <- ordinate(ps.prop, method=\"NMDS\", distance=\"bray\")\n",
    "\n",
    "# Plot beta diversity                                   \n",
    "plot_ordination(ps.prop, ord.nmds.bray, color=\"Sex\", title=\"Bray NMDS\") + theme_pubr(legend=\"right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also view by BMI and protein intake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot_ordination(ps.prop, ord.nmds.bray, color=\"BMI_ordinal\", title=\"Bray NMDS\") + scale_color_manual(values = c('Healthy weight' = '#f7b801', 'Overweight' = '#f18701', \n",
    "                                'Obese' = '#780116')) + theme_pubr(legend=\"right\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot_ordination(ps.prop, ord.nmds.bray, color=\"Protein_Level\", title=\"Bray NMDS\") +\n",
    "scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  theme_pubr(legend=\"right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will explore the distribution of taxa at the phylum level across our samples using three different datasets: the full dataset, the proportion-transformed dataset, and the subset of the top 20 most abundant taxa. Visualizing these datasets will allow us to examine the overall microbial composition, highlight differences in relative abundance, and focus on the most dominant taxa within the community."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Creating top 20 most abundant\n",
    "top20    <- names(sort(taxa_sums(ps), decreasing=TRUE))[1:20]\n",
    "ps.top20 <- prune_taxa(top20, ps.prop)\n",
    "\n",
    "# Plotting by Sex\n",
    "plot_bar(ps, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Sex, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.prop, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Sex, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.top20, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Sex, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Plotting by BMI\n",
    "plot_bar(ps, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~BMI_ordinal, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.prop, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~BMI_ordinal, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.top20, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~BMI_ordinal, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Plotting by protein intake\n",
    "plot_bar(ps, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Protein_Level, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.prop, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Protein_Level, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))\n",
    "plot_bar(ps.top20, x=\"SampleID\", fill=\"Phylum\") + facet_wrap(~Protein_Level, scales=\"free_x\") + theme_bw() +  theme(axis.text.x =element_text(angle=90))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will focus on a single phylum of interest, *Firmicutes*, to examine its distribution across our study participants. By extracting this specific phylum, we can explore how its relative abundance varies across key categories such as sex, BMI levels, and protein intake. This detailed investigation will provide insights into potential associations between this microbial group and individual characteristics, helping us better understand the role this phylum may play in the overall microbiome composition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Subset data\n",
    "ps.phylum <- tax_glom(ps.prop, taxrank = \"Phylum\")\n",
    "dat       <- psmelt(ps.phylum)\n",
    "idx.keep  <- which(dat$Phylum == \"Firmicutes\") \n",
    "Firmicutes      <- dat[idx.keep,]\n",
    "table(Firmicutes$Phylum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "ggplot(Firmicutes)+\n",
    "  geom_col(aes( x=Sample,y=Abundance, fill=Phylum), position=\"dodge\")+\n",
    "  facet_grid(col=vars(Protein_Level), space=\"free_x\",scales=\"free\")+\n",
    "  scale_fill_manual(values =  c('#780116')) +\n",
    "  labs(title = \"Abundance by Protein Level\") +\n",
    "  theme_bw() +\n",
    "  theme(axis.text.x =element_text(angle=90)\n",
    "  ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use a Kruskal-Wallis test to determine if there are different amounts of *Firmicutes* by protein intake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Create comparison groups \n",
    "my_comparisons <- list( c(\"Low\", \"Moderate\"), c(\"High\", \"Low\"), c(\"Moderate\", \"High\") )\n",
    "\n",
    "# Plot\n",
    "ggplot(Firmicutes, aes(x=Protein_Level,y=Abundance, fill = Phylum, color = Protein_Level))+\n",
    "  geom_boxplot(outlier.shape = NA) +\n",
    "  geom_jitter(width = 0.2, size = 2, alpha = 0.5) +\n",
    "  scale_color_manual(values = c('Low' = '#f7b801', 'Moderate' = '#f18701', \n",
    "                                'High' = '#780116')) +\n",
    "  stat_compare_means(comparisons = my_comparisons, label.y = c(1, 1.2, 1.1))+\n",
    "  stat_compare_means(label.y = 1.3) +\n",
    "  labs(title = \"Firmicutes Abundance by Protein Intake\", x = \"Protein Intake\") +\n",
    "  theme_classic() +\n",
    "  guides(color = \"none\") +\n",
    "  theme(legend.position = \"none\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems there is no significant difference in *Firmicutes* abundance between levels of protein intake. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "In this module, we successfully explored taxonomic classification and diversity analysis using microbiome sequencing data. By integrating dietary and lifestyle metadata with sequencing results, we evaluated the diversity within and between microbial communities. We analyzed alpha diversity using indices such as Shannon and Inverse Simpson to understand richness and evenness, and beta diversity through Non-metric Multidimensional Scaling to assess community composition differences. Additionally, we visualized microbial composition at various taxonomic levels and identified significant associations, such as the relationship between protein intake and *Firmicutes* abundance. These analyses provide a foundation for understanding how diet and other factors influence microbiome composition and diversity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clean up\n",
    "Please remember to shut down your Sagemaker instance. <br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<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.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}