{"cells":[{"metadata":{"_uuid":"822ff03139a0e00af3b834e6ad8d142f8a10c5bd"},"cell_type":"markdown","source":"# Which proteins come together?\n\nTake a look at this wonderful image of an animal cell provided by LadyofHats (Mariana Ruiz) in the public domain via Wikimedia Commons. We can see some of our target proteins and may conclude that some of them are likely to occur together. In our images all of them should be present but only one of them are stained in the green channel. Now let's assume that staining is sometimes not that easy and it happens that multible organelles are likely to be stained together. If this is true we could expect groupings - some targets may be likely to be one-hot at the same time over a broader range of image samples in our data set.\n\nAs the target distribution and dependencies are an entry point to setup an objective or loss function, it could be worth it to dive into dive with me into target group analysis using a latent variable model. \n\n**If you like my kernel** you can make be very happy with an **upvote and/or comment** ;-)! The motivation I gain out of your feedback pushes me to share my ideas instead of hiding them. Thank you!"},{"metadata":{"_uuid":"4cdf45df629f9e62d0ec750048e3da692ea20f32"},"cell_type":"markdown","source":"![Animal cell organelles](https://upload.wikimedia.org/wikipedia/commons/4/48/Animal_cell_structure_en.svg)"},{"metadata":{"_uuid":"316cce53e8ca13c185774a72f29ac65435be3f0c"},"cell_type":"markdown","source":"## What can you find within this kernel?\n\n1. A motivation why there is a hint that target groups exist in our data. \n2. A short explanation what is a bernoulli mixture model, how it learns and its implementation.\n3. An overview of target groups found by clustering with this model. \n4. Some explorations related to certainty of cluster assignment and target combination anomalies.\n5. A conclusion \n\nLet's go! :-)"},{"metadata":{"_uuid":"06923aa2089236f71c3eed1037d8f9676fb0fa7d"},"cell_type":"markdown","source":"# Preliminary Work\n\nAs usual - loading packages:"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport seaborn as sns\nsns.set()\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom PIL import Image\nfrom scipy.misc import imread\n\nimport tensorflow as tf\n\n\nimport os\nprint(os.listdir(\"../input\"))\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=RuntimeWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e2d56e08bfb1940cc0dcbb7c27977acf1c0b2dfa"},"cell_type":"markdown","source":"and target information data given by train.csv:"},{"metadata":{"trusted":true,"_uuid":"1ac02466bf8f2d58f67f5d1067093da45ef49f12"},"cell_type":"code","source":"train_labels = pd.read_csv(\"../input/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dd8b708a22e7ac156b24449f0544b4c8bf7a239c"},"cell_type":"markdown","source":"After extracting the labels per id with zero-hot-encoding:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"b545a6b119ebcd1f4b326b521fa785b3e58bc468"},"cell_type":"code","source":"label_names = {\n    0:  \"Nucleoplasmn\",  \n    1:  \"Nuclear membrane\",   \n    2:  \"Nucleoli\",   \n    3:  \"Nucleoli fibrillar center\",   \n    4:  \"Nuclear speckles\",\n    5:  \"Nuclear bodies\",   \n    6:  \"Endoplasmic reticulum\",   \n    7:  \"Golgi apparatus\",   \n    8:  \"Peroxisomes\",   \n    9:  \"Endosomes\",   \n    10:  \"Lysosomes\", \n    11:  \"Intermediate filaments\",   \n    12:  \"Actin filaments\",   \n    13:  \"Focal adhesion sites\",   \n    14:  \"Microtubules\",   \n    15:  \"Microtubule ends\",   \n    16:  \"Cytokinetic bridge\",   \n    17:  \"Mitotic spindle\",   \n    18:  \"Microtubule organizing center\",   \n    19:  \"Centrosome\",   \n    20:  \"Lipid droplets\",   \n    21:  \"Plasma membrane\",   \n    22:  \"Cell junctions\",   \n    23:  \"Mitochondria\",   \n    24:  \"Aggresome\",   \n    25:  \"Cytosol\",   \n    26:  \"Cytoplasmic bodies\",   \n    27:  \"Rods & rings\"\n}\n\nreverse_train_labels = dict((v,k) for k,v in label_names.items())\n\ndef fill_targets(row):\n    row.Target = np.array(row.Target.split(\" \")).astype(np.int)\n    for num in row.Target:\n        name = label_names[int(num)]\n        row.loc[name] = 1\n    return row","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0c367f67bf61b62b57bef9712968fe3064fe63ea","_kg_hide-input":true},"cell_type":"code","source":"for key in label_names.keys():\n    train_labels[label_names[key]] = 0\n\ntrain_labels = train_labels.apply(fill_targets, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c05eecb05fd15fb162767165e31c5726595b7445"},"cell_type":"code","source":"train_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"812141fa937c6b5f82e1ee94ef5f4776129a408d"},"cell_type":"markdown","source":"We are ready to start! :-D"},{"metadata":{"_uuid":"92329908fdfd918b7780bb0afd53226ecc588d86"},"cell_type":"markdown","source":"# Why do proteins come together?"},{"metadata":{"trusted":true,"_uuid":"5b061b629d0d919218e6114ab2cdc7eec1d2cc41","_kg_hide-input":true},"cell_type":"code","source":"target_counts = train_labels.drop([\"Id\", \"Target\"],axis=1).sum(axis=0).sort_values(ascending=False)\nplt.figure(figsize=(15,15))\nsns.barplot(y=target_counts.index.values, x=target_counts.values, order=target_counts.index);","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false,"_uuid":"19eb4d692bdb7ab0bd74f217cfb1cf9b7fd75a7e"},"cell_type":"markdown","source":"This is already known! We have some very seldom targets like rods & rings and very common ones like nucleoplasmn, cytosol and plasma membrane. But wouldn't it be nice to know if some proteins are likely to come together? We have already seen that lysosomes, endosomes and endoplasmatic reticulum have target correlations:\n\n"},{"metadata":{"trusted":true,"_uuid":"dcc8f3acd60f4d07bb620dd8b2f66150c48e77db"},"cell_type":"code","source":"train_labels[\"number_of_targets\"] = train_labels.drop([\"Id\", \"Target\"],axis=1).sum(axis=1)\n\ndef find_counts(special_target, labels):\n    counts = labels[labels[special_target] == 1].drop(\n        [\"Id\", \"Target\", \"number_of_targets\"],axis=1\n    ).sum(axis=0)\n    counts = counts[counts > 0]\n    counts = counts.sort_values()\n    return counts\n\nlyso_endo_counts = find_counts(\"Lysosomes\", train_labels)\n\nplt.figure(figsize=(15,5))\nsns.barplot(x=lyso_endo_counts.index.values, y=lyso_endo_counts.values, palette=\"Blues\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5ee2618bf024a9e31069ff447416f223ceb45931"},"cell_type":"markdown","source":"**You can see that lysosomes and endosomes always come together and that it is likely that the endoplasmatic reticulum is present as well**. I have read about lysosomes that their membrane is produced by ribosomes that are located in the rough endoplasmatic reticulum and transferred to the golgi apparatus afterwards. Later they are shipped to endosomes, small bubbles full of stuff, and seem to fuse with them to digest all the stuff inside. What if staining is done with molecules that are found in all three participants: the rough ER, endosomes and lysosomes? Then we will always see that they often come together. \n\nBut isn't it nice to know? If our model is sure that lysosomes are present, we can automatically say endosomes is hot as well. Great! Hence instead of predicting both target classes we can reduce to both to one single lyso-endo class."},{"metadata":{"trusted":true,"_uuid":"45b7c2d87fd3fa439f46b7469a10e12252f7c6d4","_kg_hide-input":true},"cell_type":"code","source":"count_perc = np.round(100 * train_labels[\"number_of_targets\"].value_counts() / train_labels.shape[0], 2)\nplt.figure(figsize=(20,5))\nsns.barplot(x=count_perc.index.values, y=count_perc.values, palette=\"Reds\")\nplt.xlabel(\"Number of targets per image\")\nplt.ylabel(\"% of data\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a5591d0c882150bae5696a59c9cf74a3724e56e4"},"cell_type":"markdown","source":"Here you can see that most of the target proteins per image come alone or as a pair. Hopefully we can find nice pair-structures in our data. We will see... "},{"metadata":{"trusted":true,"_uuid":"894c3d1273925dc0d44dd9f3bf91cc70e6b43ec2"},"cell_type":"code","source":"targets = train_labels.drop([\"Id\", \"Target\", \"number_of_targets\"], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"49adbc20ebd6a0c2e8530b1d18e40bfe47397d83"},"cell_type":"markdown","source":"# How can we uncover hidden protein groups?\n\nThe presence of target proteins is given by a binary, discrete value:\n\n* 0 for absence\n* 1 for presence \n\nAs I like to find groupings or clusters of these **discrete quantities** I liked to use a discrete clustering algorithm. As I'm currently on the way of my learning path to discover mixture models the choice using a Bernoulli mixture model was an easy one.  What's the idea behind it?\n\n*If you are not interested in mixture models and math etc. skip this chapter and jump to results analysis ;-) *\n\nMixture models assume that the data we **observed $X$ was generated by some latent variable $Z$.** In our case this could be the **staining process** of microscope preparates used to obtain the 4 channel images. Perhaps this process is limited and this limitation has caused our target groupings. You have already seen that we were able to reduce the number of target classes by fusing endosomes and lysosomes. Consequently we can say that there are fewer configurations of the latent variable than target classes. \n\n\n### Model description\n\nThese configurations can be seen as $K$ components of the mixture model. We don't know how many of them are actually there and we will have to estimate them during the analysis. Each component tries to explain one target group we are seeking for. And for each sample $x_{n}$ of our $N$ data spots there exists a related latent or hidden variable $z_{n}$ that holds 1 for the component $k$ that generated $x_{n}$ and 0 for all others. Imagine you would already know them, then we could describe the probability density our data as follows:\n\n$$ P(X) = \\sum_{Z} P(X, Z|\\theta |\\theta) = \\sum_{Z} P(Z|\\theta) \\cdot P(X|Z, \\theta)$$\n\nIf you are not very familiar with probabilities: To obtain this equation you need to know the sum and product rules of probabilites. We are now summing over all $K$ configurations of the latent variable Z, hence over $K$ assumed clusters. This formulation is very general. Now image you take a look at one of the components, on one single cluster. How is the target protein data distributed in that group?\n\n<a title=\"Classical Numismatic Group, Inc. http://www.cngcoins.com [GFDL (http://www.gnu.org/copyleft/fdl.html), CC-BY-SA-3.0 (http://creativecommons.org/licenses/by-sa/3.0/) or CC BY-SA 2.5 (https://creativecommons.org/licenses/by-sa/2.5)]\" href=\"https://commons.wikimedia.org/wiki/File:Ephesos_620-600_BC.jpg\"><img width=\"256\" alt=\"Ephesos 620-600 BC\" src=\"https://upload.wikimedia.org/wikipedia/commons/4/4f/Ephesos_620-600_BC.jpg\"></a>\n\nEach target protein $x_{d}$ itself is binary. If we like to describe the distribution of nucleoplasmn it is like one we would obtain by tossing a coin. And for each of the $D=28$ target protein we could flip such a coin, with zero on one side and one on the other. If we assume that  the target proteins are independent within one component $k$ than we can write:\n\n$$ p(x_{n}|\\mu_{k}) = \\prod_{d=1}^{D} \\mu_{k,d}^{x_{n,d}} (1-\\mu_{k,d})^{(1-x_{n,d})}$$\n\nWhereas $\\mu_{k,d}$ stands for the probability to observe the current target protein $x_{n,d}$. This sounds difficult first but let's try to understand it with lysosomes and endosomes again: If we know that they belong to group component $k=2$, for example. Then the probability to observe lysosomes and endosomes within this group could perhaps be $\\mu_{2,lyso} = \\mu_{2,endo} =  0.99$ for both and $\\mu_{2,ER} = 0.4$ for endoplasmatic reticulum. For all other proteins the probability to observe them should be very low, for example $\\mu_{2, k} = 0.01.$ If we now have a sample that fits to this target combination it would yield a high value for $p(x_{n}| \\mu_{2})$. Hence it's probable that it belongs to that group.\n\nOk, we are close to fullfill our model decription. We try to explain our data by the probability density $P(X)$. Imagine that all samples were independently drawn from this distribution. With this assumption we can split into a product over all $N$ samples. \n\n$$ P(X) = \\sum_{Z} P(X, Z|\\theta) = \\sum_{Z} P(Z|\\theta) \\cdot P(X|Z, \\theta) = \\prod_{n} \\sum_{z} p(z|\\pi) \\cdot p(x|z, \\mu)$$\n\nAs there exists one true component $z_{n,k}$ for each sample, we can say that $z_{n}$ is one-hot-encoded and can be described by a multinomial distribution:\n\n$$p(z|\\pi) = \\prod_{k=1}^{K} \\pi_{k}^{z_{k}}$$\n\nThe same holds for the conditional probability $p(x|z, \\mu)$:\n\n$$p(x|z, \\mu) = \\prod_{k=1}^{K} p(x|\\mu_{k})^{z_{k}}$$\n\nand finally we obtain: \n\n$$ P(X) = \\prod_{n} \\sum_{z} p(z|\\pi) \\cdot p(x|z, \\mu)  = \\prod_{n} \\sum_{k}^{K} \\pi_{k} \\cdot p(x_{n}|\\mu_{k}) = \\prod_{n} \\sum_{k=1}^{K} \\pi_{k} \\cdot \\prod_{d=1}^{D} \\mu_{k,d}^{x_{n,d}} (1-\\mu_{k,d})^{(1-x_{n,d})} $$\n\nWith this model we like to describe our target data. And to make this model fit to what we observe we will maximize this probability density $p(X)$ with respect to our model parameters $\\pi$ and $\\mu$. This maximization is usually performed by taking the log first as it often makes it simpler to take derivarives. But in our case we will stuck....\n\n$$ \\ln p(X|\\mu, \\pi) = \\sum_{n=1}^{N} \\ln \\left( \\sum_{k=1}^{K} \\pi_{k} \\cdot p(x_{n}|\\mu_{k}) \\right) $$ \n\nDo you see the problem? The sum over all $K$ components prevents the log to act on $\\pi_{k} \\cdot p(x_{n}|\\mu_{k})$ and consequently we can't make the bernoullis tractable for taking derivatives. :-(  "},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"### Expectation maximization\n\n$$ \\ln p(X,Z|\\mu, \\pi) = \\sum_{n=1}^{N} \\sum_{k=1}^{K} \\gamma (z_{n,k}) \\left( \\ln \\pi_{k} + \\sum_{d=1}^{D} x_{n,d} \\ln \\mu_{k,d} + (1- x_{n,d}) \\ln (1-\\mu_{k,d}) \\right)$$\n\nNow the hero comes: Instead of maximizing the log-likelihood we are going to act as if we already know our parameters $\\pi_{k}$ and $\\mu_{k}$ by initializing them randomly.  \n\n#### E-Step: responsibilites\n\nThen we compute how responsible each component of our mixture $k$ was to generate the data spot $x_{n}$ and its targets. This is like a soft assigment to cluster components. The one cluster with the highest responsibility would be the winner for cluster assignment. \n\n$$ \\gamma( z_{n,k} ) = \\frac{\\pi_{k} p(x_{n}|\\mu_{k})} {\\sum_{j=1}^{K} \\pi_{j} p(x_{n}|\\mu_{j})} $$\n\n#### M-Step: maximization\n\nNow we know how responsible each component is for generating data point $x_{n}$ and given this information we can recalculate the parameters we initialized randomly at the starting point.\n\n$$ N_{k} = \\sum_{n=1}^{N} \\gamma(z_{n,k})$$\n\n$$\\mu_{k} = \\frac{1}{N_{k}} \\sum_{n=1}^{N} \\gamma(z_{n,k}) x_{n}$$\n\n$$\\pi_{k} = \\frac{N_{k}}{N} $$\n\nYou may wonder where do these nice equations come from?! To understand it, you need to know more about expectation maximization and I don't like to blow up this kernel even more with math. So if you like, you should dive deeper into this topic by reading books or watching videos in the orbit out there. ;-)\n\nPerforming E- and M-Step iteratively one step after another we will end up with a nice fit of our model to the target protein data.\n"},{"metadata":{"_uuid":"046d4d1cd26caa98e78ba67c6455c57431dd5d37"},"cell_type":"markdown","source":"### Implementation of the model"},{"metadata":{"trusted":true,"_uuid":"6af87f8e5b43b25410c6b4d11c7fb5cf0f3cde33"},"cell_type":"code","source":"from scipy.special import logsumexp\n\nclass BernoulliMixture:\n    \n    def __init__(self, n_components, max_iter, tol=1e-3):\n        self.n_components = n_components\n        self.max_iter = max_iter\n        self.tol = tol\n    \n    def fit(self,x):\n        self.x = x\n        self.init_params()\n        log_bernoullis = self.get_log_bernoullis(self.x)\n        self.old_logL = self.get_log_likelihood(log_bernoullis)\n        for step in range(self.max_iter):\n            if step > 0:\n                self.old_logL = self.logL\n            # E-Step\n            self.gamma = self.get_responsibilities(log_bernoullis)\n            self.remember_params()\n            # M-Step\n            self.get_Neff()\n            self.get_mu()\n            self.get_pi()\n            # Compute new log_likelihood:\n            log_bernoullis = self.get_log_bernoullis(self.x)\n            self.logL = self.get_log_likelihood(log_bernoullis)\n            if np.isnan(self.logL):\n                self.reset_params()\n                print(self.logL)\n                break\n\n    def reset_params(self):\n        self.mu = self.old_mu.copy()\n        self.pi = self.old_pi.copy()\n        self.gamma = self.old_gamma.copy()\n        self.get_Neff()\n        log_bernoullis = self.get_log_bernoullis(self.x)\n        self.logL = self.get_log_likelihood(log_bernoullis)\n        \n    def remember_params(self):\n        self.old_mu = self.mu.copy()\n        self.old_pi = self.pi.copy()\n        self.old_gamma = self.gamma.copy()\n    \n    def init_params(self):\n        self.n_samples = self.x.shape[0]\n        self.n_features = self.x.shape[1]\n        #self.gamma = np.zeros(shape=(self.n_samples, self.n_components))\n        self.pi = 1/self.n_components * np.ones(self.n_components)\n        self.mu = np.random.RandomState(seed=0).uniform(low=0.25, high=0.75, size=(self.n_components, self.n_features))\n        self.normalize_mu()\n    \n    def normalize_mu(self):\n        sum_over_features = np.sum(self.mu, axis=1)\n        for k in range(self.n_components):\n            self.mu[k,:] /= sum_over_features[k]\n            \n    def get_responsibilities(self, log_bernoullis):\n        gamma = np.zeros(shape=(log_bernoullis.shape[0], self.n_components))\n        Z =  logsumexp(np.log(self.pi[None,:]) + log_bernoullis, axis=1)\n        for k in range(self.n_components):\n            gamma[:, k] = np.exp(np.log(self.pi[k]) + log_bernoullis[:,k] - Z)\n        return gamma\n        \n    def get_log_bernoullis(self, x):\n        log_bernoullis = self.get_save_single(x, self.mu)\n        log_bernoullis += self.get_save_single(1-x, 1-self.mu)\n        return log_bernoullis\n    \n    def get_save_single(self, x, mu):\n        mu_place = np.where(np.max(mu, axis=0) <= 1e-15, 1e-15, mu)\n        return np.tensordot(x, np.log(mu_place), (1,1))\n        \n    def get_Neff(self):\n        self.Neff = np.sum(self.gamma, axis=0)\n    \n    def get_mu(self):\n        self.mu = np.einsum('ik,id -> kd', self.gamma, self.x) / self.Neff[:,None] \n        \n    def get_pi(self):\n        self.pi = self.Neff / self.n_samples\n    \n    def predict(self, x):\n        log_bernoullis = self.get_log_bernoullis(x)\n        gamma = self.get_responsibilities(log_bernoullis)\n        return np.argmax(gamma, axis=1)\n        \n    def get_sample_log_likelihood(self, log_bernoullis):\n        return logsumexp(np.log(self.pi[None,:]) + log_bernoullis, axis=1)\n    \n    def get_log_likelihood(self, log_bernoullis):\n        return np.mean(self.get_sample_log_likelihood(log_bernoullis))\n        \n    def score(self, x):\n        log_bernoullis = self.get_log_bernoullis(x)\n        return self.get_log_likelihood(log_bernoullis)\n    \n    def score_samples(self, x):\n        log_bernoullis = self.get_log_bernoullis(x)\n        return self.get_sample_log_likelihood(log_bernoullis)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"56e4c1cd3dcda777c1a5da3c62dcb4a1ef940cd5"},"cell_type":"markdown","source":"### Estimating the number of components\n\nI don't know how many target groups are there in advance. But in contrast to hard clustering algorithms like k-means we can use a test set to tune the number of components to choose as a hyperparameter. "},{"metadata":{"trusted":true,"_uuid":"d1bfa7d6d31324655750fa1383aab5d975c4b627"},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX = targets.values\nx_train, x_test = train_test_split(X, shuffle=True, random_state=0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3766b15af4dfc74a03cc57a1849f514165a536f5"},"cell_type":"markdown","source":"Let's try out some values:"},{"metadata":{"trusted":true,"_uuid":"283a371abca76f8b7e8c63fe856239769d52eca3"},"cell_type":"code","source":"components_to_test = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc4a324536b7fc94e7ca5e908a3d103a4682d38a"},"cell_type":"markdown","source":"And fit multiple models... "},{"metadata":{"trusted":true,"_uuid":"b23a6f0a6e6f9afb301ebc3e08c990ca05699ac9","_kg_hide-input":true},"cell_type":"code","source":"scores = []\n\n\nfor n in range(len(components_to_test)):\n    if n > 0:\n        old_score = score\n    model = BernoulliMixture(components_to_test[n], 200)\n    model.fit(x_train)\n    score = model.score(x_test)\n    scores.append(score)\n    if n > 0: \n        if score < old_score:\n            estimated_components = components_to_test[n-1]\n            break\n        ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"03633974426986dbdcb9463abbafac8b5f49a0c3"},"cell_type":"markdown","source":"In the end, we obtain that these number of components was nice for train and test:"},{"metadata":{"trusted":true,"_uuid":"852a2fd0ee9031e3b861fc1bdce5ae35bfff19e1","_kg_hide-input":true},"cell_type":"code","source":"estimated_components","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9d42c95b213ca5c5992f66c4a736bfcdc90089ec"},"cell_type":"markdown","source":"To obtain results to work with, let's refit out model with this estimated number of components:"},{"metadata":{"trusted":true,"_uuid":"fb21a2ce3995a778afa0fd0500f1c494ff9ab31e"},"cell_type":"code","source":"model = BernoulliMixture(estimated_components, 200)\nmodel.fit(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7573962d1b4a63de6868e72b3a9807307f59c9ad"},"cell_type":"code","source":"results = targets.copy()\nresults[\"cluster\"] = np.argmax(model.gamma, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a429ae73a2ac754689c18bb00b0f8c6c0e504831"},"cell_type":"markdown","source":"## What kind of target groups are found?"},{"metadata":{"_uuid":"4cdbeebe4e5d941031efb7d5e0ba32eb274154d6"},"cell_type":"markdown","source":"### How are specific proteins distributed over all clusters in percent?"},{"metadata":{"trusted":true,"_uuid":"3791c585f38446229dc5bf41faf8d5a1ba66d061","_kg_hide-input":true},"cell_type":"code","source":"grouped_targets = results.groupby(\"cluster\").sum() / results.drop(\"cluster\", axis=1).sum(axis=0) * 100\ngrouped_targets = grouped_targets.apply(np.round).astype(np.int32)\n\nplt.figure(figsize=(20,15))\nsns.heatmap(grouped_targets, cmap=\"Blues\", annot=True, fmt=\"g\", cbar=False);\nplt.title(\"How are specific proteins distributed over clusters in percent?\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9e8bbba84c05e8990443dc69c2c59b28311309b6"},"cell_type":"markdown","source":"### Take-Away\n\n* This looks great! You can see that several clusters only hold one specific target protein!\n* For each target protein you can see the percentage of its occurences that are placed into specific clusters.\n* One example: 97 % of Mitochondria target proteins are located in cluster 18. Only a few percents are hold by cluster 12 and 8. There is one percent missing to fill up to 100 % but this is caused by rounding errors and should not worry you.\n* We find that **a lot of cellular components are roughly hold by their own clusters**: \n    * Nuclear speckles\n    * Endoplasmatic reticulum, Lysosomes and Endosomes\n    * Golgi apparatus\n    * Intermediate filaments\n    * Microtubules\n    * Lipid Droplets\n    * Centrosomes\n    * Plasma membrane\n    * Mitochondria\n* In addition we find **meaningful combinations of cellular components**:\n   * Focal adhesion sides come together with actin fliaments and are spread over 3 clusters that are very similar up to absence or presence of plasma membrane, nuclear membrane and cell junctions. \n    * Aggresomes come alone or sometimes with cytoplasmic bodies. This makes sense: If bodies consist of viral capsids then there could be a lot of garbage or clutted proteins as well either due to immune response or to the viral building process. \n    * Cellular devision often comes in groups of either microtubules & mitotic spindle and cytokinetic bridge or as microtubule organizing center & mitotic spindle and cytokinetic bridge.\n    * As biologist one may see much more!\n* **Some components are spread over several clusters** like nucleoplasmn.\n "},{"metadata":{"_uuid":"007c8bedb0e2c6d76ac9439c19eb0b1262b0afa7"},"cell_type":"markdown","source":"### Give clusters a name\n\nNaming clusters will make it easier to understand the patterns. I like to do this given the information of the blue map that show us in which cluster one can find a specific target protein. This information does not suffer under the target imbalance problem and show nice couplings between targets.  "},{"metadata":{"trusted":true,"_uuid":"e499bfffaf3e26078bf2956c965dffe2bfc1ce02"},"cell_type":"code","source":"cluster_names = {\n    0: \"Actin filaments & Focal adhesion sites\",\n    1: \"Aggresomes\",\n    2: \"Microtubules, Mitotic spindle, Cytokinetic Bridge\",\n    3: \"RodsRings, Microtubule ends, Nuclear bodies\",\n    4: \"Some nuclear membranes\",\n    5: \"various - RodsRings\",\n    6: \"various - Mitotic spindle, Organizing center, Cytosol\",\n    7: \"Nuclear bodies & Aggresomes\",\n    8: \"Nuclear speckles\",\n    9: \"Nuclear membrane & Actin filaments & Focal adhesion sites\",\n    10: \"Endoplasmatic reticulum & Endosomes & Lysosomes\",\n    11: \"Low dense 1\",\n    12: \"Mitotic spindle, Organizing center\",\n    13: \"Plasma membrane & various\",\n    14: \"Nucleoli fibrillar center & Peroxisomes\",\n    15: \"Low dense 2\",\n    16: \"Nucleoli fibrillar center & Cytoplasmic bodies\",\n    17: \"Nucleoli & Microtubule ends & Peroxisomes & Rods Rings\",\n    18: \"Mitochondria & Lipid droplets & RodsRings & Nucleoli\",\n    19: \"Low dense 3\",\n    20: \"Golgi apparatus\",\n    21: \"Intermediate filaments\",\n    22: \"Centrosome\",\n    23: \"Cytoplasmic bodies & Aggresomes\",\n    24: \"Lipid droplets & Peroxisomes & Cell junctions\"\n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d85790713e29acb4a18627ce428f2823fde03cd2"},"cell_type":"markdown","source":"### How many targets are hot within one cluster?\n\nLet's go one step further. We already know that cluster 10 holds all samples of endoplasmatic reticulum, lysosomes and endosomes. But we don't know how present each target is given only the targets of one cluster. "},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7385d53881435696c50ff1ed612c8b15cc7fa2ee"},"cell_type":"code","source":"cluster_size = results.groupby(\"cluster\").Nucleoplasmn.count()\ncluster_composition = results.groupby(\"cluster\").sum().apply(lambda l: l/cluster_size, axis=0) * 100\ncluster_composition = cluster_composition.apply(np.round).astype(np.int)\n\ncluster_composition = cluster_composition.reset_index()\ncluster_composition.cluster = cluster_composition.cluster.apply(lambda l: cluster_names[l])\ncluster_composition = cluster_composition.set_index(\"cluster\")\n\nplt.figure(figsize=(20,20))\nsns.heatmap(cluster_composition, cmap=\"Oranges\", annot=True, fmt=\"g\", cbar=False);\nplt.title(\"How present alias hot are specific targets within one cluster?\");\nplt.ylabel(\"\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4c6d631c1ab83c8d1f78f735399d899f9c6619f6"},"cell_type":"markdown","source":"### Take-Away\n\n* That's great! This may yields insights what kind of targets we will always find per cluster. Take a look at the endoplasmatic reticulum, lysosome, endosome cluster. You can see that almost all samples of that cluster hold the endoplasmatic reticulum with 1. In contrast due to the seldomness of lysosomes and endosomes the hotness of them is very low. A lot of samples do not have them even though all targets of lysosomes and endosomes can be find within that cluster!!!\n* Well you might ask yourself now, why this map helps you: Take a look at the aggresomes for example. It's a seldom target but once you have found the cluster 1 you can definitely say that all of its samples have an aggresome present. :-) Hence if seldom targets have their own nice, little cluster and you have found it, you are done with its targets. "},{"metadata":{"_uuid":"2c9a856f82dae623fb90b6397d2f177cf77c419b"},"cell_type":"markdown","source":"### Take away\n\n* This map seems to be more difficult to understand than the one before.\n* It's highly influenced by the frequency of target proteins and **yields more insights of the imbalance of target proteins per cluster**. Let's try to collect important insights! :-)\n* First of all it comes out very clear that **nucleoplasmn and cytosol are the dominating targets** that influence the composition of many clusters. \n* Now let's come to our example of endoplasmatic reticulum, lysosomes and endosomes: You can see that the cluster that hold these target proteins is mainly occupied by endoplasmatic reticulum followed by cytosol. Lysosomes and endosomes in contrast only appear with very small percentages of 3 and 2 % due to their seldomness. **Hence even though they are all located in this cluster, they are seldom in this cluster as well and not only in the overall data!**\n* This situation is even more worse for Rods and Rings and Microtubule ends. Due to their seldomness they are overwhelmed by more common classes like nucleoplasm, nucleoli and nuclear bodies. \n* Luckily we can see some **nice clusters that are more specific to seldom targets** like that with lipid droplets, cell junctions and peroxisomes. \n* Same holds for aggresomes and aggresomes with cytoplasmic bodies. "},{"metadata":{"trusted":true,"_uuid":"3b1534abfd383f4870257aa1457b1e26ea8a2235"},"cell_type":"code","source":"results[\"cluster_names\"] = results.cluster.apply(lambda l: cluster_names[l])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e7eecc7417024654645da1183a6e936485a79ea4"},"cell_type":"markdown","source":"Both maps already yielded some insights which targets are likely to come together and how seldom a target is given a cluster. We can go one step further. The model was trained with $\\mu$ and $\\pi$. The latter shows us the probability per cluster: "},{"metadata":{"_uuid":"d84603bef889884af8d842b911ba05a7b1b1e3bb"},"cell_type":"markdown","source":"### How does the prior probability $\\pi_{k}$ per cluster look like?"},{"metadata":{"trusted":true,"_uuid":"0d2c5bd4045a4d590f631339d517093dad90036a"},"cell_type":"code","source":"cluster_ids = np.arange(0, estimated_components)\nnames = [cluster_names[l] for l in cluster_ids]\n\npi = pd.Series(data=model.pi, index=names).sort_values(ascending=False)\nplt.figure(figsize=(20,5))\nsns.barplot(x=pi.index, y=pi.values, palette=\"Reds_r\", order=pi.index)\nplt.xticks(rotation=90);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9f7bd203cdf14c46ac0de1494b510422f7f10a52"},"cell_type":"markdown","source":"### Take-away\n\nUhh, that's interesting! I haven't thought that the cluster with mitochondria has the highest prior probability. I expected cluster 5 and 6 to be most probable:\n\n* 5: \"various - RodsRings\"\n* 6: \"various - Mitotic spindle, Organizing center, Cytosol\"\n\nBoth have high counts of nucleoplasmn and cytosol, hence most common targets. Even though Mitochondria with Nucleoli are at the top we can see that clusters that hold most common targets have higher prior probabilities."},{"metadata":{"_uuid":"8d5d3454d38f27cb218967ebe99abf85819ff9de"},"cell_type":"markdown","source":"## How do the target probabilities $\\mu_{kd}$ look like given the cluster?"},{"metadata":{"trusted":true,"_uuid":"19ea4bd8f935867f9e2c70f2d49fe1b3a1f7814e"},"cell_type":"code","source":"model.mu.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f4673756c8846a7b5fd30c3f6770b7e7f8aba2e0","_kg_hide-input":true},"cell_type":"code","source":"mu = pd.DataFrame(data=model.mu * 100, index=names, columns=results.drop([\"cluster\", \"cluster_names\"], axis=1).columns.values)\nmu = mu.apply(np.round)\n\nplt.figure(figsize=(20,20))\nsns.heatmap(mu, cmap=\"Purples\", annot=True, fmt=\"g\", cbar=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5df75b5653b6ab8ca3c4979088e8cdb3b875c422"},"cell_type":"markdown","source":"### Take-Away\n\n* We can see that the probability of each target to be present yields a lot of insights. Let's consider the golgi apparatus for example. If we choose the corresponding cluster named the same we can observe a probability to occur of 100%. Hence if we now that a sample is located within that cluster, we can say it has a golgi apparatus!\n* There a lot of target proteins that have a very high probability to occur in their cluster!\n* But there are some strange probabilities as well. Think of cluster 1 again: It's nearly fully occupied by aggresomes but the mu-proability for this target given this cluster is very low with 18 %. That's strange! It's related to the seldomness of the target and could be caused by a shift due to weighted sum during mu calculation. As the sum is taken over all samples, nearby targets (of other components) that still have some higher responsibility for the aggresome cluster can cause a shift of the mu-center of the cluster. "},{"metadata":{"_uuid":"7a5f1dea325ec0ce5d24b9d59fd63519e033f544"},"cell_type":"markdown","source":"## How many samples do the clusters hold?"},{"metadata":{"trusted":true,"_uuid":"59c28bbed372b3b24ff8ac7f50d863da8323ed52","_kg_hide-input":true},"cell_type":"code","source":"cluster_counts = results.groupby(\"cluster\").cluster.count()\ncluster_counts = cluster_counts.sort_values()\nnames = [cluster_names[num] for num in cluster_counts.index]\n\nplt.figure(figsize=(20,5))\nsns.barplot(x=names, y=cluster_counts.values, order=names)\nplt.xticks(rotation=90);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b88ca92f080f9e20bbbc6e79fb7d594b4cbfd4ac"},"cell_type":"markdown","source":"### Take-Away\n\n* We can see that the clusters are occupied by very different amounts of samples.\n* This fits to what we have observed: Small clusters are mainly responsible for seldom targets and do from nice little groups.\n* In contrast big clusters hold various of different target proteins and perhaps clustering was not well for these samples. "},{"metadata":{"_uuid":"acf4ec40b839e5bb386780d153b326be1c4646be"},"cell_type":"markdown","source":"## How many multilabel target  do the clusters hold?"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"3d5e3567d1cafd94222281fae54ebbe180b91be9"},"cell_type":"code","source":"results[\"number_of_targets\"] = results.drop(\"cluster\", axis=1).sum(axis=1)\n\nmultilabel_stats = results.groupby(\"cluster_names\").number_of_targets.value_counts() \nmultilabel_stats /= results.groupby(\"cluster_names\").number_of_targets.count()\nmultilabel_stats = multilabel_stats.unstack()\nmultilabel_stats.fillna(0, inplace=True)\nmultilabel_stats = 100 * multilabel_stats\nmultilabel_stats = multilabel_stats.apply(np.round)\nmultilabel_stats = multilabel_stats.astype(np.int)\n\nplt.figure(figsize=(20,5))\nsns.heatmap(multilabel_stats.transpose(),\n            square=True,\n            cbar=False,\n            cmap=\"Greens\",\n            annot=True);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"595f34f4d686d753032216ed4c5febd717a3dca2"},"cell_type":"markdown","source":"### Take-Away\n\n* First of all we can see that the low dense clusters all have multiple targets. This can be the reason why their samples were not assigned to other clusters that would fit but are different in the number of total targets.\n* In addition we can confirm what we already know: Many samples only have one or two target proteins present, not more! This is especially interesting for clusters that hold various components. Often these clusters have relations to the cellular devision process.\n* There are only two clusters that have at least two targets present:\n    * Nucleoli fibrillar center & Cytoplasmic bodies\n    * Rods & Rings, Microtubule ends, Nuclear bodies "},{"metadata":{"_uuid":"0570a6be51f203c2ba509f0c17832cb10c2cec83"},"cell_type":"markdown","source":"## Which clusters are anomalistic?"},{"metadata":{"trusted":true,"_uuid":"f5327fa483007516a8b69050a88d92fa8a53d92d"},"cell_type":"code","source":"sample_logLs = model.score_samples(X)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0044023b047f4e628c0fc0934397d8f92ec1d6ae"},"cell_type":"markdown","source":"For each sample we can try to find out how dense the region of the target feature space is, hence how many other samples are around it with the same kind of target structure. This density is given by the sample log-likelihood. "},{"metadata":{"trusted":true,"_uuid":"5f260e3dda24c18d0b0368d054dee706dd41d0ca"},"cell_type":"code","source":"my_threshold = np.quantile(sample_logLs, 0.05)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"093ac32d1897b06a26ede300a362184489240345"},"cell_type":"markdown","source":"To define which samples are anomal, I will setup a threshold of 5 %. If we compare this threshold with the distribution of the sample log-likelihoods, you can see..."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"3134785a2a57d0261c3628b85eb1b223c3685e5c"},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(sample_logLs)\nplt.axvline(my_threshold, color=\"Red\")\nplt.xlabel(\"Sample log likelihood of bernoulli mixture\")\nplt.title(\"Choosing a threshold to detect anomalies\")\nplt.ylabel(\"Density\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5455a083d9cd9d57f33703fcfb2b55a17feb1161"},"cell_type":"markdown","source":"... that below this threshold there are some samples in the negative regim (< -8)  that are truely outliers compared to the other ones (> -8)."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"6e58ceedfb851bb91a896302810a65cd4ef3e09d"},"cell_type":"code","source":"results[\"anomaly\"] = np.where(sample_logLs <= my_threshold, 1, 0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8c69b9339ae1db6117fbf9be556c3ffaa204784c"},"cell_type":"markdown","source":"### How anomalistic are the clusters?"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"9b35d1a8abd573a67acc59377c828a6892fa6425"},"cell_type":"code","source":"anomalies = results.groupby(\"cluster_names\").anomaly.value_counts() / results.groupby(\"cluster_names\").cluster.count() * 100\nanomalies = anomalies.unstack()\nanomalies.fillna(0, inplace=True)\nanomalies = anomalies.apply(np.round)\nanomalies = anomalies.astype(np.int)\n\nplt.figure(figsize=(20,5))\nsns.heatmap(anomalies.transpose(), cmap=\"Reds\", annot=True, square=True, cbar=False);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d50af42c5c86adf48219ff7a94f52ae93306c6be"},"cell_type":"markdown","source":"### Take-Away\n\n* I expected the low-dense clusters to be more anomalistic than the others but this seems to be only true for low-dense 3 and a bit for low dense 1.\n* Interestingly some other clusters have anomalies as well. This is especially the case for clusters that tend to have only one target protein but have some with 2 target proteins as well. How anomal a cluster may depend on the total number of targets of the samples in that cluster.\n* In contrast the various Rods&Rings cluster that have uncertain cluster assignments (look at the chapter below) is not located in a low dense region. Hence there are other clusters and targets around that could suite as well. "},{"metadata":{"_uuid":"c9de421db0286cffd5321678f52d588bdc17d0ae"},"cell_type":"markdown","source":"## How certain are the cluster assignments?"},{"metadata":{"trusted":true,"_uuid":"42403d4ed3b94a1059af9258cc5fd1f7100c6656"},"cell_type":"code","source":"results[\"certainty\"] = np.sort(model.gamma, axis=1)[:,-1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e5f9a9b672608695d626d81521b776b9b721b95d","_kg_hide-input":true},"cell_type":"code","source":"certainties = results.certainty.values\n\nplt.figure(figsize=(20,5))\nsns.distplot(certainties, color=\"Orange\")\nplt.xlabel(\"Certainty of cluster assignment\")\nplt.ylabel(\"Density\")\nplt.title(\"How sure was the model in predicting the winner?\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"75032ecaff025a7bfee073c80cfd410a39c1820b"},"cell_type":"markdown","source":"### Take-Away\n\n* The certainties vary on a broad range. This suggests that we might have further cluster interactions or outlier target samples.\n* Let's have a look at the certainty distribution and statistics per cluster:"},{"metadata":{"trusted":true,"_uuid":"900727bdfabcdcab4e95815d191589d0dc90cbee"},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.boxplot(x=\"cluster_names\", y=\"certainty\", data=results)\nplt.ylim([0,1])\nplt.xticks(rotation=90)\nplt.xlabel(\"\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"88caa972636f7773b43a036e8a6fee87ccd51850"},"cell_type":"markdown","source":"### Take-Away\n\n* Ah! Cool! Take a look at the various - RodsRings cluster. The model was very uncertain in its cluster assignments. If you like to take a look at the blue map, you can see that a lot of targets are present which correspond to the cellular devision process. I think we should consider the alternative clusters for this cluster, this are those that yielded a high probability to be responsible as well (that's gamma in the bernoulli model). The same should be done for the cluster that solely holds nuclear membranes.\n* In addition the occurence of seldom target proteins seems to make a cluster assignment more uncertain. The aggresome cluster for example holds 50 % of all aggresome proteins in the competition, but this cluster there is a more common target - nucleoplasmn - as well. If you take a look at the probability $\\mu$ for Aggresomes in that cluster it's quite low with 18. Perhaps this cluster is not nice shaped and occupied. We don't know, how many nucleoplasmn targets come alone in that cluster or how many are coupled with aggresomes. Let's make a map for this case: "},{"metadata":{"trusted":true,"_uuid":"d32aea578cf8c0e20b0d7b124c64fe7b394b43eb","_kg_hide-input":true},"cell_type":"code","source":"aggresome_cluster = results.loc[results.cluster==1].drop([\"cluster\",\n                                 \"cluster_names\",\n                                 \"number_of_targets\", \n                                  \"anomaly\", \n                                 \"certainty\"], axis=1).copy()\n\ncounts = aggresome_cluster.sum()\ncolumns_of_interest = list(counts[counts>0].index.values)\naggresome_cluster = aggresome_cluster.loc[:,columns_of_interest]\n\naggresome_combinations = pd.DataFrame(index=aggresome_cluster.columns.values,\n                                      columns=aggresome_cluster.columns.values)\n\nfor col in aggresome_combinations.columns:\n    aggresome_combinations.loc[col,:] = aggresome_cluster[aggresome_cluster[col] == 1].sum()\n\nmask = np.zeros_like(aggresome_combinations, dtype=np.bool)\nmask[np.triu_indices_from(mask, k=1)] = True\n\n\nplt.figure(figsize=(8,8))\nsns.set(style=\"white\")\nsns.heatmap(aggresome_combinations, mask=mask, cmap=\"Reds\",\n            square=True, linewidths=.5, cbar_kws={\"shrink\": .5}, vmin=0, vmax=50, annot=True, fmt=\"g\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"82340fcfe22eb27aec94a2ce1abbf91c590cf1b8"},"cell_type":"markdown","source":"### Take-Away\n\n* Ok, a lot of aggresomes seem to come alone as they have a lot of counts with themselves but only few counts for other target proteins.\n* Samples with two targets present are probably given by Nucleoplasmn-Aggresome combinations as this is the second highest count in the map.\n* By only considering the aggresome row we can see that there are some samples that are coupled with cytokinetic bridge or cell junctions with have themselves some couplings with the nucleoplasm. This cases are the 3-couped-targets of that cluster. \n* As there are no 4-coupled-targets we can conclude that there are some combinations that do not really fit to the aggresome cluster. One example: Take a look at the aggresome row again and the microtubule organizing center column. You can see that there are zero couplings of aggresomes with this target. But you can see that the organizing center has couplings with cell junctions and nucleoplasmn. Consequently some of the 2- or 3-couplings could be given by such combinations without aggresomes. "},{"metadata":{"trusted":true,"_uuid":"895be19bea95eda0f77c464ab3b023f1830129b2"},"cell_type":"code","source":"aggresome_cluster.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"64b5e264273f46acbb9221bbfebb3d62a9c5f148"},"cell_type":"markdown","source":"Considering the total number of samples in the cluster given by 183, we can say, that there are 7 samples that do not have an aggresome!"},{"metadata":{"_uuid":"b7ced05fc6a33ef4a965b9029b4a6813db04ace1"},"cell_type":"markdown","source":"## Which kind of cluster interactions do we have?"},{"metadata":{"trusted":true,"_uuid":"67e92d15e656e8c03f3f9a0392d764fae170775d","_kg_hide-input":true},"cell_type":"code","source":"results[\"alternative_cluster\"] = np.argsort(model.gamma, axis=1)[:,-2]\nresults[\"alternative_names\"] = results.alternative_cluster.apply(lambda l: cluster_names[l])\nresults[\"alternative_certainties\"] = np.sort(model.gamma, axis=1)[:,-2]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aa04cf608a2deb42afea357508c4704072a0616f","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"competition = np.round(100 * results.groupby(\n    \"cluster_names\").alternative_names.value_counts() / results.groupby(\n    \"cluster_names\").alternative_names.count())\ncompetition = competition.unstack()\ncompetition.fillna(0, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fd9d95ff639038634a666be54019f2c8593b861c","_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(20,15))\nsns.heatmap(competition, cmap=\"Greens\", annot=True, fmt=\"g\", square=True, cbar=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"106385b524965efdc2ca0b9dbca6c73285c9e665"},"cell_type":"markdown","source":"### Take-Away\n\n* Let's only consider some example clusters that have low cluster assignment certainties. \n    * The various Rods & Rings cluster has a tendency to join with the nucleoli fibrillar center & cytoplasmic bodies cluster.\n    * Aggresomes like to join with the mitotic spindle & organizing center cluster."},{"metadata":{"_uuid":"5109a227a474a7a6f1531c5d067917f129e81d31"},"cell_type":"markdown","source":"## Conclusion\n\nWithin this kernel you can find groups of targets that are likely to occur together. It covers many different topics that can help you to solve parts of the imbalance class problem or to improve target predictions. It's your turn now to figure out how to do this! ;-)\n\nHappy coding!"},{"metadata":{"_uuid":"230f740a1871d6806c1ac137f96e1a875c11983f","trusted":true},"cell_type":"code","source":"results.to_csv(\"target_group_analysis.csv\")","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}