3-CC_Distribution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "e4a7dc6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import itertools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import pearsonr\n",
    "from scipy.signal import chirp, find_peaks, peak_widths\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "194e7a90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# prep info\n",
    "# contact: c\n",
    "# status inlcude: r, s, a\n",
    "\n",
    "root_path = r'/Users/xiaoqiansun/Desktop/MedLu/TubeTest/Data'\n",
    "\n",
    "\n",
    "days = 6\n",
    "fpsN = 15\n",
    "fpsB = 30\n",
    "threshold = 0.1\n",
    "baselinePeriod = 180 #(3min(180s) baseline)\n",
    "redundentTime = 190 #10s+3min\n",
    "\n",
    "\n",
    "M32_list = ['32-1_28', '32-2_18', '32-3_15', '32-4_20', '32-5_12', '32-6_15']\n",
    "M33_list = ['33-1_28', '33-2_35', '33-3_32', '33-4_30', '33-5_30', '33-6_26']\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b34b869",
   "metadata": {},
   "source": [
    "# Define function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "e6039a13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to keep certain status frrame interval\n",
    "#----------------------------------------------------------------------------------\n",
    "def pick_BInterval(m1_tube, m2_tube, s1, s2, c=None):      \n",
    "      \n",
    "    df1 = m1_tube.copy()\n",
    "    df2 = m2_tube.copy()\n",
    "    \n",
    "    if c:\n",
    "        df1 = df1[(df1['mouse1']== s1) & (df1['mouse2']== s2) & (df1['contact']== c)]\n",
    "        df2 = df2[(df2['mouse1']== s1) & (df2['mouse2']== s2) & (df2['contact']== c)]\n",
    "    \n",
    "    else:\n",
    "        df1 = df1[(df1['mouse1']== s1) & (df1['mouse2']== s2)]\n",
    "        df2 = df2[(df2['mouse1']== s1) & (df2['mouse2']== s2)]\n",
    "    \n",
    "        \n",
    "    df1 = df1.drop(['mouse1', 'mouse2', 'contact'], axis=1)\n",
    "    df2 = df2.drop(['mouse1', 'mouse2', 'contact'], axis=1)\n",
    "    \n",
    "    return(df1, df2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "29f10d33",
   "metadata": {},
   "outputs": [],
   "source": [
    "#function to get corr between each two cell of a dataframe\n",
    "#----------------------------------------------------------------------------------\n",
    "def corr_list(df, threshold):\n",
    "    end_range = df.shape[1]\n",
    "    \n",
    "    corr_list = []\n",
    "    corrP_list = []\n",
    "    corrN_list = []\n",
    "    \n",
    "    if df.shape[0] >2:\n",
    "        for i in range(end_range):\n",
    "            a = i+1\n",
    "            for a in range(a,end_range):\n",
    "                corr, _ = pearsonr(df.iloc[:,i], df.iloc[:,a])\n",
    "                \n",
    "                if corr >= threshold:\n",
    "                    corrP_list.append(corr)\n",
    "                elif corr <= -threshold:\n",
    "                    corrN_list.append(corr)\n",
    "                \n",
    "                corr_list.append(corr)\n",
    "    \n",
    "    return(corr_list, corrP_list, corrN_list)\n",
    "\n",
    "\n",
    "#function to calculate FWHM\n",
    "#----------------------------------------------------------------------------------\n",
    "def FWHM(distribute):\n",
    "    \n",
    "    # calculate sigma and mu for best fit\n",
    "    sigma = np.nanstd(distribute) # standard deviation of distribution\n",
    "    mu = np.nanmean(distribute) # mean of the distribution\n",
    "    \n",
    "    fig, ax = plt.subplots()\n",
    "    n, bins, patches = ax.hist(distribute, bins=100, density=1)\n",
    "    \n",
    "    # generate best fit curve for the histgram\n",
    "    y = ((1 / (np.sqrt(2 * np.pi) * sigma)) *np.exp(-0.5 * (1 / sigma * (bins - mu))**2))\n",
    "\n",
    "    # calculate FWHM \n",
    "    peaks, _ = find_peaks(y)\n",
    "    results_half = peak_widths(y, peaks, rel_height=0.5)\n",
    "    FWHM = results_half[0]\n",
    "    \n",
    "    plt.close()\n",
    "    \n",
    "    return(FWHM, sigma, mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "502dc9da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def CC_inervalDays(df1, df2, behavior1, behavior2, chebavior, CC_32, CC_33, day, threshold):\n",
    "    \n",
    "    bdf1, bdf2 = pick_BInterval(df1, df2, behavior1, behavior2, c=chebavior)\n",
    "    \n",
    "    bdf1_cc, bdf1_ccP, bdf1_ccN = corr_list(bdf1,threshold)\n",
    "    bdf2_cc, bdf2_ccP, bdf2_ccN = corr_list(bdf2,threshold)\n",
    "    \n",
    "    CC_32.loc[day, CC_32.columns[1:]] = [np.mean(bdf1_ccP), len(bdf1_ccP)/(len(bdf1_cc)+0.0001), \n",
    "                                         np.mean(bdf1_ccN), len(bdf1_ccN)/(len(bdf1_cc)+0.0001)]\n",
    "    CC_33.loc[day, CC_33.columns[1:]] = [np.mean(bdf2_ccP), len(bdf2_ccP)/(len(bdf2_cc)+0.0001), \n",
    "                                         np.mean(bdf2_ccN), len(bdf2_ccN)/(len(bdf2_cc)+0.0001)]\n",
    "    \n",
    "\n",
    "    return(CC_32, CC_33, bdf1_cc, bdf2_cc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d44b0637",
   "metadata": {},
   "source": [
    "# CC Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "8f58c980",
   "metadata": {},
   "outputs": [],
   "source": [
    "c32_list = []\n",
    "c33_list = []\n",
    "sheetNames = ['as', 'sa', 'rs', 'sr', 'ar', 'ra', 'ssnc', 'ssc']\n",
    "for i in range(len(sheetNames)):\n",
    "    c32_list.append(pd.DataFrame({'Days':np.arange(1,7), '01P':np.nan, '01P%': np.nan, '01N':np.nan, '01N%': np.nan}))\n",
    "    c33_list.append(pd.DataFrame({'Days':np.arange(1,7), '01P':np.nan, '01P%': np.nan, '01N':np.nan, '01N%': np.nan}))\n",
    "    \n",
    "cc32_chain = [ [] for _ in range(len(sheetNames)) ]\n",
    "cc33_chain = [ [] for _ in range(len(sheetNames)) ]\n",
    "\n",
    "# update dataframe \n",
    "for bb in range(len(sheetNames)):\n",
    "\n",
    "    b = sheetNames[bb]\n",
    "    behavior1 = b[0]\n",
    "    behavior2 = b[1]\n",
    "    if len(b) != 2:\n",
    "        chebavior = b[2:]\n",
    "    else:\n",
    "        chebavior = None\n",
    "\n",
    "\n",
    "    # each day in this behavior combination\n",
    "    for day in range(days):\n",
    "\n",
    "        m32_folder = os.path.join(root_path, M32_list[day])\n",
    "        m33_folder = os.path.join(root_path, M33_list[day])\n",
    "\n",
    "\n",
    "        m32_tube = pd.read_csv(os.path.join(m32_folder,'m32_tube.csv'), index_col = 'Frame')\n",
    "        m33_tube = pd.read_csv(os.path.join(m33_folder,'m33_tube.csv'), index_col = 'Frame') \n",
    "\n",
    "        # reset\n",
    "        c32_list[bb],c33_list[bb],c32,c33 = CC_inervalDays(m32_tube, m33_tube, \n",
    "                                                           behavior1, behavior2, chebavior, \n",
    "                                                           c32_list[bb], c33_list[bb], day, threshold)\n",
    "        cc32_chain[bb] = itertools.chain(cc32_chain[bb],c32); cc33_chain[bb] = itertools.chain(cc33_chain[bb],c33)\n",
    "\n",
    "\n",
    "        \n",
    "# write to summary\n",
    "writer = pd.ExcelWriter('CC.xlsx')\n",
    "for k in range(len(sheetNames)):\n",
    "    c32_list[k].to_excel(writer, sheet_name='Mouse 32 '+sheetNames[k])\n",
    "    c33_list[k].to_excel(writer, sheet_name='Mouse 33 '+sheetNames[k])\n",
    "writer.save()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce74930f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1a9ddc10",
   "metadata": {},
   "source": [
    "# Distribution, FWHM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "965a8bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Overlay 2 histograms to compare them\n",
    "#------------------------\n",
    "\n",
    "def overlaid_histogram(figsize, data1, data2, n_bins = 0, data1_name=\"\", data1_color=\"#539caf\", \n",
    "                       data2_name=\"\", data2_color=\"#7663b0\", x_label=\"\", y_label=\"\", title=\"\", y_limit=None):\n",
    "    # Set the bounds for the bins so that the two distributions are fairly compared\n",
    "\n",
    "    if len(data1) > 0 and len(data2) > 0:\n",
    "        max_nbins = 10\n",
    "        data_range = [min(min(data1), min(data2)), max(max(data1), max(data2))]\n",
    "        binwidth = (data_range[1] - data_range[0]) / max_nbins\n",
    "\n",
    "\n",
    "        if n_bins == 0:\n",
    "            bins = np.arange(data_range[0], data_range[1] + binwidth, binwidth)\n",
    "        else: \n",
    "            bins = n_bins\n",
    "\n",
    "        weights1 = np.ones_like(data1)/float(len(data1))\n",
    "        weights2 = np.ones_like(data2)/float(len(data2))\n",
    "\n",
    "        # Create the plot\n",
    "        _, ax = plt.subplots(figsize = figsize)\n",
    "        ax.hist(data1, bins = bins, color = data1_color, alpha = 1, \n",
    "                label = data1_name, weights=weights1)\n",
    "        ax.hist(data2, bins = bins, color = data2_color, alpha = 0.5, \n",
    "                label = data2_name ,weights=weights2)\n",
    "\n",
    "        if y_limit:\n",
    "            ax.set_ylim(0, y_limit)\n",
    "\n",
    "        # https://stackoverflow.com/questions/925024/how-can-i-remove-the-top-and-right-axis-in-matplotlib\n",
    "        ax.spines['right'].set_visible(False)\n",
    "        ax.spines['top'].set_visible(False)\n",
    "\n",
    "        # https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.axes.Axes.tick_params.html\n",
    "        # https://matplotlib.org/3.1.1/gallery/text_labels_and_annotations/fonts_demo.html\n",
    "        ax.tick_params(axis='both', which='major', labelsize=20)\n",
    "        labels = ax.get_yticklabels()\n",
    "        [label.set_fontname('Arial') for label in labels]\n",
    "        labels = ax.get_xticklabels()\n",
    "        [label.set_fontname('Arial') for label in labels]\n",
    "        #ax.set_ylabel(y_label,fontname=\"Arial\",fontsize = 20)\n",
    "        #ax.set_xlabel(x_label,fontname=\"Arial\",fontsize = 20)\n",
    "        #ax.set_title(title, fontsize = 16, fontweight = 'bold')\n",
    "        #ax.legend(loc = 'best')\n",
    "\n",
    "        plt.savefig(title)\n",
    "        \n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "52f783e7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "at as mouse 32 has len= 556 and mouse 33 has len= 2286\n",
      "at sa mouse 32 has len= 771 and mouse 33 has len= 2069\n",
      "at rs mouse 32 has len= 907 and mouse 33 has len= 2664\n",
      "at sr mouse 32 has len= 171 and mouse 33 has len= 435\n",
      "at ar mouse 32 has len= 0 and mouse 33 has len= 0\n",
      "at ra mouse 32 has len= 534 and mouse 33 has len= 1199\n",
      "at ssnc mouse 32 has len= 907 and mouse 33 has len= 2664\n",
      "at ssc mouse 32 has len= 841 and mouse 33 has len= 2229\n"
     ]
    }
   ],
   "source": [
    "FWHM_32 = []\n",
    "FWHM_33 = []\n",
    "figsize = (6,6)\n",
    "\n",
    "for i in range(len(sheetNames)):\n",
    "    behavior = sheetNames[i]\n",
    "    cc32L = list(cc32_chain[i])\n",
    "    cc33L = list(cc33_chain[i])\n",
    "    print('at', behavior, 'mouse 32 has len=', len(cc32L), 'and mouse 33 has len=', len(cc33L))\n",
    "    overlaid_histogram(figsize,cc32L, cc33L, n_bins = 200, \n",
    "                       data1_name=\"Mouse 32\", data1_color=\"orange\", \n",
    "                       data2_name=\"Mouse 33\", data2_color=\"blue\", \n",
    "                       x_label=\"Pearson Correlation Coefficient Distribution\", y_label=\"Probability\",\n",
    "                       title=\"Distribution At \" + behavior)\n",
    "    \n",
    "    FWHM_32.append(FWHM(cc32L)[0])\n",
    "    FWHM_33.append(FWHM(cc33L)[0])\n",
    "    \n",
    "FWHM_Summary = pd.DataFrame({'behavior':sheetNames, 'Mouse 32 FWHM':FWHM_32, 'Mouse 33 FWHM':FWHM_33})\n",
    "FWHM_Summary.to_excel('FWHM_Summary.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "411cfacc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d1cba5",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}