bb22_sf_microbial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cd872fb8",
   "metadata": {},
   "source": [
    "## Upload functions and necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f12e05f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bb2022_functions import *\n",
    "%matplotlib inline\n",
    "from Bio.SeqIO.FastaIO import SimpleFastaParser\n",
    "from Bio import SeqIO\n",
    "pd.options.mode.chained_assignment = None  # default='warn'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef8c5429",
   "metadata": {},
   "source": [
    "## Import and format metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fb6bbc2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "md = pd.read_csv(\"metadata_merged.csv\")\n",
    "merged = pd.read_csv(\"metadata_niskin.csv\")\n",
    "all_md = pd.read_csv(\"allmetadata.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "402ff79a",
   "metadata": {},
   "source": [
    "### Visualize metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da41acb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "maxvals = plot_nutrients(all_md, 'euphotic')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc70d5e",
   "metadata": {},
   "source": [
    "## Add microbial communities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "479815e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Community is 02-PROKs\n",
      "Found all 16S tables.\n",
      "Successfully saved all tables.\n",
      "Set up metadata ...\n",
      "Saved merged_asvs_metadata.tsv\n",
      "Appended all taxonomies to taxos\n",
      "Saved separated by metadata dataframe.\n",
      "Community is 02-PROKs\n",
      "Removed cyanobacteria and chloroplast from 02-PROKs\n"
     ]
    }
   ],
   "source": [
    "#generate a dataframe from all specified amplicon\n",
    "df, comm = consolidate_tables('16S') #frac='pooled') #16S, chloroplast, or 18S\n",
    "merged = merge_metadata(df, all_md)\n",
    "separated, contaminants = pick_metadata(comm, merged)\n",
    "newseparated = make_defract(all_md, separated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "64ad45e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Values were updated\n"
     ]
    }
   ],
   "source": [
    "#apply changes to taxonomy according to NCBI identified ASVs\n",
    "newdf = apply_replacement(newseparated, \"feature_id\", \"Genus\") \n",
    "# or replace Genus with PRSpecies if dealing with phytoref"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "631b6453",
   "metadata": {},
   "source": [
    "To generate MIMARKS file for NCBI sequence submission; the output is a .csv file for the samples and their metadata for submission (i.e sampleid, size fraction, date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5782332e",
   "metadata": {},
   "outputs": [],
   "source": [
    "make_MIMARKS(newseparated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82cec6f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot of dna concentrations per sample\n",
    "dnacon(newseparated, depth='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3738433b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#rarefaction curves per community\n",
    "rarefy_curve(comm, newseparated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a14289f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot sampling depth histogram to pick library depth for rarefying\n",
    "plt.hist(newdf['Total'], bins=10, edgecolor='black')\n",
    "\n",
    "# Add lines for quantiles\n",
    "quantiles = [0.1, 0.2, 0.5]\n",
    "colors = ['red', 'blue', 'green']\n",
    "labels = ['10th Percentile', '20th Percentile', '50th Percentile (Median)']\n",
    "\n",
    "for q, color, label in zip(quantiles, colors, labels):\n",
    "    plt.axvline(x=newdf['Total'].quantile(q), color=color, linestyle='--', label=label)\n",
    "plt.legend()\n",
    "\n",
    "# Save and show the plot\n",
    "plt.savefig('sampling_depths.png', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "# Suggest a rarefaction depth\n",
    "suggested_depth = int(newdf['Total'].quantile(0.5))\n",
    "print(f\"Suggested rarefaction depth: {suggested_depth}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d5a8df5",
   "metadata": {},
   "source": [
    "We can also sort the samples by their library size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75273c16",
   "metadata": {},
   "outputs": [],
   "source": [
    "newseparated[['Total','sampleid']].sort_values('Total').drop_duplicates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebab5e55",
   "metadata": {},
   "source": [
    "#### Explore the taxonomy in the samples and compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "945f9da7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Produce interactive taxonomic barplots with plotly\n",
    "phyld, top10d = taxbarplot(comm, newseparated, 'Genus', 5, 10, 'size_code')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5be826b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize the static barplots with seaborn, and each size fraction separately\n",
    "taxonomic_barplots(comm, newdf, [5,60], 'Genus', 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c534508d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate the heatmap for the top genus from each sample\n",
    "heatmap_top1(comm, newseparated, 'Genus')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baae7dca",
   "metadata": {},
   "source": [
    "The above plot uses taxonomy, but we can generate the same plot but by comparing whether 80% of the features in each samples are also found in the whole (unfractionated samples). This was quantified by dividing the number of shared features by the toal number of features. If a square has a red color (closer to 1), it's very similar to the unfractionated sample, and the bluer the square, the more different it is from the unfractionated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8afbafc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "grab_80perc(comm, newseparated, 0.8, 'feature_id')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f2a2e05",
   "metadata": {},
   "source": [
    "We can plot alpha diversity measurements, whether as 'shannon_diversity' or 'nASVs' which is the richness quantified by the total number of ASVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "155a887c",
   "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(newdf, comm, 'all', 'shannon_diversity', 'Shannon Diversity Index')\n",
    "#results gives the corrected p-values for pairwise comparisons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ac93b5",
   "metadata": {},
   "source": [
    "Compare the slopes of linear regressions of the richness change over time. Each value represents how much a size fraction (column) differs in comparison to the average slope (averaged between all size fractions) for each depth (rows)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4d80b53",
   "metadata": {},
   "outputs": [],
   "source": [
    "tohm, z_sc_df = get_slopes(comm, separated)\n",
    "#a zscore of 1= 1 std away from the mean,\n",
    "#positive values=higher than mean, neg= smaller than mean"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efa8fadd",
   "metadata": {},
   "source": [
    "### Beta diversity and ANCOM analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d19b1a82",
   "metadata": {},
   "source": [
    "Optionally we can run ANCOM with removed low abundance features with a given threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc00b39d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#only if we want to run ANCOM pairwise\n",
    "news2 = newseparated[newseparated.size_code != 'L']\n",
    "news2 = news2[news2.size_code != 'SL']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70cabc53",
   "metadata": {},
   "outputs": [],
   "source": [
    "depths = [1,5,10,30,60]\n",
    "for depth in depths:\n",
    "    pca, pca_features, sfdclr, dm = pcaplot(newseparated, depth, comm, 'size_code', 'DFr', 'week')\n",
    "    DAresults, DARejected_SC_taxonomy, prcentile = run_ancom(comm, newseparated, sfdclr, depth, 'size_code', threshold=0)\n",
    "\n",
    "    #save outputs\n",
    "    DAresults.to_csv('outputs/ANCOM/chloroplast/none/'+comm+'_D'+str(depth)+'_WSLSL.csv')\n",
    "    DARejected_SC_taxonomy.to_csv('outputs/ANCOM/chloroplast/none/'+comm+'_D'+str(depth)+'_Trueonly_WSLSL.csv')\n",
    "\n",
    "    notify()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d412eea7",
   "metadata": {},
   "source": [
    "Depending on ancom results, we can investigate single features temporal dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd391fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_id = '5a94578dd1d7cdd039a52f1c7079f874'\n",
    "newdf.loc[newdf['feature_id'] == f_id, 'Taxon'].tolist()[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0abfc9",
   "metadata": {},
   "source": [
    "Visualize the time series of a single feature in each size fraction over the 16 weeks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e5bfdc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "timeseries_fid(comm, newseparated, f_id, 'Oligotrichia', 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56cc9acd",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "feature_id_summary = count_feature_id_presence_with_depth_and_W('outputs', comm)\n",
    "top_asvs_summary = filter_top_asvs(feature_id_summary, method=\"top_W_sum\", n=50)\n",
    "plot_asv_heatmap(comm, feature_id_summary, file_filter=\"WSLSL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb8900a1",
   "metadata": {},
   "source": [
    "### Export files to R analyses"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8bdf7ac",
   "metadata": {},
   "source": [
    "To create a phyloseq object, you need an ASV table, taxonomy file and metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fef52ea",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "238fc456",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffd74b7b",
   "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": 5
}