CSS844_SVA_Group_PlottingCleanedMatrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSS844"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code in this notebook takes the \"CleanSV_X.csv\" files from an R script included in the GitHub repository and generates PCA, MDS, and t-SNE plots of these datasets with color-coding for the family, tissue, and stress factors. Note that the visualization code is adapted from Rei Doko's code from HRT841. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#imports\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.manifold import MDS\n",
    "from scipy import stats\n",
    "from sklearn.manifold import TSNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading in data\n",
    "# Note that the CleanSV datasets come from an R script included with this notebook\n",
    "\n",
    "master_data_list = []\n",
    "for i in range(16): \n",
    "    filename = \"CleanSV_{0}.csv\".format(i+1)\n",
    "    to_import = pd.read_csv(filename)\n",
    "    master_data_list.append(to_import)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Necessary functions for visualization, adapted from Rei Doko's code\n",
    "\n",
    "def apply_basic_filter(df):\n",
    "    filtered_df = df.copy() # so we don't modify the original dataframe if we didn't want to\n",
    "   \n",
    "    # (Lifted from RNAseq_Meta_analysis.ipynb)\n",
    "    # Filter out non-expressed genes\n",
    "    filtered_df = filtered_df.loc[filtered_df.sum(axis=1) > 10, :]\n",
    "    # Filter out lowly expressed genes\n",
    "    mask_low_vals = (filtered_df > 1).sum(axis=1) > 10\n",
    "    filtered_df = filtered_df.loc[mask_low_vals, :]\n",
    "    \n",
    "    return filtered_df\n",
    "\n",
    "def log_transform(df):\n",
    "    log_df = df.copy() # so we don't modify the original dataframe if we didn't want to\n",
    "    \n",
    "    # (Lifted from RNAseq_Meta_analysis.ipynb)\n",
    "    for c in [c for c in log_df.columns if np.issubdtype(log_df[c].dtype , np.number)]:\n",
    "        log_df[c] += 1\n",
    "    for c in [c for c in log_df.columns if np.issubdtype(log_df[c].dtype , np.number)]:\n",
    "        log_df[c] = np.log(log_df[c])\n",
    "    return log_df\n",
    "\n",
    "def apply_pca(data, num_c=2, return_pca=False):\n",
    "    \"\"\"\n",
    "    data : data to fit PCA on\n",
    "    num_c : Number of principal components for PCA\n",
    "    return_pca : should the function also return the PCA object (for looking at things like explained variance ratio) \n",
    "    \"\"\"\n",
    "    pca = PCA(n_components=num_c)\n",
    "    real_PCs = pca.fit_transform(data)\n",
    "    real_PCs_df = pd.DataFrame(data=real_PCs, columns=[f'PC{i}' for i in range(1, num_c+1)])\n",
    "    if return_pca:\n",
    "        return real_PCs_df, pca\n",
    "    else:\n",
    "        return real_PCs_df\n",
    "\n",
    "def zscore_dataframe(df, axis=0):\n",
    "    \"\"\"\n",
    "    df : pandas dataframe\n",
    "    axis : which axis to zscore on (0 for zscore across each row, 1 for across each column)\n",
    "    \"\"\"\n",
    "    zscored_data = stats.zscore(df, axis=axis) # find the zscore of the data (returns a numpy array)\n",
    "    zscore_df = pd.DataFrame(zscored_data, columns=df.columns) # turn it back into a dataframe\n",
    "    zscore_df = zscore_df.set_index(df.index) # set the orthogroup as the index\n",
    "    return zscore_df\n",
    "\n",
    "def scatter_plot(df, x, y, ax, labels=None, title=\"\", palette='tab10'):\n",
    "    sns.scatterplot(x=x, y=y, data=df,\n",
    "                    ax=ax,\n",
    "                    alpha=0.7, # opacity of the points\n",
    "                    linewidth=0, # outline of the data points\n",
    "                    hue=labels,\n",
    "                    palette=palette)\n",
    "    ax.set_title(title)\n",
    "    \n",
    "def apply_mds(data, num_dim=2, return_object=False):\n",
    "    \"\"\"\n",
    "    data : data to fit MDS on\n",
    "    num_dim : Number of dimensions of reduced data\n",
    "    return_object : should the function also return the MDS object\n",
    "    \"\"\"\n",
    "    mds = MDS(n_components=num_dim)\n",
    "    real_MDS = mds.fit_transform(data)\n",
    "    real_MDS_df = pd.DataFrame(data=real_MDS, columns=[f'MDS{i}' for i in range(1, num_dim+1)])\n",
    "    if return_object:\n",
    "        return real_MDS_df, mds\n",
    "    else:\n",
    "        return real_MDS_df\n",
    "    \n",
    "\n",
    "    \n",
    "def apply_tsne(data, num_dim=2, return_object=False):\n",
    "    \"\"\"\n",
    "    data : data to fit t-SNE on\n",
    "    num_dim : Number of dimensions of reduced data\n",
    "    return_object : should the function also return the t-SNE object \n",
    "    \"\"\"\n",
    "    tsne = TSNE(n_components=num_dim)\n",
    "    real_TSNE = tsne.fit_transform(data)\n",
    "    real_TSNE_df = pd.DataFrame(data=real_TSNE, columns=[f'tSNE{i}' for i in range(1, num_dim+1)])\n",
    "    if return_object:\n",
    "        return real_TSNE_df, tsne\n",
    "    else:\n",
    "        return real_TSNE_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-f68627a5b908>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     44\u001b[0m     \u001b[0mlabel_df\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     45\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 46\u001b[1;33m     \u001b[0mtsne_df\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mapply_tsne\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mexpr_df_zscore\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     47\u001b[0m     \u001b[0mmds_df\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mapply_mds\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mexpr_df_zscore\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     48\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-6-08fad627d85a>\u001b[0m in \u001b[0;36mapply_tsne\u001b[1;34m(data, num_dim, return_object)\u001b[0m\n\u001b[0;32m     79\u001b[0m     \"\"\"\n\u001b[0;32m     80\u001b[0m     \u001b[0mtsne\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTSNE\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_components\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnum_dim\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 81\u001b[1;33m     \u001b[0mreal_TSNE\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtsne\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     82\u001b[0m     \u001b[0mreal_TSNE_df\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mreal_TSNE\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34mf'tSNE{i}'\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum_dim\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     83\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mreturn_object\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\sklearn\\manifold\\t_sne.py\u001b[0m in \u001b[0;36mfit_transform\u001b[1;34m(self, X, y)\u001b[0m\n\u001b[0;32m    893\u001b[0m             \u001b[0mEmbedding\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mtraining\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mlow\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mdimensional\u001b[0m \u001b[0mspace\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    894\u001b[0m         \"\"\"\n\u001b[1;32m--> 895\u001b[1;33m         \u001b[0membedding\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_fit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    896\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0membedding_\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0membedding\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    897\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0membedding_\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\sklearn\\manifold\\t_sne.py\u001b[0m in \u001b[0;36m_fit\u001b[1;34m(self, X, skip_num_points)\u001b[0m\n\u001b[0;32m    766\u001b[0m             \u001b[0mt0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    767\u001b[0m             distances_nn, neighbors_nn = knn.kneighbors(\n\u001b[1;32m--> 768\u001b[1;33m                 None, n_neighbors=k)\n\u001b[0m\u001b[0;32m    769\u001b[0m             \u001b[0mduration\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mt0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    770\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mverbose\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\sklearn\\neighbors\\base.py\u001b[0m in \u001b[0;36mkneighbors\u001b[1;34m(self, X, n_neighbors, return_distance)\u001b[0m\n\u001b[0;32m    452\u001b[0m                 delayed_query(\n\u001b[0;32m    453\u001b[0m                     self._tree, X[s], n_neighbors, return_distance)\n\u001b[1;32m--> 454\u001b[1;33m                 \u001b[1;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mgen_even_slices\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_jobs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    455\u001b[0m             )\n\u001b[0;32m    456\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, iterable)\u001b[0m\n\u001b[0;32m    919\u001b[0m             \u001b[1;31m# remaining jobs.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    920\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_iterating\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 921\u001b[1;33m             \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdispatch_one_batch\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0miterator\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    922\u001b[0m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_iterating\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_original_iterator\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    923\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36mdispatch_one_batch\u001b[1;34m(self, iterator)\u001b[0m\n\u001b[0;32m    757\u001b[0m                 \u001b[1;32mreturn\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    758\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 759\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_dispatch\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtasks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    760\u001b[0m                 \u001b[1;32mreturn\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    761\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36m_dispatch\u001b[1;34m(self, batch)\u001b[0m\n\u001b[0;32m    714\u001b[0m         \u001b[1;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    715\u001b[0m             \u001b[0mjob_idx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_jobs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 716\u001b[1;33m             \u001b[0mjob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mapply_async\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbatch\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    717\u001b[0m             \u001b[1;31m# A job can complete so quickly than its callback is\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    718\u001b[0m             \u001b[1;31m# called before we get here, causing self._jobs to\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\_parallel_backends.py\u001b[0m in \u001b[0;36mapply_async\u001b[1;34m(self, func, callback)\u001b[0m\n\u001b[0;32m    180\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mapply_async\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    181\u001b[0m         \u001b[1;34m\"\"\"Schedule a func to be run\"\"\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 182\u001b[1;33m         \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mImmediateResult\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    183\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    184\u001b[0m             \u001b[0mcallback\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\_parallel_backends.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, batch)\u001b[0m\n\u001b[0;32m    547\u001b[0m         \u001b[1;31m# Don't delay the application, to avoid keeping the input\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    548\u001b[0m         \u001b[1;31m# arguments in memory\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 549\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbatch\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    550\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    551\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    223\u001b[0m         \u001b[1;32mwith\u001b[0m \u001b[0mparallel_backend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_jobs\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_n_jobs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    224\u001b[0m             return [func(*args, **kwargs)\n\u001b[1;32m--> 225\u001b[1;33m                     for func, args, kwargs in self.items]\n\u001b[0m\u001b[0;32m    226\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    227\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__len__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m    223\u001b[0m         \u001b[1;32mwith\u001b[0m \u001b[0mparallel_backend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_jobs\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_n_jobs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    224\u001b[0m             return [func(*args, **kwargs)\n\u001b[1;32m--> 225\u001b[1;33m                     for func, args, kwargs in self.items]\n\u001b[0m\u001b[0;32m    226\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    227\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__len__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\sklearn\\neighbors\\base.py\u001b[0m in \u001b[0;36m_tree_query_parallel_helper\u001b[1;34m(tree, data, n_neighbors, return_distance)\u001b[0m\n\u001b[0;32m    289\u001b[0m     \u001b[0munder\u001b[0m \u001b[0mPyPy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    290\u001b[0m     \"\"\"\n\u001b[1;32m--> 291\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mtree\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mquery\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_neighbors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreturn_distance\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    292\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    293\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Loop to generate PCA, MDS, and tSNE plots with increasing # of surrogate variables\n",
    "\n",
    "# (1) Do the filtering and all that stuff up until the PCA, MDS, and tSNE calculations, then\n",
    "# (2) loop across family, tissue, stress to generate three 1x3 plots with the PCA, MDS, and tSNE results \n",
    "\n",
    "for i in range(15):\n",
    "    \n",
    "    print('Iteration',i)\n",
    "    raw_expr_df = master_data_list[i] # Change as needed to where you have your copy of the dataset\n",
    "    raw_expr_df = raw_expr_df.set_index('Orthogroup')\n",
    "\n",
    "    expr_df = apply_basic_filter(raw_expr_df)\n",
    "\n",
    "    expr_df_trans = expr_df.transpose()\n",
    "\n",
    "    real_PCs_df = apply_pca(expr_df_trans)\n",
    "\n",
    "    expr_df_zscore = zscore_dataframe(expr_df_trans, axis=1)\n",
    "\n",
    "    z_real_PCs_df = apply_pca(expr_df_zscore)\n",
    "\n",
    "    raw_labels = pd.read_csv(\"factors_v1_wPloidyInfo_Final.csv\")\n",
    "    raw_labels = raw_labels.set_index('sra')\n",
    "\n",
    "    labeling_rows = []\n",
    "\n",
    "    columns = [('Ploidy_low', 'Ploidy_low'), \n",
    "               ('Ploidy_med', 'Ploidy_med'), \n",
    "               ('Ploidy_med_high', 'Ploidy_med_high'),\n",
    "               ('Ploidy_high','Ploidy_high'),\n",
    "               ('Ploidy_AncientPolyploids','Ploidy_AncientPolyploids'),\n",
    "              ('tissue','tissue'),\n",
    "              ('family','family'),\n",
    "              ('stress','stress'),]\n",
    "\n",
    "    for srr in list(expr_df_zscore.index):\n",
    "        entry = raw_labels.loc[srr]\n",
    "        row = {'SRA' : srr}\n",
    "        for old_colname, new_colname in columns:\n",
    "            row[new_colname] = entry[old_colname]\n",
    "        labeling_rows.append(row)\n",
    "\n",
    "    label_df = pd.DataFrame(labeling_rows)\n",
    "    label_df\n",
    "\n",
    "    tsne_df = apply_tsne(expr_df_zscore)\n",
    "    mds_df = apply_mds(expr_df_zscore)\n",
    "\n",
    "    #fig loop\n",
    "    plot_types = ['family','tissue','stress']\n",
    "    \n",
    "    for j in range(3):\n",
    "        plot_type = plot_types[j]\n",
    "        labels = label_df[plot_type]\n",
    "\n",
    "        fig, axs = plt.subplots(1,3, figsize=(24,8), dpi=100)\n",
    "        axs=axs.flatten()\n",
    "\n",
    "        scatter_plot(z_real_PCs_df, 'PC1', 'PC2', axs[0], labels, \"PCA on expression z-scores\")\n",
    "        scatter_plot(mds_df, 'MDS1', 'MDS2', axs[1], labels, \"MDS on expression z-scores\")\n",
    "        scatter_plot(tsne_df, 'tSNE1', 'tSNE2', axs[2], labels, \"t-SNE on expression z-scores\")\n",
    "\n",
    "        for ax in axs:\n",
    "            if ax.get_legend(): # legend could be Nonetype\n",
    "                ax.get_legend().set_visible(False)\n",
    "        axs[2].legend(bbox_to_anchor=(1.01, 0.55), loc=2, borderaxespad=0., frameon=False)\n",
    "        #fig.suptitle('Factorized by Family',size=20)\n",
    "        filename = 'CleanSV_Dataset_{0}_plot_type_{1}_afterOutlierRemoval'.format(i,plot_type)\n",
    "        plt.savefig(filename)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}