exploration/Notes_on_PCA_with_prcomp_and_factoextra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.utilities import *\n",
    "%run helpers/notebook_setup.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%R library(factoextra);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While analyzing the PCA plots for protein levels two problems have arisen:\n",
    "- the loadings plot did not show any variables that define PC2\n",
    "  - in order to select variables that explain both PC1 and PC2, a way to extract the contributions of variables of each PC (but not of joint PC1 and PC2) was needed\n",
    "- it was not obvious how are the contributions defined, given different terminologies in use\n",
    "\n",
    "This notebook addresses these questions using the protein levels data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "path = 'data/clean/protein/levels.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_levels = read_csv(path, index_col=2).drop(['soma_id', 'target_full_name', 'entrez_gene_symbol'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def z_score(v):\n",
    "    return (v - v.mean()) / v.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import log10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_matrix = protein_levels.applymap(log10)\n",
    "z_log_matrix = log_matrix.apply(z_score).apply(z_score, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i z_log_matrix -o pca\n",
    "pca = prcomp(t(z_log_matrix), scale=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A quick look at the distribution of contributions to PC1/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "pc_contributions_distribution = function(pca_result, axes, choice='var', labels=FALSE) {\n",
    "    g = fviz_pca_contrib(pca, choice=choice, axes=axes)\n",
    "    if(!labels)\n",
    "        g = g + theme(\n",
    "            axis.text.x=element_blank(),\n",
    "            axis.ticks.x=element_blank()\n",
    "        )\n",
    "    g\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAEsCAIAAAD2HxkiAAAgAElEQVR4nO3deVwT19oH8DMJKERQWQSJgBRxQ6+WWi3U3loVelu13lZ9q7d6xe2qVWttxavWYktrX7Sor1KuW3HHrSougF7rUktFAdmhCkSRfU1YQ0hCMuf9YzClYdhiwsnyfP/oZzjMmXkm9GdmPxTGGAEAyOGQLgAAUwchBIAwCCEAhEEIASCMfAizs7M//PBDV1dXOzu76dOnZ2Zman0VMTEx48eP1/piNXPy5ElXV9eIiAhdr4h1q0tKSszNzXW9apW//OUv1HNDhgxZtWpVY2Njj63dUBAOYUJCgre395AhQ6Kion755RdnZ+cpU6aUlpZ2dzkuLi7JycltF05RVE5OjmaL0qx7p3bv3h0UFLRgwYIXX5SOKuwY60fdgV27dpWXl5eUlBw8ePDhw4dTpkxhTsi/YPFHjhz55JNPNOurbwiHcP369StXrgwODh47duyYMWMOHjzo7u7+448/vviSq6qq+Hz+sWPHnJycNFvCkCFDXqR7e6RSqbu7+4sv58U3sGf069fP0dGRz+f7+vrevn07Jyfn0qVL6AU+3rS0tK+//nrjxo06KJYQTE5ubi5CqKysrHVjYWGhQCBgph88eODj42NlZTVq1Khz585hjIVCYf/+/W/evDlr1qxBgwbNnz9fLpdPmjSJw+HY2dldvnxZKBTa29vv2bOHy+VeuXIFIaRQKKKjo8eNG/f11187OTm5uLh89dVXSqUyOzvbyspKtV5fX9+IiAi1RTHdu16J2ga27TVz5kwul2tra3vq1KnWc06fPn3Lli3MdGNjI4/HS0hIwBhHRUWNGjWqd+/egwYNCg4OZtbbdgNZ52Td6uLiYjMzM9Xn7+fnZ21tzfzDxzQGBwc7OzvzeLzJkyfn5eWpbVHrz4d1A9WMHj06PDy8dcuiRYuWL1/ObAhTPPNJRkRETJ8+3cXFJTw8PDQ01M/Pb8CAAfv27Wu7zPDw8JUrV3p6eq5Zs6btbw0RyRDGxMTY2Ni099uqqqp+/frt3bu3srLy8uXLffr0efjwoVAoNDc3DwwMxBg3NjY6Ojoyf3tnZ+ekpCSMsVAoNDMzmzdvXnFxserPHB0dzeFw1qxZU15efufOHUdHxwMHDrCGUG1RTPduVdJx/RjjESNG3L17V21jIyIihg0bxkyfPXvW09MTYyyTyfr06bN3797y8vKYmBgulysQCFg3kHVO1q1WhVAqlQ4ePHjbtm2VlZU3btywsbGJjo6+c+eOlZVVXFxcQUHBtGnT5s+f3/bvovp82tvA1tqGMCgoyM/PD/85hBwOZ/v27RjjS5cuURS1bds2Zrpv377t/e+xcuVKCKEWREZGOjg4tPfb8PBwb29v1Y8rVqz49NNPmb9ccXEx0zhr1qywsDDcJjnPnj3Drf7M0dHRffv2lUqlTK/du3dPmjSp6yHsViUd14/bCWFDQ4OlpWVmZibG+P333w8JCcEYy2SytLQ0ZoZHjx5ZWVnFx8ezbiDrnKxbrQrhpUuXxo4dqypg69ata9asiY6Otra2/vnnn5llNjQ0tP27qD6f9jawtS6GECFUWlqKMa6srFSbZr7n2zKmEJI8JvTw8KisrKytrW3dePXq1dWrVyOEiouLPTw8Ws9cWFjITDs6OjITHA57/YMGDVJrcXFx6d27NzM9YsSItud+aJpur07NKumgV1tWVlYzZsy4cOFCXV3dzZs3mdM2vXr1SkhImDp16iuvvBIYGNh6FWob2N6cHWx1Xl5edna203OHDh2qra2dNm3azp07N23aZGNj4+/vzxwvaPCxdCA/P/+ll15q225ra4uef4ytpxFChw8ftra2tra2fu+99zpdviEiGUJPT8/hw4f/8MMPrRv3798vlUoRQoMGDXr69KmqPS8vz9nZmZmmKKrjJbeNRFFRkUwmY6YfP37s5uaGEFIqlfj5rbMFBQXtLU2zSjroxeof//jHxYsXL1++PGXKlIEDByKEbty4sWXLltDQ0JSUlAsXLlhaWra3ge3NybrVDCcnJy8vr7LnMjIydu7cWVRUNHny5OTk5NzcXHd396lTp3bwb1N3NxAhJBaLL1265Ofn1/FsapYuXdrQ0NDQ0BAVFdWtjoaCZAi5XO6ePXuCg4O3bt2alZWVm5u7efPmX3/99d///jdCaObMmb///vu+fftqa2uvXbt28uTJ+fPnd7C0urq6Dn7b2Nj4+eefV1VVxcbGfv/990uXLrW1tW1qajpx4oRcLv/hhx9KSkraW1R3K9Gs17vvvltYWLh9+/bFixczLZWVlTwez8HBobKyMjg4uKKioqysjLVve3O23WpVl3feeefJkycHDhyoq6uLj48fM2ZMbGzs9evXvb2909PTzc3NrayseDwe674G8/l0cQMbGhpEIlFlZWVsbKyfn9+wYcNmz57d6adnWkjvD+PExMRp06Y5Ojra2dn5+fklJiaqfhUXF/faa6/16dPH09Pz7NmzuNWBBDPDnDlzmCOxgIAAa2vryMjI1jOopmNjY9euXRsQEODg4ODm5rZr1y6me0hIiL29/aBBgzZv3rx58+aYmJj2FtX1Slpr2wu3c0zI8Pf3HzBggOosq1QqnTt3rrW19ZAhQ3bs2LFz585hw4YVFRW13UDWOW/cuNF2q1ufHU1KSnrjjTf69Onz0ksvhYaGYoxlMtnSpUvt7Ox4PJ63t/f9+/fbFqn6fNrbwNZGjx6t+j/Nzc1txYoVquNMtWNC5ti17bQpHBNS2JAfZZLL5QqFgsfjkS5EjyoxJjU1NTY2NqSr0DnDDiEARoD8vaMAmDgIIQCEQQgBIAxCCABhEEIACIMQAkCYGcF1Nzc30zRNUeqXSZh7wTpt7OJsTKPerqJbfY1jFcgY/45dX4WZmRmXy/1TCyKn4xvNADBW9vb2rX+E3VEACIMQAkAYhBAAwiCEABAGIQSAMAghAIRBCAEgDEIIAGEQQgAIgxACQBjJEAoqxATXDoCeIBnCyIfFBNcOgJ4gGcKs4vqMIriHG5g6kiGkMT59v4hgAQDoA8InZoprJM8qYehWYNI0D6FCodi9e/c333xz7NgxtV9t27aNGU9CKBT6+/tv2LBhw4YNrOPvYoz/N+qxxjUAYAQ0f6g3Pj6ez+fPmzcvODi4uLiYGQxELBYHBQWpxkCuqqqaNm3a3LlzO1hOU7OyViLvz+ulcSUAGDTNvwkFAgEz7LO7u7tAIGAaraysduzYMWbMGObHioqK0tLSsLCwu3fvtrccjFHQpWyNywDA0GkeQolEwjylb2dnJxb/ccWPw+GoBgzj8XijR4+eN2/e7du3MzIymMZVq1ZRFEVRVKOw5RJFZX3T9mjIITBRmoeQx+OJRCKEkEgksrKyYp1nwoQJfn5+9vb2U6ZMUY04uXPnzurq6urqap4tXzVnZlHd+tMZMCwGMEGah3Do0KH5+fkIocLCwtYjtrZ25syZ9PR0Zh4nJyemkcfj2djY2NjYUH8e+668rumLc5ka1wOAgdI8hN7e3iUlJSEhIQ4ODi4uLrm5uWFhYWrz+Pr6nj17NjAwsKamxsfHp9NlFtZIfnlUpXFJABgikkOj/e3bGNb2HfPGOttY9HAxAPQYA3jl4deRmTBqIjAd+hjCJjm94GBiwJl0iVxBuhYAdE4fQ4gQQhiV1UqXH07eHp2tUNKkqwFAh/Q1hAghhDBCmUV1S8KTKuplpGsBQFf0OoQMJY3Xn0pLLqgmXQgAOmEAIUQIYYR2XxccvP2UdCEAaJ9hhBAhhDCKzRUuPJhY2SAnXQoA2mQ4IUQIIaSk8eenUivhEBEYEQMLIUIIY/TZqbTQGwLShQCgHYYXQkZCXvWW85lwRR8YAUMNIUIoXyhZeSSJhntrgIEz4BAihMQy5b8OJzXBjTXAkBl2CBFC0mb6X4dTCoQS0oUAoCGDDyFCCCO85XzW0dgC0oUAoAnNX/SkVzDCt34vv/Wo3KGv5efvDnWxsSRdEQBdZQzfhH/AqLKuafPZjHOJJaRLAaCrjCuECCGEMEJXk4uXH00ur5eSrgWAzhlhCBmNUsX6U+kbzmY0NStJ1wJAR4w2hIzSmqZl4Ulrj6dWNsCdbkBPGcmJmY6JJPLPItLsrXv5v+H2ipsN6XIA+BOTCCFD2CDfdT23Ty/uZ9OGjXTqS7ocAFoY+e5oW41y5bbLj7ddhVFogL4wuRAyHpfUR9yHi/tAL5hoCBFC19PLv/gJXvgNyDPdECKECkSStRGpzTQ8hwFIMukQIoREDfIlhxJ/zYF37wNiTD2ECCEao0N38uYfTDx+Lx+eTgQ9D0L4HI1/zqzwP/jwv+nlpEsBpgVC+Cc0xifvFyz+8WFelbjzuQHQBgghC7mCDrz46Ni9Z6QLASYBQtgOjG9mVq44lqyAc6dAxyCEHRE3KdafSSddBTByEMJOCOtlq08kk64CGDMIYedqGxWLfkysaWomXQgwThDCLmlW4LXHU76+9LhRBq9XBFoGIewqGiNBef3Hx1IE5fWkawFGBULYPUoaf33pcTMNgwcDrYEQamLFkeS8ykbSVQAjASHUhKyZDryYFXgxi3QhwBhACDWXV9m48EDij7/AjTXghUAIX4gS47vZlfCVCF4EhFAL8iobr6XBsxdAQxBC7TgdX5AnhFM1QBMQQu3AGH114Xe4lA80ACHUGhrjLRfg4BB0G4RQm6rqZT/cfALPPoFugRBqWfwT0fIjyZBD0HUQQu2TyBQfH02GG9tAF0EIdaJBqth8LoN0FcAwQAh1pbi6aemPSTIFjI4IOgEh1CGpQvm/kY/lCtgzBR2BEOrWE1HjkvCH+XAdH7QPQqhzGKOtF7KS82pIFwL0FISwJygx2n0j94dbTzC8Zh+0ASHsOfECkf+hh5FJpaQLAfoFQtijlDS++LBo4cGH0SkQRdACQkiAkqbPJBSFxOSQLgToBQghMWmFtV9F/k66CkAehJCkJxXiNSfTlDDchWmDEBJWI5Z9ciK1tKaJdCGAGAgheXVNzduj4PjQdEEI9YKoUfbvs3DDt4mCEOqLkpqm3f/NJV0FIABCqEeSn9WEXIccmhwIoX5Jy6/554GE7LIG0oWAngMh1Ds0Rt9efrQVLiGaDAihnnpaIc4tgwegTAKEUH+diS8gXQLoCWZtm2iajoqKSkhISEtLo2nay8trwoQJf//73zkcSGyPEpQ3lNVKnfpbkC4E6JZ6rq5fvz5z5sz09PTXXntt9+7de/fu9fHxycrKmjlz5vXr14mUaLIwQhvPpacV1pIuBOiW+jehUCi8dOmSubm5qmX48OEzZ85UKBRnzpzp2doAUtIoJCbn4JJxVr1Z9lmAcaA6ftZboVDcvn178uTJvXr1avur0NBQsVjs6uq6aNGi1r/atm1bQECAhYVFB/MghP72bYwWtsAEcBA1yrXfxunDKdKVAK2wt7dv/WMnh3lz5sxJSUmZO3du21/Fx8fz+fytW7eWlZUVFxczjWKxeMOGDYmJiR3MA7qLRjizsPbAnTzShQCdUA/h48ePW/9YVFS0efPm0lKWx8AFAoG7uztCyN3dXSAQMI1WVlY7duwYM2ZMB/OkpKScP3/+/PnzChk8OtAN93Kr6qUw6pMRUg/h/v37V69erUrd+PHjly1b5unp2banRCJhvlXt7OzEYvEfS+RwKIrqYJ579+4dOnTo0KFDChlcB+sOjD49mVYnaSZdB9Ay9RCGhoauW7du06ZNgYGB9fX1Bw4cWLdu3ZEjR9r25PF4IpEIISQSiaysrFiXzjrP2rVrb968efPmTYu+9qy9QHvkCuX3MdkyeJuwcWE5Jhw6dOjx48dnzJixePHi0NDQYcOGqb7Z1GbLz89HCBUWFnp4eLAuvSvzgG7JF0o+PpZCugqgTSwhPHDgwOuvv56amnrx4sWXXnppzpw5Z8+ebTubt7d3SUlJSEiIg4ODi4tLbm5uWFhYx/PoZAtMj6xZOX9/whcwIKmxYLlE4eXllZqaOmHCBOYkp1KpPHbs2NKlS7W+brhE8YI8B/XdMnMk6SpAt3V+icLS0jIuLq53797Mj1wuVxcJBC/uUUl9Vb2UdBXgRbGE8Pjx45GRkawnY4C+Of2giHQJ4EWp7442NTVZWlqyztrBrzQDu6NaQVHoH96u0192Il0I6KpOdkfXrVsXFBRUUVHRurGioiIoKOjTTz/VeXWg+zBGZ+ILqxvlpAsBGmI5MRMXF/fdd9/l5+czeRUKhW5ublu2bJk4caJ21w3fhFpk3dvswJJxpKsAXaL2TdjRDdzM3Z7Ozs46KgVCqF2jnftufg9OlhqAbtzA7ezsrLsEAq3LKq6/LxCSrgJ0Gzwsb1SOxuYrYGQLQwMhNCoSufKHm09IVwG6p/MQwgjPhiW9sIZ0CaB72EN4+/btdevWIYQ++OCDvn37Hjt2rEeLAi+gWYHrJHC5wpCwh3DFihUzZsx48OCBVCrNycn57rvvergs8CJWH0+NSikhXQXoKvYQymQyX1/fmJiY999/n8/nNzfDg6SGBCN0LqH4XDzc0WYY2N/hNXbsWH9//9jY2NTU1K+++srJCW6JMjAYoatppZNHOThY9yZdC+gE+zfhqVOnvL29L1++3L9/fy6XCy87NEgYBZ6HAS0MAHsIuVwun8/Pzs4+d+7c8OHDk5OTe7gsoBViWfP2qMedzweIYt8dXbJkiUwmY16UhhDicDizZ8/uwaqA1mQW14dcy94wbQTpQkC72EOoUCiuXLnSw6UAHUkrqDt4++mKqUNIFwLYse+Oenh41NbCEAjGIzZXuOTHxN+L60gXAliwfxOWlJTw+XwfHx9bW1um5fz58z1YFdA+mQLvui44vGwc67vzAEHsjzLFxcWptWj9YUIEjzKR0J9n/h//V0hXYeq69CjTxIkTxWJxdHT0lStX6urqdJFAQEStpPmHm09JVwH+hD2EQUFBW7du5fP5rq6u33zzTVBQUA+XBXQn/okwqwgODvUI++6oh4dHRkYGj8dDCEkkEi8vr5ycHK2vG3ZHSeltzj2y7FXSVZiurj5Zrzp8h1GyjY+sWbnqeAo8oqYn2M+Ozp8/f+rUqQsXLqQo6uTJk6zjEwKDVidpPvlb/sK/upEuBLT/oqdr167dvXsXIfTmm2/OmDFDF+uG3VHivv5g1NCB7CNqAd3p5G1r48ePv3DhQkBAgFo3XVwnhBASZ8blHFs+Hq4b9jC1EKrvjgYEBNja2jKP1QOjp1DS0all73nBo2okqYeQOfw7ffr0f/7zH1XjokWL4FKhsTqbUGRlwZ080oF0IaZLPYQjR45ECBUWFt65c4dpUSgU/fv37+m6QI/BOPzus/Ja6T98XEmXYqLUQ3jv3j2E0KpVq/bt26dq1FEIZ8dd6i9puWp8f4T3Y5eWx20GVxb6prf8E9DMNTsxZYGqy3uJMQ51Vcx0yhCvVPexzLRjbcWMh9dVsx2b+k8lh8tMv516y0VYzEz/7uoZP3wCM20jrp314LKqy+k35zb1bhnu5s2s34aWtdxW8sRpyK+j/8pMW8qaPoo9p+oS6fN+jVXLJ+Odkziq8BEzXWTv/LOXLzPNoenFt0+ousS8+m65jSMz/XJe+rinqcx0Zb8BUROmq2ZbeCfCXKlgpm+NnVLg0BKPkUXZr2fHM9N1vL4XJs5SdZn723kraSMz/ZvnxNxBQ5lp9/JnkzN/Zaal5han3pqn6vJ+/FW7hmpmOqloPPZZwhwccvPyLE6dapmJohq//FLVxSIigvvsGTPdPHGifMqUls0sLbU8fFg1m2TDBmxhwUz3vnDB7HHLM43NXl7y5yf5qJoaXqtRZSWrV+PnNyr3iokxT2kZjVgxcqRszpyWLjIZ7/vvVV2ali6l+fyWLr/8Yn7vHjOtdHOT/vOfqtn6fPcdolsGGJd+9JFySMvTJOYPHvS6dYuZph0dm5YvV3Xh7dpFSSTMtOyDDxSjRzPTZmlpvaOiWrpYWze1Omqz3LePIxK1dHnnHcX48cw0NyfH4qefmGlsZibZvBm1oR5COzs7hNDChQvj4+P/2Iw+fd566622nYHRoJX4xG/5/nDFggT2SxSLFy9mJioqKu7evbt27drt27drfd1wdlTfLP6rm+9oR9JVGL9uDAjDqKio+PjjjyMjI7VeCoRQD/Uy4+xf9IqFOZd0IcasGwPCMAYMGKCLG0eBfpIr6F3XcklXYVrYb1tT7Y4ihLKysl5++eWeqgeQ96i0/npG+btjBpIuxFSwh3DO8/NRCKEFCxbARUJTE3G/4NWXbAdY9yJdiElg3x2dOnVqTk5OeHj4iRMnnj17ZmbGnlVgtDD68mImPGbRM9odi+LOnTvLli3z9/e/cuXK+vXre7gsQJy4SfHt5UekqzAJ7GdHXV1dnz59am5ujhBSKBQjRox48kT7o97B2VH9N/vVQbPGw2jNWtals6ODBw+uq2u5l6WmpkatDzAdkcmljTIF6SqMnPrBHvM6mcGDB7/yyiuzZs1CCF28eHHhwoUESgN6AGN8LqF4yZtupAsxZuohHDhwIPPfv/615YbJL7/8Em7gNmXphTUIuZGuwpiphzA8PJz1oV54w4XJqhY355TWD+f3JV2I0YKHekEnaIz33nyyD14ZrDMsZ0erq6sPHTq0adMmXa8bzo4akH6WvfYsGNPLDO4p1YLOz47269cvIiLi2fPHxgBACNU1yQNOZ5Cuwjix3ArD5XJHjRrl5eX1+uuv9+nTh2mEAWGAqFG+83pOwLvDSRdibGBAGNA9018e+JHPYNJVGLZO3rbGgBc9gfbEpFUUVDVuem8kDLGmLfCiJ9BdOKukIT6v2meIHelKjATJFz0BwxX281O7D3oPg7d3awP7MaFYLL59+7ZUKmV+NDMzmz17ttbXDceEBo1DUceWj+dyYKe027p0TLhkyRKZTObu7s78yOFwdBFCYNBojEN/Fnz2zjDShRg89hAqFIorV670cCnA4CQ9q9n/S97Hk91JF2LY2B9l8vDwqK2t7eFSgCG6l131U2IR6SoMG/s3YUlJCZ/P9/HxsX3+UmS4WA/acyW5lOJQ//MqPPurIbhYD7SD19ts/6JXzOA8TRd06cn6iRMnisXi6OjoK1eu1NXVwZV60CmJTPH5qTTSVRgk9hAGBQVt3bqVz+e7urp+8803zOP2AHRMJJb/33/hxcHdxr476uHhkZGRwePxEEISicTLy0sXL+GG3VGjZNen166PxpjDQ0/t6+pr8FV3BnI4nb8qHwAVUaP8X+EpeRVi0oUYDPazo/Pnz586derChQspijp58iS82wJ0SzOmAyN/X/bWkMkj4T19nWt3VKZr167dvXsXIfTmm2/OeD6wo3bB7qhxoxAK/MBz+EBr0oXonc6HRqNpury8nM/nI4SioqImT55sZaWT+3QhhEbPysLs4OJxpKvQO50cE+bm5o4cOfLgwYPMj3v27PH09MzIgPcaAE2IpYrMonrSVeg79W/Ct99+e9asWStXrlS1/PTTT4cOHbr1fHRvLYJvQlNAUSho9qghA+Chpz90sjtqZ2dXUFCgtv/p4uJSVKT9+wMhhCbCxd5y+/+MIV2FHulkd3TAgAFq1ySamposLCx0XhcwXsXCppiMMtJV6C/1SxTvvvvuF198sXv3biaKGOPAwEBfX9+2PRUKRWhoqFgsdnV1XbRoEWujUChcv369g4MDQuizzz5jTvYAU4MRikwonj7GiXQheko9hMHBwR9++OGIESPeeOMNLpf74MEDOzu76Ojotj3j4+P5fP68efOCg4OLi4udnZ3bNjY0NEybNg0uMwKpgj4RV7RwogvpQvSReggtLCyuXr2anJycnJwsl8sXLFgwadIk1p4CgWDUqFEIIXd3d4FAwIRQrZGiqNLS0rCwsNGjR7/11lu63RSg325klEpk8pVThpAuRO+w3zEzbty4ceM6ubwjkUiY40s7OzuxWMza6OjoOHr0aC8vr71799ra2o4ZMwYhdPr06V9//RUhJO//Zi8reIWUCYnLFS2f7M6BdyX+meb3hfJ4PJFIhBASiUSqs6lqjRMmTPDz87O3t58yZUpubsv99fb29u7u7u7u7hQX7vE1LTTGB27lka5C72gewqFDh+bn5yOECgsLPTw8WBvPnDmTnp7O/Ojk1HJc/vbbb2/cuHHjxo3mlnBDk8m5/0QUelNAugr9onkIvb29S0pKQkJCHBwcXFxccnNzw8LC1Bp9fX3Pnj0bGBhYU1Pj4+OjxbqBgcIIJzypPptQTLoQPdLuDdw9AC7Wm7KTKyeY7MFhV58nBECncssbSJegLyCEgIwd0TnE9sH0DIQQkCFX0MIGKekq9AKEEBCTXgAvmEYIQggIOv0AXt2NEIQQECRT0Jt/yiRdBXkQQkBSoUiy8miyTEGTLoQkCCEgrEGq+OriI9JVkAQhBOQVVTfe+r2SdBXEQAiBXjga+yyzqI50FWRACIG++M+tJzS5mygJghACfdEgVSw/mmSCJ2kghECPNMnoNcdTlbRp5RBCCPSLRK5YEp5kUvulEEKgdxRKvOJICukqeg6EEOgjiVwRcDa1TtJMupCeACEEeqqsRr7qeMqKIynSZiM/RIQQAr0mljUHXshS0sZ8iAghBPqutLbp4bMa0lXoEIQQGICLiUVG/FUIIQQGoLRWevp+AekqdAVCCAzDtfTy0BtPSFehExBCYDAS8kQ/PTTCF5ZCCIEhuZJUUttkbBcPIYTAwGw+l2lkJ2kghMDA1Dc1l9Q0ka5CmyCEwPBsPJeRYURPAEMIgQHC6Fy88bwuEUIIDFK+sPHYb0Zy5RBCCAzVzayKOomcdBVaACEEhguvPpFaWmvwA1pACIEBwxhtOJP+f9cNe+hfCCEweEn51dEZ5aSr0ByEEBiDM3EFRdWGevEQQgiMxDeXHxEc+/1FQAiBkZDIFNsuPyZdhSYghMB4ZJc3nEswvMcsIITAqFxNKYlKKVHEs4YAAAR2SURBVCVdRfdACIGxiUwqbmpWkq6iGyCEwNjIlfgng7qzFEIIjNDPWRWfn0knXUVXQQiBcaqolVbUGcYdbRBCYLS+i8omXUKXQAiB0RI1yMJuGcBtpRBCYMweCKqziupJV9EJCCEwct9f1/edUgghMHJKJV525CHpKjoCIQTGr0lG33lcSbqKdkEIgUk4GptPuoR2QQiBSaBp/HNmBekq2EEIgamISS0jXQI7CCEwFcJG2Y7oHNJVsIAQAhOSUVS74nBKvVS/hpSBEALTIpY3rzqaoldDO0EIgcnBCH1yPPVqqr48+wshBKaIxvhcfNHyI0kSOfnHfyGEwHQ1ypT/Opz0yYmUBpmCYBkQQmDqqhubVx9NvpZB7AIGhBAApMToVFzhggOJjTICJ2wghAC0wBivPp4qU9A9vF4IIQB/aFbipeEP99960pMrhRAC8CcYo3tPRHECUY+tEUIIQBsYHYvN77GRLSCEALCQyBWLf3woEvfESMAQQgDYNSvx2pOpHx9NqZXo9pQphBCAjtRLmz+NSCut0+HghxBCADqhUNIbTmf8mqOrF2RACAHokkO/PPvkZFpBtUTrSzbT+hIBME4YVYtlX5zLdLKxePsvjr6ejhyK0sqCKYIjDP/t2xhSqwbgBXE4aNKIAcsmuWvQ197e/k+L0lJJAJgWmka/PKpaGp7U+MIPQ2m+O6pQKEJDQ8Visaur66JFi1gbWecBwGhIm5XLDydZ9uJ4udkunOhqbWGuwUI0/yaMj4/n8/lbt24tKysrLi5mbWSdBwAj0ySn7+cKVx5NmX8gYfGPDx8+q+lWd82/CQUCwahRoxBC7u7uAoHA2dm5bWN+fn7beQAwWhjJFfSe/+YihDhcTt/eHNs+Fnxbi6GOfUc7Ww3sZ8naSfMQSiQS5vjSzs5OLBazNrLOExAQcOTIEYTQKP/vebZOGhcAgD6jlXSthK6ViPOqxPdyhKp2Loe6tmVa6zk1DyGPxxOJRO7u7iKRyMHBgbWRdZ5NmzatWrUKIfRU1CxFmKIpTCEOjTCFEUI0QhwKIcyhMeZSCFGYxghTiMIURSFEY0RRiEIY04hDIZqDEeYihBHGCOGWvhSNUeu+HNyyCkQhjDCmEIUQwhwlps0oCiFMI0QjzEEcCiGEMUVRmEIYY8RBFE3RCHGxWnldXYUC0+YUhTBNUxSNMJeiEKZ0twolwmYUxSyIgxCiEK3RKjDzQSBMURTGHBrTZhTFfMg0xlxO11fB8ndEaqvgUPjPf0etrILS9f8qFEI029+xs1UM6sdTi5LmIRw6dGh+fv748eMLCwtff/111kZzc/O289jb2zNfj337CttdOgAmQ/MTM97e3iUlJSEhIQ4ODi4uLrm5uWFhYWqNaj9qsW4AjAbJi/VCIXwTAlMEF+sB0C8QQgAIgxACQBiEEADCIIQAEAYhBIAwCCEAhBlhCCktPe9sKOsFGiPyJ2u7UpIX6wEAyCi/CQEwLBBCAAiDEAJAGIQQAMIghAAQBiEEgDAIIQCE/T923mV/DtE9sgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 300 -h 300 -u px\n",
    "# \"Reference line corresponds to the expected value if the contribution where uniform\"\n",
    "pc_contributions_distribution(pca, axes=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cut-off point for PC1 is much more to the right (PC1 has a heavy tail), thus more variables contribute to this dimension (compared to PC2, below). Dr Tim noted that similar pattern appears when the data are not centered, as then the first PC is essentially swapping the dataset into the origin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the z_log_matrix might not be centered on zero (due to the order of z-scores) the `prcomp` function does center it anyway (when scale=T)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAEsCAIAAAD2HxkiAAAgAElEQVR4nO3deVwTZ/4H8O8k4UoCKIEg4RC5PLC1WK14ra1ID7XW9ar7s151Vy26tru1/alVWmz7w3qtolW01tZqay9FV7DF1qNaKorWW4GoIIdyRUBCSCCZ5/fHYMAQDgPME8j3/Yevh2Geme/EfJh7hiGEAEKIHgHtAhCydRhChCjDECJEGYYQIcoohDA9PX3KlCl+fn4ymWzMmDFXrlxp81kkJSUNHDiwzSdrmd27d/v5+e3Zs6e9Z2R2qfPz8+3s7Np71kZPPPEE81BgYGBUVFRlZSVvc++g+A7hmTNnwsPDAwMDDx06dPz4cR8fn5EjR969e/dxp+Pr63v+/PmGE2cYJiMjw7JJWda9WevXr4+JiXnttddaP6l2qrBpZj/qJqxbt66goCA/P3/btm1paWkjR47kjsBbXHxaWtqwYcOcnZ0DAgJWr179uN2tH98hfPvtt+fPnx8bG9uvX78nn3xy27ZtAQEBn332WeunXFxcrFAovvzySy8vL8umEBgY2JrujdFqtQEBAa2fTusXkB+urq6enp4KhWLUqFFHjx7NyMhISEgASz/eysrKV155ZdSoUVlZWbt27dqwYcPnn3/ePoXTQ3iUmZkJAPfu3as/MCcnR6lUcu3Tp08PHjxYKpWGhoZ+9913hJCSkpIuXbr88ssvEyZM8Pb2njZtWnV19YgRIwQCgUwmO3DgQElJibu7+4YNG4RC4cGDBwFAr9cnJiY+/fTTH3zwgZeXl6+v7/vvv28wGNLT06VSqXG+o0aN2rNnj8mkuO4tr8RkARv2GjdunFAodHNz+/rrr+uPOWbMmPfee49rV1ZWisXiM2fOEEIOHToUGhrq4ODg7e0dGxvLzbfhApod0+xS5+XliUQi4+cfGRnJrVI+++wzbmBsbKyPj49YLH7uuedu375tskT1Px+zC2iib9++O3bsqD9k1qxZc+fO5RaEK577JPfs2TNmzBhfX98dO3bExcVFRkZ6eHhs2bLFZIKnTp2SyWQsy3I/Llu2bPLkyQ3n26HxGsKkpKSuXbs29tvi4mJXV9eNGzcWFRUdOHBAIpGkpaWVlJTY2dmtWLGCEFJZWenp6cn93/v4+Jw7d44QUlJSIhKJpk6dmpeXZ/xvTkxMFAgECxcuLCgoOHbsmKenZ3x8vNkQmkyK6/5YlTRdPyGkV69eJ06cMFnYPXv2hISEcO1vv/22T58+hBCdTieRSDZu3FhQUJCUlCQUCpVKpdkFNDum2aU2hlCr1Xbv3v2jjz4qKipKTk7u2rVrYmLisWPHpFJpSkrKnTt3Ro8ePW3atIb/L8bPp7EFrK9hCGNiYiIjI8mjIRQIBKtWrSKEJCQkMAzz0UcfcW0XFxeTCarV6pycHK6t1+uHDx++Zs2axr5CHRSvIdy/f79cLm/stzt27AgPDzf+OG/evDfffJP7n8vLy+MGTpgwYfPmzaRBcrKyski9/+bExEQXFxetVsv1Wr9+/YgRI1oewseqpOn6SSMhrKiocHJyunLlCiFk/Pjx3BdLp9NdvHiRG+H69etSqTQ1NdXsApod0+xSG0OYkJDQr18/YwHR0dELFy5MTEx0dnY+cuQIN82KioqG/y/Gz6exBayvhSEEgLt37xJCioqKTNrcer6hjIyMyMjIkSNH3r9/3+wIHRev+4RBQUFFRUVlZWX1B/73v/9dsGAB9/0OCgqqP3JOTg7X9vT05BoCgfmCvb29TYb4+vo6ODhw7V69ejU89sOybGN1WlZJE70akkqlY8eO/fHHH8vLy3/55RfusI29vf2ZM2ciIiL69++/YsWK+rMwWcDGxmxiqW/fvp2enu710Pbt28vKykaPHr127dolS5Z07dp15syZ3P6CBR9LE7Kzs3v06NFwuJubGzz8GOu3AeDzzz93dnZ2dnZ++eWXAaCqquqdd94ZOnRoRETEkSNHunbt2uxMOxYRnzPr06dPz549N23atGLFCuPArVu3KhQKAPD29k5OTjYOv337to+PD9dmGKbpKTeMRG5urk6n476RN27c8Pf3BwCDwUAI4aZ2586dxqZmWSVN9DLrb3/7W3R0dI8ePUaOHNmtWzcASE5Ofu+9906cOBEaGgoA3ECzC9jYmGaXmuPl5RUWFnb69Gnux+LiYpZlc3Nzn3vuublz5xYUFGzatCkiIkKlUjXxl+6xFhAA1Gp1QkLC4x54mzNnzpw5c7g2IWT8+PH29vbXr1/38PB4rOl0FLyuCYVC4YYNG2JjY6Ojo69evZqZmbl06dLffvvt3XffBYBx48Zdu3Zty5YtZWVlhw8f3r1797Rp05qYWnl5eRO/rays/Pe//11cXHzy5MnVq1fPmTPHzc2tqqrqq6++qq6u3rRpU35+fmOTetxKLOv10ksv5eTkrFq1avbs2dyQoqIisVgsl8uLiopiY2MLCwvv3btntm9jYzZcamOXF1988ebNm/Hx8eXl5ampqU8++eTJkyd/+umn8PDwS5cu2dnZSaVSsVhsNoHc59PCBayoqFCpVEVFRSdPnoyMjAwJCZk4cWKzn15jkpOTz549+8knnzx48ODWrVu3bt0qKCiweGpWiv8t4LNnz44ePdrT01Mmk0VGRp49e9b4q5SUlEGDBkkkkj59+nz77bek3o4EN8KkSZO4PbHFixc7Ozvv37+//gjG9smTJxctWrR48WK5XO7v779u3Tqu+5o1a9zd3b29vZcuXbp06dKkpKTGJtXySupr2Is0sk/ImTlzpoeHh/Eoq1arffXVV52dnQMDAz/55JO1a9eGhITk5uY2XECzYyYnJzdc6vpHR8+dOzds2DCJRNKjR4+4uDhCiE6nmzNnjkwmE4vF4eHhf/zxR8MijZ9PYwtYX9++fY1fLX9//3nz5hn3M032Cbl914Ztk33CmJgYk2/spEmTzH6YHRdDOtStTNXV1Xq9XiwW0y7EiirpTEpLSzvfLl+zOlgIEep88AJuhCjDECJEGYYQIcowhAhRhiFEiDIMIUKU8XrZWk1NDXfFZmMXf7XkfInZvvx3tNqZYrWWdWxJ3zapVigUCoXC+r/iNYRNX2iGkI1wd3ev/yNujiJEGYYQIcowhAhRhiFEiDIMIUKUYQgRogxDiBBlGEKEKMMQIkQZhhAhyngNobJQzefsEOoQeA3h/rQ8PmeHUIeAm6MIUYYhRIgyXkP4oKoGn6+IkAleQ3inRFOp1fM5R4SsH26OIkQZhhAhyix/vIVer4+Li1Or1X5+frNmzeIGEkI+++yzoqIiFxeXf/7znw2fyaHW6aWOvD5TAyErZ/maMDU1VaFQREdH37t3Ly+v9gTguXPnJBLJ8uXLw8LCCgsLTboQgEt3yhpMCSGbZnkIlUplQEAAAAQEBCiVSm7g9evXDQZDXFxcVVVV/XdcGt0urrR4jgh1SpaHUKPRcA+NkslkanXt9WhqtbqgoGDq1Klnzpy5cOECNzAuLi4yMjIyMlL7oKT1FSPUyVi+eyYWi1UqVUBAgEqlksvl3ECJRDJ48GC5XD58+PBbt26FhYUBwLBhw7y8vAAg/lIzb71GyAZZviYMDg7Ozs4GgJycnKCgIG5gUFDQzZs3ASArK8vT05Mb2L9//8mTJ0+ePFnk4NTaehHqdCwPYXh4eH5+/po1a+Ryua+vb2Zm5ubNm8PDw7OyspYtW1ZSUjJkyJA2LBShzorXN/W+8GHSsBD3NyICeZsjQlYIn8CNkHXBECJEGd8hfFCFF3Aj9Ai+Q3hHhSfrEXoE3yFU461MCD2K7xDiTb0ImcADMwhRhiFEiDLeN0dZ0BtYnmeKkDXjf5+QVBtwxxChOrg5ihBlGEKEKKMQQp0e9wkRqkMhhNn4hAuE6qEQwpv4biaE6qEQwvuVNfzPFCGrhQdmEKIMQ4gQZRhChCjDECJEGYYQIcowhAhRRuU8YQX/M0XIalEIYfEDLf8zRchqUQghj08bRqgDwH1ChCijEEIW14QI1UMlhERbY+B/vghZJ9wcRYgyOiHU4yYpQg/RCeGp9GIq80XICtEJYX5pFZX5ImSF6ITwbhmer0eoFp0Qlmnw5nqEauHRUYQooxPConKdWofvSEMIgFYICRB88CFCHNwcRYgyDCFClGEIEaIMQ4gQZRhChCijFsIUZQmtWSNkVaiF8FreA1qzRsiqUAtheRVeuYYQAMUQGgzEgHcVIkQxhAQAI4gQ4NFRhKjDECJEGYYQIcpohhDvZkII6IYw9uANinNHyErQDGF+mRZXhgjRDCEh5GxWGcUCELIGlA/M3L2voVsAQtTh0VGEKKMeQrxsBtk6yiH89WoR3QIQoo5yCGsMLME39yLbRn1zFKqq8V2FyKbRD+GRKwW0S0CIJvohPJtVSrsEhGiiH0KEbBz9EOpZlnYJCNFEP4RF5TraJSBEk6jhIJZlDx06dObMmYsXL7IsGxYW9swzz7zyyisCQbsktsZAruWVh/q4tsfEEbJ+prn66aefxo0bd+nSpUGDBq1fv37jxo2DBw++evXquHHjfvrpp3Yq4kY+Pv4Q2S7TNWFJSUlCQoKdnZ1xSM+ePceNG6fX6/fu3dtORaRllU4a5NtOE0fIypmuCadPn14/gXq9Pjk5ubq6WiQSTZ8+vZ2KyC/VavCUPbJVzezmTZo06c8//3z11VfbtQgC5OD5u+06C4SslmkIb9x45JETubm5S5cuvXvXTEL0ev369etXrlz55Zdfmvzq0qVL//nPfx6rjhv55Y81PkKdhmkIt27dumDBAmPqBg4c+Pe//71Pnz4Ne6ampioUiujo6Hv37uXl5RmHazSavXv34mXZCLWQ6YGZuLg4pVK5ZMmS7t27v/POO/Hx8VevXg0NDW3YU6lUcsMDAgKUSqWPjw83fNeuXZMmTTp58qRxzJKSkgcPHgAAMegZoZmTIgBQocWHzSAbZWafMDg4eNeuXWPHjp09e3ZcXFxISAjDMA1H02g07u7uACCTydRqNTcwLS3N09PTz8+v/pirVq0aMGDAgAEDqsqLG6ujRK3DlSeyTWZCGB8fP2TIkAsXLuzbt69Hjx6TJk369ttvG44mFotVKhUAqFQqqVTKDTxw4MC1a9c+/fTTq1evJicncwPXrl17//79+/fvi928GquDJZBfWtUGC4RQR2MmhNu2bTt9+vTOnTsB4OWXX05ISKisrGw4WnBwcHZ2NgDk5OQEBQVxAz/++OMVK1YsWLCgb9++L7zwwmOVciUXT9kjW2QmhE5OTikpKQ4ODtyPQqFwzpw5DUcLDw/Pz89fs2aNXC739fXNzMzcvHlza0opVle3pjtCHRTTcE9MqVTGx8fPnz8/ODi4bWf2wodJTfw2yFMaM8HMESCEOhnuYIqR6bHKqqqq4ODgdevWNexZVVXl5OTUfpXdUZnZ6EWo0zMN4VtvvaVQKObPn+/p6WkcWFhYGB8fn5+fv3379vYrpUaPR0eRLTIN4bZt21JSUmbPnp2dnc2tNEtKSvz9/d97772hQ4fSqBChTs7MqfOhQ4cePnwYALjrYIxn4RFC7cH89Ssc/uNHCDF7YQBCnRj9x1vUd69MS7sEhPhmXSEsxVOFyPY0H0I+L+lMuoQPAkY2x3wIjx49+tZbbwHAX//6VxcXl4Z3DLaTsipcEyKbYz6E8+bNGzt27OnTp7VabUZGxscff8xPNcUPMITI5pgPoU6nGzVqVFJS0vjx4xUKRU1NDT/V4N1MyAaZP0XRr1+/mTNnnjx58sKFC++//76XV6O3ILUtbY0hV6XxlYn5mR1C1sD8mvDrr78ODw8/cOBAly5dhEJh+z3s0AQBOHThHj/zQshKmF8TCoVChUKRnp6enp7es2fP8+fP+/v781OQCs9SIBtjPoSvv/66TqcLCAjgfhQIBBMnTuSnoMJyvL8e2RbzIdTr9QcPHuS5FE6ppuZybtmTvl2ozB0h/pnfJwwKCiorK+O5FKN1hzNozRoh/plfE+bn5ysUisGDB7u5uXFDfvjhB95q0rNws7AiyNOZtzkiRJH5EEZFRUVFRfFcSn23izQYQmQjzG+ODh06VK1WJyYmHjx4sLy8HG/nRaj9mA9hTExMdHS0QqHw8/NbuXJlTEwMz2WlKEt4niNCtJh52hoABAUFXb58WSwWA4BGowkLC8vIaIODJU0/ba0+O6Hgy7kDWz9HhKyQydPWGr2VyXiHezu9JbtpegNeRIpshfkDM9OmTYuIiJgxYwbDMLt3727v9xMiZMvMhzAmJmbQoEEnTpwAgCVLlowdO5bXohCyJaYhHDhw4I8//rh48WLjkKysrF27dvF5nhAhm2IawsWLF7u5uXG31SOEeGAaQm7375tvvvn000+NA2fNmsXzqUICpFKnlzg09URGhDoH02957969ASAnJ+fYsWPcEL1e36ULhcupD/55938G+zU/HkIdnGkIf//9dwCIiorasmWLcSCVEP6WXowhRLbANIQymQwAZsyYkZqaahwokUieffZZPssCgEqtgec5IkSF+Z2uH3/8kWsUFhaeOHFi0aJF/IeQACmq0MqdHXmeL0I8Mx/CL774wtguLCx84403+KrnEedvl77Uj6dnTCFES/OXpHl4eLTJhaMWuHYX32KPOj/za8LZs2cb21evXn3qqaf4qudReAEpsgHmQzhp0iRj+7XXXqN1P+HFnLKCcm03V9wtRJ2Z+c3RiIiIjIyMHTt2fPXVV1lZWSIRnZPmhMC+tHwqs0aIN+bTNW/ePJVK9cYbbzAMs3Xr1itXrmzcuJHnyjjX88upzBch3pgP4fHjx2/dumVnZwcAzz//fK9evWiFUKfHs4WokzO/Odq9e/fy8tpVUGlpqcmNwHyqqmYfVOlpzR0hHpiuCbnHyXTv3r1///4TJkwAgH379s2YMYNCaQ8t+fbK8r/2UnRxolgDQu3HNITdunXj/h0+fDg3ZPny5VSuHTUq11bv+v3O0rG9KNaAUPsxDeGOHTtMburl0H3CRXo+nrVHnVbHuKlXz5J7ZVqvLnjCEHVCZm7qvX///qlTp5YsWUKloMb8cq1gxlB/2lUg1PbMHB11dXXds2dPVlYW/9U04di1YtolINQuzJwnFAqFoaGhYWFhQ4YMkUgk3EDqD3qqMbA6PesgovAQVITalfmT9YsWLVq0aBHPpTTryJWCl8MUtKtAqI2ZDyH1Bz2ZlVeKL/FFnZD1PuipoT+zqb23FKH2Y70PempIU63PKHjQs5sL7UIQakvm38qkVquPHj2q1Wq5H0Ui0cSJE1s/s5a/lakxQgHz+Zyn7UTC1heDEC0mF2Ob3yd8/fXXdTpdQEAA96NAIGiTELaegSVX76nDfF1pF4JQmzEfQr1ef/DgQZ5LaaG4nzN3/mMgQ7sMhNqK+dNuQUFBZWVWehSkWs8euVJIuwqE2oz5NWF+fr5CoRg8eLCbmxs3hPrJ+vq+P5P7fF+58TWmCHVo5kMYFRUVFRXFcyktp60xHLteFBHqSbsQhNqA+c3RoUOHqtXqxMTEgwcPlpeXW8OZehN/5ljp1jJCj8t8CGNiYqKjoxUKhZ+f38qVK7nb7a1KyYNq2iUg1DbMnycMCgq6fPmyWCwGAI1GExYW1iYP4W79eUIj765Oq6c+2VZTQ4hPJucJG70pwXjYQyCwxhsXSh7oaJeAUNswf2Bm2rRpERERM2bMYBhm9+7ddJ9tYZbOwFZV653s8VW+qMMzvzkKAIcPHz5x4gQA/OUvfxk7dmybzKwNN0cBYMozPq887d2GE0SIH81vjrIse/fu3dGjR69evXr48OH8v5mwhY5cxVP2qDMwDWFmZmbv3r23bdvG/bhhw4Y+ffpcvnyZ98KaV6apuVWspl0FQq1lujn6/PPPT5gwYf78+cYh33///fbt23/99dfWz6xtN0cBQCBgds97pm2niVB7a2Zz9Pz586+99lr9IVOmTKH1ktBmsSz59To+AAp1bKYh9PDwMDknUVVV5ehovQ/8/D41h3YJCLWKaQhfeumlZcuWsSzL/UgIWbFixahRo3gvrKX0D0tFqIMyPc8WGxs7ZcqUXr16DRs2TCgUnj59WiaTJSYmUimuJXQ1bKmmuqvYnnYhCFnI/HnC8+fPnz9/vrq6+oknnhgxYoTZnnq9Pi4uTq1W+/n5zZo1ixuo0WhiY2P1er1EInn33Xft7R/JRpsfmOHYCwXxs/s72OEzL1DHYHJgptGT9c36/fff8/Lypk6dGhsbO336dB8fHwBISkrSarUTJ07cu3evh4eHyXZsO4UQAPxkktgpfdtp4gi1rRY9Y6YllEplaGgoAAQEBCiVSi6EwcHBMpkMAKRSKfeiXwA4cOBAamoqAFTDk/aSdnk8TOEDfCQp6qgsvzhbo9FwgZbJZGp17UnzkJAQmUyWlpb2xx9/DBw4sG1qbAFdDfvrNbyABnVIlq8JxWKxSqUKCAhQqVRyudw4/IcffsjLy1u+fDl3JxQAjB8/fvz48dCem6MA8FXKnVF4rz3qgCxfEwYHB2dnZwNATk5OUFAQN/D06dNqtfpf//qX8U0yvDEYyKEL93ieKUKtZ/mBGb1ev3nz5pqaGrlcPnPmzMzMzCNHjohEohs3bjg5OQHAmDFjjO/c5rTrmhAApA6iba8/3a6zQKj12uzoqAXaO4QAEDcjTCbBc4bIqrX0zvoOaun3V1ge/6wg1HqdLYSVWv2Xp7JpV4HQY+hsIQSAE9eLK3V62lUg1FKdMIQGQuZ98WfixQLahSDUIp0whABACNl7+k6uSkO7EISa1zlDyMGHX6AOoTOHcE/KHT5PwCBkmc4cwqpqdu4X50s1+MB8ZNU6cwgBQKMzvLn7Yo0B775H1quThxAADCz5LjWXdhUINarzhxAAfr5SoGdx5xBZKZsIISEQ9cWfJWp8hwyyRjYRQgCorNYv/uZyeVUN7UIQMmUrIQSAGgN74M+7tKtAyJQNhRAAfrteRLsEhEzZVgh1evZaXjntKhB6hG2FEABWJ2XicVJkVWwuhHqW/SrlDu0qEKpjcyEEgCOXC5Z8d7WwXEu7EIQAbDOEAJB7v/LtvZeOp+NxGkSfjYYQAAiBHcezViWl0y4E2TrbDSHnSk75zUK87RDRZOshBIBPEjPwtkNEEYYQNNX6mdvTDHiFN6IEQwgAYGDJvnN5tKtANgpDWOvg+bvbjt2mXQWyRRjCOiczilcnZuBWKeIZhvARl3LLdv+eTbsKZFswhKZ+uVpUXIG3/yL+YAhNsYS8/c3l/DJ8/zbiCYbQDAPLvrv38tt7L+nxMW2o/WEIG1VQpp25PW3/ubyqagPtWlBnhiFsxr60/Hk7z/2hLKZdCOq0MITNMxD49Nfb+9PyaReCOicMYUvtO5c3Pf7M5qO3aBeCOhsM4WNgCZzOLJm1PU2tw71E1GYwhI+txsDO/+Lcn1mltAtBnYSIdgEdEiGw7udMRsA829Pj78/2oF0O6thwTWg5wpLjN4rWHM6kXQjq2DCErXXxTumawxm0q0AdGIawDVy8U7Ydb4NClsIQto3fMkq+ScmhXQXqkBg+H6+SMOQV2YP7XDu5f+S54Ke5dki+cvLv+7h2jchu9cS3jV1mHd3tXVJ7lvxk3+GnQodybb/i3OnHvjaO9smkxXph7UGmV0/9EHS39mxeWsiAI2GjuLZHecncn3cYu2wct1DtJOXa484kPpF9lWtf8e/730Fjuba0Sv3mfzcbu2x/8e/Fru5c+/kLvw7MPMe1byoCvxs+GQDshExoN8kHB9YJGIb7leZ//9cQEsK1HRIT7RMSuLYhMFCzbJlxys4LFoC29jmoVQsX6sPCuLb98eMOe/ZwbdbTs/L//s/YRfruu4xKxbW1s2bVDB/Ote3OnnXcto1rE6lUvXGjsYtk5UrBndoHH+smT65+8UWuLbp+3WndutqRGKZiR92nJF67VnjjBteuHj1aN3Ei1xZmZ4s//NA4mnrTJiIWc22nrVtF52o/mZoRI7QzZnBtQXGxZMkSY5fKVatYDw+u7bh7t92JE1xbP2BA1Rtv1NZSVSVduNDYRbNihcHfv/bD3L/fPimp9sPs1Uvzzjt1H+Y//gFs7UW/Vf/6l75v39oPMznZ4fvvaz9MX9/KDz4wdpG+9RZTUcG1tXPn1gwaxLXtUlIcd+7k2sTNTb1mjbGLZPlywb17XFs3bVr1yJFcW3TxotOmTbUjOThUbNli7CKOjRXevFk9dqzzP/4B9fB6dPS6b2+ptpJrF7t6GIeXS1wuBD7FtQ0CYf0uGd7BRQ/HvOfWzTi8wklq7AIALFO3Sr/VLaDCyZlr57r7GIdX2TvW71IjsjO278j9jBnOlynqj1O/i9bewdjOdfcRGfRc21hhjYFcyqvYznTvIrGbOMAHAIirq7GLwc+vZsSI2oI96hYfAGqGDYOa2te2EXf3ui5eXsYuxMXlkS6DBjGVtR8m6+VV91HI5XVdHBwe6dK/v+Dhl5j19a3r0qWLsQs8/PPB0ffrx8rltcX0qDsOTJyd67oAgKjui6Tv04dIJLVdHv4BAgDi6Fi/C3F0rFvM4GB4uDIw+PnVjSMUPtLF2bmuS48edR9mt7ovBgDU/OUvxqmxXbvWLaavb90n4+b2SJfBgxld7f1rxuUFANbTs67Lw4WqXcyBA5ny2veaGBR13xlWJqurWfRIvvRPPcV6exv/jhjxuiZ84cMk3uZFHcPAosjgZwLdmh8V2Rj3en9nAfcJ2w8hsPEXZfT+a/g8RdQ0DGF7InCrUD1j29kNP9+s1OlpV4OsFF4x0+5YAmlZqrRslVDADPDvOn9koL0I//ahOhhCvhAwGMiZW/fP3LpvJxBE9PWYPtSfdk3IKuCfZApqWPbny4Vzd54z4O4iwjUhRZU6w4z4syKGCVQ4LxoV2EVsT7siRAeeorAKDECIl3TyIN/eXi7Nj406OJNTFLgmtAoEIOOe+qMDNxgAgUDgKxN/8NfedkLcWbAJGGC7qTEAAASJSURBVELrQgAMLJtdrJ61PU3IMK5iu97erv8zxK+LE/5PdVr4X2u9DITcr6xOySxOySwWCQTvvdI7pJuUdlGo7WEIOwY9y8YkXBMKmGBPyT+fD+kitmu+D+og8MBMR+UgEvT1dXnz+RChgGl+bGRNTA7MYAg7PDuRYHx/xfinvWkXgloKQ9h5MWAvEHi42kX09hoZ6oEHV60WhtBWOIiYHu6S14b16OEhpl0LegSG0PYwIGIYmYtDRB+PiD6ejnbC5rug9oQhRCASCv3lTiN7yUf08mh+bNTWMITIBCMQMvYMODvZ9/J2jnxCHiCTMAwecW1HeNkaMkFYA9ECaCu0xenaU+nFAMAACIVCRReHPt7Og4NkPeRSIcay3WAIkRkEQG8w5Kg0OSrNz5cLuYECAWMvEjg72D3h22X8AC+ZBG/7aBsYQtRSLEu01QZtteHY9YJj1wsAQMCASCCQOgp7K6QRfbv17Ia3gFgCQ4gsxxKoNrD3K9kUZWmKshSg9kisvUgodhB2kdg928tjRG85nq9sGh6YQTwRMCAQMM72Ij93p9FPeff1sd3VJh6YQXSwBFgDKa2qKc2tuZT7AAAYAAIgEDAOIqGjg7Cro6ibq+OL/boFym3rZhHLQ6jX6+Pi4tRqtZ+f36xZs5oYiJBZ3DYYy5Kqan1Vtb60Qne7uPKPm7XP9gcGGAZEjFAkYhzsBGJ7gcRRJHNy7Ops383ZYWBQV1fHTnIrieUhTE1NVSgUU6dOjY2NzcvL8/HxaWwgQpYgQAjUgKHGAFU6KAMAACWouV9+cSqLazAADAOEYQSECAWMQCgQ2wvF9iKJvcBVau/sJFR0cfJzl3ZzdXBzstLDuZaHUKlUhoaGAkBAQIBSqeTyZnYgQu2HABACQIgBwGAgYDBoqw33obrZjgwDQBiBABiGETJEKBCIhIyIYRwchU4iob29UGIntBMKRAwjdhS5OAmdHURiRzuJvcDBQSh2EDmKRPYCAsA4O4qEglYde7I8hBqNhtu/lMlkarW6iYGLFy/euXMnAITOXC1282pkegjxihAAIAYWAIgeAMBQ+wtNTZtMn3m4vW3yo0jAJL03uv6YlodQLBarVKqAgACVSiV/+CIbswOXLFkSFRUFALdUNVrCbeqzhAEAEBAghCEMMACEMFyhDNdiGJYhQGoP4JKHV1IxQAgDDOG61A6v6wuEASACII/bscFMWQBgGAHDACGEISbVAmHAmqo19mUAjNUCAAsMADBtUq2xL3C926ZaqNe3tloAIAICpK2q5cYlAFxXrmMbfP0ACPcWq4fVAgNMcx29XE1varE8hMHBwdnZ2QMHDszJyRkyZEgTA93d3bnVo4tLicWzQ6izsnxbNjw8PD8/f82aNXK53NfXNzMzc/PmzSYD27BQhDorXk/Wl5TgmhAhfD8hQlYGQ4gQZRhChCjDECJEGYYQIcowhAhRhiFEiDJrCaHFj/dqzXPBqMzUYlitFc7Ugo4Nu/B6sh4h1JC1rAkRslkYQoQowxAiRBmGECHKMIQIUYYhRIgyDCFClP0/V7VoNchLofIAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 300 -h 300 -u px\n",
    "pc_contributions_distribution(pca, axes=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While PC2 contributions are much more clearly defined by a small number of proteins, these dimensions do not account for enough variation to be shown on the original biplot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting the top contributors for each of the PCs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loadings are named \"rotations\" in `prcomp`; columns correspond to the eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -o loadings -o stdev\n",
    "loadings = as.data.frame(pca$rotation)\n",
    "stdev = as.data.frame(pca$sdev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "      <th>PC4</th>\n",
       "      <th>PC5</th>\n",
       "      <th>...</th>\n",
       "      <th>PC78</th>\n",
       "      <th>PC79</th>\n",
       "      <th>PC80</th>\n",
       "      <th>PC81</th>\n",
       "      <th>PC82</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>CHIP</th>\n",
       "      <td>-0.030286</td>\n",
       "      <td>-0.010796</td>\n",
       "      <td>0.019527</td>\n",
       "      <td>-0.018871</td>\n",
       "      <td>-0.001675</td>\n",
       "      <td>...</td>\n",
       "      <td>0.019079</td>\n",
       "      <td>0.002251</td>\n",
       "      <td>0.022151</td>\n",
       "      <td>-0.027124</td>\n",
       "      <td>-0.630446</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CEBPB</th>\n",
       "      <td>-0.013677</td>\n",
       "      <td>-0.053114</td>\n",
       "      <td>-0.035227</td>\n",
       "      <td>0.040743</td>\n",
       "      <td>0.016014</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.048504</td>\n",
       "      <td>0.010487</td>\n",
       "      <td>-0.029901</td>\n",
       "      <td>-0.013329</td>\n",
       "      <td>0.216807</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSE</th>\n",
       "      <td>0.033534</td>\n",
       "      <td>-0.011904</td>\n",
       "      <td>0.032665</td>\n",
       "      <td>0.024509</td>\n",
       "      <td>-0.001668</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.038495</td>\n",
       "      <td>-0.003225</td>\n",
       "      <td>0.016193</td>\n",
       "      <td>-0.005210</td>\n",
       "      <td>0.366605</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PIAS4</th>\n",
       "      <td>0.001750</td>\n",
       "      <td>-0.052516</td>\n",
       "      <td>-0.044091</td>\n",
       "      <td>0.049317</td>\n",
       "      <td>0.019133</td>\n",
       "      <td>...</td>\n",
       "      <td>0.009833</td>\n",
       "      <td>-0.070958</td>\n",
       "      <td>0.014760</td>\n",
       "      <td>-0.022986</td>\n",
       "      <td>0.276844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IL-10 Ra</th>\n",
       "      <td>0.032442</td>\n",
       "      <td>0.015214</td>\n",
       "      <td>0.002669</td>\n",
       "      <td>-0.007186</td>\n",
       "      <td>0.022590</td>\n",
       "      <td>...</td>\n",
       "      <td>0.008955</td>\n",
       "      <td>0.064745</td>\n",
       "      <td>0.062419</td>\n",
       "      <td>-0.017625</td>\n",
       "      <td>0.089330</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 82 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               PC1       PC2       PC3       PC4       PC5  ...      PC78  \\\n",
       "CHIP     -0.030286 -0.010796  0.019527 -0.018871 -0.001675  ...  0.019079   \n",
       "CEBPB    -0.013677 -0.053114 -0.035227  0.040743  0.016014  ... -0.048504   \n",
       "NSE       0.033534 -0.011904  0.032665  0.024509 -0.001668  ... -0.038495   \n",
       "PIAS4     0.001750 -0.052516 -0.044091  0.049317  0.019133  ...  0.009833   \n",
       "IL-10 Ra  0.032442  0.015214  0.002669 -0.007186  0.022590  ...  0.008955   \n",
       "\n",
       "              PC79      PC80      PC81      PC82  \n",
       "CHIP      0.002251  0.022151 -0.027124 -0.630446  \n",
       "CEBPB     0.010487 -0.029901 -0.013329  0.216807  \n",
       "NSE      -0.003225  0.016193 -0.005210  0.366605  \n",
       "PIAS4    -0.070958  0.014760 -0.022986  0.276844  \n",
       "IL-10 Ra  0.064745  0.062419 -0.017625  0.089330  \n",
       "\n",
       "[5 rows x 82 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadings.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helpers.unsupervised import top_variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"highlight\"><pre><span></span><span class=\"k\">def</span> <span class=\"nf\">top_variables</span><span class=\"p\">(</span><span class=\"n\">loadings</span><span class=\"p\">,</span> <span class=\"n\">pc</span><span class=\"p\">,</span> <span class=\"n\">n</span><span class=\"o\">=</span><span class=\"mi\">20</span><span class=\"p\">):</span>\n",
       "    <span class=\"k\">return</span> <span class=\"n\">loadings</span><span class=\"p\">[</span><span class=\"n\">pc</span><span class=\"p\">]</span><span class=\"o\">.</span><span class=\"n\">abs</span><span class=\"p\">()</span><span class=\"o\">.</span><span class=\"n\">nlargest</span><span class=\"p\">(</span><span class=\"n\">n</span><span class=\"p\">)</span>\n",
       "</pre></div>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_source(top_variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def top_variables(loadings, pc, n=20):\n",
    "    return loadings[pc].abs().nlargest(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick check to verify that my results are same as of the R-package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hnRNP A/B    0.039330\n",
       "IMB1         0.039219\n",
       "STAT3        0.039196\n",
       "PTK6         0.039172\n",
       "FGF9         0.039125\n",
       "Name: PC1, dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top_variables(loadings, 'PC1', n=5) # I do not scale by variance - do not need it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAADICAIAAABOPGvMAAAgAElEQVR4nO2deVxTx/r/JwlbQhBkLWEVZBGpCBQE1IvFpRWr9lqq9qe1Krbutm794l6sXmxdQatogYtr8V4qUoGIilAsIpsKKvsaliAEJSRAAiHn98fBNDecBMRzEjHz/oPXZM6c88wMn0zmzDnPPCQEQQAEoh6QVV0BCER5QLlD1Agod4gaAeUOUSNUL/fS0tKFCxdaW1sbGRnNmTPnyZMnuJtITk728vLC/bLD4+LFi9bW1pcuXSLaEGarGxsbNTU1iTYt4f333ye9wt7eft26dZ2dnUqzPhAVyz0nJ8fHx8fe3v7GjRvp6emWlpYBAQFNTU2vex0rK6uCgoKBFyeRSGVlZcO71PBOH5Rjx46FhoYuXbr0zS9FUA0Vg9nVCjh69Ghzc3NjY+PZs2fz8vICAgLQxcA3rHxMTMzGjRtf9ywVy33r1q1r1qwJCwtzc3ObMGHC2bNn7ezsfv311ze/cmtrK4PBiI2NNTc3H94V7O3t3+R0eQgEAjs7uze/zps3UDno6+ubmZkxGIwZM2akpaWVlZUlJCSAN+jex48f//DDD//3f/83nNogqqO8vBwAwGazpTNZLFZFRQWazs7O9vX1pdPp48ePv3r1KoIgHA7HwMDg9u3bCxYssLCwWLJkSU9Pj7+/P5lMNjIyun79OofDMTY2PnHiBIVCSUxMBACIRKKkpCRPT88ffvjB3Nzcyspq3759fX19paWldDpdYnfGjBmXLl2SuRR6+tBrItPAgWfNmzePQqEYGhpevnxZuuScOXN27dqFpjs7O2k0Wk5ODoIgN27cGD9+vLa2toWFRVhYGGp3YAMxS2K2uqGhQUNDQ9L/M2fO1NPTQ4cYNDMsLMzS0pJGo3344YfV1dUyLZLuH8wGyuDq6hoVFSWds3z58m+++QZtCFp5tCcvXbo0Z84cKyurqKioiIiImTNnmpiYnD59euA1o6Ki1qxZ4+LismHDhoFHFaNKuScnJ48ePVre0dbWVn19/fDw8JaWluvXr+vq6ubl5XE4HE1NzT179iAI0tnZaWZmhvaypaVlfn4+giAcDkdDQ2Px4sUNDQ2SDk1KSiKTyRs2bGhubr57966ZmVlkZCSm3GUuhZ7+WjVRXH8EQZydnTMyMmQae+nSJUdHRzQdFxfn4uKCIIhQKNTV1Q0PD29ubk5OTqZQKBUVFZgNxCyJ2WqJ3AUCgY2NzYEDB1paWlJTU0ePHp2UlHT37l06nZ6VlVVXVxcYGLhkyZKB/xdJ/8hroDQD5R4aGjpz5kzkf+VOJpMPHTqEIEhCQgKJRDpw4ACaHjVqlDx5rFmzZoTJ/dq1a6ampvKORkVF+fj4SD6uXr3622+/RfuooaEBzVywYMGpU6eQARqtqalBpDo0KSlp1KhRAoEAPevYsWP+/v5Dl/tr1URx/RE5cufxeFQq9cmTJwiCfPrpp4cPH0YQRCgUPn78GC1QXFxMp9MfPHiA2UDMkpitlsg9ISHBzc1NUoG9e/du2LAhKSlJT0/v1q1b6DV5PN7A/4ukf+Q1UJohyh0A0NTUhCBIS0uLTBr97RrI8OSuyrn72LFjW1pa2tvbpTP/+OOP9evXAwAaGhrGjh0rXZjFYqFpMzMzNEEmY9ffwsJCJsfKykpbWxtNOzs7D7wbFovF8uo5vJooOGsgdDr9k08+iY+P53K5t2/fRm9ktbS0cnJypk+f7uHhsWfPHmkTMg2UV1JBq6urq0tLS81fce7cufb29sDAwCNHjoSEhIwePfqrr75CZ5vD6BYF1NbWjhkzZmC+oaEheNWN0mkAQHR0tJ6enp6e3ty5cwe9vmJUKXcXFxcnJ6eTJ09KZ545c0YgEAAALCwsqqqqJPnV1dWWlpZomkQiKb7yQPHV19cLhUI0XVJSYmtrCwDo6+tDXr0yVFdXJ+9qw6uJgrMw+eKLL37//ffr168HBAS89957AIDU1NRdu3ZFREQ8fPgwPj6eSqXKa6C8kpitRjE3N3d3d2e/oqio6MiRI/X19R9++GFBQUF5ebmdnd306dMVjAKv20AAAJ/PT0hImDlzpuJiMgQHB/N4PB6Pd+PGjdc6cSCqlDuFQjlx4kRYWNjevXufPn1aXl6+Y8eOP//88/vvvwcAzJs379mzZ6dPn25vb09JSbl48eKSJUsUXI3L5So42tnZuWXLltbW1szMzJ9//jk4ONjQ0LC7u/vChQs9PT0nT55sbGyUd6nXrcnwzpo9ezaLxTp06NCKFSvQnJaWFhqNZmpq2tLSEhYW9vz5czabjXmuvJIDWy055eOPP66srIyMjORyuQ8ePJgwYUJmZiaTyfTx8SksLNTU1KTT6TQaDfP3E+2fITaQx+O1tbW1tLRkZmbOnDnT0dHxs88+G7T3iOJ1Zz+4k5ubGxgYaGZmZmRkNHPmzNzcXMmhrKysSZMm6erquri4xMXFIVITPrRAUFAQOmPetm2bnp7etWvXpAtI0pmZmZs2bdq2bZupqamtre3Ro0fR0w8fPmxsbGxhYbFjx44dO3YkJyfLu9TQayLNwLMQOXN3lK+++srExESywiMQCBYtWqSnp2dvb//TTz8dOXLE0dGxvr5+YAMxS6ampg5stfTKTH5+/pQpU3R1dceMGRMREYEgiFAoDA4ONjIyotFoPj4+9+/fH1hJSf/Ia6A0rq6uEqXZ2tquXr1acj8gM3dH7zEGpvGdu5OQkfwCcE9Pj0gkotFoqq7IW1STd4mXL1+OHj0axwuObLlDIK+F6t+ZgUCUBpQ7RI2AcoeoEVDuEDUCyh2iRkC5Q9SIES93lSykqpVR5dslzuiIl/vLly/7+vqUbJTL5fb29irZKI/Hk7wAozQ6Ozu7u7uVbLS7u5sgH78RL3cIZOhAuUPUCCh3iBoB5Q5RIzRUXQFsUh6yUh4O7hoDABCJRBQKZVCHDwCAE8NgY6CrggIiMdLKHdJtWUeHgNZL1tAY/G6Vpq2hT9MayjUhSmD4cheJRBEREXw+39raevny5dKHDhw4sG3bNh0dHQ6Hs3XrVlNTUwDA5s2bGQzGEC9ezuZWsBW5awwDYe8gCzjXHlRHp5Xia1SfpvWfrYqcd+6XNe//TwG+q25ULY3/bJ2ppSH3p7u2hbf/vwXiISz2icVidFOkQUtqUsg/f+kzmq4tr0B3jyg2vaxHJNc9SgL6jvtQtn+ikEn/b6qDoXyjMgxf7g8ePGAwGIsXLw4LC2toaEAdt/h8fmhoqGSvnNbW1sDAwEWLFg3bijIZyn8C92t2dPXivsLc3SPq7RMrkDuHJ2h8gf9KX0d3jwK5t3C7r+fW4m7Ux9HMkG4yxMLDl3tFRcX48eMBAHZ2dhUVFajc6XT6Tz/9tHfvXrTM8+fPm5qaTp065erqOm3aNDSztLS0oaEBAODl5SXPH0KBi+SwQRBE8WI5Qev3qjLaS5b7PRKJREQYFYlEChqrfKNkMplCoUjnDF/uXV1dxsbGAAAjIyM+ny9tQ/LbR6PRXF1d3d3dw8PDDQ0NJ0yYAAC4fv16fHw8AOD8+fMDtwyQtGHYFZOHWCxW/PCCiCdHCIIoNkrQk6Ouri4gosg7ijq/4053d3dnp9xpD0GPqwQCgbwe1tHRwU3uNBqtra3Nzs6ura0NnZ0PxNvbG00EBASUl5ejcg8JCQkJCVF8cS0t/G/vKBSKgYGBggI6Oq24GyWRSIqN0mg83I0CAPT19XW15f5z6W2EPBLW09MzMNCTd5TbK/fr9ybQ6XTFPSzN8BciHRwcamtrAQAsFkt6vxFpfvvtt8LCQrTMW76VIUQdGL7cfXx8GhsbDx8+bGpqamVlVV5efurUKZkyM2bMiIuL27Nnz8uXL319fd+sqhDImzL8yYyGhsZ3330n+ejo6Ojo6Iimf/zxRzRhYmISFhb2JvWDQHAEPlWFqBFQ7hA1AsodokZAuUPUCCh3iBoB5Q5RI6DcIWoElDtEjcB4zCQWi2/cuJGTk/P48WOxWOzu7u7t7T1//nx5kWEgkJGCrIKZTOa8efMKCwsnTZp07Nix8PBwX1/fp0+fzps3j8lkqqSKEAheyI7uHA4nISFB2pHEyclp3rx5IpHot99+U27dIBCckR3dv/zyS2mti0Si1NTUnp4eDQ2NL7/8Url1g0BwZpDpeFBQ0MOHD0eK9x0EohjZyUxJScm4ceMkH+vr669fvz5p0iQibCvwHoLOe29uFDrvDXTekx3dz5w5s379ekm8WS8vr1WrVrm4uBBRUQhEyciO7hERERUVFSEhITY2Ntu3b4+MjHz69Cnqgo07CnZWIGLRk0QiKd7LQWYkwAtVGdXUlOvMoKFByP5CGhoaChqrEqOyhQdmOTg4nD9/Pjc3d8WKFf7+/mvWrBnKNiMQyNsPxiAaGRnp5+f36NGj33//fcyYMUFBQXFxccqvGQSCOxhyP3v2bHZ2dkxMDABg7ty5CQkJBG22DYEoGQy5U6nUrKwsbe3+7aAoFEpwcLByawWBEAKG3M+fP3/t2jV0dIdA3iVkb1W7u7sdHByOHj06sGh3dzeVSlVKrSAQQpCV+3fffcdgMNasWWNmZibJfP78eWRkZGNj47lz55RbPQgET2Tlfvbs2aysrBUrVtTW1qJbQHI4HFtb2127dk2ePFkVNYRAcANj3X3y5MkpKSkAAHSfXnRrXwjkHUDRgy4odMg7BnRQgqgRUO4QNWJwuaskNjkEQgTYck9LS0N39/3nP/85atSo2NjYgWVEItGxY8f2798/8OiBAwfQ+BAKykAgygdb7qtXr/7kk0+ys7MFAkFZWdnBgwcHlkFDke3du5fNZqNrOAAAPp+/ffv23NxcBWUgEFWBLXehUDhjxozk5ORPP/2UwWBgeotUVFTY2dmBV6HI0Ew0FBkalEZeGQhEVWAvRLq5uX311VeZmZmPHj3at28fZpyZoYQiwyxz6NChQUOR9fT0DL9Ncujr62tvb1dQgIgAXQiCKDba1dWFu1EAAJfL7dWS6zgi/f/CER6P164p1xeRxyOkpXw+v70d271DR0dHR0dHOgdb7pcvX75y5cqWLVsMDAwoFArmlhtDCUWGWebTTz/94IMPAACWlpbyAk0S4flCJpN1dXUVFBi6U8zQIZFIio1K3jzFFxqNpiAUmY4OIUHwqFSqgsZSuwlZ89DR0ZFndKBPHHaPUCgUBoNRWlpaWlrq5ORUUFBga2srUwYNRebl5cVisfz8/DCvg1nG2dnZ2dlZcRug896bG4XOewPBVtXKlStjYmIevOL+/fsDywwlFJlMmddqBgSCO9hfOJFIlJiYOMiZQwhFJlMGAlEt2KP72LFjFd9gQSAjEezRvbGxkcFg+Pr6Ghoaojn//e9/lVgrCIQQsOW+bt26devWKbkqEAjRYE9mJk+ezOfzk5KSEhMTuVwudOyAvBtgyz00NHTv3r0MBsPa2nr//v2hoaFKrhYEQgTYk5mLFy8WFRWhz4BWrlzp7u6+b98+5VYMAsEfuU9zJC8CwBg1kHcG7NF9yZIl06dPX7ZsGYlEunjxItzfHfJugC330NDQSZMmZWRkAABCQkI++eQTpVYKAiEGWbl7eXnFx8dv27ZNklNTU3P+/Hm47g55B5CV+7Zt2wwNDeGTf8g7iazc0Wn6lStXfvnlF0nm8uXL4dI75B1AVu5oYCYWi3X37l00RyQSGRgYEGL85k3A4/Wn3d3B2LH96cZGh79uddZwAAB9ZErWOF/JGROrC0d1959Sa2rNMrFG0/pdHW41RZJi91wmI69Wllzrig35LwAAxno64LE+mDixv1BHB0hN/bsyc+dKkk4N5WbcFjT9XN+0zLL/1TctUY9PWa6kWIG9e6dO/5vW9s3VFm39EX5e0A2f2vSH9yEhCJCeB06bBkxM+tPFxeDZs/dq2v7xrK6DqvfYzk1SanJJNkXc7ydRbDWOM8oITTNesMeyq9B0lxY138FTcop3eZ5OrxBNlzMc/rbIYoGcnP60tjaYN09yxKPqEV3Qv5V5tZltg3H/tkKj+e3v1z1F0wgg3Rv/90g3ofapQWf/y1T1xpY1ZrZomt7N96h+TPujDYyiAgDA/PlAS6v/nOxs0NAAAKB1dDs0vaxg9P+XdXoE3hX5kivnOXzQrdXvijGWXcV4wUbTraOMS6z6XxfX6BP5lT6QnPLIzo1H1ev/8PQpKCnpT5uaAn9/gIWs3P/66y8AwLp1606fPi3JJEruV68CiQPrhg1/y72u7v2b8QyeAADQo6ElLXff0hybVhaavj1xukTuhrwXgfk3JcXuO/uIKP1N86x8OK6hFACgo0kBHqZ/y/3lSyC95eX06ZKkW+0Tj6pHaPqh3USJ3LV7hdJWyiwcJHIfxyqdUtL/mnSppZNE7mRE/D9WnJ3/lvujRyA21oonCOTw60yspeU+6+FtbVG/PxdHz0gid5sWlqQCraOMpeUeUPSnRIhdflI711ZV/V0BfX1puU8pvi9RVcoHH0vkbtzBkViRkbt3eZ7k+5bx/j8kcjfo4gbm39Rj5wDUhWr27L/lfucOyMwEAOj19LnqjpXInSbsku7MZ9YuErm/X/tsUnn/sPLEZry03KVPqTaz/VvuublA4oTk5iZP7iTMfTWSk5OlP+rq6k6bNg3z/Dehu7tbXoS9M3fKbxc14WvO0pAWsdxbQYGr2bVXs2vxNaqjSbmycaqCAneesE/fLsPXKADg0oapNPnOe49qX/x4rUje0WET/pWXlZFcb6b6ts5vz+fhbnTPggnutoaYhzQ1NbUk3zoAgLyFSNSXFADw/PnzjIyMTZs2ESF3Bdtnq8R5T6ZrcOFtdd4jxG10EOe9LvyDhwKFznsDwe6Rf//735L08+fP165di0O9IBBVM/gLAiYmJmVl+P/aQiDKB3t0X7FihST99OnTiZLbOwhkJIMt96CgIEl66dKlcNEd8m6APZmZPn16WVlZVFTUhQsXampqCNoyAQJRMtg6Xr16dVtb29q1a0kk0pkzZ548eRIeHq7kmkEguIMt9/T09KqqKnS3mlmzZjk7O0O5Q94BsCczNjY2XC4XTb98+RLd5xECGenIju6oW6qNjY2Hh8eCBQsAAL///vuyZctUUDUIBG9k5f7ee++hf6dO7X/0vXv3bqLemYFAlIus3KOiomTcO1Cg/x7kHQC6d0DUCAz3jhcvXty7dy8kJEQlFYJAiANjZUZfX//SpUs1NTXKrw0EQigY6+4UCmX8+PHu7u5+fn6SVysHumaLRKKIiAg+n29tbb18+XLMTA6Hs3XrVjRux+bNmxkMBoFNgUAGA/sx06ZNmzZt2qT4TDSq3uLFi8PCwhoaGtCI8jKZPB4vMDAQ3uZC3hKw5T4U1+yKiorx48eDV1H1ULnLZJJIpKamplOnTrm6ukocRP7666/i4mIAwPz58/X19TEr0NcnN6LVsEEQRHGwMZFIhLtRMFiEM8yohrgYpSBy33QiItIbAEAoFAoEcuPGCIVCIoz29PTI62ENDQ2Z172G75qNGVVPJtPMzMzV1dXd3T08PNzQ0BANQFlaWnrnzh0AQEBAgLxQZPKc+t4EBEEUa4ug75jyjQIARCJRL1lu6C8CjcpvLEGjSV9fnzyjJBJpELkP3TUbM6qeTKa3d79vaEBAQHl5OSr3VatWrVq1SnEbiAiCRyaT9fT0FBQgwo+ORCIpNqqjQ0iUFDqdrsB5j0rFP6QmAEBXV1dBY3UJsQmoVKriHpZGdmXGyMjIyMgoOjr6r7/+uvOK69evDzwTjaoHAGCxWGNfbSIgk/nbb78VFhaiHzGDs0IgygR7AFi5cqVQKEQDXgMAyGTyZ599JlPGx8fn1KlT0pH3bt26tWbNGunMGTNmHDt2LD4+3sjIyNfXd4AdCESp4B95TzrTxMQkLCwMj3pCIDgAI+9B1AgYeQ+iRsDIexA1Akbeg6gRMPIeRI2AkfcgagSMvAdRI2DkPYgaASPvQdQIDLmLxeLm5ubAwMDAwMAbN24QsbM7BKISZOfl5eXl48aNO3v2LPrxxIkTLi4uRUX4R3qAQJSPrNw3bNiwefNmycpjWlrakSNHtmzZovSKQSD4Iyv3goKCpUuXSucsXLgQhjOAvBvIzt1NTExkVh67u7t1dHSIsC0UCjEDoQHovIeHUei8N7jz3uzZs3fu3Hns2DFU9AiC7NmzZ8aMGURUVFNTU57ciVjsJ5FIip2kCHrCoNgohSI3Pt4bGtXUlCt3gjbs19TUVNBYIjzUAAAaGhryrjzwHyrb7LCwsIULFzo7O0+ZMoVCoWRnZxsZGSUlJRFRUQXykjzkwhfF2iJI7qoyqsAuQUbJZPLbZlQGWbnr6Oj88ccfBQUFBQUFPT09S5cu9ZcTkRUCGXFg/6h5enp6enpiHoJARi7wfRiIGgHlDlEjoNwhagSUO0SNgHKHqBFQ7hA1AsodokZAuUPUCCh3iBoB5Q5RI6DcIWrE8F8EHUooMswyEIiqGP7ojkYd27t3L5vNbmhowMzELAOBqIrhy72iogKNd4BGHcPMxCwDgaiK4U9mhhKKDLPMtm3bYmJiAAB37tyxsbHBvHhfDwGOXoi4ra1NwfFeIf7RgzQpJMVGhYIu3I0CALjtLwUacscyQRefCKOdvI42sly3wE4eIcGZBF38tjbsllKpVJlYd8OX+1BCkWGWCQkJQXfTtrCwkOdFtnb2hEVTh6R4Ho+nq6s7FE8Zuo6mHlWR/9gX/qNmutsOxSifz6dSqUNxoqFqaRjoaikoEOil7+kwpOjKXV1dir3jJGhpkI30FLkX+xgYxDCM5ThO/g8CgYBEIg0lSBuFTDLVpyooYGAAYtf794kHt4o6MQ/FQ5pEAmYGNHmeb4M77w0dNOqYl5cXi8Xy8/PDzNTU1BxYxtjYGB3yFUCjUGg6ilQi4QW5V1+fjovHJ4VCsTQekj9lu0afrq4uLs6XFAAsjYcUOK6jA9HW1sYrPKCF0ZCM8vkkMpksLx7o62JuSB9Ksa4uilgsptOHVPi1GP7c3cfHp7GxUToU2alTp2QyZT7iWG8IZBiQ5O0FMFJ48eKFvr4+Qf788mhvb8drdB86HR0dOI7uQ4TP5+M4ug+Rrq6ut250h0BGHCNe7j/++GNra6uSjYaHhyt/XTU2NjY/P1/JRhMSEtLS0pRsNC0tLSEhgYgrj3i5X758mcvlKtloYmIim81WstE7d+5UVlYq2Wh2drbyN8QtKirKzs4m4sojXu4QyNAZ8beqdXV1DAZDyXeNjY2NhoaGVKqiZWbcaW5u1tXV1dMb0gIiXnA4HA0NDQMDA2UabW9vF4lEg65WD4MRL3cIZOjAyQxEjYByf23EYrHyjba2tjY2NirZaFtbW11dnZKN8ng8Ho9H0MUJ2fhYOWRkZFy5csXe3n7Lli1Km7sLBIKDBw/u2rWLoD3vMbl3796tW7d6e3unTJmitLBw+fn5aWlp5ubmn3/+udLuUrKyspKSkjQ0NCwtLVetWoX/00NkZJKSknL8+HE2m33x4sXU1FSl2RWLxcHBwTt37uzu7laOxczMzNWrV/f29vb19e3evZvNZivH7pEjRzo7O9lsdkRExOXLlwUCAdEW8/LyNm/ejBqKjo7etWtXb28vviYoP/zwA85fIOKJj49PTEw8ePCgnp5eS0tLYmJiTk6OlZWVoaEh0abZbLa+vr6np2dUVJSfnx9BcQEkMJnM4uJiOp1uYGDQ3t5+48aNwsLC+/fv+/v7E7QFPgDg7t27Ojo6XV1dtbW1f/zxx/z587lcblxc3LRp04gzymQyi4qKqFQqnU43NzfX0tJ68OBBfX39Bx98gKcZfL89yqGtrS0oKCg8PLy0tHTnzp21tbWVlZUhISHEWUxPT//6668PHTrU09OD5uTk5OzevZvQMT4lJWXx4sXt7e0dHR09PT1CoRB9MzYyMjIvL484u2fOnKmqqhIIBCEhIfv370czw8PDnz59SpDFlJSUiIgIsVjc1ta2fv36yMjI7du319fX7927l8/n42hoRN6qGhoafvHFFxQKZf/+/SEhITY2Nvb29ubm5gSZYzKZjx8/3r9/v4WFRXp6Oprp7e09d+7cgwcPEhTR6ebNm5WVlVu3bkUFp6mpqaXV/0Z0V1eXiYkJEUZRqFSqUCjU1tbevn17W1tbRkZGT09Pd3c3Ee9sAQDq6+vj4uK+/vprEonE4XBMTU29vLxCQkIsLS319fURXBfKR9JkpqSkxMDAAL19sbGxYbPZ1tbWDx488Pb2vnjxoqWlpYODA+5Gy8rKYmNj9+/fP3DiZGFh4eTkRNAjGAqFMmvWLAsLCxMTk5MnT/r6+gIA9uzZw2QyPT09vb29cbdYVlYWFxfn5eVVXl6upaVlZWVFpVKnTZv24MGD+Ph4Pz8/Dw8P3I0CAPT19bW1ta9cuWJiYnL58uXvv//exsYmKysrNjZ20qRJzs7OONoaMY+ZOjo6VqxYYW5uPnfu3BkzZlAoFCaT6eXllZycnJ+f7+fn98UXX+BuND09fcKECUlJSVwu96OPPrpw4cI333wjEomioqLCwsJwN4dSUlIyduxY6bWmhw8fXr58ed++faNGjUIQhKAJ9L179549e9bT0zNmzJhbt26RyWRLS0snJycnJydHR0cijPJ4vIaGhjFjxujo6KSkpMTExJw9e9bIyAh3QxJGzOiura3d2NgYEBBQU1MTHR1No9H8/f2Tk5OXLl1qYWFBUGzAiIiItLS01atXs1ismJiYQ4cOmZqaGhoaFhcXT5o0iQiLHR0d3377bVZWFoVCGTNmDOp+Zm5ujo7x06ZNI+jmGL0nXrduXVlZWWpq6sqVK4ODg62srNrb29vb28eOHYu7xczMzJMnT6JL+xMmTHBwcKDT6fHx8WgMPNzNoYwYuQMArK2tL1y4sGPHDgRBbt68yWQyhULh+PHj0c0OiMDU1LSlpeXWrVvBwcF9fX1ET5zu3btnZWX1/OuOm7cAAAZwSURBVPlz6W+1jY0NiUQyNzefMGECQRMnJpN54cIF9GHCxIkT29vby8rKfH19DQwM7O3tidA6h8O5evXqoUOH/Pz8xo0b969//UtTUzMgIEAgEFy+fNnf35+ogIQjSO6jRo0qLi4uLCzkcDh79uzx9PSkUqmurq64G+Lz+eh9obm5eWZm5ueff/7LL78EBwezWKzo6GhbW9sFCxbgbhQdX/38/GxsbKS/1UlJSXp6ejY2NgS9HIbeEy9cuPCXX37x9fXV1taeOHFieXl5ZmYmQb9gAICHDx8aGho6OTmJxeKff/75o48+io+Pp9PpAQEBVlZWRLwc1g+Oqzy409fXJ5PT1NS0bNkykUhEnNGioqIFCxacPn26ra0NQZDS0tLY2Fh0kYTD4Tx+/JgIo5I1R/Tj8ePHIyMjw8PD+/r6mpqa7ty5Q4RRtHurq6vFYjGCIAUFBVu2bOFyuehRglqal5cnEolKSkpOnjyJ5jQ3NyMIcvXq1cOHD0t6gCDe3oXItLS0b775Zu3atTU1NZJMc3NzT0/P3Nxcgoy2tbVdunRp6tSpHA7n7NmzsbGxDAajtbV19OjRa9euLSgocHNzw92o9JpjR0cHAGDRokX3799fv349mUw2NzefPn067kYl3QteBW328PBYsmRJaGgoWgciWspkMrOzs8lkspOTE5vNfvz4MQDAzMwMAFBZWblx40Z9fX3cjUrzlk5mmExmaWnpvn37nJ2dT548OXv2bMkhW1vbyspKR0dHIuzSaDQjI6Ps7GxjY2NLS0sPD4+YmBgNDY2KiopZs2bZ29sTYXTgmqOxsXF9fT0AgKDtG+R1L6H3xEwms7KycsOGDSQSCUEQT0/Pc+fOVVVVcbncq1evuru747vmiA2hvx3Do6GhYdmyZegDy46OjnXr1im5AgUFBdu3bz9+/PiVK1cQBMnPz79x4wYRhnJycnbs2HH06FHJyyGSGUVzczNBRgft3qamJtyNslgsidGysrKffvoJQRCRSJSbm5uYmFhZWYm7RUzeOrkXFxefP3/+9u3bu3fv5vF4Bw8eJGgSqRgZxRPBo0ePfv75546ODhaLJWN6y5YthL6SpZLuTU5O3r17d1FR0c6dO3k8nnKMyvDWyR1BkHPnzqGKR6ewqqoGqvjs7GyCrn/kyJGGhgYEQerq6o4fP75z587S0lL0EBHjqwwq6d7k5OTPPvuMw+Eo06g0b9FTVUTqeeGvv/6qra3NYDD+/PPPPXv2SF4XUTJsNpu4V3ESExPLyso0NTVfvHgxb948LS2tyMjIM2fOEGRuIHfu3FF+96akpGRnZ6vqf/q23KqWlJRs2rSprKyso6MDfVUjLy+vu7v7/fffv3DhAnHPHRRDqB+0s7Oztra2s7PzokWL0PvUjIyMWbNmEfeSrQx2dnZisVjJ3evg4ED0syQFvBWje3NzM5/PT09P7+zsdHNzKywsbGhoYDAYxcXF//jHP6ZMmULcc1OVIxaLq6qqrKysoqOj7e3tP/74YyVXoLq6WvndW1FRQcRj6UFRvdxTU1NTU1NJJJK3tze6DeKyZcsAAI2NjYWFhQiCzJkzR7U1JJSurq7o6Ojm5uaPP/546tSpqq7OO46K5S5ZiwUAbN++/ciRI+isHVU8BIIvqnyq2tzcHBcXFxwcTCKRent7uVzu7t273dzchELhhQsXVFgxyLuKKuX+3nvvrVix4uDBgx0dHYcPH161alVISMjt27cBAES8hQeBqH7unpGRERkZuWnTJjS8B4IghYWFEydOVG2tIO8kqt9nZtq0aQCA5ORkDw8PHR0dEokEtQ4hCJWN7ikpKffv3x83btySJUsAABkZGbdv3w4NDSV6KwuIOqMauVdXVz98+DAoKCg+Pp7L5QYHBwMA6urq5MWdhEBwQQW3qgiCMJlMdPPBoKAgfX19dB0Gah1CNMoY3Z88edLQ0DB79uzExMT09HRfX9+FCxeePHnSwcEBfdMajusQ5aCM0d3V1bWioiI5OZnFYp04cYJCoVy8eHHjxo0VFRVMJhPAcR2iLJQhdxKJtHHjxqqqKnQj46CgIBqNhiqe0E1FIBAZlHeriiCI9ASmuLjYxcVFOaYhEBSlrszIKB4CUTLKXohEECQvL4+ILQ4hkEFR/UsEEIjSeHv3mYFAcAfKHaJGQLlD1Agod4gaAeUOUSOg3CFqBJQ7RI34/4imu93dqB2MAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 250 -h 200 -u px\n",
    "fviz_pca_contrib(pca, choice=\"var\", axes=1, top=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How are the contributions are defined?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I recreate the calculation of contributions step-by-step for benefit of the reader (including me during the write-up):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) The loadings are normalized (scaled to unit vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import isclose\n",
    "\n",
    "assert all(\n",
    "    isclose(pc_squared, 1) # equal 1 or very, very close\n",
    "    # (rtol=1e-05, atol=1e-08) - to allow for rounding errors\n",
    "    for pc_squared in loadings.pow(2).sum()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) \"sdev\" gives us the square root of eigenvalues, which can be used to calculate the variance explained by each principal component:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>variance_explained</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PC1</th>\n",
       "      <td>0.479502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PC2</th>\n",
       "      <td>0.106114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PC3</th>\n",
       "      <td>0.045693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PC4</th>\n",
       "      <td>0.034823</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PC5</th>\n",
       "      <td>0.030572</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     variance_explained\n",
       "PC1            0.479502\n",
       "PC2            0.106114\n",
       "PC3            0.045693\n",
       "PC4            0.034823\n",
       "PC5            0.030572"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "variance_by_pc = stdev.pow(2)\n",
    "variance_explained = variance_by_pc / variance_by_pc.sum()\n",
    "variance_explained.columns = ['variance_explained']\n",
    "variance_explained.index = 'PC' + variance_explained.index\n",
    "variance_explained.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In documentation the contributions are defined as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\text{contrib} = \\frac{\\text{variable.cos}^2 \\cdot 100}{\\text{total cos}^2 \\text{of the component}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the loadings are cosines of rotation of variables required to form the PCs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "absolute_contributions = loadings.pow(2) * variance_explained.variance_explained\n",
    "relative_contributions = 100 * absolute_contributions / absolute_contributions.sum() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hnRNP A/B    0.154683\n",
       "IMB1         0.153810\n",
       "STAT3        0.153632\n",
       "PTK6         0.153442\n",
       "FGF9         0.153075\n",
       "Name: PC1, dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relative_contributions['PC1'].nlargest()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which agrees with the values on the plot above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**But how are the contributions for more than one PC calculated?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are just means of contributions (i.e. means of squared loadings, weighted by variance explained by each of the principal components):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KREM2        0.127032\n",
       "hnRNP A/B    0.126739\n",
       "SPARCL1      0.126043\n",
       "IMB1         0.126038\n",
       "STAT3        0.125892\n",
       "dtype: float64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pc_1_2_contrib = absolute_contributions[['PC1', 'PC2']].mean(axis=1)\n",
    "(100 * pc_1_2_contrib / pc_1_2_contrib.sum()).nlargest()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please, compare with this plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAADICAIAAABOPGvMAAAgAElEQVR4nO2deVwTR//HJxeQcBMOCTcIUkCtUhClFsvVimct9T7L81MqlMd6ggrWE285PPDAIh7leUpFKoeoCGKRQ1E8ilwqRpAbCQRyELK/P1ZjnrCJiNmgZt5/8JrMzu5nZvNhMju7810CgiAAAlEOiINdAQhEcUC7Q5QIaHeIEgHtDlEiBt/u5eXlM2fONDc3p9PpkyZNevDggdwl0tPTXVxc5H7YgXH69Glzc/MzZ87gLYTZ6rq6OgqFgre0iOHDhxNeY2Njs3z58q6uLoWp92WQ7V5UVOTm5mZjY3Px4sWcnBxTU1NPT88XL16863HMzMxKSkr6HpxAIFRUVAzsUAPb/a3s379/8+bN8+fPf/9D4VRD2WCeahns27evoaGhrq7u6NGjt27d8vT0RCcD37PyJ0+e/PnnnzE33bp168svv9TU1LS2tt69e7f4pkG2+6pVqwIDAyMjI0eOHDlixIijR49aW1sfP378/Y/c3NzMYDASEhKMjY0HdgQbG5v32V0aXC7X2tr6/Y/z/g1UDNra2kZGRgwGw9vbOzs7u6KiIiUlBbzH6S0tLf3111/XrVuHubWrq2vatGne3t5Pnz49depUVFRUfHz8m83I4FFZWQkAqK+vF89kMplVVVVouqCgYOzYsRoaGo6Ojv/5z38QBGlpadHR0bly5cqMGTNMTEzmzZvH5/M9PDyIRCKdTr9w4UJLS4u+vn5UVBSJREpNTQUACASCtLQ0Z2fnX3/91djY2MzMbNOmTb29veXl5RoaGiJdb2/vM2fOSBwK3b3/NZFoYN+9pk6dSiKR9PT0zp49K15y0qRJGzZsQNNdXV00Gq2oqAhBkIsXLzo6OqqqqpqYmERGRqK6fRuIWRKz1bW1tWQyWXT+fXx80F7w+PHjaGZkZKSpqSmNRvv666+fPHki0SLx84PZQAmcnJxOnDghnrN48eKlS5eiDUErj57JM2fOTJo0yczM7MSJEzExMT4+PgYGBocPH+57zBMnTgQGBjo4OAQHB/fdeuPGDTqdLhQK0Y/r16//4YcfRFsH0+7p6em6urrStjY3N2tra0dHRzc1NV24cEFdXf3WrVstLS0UCiU8PBxBkK6uLiMjI/Qsm5qa3r59G0GQlpYWMpk8e/bs2tpa0QlNS0sjEonBwcENDQ3Xrl0zMjKKi4vDtLvEodDd36kmsuuPIIi9vX1ubq5EY8+cOWNnZ4emk5KSHBwcEATh8Xjq6urR0dENDQ3p6ekkEqmqqgqzgZglMVstsjuXy7WwsNi2bVtTU1NWVpaurm5aWtq1a9c0NDTy8/OfPXvm5+c3b968vt+L6PxIa6A4fe2+efNmHx8f5H/tTiQSd+7ciSBISkoKgUDYtm0bmtbS0pJmj8DAQEy7s9lsJpOJpgUCwfjx4/fs2SPaOph2P3/+vKGhobStJ06ccHNzE31ctmzZv//9b/Qc1dbWopkzZsw4ePAg0sejT58+RcROaFpampaWFpfLRffav3+/h4dH/+3+TjWRXX9Eit07OzupVOqDBw8QBJk+fTr6DfF4vNLSUrRAWVmZhoZGYWEhZgMxS2K2WmT3lJSUkSNHiioQERERHByclpamqal5+fJl9JidnZ19vxfR+ZHWQHH6aXcAwIsXLxAEaWpqkkijv119kWZ3ERUVFT4+Pp6enm1tbaLMwRy7Dx06tKmpqb29XTzzr7/+CgoKAgDU1tYOHTpUvDCTyUTTRkZGaIJIxK6/iYmJRI6ZmZmqqiqatre373s1LBQKpdVzYDWRsVdfNDQ0Jk+enJyczGKxrly5gl7IqqioFBUVeXl5jR49Ojw8XFxCooHSSspo9ZMnT8rLy41fc+zYsfb2dj8/v71794aGhurq6i5atAgdbQ7gtMigpqbGysqqb76enh54fRrF0wCA+Ph4TU1NTU3NKVOmYB5TogCHw1mzZo27u7uXl9fly5d1dXVFJclvrR9+ODg4DBs2LDY2Njw8XJR55MgRBoMBADAxMcnKyhLlP3nyxNTUFE0TCATZR+5rvufPn/N4PPS7f/TokaWlJQCgt7cXQRD0aM+ePZN2tIHVRMZemMyZMyciIsLKysrT03PIkCEAgKysrA0bNuTm5jo6OgIA0EzMBkoridlqFGNj41GjRhUUFKAfm5ubhULh8+fPv/7666VLlzY0NMTGxnp5ebW2tsroU96pgQAANpudkpLyrlMRAQEBAQEB/SyA/jyqqKiUlZUZGBhIlBzM3p1EIkVFRUVGRkZERDx8+LCysjIsLOz69etr164FAEydOvWff/45fPhwe3t7RkbG6dOn582bJ+NoLBZLxtaurq6VK1c2Nzfn5eXt3r07ICBAT0+Pw+EkJiby+fzY2Ni6ujpph3rXmgxsr4kTJzKZzJ07dy5ZsgTNaWpqotFohoaGTU1NkZGRjY2N9fX1mPtKK9m31aJdvv322+rq6ri4OBaLVVhYOGLEiLy8vMzMTDc3t3v37lEoFA0NDRqNhul19Pz0s4GdnZ2tra1NTU15eXk+Pj52dnbff//9W8/egMnKyiouLt61a1dHR8fjx48fP37c0NDwZrOM0Y9iKC4u9vPzMzIyotPpPj4+xcXFok35+fljxoxRV1d3cHBISkpCxAZ8aAF/f390xLx69WpNTc3z58+LFxCl8/LyQkJCVq9ebWhoaGlpuW/fPnT3PXv26Ovrm5iYhIWFhYWFpaenSztU/2siTt+9ECljd5RFixYZGBiIZni4XO6sWbM0NTVtbGx27dq1d+9eOzu758+f920gZsmsrKy+rRafmbl9+/aXX36prq5uZWUVExODIAiPxwsICKDT6TQazc3N7ebNm30rKTo/0hoojpOTk8hplpaWy5YtE10PSIzd0WuMvul3Hbtv3rxZ4h/A399ftJWAfMwPAPP5fIFAQKPRBrsiH1BNPiVevnwpPvJ+fz5uu0Mg78TgPzMDgSgMaHeIEgHtDlEioN0hSgS0O0SJgHaHKBEfvd0HZSJVqUQVr4uf6Edv95cvX/b29ipYlMVi9fT0KFi0s7OTx+MpWLSrq4vD4ShYlMPh4LTG76O3OwTSf6DdIUoEtDtEiYB2hygRg7m8QwYZd5gZd96+NAYAIBAISCTSWxd8AACGMXR+9nOSUeBpY8e1h/2K+cHj8SgUirR1D+IYaKlNdbHszzEhCuADtXtlPauqXtZyjQHA63nLBE5+ReN/bz6WryhVhSzb7neftkSn9yuSlFAoROMTvbWkGoUU+68vKSSp/41PGzu2JN/pz2Rf/0VJRMLehWN1NVSlFejg8HdfKBX0vl0UnWojkUhvLUkkEIInOjL01N9aEmUw7c7hcKStEBUIBHKXEwqFsqe3+Hy+3EURBJEtymxsr3/ZLXfd9g42TUWqXepaWC/a5D/T19jGUiFIdd6L1q5b1c1yF31S36Yt5V+MQqGoqKiI5wym3alUqrRNZLL8K0YkEtXVZXUDEqdGLhAIBNmioqXT8oVGo6mrSj2Hamry/wcDAFCpVBmNpXZLXfz+Pqipqck+w+LAS1WIEgHtDlEioN0hSgS0O0SJgHaHKBHQ7hAlAtodokRAu0OUCGh3iBIB7Q5RIqDdIUoEtDtEiYB2hygRGA/NCYXCixcvFhUVlZaWCoXCUaNGubq6Tps2rT+rGSCQDxlJB2dmZk6dOvXevXtjxozZv39/dHT02LFjHz58OHXq1MzMzEGpIgQiLyR795aWlpSUFPH3iA8bNmzq1KkCgeD3339XbN0gEDkj2bsvWLBA3OsCgSArK4vP55PJ5AULFii2bhCInHnLcNzf3//OnTuzZs1STG0gEFyRHMw8evTos88+E318/vz5hQsXxowZ03dPgUAQExPDZrPNzc0XL14svmnbtm2rV69WU1OTUQYCUTySvfuRI0eCgoJE75t1cXH517/+5eDg0HfPwsJCBoMRERFRX19fW1uLZrLZ7DVr1hQXF8soA4EMFpK9e0xMTFVVVWhoqIWFxZo1a+Li4h4+fIi+nFaCqqoqNN/a2rqqqgp9hayGhsauXbsiIiJklGlpaeno6AAAmJiYSFuCjVMAWNnBU2W8OPtjFO3tlRotAydRoVAoo7GKFyUSiRIhQzDcZmtre+rUqeLi4iVLlnh4eAQGBmKGGenu7tbX1wcA0Ol0NpuNqYFZZufOnSdPngQAXL161cLCArOieIS67e3tlXghvQRcLlfuogiCyBbt7sYlKACLxeqRHnhD/PuSIx0dHe1kqRFTOjpwaSmbzW5vx+40qVSqxKs/McrFxcWdOnVq0aJFf/7558WLF/39/efOnTt79myJYjQarbW11draurW11dDQEFMPs8zevXv37t0ruw1qav2K5vVOkEgkOp0uowCN1iZ3UQKBIFtUQwMXE+jp6ckIvKHVjktHq6OjQ6drStvKFso/rgkAQEtLS/YZFgdjZubo0aMFBQVoBzxlypSUlBTMwEC2trY1NTUAACaTOXToUMyj96cMBKIwMOxOpVLz8/NF4X5IJFJAQEDfYm5ubnV1dXv27DE0NDQzM6usrDx48KDsMnKvPQTyTmD83p06dSouLg7t3WXtSSavWLFC9NHOzs7Ozg5Nb926FbMMBDK4SNqdw+HY2tru27evb1EOhyMjzB0E8uEjafcVK1YwGIzAwEAjIyNRZmNjY1xcXF1d3bFjxxRbPQhEnkja/ejRo/n5+UuWLKmpqUHnEFtaWiwtLTds2ODu7j4YNYRA5AbG2N3d3T0jIwMAgN4HRe8NQSCfALLiSkOjQz4x4AIliBIB7Q5RIt5u90F5NzkEggfYds/OzkZvD3333XdaWloJCQkKrRQEgg/Ydl+2bNnkyZMLCgq4XG5FRcX27dsVXC0IBA+w7c7j8by9vdPT06dPn85gMHp6ehRcLQgED7AnIkeOHLlo0aK8vLy7d+9u2rTJ2NhYwdWCQPAAu3c/e/asm5vbhQsXdHR0SCQSDLkB+TTA7t1JJBKDwSgvLy8vLx82bFhJSYmlpaViKwaByB9su//44488Hs/a2hr9SCQSv//+ewXWCgLBBWy7CwSC1NRUBVcFAsEbbLsPHTq0vb1dR0cHX/HvvwfPnr1Kr10LZs58lc7NnbsiaFI3HwDAo6isCtgt2mNVSpRVYw2aTneZmOn8zasK1z9ekRorKrbi//YKSK+aFph53OnZPwAAFTIRdC8Fq1a9KvTkyRtFAMClS6LkwmtnXStvoeliuy8SPeejae3uju2JEaJi22euq9d7dRH/fX7K1w+uo+mHFg5xE5eiaVKvAHzxxRuV48fBqFGv0idOgLi4cd18GxbnqZHlvu/erIPZF79OtefV4vRj3wTctxqOpr/65++ZN5LRdJO2wZY5G0S7bD2zWZf9Ek2fnTAbgFdnBmRlgQ2vi+nogKtXRbuEJu81bXkVDSVl7NTskZ5o+rPn5UHpcWgaIRB+XnZAtEvIxUN2dVVo+vIo77/GTEbT5s3P1/65z/hPDUAmAQDAtWtAS+u1TCgqaizonaM//HePVyG69Dtafz23VXTkX+dubNHSR9Oz8/74siwfTd+1HhnvuwRNU/ncPSdDRbvs/n4l08D81YdDh8Bvv71Kf/EFiIsDWGDbva6ujsFgjB07Vk9PD835448/MEu+D8jKlUC0Dn/YMCC6fevg8PfiFcVVTQAA4f+GHb7oOkmd+2rhrMhqAIAG3SHxPotFH3uJbxbhXxnlVWTnAgAw1Kb+MnnSGxVDQxAZ+ebQGhoI8mpp9nWn8Q8sXsUaadPUExXpVqGKq7Rr6IrSRfauT4ZYoekOdS1RvpBIQsRVLC3fVMDTE1hYlFc3ny962qWmLt7MRM95xNdhKp4Zmovyy8zsRRXgUVTFd/nP+B8oAj6aZhqaIwjy6nb4iBFvmkmhoOropgtuU6g8Drqljs4QHaqObvKmmf8bhCLji29vOLx6DrxR982KiBYterzP4l8mDzfUpgIAgJram2bOnQu8vAAAL1mc/Lw3K+47qZriJ7OT+uak5TuMrTCxRdPtGm/6XD6ZIr4L+u/xqqW+vuD1YjqgqytSlwiiQcB8RiA/P18iB4+H3blcrrQnFA5dLr9yX87BCMzo6rFLMMKhiUi6+TTp5lP5ilJVyL+HfCWjwJX7Lw5dLpevKADgXIgHTXrgjbs1bZuTS+UuGrtkjBldXdrW561dP/9WJHfRTf6fj7LUw9xEJpPFA54Cab27u7t7VlZWbm5ub2/vhAkT/Pz85F5LAICampq0TdLCLb0PBAJB9uJDiVMjL2SLqqjgEo5CTU2NKj3wBk6iqqqqMhqrqio1BM37oKKi0v81pdjz7ps3b46IiGAwGObm5lu2bNm8ebP8qgeBDBrYHcDp06fv37+PRmD68ccfR40atWnTJsVWDAKRP1IfABaN8eE7aiCfDNi9+7x587y8vBYuXEggEE6fPg3ju0M+DbDtvnnz5jFjxuTm5gIAQkNDJ0+erNBKQSD4IGl3FxeX5OTk1atXi3KePn166tQpPObdIRAFI2n31atX6+npwUh3kE8SSbujw/Rz584dOnRIlLl48WIYUwnyCSBpd/TFTEwm89q1a2iOQCDA/eEZCEQhSNr977//BgAsX7788OHDokxod8ingaTd0RchLFy4sLCwUJSprq4+YcIERVYLAsED7InI5ORXT5k2Njbm5uaGhIRAu0M+AbDt/pvo0WEAGhsbf/rpJ0XVBwLBkbc/IGBgYFBRUaGAqkAgeIPduy9ZskSUfvjw4eeff66o+kAgOIJtd39/f1F6/vz5cNId8mmAPZjx8vKqqKg4ceJEYmLi06dP8VhsAYEoHmwfL1u2rLW19aeffiIQCEeOHHnw4EF0dLSCawaByB1su+fk5Dx+/BhdzObr62tvb9/X7gKBICYmhs1mm5ubL168GDOzpaVl1apV6Puyf/nlFwaDASCQwQPb7hYWFiwWC30V2cuXL9GEBIWFhQwGY/bs2ZGRkbW1teibbSQyOzs7/fz84OPykA8ESbujy1ItLCxGjx49Y8YMAMCff/65cOHCvntWVVU5OjoCAKytrauqqlC7S2QSCIQXL14cPHjQyclJdKPqwoUL6C3boKAgUWAPCQQC+S/jFQqFmG+7F8Hn8+UuiiCIbFEejyd3UQBAd3c3EEiNRMDlcvEQ5XA4XV1Sp7Y5HA4eolwuV9oZplAoEovQJe0+ZMgQ9O/48ePRnI0bN2I+M9Pd3Y32+nQ6nc1mY2YaGRk5OTmNGjUqOjpaT09vxIgR79UyCOT9kLT7iRMnJJZ3oPQdkNBotNbWVmtr69bWVnR03jfT1dUVzff09KysrETtPn369OnTp7+lWjjMBRGJRHV1qVFQAD7hKAgEgmxRVVVVGVsHDI1GU5ceeENNrVvapveBSqXKaCy1W4iHqJqamuwzLM7Al3fY2trW1NS4uLgwmcxx48ZhZv7+++8ODg4jR45kMplDhw59p2ZAIHIHY3lHW1vbjRs3QkNDMXcQ4ebmdvDgwT179hgaGpqZmVVWVl6+fDkwMFA809vbe//+/cnJyXQ6fezYsbi1AgLpFxi/d9ra2mfOnJk1a5aVlZWsPclk8R8BOzs7Ozs7AIB4poGBQaR4hEQIZFDBsDuJRHJ0dBw1atS4ceNEoyK4NBvyCYB9NRMSEhISEqLgqkAgeINtd7g0G/JJApdmQ5QIuDQbokRgL82Oj4/Pzs4W3Womk8nwVWSQTwD45j2IEgHfvAdRIrCfX0PfvKfgqkAgeDOYb96DQBQMtt2XL1++fPlyBVcFAsEb7MGMu7s7m81OS0tLTU1lsVjwHhPk0wC+eQ+iRMA370GUCPjmPYgSAd+8B1Ei4Jv3IEoEht2FQmFDQ4Ofn5+fn9/Fixfxi+zO5XIRBMHchEfgDQRBZMd+6OnpkbsoeFvACTyifQAAuFwuUSg18AZOojwej8ORuh4cpxAjfD5f2hkmk8loaLA3ORIlKisrp0yZMnv2bHQ2JioqKigoKC0tDY+YGTIW4eMRiYBAIKipqckogFMoTNmiEt+HvFBVVVWTHokAJ1EVFRUZjVVRwaU3oVAoss+wOJJnJDg4+JdffgkMDEQ/Zmdn//e//125cuXVq1flWUcAgNjVsMKQrYhTfQZLVMaRlUdUAslZl5KSkvnz54vnzJw5E77OAPJpIGl3AwMDiZlHDofT/x8LCORDRtLuEydOXL9+vVD4KuATgiDh4eHe3t4KrxgEIn8kx+6RkZEzZ860t7f/8ssvSSRSQUEBnU5PS0sblMpBIPJF0u5qamp//fVXSUlJSUkJn8+fP3++h4fHoNQMApE72HNVzs7Ozs7OCq4KBII38HkYiBIB7Q5RIqDdIUoEtDtEiYB2hygR0O4QJQLaHaJEQLtDlAhod4gSAe0OUSKg3SFKxMCXqwkEgpiYGDabbW5uvnjxYsxMzDIQyGAx8N69sLCQwWBERETU19fX1tZiZmKWgUAGi4HbvaqqCn3fgbW1dVVVFWYmZhkIZLAY+GCmu7tbX18fAECn09lsNmYmZpnVq1efPHkSAHD16lULCwvMg/fycQjSgAhbW1tlbO/hceWuSSERZIvyuN1yFwUAsNpfcslS+zJuNxsP0a7Ojlai1JAeXZ3yP70AAG43u7UVu6VUKhUN/Chi4Han0Witra3W1tatra2GhoaYmZhlQkND0WjaJiYm0mJd/DRxxKzx/XJ8Z2enurp6fyL7aahRNKmyAk7M8dDyGWXZH1E2m02lUkkkqYFcRFBVyDrqKjIK+LloO9sy+iPa3d1NoVD6EzNDhUyka8paXuymo3OSoS8lxs//wOVyCQSCjBApIkhEgqE2VUYBHR2QEOTRK3y7Ko/HQxCkPyukCQRgpEOTFoegrysGbndbW9uamhoXFxcmkzlu3DjMTAqF0reMvr4+2uXLgEYi0dRkuUREG7FHW1ujP857KyQSyVS/XwFY2sm96urqconWQgLAVF+zPyU7OhBVVdX+OK8/mND7JcpmE4hEokQfOWCM9TT6U6y7myQUCjU0+lX4nRj42N3Nza2urm7Pnj2GhoZmZmaVlZUHDx6UyJT4KMd6QyADgCAtbN3HQltbm7a2tlx69/7T3t4ur969/3R0dMixd+8nbDZbjr17P+nu7v7gencI5KPjo7f71q1bm5ubFSwaHR2t+HnVhISE27dvK1g0JSUlOztbwaLZ2dkpKSl4HPmjt/vZs2dZLJaCRVNTU+vr6xUsevXq1erqagWLFhQU3L9/X8Gi9+/fLygowOPIH73dIZD+89Ffqj579ozBYCj4qrGurk5PT49KlTXNLHcaGhrU1dU1Nfs1gSgvWlpayGSyjo6OIkXb29sFAsFbZ6sHwEdvdwik/8DBDESJ+LjtXldXV1paKopX/AkzKG0cFNHm5ua6ujqcDv4R2/3vv/8+fvz4zZs3d+zYoTDR3NzcpUuX7tq1C6cXOWHC5XI3bdrE5eLyiNUHJXrjxo2YmJjY2Ficgk5/xHa/cOHC2rVrly9fbmBg8OjRIwUoZmZmlpaWbtmyxcTEJCcnRwGKKKqqqvX19Vu3blWk+RQveuPGjbNnz27atGnHjh1FRUUNDQ1ylyD9+uuvcj8o3uTk5GhpadXX1/P5fEtLy4qKCjqdPmTIEFxFk5OTU1NTt2/frqmp2dTUlJqaWlRUZGZmpqenh6suAKC+vl5bW9vZ2fnEiRPjxo3D6ZVpgyuamZlZVlamoaGho6PT3t5+8eLFe/fu3bx508PDQ54vdUI+QlasWBEUFFRbW9vU1IQgyK5du9ra2oRCYWpqak9PD06ira2t/v7+0dHR5eXl69evr6mpqa6uDg0NxUkOk6Kioo0bN3I4HPwkcnJy/u///m/nzp18Pl9hohkZGbNnz25vb+/o6ODz+TweD30GOC4u7tatW3IU+igHMwsWLDAxMTlw4AD6QHNPT4+Ojk5cXFxXVxd+nZCent6cOXNIJNKWLVtCQ0MtLCxsbGyMjY1xkquurj5w4ADyv9PErq6uU6ZM2b59Ox7vnQVSRmt4i166dKm6unrVqlVbtmxBEIRCoaiovHr2u7u728DAQI5aH9Nghs1moyfC2Ng4Ly/vhx9+OHTokLOzM5PJTEtLMzExmTNnjtxFHz16pKOjgz5xaWFhUV9fb25uXlhY6Orqevr0aVNTU1tbW7mLAgD09PRqamquXLni5uYm/mtuYmIybNgwPO77VFRUJCQkbNmype9oDT9RAACJRPL19TUxMTEwMIiNjR07diwAIDw8PDMz09nZ2dXVVZ5icvylwJWysrIZM2YcPny4tbUVQZDy8vKEhAS0V6ipqUlJScFDlMVizZgxIygo6NKlSwKBAEGQjIyM5ubmhISE4ODgc+fO4SEqEAjQn3IEQc6dO7dv3z6hUIh+TEpKamtrw0P02rVrLS0tCQkJihytlZWViYZMKCUlJStXrmSxWAiCiFotRz6O3j0vL8/BwYHD4TQ1Nd29e/fx48djxozJy8sbN26ck5NTVVXVpEmT8NBVVVWtq6vz9PR8+vRpfHw8jUbz8PBIT0+fP3++iYkJHi8kvHHjxu7duy9evKiqqjp06NDhw4ejffyYMWNOnz5NJBJdXFzkLgoAiImJyc7OXrZsGZPJPHny5M6dOw0NDfX09MrKysaMGYOHYkdHx7///e/8/HwSiWRlZYWOS42NjdE+fsKECXiMSz8Cu2dmZj569Mjd3d3Z2bmystLQ0HD06NEnT54kk8lVVVW+vr42Njb4qZubmycmJoaFhSEIcunSpczMTB6P5+joiEZYkC+NjY3nz5/funXrd999Z21tjY5hUMcfPXrUysoKj9EaiqGhYVNT0+XLlwMCAnp7exUwWisuLubxeN98842oK7GwsCAQCMbGxiNGjMBr4PSB2z0zM7O6ujo4OJhAIDx//vyrr74qKip6+fJlcHCwmpoakUi0s7PDtQJaWlplZWX37t1raWkJD64i2KEAAAg1SURBVA93dnamUqlOTk54aJWUlAwZMgRtEYFA6OjoOHTo0NixY4cPH25ubu7l5SV3RczLoYCAACaTGR8fb2lpOWPGDLmLgtfTjnPnzhXvStLS0jQ1NS0sLHB8DE7uwyM5wmQyFy5ciE6BVVRUrF+/ns1mIwhy7NgxnMbNCIL09vZK5Lx48WLhwoXo2B0n6uvrq6qqKioqoqKixPN//vln9FoFD+7fvy/tcqilpaW0tBQnXdG0I4IgBw4ciIuLi46O7u3tffHixdWrV3ESRflwe/eGhgYOh2NsbHzu3DkDA4OzZ8+GhYWh6xctLS3JZDIek4DZ2dmRkZGZmZmOjo66urpopqam5vPnzwEAOK0uz8rKSkxMvH37NolEqq2tVVVVFQldv37dy8sLj8ebW1tbjxw54ujo2N7e3vdy6J9//nF3d5e7KHg97Thz5kz0h8vOzu63337bvn07iUTS1NTEY4gozgdqd5EDDAwMbGxsoqKiduzYgY7nEhMT6+vrPT095S6amZlZXl6+adMme3v72NjYiRMnijZZWlpWV1fjMXDKzMysqqoKDw/39fVNSEgICwtLSEioqqri8Xh//PGHs7Ozo6Oj3EUBADQajU6nFxQU6Ovrm5qaKuxySGLa0cfHp6mpCeDWlUigiNvR7wo6Xt+3bx8AYM2aNXv37gUAREVFhYeHJyUlUSgUPAaUdXV1SUlJR48eVVFRYTAYErdUjIyMJk+eLHfRhoaGpKSkI0eOEAgEPp/PYrF27do1adIkCoVSW1s7bdq0oUOHyl1UxOjRowEASUlJAID79+9HRESUlJTgtCixuLj4woUL+vr6ISEh6CU4qr558+bAwMCKigo8RPvywS3vaGhoWLdu3ZEjR2g0Gp/PDwoKMjIymjp1aktLy/nz5728vHC6l3Tr1i0Gg3H9+vV169bFxMRMmjRp5MiRchfqS25u7pUrV9atWxcbG+vt7e3o6BgdHW1oaLhgwYL+hNF6f+7cuZOUlMRgMIyMjHCa+SktLb1y5UpgYGB7e7tEL37nzp2zZ8/u2LFDQQFFcL0yeCeys7NfvHiBIEhOTs769etZLNa2bdsKCws7Ozu3bdt27Nix+/fv46d+7NixU6dOXblyZdasWTdv3sRPCEX8DktOTs6sWbPy8/PRj0Kh8O7du3hXQJySkpI1a9YUFBTgdPy9e/fW1tYiCPLs2bMDBw6sX7++vLxctBX90hXDB9S7x8XF+fr6ohcrubm5cXFxISEhaKg9BEHu3bv3+eefy10UQRDRLfrjx4+rqqqifXx4eLjoyQ2509HRsWTJEmNj4ylTpnh7e5NIJLSPDw8PV0yP3pf6+nr8nv9JTU2tqKigUChtbW1Tp05VUVGJi4s7cuQITnIy+IAuVf/5558hQ4agjwRZWlrq6+tfunTJ3d2dTCYTCAQ8nu999OhRSEhIRUVFR0cHlUqdMGHCrVu3OBzO8OHDExMTPTw8+hNpdQBg3qwlk8m//fYbfqKywXXFt729vaqqqr29/axZs9CL1NzcXF9fX3k+2ds/Bt/uFRUVSUlJLi4ulZWVKioqorGdpaUliUTCzwENDQ29vb3oGTczM8vJyfnzzz9VVVVv376tp6c3bdo0Op0ud1GhUIgqStysTUtLGz58uJeXlwKenh8UGAyGvr5+dXW1urr68ePHv/jiC5xu1srmgxjMnDp1isViWVlZXb58mUgkmpqaDhs2bNiwYXZ2dkwmU1oM+PchKysrKyuLQCC4urqisRcXLlwIAKirq7t37x6CIHg8hJOdnf37779TKJS1a9daWVlFRUVRqVQejxccHNzY2FhWVobHfdMPh+7u7vj4+IaGhm+//Xb8+PGDUwmFXSXIJiEh4ccffywqKhIKhc+ePcvKyrp06RJOWhkZGTExMUKhUCgUrlq1Cnl9nYqTnEg0Ojqax+OVl5eHhIQgCrlZC5HgQ1nesWjRoq+++urvv/8GAJibm/v6+n7zzTd4CKFT3QEBAQQCoaenh8Vibdy4ceTIkTweLzExEQ9F8HpSf+nSpeKT+sbGxs7OzsXFxTiJQvoy+GN3EZ9//nllZWVeXh5OT5yiaGho6OrqJiQkuLi4HDhwwN/ff8aMGefOnVNXV7e3t8fp3p6Wlpampubp06ddXV1jY2PnzJmDXnnjd7MWgs1g/7xIgt+TSeIobKq7qKgoLCxs3759PT09CpvUh0jjg3uIQDH3MidMmAAASE9PHz16tJqaGoFAwGNSv7S09Pr162FhYe3t7WQyGV0OkpGR4ezsjN+kPkQGH8rYXTFkZGRs3Ljx7NmzAIAJEyb4+Phs3boVpxXHAIDs7Oy5c+dqamoiCBIVFbVhwwYzMzMPD4+tW7f29vbiJAqRwQcxEakYnjx5cufOHX9//+TkZBaLFRAQAAB49uwZHhOdKNLuJj558gTvJ10hmHxwgxmcQBAkMzMTDYiFOj4xMXHhwoX4eR0AMG3atNu3b1OpVPQh3t7eXhqNhiAI9Ppg8Sn37g8ePKitrZ04cWJqampOTo6bm1tTU5OtrS36IDuu/bo4QqHw8ePHZmZm8fHxNjY23377rQJEIZh8yr27k5NTTk6OUChkMplRUVHJycna2troO5UmTpyoGK8DALhc7qVLlwb5biIEAPBp9+4AAARBYmNj2Wz2+vXrAQDJycldXV2fffaZnIP1QD4SPnG7g9eOF41hysrKHBwcBrtSkMHh07c76ON4iNKiFHYHACAIcuvWLTiGUXKUxe4QCFC2u6oQJQfaHaJEQLtDlAhod4gSAe0OUSKg3SFKBLQ7RIn4f+qNoOvABy4fAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 250 -h 200 -u px\n",
    "fviz_pca_contrib(pca, choice=\"var\", axes=1:2, top=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is of course different from getting a mean of relative contributions (which would assume equal importance of each of the principal components):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plasmin    0.315266\n",
       "ANP        0.311123\n",
       "STAB2      0.303710\n",
       "FGF-8B     0.302288\n",
       "PECAM-1    0.293163\n",
       "dtype: float64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relative_contributions[['PC1', 'PC2']].mean(axis=1).nlargest()"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}