ANALYSIS_atlas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=top></a>\n",
    "\n",
    "# Multi-Channel Atlas Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "\n",
    "----\n",
    "\n",
    "1. [Preparations](#prep)\n",
    "2. [PCA & tSNE](#pca)\n",
    "3. [Covariate bi-graph](#graph)\n",
    "4. [Tissue Consensus Map](#map)\n",
    "5. [Archetype Overlay](#archo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=prep></a>\n",
    "\n",
    "## 1. Preparations\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import modules\n",
    "\n",
    "# External, general\n",
    "from __future__ import division\n",
    "import os, sys, pickle\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# External, specific\n",
    "import pandas as pd\n",
    "import ipywidgets as widgets\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.manifold.t_sne import TSNE\n",
    "import networkx as nx\n",
    "import scipy.stats as stats\n",
    "\n",
    "# Internal\n",
    "import katachi.utilities.loading as ld\n",
    "import katachi.utilities.plotting as kp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load basic data\n",
    "\n",
    "# Prep loader\n",
    "loader = ld.DataLoaderIDR()\n",
    "loader.find_imports(r\"data/experimentA/extracted_measurements/\", recurse=True, verbose=True)\n",
    "\n",
    "# Import shape spaces\n",
    "fspace_TFOR_pca, prim_IDs, fspace_idx = loader.load_dataset(\"shape_TFOR_pca_measured.tsv\")\n",
    "fspace_CFOR_pca, _, _ = loader.load_dataset(\"shape_CFOR_pca_measured.tsv\", IDs=prim_IDs)\n",
    "print \"Imported TFOR shape space of shape:\", fspace_TFOR_pca.shape\n",
    "print \"Imported CFOR shape space of shape:\", fspace_CFOR_pca.shape\n",
    "\n",
    "# Import TFOR centroid locations\n",
    "centroids = loader.load_dataset(\"_other_measurements.tsv\", IDs=prim_IDs)[0][:,3:6][:,::-1]\n",
    "print \"Imported TFOR centroids of shape:\", centroids.shape\n",
    "\n",
    "# Import & standardize engineered features\n",
    "covar_df, _, _ = loader.load_dataset(\"_other_measurements.tsv\", IDs=prim_IDs, force_df=True)\n",
    "del covar_df['Centroids RAW X']; del covar_df['Centroids RAW Y']; del covar_df['Centroids RAW Z']\n",
    "covar_names = list(covar_df.columns)\n",
    "covar_df_z = (covar_df - covar_df.mean()) / covar_df.std()\n",
    "print \"Imported engineered features of shape:\", covar_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load atlas data\n",
    "\n",
    "# Set embedding type for atlas data\n",
    "embedding_type = 'TFOR'\n",
    "#embedding_type = 'CFOR'\n",
    "\n",
    "# Specify channels to load\n",
    "channels = ['CDMPRtagRFPt', 'mKate2GM130', 'tagRFPtUtrCH', 'b4galT1tagRFPt', \n",
    "            'NLStdTomato', 'mKate2rab11', 'mKate2rab5', 'lysotrackerdeepred']\n",
    "\n",
    "# Load channel data\n",
    "obs_fspace = {}\n",
    "atlas_fspace = {}\n",
    "for channel in channels:\n",
    "    \n",
    "    # Load observed channel data (and map to full space)\n",
    "    obs, obs_IDs, obs_idx = loader.load_dataset(channel+'_'+embedding_type+'_'+'pca_measured.tsv', IDs=prim_IDs)\n",
    "    obs_full = np.empty((fspace_TFOR_pca.shape[0], obs.shape[1]))\n",
    "    obs_full[:] = np.nan\n",
    "    obs_full[np.in1d(prim_IDs, obs_IDs)[fspace_idx]] = obs\n",
    "    obs_fspace[channel] = obs_full\n",
    "    \n",
    "    # Load predicted channel data\n",
    "    atl, _, _ = loader.load_dataset(channel+'_'+embedding_type+'_'+'pca_predicted.tsv', IDs=prim_IDs)\n",
    "    atlas_fspace[channel] = atl\n",
    "    \n",
    "    # Report\n",
    "    print channel+\":\\tobs.shape is\", obs.shape, \"-- obs_full.shape is\", obs_full.shape, \"-- atl.shape is\", atl.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load archetype data\n",
    "\n",
    "# Source fspace\n",
    "#archetype_type = 'TFOR'\n",
    "archetype_type = 'CFOR'\n",
    "\n",
    "# Load archetype predictions\n",
    "archetype_classes, _, _ = loader.load_dataset(\"_archetype_\"+archetype_type+\"_classifications.tsv\", \n",
    "                                              IDs=prim_IDs)\n",
    "print \"Imported archetype classifications of shape:\", archetype_classes.shape\n",
    "archetype_probas, _, _  = loader.load_dataset(\"_archetype_\"+archetype_type+\"_probabilities.tsv\", \n",
    "                                              IDs=prim_IDs)\n",
    "print \"Imported archetype probabilities of shape:\", archetype_probas.shape\n",
    "\n",
    "# PCA of archetype space\n",
    "archetype_pca = PCA().fit_transform(archetype_probas)\n",
    "\n",
    "# Archetype annotation and visualization\n",
    "archetype_decodedict = {0 : 'unclassified', 1 : 'central', 2 : 'peri',\n",
    "                        3 : 'inter', 4 : 'leader'}\n",
    "archetype_encodedict = {name:key for key,name in archetype_decodedict.iteritems()}\n",
    "archetype_colors = {0 : 'lightgray', 1 : 'royalblue', 2 : 'limegreen',\n",
    "                    3 : 'cyan', 4 : 'orangered'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=pca></a>\n",
    "\n",
    "## 3. PCA & tSNE\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot PCA\n",
    "\n",
    "# Set interactions\n",
    "from ipywidgets import interact\n",
    "@interact(channel=atlas_fspace.keys(),\n",
    "          PCx=(1, atlas_fspace[atlas_fspace.keys()[0]].shape[1], 1),\n",
    "          PCy=(1, atlas_fspace[atlas_fspace.keys()[0]].shape[1], 1),\n",
    "          highlight_obs=False)\n",
    "\n",
    "# Show \n",
    "def show_PCs(channel='mKate2GM130', PCx=1, PCy=2, highlight_obs=False): \n",
    "    \n",
    "    # Prep\n",
    "    plt.figure(figsize=(9,7))\n",
    "\n",
    "    # Plot\n",
    "    rnd = np.random.permutation(np.arange(atlas_fspace[channel].shape[0]))\n",
    "    if not highlight_obs:\n",
    "        plt.scatter(atlas_fspace[channel][rnd,PCx-1], atlas_fspace[channel][rnd,PCy-1],\n",
    "                    c=fspace_idx[rnd], cmap=plt.cm.plasma,\n",
    "                    s=10, edgecolor='', alpha=0.5)\n",
    "    else:\n",
    "        plt.scatter(atlas_fspace[channel][rnd,PCx-1], atlas_fspace[channel][rnd,PCy-1],\n",
    "                    c='gray', s=10, edgecolor='', alpha=0.5, label='predicted')\n",
    "        plt.scatter(obs_fspace[channel][rnd,PCx-1], obs_fspace[channel][rnd,PCy-1],\n",
    "                    c='darkcyan', s=10, edgecolor='', label='observed')\n",
    "\n",
    "    # Cosmetics \n",
    "    if not highlight_obs:\n",
    "        cbar = plt.colorbar()\n",
    "        cbar.set_label('sample', rotation=270, labelpad=15)\n",
    "    else:\n",
    "        plt.legend(frameon=False)\n",
    "    plt.xlabel(\"PC \"+str(PCx))\n",
    "    plt.ylabel(\"PC \"+str(PCy))\n",
    "    plt.title(channel)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Perform tSNE on data from one channel\n",
    "### Warning: Running the tSNE on the full dataset can take a while!\n",
    "\n",
    "# Settings\n",
    "channel   = 'mKate2GM130'\n",
    "num_PCs   = 10  # Number of PCs to include\n",
    "subsample = 10  # Subsampling so it doesn't take forever\n",
    "\n",
    "# Running the tSNE\n",
    "tsne = TSNE(n_components=2, random_state=42, \n",
    "            perplexity=30, learning_rate=200)\n",
    "fspace_tsne = tsne.fit_transform(atlas_fspace[channel][::subsample, :num_PCs])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Perform tSNE on data from all channels\n",
    "### Warning: Running the tSNE on the full dataset can take a while!\n",
    "\n",
    "# Settings\n",
    "num_PCs   = 20  # Number of PCs to include\n",
    "subsample = 10  # Subsampling so it doesn't take forever\n",
    "\n",
    "# Assembling the channels + new PCA\n",
    "fspace_concat = np.concatenate([atlas_fspace[ch] for ch in atlas_fspace.keys()]\n",
    "                               +[fspace_TFOR_pca, fspace_CFOR_pca], \n",
    "                               axis=1)\n",
    "print 'Concatenated shape:', fspace_concat.shape\n",
    "fspace_concat_pca = PCA(n_components=num_PCs).fit_transform(fspace_concat)\n",
    "print 'Concat. PCA shape: ', fspace_concat_pca.shape\n",
    "\n",
    "# Running the tSNE\n",
    "tsne = TSNE(n_components=2, random_state=42, \n",
    "            perplexity=30, learning_rate=200)\n",
    "fspace_tsne = tsne.fit_transform(fspace_concat_pca[::subsample, :])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Plot tSNE\n",
    "\n",
    "# Prep\n",
    "plt.figure(figsize=(4,4))\n",
    "\n",
    "# Plot\n",
    "rnd = np.random.permutation(np.arange(atlas_fspace[channel].shape[0] // subsample))\n",
    "plt.scatter(fspace_tsne[rnd,0], fspace_tsne[rnd,1],\n",
    "            c=fspace_idx[rnd], cmap=plt.cm.plasma,\n",
    "            s=4, edgecolor='', alpha=0.8)\n",
    "\n",
    "# Labels\n",
    "plt.xlabel(\"tSNE 1\", fontsize=16)\n",
    "plt.ylabel(\"tSNE 2\", fontsize=16, labelpad=-5)\n",
    "plt.title(\"tSNE\", fontsize=16)\n",
    "\n",
    "# Done\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=graph></a>\n",
    "\n",
    "## 3. Covariate Bi-Graph\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Compute covariate correlations\n",
    "\n",
    "# USER INPUT: Parameters\n",
    "channel1     = 'covars'       # Upper channel\n",
    "channel2     = 'mKate2GM130'  # Lower channel\n",
    "num_PCs      = 8              # Number of PCs to include\n",
    "corr_measure = 'pearsonr'     # Correlation measure to use (currently only pearsonr is available!)\n",
    "threshold    = 0.30           # Threshold to include a correlation as relevant\n",
    "\n",
    "# Select the datasets to compare\n",
    "if channel1 == 'covars':\n",
    "    ch1_data = covar_df_z\n",
    "elif channel1 == 'TFOR':\n",
    "    ch1_data = pd.DataFrame(fspace_TFOR_pca[:,:num_PCs])\n",
    "elif channel1 == 'CFOR':\n",
    "    ch1_data = pd.DataFrame(fspace_CFOR_pca[:,:num_PCs])\n",
    "else:\n",
    "    ch1_data = pd.DataFrame(StandardScaler().fit_transform(atlas_fspace[channel1][:,:num_PCs]))\n",
    "ch2_data = pd.DataFrame(StandardScaler().fit_transform(atlas_fspace[channel2][:,:num_PCs]))\n",
    "ch1_names = covar_names if channel1=='covars' else [str(num) for num in range(1, num_PCs+1)]\n",
    "    \n",
    "# Compute correlations\n",
    "print \"Preparing correlation data...\"\n",
    "if corr_measure == 'pearsonr':\n",
    "    dists = ch1_data.expanding(axis=1).corr(ch2_data, pairwise=True).iloc[-1, :, :].get_values()\n",
    "else:\n",
    "    raise NotImplementedError()\n",
    "\n",
    "# Generate the plot\n",
    "print \"Generating bigraph...\\n\"\n",
    "kp.covar_pc_bigraph(dists, threshold, ch1_names, \n",
    "                    height=0.6, verbose=True, show=False)\n",
    "\n",
    "# Done\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=map></a>\n",
    "\n",
    "## 4. Tissue Consensus Map\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot backmapping\n",
    "\n",
    "# PARAMETERS\n",
    "PC = 4\n",
    "interesting = 'leaders'\n",
    "save_pdf = False\n",
    "### Contour plot backmapping plot for publication\n",
    "\n",
    "# Set interactions\n",
    "@widgets.interact(channel=channels+['TFOR','CFOR'],\n",
    "                  standardized=['no','z'],\n",
    "                  PC=(1, 20, 1))\n",
    "\n",
    "# Plot\n",
    "def contour_backmap(channel='mKate2GM130',\n",
    "                    standardized='no', PC=1):  \n",
    "\n",
    "    # Settings\n",
    "    xlim = (-130, 8)\n",
    "    ylim = ( -19, 19)\n",
    "\n",
    "    # Select data\n",
    "    if channel=='TFOR':\n",
    "        plot_values = fspace_TFOR_pca[:,PC-1]\n",
    "    elif channel=='CFOR':\n",
    "        plot_values = fspace_CFOR_pca[:,PC-1]\n",
    "    else:\n",
    "        plot_values = atlas_fspace[channel][:,PC-1]\n",
    "    \n",
    "    # Standardization\n",
    "    if standardized=='z':\n",
    "        plot_values = StandardScaler().fit_transform(plot_values[:,np.newaxis]).squeeze()\n",
    "\n",
    "    # Tools for smoothing on scatter\n",
    "    from katachi.utilities.pcl_helpers import pcl_gaussian_smooth\n",
    "    from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "    # Cut off at prim contour outline\n",
    "    kernel_prim = stats.gaussian_kde(centroids[:,1:].T)\n",
    "    f_prim = kernel_prim(centroids[:,1:].T)\n",
    "    f_prim_mask = f_prim > f_prim.min() + (f_prim.max()-f_prim.min())*0.1\n",
    "    plot_values    = plot_values[f_prim_mask]\n",
    "    plot_centroids = centroids[f_prim_mask]\n",
    "\n",
    "    # Smoothen\n",
    "    pdists = squareform(pdist(plot_centroids[:,1:]))\n",
    "    plot_values = pcl_gaussian_smooth(pdists, plot_values[:,np.newaxis], sg_percentile=0.5)[:,0]\n",
    "\n",
    "    # Initialize figure\n",
    "    fig, ax = plt.subplots(1, figsize=(8, 3.25))\n",
    "\n",
    "    # Contourf plot\n",
    "    cfset = ax.tricontourf(plot_centroids[:,2], plot_centroids[:,1], plot_values, 20, \n",
    "                           cmap='plasma')\n",
    "\n",
    "    # Illustrative centroids from a single prim\n",
    "    plt.scatter(centroids[fspace_idx==prim_IDs.index(prim_IDs[0]), 2], \n",
    "                centroids[fspace_idx==prim_IDs.index(prim_IDs[0]), 1],\n",
    "                c='', alpha=0.5)\n",
    "\n",
    "    # Cosmetics\n",
    "    ax.set_xlabel('TFOR x', fontsize=16)\n",
    "    ax.set_ylabel('TFOR y', fontsize=16)\n",
    "    plt.tick_params(axis='both', which='major', labelsize=13)\n",
    "    plt.xlim(xlim); plt.ylim(ylim)\n",
    "    ax.invert_yaxis()  # To match images\n",
    "\n",
    "    # Colorbar\n",
    "    cbar = plt.colorbar(cfset, ax=ax, pad=0.01)\n",
    "    cbar.set_label(channel+' (PC'+str(PC)+')', rotation=270, labelpad=10, fontsize=16)\n",
    "    cbar.ax.tick_params(labelsize=13)\n",
    "\n",
    "    # Done\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=archov></a>\n",
    "\n",
    "## 5. Archetype Overlay\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot archetype overlay\n",
    "\n",
    "# Set interactions\n",
    "@interact(channel=atlas_fspace.keys(),\n",
    "          dim=(1, atlas_fspace[atlas_fspace.keys()[0]].shape[1], 1), \n",
    "          PCx=(1, archetype_pca.shape[1], 1),\n",
    "          PCy=(1, archetype_pca.shape[1], 1),\n",
    "          vmax_factor=(0.0, 1.0, 0.1),\n",
    "          show_classes=False)\n",
    "\n",
    "# Show \n",
    "def show_proba_fspace_overlay(channel='mKate2GM130', dim=1, PCx=1, PCy=2, \n",
    "                              vmax_factor=1.0, show_classes=False): \n",
    "    \n",
    "    # Prep plot\n",
    "    plt.figure(figsize=(9,7))\n",
    "\n",
    "    # Grab overlay data\n",
    "    overlay_data = atlas_fspace[channel][:, dim-1]\n",
    "    \n",
    "    # Create scatter\n",
    "    if not show_classes:\n",
    "        scat = plt.scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                           color=overlay_data, cmap='viridis',\n",
    "                           vmax=vmax_factor*np.max(overlay_data),\n",
    "                           edgecolor='', s=10, alpha=0.7)\n",
    "        \n",
    "    # Show the archetypes instead (for reference)\n",
    "    if show_classes:\n",
    "        for key in archetype_decodedict.keys():\n",
    "            mask = archetype_classes==key\n",
    "            if np.any(mask):\n",
    "                scat = plt.scatter(archetype_pca[mask, PCx-1], archetype_pca[mask, PCy-1],\n",
    "                                   color=archetype_colors[key], edgecolor='', \n",
    "                                   s=10, alpha=0.5, label=archetype_decodedict[key])\n",
    "        plt.legend(frameon=False, fontsize=16)\n",
    "\n",
    "    # Cosmetics  \n",
    "    plt.xlabel(\"PC \"+str(PCx), fontsize=18)\n",
    "    plt.ylabel(\"PC \"+str(PCy), fontsize=18)\n",
    "    if not show_classes: plt.title(channel + ', dim '+str(dim), fontsize=20)\n",
    "    plt.gca().tick_params(axis='both', which='major', labelsize=16)\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    # Show\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Archetype overlay panel\n",
    "\n",
    "# Prep\n",
    "fig, ax = plt.subplots(1, 4, figsize=(12,3), sharex=True, sharey=True)\n",
    "PCx = 1; PCy=2\n",
    "\n",
    "# Plot archetypes\n",
    "for key in archetype_decodedict.keys()[::-1]:\n",
    "    mask = archetype_classes==key\n",
    "    if np.any(mask):\n",
    "        scat = ax[0].scatter(archetype_pca[mask, PCx-1], archetype_pca[mask, PCy-1],\n",
    "                             color=archetype_colors[key], edgecolor='', s=7, alpha=0.5, \n",
    "                             label=archetype_decodedict[key].replace('Cells','').replace('Rosette',''))\n",
    "        \n",
    "# Add archetype legend\n",
    "legend = ax[0].legend(frameon=False, fontsize=9)\n",
    "for i in range(len(legend.legendHandles)):\n",
    "    legend.legendHandles[i]._sizes = [30 for s in legend.legendHandles[i]._sizes]\n",
    "\n",
    "# Plot GM130 TFOR PC 3\n",
    "channel = 'mKate2GM130'\n",
    "dim = 3\n",
    "vmax_factor=1.0\n",
    "overlay_data = atlas_fspace[channel][:, dim-1]\n",
    "scat = ax[1].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                     color=overlay_data, cmap='viridis',\n",
    "                     vmax=vmax_factor*np.max(overlay_data),\n",
    "                     edgecolor='', s=7, alpha=0.7)\n",
    "\n",
    "# Plot UtrCH TFOR PC 3\n",
    "channel = 'tagRFPtUtrCH'\n",
    "dim = 3\n",
    "vmax_factor=1.0\n",
    "overlay_data = atlas_fspace[channel][:, dim-1]\n",
    "scat = ax[2].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                     color=overlay_data, cmap='viridis',\n",
    "                     vmax=vmax_factor*np.max(overlay_data),\n",
    "                     edgecolor='', s=7, alpha=0.7)\n",
    "\n",
    "# Plot NLS TFOR PC 3\n",
    "channel = 'NLStdTomato'\n",
    "dim = 3\n",
    "vmax_factor=1.0\n",
    "overlay_data = atlas_fspace[channel][:, dim-1]\n",
    "scat = ax[3].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                     color=overlay_data, cmap='viridis',\n",
    "                     vmax=vmax_factor*np.max(overlay_data),\n",
    "                     edgecolor='', s=7, alpha=0.7)\n",
    "\n",
    "  \n",
    "# Cosmetics\n",
    "for axx in ax.flatten():\n",
    "    axx.set_xticks(np.arange(-0.5, 1.1, 0.5))\n",
    "    axx.set_yticks(np.arange(-0.5, 1.1, 0.5))\n",
    "    axx.set_xlim([-0.7, 0.9]); axx.set_ylim([-0.5, 1.0])\n",
    "    axx.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "# Show\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Archetype overlay panel for publication\n",
    "\n",
    "# Prep\n",
    "fig, ax = plt.subplots(2, 2, figsize=(12,12), sharex=True, sharey=True)\n",
    "PCx = 1; PCy=2\n",
    "\n",
    "# Plot TFOR-PC1\n",
    "dim = 1\n",
    "vmax_factor=1.0\n",
    "scat = ax[0,0].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                       color=fspace_TFOR_pca[:, dim-1], cmap='viridis',\n",
    "                       vmax=vmax_factor*np.max(fspace_TFOR_pca[:, dim-1]),\n",
    "                       edgecolor='', s=15, alpha=0.7)\n",
    "\n",
    "# Plot TFOR-PC3\n",
    "dim = 3\n",
    "vmax_factor=1.0\n",
    "scat = ax[0,1].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                       color=fspace_TFOR_pca[:, dim-1], cmap='viridis',\n",
    "                       vmax=vmax_factor*np.max(fspace_TFOR_pca[:, dim-1]),\n",
    "                       edgecolor='', s=15, alpha=0.7)\n",
    "\n",
    "# Plot UtrCH TFOR-PC3\n",
    "channel = 'tagRFPtUtrCH'\n",
    "dim = 3\n",
    "vmax_factor=1.0\n",
    "overlay_data = atlas_fspace[channel][:, dim-1]\n",
    "scat = ax[1,0].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                       color=overlay_data, cmap='viridis',\n",
    "                       vmax=vmax_factor*np.max(overlay_data),\n",
    "                       edgecolor='', s=15, alpha=0.7)\n",
    "\n",
    "# Plot CFOR-PC1\n",
    "dim = 1\n",
    "vmax_factor=1.0\n",
    "scat = ax[1,1].scatter(archetype_pca[:, PCx-1], archetype_pca[:, PCy-1],\n",
    "                       color=fspace_CFOR_pca[:, dim-1], cmap='viridis',\n",
    "                       vmax=vmax_factor*np.max(fspace_CFOR_pca[:, dim-1]),\n",
    "                       edgecolor='', s=15, alpha=0.7)\n",
    "\n",
    "# Cosmetics\n",
    "for axx in ax.flatten():\n",
    "    axx.set_xticks(np.arange(-0.5, 1.1, 0.5))\n",
    "    axx.set_yticks(np.arange(-0.5, 1.1, 0.5))\n",
    "    axx.set_xlim([-0.7, 0.9]); axx.set_ylim([-0.5, 1.0])\n",
    "    axx.tick_params(axis='both', which='major', labelsize=24)\n",
    "plt.tight_layout()\n",
    "\n",
    "# Show\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "[back to top](#top)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}