4_permutation_testing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/quanti/Desktop/M_HSC/\n"
     ]
    }
   ],
   "source": [
    "import itertools as itt\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import scipy.optimize as spo\n",
    "import scipy.special as spsp\n",
    "import scipy.stats as sps\n",
    "import seaborn as sns\n",
    "\n",
    "from matplotlib.ticker import FixedLocator\n",
    "\n",
    "mpl.rcParams['axes.titlesize'] = 'xx-large'\n",
    "mpl.rcParams['axes.labelsize'] = 'xx-large'\n",
    "mpl.rcParams['xtick.labelsize'] = 'x-large'\n",
    "mpl.rcParams['ytick.labelsize'] = 'x-large'\n",
    "mpl.rcParams['xtick.direction'] = 'out'\n",
    "mpl.rcParams['ytick.direction'] = 'out'\n",
    "mpl.rcParams['legend.frameon'] = True\n",
    "mpl.rcParams['legend.framealpha'] = 0.5\n",
    "mpl.rcParams['legend.fontsize'] = 'large'\n",
    "\n",
    "path_proj = os.getcwd()+ '/'\n",
    "path_sc = path_proj+'csv/Single_cell/'\n",
    "path_sort = path_proj+'csv/Sort/'\n",
    "print(path_proj)\n",
    "\n",
    "#LOAD THE PROCESSED DATA\n",
    "df = pd.read_csv(path_sc+'Pooled_data.csv', sep=';', decimal=',')\n",
    "\n",
    "#GLOBAL VARIABLES FOR ITERATION\n",
    "or_cells = ['SLAM-HSC', 'ST-HSC', 'MPP']\n",
    "conds = ['P', 'P+ILs']\n",
    "times = ['24h', '48h']\n",
    "\n",
    "#SET A RANDOM SEED TO REPLICATE THE P-VALUES\n",
    "#FROM STATISTICAL TESTS BASED ON RANDOM PERMUTATIONS\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#MAP THE VECTORS OF counts OF classes TO A DIFFERENT SUPPORT (USUALLY CONTAINING classes)\n",
    "def spread_to_support(classes, counts, support):\n",
    "    counts_sum = counts.sum()\n",
    "    if counts_sum == 0:\n",
    "        return np.zeros(len(support))\n",
    "    else:\n",
    "        return np.hstack([counts[classes==cl] if any(classes==cl) else np.zeros(1) for cl in support])\n",
    "\n",
    "#CALCULATE THE G-TEST STATISTIC ON data WITH labels,\n",
    "#ON THE CONTINGENCY TABLE OF support (DATA SUPPORT) X all_perm_labels (LABELS SUPPORT)\n",
    "def chi_squared_stat(data, labels, support, all_perm_labels, bool_stack=True):\n",
    "    if bool_stack:\n",
    "        table = np.array([\n",
    "            spread_to_support(*np.unique(np.hstack(data[labels==k]), return_counts=True),\n",
    "                              support=support)\n",
    "            for k in all_perm_labels\n",
    "        ])\n",
    "    else:\n",
    "        table = np.array([\n",
    "            spread_to_support(*np.unique(data[labels==k], return_counts=True),\n",
    "                              support=support)\n",
    "            for k in all_perm_labels\n",
    "        ])\n",
    "    return sps.chi2_contingency(table, correction=False, lambda_='log-likelihood')[0]\n",
    "\n",
    "####PERMUTATION TESTING\n",
    "#INPUT\n",
    "#    stat_func: the statistic to be calculated on data and permuted data\n",
    "#    args: fixed input variables for stat_func\n",
    "#    data: array of data to be permuted \n",
    "#    perm_label: array of labels, with same lenght as data, used to permute data\n",
    "#    n_iter: number of permuted data to be drawn\n",
    "#    boold_p_dist: if True, also return the stat_func calculated on the permuted data.\n",
    "#        Useful to plot the statistic' distribution. \n",
    "#OUTPUT: left, right and two-sided pvalues;\n",
    "#        if boold_p_dist=True, also return the permuted versions of the statistics\n",
    "\n",
    "def permutation_testing(stat_func, args, data, perm_label, n_iter=250000, boold_p_dist=False):\n",
    "    stat = stat_func(data, perm_label, *args)\n",
    "    perm_stat = np.array([\n",
    "        stat_func(data, perm_label[np.random.permutation(len(perm_label))], *args)\n",
    "        for k in np.arange(n_iter)])\n",
    "    \n",
    "    left_pval = ((perm_stat <= stat).sum()+1.)/(n_iter+1.)\n",
    "    right_pval = ((perm_stat >= stat).sum()+1.)/(n_iter+1.)\n",
    "    two_sided_pval = min(2 * min(left_pval, right_pval), 1)\n",
    "    if boold_p_dist:\n",
    "        return {'left':left_pval, 'right':right_pval, 'two-sided':two_sided_pval}, stat, perm_stat\n",
    "    else:\n",
    "        return {'left':left_pval, 'right':right_pval, 'two-sided':two_sided_pval}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC\n",
      "ST-HSC\n",
      "MPP\n",
      "{'SLAM-HSC': {'left': 0.842688629245483, 'right': 0.15731537073851704, 'two-sided': 0.31463074147703407}, 'ST-HSC': {'left': 0.999996000016, 'right': 7.999968000127999e-06, 'two-sided': 1.5999936000255998e-05}, 'MPP': {'left': 1.0, 'right': 3.9999840000639995e-06, 'two-sided': 7.999968000127999e-06}}\n"
     ]
    }
   ],
   "source": [
    "#FIG 2B\n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    print(oc)\n",
    "    df_temp = df[(df.Original_cell==oc)]\n",
    "\n",
    "    data = np.array([df_temp[(df_temp.Family==fam)].Class.values\n",
    "                     for fam in np.unique(df_temp.Family)])\n",
    "    perm_label = np.array([df_temp[(df_temp.Family==fam)].Culture_condition.iloc[0]\n",
    "                           for fam in np.unique(df_temp.Family)])\n",
    "    \n",
    "    table_support = np.unique(np.hstack(data))\n",
    "    all_perm_labels = np.unique(perm_label)\n",
    "    bool_stack = True\n",
    "    args = (table_support, all_perm_labels, bool_stack,)\n",
    "    \n",
    "    pvals_dct[oc] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "print(pvals_dct)\n",
    "with open('./pickled_data/Fig2B_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC\n",
      "ST-HSC\n",
      "MPP\n",
      "{'SLAM-HSC': {'left': 0.9653961384154464, 'right': 0.03461186155255379, 'two-sided': 0.06922372310510758}, 'ST-HSC': {'left': 0.9923960304158783, 'right': 0.007607969568121728, 'two-sided': 0.015215939136243455}, 'MPP': {'left': 0.9810240759036963, 'right': 0.018983924064303742, 'two-sided': 0.037967848128607484}}\n"
     ]
    }
   ],
   "source": [
    "###FIG 2C\n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    print(oc)\n",
    "    df_temp = df[(df.Original_cell==oc)]\n",
    "\n",
    "    data = np.array([df_temp[(df_temp.Family==fam)].Generation.max()\n",
    "                     for fam in np.unique(df_temp.Family)])\n",
    "    perm_label = np.array([df_temp[(df_temp.Family==fam)].Culture_condition.iloc[0]\n",
    "                           for fam in np.unique(df_temp.Family)])\n",
    "    \n",
    "    table_support = np.unique(data)\n",
    "    all_perm_labels = np.unique(perm_label)\n",
    "    bool_stack = False\n",
    "    args = (table_support, all_perm_labels, bool_stack,)\n",
    "    \n",
    "    pvals_dct[oc] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "print(pvals_dct)\n",
    "        \n",
    "with open('./pickled_data/Fig2C_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'MPP': {'left': 0.9810240759036963,\n",
       "  'right': 0.018983924064303742,\n",
       "  'two-sided': 0.037967848128607484},\n",
       " 'SLAM-HSC': {'left': 0.9653961384154464,\n",
       "  'right': 0.03461186155255379,\n",
       "  'two-sided': 0.06922372310510758},\n",
       " 'ST-HSC': {'left': 0.9923960304158783,\n",
       "  'right': 0.007607969568121728,\n",
       "  'two-sided': 0.015215939136243455}}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with open('./pickled_data/Fig2C_pvals', 'rb') as fp:\n",
    "    pvals_dct = pickle.load(fp)\n",
    "pvals_dct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC\n",
      "ST-HSC\n",
      "MPP\n",
      "{'SLAM-HSC': {'left': 0.646501413994344, 'right': 0.35420258318966724, 'two-sided': 0.7084051663793345}, 'ST-HSC': {'left': 0.8933844264622941, 'right': 0.10661957352170591, 'two-sided': 0.21323914704341182}, 'MPP': {'left': 0.9981200075199699, 'right': 0.0018839924640301439, 'two-sided': 0.0037679849280602877}}\n"
     ]
    }
   ],
   "source": [
    "###FIG 2D\n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    print(oc)\n",
    "    df_temp = df[(df.Original_cell==oc)&(df.Generation==0)]\n",
    "\n",
    "    data = np.array([df_temp[(df_temp.Family==fam)].Class.iloc[0]\n",
    "                     for fam in np.unique(df_temp.Family)])\n",
    "    perm_label = np.array([df_temp[(df_temp.Family==fam)].Culture_condition.iloc[0]\n",
    "                           for fam in np.unique(df_temp.Family)])\n",
    "    \n",
    "    table_support = np.unique(data)\n",
    "    all_perm_labels = np.unique(perm_label)\n",
    "    bool_stack = False\n",
    "    args = (table_support, all_perm_labels, bool_stack,)\n",
    "    \n",
    "    pvals_dct[oc] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "print(pvals_dct)\n",
    "        \n",
    "with open('./pickled_data/Fig2D_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC\n",
      "ST-HSC\n",
      "MPP\n",
      "{'SLAM-HSC': {'left': 0.23293906824372704, 'right': 0.7787488850044599, 'two-sided': 0.4658781364874541}, 'ST-HSC': {'left': 0.6389534441862232, 'right': 0.38085447658209365, 'two-sided': 0.7617089531641873}, 'MPP': {'left': 0.8721325114699541, 'right': 0.1326434694261223, 'two-sided': 0.2652869388522446}}\n"
     ]
    }
   ],
   "source": [
    "#FIG 2E\n",
    "sym_labs = ['SYM UNDIF', 'SYM DIF', 'ASYM UNDIF', 'ASYM DIF']\n",
    "\n",
    "def first_div_class(df):#sis_type, prog_type):\n",
    "    daugher1, daugher2 = df.Class.values\n",
    "    progenitor = df.Original_cell.values[0]\n",
    "    if daugher1 == daugher2:\n",
    "        if daugher1 == progenitor:\n",
    "            return sym_labs[0]\n",
    "        else:\n",
    "            return sym_labs[1]\n",
    "    else:\n",
    "        if (daugher1 == progenitor) or (daugher2 == progenitor):\n",
    "            return sym_labs[2]\n",
    "        else:\n",
    "            return sym_labs[3]\n",
    "    \n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    print(oc)\n",
    "    df_temp = df[(df.Original_cell==oc)&(df.Generation==1)]\n",
    "\n",
    "    data = np.array([first_div_class(df_temp[df_temp.Family==fam])\n",
    "                     for fam in np.unique(df_temp.Family) if len(df_temp[df_temp.Family==fam])==2])\n",
    "    perm_label = np.array([df_temp[(df_temp.Family==fam)].Culture_condition.iloc[0]\n",
    "                           for fam in np.unique(df_temp.Family) if len(df_temp[df_temp.Family==fam])==2])\n",
    "    \n",
    "    table_support = np.unique(data)\n",
    "    all_perm_labels = np.unique(perm_label)\n",
    "    bool_stack = False\n",
    "    args = (table_support, all_perm_labels, bool_stack,)\n",
    "    \n",
    "    pvals_dct[oc] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "print(pvals_dct)\n",
    "        \n",
    "with open('./pickled_data/Fig2E_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC\n",
      "ST-HSC\n",
      "MPP\n",
      "{'SLAM-HSC': {'left': 0.885312458750165, 'right': 0.11469154123383507, 'two-sided': 0.22938308246767014}, 'ST-HSC': {'left': 0.9996360014559942, 'right': 0.00036799852800588795, 'two-sided': 0.0007359970560117759}, 'MPP': {'left': 0.9999880000479998, 'right': 1.5999936000255998e-05, 'two-sided': 3.1999872000511996e-05}}\n"
     ]
    }
   ],
   "source": [
    "#SUPPLEMENTARY FIG 2B\n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    print(oc)\n",
    "    df_temp = df[(df.Original_cell==oc)&(df.Generation==1)]\n",
    "\n",
    "    data = np.array([df_temp[(df_temp.Family==fam)].Class.values\n",
    "                     for fam in np.unique(df_temp.Family)])\n",
    "    perm_label = np.array([df_temp[(df_temp.Family==fam)].Culture_condition.iloc[0]\n",
    "                           for fam in np.unique(df_temp.Family)])\n",
    "    \n",
    "    table_support = np.unique(np.hstack(data))\n",
    "    all_perm_labels = np.unique(perm_label)\n",
    "    bool_stack = True\n",
    "    args = (table_support, all_perm_labels, bool_stack,)\n",
    "    \n",
    "    pvals_dct[oc] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "print(pvals_dct)\n",
    "with open('./pickled_data/SuppFig2B_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC P 24h\n",
      "SLAM-HSC P 48h\n",
      "SLAM-HSC P+ILs 24h\n",
      "SLAM-HSC P+ILs 48h\n",
      "ST-HSC P 24h\n",
      "ST-HSC P 48h\n",
      "ST-HSC P+ILs 24h\n",
      "ST-HSC P+ILs 48h\n",
      "MPP P 24h\n",
      "MPP P 48h\n",
      "MPP P+ILs 24h\n",
      "MPP P+ILs 48h\n"
     ]
    }
   ],
   "source": [
    "#SANITY CHECK: TEST FOR WELL AND EXPERIMENT EFFECT TO FAMILIAL CELL CLASS DISTIRBUTION\n",
    "#(IN FIXED CONDITION, TIME AND PROGENITOR CLASS)\n",
    "#FAMILIES ARE PERMUTED ACROSS WELLS AND EXPERIMENTS AFTER THE WELL_EXPERIMENT LABEL\n",
    "#THE STATISTIC USED IS THE G-TEST STATISTIC ON THE CONTINGENCY TABLE OF CELL COUNTS\n",
    "#SORTED BY CELL CLASS VS WELL_EXPERIMENT LABEL\n",
    "\n",
    "pvals_dct = {}\n",
    "for oc in or_cells:\n",
    "    pvals_dct[oc] = {}\n",
    "    for cnd in conds:\n",
    "        pvals_dct[oc][cnd] = {}\n",
    "        for t in times:\n",
    "            print(oc, cnd, t)\n",
    "            df_temp = df[(df.Original_cell==oc)&(df.Culture_condition==cnd)&(df.Culture_time==t)]\n",
    "\n",
    "            data = np.array([df_temp[(df_temp.Family==fam)].Class.values\n",
    "                             for fam in np.unique(df_temp.Family)])\n",
    "            perm_label = np.array([df_temp[(df_temp.Family==fam)].Well_experiment.iloc[0]\n",
    "                                   for fam in np.unique(df_temp.Family)])\n",
    "\n",
    "            table_support = np.unique(np.hstack(data))\n",
    "            all_perm_labels = np.unique(perm_label)\n",
    "            bool_stack = True\n",
    "            args = (table_support, all_perm_labels, bool_stack,)\n",
    "\n",
    "            pvals_dct[oc][cnd][t] = permutation_testing(chi_squared_stat, args, data, perm_label)\n",
    "with open('./pickled_data/WellExperiment_test_pvals', 'wb') as fp:\n",
    "    pickle.dump(pvals_dct, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLAM-HSC P 24h \t 0.008683965264138943\n",
      "SLAM-HSC P 48h \t 0.5123499506001976\n",
      "SLAM-HSC P+ILs 24h \t 0.2304630781476874\n",
      "SLAM-HSC P+ILs 48h \t 0.6260374958500166\n",
      "ST-HSC P 24h \t 0.2593549625801497\n",
      "ST-HSC P 48h \t 0.0967436130255479\n",
      "ST-HSC P+ILs 24h \t 0.6066975732097072\n",
      "ST-HSC P+ILs 48h \t 0.7582449670201319\n",
      "MPP P 24h \t 0.24542301830792676\n",
      "MPP P 48h \t 0.6022295910816357\n",
      "MPP P+ILs 24h \t 0.37274650901396394\n",
      "MPP P+ILs 48h \t 0.3910064359742561\n"
     ]
    }
   ],
   "source": [
    "with open('./pickled_data/WellExperiment_test_pvals', 'rb') as fp:\n",
    "    pvals_dct = pickle.load(fp)\n",
    "for oc in or_cells:\n",
    "    for cnd in conds:\n",
    "        for t in times:\n",
    "            print(oc, cnd, t, '\\t', pvals_dct[oc][cnd][t]['right'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}