C_Elegans_Lab-Fall2024.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C Elegans Analysis notebook\n",
    "This Jupyter notebook was written to aid you in a statistical analysis of the morphology of DD/VD motor neurons from the four C elegans genotypes. This notebook will load in the morphology data (.cvs file) collected from the whole class that you contribute to in lab this week. \n",
    "\n",
    "Jupyter is an interactive, web-based python interface. Jupyter is a great tool for data-science, since you can iteratively develop your project and debug your code without having to continually reload data into memory. To run each cell of the jupyter notebook, use  𝚌𝚝𝚛𝚕+𝚎𝚗𝚝𝚎𝚛 . In addition to jupyter, three python packages are used in this notebook: numpy (to perform computations + sampling), matplotlib (to plot the data, pandas (data format and statistic methods), and scipy (statistics toolbox).\n",
    "\n",
    "#### Learning Objectives\n",
    "- Use descriptive statistics to quantify the basic features of the data \n",
    "- Use an equivalent to a t-test (Mann-Whitney U test) to compare morphology between wildtype and the various mutant genotypes.\n",
    "- Interprete what the p-value means with respect to the tested hypothesis \n",
    "- Interprete the obtained p-values in the context of correction for mulitple comparisons \n",
    "\n",
    "\n",
    "#### Step 0: Load the packages to be used later in the analysis\n",
    "(Don't edit this cell, just run it using 𝚌𝚝𝚛𝚕+𝚎𝚗𝚝𝚎𝚛 and then move to the next cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "#### Step 1: Load in the data from your classmates and calculate precent defective for each genotype\n",
    "(No need to edit this cell, run it using 𝚌𝚝𝚛𝚕+𝚎𝚗𝚝𝚎𝚛 and examine the loaded data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load in the class data from the .csv file\n",
    "class_data = pd.read_csv(\"MCB_160L_F21_Lab_11_class_data.csv\",\n",
    "                         skiprows=1,\n",
    "                         header=[0, 1]).dropna()\n",
    "# list of genes\n",
    "genes = [\"juIs76\", \"juIs76; unc-40\", \"juIs76; unc-6\", \"juIs76; unc-5\"]\n",
    "\n",
    "# loop over each gene and extract the data section of class_data from that gene\n",
    "for gene in genes:\n",
    "    # for each gene, calculate the percent of defective worms\n",
    "    pct_defect = 100*class_data.loc[:, (gene, \"Abnormal DD/VD\")].values/class_data.loc[:, gene].values.sum(axis=1)\n",
    "    \n",
    "    # create new field called \"Percent Defective\" within class_data to be used for the analysis\n",
    "    class_data[(gene, \"Percent Defective\")] = pct_defect    \n",
    "class_data.sort_index(axis=1, inplace=True)\n",
    "\n",
    "#print the data\n",
    "print(\"DD/VD class morphology data:\")\n",
    "print(class_data, \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Calculate the mean % pathfinding defects for each genotype.\n",
    "\n",
    "Calculate the **mean** of \"Percent Defective\" for each of the four genotypes. To show an example the <code>mean</code> method has already been applied to the column \"Percent Defective\" for the first three genotypes.  Note that the <code>%.3f</code> in the print line tells Python to include the listed number (floating point number) rounded to 3 decimal points.\n",
    "\n",
    "You should extend this code using the same **mean** method to calculate the mean for the final \"juIs76; unc-5\" mutant genotype and print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate juIs76 mean of percent defective\n",
    "juIs76_percent = class_data[\"juIs76\"][\"Percent Defective\"]\n",
    "juIs76_mean = juIs76_percent.mean()\n",
    "print(\"mean juIs76 pct pathfinding defects: %.3f \" % juIs76_mean)\n",
    "\n",
    "# Calculate juIs76; unc-40 mean of percent defective\n",
    "juIs76_unc40_percent = class_data[\"juIs76; unc-40\"][\"Percent Defective\"]\n",
    "juIs76_unc40_mean = juIs76_unc40_percent.mean()\n",
    "print(\"mean juIs76; unc-40 pct pathfinding defects: %.3f \" % juIs76_unc40_mean)\n",
    "\n",
    "# Calculate juIs76; unc-6 mean of percent defective \n",
    "juIs76_unc6_percent = class_data[\"juIs76; unc-6\"][\"Percent Defective\"]\n",
    "juIs76_unc6_mean = juIs76_unc6_percent.mean()\n",
    "print(\"mean juIs76; unc-6 pct pathfinding defects: %.3f \" % juIs76_unc6_mean)\n",
    "\n",
    "# Calculate juIs76; unc-5 mean of percent defective\n",
    "# *add code here that calculates juIs76_unc5_mean*\n",
    "\n",
    "\n",
    "juIs76_unc5_mean = ...\n",
    "print(\"mean juIs76; unc-5 pct pathfinding defects: %.3f \" % juIs76_unc5_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Calculate the standard deviation of % pathfinding defects for each genotype.\n",
    "\n",
    "Calculate the **standard devation** of \"Percent Defective\" for each of the four genotypes. To show an example the <code>std</code> method has already been applied to the column \"Percent Defective\" for the first three genotypes. Note that the <code>%.3f</code> in the print line tells Python to include the listed number (floating point number) rounded to 3 decimal points.\n",
    "\n",
    "You should extend this code using the same **std** method to calculate the mean for the final \"juIs76; unc-5\" mutant genotype and print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate juIs76 standard deviation of percent defective\n",
    "juIs76_percent = class_data[\"juIs76\"][\"Percent Defective\"]\n",
    "juIs76_std = juIs76_percent.std()\n",
    "print(\"std juIs76 pct pathfinding defects: %.3f \" % juIs76_std)\n",
    "\n",
    "# Calculate juIs76; unc-40 std of percent defective\n",
    "juIs76_unc40_percent = class_data[\"juIs76; unc-40\"][\"Percent Defective\"]\n",
    "juIs76_unc40_std = juIs76_unc40_percent.std()\n",
    "print(\"std juIs76 unc-40 pct pathfinding defects: %.3f \" % juIs76_unc40_std)\n",
    "\n",
    "# Calculate juIs76; unc-6 std of percent defective \n",
    "juIs76_unc6_percent = class_data[\"juIs76; unc-6\"][\"Percent Defective\"]\n",
    "juIs76_unc6_std = juIs76_unc6_percent.std()\n",
    "print(\"std juIs76 unc-6 pct pathfinding defects: %.3f \" % juIs76_unc6_std)\n",
    "\n",
    "# Calculate juIs76; unc-5 std of percent defective\n",
    "# *add code here that calculates juIs76_unc5_std*\n",
    "\n",
    "juIs76_unc5_std = ...\n",
    "print(\"std juIs76; unc-5 pct pathfinding defects: %.3f \" % juIs76_unc5_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "#### Step 4: Make a bar graph of the means with error bars (include this in your lab report)\n",
    "(No need to edit this cell. Run it using 𝚌𝚝𝚛𝚕+𝚎𝚗𝚝𝚎𝚛, examine the plotted data and include this plot in your lab report)\n",
    "\n",
    "Note: this code uses variables (e.g. juIs76_unc5_mean and juIs76_unc5_std) that you will calculate above so make sure to complete Step 3 before running Step 4. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "# plot bar of the average percent defective for the juIs\n",
    "ax.bar(0, juIs76_mean)\n",
    "# plot errorbar for the standard deviation of percent defective for the gene\n",
    "ax.errorbar(0, juIs76_mean ,juIs76_std,capsize=5)\n",
    "\n",
    "# plot bar of the average percent defective for the juIs unc-40\n",
    "ax.bar(1, juIs76_unc40_mean)\n",
    "# plot errorbar for the standard deviation of percent defective for the gene\n",
    "ax.errorbar(1, juIs76_unc40_mean ,juIs76_unc40_std,capsize=5)\n",
    "\n",
    "# plot bar of the average percent defective for the juIs unc-6\n",
    "ax.bar(2, juIs76_unc6_mean)\n",
    "# plot errorbar for the standard deviation of percent defective for the gene\n",
    "ax.errorbar(2, juIs76_unc6_mean ,juIs76_unc6_std,capsize=5)\n",
    "\n",
    "# plot bar of the average percent defective for the juIs unc-5\n",
    "ax.bar(3, juIs76_unc5_mean)\n",
    "# plot errorbar for the standard deviation of percent defective for the gene\n",
    "ax.errorbar(3, juIs76_unc5_mean ,juIs76_unc5_std,capsize=5)\n",
    "\n",
    "# format and label bar plot\n",
    "ax.set_ylim(bottom=0)\n",
    "ax.set_xticks(range(len(genes))) # list of genes \"juIs76\", \"juIs76; unc-40\", \"juIs76; unc-6\", \"juIs76; unc-5\"\n",
    "ax.set_xticklabels(genes)\n",
    "ax.set_title(\"% Pathfinding Defects by Gene\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "source": [
    "#### Step 5: Run a Mann-Whitney U two-sided rank test to compare the % pathfinding effects between wildtype and each mutant genotype.\n",
    "\n",
    "The Mann-Whitney U test is a nonparametric equivalent test to a t-test.  However, unlike the t-test is does not require the assumption that the distributions follows a normal distributions. The Mann-Whitney U test tests the null hypothesis that the distribution underlying sample x is the same as the distribution underlying sample y. Therefore lower p-values indicate that the two samples are likely not from the same distribution. \n",
    "\n",
    "(No need to edit this cell. Run it using 𝚌𝚝𝚛𝚕+𝚎𝚗𝚝𝚎𝚛, examine and interprete the resulting p-values for your analysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# wildtype genes and mutant gene lists\n",
    "wildtype_gene = \"juIs76\"\n",
    "mutant_genes = [\"juIs76; unc-40\", \"juIs76; unc-6\", \"juIs76; unc-5\"]\n",
    "\n",
    "# loop over mutant list and compare each to wildtype\n",
    "for curr_mutant_gene in mutant_genes:\n",
    "    # for each mutant genes, compute a mannwhitneyu 2-tailed t-test, assuming pooled variance.\n",
    "    test_results = st.mannwhitneyu(x=class_data[(wildtype_gene, \"Percent Defective\")],\n",
    "                        y=class_data[(curr_mutant_gene, \"Percent Defective\")],\n",
    "                        alternative='two-sided')\n",
    "    \n",
    "    print(\"Comparison of % Pathfinding Defects in [\" + wildtype_gene + \"] vs [\" + curr_mutant_gene + \"]\")\n",
    "    print(\"\\t t-statistic: %.3f\" % test_results[0])\n",
    "    print(\"\\t p-value: %.3e\\n\" % test_results[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6: Calculate the appropriate cut-off (alpha) value for signficance of your Mann-Whitney analysis\n",
    "In statistics when you perform muiltiple statisitcal tests looking for a positive result it is important to correct for multiple comparisons. The intuition behind this issue is that as you perform more and more tests the chances increase that a false positive result will occur just by chance. \n",
    "\n",
    "To adjust for the mulitple comparisons we will use the Bonferroni correction method. Bonferroni correction can be used when several statistical tests are being performed simultaneously (while a given alpha (e.g. 0.05) value may be appropriate for each individual comparison, it is not for the set of all comparisons). The way this correction work is the alpha value needs to be lowered to account for the number of comparisons being performed by dividing the original alpha value (typically 0.05) by the number of test being performed.  By lowering the alpha value less signficant results are found thus avoiding a higher rate of  false positives. \n",
    "\n",
    "More details about Bonerroni correction:\n",
    "https://mathworld.wolfram.com/BonferroniCorrection.html\n",
    "\n",
    "To be concrete: since we performed 3 different statstical tests (one for each mutant genotype compared to wildtype) the alpha should be divided by 3. \n",
    "Hint: In python the \"/\" operator allows for division. \n",
    "\n",
    "Below finish this segement of code to calculate your adjusted alpha so you can use it for interpreting your results in you lab assignment and lab report. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_tests = len(mutant_genes)\n",
    "alpha = 0.05\n",
    "\n",
    "# *add code here that calculates adjusted_alpha*\n",
    "adjusted_alpha = ...\n",
    "\n",
    "print(\"The adjusted cut-off for signficance (alpha) is: %.3f\" % adjusted_alpha )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}