bb2022_sf-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for importing, formatting and data manipulation\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import glob\n",
    "import datetime\n",
    "#from time import time\n",
    "#from datetime import datetime\n",
    "#from datetime import timedelta\n",
    "import tempfile\n",
    "from qiime2 import Artifact\n",
    "import zipfile\n",
    "import yaml\n",
    "import json\n",
    "\n",
    "#for plotting\n",
    "import matplotlib, random\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib_venn import venn3, venn3_circles\n",
    "from matplotlib.patches import Patch\n",
    "import seaborn as sns\n",
    "#sns.set(style=\"whitegrid\")\n",
    "import plotly.express as px\n",
    "%matplotlib inline\n",
    "from IPython.display import display\n",
    "from upsetplot import plot\n",
    "#import pyupset as pyu\n",
    "pd.set_option('display.max_rows', 15)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "pd.set_option('display.width', 1000)\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import plotly.graph_objects as go\n",
    "from statsmodels.tsa.stattools import acf, pacf\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "register_matplotlib_converters()\n",
    "import scipy as sp\n",
    "import statsmodels.api as sm\n",
    "\n",
    "#for statistical analyses\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.covariance import EllipticEnvelope\n",
    "from sklearn.ensemble import IsolationForest\n",
    "from skbio.diversity import alpha_diversity\n",
    "from skbio.stats.distance import permanova\n",
    "from skbio import DistanceMatrix\n",
    "from scipy.spatial.distance import cdist\n",
    "from skbio.stats.composition import clr\n",
    "from skbio.stats.composition import alr\n",
    "from skbio.stats.composition import ilr\n",
    "from skbio.diversity.alpha import chao1\n",
    "from skbio.stats.composition import ancom\n",
    "import scipy.stats as stats\n",
    "import statsmodels.api as sa\n",
    "import statsmodels.formula.api as sfa\n",
    "import scikit_posthocs as spPH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import and format metadata from lab, and BBMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import lab metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_df():\n",
    "    \n",
    "    filenames = glob.glob('/Users/Diana/Documents/escuela/phd/ch2/bb_data/2022/METADATA.txt')\n",
    "    #load 2022 metadata and concatenate it into one dataframe\n",
    "    md = []\n",
    "    for filename in filenames:\n",
    "        df = pd.read_csv(filename, sep='\\t')\n",
    "        md.append(df)\n",
    "        print (filename)\n",
    "    \n",
    "    md = pd.concat(md)\n",
    "    \n",
    "    #drop empty columns and rows\n",
    "    md.dropna(how='all', axis=1, inplace=True) #empty cols\n",
    "    md.dropna(how='all', inplace=True) #empty rows\n",
    "    \n",
    "    return md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = load_df()\n",
    "#md = md[[\"sampleid\", \"[DNA]ng/ul\", \"A260/280\"]].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from platform import python_version\n",
    "\n",
    "print(python_version())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modified version of load_df() to accomodate different metadata files specific to the size fractioning samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_df():\n",
    "    \n",
    "    filenames = glob.glob('/Users/Diana/Documents/escuela/phd/ch2/bb_data/2022/*.csv')\n",
    "    #load all metadata from 2022 folder and concatenate them.\n",
    "    md = []\n",
    "    for filename in filenames:\n",
    "        df = pd.read_csv(filename)\n",
    "        df = df.rename(columns={'Depth Code 1-A, 5-B, 10-C, 60-D': 'depth_code',\n",
    "                            'Size Code 3um - L 0.2um - S': 'size_code'}) \n",
    "        md.append(df)\n",
    "        print (filename)\n",
    "    \n",
    "    md = pd.concat(md)\n",
    "    \n",
    "    #drop empty columns and rows\n",
    "    md.dropna(how='all', axis=1, inplace=True) #empty cols\n",
    "    md.dropna(how='all', inplace=True) #empty rows\n",
    "    \n",
    "    return md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initial version of code for whole BB dataset uses this function but for the size fractions the metadata was slightly differently formatted so we manipulated each \n",
    "#mdsf = load_df()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#upload the absorbance per sample data\n",
    "a260230 = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/ch2/bb_data/2022/a260230.csv\")\n",
    "a260230 = a260230.dropna(how='all') #drop null rows and columns\n",
    "a260230.dropna(how='all', axis=1, inplace=True)\n",
    "a260230 = a260230.replace({pd.NA: np.nan})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#upload metadata of non size fractionated samples\n",
    "noSF = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/ch2/bb_data/2022/noSF.csv\")\n",
    "noSF = noSF.dropna(how='all')\n",
    "noSF.dropna(how='all', axis=1, inplace=True)\n",
    "noSF = noSF.replace({pd.NA: np.nan})\n",
    "noSF = noSF.dropna(how='all') #drop null rows and columns\n",
    "#uncomment the line below to remove metadata columns\n",
    "#noSF = noSF[[\"sampleid\", \"[DNA]ng/ul\", \"A260/280\", \"date\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#upload metadata of size fractionated samples\n",
    "SF = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/ch2/bb_data/2022/SF.csv\")\n",
    "SF = SF.dropna(how='all')\n",
    "SF.dropna(how='all', axis=1, inplace=True)\n",
    "SF = SF.replace({pd.NA: np.nan})\n",
    "SF = SF.rename(columns={'Depth Code 1-A, 5-B, 10-C, 60-D': 'depth_code',\n",
    "                            'Size Code 3um - L 0.2um - S': 'size_code'}) \n",
    "#uncomment the line below to remove metadata columns\n",
    "#SF = SF[[\"sampleid\", \"[DNA]ng/ul\", \"A260/280\", \"date\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Renumber dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a dictionary for months\n",
    "month_dic = {\n",
    "    \"Jan\": 1,\n",
    "    \"Feb\": 2,\n",
    "    \"Mar\": 3,\n",
    "    \"Apr\": 4,\n",
    "    \"May\": 5,\n",
    "    \"Jun\": 6,\n",
    "    \"Jul\": 7,\n",
    "    \"Aug\": 8,\n",
    "    \"Sep\": 9,\n",
    "    \"Oct\": 10,\n",
    "    \"Nov\": 11,\n",
    "    \"Dec\": 12\n",
    "}\n",
    "month_season = {\n",
    "    \"Jan\": \"Winter\",\n",
    "    \"Feb\": \"Winter\",\n",
    "    \"Mar\": \"Spring\",\n",
    "    \"Apr\": \"Spring\",\n",
    "    \"May\": \"Spring\",\n",
    "    \"Jun\": \"Summer\",\n",
    "    \"Jul\": \"Summer\",\n",
    "    \"Aug\": \"Summer\",\n",
    "    \"Sep\": \"Autumn\",\n",
    "    \"Oct\": \"Autumn\",\n",
    "    \"Nov\": \"Autumn\",\n",
    "    \"Dec\": \"Winter\"\n",
    "}\n",
    "depth_num = {\n",
    "    \"A\": 1,\n",
    "    \"B\": 5,\n",
    "    \"C\": 10,\n",
    "    \"D\": 60,\n",
    "    \"E\": 30\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fill_dates(md):\n",
    "    if 'weekn' not in md:\n",
    "        md[\"weekn\"] = md[\"sampleid\"].str.extract(r'\\.([1-9][0-9]?)[A-E]')\n",
    "    md['weekn'] = pd.to_numeric(md['weekn'])\n",
    "    md['date'] = md.groupby(['sampleid','weekn'], sort=False)['date'].apply(lambda x: x.ffill().bfill())\n",
    "\n",
    "    #add month to a new column\n",
    "    md['month_name'] = md['date'].str.split('-').str[1]\n",
    "    md['year'] = 2022\n",
    "    md=md[md.year==2022]\n",
    "\n",
    "    #add month number\n",
    "    md['month']= md['month_name'].map(month_dic)\n",
    "\n",
    "    #add day number\n",
    "    md['day'] = md['date'].str.split('-').str[0]\n",
    "    md[[\"year\", \"month\", \"day\"]] = md[[\"year\", \"month\", \"day\"]].apply(pd.to_numeric)\n",
    "\n",
    "    #remove symbol for better handling of data\n",
    "    #md.rename(columns={\"Week#\": \"Weekn\"}, inplace=True)\n",
    "    #md.rename(columns={\"Depth\": \"depth\"}, inplace=True) #to match dfo\n",
    "\n",
    "    #change to int to remove decimals from date columns\n",
    "    md.year = md.year.apply(int)\n",
    "    md.day = md.day.apply(int)\n",
    "    md.month = md.month.apply(int)\n",
    "    #md.depth = md.depth.apply(int)\n",
    "    #md.weekn = md.weekn.apply(int)\n",
    "\n",
    "    #change to str to aggregate them into time_string to match dfos formatting of the date\n",
    "    md.year = md.year.apply(str)\n",
    "    md.month = md.month.apply(str)\n",
    "    md.day = md.day.apply(str)\n",
    "\n",
    "    md[\"depth_code\"] = md[\"sampleid\"].str.extract(r'[1-9][0-9]?([A-E])')\n",
    "    md['depth']= md['depth_code'].map(depth_num)\n",
    "    md['depth'] = pd.to_numeric(md['depth'])\n",
    "\n",
    "    #add leading zero to match date format in dfo metadata\n",
    "    md['month'] = md['month'].str.zfill(2)\n",
    "    md['day'] = md['day'].str.zfill(2)\n",
    "\n",
    "    md['time_string'] = md[['year', 'month', 'day']].agg('-'.join, axis=1)\n",
    "    \n",
    "    md[\"size_code\"] = md[\"sampleid\"].str.extract(r'[1-9][0-9]?[A-E]([L-S])')\n",
    "    md[\"size_code\"] = md[\"size_code\"].fillna('W')\n",
    "    \n",
    "    return md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SF = fill_dates(SF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noSF = noSF[noSF['weekn'] < 17]\n",
    "noSF = fill_dates(noSF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#verify which columns are shared between two dataframes\n",
    "a = np.intersect1d(SF.columns, noSF.columns)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colsnum = ['A260/280', '[DNA]ng/ul', 'depth', 'elution_volume', 'filtration_volume ', 'weekn', 'year']\n",
    "noSF[colsnum] = noSF[colsnum].apply(pd.to_numeric, errors='coerce', axis=1)\n",
    "SF[colsnum] = SF[colsnum].apply(pd.to_numeric, errors='coerce', axis=1)\n",
    "\n",
    "colstr = ['sampleid']\n",
    "SF[colstr] = SF[colstr].astype(\"string\")\n",
    "\n",
    "SF = SF.replace({pd.NA: np.nan})\n",
    "noSF = noSF.replace({pd.NA: np.nan})\n",
    "\n",
    "mdsf = noSF.merge(SF, on=['A260/280', 'Extracted_By', 'Notes', '[DNA]ng/ul', 'date', 'day',\n",
    "                           'depth', 'depth_code', 'elution_volume', 'extraction_date',\n",
    "                           'filtration_volume ', 'month', 'month_name', 'sampleid',\n",
    "                           'size_code', 'time_string', 'weekn', 'year'], how='outer')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.intersect1d(mdsf.columns, a260230.columns)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fill missing cell values with matching column name from other dataframe\n",
    "\n",
    "a260230[\"weekn\"] = a260230[\"sampleid\"].str.extract(r'\\.([1-9][0-9]?)[A-E]')\n",
    "a260230['weekn'] = pd.to_numeric(a260230['weekn'])\n",
    "    \n",
    "a260230 = a260230.fillna(noSF)\n",
    "mdsf = mdsf.fillna(a260230)\n",
    "\n",
    "a260230[\"depth_code\"] = a260230[\"sampleid\"].str.extract(r'[1-9][0-9]?([A-E])')\n",
    "a260230['depth']= a260230['depth_code'].map(depth_num)\n",
    "a260230['depth'] = pd.to_numeric(a260230['depth'])\n",
    "\n",
    "\n",
    "mdsf2 = mdsf.merge(a260230, on=['A260/280', '[DNA]ng/ul', 'extraction_date', 'sampleid', 'weekn', 'depth_code','depth'], how='outer')\n",
    "\n",
    "mdsf[\"weekn\"] = mdsf[\"sampleid\"].str.extract(r'\\.([1-9][0-9]?)[A-E]')\n",
    "mdsf['weekn'] = pd.to_numeric(mdsf['weekn'])\n",
    "mdsf['date'] = mdsf.groupby(['weekn'], sort=False)['date'].apply(lambda x: x.ffill().bfill())\n",
    "\n",
    "mdsf = mdsf[mdsf['weekn'] < 17]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = mdsf2.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md.sort_values(by=['weekn', 'depth'],inplace=True)\n",
    "md = md.fillna(method='ffill')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md.to_csv('regardons.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import and manage BBMP data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Metadata __md__ is formatted. It contains 38 columns.\n",
    "__md__ is the lab's metadata for sampling, extraction and sequencing. \\\n",
    "__dfo_md__ is BBMP remote sensing data (salinity, pH, temperature, density..) \\\n",
    "__bio_niskin__ is nutrient data \\\n",
    "Format __bio_niskin__ data to merge with __md__. __bio_niskin__ is 32 columns, including year, month, day, and depth. __dfo_md__ also has 32 columns, including year_time, month_time, day_time. To merge these data with __md__, we will change the time stamps columns to the same name, and generate a time_string column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfo_md = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/ch2/bb_data/bbmp_aggregated_profiles.csv\")\n",
    "bio_niskin = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/ch2/bb_data/BBMP_Data_2022.csv\")#\n",
    "#dfo_metadata_y14 = pd.read_csv(\"/Users/Diana/Documents/escuela/phd/bb_data/2019/data_export/trim-analysis/dfo_metadata_y14.tsv\", sep='\\t')\n",
    "\n",
    "#change to str to aggregate them into time_string\n",
    "bio_niskin = bio_niskin[bio_niskin.year==2022]\n",
    "bio_niskin.year = bio_niskin.year.apply(str)\n",
    "bio_niskin.month = bio_niskin.month.apply(str)\n",
    "bio_niskin.day = bio_niskin.day.apply(str)\n",
    "#add leading zero to match date format in dfo metadata\n",
    "bio_niskin['month'] = bio_niskin['month'].str.zfill(2)\n",
    "bio_niskin['day'] = bio_niskin['day'].str.zfill(2)\n",
    "\n",
    "bio_niskin['time_string'] = bio_niskin[['year', 'month', 'day']].agg('-'.join, axis=1)\n",
    "\n",
    "#make a new column for time_string without the time\n",
    "dfo_md=dfo_md[dfo_md.year_time==2022]\n",
    "dfo_md['time_string_time'] = dfo_md['time_string']\n",
    "dfo_md['time_string'] = dfo_md['time_string'].str.split(' ').str[0]\n",
    "\n",
    "#renaming columns to ensure correct merging\n",
    "dfo_md.rename(columns={\"depth\":\"bbmpdepth\",\"pressure\": \"depth\", \"year_time\": \"year\", \"month_time\": \"month\", \"day_time\": \"day\"}, inplace=True)\n",
    "\n",
    "#change to int to remove decimals from date columns\n",
    "cols = ['year', 'depth', 'month', 'day']\n",
    "md[cols] = md[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "dfo_md[cols] = dfo_md[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "bio_niskin[cols] = bio_niskin[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "\n",
    "#drop empty columns and rows\n",
    "dfo_md.dropna(how='all', axis=1, inplace=True) #empty cols\n",
    "dfo_md.dropna(how='all', inplace=True) #empty rows\n",
    "\n",
    "bio_niskin.dropna(how='all', axis=1, inplace=True) #empty cols\n",
    "bio_niskin.dropna(how='all', inplace=True) #empty rows\n",
    "\n",
    "#make a season column\n",
    "md['season'] = ''\n",
    "\n",
    "for month, season in month_season.items():\n",
    "    md.loc[md['month_name'] == month, 'season'] = season"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bio_niskin data has exact recorded depths, whereas BB sample data is restricted to categories: make a new column to allow for data merging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depths = np.array([1,5,10,60])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bio_niskin2= bio_niskin.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length = 48 #number of weeks for the tile repeat\n",
    "bio_niskin2['NewDepth'] = pd.DataFrame({'NewDepth': np.tile(depths, length)}) #tile depth categorical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bio_niskin2=bio_niskin2.assign(NewDepth=depths[np.arange(len(bio_niskin2)) % len(depths)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#order the two depth columns at the end of the dataframe to visually examine\n",
    "cols_at_end = ['depth', 'NewDepth']\n",
    "bio_niskin3 = bio_niskin2[[c for c in bio_niskin2 if c not in cols_at_end] \n",
    "        + [c for c in cols_at_end if c in bio_niskin2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rename columns to ensure correct merging\n",
    "bio_niskin3.rename(columns={'depth': 'truedepth', 'NewDepth': 'depth'}, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make merging columns to same type\n",
    "bio_niskin3[cols] = bio_niskin3[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "md[cols] = md[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "dfo_md[cols] = md[cols].apply(pd.to_numeric, errors='ignore', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.intersect1d(md.columns, dfo_md.columns)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#preview column types to allow for merging\n",
    "#pd.set_option('display.max_rows', 35)\n",
    "#md.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert merging columns to same type\n",
    "colsnum = ['day', 'month']\n",
    "dfo_md[colsnum] = dfo_md[colsnum].apply(pd.to_numeric, errors='coerce', axis=1)\n",
    "md[colsnum] = md[colsnum].apply(pd.to_numeric, errors='coerce', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#merging party\n",
    "merged = pd.merge(md, dfo_md, on=['day', 'depth', 'month', 'pH', 'sigmaTheta', 'theta',\n",
    "                                  'time_string', 'year'], how=\"left\")\n",
    "allyears = pd.merge(md, dfo_md, on=['day', 'depth', 'month', 'pH', 'sigmaTheta', 'theta',\n",
    "                                    'time_string', 'year'], how=\"outer\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = merged[merged['weekn'] < 17]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged.to_csv('regardons2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert column type to numeric for merging\n",
    "allyears[cols] = allyears[cols].apply(pd.to_numeric, errors='ignore', axis=1)\n",
    "\n",
    "#merged = merged.drop(index=237) #delete a row with missing information\n",
    "merged[cols] = merged[cols].apply(pd.to_numeric, axis=1)\n",
    "bio_niskin3[cols] = bio_niskin3[cols].apply(pd.to_numeric, axis=1)\n",
    "\n",
    "#add nutrient data\n",
    "#uncomment the line below if  need access to metadata outside the 16weeks samples in 2022\n",
    "#preall_md= pd.merge(allyears, bio_niskin3, on=[\"day\", \"month\", \"year\", 'depth'], how=\"outer\")\n",
    "all_md = pd.merge(merged, bio_niskin3, on=[\"day\", \"month\", \"year\", 'depth'], how=\"left\")\n",
    "\n",
    "#split dfs by depth\n",
    "shallow_depths = [1, 5, 10]\n",
    "shallow = all_md[all_md[\"depth\"] < 30]\n",
    "#shallow = shallow.groupby(['year', 'month', \"day\"]).mean().reset_index()\n",
    "deep = all_md[all_md.depth == 60]\n",
    "\n",
    "#split dfs by season\n",
    "year_season = all_md.groupby(by = ['year','season']).mean().reset_index()\n",
    "\n",
    "Winter = year_season.loc[year_season['season'] == 'Winter',:]\n",
    "Spring = year_season.loc[year_season['season'] == 'Spring',:]\n",
    "Summer = year_season.loc[year_season['season'] == 'Summer',:]\n",
    "Autumn = year_season.loc[year_season['season'] == 'Autumn',:]\n",
    "\n",
    "#save output as csv\n",
    "all_md.to_csv('allmetadata.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d1 = all_md[all_md.depth_code == 'A']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(data=d1, x=\"weekn\", y=\"Chlorophyll A\", color=\"0.8\", linewidth=.75, kind=\"line\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.relplot(data=d1, x=\"weekn\", y=\"Phosphate\", color=\"0.8\", linewidth=.75, kind=\"line\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find rows with null values at given column\n",
    "emptynit = merged[merged['depth'].isna()]\n",
    "emptynit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emptynit = merged[merged['temperature'].isna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emptynit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting party"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plotly seasonal averages figure\n",
    "fig2 = go.Figure()\n",
    "for template in [\"plotly_white\"]:\n",
    "    fig2.add_trace(go.Scatter(x=Winter['year'], y=Winter['temperature'],\n",
    "                    mode='lines',\n",
    "                    name='Winter',\n",
    "                    marker_color='#838B8B'))\n",
    "    fig2.add_trace(go.Scatter(x=Spring['year'], y=Spring['temperature'],\n",
    "                    mode='lines',\n",
    "                    name='Spring',\n",
    "                    marker_color='#FFB5C5'))\n",
    "    fig2.add_trace(go.Scatter(x=Summer['year'], y=Summer['temperature'],\n",
    "                    mode='lines',\n",
    "                    name='Summer',\n",
    "                    marker_color='#87CEFF'))\n",
    "    fig2.add_trace(go.Scatter(x=Autumn['year'], y=Autumn['temperature'],\n",
    "                    mode='lines',\n",
    "                    name='Autumn',\n",
    "                    marker_color='#FF8000'))\n",
    "    fig2.update_layout(\n",
    "    height=800,\n",
    "    xaxis_title=\"Years\",\n",
    "    yaxis_title='Temperature in degree',\n",
    "    title_text='Average Temperature seasonwise over the years',\n",
    "    template=template)\n",
    "\n",
    "fig2.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#seaborn season averages plot\n",
    "sns.lineplot(year_season['year'],year_season['temperature'], hue =year_season[\"season\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detect and plot anomalies in variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def detect_anomalies(metadata, df, dpt, yr=all, month=all):\n",
    "    \n",
    "    sfd=df[df.depth==dpt]\n",
    "    \n",
    "    md_col = sfd[['event_id', metadata, \"year\", \"month\"]].copy()\n",
    "    md_col = md_col[md_col[metadata].notna()]\n",
    "    if yr != all:\n",
    "        #mdcol_yr = md_col[md_col.Year == yr]\n",
    "        mdcol_yr = md_col[md_col['year'].isin(yr)]\n",
    "    else: \n",
    "        mdcol_yr = md_col\n",
    "        \n",
    "    if month != all:\n",
    "        #mdcol_yr = mdcol_yr[mdcol_yr.Month == month]\n",
    "        mdcol_yr = mdcol_yr[mdcol_yr['month'].isin(month)]\n",
    "    \n",
    "    mdcol_yr = mdcol_yr.drop(columns=['year', \"month\"])\n",
    "    mdcol_yr = mdcol_yr.set_index(['event_id'])\n",
    "    \n",
    "    #modelling time\n",
    "    outliers_fraction = float(.01)\n",
    "    scaler = StandardScaler()\n",
    "    np_scaled = scaler.fit_transform(mdcol_yr.values.reshape(-1, 1))\n",
    "    data = pd.DataFrame(np_scaled)\n",
    "    # train isolation forest\n",
    "    model =  IsolationForest(contamination=outliers_fraction)\n",
    "    model.fit(data)\n",
    "    \n",
    "    #predict data\n",
    "    mdcol_yr['anomaly'] = model.predict(data)\n",
    "    \n",
    "    # visualization\n",
    "    fig, ax = plt.subplots(figsize=(10,6))\n",
    "    a = mdcol_yr.loc[mdcol_yr['anomaly'] == -1, [metadata]] #anomaly\n",
    "    ax.plot(mdcol_yr.index, mdcol_yr[metadata], color='black', label = 'Normal')\n",
    "    ax.scatter(a.index,a[metadata], color='red', label = 'Anomaly')\n",
    "    #plt.axvline(36, ls='--')\n",
    "    plt.legend()\n",
    "    plt.show();\n",
    "    #add axes names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detect_anomalies('Temperature', bio_niskin3, 1, yr={2022}, month={1,2,3,4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detect_anomalies('Temperature', bio_niskin3, 5, yr={2022}, month={1,2,3,4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detect_anomalies('Temperature', bio_niskin3, 10, yr={2022}, month={1,2,3,4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detect_anomalies('Temperature', bio_niskin3, 60, yr={2022}, month={1,2,3,4})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add prokaryotic community"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Upload functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Special thanks to Alex Manuele https://github.com/alexmanuele\n",
    "def consolidate_tables(MG):\n",
    "    if MG == '16S':\n",
    "        comm = '02-PROKs'\n",
    "    else :\n",
    "        comm = '02-EUKs'\n",
    "        \n",
    "    table_list = glob.glob('{0}/table.qza'.format('/Users/Diana/Documents/escuela/phd/size_fractions/BB22_size-fraction-comparison-analysed/to_transfer/'+comm))\n",
    "    print(\"Found all \"+MG+\" tables.\")\n",
    "\n",
    "        \n",
    "    dataframes = []  \n",
    "    for table_path in table_list:\n",
    "        with tempfile.TemporaryDirectory() as tempdir:\n",
    "            #load table, dump contents to tempdir\n",
    "            table = Artifact.load(table_path)\n",
    "            #Make sure the tables are all FeatureFrequency type\n",
    "            assert str(table.type) == 'FeatureTable[Frequency]', \"{0}: Expected FeatureTable[Frequency], got {1}\".format(table_path, table.type)\n",
    "            Artifact.extract(table_path, tempdir)\n",
    "            #get the provenance form the tempdir and format it for DF\n",
    "            prov = '{0}/{1}/provenance/'.format(tempdir, table.uuid)\n",
    "            action = yaml.load(open(\"{0}action/action.yaml\".format(prov), 'r'), Loader=yaml.BaseLoader)\n",
    "            paramlist = action['action']['parameters']\n",
    "            paramlist.append({'table_uuid': \"{}\".format(table.uuid)})\n",
    "            paramdict = {}\n",
    "            for record in paramlist:\n",
    "                paramdict.update(record)\n",
    "\n",
    "            # Get the data into a dataframe\n",
    "              #Biom data\n",
    "            df = table.view(pd.DataFrame).unstack().reset_index()\n",
    "            df.columns = ['feature_id', 'sample_name', 'feature_frequency']\n",
    "            df['table_uuid'] = [\"{}\".format(table.uuid)] * df.shape[0]\n",
    "              #param data\n",
    "            pdf = pd.DataFrame.from_records([paramdict])\n",
    "              #merge params into main df\n",
    "            df = df.merge(pdf, on='table_uuid')\n",
    "            \n",
    "\n",
    "            #I like having these columns as the last three. Makes it more readable\n",
    "            cols = df.columns.tolist()\n",
    "            reorder = ['sample_name', 'feature_id', 'feature_frequency']\n",
    "            for val in reorder:\n",
    "                cols.append(cols.pop(cols.index(val)))\n",
    "            df = df[cols]\n",
    "            df['table_path'] = [table_path] * df.shape[0]\n",
    "            df['sample_name'] = df['sample_name'].str.replace('-', '.')\n",
    "            dataframes.append(df)\n",
    "            \n",
    "            # Adding table_id, forward and reverse trim columns\n",
    "            #df['table_id'] = str(table_path.split('/')[-3]) #add a table_id column\n",
    "            #df['forward_trim'], df['reverse_trim'] = df['table_id'].str.split('R', 1).str\n",
    "            #df['forward_trim'] = df['forward_trim'].map(lambda x: x.lstrip('F'))\n",
    "            #df[\"forward_trim\"] = pd.to_numeric(df[\"forward_trim\"])\n",
    "            #df[\"reverse_trim\"] = pd.to_numeric(df[\"reverse_trim\"])\n",
    "\n",
    "    #Stick all the dataframes together\n",
    "    #outputfile=\"merged_all_tables.tsv\"\n",
    "    df = pd.concat(dataframes)\n",
    "    df['sample_name'] = df['sample_name'].str.replace(r'\\.S([1-9]|[1-9][0-9]|[1-9][0-9][0-9]).L001\\.','', regex=True)\n",
    "    \n",
    "    #df.to_csv(comm+'/merged_all_tables.tsv', sep='\\t', index=False)\n",
    "    print(\"Successfully saved all tables.\")\n",
    "    return df, MG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_metadata(df):\n",
    "    #df = pd.read_csv('02-PROKs/'+'/merged_all_tables.tsv', sep='\\t')\n",
    "\n",
    "    tables = df[['sample_name', 'feature_id', 'feature_frequency']].copy()\n",
    "    tables.rename(columns={'sample_name':'sampleid'}, inplace=True)\n",
    "\n",
    "    all_md['sampleid'] = all_md['sampleid'].str.replace('_', '.')\n",
    "    merged = pd.merge(tables,all_md, on='sampleid', how='left') #all_md is the metadata file\n",
    "    merged = merged[merged.feature_frequency != 0]\n",
    "    \n",
    "    merged['year'] = 2022\n",
    "\n",
    "    merged[\"size_code\"] = merged[\"sampleid\"].str.extract(r'[1-9][0-9]?[A-E]([L-S])')\n",
    "    merged[\"size_code\"] = merged[\"size_code\"].fillna('W')\n",
    "    merged[\"depth_code\"] = merged[\"sampleid\"].str.extract(r'[1-9][0-9]?([A-E])')\n",
    "    merged['depth']= merged['depth_code'].map(depth_num)\n",
    "    merged[\"weekn\"] = merged[\"sampleid\"].str.extract(r'\\.([1-9][0-9]?)[A-E]')\n",
    "    merged['weekn'] = pd.to_numeric(merged['weekn'])\n",
    "    merged['depth'] = pd.to_numeric(merged['depth'])\n",
    "    merged['date'] = merged.groupby('weekn', as_index=False)['date'].transform('first')\n",
    "    \n",
    "    merged['Total'] = merged['feature_frequency'].groupby(merged['sampleid']).transform('sum')\n",
    "    merged['ratio'] = merged['feature_frequency']/merged['Total']\n",
    "    merged['nASVs'] = merged['feature_id'].groupby(merged['sampleid']).transform('count')\n",
    "    merged['weekdepth'] = merged[\"weekn\"].astype(str) + merged[\"depth\"].astype(str)\n",
    "    merged['avg'] = merged['nASVs'].groupby(merged['weekdepth']).transform('mean')\n",
    "    merged['diff'] = merged['nASVs'] - merged['avg']\n",
    "\n",
    "    print('Set up metadata ...')\n",
    "    \n",
    "    #merged.to_csv(comm+'/merged_asvs_metadata.tsv', sep = '\\t')\n",
    "    print('Saved merged_asvs_metadata.tsv')\n",
    "    \n",
    "    return merged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pick_metadata(merged, depth='all', size_fraction='both', year='all', R='all', F='all', txsubset = 'all'):\n",
    "#make df of features/composition+run+comm\n",
    "\n",
    "    depth = depth\n",
    "    year = year\n",
    "    size_fraction = size_fraction\n",
    "    txsubset = txsubset\n",
    "        \n",
    "    files = glob.glob('{0}/*/class/*/data/taxonomy.tsv'.format('/Users/Diana/Documents/escuela/phd/size_fractions/BB22_size-fraction-comparison-analysed/to_transfer'))\n",
    "    taxos = []\n",
    "#    if not os.path.exists(path+composition):\n",
    "#        os.mkdir(path+composition)\n",
    "    for filename in files:\n",
    "        tax = pd.read_csv(filename, sep='\\t')\n",
    "        taxos.append(tax)\n",
    "        \n",
    "    print('Appended all taxonomies to taxos')\n",
    "    taxos = pd.concat(taxos)\n",
    "    taxos = taxos.rename(columns={\"Feature ID\": \"feature_id\"}, errors=\"raise\")\n",
    "    taxos = taxos.drop_duplicates()\n",
    "\n",
    "    separated = merged.merge(taxos, how='left', on='feature_id') #merged excludes features of frequency = 0\n",
    "    separated = separated.drop_duplicates()\n",
    "    \n",
    "    if depth != 'all':\n",
    "        separated = separated[separated[\"depth\"] == depth]\n",
    "    if size_fraction != 'both':\n",
    "        separated = separated[separated[\"size_fraction\"] == size_fraction]\n",
    "\n",
    "    separated[['Domain', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species']] = separated['Taxon'].str.split('; ', expand=True)\n",
    "    cols = ['Domain', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species']\n",
    "    for col in cols:\n",
    "        separated[col] = separated[col].fillna('Unassigned')\n",
    "        \n",
    "    separated['Month'] = separated['date'].str.split('-').str[1]\n",
    "    \n",
    "    #separated['total'] = separated.groupby(['table_id','sample-id'])['feature_frequency'].transform('sum')\n",
    "    #separated['ratio'] = separated['feature_frequency']/(separated['total'])\n",
    "    #separated_taxonomies = separated.copy()\n",
    "    \n",
    "    #make a dictionary with keys for id-ing the taxon belonging to this sub-community\n",
    "    #separated_dic = pd.Series(separated.Taxon.values,separated.feature_id.values).to_dict()\n",
    "    print('Saved separated by metadata dataframe.')\n",
    "    \n",
    "    return separated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def taxbarplot(separated, level, depth, topn): #separated is the df, #level is a string of taxonomic level column name, depth is an integer\n",
    "    sfd=separated[separated.depth==depth]\n",
    "    toptaxa = sfd[['feature_id', 'feature_frequency', 'Taxon', 'size_code', 'depth','weekn', level]].copy()\n",
    "    toptaxa = toptaxa.drop_duplicates()\n",
    "    df_agg = toptaxa.groupby(['size_code',level, 'depth']).agg({'feature_frequency':sum})\n",
    "    topd = df_agg['feature_frequency'].groupby('size_code', group_keys=False).nlargest(topn)\n",
    "    topd = topd.to_frame()\n",
    "    topd = topd.reset_index()\n",
    "\n",
    "\n",
    "    df_agg = df_agg.reset_index()\n",
    "    df_agg['set_name'] = df_agg['size_code']+df_agg['depth'].astype(str)\n",
    "    \n",
    "    cumulab = separated[['feature_frequency', 'depth', 'size_code', 'Genus']].copy()\n",
    "    cumulab1 = cumulab.groupby(['Genus']).agg({'feature_frequency':sum})\n",
    "\n",
    "    resultpivot = df_agg.pivot_table(index=level, columns='set_name', values='feature_frequency')\n",
    "    resultpivot = resultpivot.fillna(0)\n",
    "    resultpivot[resultpivot != 0] = 1\n",
    "    tosave = pd.merge(resultpivot, cumulab1, left_index=True, right_index=True)\n",
    "    tosave.to_csv(level+'_'+str(depth)+'16S_relab.csv')\n",
    "    \n",
    "    top10d_list = topd[level].unique()\n",
    "    top10d = sfd.copy()\n",
    "    top10d.loc[~top10d[level].isin(top10d_list), level] = 'Other' #isnot in top list\n",
    "    phyld = top10d.groupby(['size_code','weekn', level])['ratio'].sum()\n",
    "    phyld = phyld.reset_index()\n",
    "\n",
    "\n",
    "    fig = px.bar(phyld, x=\"size_code\", y=\"ratio\", facet_col=\"weekn\", color=level, labels={\n",
    "                     \"feature_frequency\": \"Relative abundance\",\n",
    "                     \"size_code\": \"\",\n",
    "                     \"weekn\": \"w\"}, color_discrete_map=palette_dict)\n",
    "    fig.update_xaxes(type='category', dtick=1)\n",
    "    fig.update_layout(\n",
    "        title=\"Relative abundance of top 10\" + level + 'observed at Depth' + str(depth),\n",
    "        yaxis_title=\"Relative abundance\",\n",
    "        xaxis_title=\"Size fraction\",\n",
    "        legend_title=level,\n",
    "        font=dict(size=8)\n",
    "    )\n",
    "\n",
    "    fig.show()\n",
    "    #fig.write_image(\"outputs/fig1.png\")\n",
    "    #fig.to_image(format=\"png\")\n",
    "    \n",
    "    return phyld, top10d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcaplot(separated, depth, comm, columnperm, spc):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        folder = '02-PROKs'\n",
    "    else:\n",
    "        folder = '02-EUKs'\n",
    "        \n",
    "    \n",
    "    if depth == 'all':\n",
    "        df = separated.copy()\n",
    "    else:\n",
    "        df=separated[separated.depth==depth]\n",
    "        \n",
    "    \n",
    "    if 'SL' in separated['size_code'].unique():\n",
    "        #sizecode palette codes\n",
    "        sizecodes = ['S', 'L', 'W', 'SL']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "        dicsc = pd.Series(df.size_code.values,index=df.sampleid).to_dict()\n",
    "        color_rows_sc = {k: palette_dict[v] for k, v in dicsc.items()}\n",
    "        seriescr = pd.Series(color_rows_sc)\n",
    "    \n",
    "    else:\n",
    "        #sizecode palette codes\n",
    "        sizecodes = ['S', 'L', 'W']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "        dicsc = pd.Series(df.size_code.values,index=df.sampleid).to_dict()\n",
    "        color_rows_sc = {k: palette_dict[v] for k, v in dicsc.items()}\n",
    "        seriescr = pd.Series(color_rows_sc)\n",
    "    \n",
    "    #month palette code\n",
    "    df['Month'] = df['date'].str.split('-').str[1]\n",
    "    months = ['Jan', 'Feb', 'Mar', 'May', 'Apr']\n",
    "    palette_colors = sns.color_palette(\"flare\")\n",
    "    palette_dict_month = {monthname: color for monthname, color in zip(months, palette_colors)}\n",
    "    dic = pd.Series(df.Month.values,index=df.sampleid).to_dict()\n",
    "    color_rows_month = {k: palette_dict_month[v] for k, v in dic.items()}\n",
    "    seriesmonthcr = pd.Series(color_rows_month)\n",
    "\n",
    "    dfcolors = pd.DataFrame({'Month': seriesmonthcr,'Size code':seriescr})\n",
    "    \n",
    "    topiv = df[['feature_id', 'feature_frequency', 'sampleid']].copy()\n",
    "    topiv = topiv.drop_duplicates()\n",
    "    \n",
    "    sfdpiv= topiv.pivot(index='sampleid', columns='feature_id', values='feature_frequency')\n",
    "    sfdpiv=sfdpiv.fillna(0)\n",
    "    sfdclr=sfdpiv.mask(sfdpiv==0).fillna(0.1)\n",
    "    clr_transformed_array = clr(sfdclr)\n",
    "    samples = sfdpiv.index\n",
    "    asvs = sfdpiv.columns\n",
    "    \n",
    "    #Creating the dataframe with the clr transformed data, and assigning the sample names\n",
    "    clr_transformed = pd.DataFrame(clr_transformed_array, columns=asvs)\n",
    "    #Assigning the asv names\n",
    "    clr_transformed['samples'] = samples\n",
    "    clr_transformed = clr_transformed.set_index('samples')\n",
    "    clr_transformed.head()\n",
    "\n",
    "    #calculate distance matrix\n",
    "    dist = cdist(clr_transformed, clr_transformed, 'euclid')\n",
    "    distance_matrix = pd.DataFrame(dist, columns=samples)\n",
    "    distance_matrix['samples'] = samples\n",
    "    distance_matrix = distance_matrix.set_index('samples')\n",
    "\n",
    "    #format for pca\n",
    "    dm = DistanceMatrix(distance_matrix)\n",
    "\n",
    "    pca = PCA(n_components=2)\n",
    "    pca_features = pca.fit_transform(distance_matrix)\n",
    "    \n",
    "    ####\n",
    "    sns.set(rc={\"figure.figsize\":(4, 3)})\n",
    "    sns.set_style(\"whitegrid\", {'axes.grid' : False})\n",
    "    plot_df = pd.DataFrame(data = pca_features, columns = ['dim1', 'dim2'], index = sfdpiv.index)\n",
    "    plot_df['dim1'] = plot_df['dim1']/1000\n",
    "    plot_df['dim2'] = plot_df['dim2']/1000\n",
    "    if depth =='all':\n",
    "        plot_df2 = pd.merge(plot_df,df[['sampleid','size_code','depth']],on='sampleid', how='left')\n",
    "    else:\n",
    "        plot_df2 = pd.merge(plot_df,df[['sampleid','size_code','weekn']],on='sampleid', how='left')\n",
    "        \n",
    "    \n",
    "    ##divide into pre-post bloom\n",
    "    def get_stage(weekNb):\n",
    "        if weekNb < 10:\n",
    "            return 'Pre-bloom'\n",
    "        elif weekNb == 10 :\n",
    "            return 'Bloom'\n",
    "        elif weekNb > 10:\n",
    "            return 'Bloom'\n",
    "    \n",
    "    if depth != 'all':\n",
    "        plot_df2['Time'] = plot_df2['weekn'].apply(get_stage)\n",
    "    \n",
    "    plot_df2 = plot_df2.rename(columns={'size_code': 'Size code'})\n",
    "    \n",
    "    pc1v = round(pca.explained_variance_ratio_[0]*100)\n",
    "    pc2v = round(pca.explained_variance_ratio_[1]*100)\n",
    "    \n",
    "    #plot_df2 = plot_df2.drop_duplicates()\n",
    "    #dfperm = plot_df2.set_index('sampleid')\n",
    "    \n",
    "    #permanova2 = permanova(dm, dfperm, columnperm)\n",
    "    #results = permanova2(999)\n",
    "    \n",
    "    #plot\n",
    "    \n",
    "    if depth == 'all':\n",
    "        var2 = 'depth'\n",
    "    else:\n",
    "        var2 = 'Time'\n",
    "    \n",
    "    sns.set_style(\"white\")\n",
    "    ax=sns.scatterplot(x = 'dim1', y = 'dim2', hue= 'Size code', style=var2, data = plot_df2, \n",
    "                       palette=palette_dict)#, size = 'Week_Group')#,palette=sns.color_palette(\"dark:salmon_r\", as_cmap=True))\n",
    "    plt.ylabel('PCo2 (' + str(pc2v) + '% variance explained)')\n",
    "    plt.xlabel('PCo1 (' + str(pc1v) +'% variance explained)')\n",
    "    ax.set_title('Depth ' + str(depth) + 'm', loc='left', weight='bold')\n",
    "    plt.legend(frameon=False)\n",
    "    sns.move_legend(ax, \"upper left\", bbox_to_anchor=(1, 1))\n",
    "    sns.despine()\n",
    "    plt.savefig('outputs/'+folder+'/D'+str(depth)+spc+'_PCAplot.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "    \n",
    "    print ( \"Components = \", pca.n_components_ , \";\\nTotal explained variance = \",\n",
    "      round(pca.explained_variance_ratio_.sum(),5)  )\n",
    "    \n",
    "    print (\"Components 1 and 2 are\", pca.explained_variance_ratio_)\n",
    "    \n",
    "    # Retrieve Loadings\n",
    "    loadings = pca.components_\n",
    "\n",
    "    # Summarize Loadings by Metadata Category\n",
    "    metadata_groups = plot_df2[var2].unique()\n",
    "    metadata_contributions = {}\n",
    "    \n",
    "    for group in metadata_groups:\n",
    "        group_variables = plot_df2.loc[plot_df2[var2] == group, 'sampleid']\n",
    "        group_loadings = np.abs(loadings[:, [list(distance_matrix.columns).index(var) for var in group_variables]]).mean(axis=1)\n",
    "        metadata_contributions[group] = group_loadings\n",
    "\n",
    "    # Visual Representation\n",
    "    for group, contributions in metadata_contributions.items():\n",
    "        plt.barh(contributions, group) #range(1, len(contributions) + 1),\n",
    "\n",
    "    plt.ylabel('Principal Component')\n",
    "    plt.xlabel('Average Loading Contribution')\n",
    "    sns.despine()\n",
    "    plt.legend(frameon=False)\n",
    "    plt.savefig('outputs/'+folder+'/D'+str(depth)+spc+'_PCAplot_brplot.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "        \n",
    "\n",
    "    ##clustermap\n",
    "    ax = sns.clustermap(distance_matrix, method=\"complete\", cmap='RdBu', annot=True,\n",
    "               yticklabels=True, row_colors = dfcolors,\n",
    "               annot_kws={\"size\": 7}, figsize=(15,12));\n",
    "\n",
    "    handles1 = [Patch(facecolor=palette_dict_month[key]) for key in palette_dict_month]\n",
    "    plt.legend(handles1, palette_dict_month, title='Month',\n",
    "               bbox_to_anchor=(1, 1), bbox_transform=plt.gcf().transFigure, loc='upper left')\n",
    "    \n",
    "    plt.savefig('outputs/'+folder+'/D'+str(depth)+spc+'_clustermap.png', dpi=200, bbox_inches=\"tight\")\n",
    "\n",
    "\n",
    "    return pca, pca_features, sfdclr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boxplot_depth(separated, comm, depth, ycolumn, yaxislabel='def'):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "    \n",
    "    if yaxislabel != 'def':\n",
    "        ycol = ycolumn\n",
    "    \n",
    "    #sfd=separated[separated.depth==depth]\n",
    "    sfd = separated.copy()\n",
    "    \n",
    "    #sfd_S = sfd[['size_code', 'nASVs', 'weekn']].copy()\n",
    "    #sfd_S = sfd_S.drop_duplicates()\n",
    "    #sdfpv = sfd_S.pivot(index='weekn', columns='size_code', values='nASVs')\n",
    "    #fvalue, pvalue = stats.f_oneway(sdfpv['L'], sdfpv['S'], sdfpv['W'])\n",
    "    \n",
    "    sfd_LM = sfd[['size_code', 'nASVs']].copy()\n",
    "    sfd_LM = sfd_LM.drop_duplicates()\n",
    "    lm = sfa.ols('nASVs ~ C(size_code)', data=sfd_LM).fit()\n",
    "    anova = sa.stats.anova_lm(lm)\n",
    "    results = spPH.posthoc_ttest(sfd_LM, val_col='nASVs', group_col='size_code', p_adjust='holm')\n",
    "    \n",
    "    if 'SL' in separated['size_code'].unique():\n",
    "        #sizecode palette codes\n",
    "        sizecodes = ['S', 'L', 'W', 'SL']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "    \n",
    "    else:\n",
    "        #define color palettes\n",
    "        sizecodes = ['S', 'L', 'W']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "    \n",
    "    \n",
    "    \n",
    "    #plot\n",
    "    sns.set(rc={\"figure.figsize\":(4, 3)})\n",
    "    sns.set_style(\"ticks\")\n",
    "    sns.boxplot(data=sfd, x=\"size_code\", y=ycolumn, palette=palette_dict, order=sizecodes)#, hue=\"size_code\")\n",
    "    sns.despine()\n",
    "    plt.ylabel(yaxislabel, fontsize=20)\n",
    "    plt.xlabel('Size fraction', fontsize=20)\n",
    "\n",
    "    #g.tick_params(labelsize=15)\n",
    "    plt.savefig('outputs/'+comm_id+'/D'+str(depth)+'_adboxplot.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.clf()\n",
    "    plt.cla()\n",
    "    plt.close()\n",
    "    \n",
    "    \n",
    "    sns.set(rc={\"figure.figsize\":(7, 3)})\n",
    "    sns.set_style(\"ticks\")\n",
    "    ax=sns.barplot(data=sfd, x=\"weekn\", y=\"diff\", hue=\"size_code\", palette=palette_dict,\n",
    "                  capsize=.15, errwidth=0.5)#, hue=\"size_code\")\n",
    "    sns.despine()\n",
    "    plt.ylabel('Number of ASVs relative to weekly average', fontsize=20)\n",
    "    plt.xlabel('Week number', fontsize=20)\n",
    "    sns.move_legend(ax, \"upper left\", bbox_to_anchor=(1, 1))\n",
    "    plt.savefig('outputs/'+comm_id+'/D'+str(depth)+'_avgbarplot.png', dpi=200, bbox_inches=\"tight\")\n",
    "\n",
    "    plt.clf() \n",
    "    \n",
    "\n",
    "#     glue = sfd[['size_code', 'weekn', 'diff']].copy()\n",
    "#     glue = glue.drop_duplicates()\n",
    "#     glue = glue.pivot(index=\"size_code\", columns=\"weekn\", values=\"diff\")\n",
    "#     floored_data = glue.apply(np.floor)\n",
    "#     sns.set_style('ticks')\n",
    "#     plt.figure(figsize=(8, 2))\n",
    "#     cmap = sns.diverging_palette(240,240, as_cmap=True)\n",
    "#     ax = sns.heatmap(floored_data, yticklabels=True, linewidths=.5, annot=True, annot_kws={\"fontsize\":8},\n",
    "#                     cmap = cmap)\n",
    "#     plt.savefig('outputs/'+comm_id+'/heatmap_nasv_change_d'+str(depth)+'_annot.png', bbox_inches='tight', dpi=300)\n",
    "#     plt.clf() \n",
    "    \n",
    "#     ax = sns.heatmap(floored_data, fmt='.1f', yticklabels=True, linewidths=.5,\n",
    "#                     cmap = cmap)\n",
    "#     plt.savefig('outputs/'+comm_id+'/heatmap_nasv_change_d'+str(depth)+'.png', bbox_inches='tight', dpi=300)\n",
    "    \n",
    "    \n",
    "    plt.clf() \n",
    "    sns.set(rc={\"figure.figsize\":(7, 3)})\n",
    "    sns.set_style(\"ticks\")\n",
    "    ax=sns.lineplot(x = \"weekn\", y = ycolumn, data=sfd, hue=\"size_code\", palette=palette_dict)\n",
    "    sns.despine()\n",
    "    plt.ylabel(yaxislabel, fontsize=20)\n",
    "    plt.xlabel('Week', fontsize=20)\n",
    "    plt.legend(title='Size fraction')\n",
    "    sns.move_legend(ax, \"upper left\", bbox_to_anchor=(1, 1))\n",
    "    plt.savefig('outputs/'+comm_id+'/D'+str(depth)+'_adlineplot.png', dpi=200, bbox_inches=\"tight\")\n",
    "    \n",
    "    return anova, results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def upsetprep(comm, level, separated):\n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    cumulab = separated[['feature_frequency', 'depth', 'size_code', level]].copy()\n",
    "    cumulab1 = cumulab.groupby([level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "    for d in depths:\n",
    "        #make csv\n",
    "        sfd=separated[separated.depth==d]\n",
    "        \n",
    "        toptaxa = sfd[['feature_id', 'feature_frequency', 'Taxon', 'size_code', 'depth','weekn', level]].copy()\n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level, 'depth']).agg({'feature_frequency':sum})\n",
    "        topd = df_agg['feature_frequency'].groupby('size_code', group_keys=False).nlargest(10)\n",
    "        topd = topd.to_frame()\n",
    "        topd = topd.reset_index()\n",
    "\n",
    "        df_agg = df_agg.reset_index()\n",
    "        df_agg['set_name'] = df_agg['size_code']+df_agg['depth'].astype(str)\n",
    "    \n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='set_name', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "        resultpivot[resultpivot != 0] = 1\n",
    "        tosave = pd.merge(resultpivot, cumulab1, left_index=True, right_index=True)\n",
    "        tosave.to_csv('csvs/'+comm_id+'/'+level+'_d'+str(d)+'_relab.csv')\n",
    "        \n",
    "        \n",
    "        #make json\n",
    "        data = {\n",
    "            \"file\": \"https://raw.githubusercontent.com/dianahaider/size_fractions/main/csvs/\"+comm_id+'/'+level+'_d'+str(d)+'_relab.csv',\n",
    "            \"name\": comm + level,\n",
    "            \"header\": 0,\n",
    "            \"separator\": \",\",\n",
    "            \"skip\": 0,\n",
    "            \"meta\":[\n",
    "                {\"type\":\"id\", \"index\":0, \"name\":\"Name\"},\n",
    "                {\"type\":\"integer\", \"index\":4, \"name\":\"Rel. ab.\"}\n",
    "            ],\n",
    "            \"sets\": [\n",
    "                {\"format\": \"binary\", \"start\":1, \"end\": 3}\n",
    "            ]\n",
    "        }\n",
    "        \n",
    "        with open('json/'+comm_id+'/'+level+'_d'+str(d)+'.json', 'w') as f:\n",
    "            json.dump(data, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_per_fid(comm, separated, depth, fid):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "    \n",
    "    if 'SL' in separated['size_code'].unique():\n",
    "        #sizecode palette codes\n",
    "        sizecodes = ['S', 'L', 'W', 'SL']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "    \n",
    "    else:\n",
    "        #sizecode palette codes\n",
    "        sizecodes = ['S', 'L', 'W']\n",
    "        palette_colors = sns.color_palette()\n",
    "        palette_dict = {sizecode: color for sizecode, color in zip(sizecodes, palette_colors)}\n",
    "    \n",
    "    sfd=separated[separated.depth==depth]\n",
    "    sfd['weekfid'] = sfd[\"weekn\"].astype(str) + sfd[\"feature_id\"].astype(str)\n",
    "    sfd['avg_p_id'] = sfd['ratio'].groupby(sfd['weekfid']).transform('mean')\n",
    "    sfd['diff_p_id'] = sfd['ratio'] - sfd['avg_p_id']\n",
    "    \n",
    "    sfd_f=sfd[sfd.feature_id==fid]\n",
    "    \n",
    "    ttl = sfd_f['Taxon'].iloc[0]\n",
    "    \n",
    "    sns.set(rc={\"figure.figsize\":(7, 3)})\n",
    "    ax=sns.barplot(data=sfd_f, x=\"weekn\", y=\"diff_p_id\", hue=\"size_code\", palette=palette_dict)#, hue=\"size_code\")\n",
    "    sns.move_legend(ax, \"upper left\", bbox_to_anchor=(1, 1))\n",
    "    plt.title(ttl)\n",
    "    plt.ylabel('Ratio difference')\n",
    "    plt.xlabel('Week number')\n",
    "    plt.savefig('outputs/'+comm_id+'/D'+str(depth)+fid+'.png', dpi=200, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_ancom(separated, sfdclr, depth, ancomcol):\n",
    "    \n",
    "    sfd=separated[separated.depth==depth]\n",
    "\n",
    "        \n",
    "    df_ancom = sfd[['sampleid', ancomcol]].copy()\n",
    "    df_ancom = df_ancom.drop_duplicates()\n",
    "    df_ancom = df_ancom.set_index('sampleid')\n",
    "    \n",
    "    results = ancom(table=sfdclr, grouping=df_ancom[ancomcol])\n",
    "    \n",
    "    DAresults = results[0].copy()\n",
    "    DARejected_SC = DAresults.loc[DAresults['Reject null hypothesis'] == True]\n",
    "    DARejected_SC.sort_values(by=['W'])\n",
    "    \n",
    "    taxonomy = sfd[['feature_id', 'Confidence', 'Taxon', 'Phylum', 'Class', 'Family', 'Genus', 'Species']].copy()\n",
    "    taxonomy = taxonomy.drop_duplicates()\n",
    "    DARejected_SC_taxonomy = pd.merge(DARejected_SC, taxonomy, on=\"feature_id\", how=\"left\")\n",
    "    DARejected_SC_taxonomy.sort_values(by='W')\n",
    "    \n",
    "    prcentile = results[1].copy()\n",
    "    \n",
    "    return DARejected_SC_taxonomy, prcentile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subtitle = 'From Jan7 2022 to Apr27 2022'\n",
    "def plot_stackedbar_(df, labels, colors, title, subtitle, level):\n",
    "        \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    fields = df.columns.tolist()\n",
    "    \n",
    "    # figure and axis\n",
    "    fig, ax = plt.subplots(1, figsize=(8, 5))\n",
    "# plot bars\n",
    "    left = len(df) * [0]\n",
    "    for idx, name in enumerate(fields):\n",
    "        plt.barh(df.index, df[name], left = left, color=colors[idx])\n",
    "        left = left + df[name]\n",
    "# title and subtitle\n",
    "    plt.title(title, loc='left')\n",
    "    #plt.text(0, ax.get_yticks()[-1] + 0.75, subtitle)\n",
    "# legend\n",
    "    plt.legend(labels, bbox_to_anchor=([1, 1, 0, 0]), ncol=1, frameon=False)\n",
    "# remove spines\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "# format x ticks\n",
    "    xticks = np.arange(0,110,10)\n",
    "    xlabels = ['{}%'.format(i) for i in np.arange(0,101,10)]\n",
    "    plt.xticks(xticks, xlabels)\n",
    "# adjust limits and draw grid lines\n",
    "    plt.ylim(-0.5, ax.get_yticks()[-1] + 0.5)\n",
    "    ax.xaxis.grid(color='gray', linestyle='dashed')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.ylabel(\"Depth (m)\")\n",
    "    plt.savefig('outputs/'+comm_id+'/'+level+'alldepths_stacked_perc_weighted.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subtitle = 'From Jan7 2022 to Apr27 2022'\n",
    "def plot_stackedbar_p(df, labels, colors, title, subtitle, level):\n",
    "        \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    fields = df.columns.tolist()\n",
    "    \n",
    "    # figure and axis\n",
    "    fig, ax = plt.subplots(1, figsize=(8, 5))\n",
    "# plot bars\n",
    "    left = len(df) * [0]\n",
    "    for idx, name in enumerate(fields):\n",
    "        plt.barh(df.index, df[name], left = left, color=colors[idx])\n",
    "        left = left + df[name]\n",
    "# title and subtitle\n",
    "    plt.title(title, loc='left')\n",
    "    #plt.text(0, ax.get_yticks()[-1] + 0.75, subtitle)\n",
    "# legend\n",
    "    plt.legend(labels, bbox_to_anchor=([1, 1, 0, 0]), ncol=1, frameon=False)\n",
    "# remove spines\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "# format x ticks\n",
    "    xticks = np.arange(0,110,10)\n",
    "    xlabels = ['{}%'.format(i) for i in np.arange(0,101,10)]\n",
    "    plt.xticks(xticks, xlabels)\n",
    "# adjust limits and draw grid lines\n",
    "    plt.ylim(-0.5, ax.get_yticks()[-1] + 0.5)\n",
    "    ax.xaxis.grid(color='gray', linestyle='dashed')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.ylabel(\"Depth (m)\")\n",
    "    plt.savefig('outputs/'+comm_id+'/'+level+'alldepths_stacked_perc_weighted.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subtitle = 'From Jan7 2022 to Apr27 2022'\n",
    "def plot_stackedbar_p(df, labels, colors, title, subtitle, level):\n",
    "        \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    fields = df.columns.tolist()\n",
    "    \n",
    "    # figure and axis\n",
    "    fig, ax = plt.subplots(1, figsize=(8, 5))\n",
    "# plot bars\n",
    "    left = len(df) * [0]\n",
    "    for idx, name in enumerate(fields):\n",
    "        plt.barh(df.index, df[name], left = left, color=colors[idx])\n",
    "        left = left + df[name]\n",
    "# title and subtitle\n",
    "    plt.title(title, loc='left')\n",
    "    #plt.text(0, ax.get_yticks()[-1] + 0.75, subtitle)\n",
    "# legend\n",
    "    plt.legend(labels, bbox_to_anchor=([1, 1, 0, 0]), ncol=1, frameon=False)\n",
    "# remove spines\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "# format x ticks\n",
    "    xticks = np.arange(0,110,10)\n",
    "    xlabels = ['{}%'.format(i) for i in np.arange(0,101,10)]\n",
    "    plt.xticks(xticks, xlabels)\n",
    "# adjust limits and draw grid lines\n",
    "    plt.ylim(-0.5, ax.get_yticks()[-1] + 0.5)\n",
    "    ax.xaxis.grid(color='gray', linestyle='dashed')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.ylabel(\"Depth (m)\")\n",
    "    plt.savefig('outputs/'+comm_id+'/'+level+'alldepths_stacked_perc_weighted.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc_defrac(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'SF', 'NSF', 'Both', 'DFr'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df1 = resultpivot.copy()\n",
    "    \n",
    "        df = resultpivot[['L', 'S', 'W']].copy()\n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "        \n",
    "        DFr = df1[(df1['SL'] != 0)]\n",
    "        DFr = DFr[['SL']].copy()\n",
    "    \n",
    "        total = resultpivot.to_numpy().sum()\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = SF.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = Both.to_numpy().sum()/total *100\n",
    "    \n",
    "        Wonly_value = Wonly.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "        Lonly_value = Lonly.to_numpy().sum()/total *100\n",
    "        LS_value = LS.to_numpy().sum()/total *100\n",
    "        DFr_value = DFr.to_numpy().sum()/total *100\n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'SF'] = SF_value\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot.loc[d,'Both'] = Both_value\n",
    "        dfplot.loc[d,'DFr'] = DFr_value\n",
    "        \n",
    "        dfplot_unweighted.loc[d,'Depth'] = depths[d]\n",
    "        dfplot_unweighted.loc[d,'SF'] = len(Lonly) + len(Sonly) + len(LS)\n",
    "        dfplot_unweighted.loc[d,'NSF'] = len(Wonly)\n",
    "        dfplot_unweighted.loc[d,'Both'] = len(LW) + len(SW) + len(LSW)\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "        \n",
    "    return dfplot, level, dfplot_unweighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc_defrac_unweighted(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'SF', 'NSF', 'Both'])\n",
    "    dfplot_unweighted = pd.DataFrame(columns=['Depth', 'SF', 'NSF', 'Both'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df1 = resultpivot.copy()\n",
    "    \n",
    "        df = resultpivot[['L', 'S', 'W']].copy()\n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "    \n",
    "        total = len(resultpivot)\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = len(SF)/total *100\n",
    "        \n",
    "        Sonly_value = len(Sonly)/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = len(Both)/total *100\n",
    "    \n",
    "        Wonly_value = len(Wonly)/total *100\n",
    "        \n",
    "        Sonly_value = len(Sonly)/total *100\n",
    "        Lonly_value = len(Lonly)/total *100\n",
    "        LS_value = len(LS)/total *100\n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'SF'] = SF_value\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot.loc[d,'Both'] = Both_value\n",
    "        \n",
    "        dfplot_unweighted.loc[d,'Depth'] = depths[d]\n",
    "        dfplot_unweighted.loc[d,'SF'] = Lonly_value + Sonly_value + LS_value\n",
    "        dfplot_unweighted.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot_unweighted.loc[d,'Both'] = 100 - (Lonly_value + Sonly_value + LS_value + Wonly_value)\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "        \n",
    "    return dfplot, level, dfplot_unweighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'SF', 'NSF', 'Both'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df = resultpivot.copy()\n",
    "    \n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "    \n",
    "        total = df.to_numpy().sum()\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = SF.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = Both.to_numpy().sum()/total *100\n",
    "    \n",
    "        Wonly_value = Wonly.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "        Lonly_value = Lonly.to_numpy().sum()/total *100\n",
    "        LS_value = LS.to_numpy().sum()/total *100\n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'SF'] = SF_value\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot.loc[d,'Both'] = Both_value\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "        \n",
    "    return dfplot, level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc_SLNSF(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'Sonly', 'Lonly', 'LS', 'NSF'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df = resultpivot.copy()\n",
    "    \n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "    \n",
    "        total = df.to_numpy().sum()\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = SF.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = Both.to_numpy().sum()/total *100\n",
    "    \n",
    "        Wonly_value = Wonly.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "        Lonly_value = Lonly.to_numpy().sum()/total *100\n",
    "        LS_value = LS.to_numpy().sum()/total *100\n",
    "        \n",
    "        NewTotal = Sonly_value + Lonly_value + LS_value + Wonly_value\n",
    "        \n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'Sonly'] = Sonly_value\n",
    "        dfplot.loc[d,'Lonly'] = Lonly_value\n",
    "        dfplot.loc[d,'LS'] = LS_value\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "    dfplot_normalized = dfplot/NewTotal *100\n",
    "        \n",
    "    return dfplot, dfplot_normalized, level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc_LSW(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'NSF', 'LW', 'SW', 'LSW'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df = resultpivot.copy()\n",
    "    \n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "    \n",
    "        total = df.to_numpy().sum()\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = SF.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = Both.to_numpy().sum()/total *100\n",
    "    \n",
    "        Wonly_value = Wonly.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "        Lonly_value = Lonly.to_numpy().sum()/total *100\n",
    "        LS_value = LS.to_numpy().sum()/total *100\n",
    "        \n",
    "        LW_value = LW.to_numpy().sum()/total *100\n",
    "        SW_value = SW.to_numpy().sum()/total *100\n",
    "        LSW_value = LSW.to_numpy().sum()/total *100\n",
    "        \n",
    "        NewTotal = Sonly_value + Lonly_value + LS_value + Wonly_value\n",
    "        \n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot.loc[d,'LW'] = LW_value\n",
    "        dfplot.loc[d,'SW'] = SW_value\n",
    "        dfplot.loc[d,'LSW'] = LSW_value\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "    dfplot_normalized = dfplot/NewTotal *100\n",
    "        \n",
    "    return dfplot, dfplot_normalized, level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcperc_LS_W(comm, separated, level):\n",
    "    \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    depths = [1, 5, 10, 30, 60]\n",
    "    \n",
    "    level = level\n",
    "    \n",
    "    dfplot = pd.DataFrame(columns=['Depth', 'NSF', 'LW', 'SW'])\n",
    "    \n",
    "    for d in range(len(depths)):\n",
    "        sfd=separated[separated.depth==depths[d]]\n",
    "        toptaxa = sfd[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "        toptaxa = toptaxa.drop_duplicates()\n",
    "        df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "        df_agg = df_agg.reset_index()\n",
    "        resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "        resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "        df = resultpivot.copy()\n",
    "    \n",
    "        Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "        Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "        Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "        LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "        LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "        SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "        LSW = df[~(df == 0).any(axis=1)]\n",
    "    \n",
    "        total = df.to_numpy().sum()\n",
    "    \n",
    "        SFdf = Lonly, LS, Sonly\n",
    "        SF = pd.concat(SFdf)\n",
    "        SF_value = SF.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "    \n",
    "    \n",
    "        Bothdf = LW, LSW, SW\n",
    "        Both = pd.concat(Bothdf)\n",
    "        Both_value = Both.to_numpy().sum()/total *100\n",
    "    \n",
    "        Wonly_value = Wonly.to_numpy().sum()/total *100\n",
    "        \n",
    "        Sonly_value = Sonly.to_numpy().sum()/total *100\n",
    "        Lonly_value = Lonly.to_numpy().sum()/total *100\n",
    "        LS_value = LS.to_numpy().sum()/total *100\n",
    "        \n",
    "        LW_value = LW.to_numpy().sum()/total *100\n",
    "        SW_value = SW.to_numpy().sum()/total *100\n",
    "        LSW_value = LSW.to_numpy().sum()/total *100\n",
    "        \n",
    "        NewTotal = Sonly_value + Lonly_value + LS_value + Wonly_value\n",
    "        \n",
    "        \n",
    "        dfplot.loc[d,'Depth'] = depths[d]\n",
    "        dfplot.loc[d,'NSF'] = Wonly_value\n",
    "        dfplot.loc[d,'LW'] = LW_value\n",
    "        dfplot.loc[d,'SW'] = SW_value\n",
    "        \n",
    "\n",
    "        venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "        plt.savefig(\"outputs/\"+comm_id+\"/D\"+str(depths[d])+level+\"_venn.png\")\n",
    "        plt.clf()\n",
    "        plt.cla()\n",
    "        plt.close()\n",
    "    \n",
    "    dfplot['Depth'] = dfplot['Depth'].astype(str)\n",
    "    dfplot = dfplot.set_index('Depth')\n",
    "    dfplot_normalized = dfplot/NewTotal *100\n",
    "        \n",
    "    return dfplot, dfplot_normalized, level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subtitle = 'From Jan7 2022 to Apr27 2022'\n",
    "def plot_stackedbar_p_SLNSF(df, labels, colors, title, subtitle, level, xmax=110, xtick=10):\n",
    "        \n",
    "    if comm == '16S':\n",
    "        comm_id = '02-PROKs'\n",
    "    else:\n",
    "        comm_id = '02-EUKs'\n",
    "        \n",
    "    fields = df.columns.tolist()\n",
    "    \n",
    "    name =[x for x in globals() if globals()[x] is df][0]\n",
    "    \n",
    "    # figure and axis\n",
    "    fig, ax = plt.subplots(1, figsize=(8, 5))\n",
    "# plot bars\n",
    "    left = len(df) * [0]\n",
    "    for idx, name in enumerate(fields):\n",
    "        plt.barh(df.index, df[name], left = left, color=colors[idx])\n",
    "        left = left + df[name]\n",
    "# title and subtitle\n",
    "    plt.title(title, loc='left')\n",
    "    #plt.text(0, ax.get_yticks()[-1] + 0.75, subtitle)\n",
    "# legend\n",
    "    plt.legend(labels, bbox_to_anchor=([1, 1, 0, 0]), ncol=1, frameon=False)\n",
    "# remove spines\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "# format x ticks\n",
    "    xticks = np.arange(0,xmax,xtick)\n",
    "    xlabels = ['{}%'.format(i) for i in np.arange(0,xmax,xtick)]\n",
    "    plt.xticks(xticks, xlabels)\n",
    "# adjust limits and draw grid lines\n",
    "    plt.ylim(-0.5, ax.get_yticks()[-1] + 0.5)\n",
    "    ax.xaxis.grid(color='gray', linestyle='dashed')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.ylabel(\"Depth (m)\")\n",
    "    plt.savefig('outputs/'+comm_id+'/'+level+'alldepths_stacked_perc_weighted'+name+'.png', dpi=200, bbox_inches=\"tight\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generate a dataframe from all specified amplicon\n",
    "df, comm = consolidate_tables('16S')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#only if needed redefine the metadata file used to create merged\n",
    "#all_md = md.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = merge_metadata(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separated = pick_metadata(merged)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove chloroplast and cyanobacteria from 16S\n",
    "chloroplast = separated[separated.Taxon.str.contains(\"Cyanobacteria\")]\n",
    "separated = separated[~separated.Taxon.str.contains(\"Cyanobacteria\")]\n",
    "separated = separated[~separated.Taxon.str.contains(\"Chloroplast\")]\n",
    "#count the number of features removed\n",
    "chloroplast['feature_id'].nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separated = separated.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following code to generate \"newseparated\" which is the union of small and large size fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make sure all size codes are indicated\n",
    "all_md[\"size_code\"] = all_md[\"sampleid\"].str.extract(r'[1-9][0-9]?[A-E]([L-S])')\n",
    "all_md[\"size_code\"] = all_md[\"size_code\"].fillna('W')\n",
    "\n",
    "#only keep values from weeks 1 to 16\n",
    "sep_SL = all_md[all_md.size_code != \"W\"]\n",
    "sep_SL = sep_SL.drop(sep_SL[sep_SL.weekn > 16].index)\n",
    "\n",
    "#sum [DNA] of small and large size fractions\n",
    "sep_SL['[DNAt]'] = sep_SL.groupby(['weekn', 'depth'])['[DNA]ng/ul'].transform('sum')\n",
    "\n",
    "#separate small and large size fraction\n",
    "sep_S = sep_SL[sep_SL.size_code == 'S']\n",
    "sep_L = sep_SL[sep_SL.size_code == 'L']\n",
    "\n",
    "#calculate DNA proportion per size fraction\n",
    "sep_SL['DNApr'] = sep_SL['[DNA]ng/ul']/sep_SL['[DNAt]']\n",
    "\n",
    "#merge with separated on common columns to get corresponding rel. abundances\n",
    "sep_SL = sep_SL[['sampleid', 'DNApr', '[DNAt]']].copy()\n",
    "sepSLRA = pd.merge(separated, sep_SL, on=['sampleid'], how='left') #all_md is the metadata file\n",
    "\n",
    "#exclude ASVs from the whole water\n",
    "sep_SLRA = sepSLRA[separated.size_code != \"W\"]\n",
    "\n",
    "#calculate corrected per sample ratio, and corrected feature frequency of de-fractionated samples\n",
    "sep_SLRA['Newfeature_frequency'] = sep_SLRA['feature_frequency'] * sep_SLRA['DNApr']\n",
    "sep_SLRA['Newff'] = sep_SLRA.groupby(['feature_id', 'weekn', 'depth'])['Newfeature_frequency'].transform('sum')\n",
    "\n",
    "\n",
    "#sep_SLRA = sep_SLRA.drop(['sampleid', 'size_code'], axis=1)\n",
    "sep_SLRA['sampleid'] = \"BB22.\" + sep_SLRA['weekn'].astype(str) + sep_SLRA['depth_code'] + \"SL\"\n",
    "\n",
    "#uncomment the line below if keeping small and large original sample\n",
    "#sep_SLRA['size_code'] = sep_SLRA['size_code'] + '-DFr'\n",
    "\n",
    "#uncomment the line above if merging smallandlarge\n",
    "sep_SLRA['size_code'] = 'SL'\n",
    "\n",
    "#drop unecessary columns which might rise merging conflicts\n",
    "sep_SLRA = sep_SLRA.drop(['feature_frequency', 'Total', 'ratio', 'nASVs', 'weekdepth', 'avg',\n",
    "                          'diff', 'extraction_date', '[DNA]ng/ul', 'A260/280', 'A260/230',\n",
    "                          'Newfeature_frequency'], axis=1)\n",
    "sep_SLRA.rename(columns={'Newff':'feature_frequency'}, inplace=True)\n",
    "sep_SLRA = sep_SLRA.drop_duplicates()\n",
    "\n",
    "#recalculate ratios\n",
    "sep_SLRA['Total'] = sep_SLRA['feature_frequency'].groupby(sep_SLRA['sampleid']).transform('sum')\n",
    "sep_SLRA['ratio'] = sep_SLRA['feature_frequency']/sep_SLRA['Total']\n",
    "sep_SLRA['nASVs'] = sep_SLRA['feature_id'].groupby(sep_SLRA['sampleid']).transform('nunique')\n",
    "\n",
    "sep_SLRA = sep_SLRA.drop_duplicates()\n",
    "\n",
    "#make new df dependingg on plotting needs\n",
    "sep_WO = separated[separated.size_code == \"W\"]\n",
    "sep_WO = sep_WO.drop_duplicates()\n",
    "\n",
    "sep_S = separated[separated.size_code == \"S\"]\n",
    "sep_L = separated[separated.size_code == \"L\"]\n",
    "\n",
    "\n",
    "sep_WO.reset_index(inplace=True, drop=True)\n",
    "sep_SLRA.reset_index(inplace=True, drop=True)\n",
    "\n",
    "#newseparated = pd.concat([sep_SLRA.reset_index(drop=True), sep_WO.reset_index(drop=True)], axis=0).reset_index(drop=True)\n",
    "newseparated = pd.concat([sep_SLRA, sep_WO, sep_L, sep_S], ignore_index=True)\n",
    "\n",
    "newseparated['weekdepth'] = newseparated[\"weekn\"].astype(str) + newseparated[\"depth\"].astype(str)\n",
    "newseparated['avg'] = newseparated['nASVs'].groupby(newseparated['weekdepth']).transform('mean')\n",
    "newseparated['diff'] = newseparated['nASVs'] - newseparated['avg']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Richness analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd1 = separated[['sampleid','size_code', 'weekn', 'nASVs', 'depth']].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#group the dataframe with all features to obtain either the mean or std of number of features per size fraction\n",
    "sfd1.groupby(['size_code']).std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run the visualisations for alpha diversity and run pairwise t-tests between size fractions for richness values\n",
    "anova, results = boxplot_depth(separated, comm, 60, 'nASVs', 'Number of ASVs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple regression analysis between pairs of size fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data prep\n",
    "if comm == '16S':\n",
    "    comm_id = '02-PROKs'\n",
    "else:\n",
    "    comm_id = '02-EUKs'\n",
    "\n",
    "depth = 1\n",
    "\n",
    "d1 = newseparated.loc[newseparated['depth'] == depth]\n",
    "forpl = d1[['ratio', 'feature_id', 'sampleid', 'weekn', 'depth', 'size_code', 'Phylum', 'Family']].copy()\n",
    "slwplot = forpl.pivot_table(index=[\"feature_id\", \"depth\", 'weekn','Phylum', 'Family'], columns=\"size_code\", values='ratio').fillna(0)\n",
    "slwplot = slwplot.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q75, q25 = np.percentile(slwplot['L'], [75 ,25])\n",
    "iqr = q75 - q25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def outlier(df):\n",
    "    new_df = df.copy()\n",
    "    numeric_cols = ['L', 'S', 'SL', 'W']\n",
    "    \n",
    "    q1 = np.percentile(new_df[numeric_cols],25, axis=0)\n",
    "    q3 = np.percentile(new_df[numeric_cols],75, axis=0)\n",
    "    IQR = q3 - q1\n",
    "    lower_limit = q1 - (1.5*IQR)\n",
    "    upper_limit = q3 + (1.5*IQR)\n",
    "    mask = (new_df[numeric_cols] < lower_limit) | (new_df[numeric_cols] > upper_limit)\n",
    "    new_df[numeric_cols] = new_df[numeric_cols].mask(mask)\n",
    "    return new_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newslw = outlier(slwplot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newslw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newslw['L'].isna().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rerun the regression by removing the 4 pathogens found in weeks 10 and 11\n",
    "slwplot.sort_values('W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slwplot.drop([35,36,2353,2200], axis=0, inplace=True)\n",
    "\n",
    "slwplot.sort_values('W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slwplot = slwplot.loc[slwplot['feature_id'] != '4cb5c79abe0a3b611671b1f356a86f19']\n",
    "slwplot = slwplot.loc[slwplot['feature_id'] != '51d96e96a3350beedece3878d6d0b3e7']\n",
    "slwplot = slwplot.loc[slwplot['feature_id'] != '015600cdcdfdc7333d251c9c9160eb72']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = slwplot['S']\n",
    "X = slwplot['W']\n",
    "\n",
    "X = sm.add_constant(X)\n",
    "\n",
    "#fit the model\n",
    "model = sm.OLS(Y, X, missing='drop')\n",
    "model_result = model.fit()\n",
    "model_result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(model_result.resid);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu, std = stats.norm.fit(model_result.resid)\n",
    "mu, std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "# plot the residuals\n",
    "sns.histplot(x=model_result.resid, ax=ax, stat=\"density\", linewidth=0, kde=True)\n",
    "ax.set(title=\"Distribution of residuals\", xlabel=\"residual\")\n",
    "\n",
    "# plot corresponding normal curve\n",
    "xmin, xmax = plt.xlim() # the maximum x values from the histogram above\n",
    "x = np.linspace(xmin, xmax, 100) # generate some x values\n",
    "p = stats.norm.pdf(x, mu, std) # calculate the y values for the normal curve\n",
    "sns.lineplot(x=x, y=p, color=\"orange\", ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x=model_result.resid, showmeans=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm.qqplot(model_result.resid, line='s');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(6, 5)\n",
    "\n",
    "sns.set_style('white')\n",
    "fig = sm.graphics.plot_fit(model_result,1, vlines=False, ax=ax)\n",
    "ax.set_ylabel(\"Defractionated\")\n",
    "ax.set_xlabel(\"Whole\")\n",
    "ax.set_title(\"Fitted values linear regression\")\n",
    "\n",
    "plt.savefig('outputs/'+comm_id+'/asv_'+str(depth)+'_scatter_RL.png', dpi=200, bbox_inches=\"tight\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_result.fittedvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_max = Y.max()\n",
    "Y_min = Y.min()\n",
    "\n",
    "plt.figure(figsize=(5, 4))\n",
    "\n",
    "\n",
    "ax = sns.scatterplot(x=model_result.fittedvalues, y=Y)\n",
    "ax.set(ylim=(Y_min, Y_max))\n",
    "ax.set(xlim=(Y_min, Y_max))\n",
    "ax.set_xlabel(\"Predicted value\")\n",
    "ax.set_ylabel(\"Observed value\")\n",
    "\n",
    "X_ref = Y_ref = np.linspace(Y_min, Y_max, 100)\n",
    "plt.plot(X_ref, Y_ref, color='red', linewidth=1)\n",
    "\n",
    "plt.savefig('outputs/'+comm_id+'/asv_'+str(depth)+'_scatterRL.png', dpi=200, bbox_inches=\"tight\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "investigate temporal pattern thorugh all depths of outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slwplot.sort_values('W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slwplot.loc[slwplot['feature_id'] == '75ceeaa937c64399438614ca3706cf2a']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "supsel = newseparated.loc[newseparated['feature_id'] == '75ceeaa937c64399438614ca3706cf2a'].sort_values('feature_frequency')\n",
    "supsel['depth'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "supsel.sort_values('weekn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(data = supsel, x = 'weekn', y = 'ratio', hue='depth')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#draw interactive plotly to identify outliers\n",
    "import plotly.express as px\n",
    "\n",
    "df = px.data.tips()\n",
    "fig = px.scatter(slwplot, x=\"W\", y=\"SL\", color=\"weekn\", trendline=\"ols\")\n",
    "fig.show()\n",
    "\n",
    "results = px.get_trendline_results(fig)\n",
    "print(results)\n",
    "\n",
    "#results.query(\"weekn ==  and Phylum == \").px_fit_results.iloc[0].summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate log2 fold change per feature at the phylum level of abundance between size fractions to identify which taxonomic group is driving the simple regression off x=y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add pseudo count for log-calculations and zero divisions\n",
    "slwplot['SL'] = slwplot['SL'] + 0.0000001\n",
    "slwplot['W'] = slwplot['W'] + 0.0000001\n",
    "\n",
    "#calculate log2 fold change\n",
    "slwplot['OR'] = (slwplot['W'] - slwplot['SL']) / slwplot['SL']\n",
    "slwplot['fold_change'] = slwplot['W']/slwplot['SL']\n",
    "slwplot['log2_fold_change'] = np.log2(slwplot['fold_change'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make dual plot of log2FC and mean relative abundances side by side including error bars in the plots\n",
    "data = slwplot[['log2_fold_change','Phylum', 'W', 'SL']].copy()\n",
    "\n",
    "data['Phylum'] = data['Phylum'].map(lambda x: x.lstrip('p__'))\n",
    "\n",
    "# Group by index labels and take the means and standard deviations for each group\n",
    "data['avg_W'] = data['W'].groupby(data['Phylum']).transform('mean')\n",
    "data['std_W'] = data['W'].groupby(data['Phylum']).transform('std')\n",
    "data['avg_SL'] = data['SL'].groupby(data['Phylum']).transform('mean')\n",
    "data['std_SL'] = data['SL'].groupby(data['Phylum']).transform('std')\n",
    "data['means'] = data['log2_fold_change'].groupby(data['Phylum']).transform('mean')\n",
    "data['stds'] = data['log2_fold_change'].groupby(data['Phylum']).transform('std')\n",
    "\n",
    "data['positive'] = data['means'] > 0\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2, sharey=True, figsize=(8, 10))\n",
    "\n",
    "axes[0].barh(data['Phylum'], data['means'],\n",
    "         xerr = data['stds'],\n",
    "         error_kw=dict(lw=0.5, capsize=1, capthick=0.5),\n",
    "         color=data.positive.map({True: 'g', False: 'r'}))\n",
    "\n",
    "axes[1].barh(data['Phylum'], data['avg_W'],\n",
    "            xerr = data['std_W'],\n",
    "         error_kw=dict(lw=0.5, capsize=1, capthick=0.5))\n",
    "\n",
    "#axes[2].barh(data['Phylum'], data['avg_SL'],\n",
    "#            xerr = data['std_SL'],\n",
    "#         error_kw=dict(lw=0.5, capsize=1, capthick=0.5))\n",
    "\n",
    "plt.gca().invert_yaxis()\n",
    "\n",
    "\n",
    "plt.savefig('outputs/'+comm_id+'/log2foldchange_d'+str(depth)+'.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#bar plot of log2FC per phylum without error bars\n",
    "data = slwplot[['log2_fold_change','Phylum']].copy()\n",
    "data['Phylum'] = data['Phylum'].map(lambda x: x.lstrip('p__'))\n",
    "\n",
    "# Group by index labels and take the means and standard deviations for each group\n",
    "#data['means'] = data['log2_fold_change'].groupby(data['Phylum']).transform('mean')\n",
    "#data['stds'] = data['log2_fold_change'].groupby(data['Phylum']).transform('std')\n",
    "\n",
    "data['positive'] = data['log2_fold_change'] > 0\n",
    "\n",
    "plt.figure(figsize=(6,11))\n",
    "plt.barh(data['Phylum'], data['log2_fold_change'],\n",
    "         color=data.positive.map({True: 'g', False: 'r'}))\n",
    "plt.gca().invert_yaxis()\n",
    "\n",
    "plt.savefig('outputs/'+comm_id+'/log2foldchange_d'+str(depth)+'_noerr.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "not_in_W = slwplot.loc[slwplot['W'] == 0]\n",
    "#not_in_W['feature_id'].nunique()\n",
    "not_in_W_from_S = not_in_W.loc[not_in_W['L'] == 0]\n",
    "#not_in_W_from_S['feature_id'].nunique()\n",
    "\n",
    "#calculate the percentage of features coming from the small size fraction that are absent from the whole  \n",
    "not_in_W_from_S['feature_id'].nunique()/not_in_W['feature_id'].nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(data=not_in_W, x='S', y='SL', hue='weekn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### run lmplot for each pair of sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loglog plot\n",
    "depths = [1,5,10,30,60]\n",
    "\n",
    "for depth in depths:\n",
    "    d1 = newseparated.loc[newseparated['depth'] == depth]\n",
    "    forpl = d1[['ratio', 'feature_id', 'sampleid', 'weekn', 'depth', 'size_code', 'Phylum']].copy()\n",
    "    slwplot = forpl.pivot_table(index=[\"feature_id\", \"depth\", 'weekn', 'Phylum'], columns=\"size_code\", values='ratio').fillna(0)\n",
    "    slwplot = slwplot.reset_index()\n",
    "    sns.set_style(\"white\")\n",
    "\n",
    "    #slwplot = slwplot.loc[slwplot['weekn'] = 10]\n",
    "    #slwplot = slwplot.loc[slwplot['weekn'] = 11]\n",
    "\n",
    "    slwplot = slwplot.rename(columns={\"depth\": \"Depth\"})\n",
    "    slwplot[\"weekn\"] = pd.to_numeric(slwplot[\"weekn\"])\n",
    "    g = sns.scatterplot(x=\"W\", y=\"SL\", data=slwplot, palette=['black'])#, hue='Phylum', alpha=0.6) #, hue=\"weekn\");\n",
    "    \n",
    "    #uncomment for log-log\n",
    "    #ax.set(xscale=\"log\", yscale=\"log\")\n",
    "    g.set_ylabel(\"Defractionated\",fontsize=15)\n",
    "    g.set_xlabel(\"Unfractionated\",fontsize=15)\n",
    "    g.tick_params(labelsize=12)\n",
    "    #g.legend(loc='center left', bbox_to_anchor=(1.25, 0.5), ncol=3)\n",
    "    plt.legend([],[], frameon=False)\n",
    "    \n",
    "    #plt.savefig('outputs/lmplot_'+comm+str(depth)+'WSL.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phylogenetic analysis: taxonomic bar plots of relative abundance per depth and size fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get a list of top taxa to provide the palette for the visualisation\n",
    "toptaxa = newseparated[['feature_id', 'feature_frequency', 'Taxon', 'size_code', 'depth','weekn', 'Phylum']].copy()\n",
    "toptaxa = toptaxa.drop_duplicates()\n",
    "df_agg = toptaxa.groupby(['size_code','Phylum', 'depth']).agg({'feature_frequency':sum})\n",
    "topd = df_agg['feature_frequency'].groupby(['size_code', 'depth'], group_keys=False).nlargest(10)\n",
    "topd = topd.to_frame()\n",
    "topd = topd.reset_index()\n",
    "listoftop = topd['Phylum'].unique()\n",
    "\n",
    "#set a palette for the toptaxa\n",
    "hex_colors_dic = {}\n",
    "rgb_colors_dic = {}\n",
    "hex_colors_only = []\n",
    "for name, hex in matplotlib.colors.cnames.items():\n",
    "    hex_colors_only.append(hex)\n",
    "    hex_colors_dic[name] = hex\n",
    "    rgb_colors_dic[name] = matplotlib.colors.to_rgb(hex)\n",
    "    \n",
    "palette_dict = {taxon: color for taxon, color in zip(listoftop, px.colors.sequential.Plasma)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phyld, top10d = taxbarplot(newseparated, 'Phylum', 1, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Top taxon longitudinal analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# newseparated.loc[(newseparated['weekn']==12)&\n",
    "#                  (newseparated['depth']==10)]\n",
    "#                & (newseparated['size_code']=='W')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd=newseparated\n",
    "toptaxa = sfd[['feature_id', 'feature_frequency', 'Taxon', 'size_code', 'depth','weekn', 'Class']].copy()\n",
    "toptaxa = toptaxa.drop_duplicates()\n",
    "df_agg = toptaxa.groupby(['size_code','Class', 'depth', 'weekn']).agg({'feature_frequency':sum})\n",
    "topd = df_agg['feature_frequency'].groupby(['size_code', 'depth','weekn'], group_keys=False).nlargest(1)\n",
    "topd = topd.to_frame()\n",
    "topd = topd.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd.loc[topd['weekn'] == 12]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list the unique top taxa\n",
    "topd['Class'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_dic = {\n",
    "    #\"c__Cyanobacteriia\": 1,\n",
    "    \"c__OM190\": 3,\n",
    "    \"c__Bacteroidia\": 2,\n",
    "    \"c__Gammaproteobacteria\": 5,\n",
    "    \"c__Alphaproteobacteria\": 4}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make a season column\n",
    "topd['comm_type'] = ''\n",
    "\n",
    "for tx, ctype in type_dic.items():\n",
    "    topd.loc[topd['Class'] == tx, 'comm_type'] = ctype\n",
    "    \n",
    "topd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd[\"sc_weekn\"] = topd[\"depth\"].astype(str) + topd[\"size_code\"]\n",
    "topd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd = topd.sort_values(['depth', 'size_code'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd.loc[topd['sc_weekn'] == '10W']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topdlist = topd['sc_weekn'].tolist()\n",
    "\n",
    "def uniqlist(seq):\n",
    "    seen = set()\n",
    "    seen_add = seen.add\n",
    "    return [x for x in seq if not (x in seen or seen_add(x))]\n",
    "\n",
    "mylist = uniqlist(topdlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glue = topd.pivot(index=\"sc_weekn\", columns=\"weekn\", values=\"comm_type\")\n",
    "glue = glue.reindex(mylist)\n",
    "glue = glue[glue.columns].astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "cmap_dict = {#1: '#77AADD',\n",
    "             2: '#EEDD88', 3: '#99DDFF', 4: '#BBCC33', 5:'#DDDDDD'}\n",
    "cmap = ListedColormap([cmap_dict[i] for i in range(2,6,1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "plt.figure(figsize=(5, 5))\n",
    "\n",
    "\n",
    "ax = sns.heatmap(glue, fmt='f', yticklabels=True, linewidths=.5, cmap=cmap)\n",
    "ax.axhline(4, ls='--')\n",
    "ax.axhline(8, ls='--')\n",
    "ax.axhline(12, ls='--')\n",
    "ax.axhline(16, ls='--')\n",
    "\n",
    "ax.set_xticks(range(1, 16, 4))\n",
    "\n",
    "plt.savefig('outputs/heatmap_'+comm+'top1clasno_chloroplast.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Longitudinal analysis of top 3 taxa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run these lines only for the 18S rRNA analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#18S has many more classes as top values to\n",
    "top518s = topd['Class'].value_counts()[:4].index.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd.loc[~topd[\"Class\"].isin(top518s), \"Class\"] = \"Other\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newdf = topd.groupby(['size_code', 'depth','weekn'])['Class'].apply(lambda x: list(set(x)))\n",
    "newdf = newdf.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = newdf.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = newdf.Class.sort_values().apply(lambda x: sorted(x))\n",
    "result = pd.DataFrame(result).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.Class.apply(tuple).unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd['Class'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_dic = {\n",
    "    'c__Cryptophyceae':1, 'c__Dinophyceae':2, 'c__Prymnesiophyceae':3,\n",
    "       'c__Mediophyceae':4, 'c__Monogononta':5, 'c__Tentaculata':6,\n",
    "       'c__Maxillopoda':7, 'c__Thecofilosea':8, 'Unassigned':9, 'c__Insecta':10,\n",
    "       'c__MAST-2':11, 'c__Hydrozoa':12, 'c__Syndiniales':13,\n",
    "       'c__Intramacronucleata':14, 'c__Pucciniomycetes':15,\n",
    "       'c__Mamiellophyceae':16, 'c__Raphidophyceae':17, 'c__MAST-1A':18,\n",
    "       'c__Bicosoecida':19, 'c__Polychaeta':20, 'c__Tremellomycetes':21,\n",
    "       'c__Embryophyta':22, 'c__Dothideomycetes':23, 'c__Incertae_Sedis':24,\n",
    "       'c__MAST-7A':25\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type_dic = {'Other':5, 'c__Dinophyceae':1, 'c__Mediophyceae':2,\n",
    "       'c__Intramacronucleata':3, 'c__Syndiniales':4}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result['liststring'] = result['Class'].apply(lambda x: ','.join(map(str, x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result['liststring'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use this dic if looking at community type by 3 top taxa\n",
    "type_dic = {'c__Bacteroidia,c__Cyanobacteriia,c__OM190':1,\n",
    "       'c__Gammaproteobacteria,c__Bacteroidia,c__Cyanobacteriia':2,\n",
    "       'c__Alphaproteobacteria,c__Bacteroidia,c__Cyanobacteriia':3,\n",
    "       'c__Bacteroidia,c__Cyanobacteriia,c__Planctomycetes':4,\n",
    "       'c__Bacteroidia,c__Alphaproteobacteria,c__OM190':5,\n",
    "       'c__Bacteroidia,c__OM190,c__Planctomycetes':6,\n",
    "       'c__Gammaproteobacteria,c__Bacteroidia,c__OM190':7,\n",
    "       'c__Gammaproteobacteria,c__Bacteroidia,c__Planctomycetes':8,\n",
    "       'c__Cyanobacteriia,c__OM190,c__Planctomycetes':9,\n",
    "       'c__Alphaproteobacteria,c__Gammaproteobacteria,c__Cyanobacteriia':10,\n",
    "       'c__Gammaproteobacteria,c__Bacteroidia,c__Alphaproteobacteria':11,\n",
    "       'c__Bacteroidia,c__Alphaproteobacteria,c__Cyanobacteriia':12,\n",
    "       'c__Nitrososphaeria,c__Gammaproteobacteria,c__Bacteroidia':13,\n",
    "       'c__Gammaproteobacteria,c__Alphaproteobacteria':14,\n",
    "       'c__Cyanobacteriia':15}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make a season column\n",
    "topd['comm_type'] = ''\n",
    "\n",
    "for tx, ctype in type_dic.items():\n",
    "    topd.loc[topd['Class'] == tx, 'comm_type'] = ctype\n",
    "    \n",
    "topd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result['comm_type'] = ''\n",
    "\n",
    "for tx, ctype in type_dic.items():\n",
    "    result.loc[result['liststring'] == tx, 'comm_type'] = ctype\n",
    "    \n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd = result.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd[\"sc_weekn\"] = topd[\"depth\"].astype(str) + topd[\"size_code\"]\n",
    "topd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topd = topd.sort_values(['depth', 'size_code'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topdlist = topd['sc_weekn'].tolist()\n",
    "\n",
    "def uniqlist(seq):\n",
    "    seen = set()\n",
    "    seen_add = seen.add\n",
    "    return [x for x in seq if not (x in seen or seen_add(x))]\n",
    "\n",
    "mylist = uniqlist(topdlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glue = topd.pivot(index=\"sc_weekn\", columns=\"weekn\", values=\"comm_type\")\n",
    "glue = glue.reindex(mylist)\n",
    "glue = glue[glue.columns].astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " '#77AADD', '#EE8866', '#EEDD88', '#FFAABB', '#99DDFF', '#44BB99', '#BBCC33', '#AAAA00', '#DDDDDD'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "cmap_dict = {1: '#77AADD', 2: '#EEDD88', 3: '#99DDFF', 4: '#BBCC33', 5:'#DDDDDD'}\n",
    "cmap = ListedColormap([cmap_dict[i] for i in range(1,6,1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap_dict ={1:'#125A56', 2:'#00767B', 3:'#238F9D', 4:'#42A7C6', 5:'#60BCE9',\n",
    "            6:'#9DCCEF', 7:'#C6DBED', 8:'#DEE6E7', 9:'#ECEADA', 10:'#F0E6B2',\n",
    "            11:'#F9D576', 12:'#FFB954', 13:'#FD9A44', 14:'#F57634', 15:'#E94C1F'}#, '#D11807', '#A01813'.\n",
    "cmap = ListedColormap([cmap_dict[i] for i in range(1,16,1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "plt.figure(figsize=(5, 5))\n",
    "\n",
    "\n",
    "ax = sns.heatmap(glue, fmt='f', yticklabels=True, linewidths=.5, cmap=cmap)\n",
    "ax.axhline(4, ls='--')\n",
    "ax.axhline(8, ls='--')\n",
    "ax.axhline(12, ls='--')\n",
    "ax.axhline(16, ls='--')\n",
    "\n",
    "ax.set_xticks(range(1, 16, 4))\n",
    "\n",
    "\n",
    "plt.savefig('outputs/heatmap_'+comm+'top1class_reducde.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    sfd=separated[separated.depth==depth]\n",
    "    toptaxa = sfd[['feature_id', 'feature_frequency', 'Taxon', 'size_code', 'depth','weekn', level]].copy()\n",
    "    toptaxa = toptaxa.drop_duplicates()\n",
    "    df_agg = toptaxa.groupby(['size_code',level, 'depth']).agg({'feature_frequency':sum})\n",
    "    topd = df_agg['feature_frequency'].groupby('size_code', group_keys=False).nlargest(topn)\n",
    "    topd = topd.to_frame()\n",
    "    topd = topd.reset_index()\n",
    "\n",
    "\n",
    "    df_agg = df_agg.reset_index()\n",
    "    df_agg['set_name'] = df_agg['size_code']+df_agg['depth'].astype(str)\n",
    "    \n",
    "    cumulab = separated[['feature_frequency', 'depth', 'size_code', 'Genus']].copy()\n",
    "    cumulab1 = cumulab.groupby(['Genus']).agg({'feature_frequency':sum})\n",
    "\n",
    "    resultpivot = df_agg.pivot_table(index=level, columns='set_name', values='feature_frequency')\n",
    "    resultpivot = resultpivot.fillna(0)\n",
    "    resultpivot[resultpivot != 0] = 1\n",
    "    tosave = pd.merge(resultpivot, cumulab1, left_index=True, right_index=True)\n",
    "    tosave.to_csv(level+'_'+str(depth)+'16S_relab.csv')\n",
    "    \n",
    "    top10d_list = topd[level].unique()\n",
    "    top10d = sfd.copy()\n",
    "    top10d.loc[~top10d[level].isin(top10d_list), level] = 'Other' #isnot in top list\n",
    "    phyld = top10d.groupby(['size_code','weekn', level])['ratio'].sum()\n",
    "    phyld = phyld.reset_index()\n",
    "\n",
    "\n",
    "    fig = px.bar(phyld, x=\"size_code\", y=\"ratio\", facet_col=\"weekn\", color=level, labels={\n",
    "                     \"feature_frequency\": \"Relative abundance\",\n",
    "                     \"size_code\": \"\",\n",
    "                     \"weekn\": \"w\"}, color_discrete_map=palette_dict)\n",
    "    fig.update_xaxes(type='category', dtick=1)\n",
    "    fig.update_layout(\n",
    "        title=\"Relative abundance of top 10\" + level + 'observed at Depth' + str(depth),\n",
    "        yaxis_title=\"Relative abundance\",\n",
    "        xaxis_title=\"Size fraction\",\n",
    "        legend_title=level,\n",
    "        font=dict(size=8)\n",
    "    )\n",
    "\n",
    "    fig.show()\n",
    "    #fig.write_image(\"outputs/fig1.png\")\n",
    "    #fig.to_image(format=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phyld, top10d = taxbarplot(newseparated, 'Class', 5, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df2 = plot_df2.drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df3 = plot_df2.set_index('sampleid')\n",
    "plot_df3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "permanova2 = permanova(dm, plot_df3, 'Size code')\n",
    "permanova2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_matrix2 = distance_matrix.reset_index()\n",
    "idedup = distance_matrix2['samples'].to_list()\n",
    "dm = DistanceMatrix(distance_matrix, ids=idedup)\n",
    "df123 = dm.to_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df123.to_csv('distance_matrix_5m16s.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df2.to_csv('METADATAtiny.txt', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca, pca_features, sfdclr = pcaplot(newseparated, 5, comm, 'Size code', 'DFr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newmetadata = newseparated[['sampleid', 'weekn', 'size_code', 'depth', 'depth_code', 'month_name']].copy()\n",
    "newmetadata = newmetadata.drop_duplicates()\n",
    "newmetadata.to_csv('METADATA.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skbio.stats.distance import permanova"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level = 'feature_id'\n",
    "if level == 'feature_id':\n",
    "    id = 'ASV'\n",
    "else:\n",
    "    id = level\n",
    "\n",
    "subtitile = 'subtitle'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplot, level = calcperc(comm, separated, level)\n",
    "# variables\n",
    "labels = ['S ∩ W', 'L ∩ W','W ∩ (S ∩ L)']\n",
    "colors = ['#1D2F6F', '#8390FA', '#C49CD3']\n",
    "title = 'Weighted proportion of shared '+id+' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p(dfplot, labels, colors, title, subtitle, level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplot, level, dfplot_unweighted = calcperc_defrac_unweighted(comm, newseparated, level)\n",
    "# variables\n",
    "labels = ['Size fractionated samples', 'Not size fractionated samples','Both']\n",
    "colors = ['#1D2F6F', '#8390FA', '#6EAF46']\n",
    "title = 'Unweighted proportion of shared '+id+' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p(dfplot, labels, colors, title, subtitle, level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplot, level = calcperc_defrac(comm, newseparated, level)\n",
    "# variables\n",
    "labels = ['SF samples', 'NSF samples','Both', 'DFr']\n",
    "colors = ['#1D2F6F', '#8390FA', '#6EAF46', '#de282e']\n",
    "title = 'Weighted proportion of shared '+id+' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p(dfplot, labels, colors, title, subtitle, level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplotSLNSF, dfplot_normalized, level = calcperc_SLNSF(comm, separated, level)\n",
    "# variables\n",
    "labels = ['S', 'L','S ∩ L', 'W']\n",
    "colors = ['#976BE5','#E56BE5','#FF96FF', '#6EAF46']\n",
    "title = 'Weighted proportion of shared '+ id +' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p_SLNSF(dfplotSLNSF, labels, colors, title, subtitle, level, 30.1, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplotLSW, dfplot_normalized, level = calcperc_LSW(comm, separated, level)\n",
    "#variables\n",
    "labels = ['NSF samples', 'Large and whole','Small and whole', 'Large, small, and whole']\n",
    "colors = ['#8390FA', '#FF0000', '#DAD746', '#E89618']\n",
    "title = 'Weighted proportion of shared '+ id +' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p_SLNSF(dfplotLSW, labels, colors, title, subtitle, level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfplotLS_W, dfplot_normalized, level = calcperc_LS_W(comm, separated, level)\n",
    "#variables\n",
    "labels = ['NSF samples', 'Large and whole','Small and whole']\n",
    "colors = ['#8390FA', '#FF0000', '#DAD746']\n",
    "title = 'Weighted proportion of shared '+ id +' between size fractionated and non size fractionated samples'\n",
    "\n",
    "plot_stackedbar_p_SLNSF(dfplotLS_W, labels, colors, title, subtitle, level, 45.1, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newseparated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level = 'feature_id'\n",
    "\n",
    "toptaxa = separated[[level, 'feature_frequency', 'Taxon', 'size_code', 'weekn']].copy()\n",
    "    \n",
    "toptaxa = toptaxa.drop_duplicates()\n",
    "df_agg = toptaxa.groupby(['size_code',level]).agg({'feature_frequency':sum})\n",
    "    \n",
    "df_agg = df_agg.reset_index()\n",
    "resultpivot = df_agg.pivot_table(index=level, columns='size_code', values='feature_frequency')\n",
    "resultpivot = resultpivot.fillna(0)\n",
    "    \n",
    "df1 = resultpivot.copy()\n",
    "    \n",
    "df = resultpivot[['L', 'S', 'W']].copy()\n",
    "Sonly = df[(df['L'] == 0) & (df['W'] == 0)]\n",
    "Wonly = df[(df['L'] == 0) & (df['S'] == 0)]\n",
    "Lonly = df[(df['S'] == 0) & (df['W'] == 0)]\n",
    "LW = df[(df['S'] == 0) & (df['W'] != 0) & (df['L'] != 0)]\n",
    "LS = df[(df['W'] == 0) & (df['S'] != 0) & (df['L'] != 0)]\n",
    "SW = df[(df['W'] != 0) & (df['S'] != 0) & (df['L'] == 0)]\n",
    "LSW = df[~(df == 0).any(axis=1)]\n",
    "        \n",
    "c = venn3(subsets = (len(Lonly), len(Sonly), len(LS), len(Wonly), len(LW), len(SW), len(LSW)), \n",
    "          set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'),\n",
    "          #set_colors=('#E56BE5', '#976BE5', '#6EAF46'),\n",
    "          set_colors=('#E56BE5', '#976BE5', '#6EAF46'), alpha = 1);\n",
    "\n",
    "\n",
    "plt.savefig(\"outputs/02-EUKs/D\"+level+\"_vennall.png\", dpi=200, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** another idea is to run ancom of sizefraction specific and compare after the categories (run ancom on ?time or month.. or some other column) and compare the number/taxonomy of differentially abundant taxa recovered;\n",
    "are we recovering the same diff ab taxa between the (1.SF samples, 2. NSF samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca, pca_features, sfdclr, dm, plot_df2, df, distance_matrix = pcaplot(newseparated, 1, comm, 'Size code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 18S ANCOM PER DEPTH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca1_18, pca_features1_18, sfdclr1_18 = pcaplot(newseparated, 1, comm, 'Size code')\n",
    "DARejected_SC_taxonomy1_18, prcentile1_18 = run_ancom(newseparated, sfdclr1_18, 1, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prcentile1_18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy1_18.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca5_18, pca_features5_18, sfdclr5_18 = pcaplot(newseparated, 5, '18S', 'Size Code')\n",
    "DARejected_SC_taxonomy5_18, prcentile5_18 = run_ancom(newseparated, sfdclr5_18, 5, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy5_18.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca10_18, pca_features10_18, sfdclr10_18 = pcaplot(newseparated, 10, '18S', 'Size code')\n",
    "DARejected_SC_taxonomy10_18, prcentile10_18 = run_ancom(newseparated, sfdclr10_18, 10, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy10_18.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca30_18, pca_features30_18, sfdclr30_18 = pcaplot(newseparated, 30, '18S', 'Size code')\n",
    "DARejected_SC_taxonomy30_18, prcentile30_18 = run_ancom(newseparated, sfdclr30_18, 30, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy30_18.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca60_18, pca_features60_18, sfdclr60_18 = pcaplot(newseparated, 60, '18S', 'size code')\n",
    "DARejected_SC_taxonomy60_18, prcentile60_18 = run_ancom(newseparated, sfdclr60_18, 60, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy60_18.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16S ANCOM PER DEPTH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newseparated.loc[newseparated['feature_id'] == '51d96e96a3350beedece3878d6d0b3e7']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = newseparated[newseparated.size_code != 'S']\n",
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = df1[df1.size_code != 'L']\n",
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newseparated = df[df.size_code != 'S']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newseparated = newseparated[newseparated.depth == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newbiom = newseparated.pivot_table(index=\"feature_id\", columns=\"sampleid\", values=\"feature_frequency\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newbiom = newbiom.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newbiom.to_csv('newbiomdepth1.tsv', sep=\"\\t\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca1_16, pca_features1_16, sfdclr1_16 = pcaplot(df1, 1, '16S', 'Size code')\n",
    "DARejected_SC_taxonomy1_16, prcentile1_16 = run_ancom(df1, sfdclr1_16, 1, 'weekn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prcentile1_16 = prcentile1_16.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prcentile1_16.loc[prcentile1_16['feature_id'] == '80c73848b68ff95bd030fccfce011294']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy1_16.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca5_16, pca_features5_16, sfdclr5_16 = pcaplot(newseparated, 5, '16S', 'Size code')\n",
    "DARejected_SC_taxonomy5_16, prcentile5_16 = run_ancom(newseparated, sfdclr5_16, 5, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy5_16.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca10_16, pca_features10_16, sfdclr10_16 = pcaplot(newseparated, 10, '16S', 'Size code')\n",
    "DARejected_SC_taxonomy10_16, prcentile10_16 = run_ancom(newseparated, sfdclr10_16, 10, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy10_16.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca30_16, pca_features30_16, sfdclr30_16 = pcaplot(newseparated, 30, '16S', 'Size code')\n",
    "DARejected_SC_taxonomy30_16, prcentile30_16 = run_ancom(newseparated, sfdclr30_16, 30, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy30_16.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca60_16, pca_features60_16, sfdclr60_16 = pcaplot(newseparated, 60, '16S', 'Size code')\n",
    "DARejected_SC_taxonomy60_16, prcentile60__16 = run_ancom(newseparated, sfdclr60_16, 60, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DARejected_SC_taxonomy60_16.sort_values(by='W')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "? run a clustermap of top10 taxa of each deapth and color rows by depth, month, size code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prcentile60__16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca, pca_features, sfdclr = pcaplot(separated, 1, '16S')\n",
    "DARejected_month_taxonomy_16_1, prcentile = run_ancom(sfdclr, 1, 'Month')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_per_fid('16S', separated, 60, '50a0c3221c68046dfc96e032aff1ccd8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_df2.sort_values('dim1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Upset plot data prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upsetprep('16S', 'Genus')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if we want to make an upset plot of all depths?\n",
    "\n",
    "frames = sfd1, sfd5, sfd10, sfd30, sfd60\n",
    "result = pd.concat(frames)\n",
    "resultpivot = result.pivot_table(index='Genus', columns='set_name', values='feature_frequency')\n",
    "resultpivot = resultpivot.fillna(0)\n",
    "resultpivot[resultpivot != 0] = 1\n",
    "tosave = pd.merge(resultpivot, cumulab1, left_index=True, right_index=True)\n",
    "tosave.to_csv('genus_all16S_relab.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Venn diagrams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Depth 1 all 16S at genus level\n",
    "#N=84\n",
    "venn3(subsets = (27, 11, 4, 7, 6, 5, 22), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "#venn3(subsets = (Lonly, Sonly, LS, Wonly, LW, SW, LSW), set_labels = ('Large >3μm', 'Small 3-02μm', 'Whole water <0.22μm'), alpha = 0.5);\n",
    "\n",
    "plt.title(\"1m depth\") \n",
    "plt.savefig(\"outputs/02-EUKs/D1_genus_venn.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm = sfa.ols('nASVs ~ C(size_code)', data=sfd_LM).fit()\n",
    "anova = sa.stats.anova_lm(lm)\n",
    "spPH.posthoc_ttest(sfd_LM, val_col='nASVs', group_col='size_code', p_adjust='holm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd_S = sfd[['size_code', 'nASVs', 'weekn']].copy()\n",
    "sfd_S = sfd_S.drop_duplicates()\n",
    "sdfpv = sfd_S.pivot(index='weekn', columns='size_code', values='nASVs')\n",
    "fvalue, pvalue = stats.f_oneway(sdfpv['L'], sdfpv['S'], sdfpv['W'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd_LM = sfd[['size_code', 'nASVs']].copy()\n",
    "sfd_LM = sfd_LM.drop_duplicates()\n",
    "lm = sfa.ols('nASVs ~ C(size_code)', data=sfd_LM).fit()\n",
    "anova = sa.stats.anova_lm(lm)\n",
    "spPH.posthoc_ttest(sfd_LM, val_col='nASVs', group_col='size_code', p_adjust='holm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dffff = spPH.posthoc_ttest(sfd_LM, val_col='nASVs', group_col='size_code', p_adjust='holm')\n",
    "dffff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd_S = sfd_S.drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd_S = sfd_S.set_index('weekn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfd_sdf = sfd_S.stack().to_frame().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdfpv = sfd_S.pivot(index='weekn', columns='size_code', values='nASVs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdfpv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fvalue, pvalue = stats.f_oneway(sdfpv['L'], sdfpv['S'], sdfpv['W'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stackedbar_p(dfplot, labels, colors, title, subtitle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.savefig('outputs/'+comm_id+'/D'+str(depth)+'_adlineplot.png', dpi=200, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_grouped.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##make new category of S+L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make sure all size codes are indicated\n",
    "all_md[\"size_code\"] = all_md[\"sampleid\"].str.extract(r'[1-9][0-9]?[A-E]([L-S])')\n",
    "all_md[\"size_code\"] = all_md[\"size_code\"].fillna('W')\n",
    "\n",
    "#only keep values from weeks 1 to 16\n",
    "sep_SL = all_md[all_md.size_code != \"W\"]\n",
    "sep_SL = sep_SL.drop(sep_SL[sep_SL.weekn > 16].index)\n",
    "\n",
    "#sum [DNA] of small and large size fractions\n",
    "sep_SL['[DNAt]'] = sep_SL.groupby(['weekn', 'depth'])['[DNA]ng/ul'].transform('sum')\n",
    "\n",
    "#separate small and size fraction\n",
    "sep_S = sep_SL[sep_SL.size_code == 'S']\n",
    "sep_L = sep_SL[sep_SL.size_code == 'L']\n",
    "\n",
    "#calculate DNA proportion per size fraction\n",
    "sep_SL['DNApr'] = sep_SL['[DNA]ng/ul']/sep_SL['[DNAt]']\n",
    "\n",
    "#merge with separated on common columns to get corresponding rel. abundances\n",
    "sep_SL = sep_SL[['sampleid', 'DNApr', '[DNAt]']].copy()\n",
    "sepSLRA = pd.merge(separated, sep_SL, on=['sampleid'], how='left') #all_md is the metadata file\n",
    "\n",
    "#exclude ASVs from the whole water\n",
    "sep_SLRA = sepSLRA[separated.size_code != \"W\"]\n",
    "\n",
    "#calculate corrected per sample ratio, and corrected feature frequency of de-fractionated samples\n",
    "sep_SLRA['Newfeature_frequency'] = sep_SLRA['feature_frequency'] * sep_SLRA['DNApr']\n",
    "sep_SLRA['Newff'] = sep_SLRA.groupby(['feature_id', 'weekn', 'depth'])['Newfeature_frequency'].transform('sum')\n",
    "\n",
    "\n",
    "#sep_SLRA = sep_SLRA.drop(['sampleid', 'size_code'], axis=1)\n",
    "#sep_SLRA['sampleid'] = \"BB22.\" + sep_SLRA['weekn'].astype(str) + sep_SLRA['depth_code'] + \"SL\"\n",
    "\n",
    "#uncomment the line below if keeping small and large original sample\n",
    "#sep_SLRA['size_code'] = sep_SLRA['size_code'] + '-DFr'\n",
    "\n",
    "#uncomment the line above if merging smallandlarge\n",
    "#sep_SLRA['size_code'] = 'SL'\n",
    "\n",
    "sep_SLRA = sep_SLRA.drop(['feature_frequency', 'Total', 'ratio', 'nASVs', 'weekdepth', 'avg',\n",
    "                          'diff', 'extraction_date', '[DNA]ng/ul', 'A260/280', 'A260/230',\n",
    "                          'Newfeature_frequency'], axis=1)\n",
    "sep_SLRA.rename(columns={'Newff':'feature_frequency'}, inplace=True)\n",
    "sep_SLRA = sep_SLRA.drop_duplicates()\n",
    "\n",
    "sep_SLRA['Total'] = sep_SLRA['feature_frequency'].groupby(sep_SLRA['sampleid']).transform('sum')\n",
    "sep_SLRA['ratio'] = sep_SLRA['feature_frequency']/sep_SLRA['Total']\n",
    "sep_SLRA['nASVs'] = sep_SLRA['feature_id'].groupby(sep_SLRA['sampleid']).transform('nunique')\n",
    "\n",
    "sep_SLRA = sep_SLRA.drop_duplicates()\n",
    "\n",
    "sep_WO = separated[separated.size_code == \"W\"]\n",
    "sep_L = separated[separated.size_code == \"L\"]\n",
    "sep_S = separated[separated.size_code == \"S\"]\n",
    "\n",
    "newseparated = pd.concat([sep_SLRA, sep_WO], ignore_index=True)\n",
    "\n",
    "newseparated['weekdepth'] = newseparated[\"weekn\"].astype(str) + newseparated[\"depth\"].astype(str)\n",
    "newseparated['avg'] = newseparated['nASVs'].groupby(newseparated['weekdepth']).transform('mean')\n",
    "newseparated['diff'] = newseparated['nASVs'] - newseparated['avg']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Permanova results from R into boxplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "permresu = pd.read_csv('R_results/post_hoc_results.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "permresu[\"depth_pairs\"] = permresu[\"depth\"].astype(str) + permresu[\"pairs\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.catplot(\n",
    "    permresu, kind=\"bar\",\n",
    "    x=\"p.adjusted\", y=\"pairs\", col=\"comm\", hue=\"depth\",\n",
    "    height=4, aspect=1.3, palette=\"Greys\", log=True\n",
    ")\n",
    "#ax.set(xlim=(0, 0.10))\n",
    "\n",
    "ax.refline(x=0.05, color='red')\n",
    "\n",
    "plt.savefig('outputs/perm_pvalues_logged.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.catplot(\n",
    "    permresu, kind=\"bar\",\n",
    "    x=\"p.adjusted\", y=\"pairs\", col=\"comm\", hue=\"depth\",\n",
    "    height=4, aspect=1.3, palette=\"Greys\", log=True\n",
    ")\n",
    "#ax.set(xlim=(0, 0.10))\n",
    "\n",
    "ax.refline(x=0.05, color='red')\n",
    "\n",
    "plt.savefig('outputs/perm_pvalues_logged.png', bbox_inches='tight', dpi=300)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qiime2-2023.5",
   "language": "python",
   "name": "qiime2-2023.5"
  },
  "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}