\n
Quick Navigation
\n\n* [0. Brief about single cell type atlas](#0)\n \n* [1. Feature Engineering](#1)\n \n* [2. EDA](#2)\n \n* [3. Image Segmentation](#3)\n \n* [4. Bernoulli Mixture](#4)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"
\n
Brief about Single Cell Type Atlas
"},{"metadata":{},"cell_type":"markdown","source":"> The **Single Cell Type Atlas** contains single cell **RNA sequencing (scRNAseq)** data from 13 different human tissues, together with in-house generated immunohistochemically stained tissue sections visualizing the corresponding spatial protein expression patterns. The scRNAseq analysis was based on publicly available genome-wide expression data and comprises all protein-coding genes in 192 individual cell type clusters corresponding to 12 different cell type groups. A specificity and distribution classification was performed to determine the number of genes elevated in these single cell types, and the number of genes detected in one, several or all cell types, respectively.\n\n"},{"metadata":{},"cell_type":"markdown","source":"
Importing dependencies
"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport cv2\nfrom PIL import Image\nfrom collections import Counter\n\nimport os\nprint(os.listdir(\"../input/hpa-single-cell-image-classification\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(\"../input/hpa-single-cell-image-classification/train.csv\")\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n
Feature Engineering
"},{"metadata":{},"cell_type":"markdown","source":"Mapping of targets in a location:\n---"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nlabel_names= {\n0: \"Nucleoplasm\",\n1: \"Nuclear membrane\",\n2: \"Nucleoli\",\n3: \"Nucleoli fibrillar center\",\n4: \"Nuclear speckles\",\n5: \"Nuclear bodies\",\n6: \"Endoplasmic reticulum\",\n7: \"Golgi apparatus\",\n8: \"Intermediate filaments\",\n9: \"Actin filaments\",\n10: \"Microtubules\",\n11: \"Mitotic spindle\",\n12: \"Centrosome\",\n13: \"Plasma membrane\",\n14: \"Mitochondria\",\n15: \"Aggresome\",\n16: \"Cytosol\",\n17: \"Vesicles and punctate cytosolic patterns\",\n18: \"Negative\"\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.Label.replace(\"|\", \" \").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":{},"cell_type":"markdown","source":"\n---\n**Each sample consists of four image files. Each file represents a different filter on the subcellular protein patterns represented by the sample (ID).**\n\n> **Red for Microtubule channels.**\n\n> **Blue for Nuclei channels.**\n\n> **Yellow for Endoplasmic Reticulum (ER) channels.**\n\n> **Green for Protein of interest.**\n\n---"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"The image with ID == 0 has the following labels:\", train_df.loc[0, \"Label\"])\nprint(\"These labels correspond to:\")\nfor location in train_df.loc[0, \"Label\"].replace(\"|\", \" \").split():\n print(\"-\", label_names[int(location)])\n\n#reset seaborn style\nsns.reset_orig()\n\n#get image id\nim_id = train_df.loc[1, \"ID\"]\n\n\ncdict1 = {'red': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0)),\n\n 'green': ((0.0, 0.0, 0.0),\n (0.75, 1.0, 1.0),\n (1.0, 1.0, 1.0)),\n\n 'blue': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0))}\n\ncdict2 = {'red': ((0.0, 0.0, 0.0),\n (0.75, 1.0, 1.0),\n (1.0, 1.0, 1.0)),\n\n 'green': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0)),\n\n 'blue': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0))}\n\ncdict3 = {'red': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0)),\n\n 'green': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0)),\n\n 'blue': ((0.0, 0.0, 0.0),\n (0.75, 1.0, 1.0),\n (1.0, 1.0, 1.0))}\n\ncdict4 = {'red': ((0.0, 0.0, 0.0),\n (0.75, 1.0, 1.0),\n (1.0, 1.0, 1.0)),\n\n 'green': ((0.0, 0.0, 0.0),\n (0.75, 1.0, 1.0),\n (1.0, 1.0, 1.0)),\n\n 'blue': ((0.0, 0.0, 0.0),\n (1.0, 0.0, 0.0))}\n\nplt.register_cmap(name='greens', data=cdict1)\nplt.register_cmap(name='reds', data=cdict2)\nplt.register_cmap(name='blues', data=cdict3)\nplt.register_cmap(name='yellows', data=cdict4)\n\n#get each image channel as a greyscale image (second argument 0 in imread)\ngreen = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_green.png'.format(im_id), 0)\nred = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_red.png'.format(im_id), 0)\nblue = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_blue.png'.format(im_id), 0)\nyellow = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_yellow.png'.format(im_id), 0)\n\n#display each channel separately\nfig, ax = plt.subplots(nrows = 2, ncols=2, figsize=(15, 15))\nax[0, 0].imshow(green, cmap=\"greens\")\nax[0, 0].set_title(\"Protein of interest\", fontsize=18)\nax[0, 1].imshow(red, cmap=\"reds\")\nax[0, 1].set_title(\"Microtubules\", fontsize=18)\nax[1, 0].imshow(blue, cmap=\"blues\")\nax[1, 0].set_title(\"Nucleus\", fontsize=18)\nax[1, 1].imshow(yellow, cmap=\"yellows\")\nax[1, 1].set_title(\"Endoplasmic reticulum\", fontsize=18)\nfor i in range(2):\n for j in range(2):\n ax[i, j].set_xticklabels([])\n ax[i, j].set_yticklabels([])\n ax[i, j].tick_params(left=False, bottom=False)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for key in label_names.keys():\n train_df[label_names[key]] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df.apply(fill_targets, axis=1)\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n
EDA
"},{"metadata":{"trusted":true},"cell_type":"code","source":"target_counts = train_df.drop([\"ID\", \"Label\"],axis=1).sum(axis=0).sort_values(ascending=False)\nsns.set_style(\"darkgrid\")\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":{},"cell_type":"markdown","source":"Take away:\n---\n> **We can see that most common protein structures belong to coarse grained cellular components like the plasma membrane, the cytosol and the nucleus.**"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[\"number_of_targets\"] = train_df.drop([\"ID\",\"Label\"],axis=1).sum(axis=1)\ncount_perc = np.round(100 * train_df[\"number_of_targets\"].value_counts() / train_df.shape[0], 2)\nplt.figure(figsize=(15,5))\nsns.set_style(\"darkgrid\")\nsns.barplot(x=count_perc.index.values, y=count_perc.values, palette=\"Reds\")\nplt.xlabel(\"Number of targets per image\")\nplt.ylabel(\"% of train data\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Take away:\n---\n> **Most train images only have 1 or two target labels.**\n\n> **More than 3 targets are sporadic!**"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,15))\n\nsns.heatmap(train_df[train_df.number_of_targets>1].drop([\"ID\", \"Label\", \"number_of_targets\"],axis=1).corr(), annot=True,cmap=\"viridis\", vmin=-1, vmax=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Take away:\n---\n\n> We can see that many targets only have ***very slight*** correlations.\n\n>In addition we find that the **Mitotic spindle** often comes together with the **Microtubules**. This makes sense as both are participants for cellular division. And in this process microtubules and thier ends are active and participate as well. Consequently we find a positive correlation between these targets.\n"},{"metadata":{},"cell_type":"markdown","source":"
\n
Image Segmentation
"},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the nucleus image\nsns.set_style(\"white\")\nret, thresh = cv2.threshold(blue, 0, 255, cv2.THRESH_BINARY)\n#display threshold image\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(thresh, cmap=\"Greys\")\nax[0].set_title(\"Threshold\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise\nkernel = np.ones((5,5),np.uint8)\nopening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)\nax[1].imshow(opening, cmap=\"Greys\")\nax[1].set_title(\"Morphological opening\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n# Marker labelling\nret, markers = cv2.connectedComponents(opening)\n# Map component labels to hue val\nlabel_hue = np.uint8(179 * markers / np.max(markers))\nblank_ch = 255 * np.ones_like(label_hue)\nlabeled_img = cv2.merge([label_hue, blank_ch, blank_ch])\n# cvt to BGR for display\nlabeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img[label_hue==0] = 0\nax[2].imshow(labeled_img)\nax[2].set_title(\"Markers\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the endoplasmic reticulum image\nsns.set_style(\"white\")\nret, thresh = cv2.threshold(yellow, 4, 255, cv2.THRESH_BINARY)\n#display threshold image\nfig, ax = plt.subplots(ncols=4, figsize=(20, 20))\nax[0].imshow(thresh, cmap=\"Greys\")\nax[0].set_title(\"Threshold\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise\nkernel = np.ones((5,5),np.uint8)\nopening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)\nax[1].imshow(opening, cmap=\"Greys\")\nax[1].set_title(\"Morphological opening\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n#morphological closing\nclosing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernel)\nax[2].imshow(closing, cmap=\"Greys\")\nax[2].set_title(\"Morphological closing\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)\n\n# Marker labelling\nret, markers = cv2.connectedComponents(closing)\n# Map component labels to hue val\nlabel_hue = np.uint8(179 * markers / np.max(markers))\nblank_ch = 255 * np.ones_like(label_hue)\nlabeled_img = cv2.merge([label_hue, blank_ch, blank_ch])\n# cvt to BGR for display\nlabeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img[label_hue==0] = 0\nax[3].imshow(labeled_img)\nax[3].set_title(\"Markers\", fontsize=15)\nax[3].set_xticklabels([])\nax[3].set_yticklabels([])\nax[3].tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the endoplasmic reticulum image\nret, thresh1 = cv2.threshold(yellow, 4, 255, cv2.THRESH_BINARY)\nret, thresh2 = cv2.threshold(yellow, 4, 255, cv2.THRESH_TRUNC)\nret, thresh3 = cv2.threshold(yellow, 4, 255, cv2.THRESH_TOZERO)\n\n#display threshold images\nsns.set_style(\"white\")\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(thresh1, cmap=\"Greys\")\nax[0].set_title(\"Binary\", fontsize=15)\n\nax[1].imshow(thresh2, cmap=\"Greys\")\nax[1].set_title(\"Trunc\", fontsize=15)\n\nax[2].imshow(thresh3, cmap=\"Greys\")\nax[2].set_title(\"To zero\", fontsize=15)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set_style(\"white\")\nfig, ax = plt.subplots(ncols=4, figsize=(20, 20))\n\n#morphological opening to remove noise after binary thresholding\nkernel = np.ones((5,5),np.uint8)\nopening1 = cv2.morphologyEx(thresh1, cv2.MORPH_OPEN, kernel)\nax[0].imshow(opening1, cmap=\"Greys\")\nax[0].set_title(\"Morphological opening (binary)\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological closing after binary thresholding\nclosing1 = cv2.morphologyEx(opening1, cv2.MORPH_CLOSE, kernel)\nax[1].imshow(closing1, cmap=\"Greys\")\nax[1].set_title(\"Morphological closing (binary)\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise after truncate thresholding\nkernel = np.ones((5,5),np.uint8)\nopening2 = cv2.morphologyEx(thresh2, cv2.MORPH_OPEN, kernel)\nax[2].imshow(opening2, cmap=\"Greys\")\nax[2].set_title(\"Morphological opening (truncate)\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)\n\n#morphological closing after truncate thresholding\nclosing2 = cv2.morphologyEx(opening2, cv2.MORPH_CLOSE, kernel)\nax[3].imshow(closing2, cmap=\"Greys\")\nax[3].set_title(\"Morphological closing (truncate)\", fontsize=15)\nax[3].set_xticklabels([])\nax[3].set_yticklabels([])\nax[3].tick_params(left=False, bottom=False)\n\nfig, ax = plt.subplots(ncols=2, figsize=(10, 10))\n# Marker labelling for binary thresholding\nret, markers1 = cv2.connectedComponents(closing1)\n# Map component labels to hue val\nlabel_hue1 = np.uint8(179 * markers1 / np.max(markers1))\nblank_ch1 = 255 * np.ones_like(label_hue1)\nlabeled_img1 = cv2.merge([label_hue1, blank_ch1, blank_ch1])\n# cvt to BGR for display\nlabeled_img1 = cv2.cvtColor(labeled_img1, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img1[label_hue1==0] = 0\nax[0].imshow(labeled_img1)\nax[0].set_title(\"Markers (binary)\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n# Marker labelling for truncate thresholding\nret, markers2 = cv2.connectedComponents(closing2)\n# Map component labels to hue val\nlabel_hue2 = np.uint8(179 * markers2 / np.max(markers2))\nblank_ch2 = 255 * np.ones_like(label_hue2)\nlabeled_img2 = cv2.merge([label_hue2, blank_ch2, blank_ch2])\n# cvt to BGR for display\nlabeled_img2 = cv2.cvtColor(labeled_img2, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img2[label_hue2==0] = 0\nax[1].imshow(labeled_img2)\nax[1].set_title(\"Markers (truncate)\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply adaptive threshold on endoplasmic reticulum image\ny_blur = cv2.medianBlur(yellow, 3)\n\n#apply adaptive thresholding\nret,th1 = cv2.threshold(y_blur, 5,255, cv2.THRESH_BINARY)\n\nth2 = cv2.adaptiveThreshold(y_blur, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 15, 3)\n\nth3 = cv2.adaptiveThreshold(y_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 15, 3)\n\n#display threshold images\nsns.set_style(\"white\")\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(th1, cmap=\"Greys\")\nax[0].set_title(\"Binary\", fontsize=15)\n\nax[1].imshow(th2, cmap=\"Greys_r\")\nax[1].set_title(\"Adaptive: mean\", fontsize=15)\n\nax[2].imshow(th3, cmap=\"Greys_r\")\nax[2].set_title(\"Adaptive: gaussian\", fontsize=15)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"
\n
Bernoulli Mixture
"},{"metadata":{},"cell_type":"markdown","source":"Bernoulli Mixture:\n---\n> **Bernoulli mixture models to solve different types of problems in pattern recognition like feature selection, classification, dimensionality reduction and rule generation.**\n\n> **Here We're going to use Bernoulli Mixture Model to see `How are specific proteins distributed over clusters` in the dataset.**\n----\n\nModel 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"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.special import logsumexp\n\n#Source: https://github.com/bayespy/bayespy\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":{"trusted":true},"cell_type":"code","source":"targets = train_df.drop([\"ID\", \"Label\", \"number_of_targets\"], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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":{"trusted":true},"cell_type":"code","source":"components_to_test = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = BernoulliMixture(estimated_components, 200)\nmodel.fit(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results = targets.copy()\nresults[\"cluster\"] = np.argmax(model.gamma, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":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":{},"cell_type":"markdown","source":"Take away:\n---\n> **This looks great! You can see that several clusters only hold one specific target protein!**\n\n> **For each target protein you can see the percentage of its occurences that are placed into specific clusters.**\n\n> **One example: 92 % of Microtubules target proteins are located in cluster 7. Only a few percents are hold by cluster 5, 14, 22 and 27. There is one percent missing to fill up to 100 % but this is caused by rounding errors and should not worry you.**"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Reference:\n---\n> https://www.kaggle.com/allunia/protein-atlas-exploration-and-baseline/\n\n> https://docs.opencv.org/\n\n> https://github.com/bayespy/bayespy\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}