02_Diff_analysis/02_Lab_journal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **CHAPTER 2. Differential analysis**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Install conda env and activate it**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "conda env create -f diff_an.yaml\n",
    "```\n",
    "\n",
    "```\n",
    "conda activate diff_an\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 0. Copy kreports from the server**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "scp -r username@host.com:\"/path/to/kreports/folder\" data/kreports/\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 1. Data parsing**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rename files. Delete `_kraken_report` from file names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_folder}\n",
    "%run scripts/rename_files.py data/kreports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install KrakenTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cloning into 'KrakenTools'...\n",
      "remote: Enumerating objects: 360, done.\u001b[K\n",
      "remote: Counting objects: 100% (110/110), done.\u001b[K\n",
      "remote: Compressing objects: 100% (53/53), done.\u001b[K\n",
      "remote: Total 360 (delta 79), reused 69 (delta 57), pack-reused 250\u001b[K\n",
      "Receiving objects: 100% (360/360), 129.37 KiB | 1003.00 KiB/s, done.\n",
      "Resolving deltas: 100% (218/218), done.\n"
     ]
    }
   ],
   "source": [
    "! git clone https://github.com/jenniferlu717/KrakenTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create folder for files groupping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir data/mpa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert kraken reports in `data/kreports` folder to MPA format and place files in `data/mpa` folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_files} {path_to_output_mpa_files}\n",
    "! ./scripts/run_kreport2mpa.sh data/kreports data/mpa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine mpa files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Number of files to parse: 10\n",
      " Number of classifications to write: 17159\n",
      "\t17159 classifications printed\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_mpa_files} {output_file_name}\n",
    "%run KrakenTools/combine_mpa.py -i data/mpa/* -o data/COMBINED.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 1.1. Getting `counts.csv` files**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare folders to place data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mkdir counts\n",
    "mkdir counts/txt\n",
    "mkdir counts/csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse `data/COMBINED.txt` file to counts files on several taxonomic levels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.1. _`Species`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"s__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__Homo_sapiens\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_species.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.2. _`Genus`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"g__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__\" \\\n",
    "| grep -v \"g__Homo\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_genus.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.3. _`Family`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"f__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__\" \\\n",
    "| grep -v \"g__\" \\\n",
    "| grep -v \"f__Hominidae\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_family.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.4. _`Order`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"o__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__\" \\\n",
    "| grep -v \"g__\" \\\n",
    "| grep -v \"f__\" \\\n",
    "| grep -v \"o__Primates\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_order.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.5. _`Class`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"c__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__\" \\\n",
    "| grep -v \"g__\" \\\n",
    "| grep -v \"f__\" \\\n",
    "| grep -v \"o__\" \\\n",
    "| grep -v \"c__Mammalia\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_class.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 1.1.6. _`Phylum`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "grep -E \"p__\" data/COMBINED.txt \\\n",
    "| grep -v \"t__\" \\\n",
    "| grep -v \"s__\" \\\n",
    "| grep -v \"g__\" \\\n",
    "| grep -v \"f__\" \\\n",
    "| grep -v \"o__\" \\\n",
    "| grep -v \"c__\" \\\n",
    "| grep -v \"p__Chordata\" \\\n",
    "| sed \"s/^.*|//g\" \\\n",
    "| sed \"s/SRS[0-9]*-//g\" \\\n",
    "> counts/txt/counts_phylum.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 1.2. Process counts files**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`processing_script` will do:\n",
    "1. Return the 1st line with sample ids from `data/COMBINED.txt`\n",
    "2. Delete '[X]__' and '_' from organisms names. [X] stands for taxonomic id (s - species, g - genus etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_species.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_genus.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_family.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_order.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_class.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_txt_file} {path_to_output_file}\n",
    "%run scripts/processing_script.py data/COMBINED.txt counts/txt/counts_phylum.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 1.3. Convert to `csv` file format**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`convert2csv` will do:\n",
    "1. Take `counts_tax_lvl.txt` as the input\n",
    "2. Give `counts_tax_lvl.csv` as the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_species.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_species.txt counts/csv/counts_species.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_genus.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_genus.txt counts/csv/counts_genus.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_family.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_family.txt counts/csv/counts_family.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_order.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_order.txt counts/csv/counts_order.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_class.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_class.txt counts/csv/counts_class.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data has been successfully converted and saved as 'counts/csv/counts_phylum.csv'.\n"
     ]
    }
   ],
   "source": [
    "# Usage\n",
    "# {path_to_script} {path_to_input_file}\n",
    "%run scripts/convert2csv.py counts/txt/counts_phylum.txt counts/csv/counts_phylum.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 1.4. Create metadata**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "metadata.csv has been created successfully.\n"
     ]
    }
   ],
   "source": [
    "# Define the data\n",
    "data = [\n",
    "    {'sample_id': 'D1', 'Group': 'Vespertilio murinus'},\n",
    "    {'sample_id': 'D2', 'Group': 'Vespertilio murinus'},\n",
    "    {'sample_id': 'D3', 'Group': 'Vespertilio murinus'},\n",
    "    {'sample_id': 'D4', 'Group': 'Vespertilio murinus'},\n",
    "    {'sample_id': 'D5', 'Group': 'Vespertilio murinus'},\n",
    "    {'sample_id': 'P1', 'Group': 'Nyctalus noctula'},\n",
    "    {'sample_id': 'P2', 'Group': 'Nyctalus noctula'},\n",
    "    {'sample_id': 'P3', 'Group': 'Nyctalus noctula'},\n",
    "    {'sample_id': 'P4', 'Group': 'Nyctalus noctula'},\n",
    "    {'sample_id': 'P5', 'Group': 'Nyctalus noctula'}\n",
    "]\n",
    "\n",
    "# Define the CSV file name\n",
    "filename = 'metadata.csv'\n",
    "\n",
    "# Write the data to the CSV file\n",
    "with open(filename, mode='w', newline='') as file:\n",
    "    writer = csv.DictWriter(file, fieldnames=['sample_id', 'Group'])\n",
    "    writer.writeheader()\n",
    "    writer.writerows(data)\n",
    "\n",
    "print(f'{filename} has been created successfully.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 2. Comparative statistics**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 2.1. Differential Microbial Abundance**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MaAsLin2` is the next generation of `MaAsLin` (Microbiome Multivariable Association with Linear Models).\n",
    "\n",
    "`MaAsLin2` is comprehensive R package for efficiently determining multivariable association between clinical metadata and microbial meta-omics features. `MaAsLin2` relies on general linear models to accommodate most modern epidemiological study designs, including cross-sectional and longitudinal, along with a variety of filtering, normalization, and transform methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir MaAsLin2_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.1. _`Species`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_species.csv MaAsLin2_results/species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_species = pd.read_csv('MaAsLin2_results/species/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.2. _`Genus`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_genus.csv MaAsLin2_results/genus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_species = pd.read_csv('MaAsLin2_results/genus/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.3. _`Family`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_family.csv MaAsLin2_results/family"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_family = pd.read_csv('MaAsLin2_results/family/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_family"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.4. _`Order`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_order.csv MaAsLin2_results/order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_order = pd.read_csv('MaAsLin2_results/order/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.5. _`Class`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_class.csv MaAsLin2_results/class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_class = pd.read_csv('MaAsLin2_results/class/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.1.6. _`Phylum`_ level**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/MaAsLin2.R metadata.csv counts/csv/counts_phylum.csv MaAsLin2_results/phylum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>feature</th>\n",
       "      <th>metadata</th>\n",
       "      <th>value</th>\n",
       "      <th>coef</th>\n",
       "      <th>stderr</th>\n",
       "      <th>N</th>\n",
       "      <th>N.not.0</th>\n",
       "      <th>pval</th>\n",
       "      <th>qval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [feature, metadata, value, coef, stderr, N, N.not.0, pval, qval]\n",
       "Index: []"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MaAsLin2_results_phylum = pd.read_csv('MaAsLin2_results/phylum/significant_results.tsv', sep='\\t')\n",
    "MaAsLin2_results_phylum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it can be seen on each taxonomic level there is no significant results in differential microbial abundance. Anyway let's visualize these results to take a closer look!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 2.2. Visualization.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please open `RStudio` and go through  `Volcano_plots_journal.R` script.<br>\n",
    "There are a lot of manual adjustments to the plots to make it executable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Part 2.3. Alpha- and Beta diversities**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir Alpha_div"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alpha diversity calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# Usage\n",
    "# {path_to_script} {path_to_metadata} {path_to_counts} {path_to_output}\n",
    "Rscript.exe scripts/Alpha_div_calculations.R counts/csv/counts_species.csv Alpha_div/alpha_div_cult.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Part 2.3.1. Alpha- and Beta diversities visualization**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please open `RStudio` and go through  `Alpha_Beta_div_journal.R` script.<br>\n",
    "There are a lot of manual adjustments to the plots to make it executable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 3. Bar-plot - Mean Relative Abundance**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please open `RStudio` and go through  `Bar_plots_journal.R` script.<br>\n",
    "There are a lot of manual adjustments to the plots to make it executable."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "diff_an",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}