family_wise_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm_notebook\n",
    "tqdm_notebook().pandas()\n",
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "import numpy as np\n",
    "from scipy.stats import ks_2samp\n",
    "from matplotlib.colors import LogNorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = \"somatosensory_converted_into_tpm\"\n",
    "families = \"clean_panther4march\"\n",
    "\n",
    "df = pd.read_csv(\"{}/{}/results/family_IC.csv\".format(dataset, families), index_col=\"family_id\")\n",
    "shuffled_df = pd.read_csv(\"{}/{}/results/shuffled_family_IC.csv\".format(dataset, families), index_col=\"family_id\")\n",
    "raw_df = pd.read_csv(\"{}.csv\".format(dataset), index_col=\"gene_id\")\n",
    "family_df = pd.read_csv(\"{}.csv\".format(families))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.family_size.hist(bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = df.family_size.value_counts()\n",
    "family_sizes = counts[counts >= 30].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summaries = []\n",
    "\n",
    "for family_size in family_sizes:\n",
    "    display(\"Family size: \"+str(family_size))\n",
    "    orig = df.loc[df.family_size == family_size]\n",
    "    shuffled = shuffled_df.loc[shuffled_df.family_size == family_size]\n",
    "\n",
    "    #KS-test\n",
    "    display(ks_2samp(orig.ic, shuffled.ic))\n",
    "    sns.distplot(orig.ic.dropna(), norm_hist=True, rug=False, kde=False, label=\"Original\")\n",
    "    sns.distplot(shuffled.ic.dropna(), norm_hist=True, rug=False, kde=False, label=\"Shuffled\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    orig_ecdf = ECDF(orig.ic)\n",
    "    shuffled_ecdf = ECDF(shuffled.ic)\n",
    "    x = np.arange(0,5,0.01)\n",
    "    sns.lineplot(x, orig_ecdf(x), label=\"Original\")\n",
    "    sns.lineplot(x, shuffled_ecdf(x), label=\"Shuffled\")\n",
    "    sns.lineplot(x, shuffled_ecdf(x) - orig_ecdf(x), label=\"Difference\")\n",
    "    plt.show()\n",
    "\n",
    "    #Permutation test\n",
    "    def get_statistics(df):\n",
    "        ics = df.ic\n",
    "        return pd.Series({\n",
    "            \"median\": ics.median(),\n",
    "            \"percentile_90\": ics.quantile(0.9),\n",
    "            \"percentile_10\": ics.quantile(0.1),\n",
    "            \"percentile_97.5\": ics.quantile(0.975),\n",
    "            \"percentile_02.5\": ics.quantile(0.025),\n",
    "            \"quantile_ratio\": ics.quantile(0.9) / ics.quantile(0.1),\n",
    "#             \"skew\": ics.skew()\n",
    "        })\n",
    "    orig_statistics = get_statistics(orig)\n",
    "    shuffled_statistics = shuffled.groupby(\"iteration\").progress_apply(get_statistics)\n",
    "#     total_shuffled_statistic = get_statistics(shuffled) #WRONG, this is not median\n",
    "    total_shuffled_statistic = shuffled_statistics.median()\n",
    "\n",
    "#     lower_count = (shuffled_statistics <= orig_statistics).sum()\n",
    "#     upper_count = (shuffled_statistics >= orig_statistics).sum()\n",
    "#     lower_pvalue = (lower_count + 1) / (shuffled_statistics.shape[0] + 1)\n",
    "#     upper_pvalue = (upper_count + 1) / (shuffled_statistics.shape[0] + 1)\n",
    "    \n",
    "    shuf_mean = shuffled_statistics.mean(axis=0)\n",
    "    orig_diff = np.abs(orig_statistics - shuf_mean)\n",
    "    shuf_diff = shuffled_statistics.subtract(shuf_mean).abs()\n",
    "    pvalue = ((shuf_diff >= orig_diff).sum(axis=0) + 1) / (shuffled_statistics.shape[0] + 1)\n",
    "    pvalues = pd.DataFrame({\n",
    "        \"orig_value\": orig_statistics,\n",
    "        \"shuffled_value\": total_shuffled_statistic,\n",
    "#         \"lower_count\": lower_count,\n",
    "#         \"lower_pvalue\": lower_pvalue,\n",
    "#         \"upper_count\": upper_count,\n",
    "#         \"upper_pvalue\": upper_pvalue,\n",
    "        \"pvalue\": pvalue,\n",
    "    })\n",
    "#     pvalues[\"significance\"] = pvalues.apply(lambda x: \"LOWER\" if x.lower_pvalue <= 0.025 else (\"HIGHER\" if x.upper_pvalue <= 0.025 else \"-----\"), axis=1)\n",
    "\n",
    "    display(pvalues)\n",
    "    _, axs = plt.subplots(3,3,figsize=(15,12))\n",
    "    for ax, statistic in zip(np.array(axs).flatten(), orig_statistics.index):\n",
    "        sns.distplot(shuffled_statistics[statistic], ax=ax, kde=False, rug=True, label=\"Shuffled\")\n",
    "        sns.distplot([orig_statistics[statistic]], ax=ax, kde=False, hist=False, rug=True, rug_kws={\"height\": 0.5}, label=\"Original\")\n",
    "        sns.distplot([total_shuffled_statistic[statistic]], ax=ax, kde=False, hist=False, rug=True, rug_kws={\"height\": 0.5}, label=\"Combined Shuffled\")\n",
    "        ax.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    summaries += [{\n",
    "        \"family_size\": family_size,\n",
    "        \"ks\": ks_2samp(orig.ic, shuffled.ic),\n",
    "        \"pvalues\": pvalues\n",
    "    }]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = summaries[0]\n",
    "def build_row(x):\n",
    "    l = [x[\"family_size\"], df.loc[df.family_size == x[\"family_size\"]].shape[0], x[\"ks\"].statistic, x[\"ks\"].pvalue]\n",
    "    l += list(x[\"pvalues\"].loc[\"quantile_ratio\", [\"orig_value\", \"shuffled_value\"]].values)\n",
    "    l += [x[\"pvalues\"].loc[\"quantile_ratio\", \"pvalue\"]]\n",
    "    l += list(x[\"pvalues\"].loc[\"percentile_90\", [\"orig_value\", \"shuffled_value\"]].values)\n",
    "    l += [x[\"pvalues\"].loc[\"percentile_90\", \"pvalue\"]]\n",
    "    l += list(x[\"pvalues\"].loc[\"percentile_10\", [\"orig_value\", \"shuffled_value\"]].values)\n",
    "    l += [x[\"pvalues\"].loc[\"percentile_10\", \"pvalue\"]]\n",
    "    l += list(x[\"pvalues\"].loc[\"median\", [\"orig_value\", \"shuffled_value\"]].values)\n",
    "    l += [x[\"pvalues\"].loc[\"median\", \"pvalue\"]]\n",
    "    return pd.Series(l, index=[\n",
    "        [\"\", \"\", \"KS\", \"KS\"] + ([\"Quantile Ratio\"] * 3) + ([\"90%\"] * 3) + ([\"10%\"] * 3) + ([\"Median\"] * 3),\n",
    "        [\"Size\", \"No. Families\", \"Statistic\", \"P-value\"] + ([\"Original\", \"Shuffled\", \"P-value\"] * 4)\n",
    "    ])\n",
    "    \n",
    "summary_table = pd.DataFrame([build_row(x) for x in summaries])#.set_index((\"\", \"Family Size\"))\n",
    "# summary_table.index.name = \"Size\"\n",
    "summary_table.loc[:,(\"\",\"Size\")] = summary_table.loc[:,(\"\",\"Size\")].astype(int)\n",
    "summary_table.loc[:,(\"\",\"No. Families\")] = summary_table.loc[:,(\"\",\"No. Families\")].astype(int)\n",
    "summary_table.to_csv(\"{}/family_stat_test_pvalues.csv\".format(dataset))\n",
    "summary_table\n",
    "# print(summary_table.to_latex(float_format=\"%.3f\", multicolumn_format=\"c\", column_format=\"|l|rr|rrr|rrr|rrr|rrr|\", index=False))\n",
    "# np.array([0, 1, 2, list(x[\"pvalues\"].loc[\"median\", [\"orig_value\", \"shuffled_value\", \"lower_pvalue\", \"upper_pvalue\"]].values)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale=1.2, style=\"ticks\", font=\"Arial\")\n",
    "shuffled_statistics = shuffled_df.groupby(\"family_size\").apply(get_statistics)\n",
    "plt.figure(figsize=(8,2.8))\n",
    "sns.scatterplot(data=df, x=\"family_size\", y=\"ic\", label=\"Original\")\n",
    "sns.lineplot(shuffled_statistics.index, shuffled_statistics[\"median\"], label=\"Shuffled (Median)\")\n",
    "sns.lineplot(shuffled_statistics.index, shuffled_statistics[\"percentile_97.5\"], label=\"Shuffled (97.5th Percentile)\")\n",
    "# sns.lineplot(shuffled_statistics.index, shuffled_statistics.percentile_90, label=\"Shuffled 90%\")\n",
    "# sns.lineplot(shuffled_statistics.index, shuffled_statistics.percentile_10, label=\"Shuffled 10%\")\n",
    "sns.lineplot(shuffled_statistics.index, shuffled_statistics[\"percentile_02.5\"], label=\"Shuffled (2.5th Percentile)\",)\n",
    "plt.ylabel(\"IC\")\n",
    "plt.xlabel(\"Family Size\")\n",
    "plt.yscale(\"log\",basey=2)\n",
    "plt.xscale(\"log\")\n",
    "plt.legend(loc=\"upper left\", bbox_to_anchor=(1, 1))\n",
    "plt.show()\n",
    "# shuffled_statistics"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.4 64-bit ('stochastic': conda)",
   "language": "python",
   "name": "python37464bitstochasticcondad3c784986b534aaebb6b91784b827746"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}