{"cells":[{"metadata":{"_uuid":"15d0a3b8abd040c61feb13b63d9f5fb9ead6a121"},"cell_type":"markdown","source":"# Exploratory data analysis of the human protein atlas image dataset\nupdate 5/10/2018: beginning of cell segmentation algorithm\n\nupdate 5/10/2018: add red + blue channels stack and whole cell identification (does not give a clean result, though)\n\nThis kernel is just the beginning of a work in progress and will be updated very often.\nWe will explore the dataset available for the human protein atlas image competition. Questions we would like to answer include:\n* what channels of the image contain the relevant information\n* how much can we reduce dimensionality of data while retaining important information"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#import modules\n\nimport 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\"))","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"## What's in the data?\n Let's import the *train.csv* data files to see what they contain. We also define a dictionary containing the map between labels of the training data (the column *target* in *train.csv*) and their biological meaning."},{"metadata":{"trusted":true,"_uuid":"399f7cb6c2467d25a8d3c2347ef4c8bfd5d85421"},"cell_type":"code","source":"#import training data\ntrain = pd.read_csv(\"../input/train.csv\")\nprint(train.head())\n\n#map of targets in a dictionary\nsubcell_locs = {\n0:  \"Nucleoplasm\", \n1:  \"Nuclear membrane\",   \n2:  \"Nucleoli\",   \n3:  \"Nucleoli fibrillar center\" ,  \n4:  \"Nuclear speckles\",\n5:  \"Nuclear bodies\",\n6:  \"Endoplasmic reticulum\",   \n7:  \"Golgi apparatus\",\n8:  \"Peroxisomes\",\n9:  \"Endosomes\",\n10:  \"Lysosomes\",\n11:  \"Intermediate filaments\",   \n12:  \"Actin filaments\",\n13:  \"Focal adhesion sites\",   \n14:  \"Microtubules\",\n15:  \"Microtubule ends\",   \n16:  \"Cytokinetic bridge\",   \n17:  \"Mitotic spindle\",\n18:  \"Microtubule organizing center\",  \n19:  \"Centrosome\",\n20:  \"Lipid droplets\",   \n21:  \"Plasma membrane\",   \n22:  \"Cell junctions\", \n23:  \"Mitochondria\",\n24:  \"Aggresome\",\n25:  \"Cytosol\",\n26:  \"Cytoplasmic bodies\",   \n27:  \"Rods & rings\" \n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a2b30702cb0b918e51a8dd5f21b5b41051559612"},"cell_type":"markdown","source":"Each image is a 4-channel image with the protein of interest in the green channel. It is the subcellular localization of this protein which is recorded in the *Target* column of the *train.csv* file. The red channel corresponds to microtubules, the blue channel to the nucleus and the yellow channel to the endoplasmid reticulum. Let's display the different channels of the image with ID == 1, since it contains several subcelullar locations for our protein of interest. Then we will overlay the green and yellow channel, as the yellow channel gives a good indication of the cell shape."},{"metadata":{"trusted":true,"_uuid":"c29862332cd7c89791044c56984a03331951fabd"},"cell_type":"code","source":"print(\"The image with ID == 1 has the following labels:\", train.loc[1, \"Target\"])\nprint(\"These labels correspond to:\")\nfor location in train.loc[1, \"Target\"].split():\n    print(\"-\", subcell_locs[int(location)])\n\n#reset seaborn style\nsns.reset_orig()\n\n#get image id\nim_id = train.loc[1, \"Id\"]\n\n#create custom color maps\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/train/{}_green.png'.format(im_id), 0)\nred = cv2.imread('../input/train/{}_red.png'.format(im_id), 0)\nblue = cv2.imread('../input/train/{}_blue.png'.format(im_id), 0)\nyellow = cv2.imread('../input/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,"_uuid":"d6f0ce4a0b1c4bba35e14b14df8a5f28cbc83e70"},"cell_type":"code","source":"#stack nucleus and microtubules images\n#create blue nucleus and red microtubule images\nnuclei = cv2.merge((np.zeros((512, 512),dtype='uint8'), np.zeros((512, 512),dtype='uint8'), blue))\nmicrotub = cv2.merge((red, np.zeros((512, 512),dtype='uint8'), np.zeros((512, 512),dtype='uint8')))\n\n#create ROI\nrows, cols, _ = nuclei.shape\nroi = microtub[:rows, :cols]\n\n#create a mask of nuclei and invert mask\nnuclei_grey = cv2.cvtColor(nuclei, cv2.COLOR_BGR2GRAY)\nret, mask = cv2.threshold(nuclei_grey, 10, 255, cv2.THRESH_BINARY)\nmask_inv = cv2.bitwise_not(mask)\n\n#make area of nuclei in ROI black\nred_bg = cv2.bitwise_and(roi, roi, mask=mask_inv)\n\n#select only region with nuclei from blue\nblue_fg = cv2.bitwise_and(nuclei, nuclei, mask=mask)\n\n#put nuclei in ROI and modify red\ndst = cv2.add(red_bg, blue_fg)\nmicrotub[:rows, :cols] = dst\n\n#show result image\nfig, ax = plt.subplots(figsize=(8, 8))\nax.imshow(microtub)\nax.set_title(\"Nuclei (blue) + microtubules (red)\", fontsize=15)\nax.set_xticklabels([])\nax.set_yticklabels([])\nax.tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d8d513f62262f86485a8a2f3a4cb2338621e126b"},"cell_type":"markdown","source":"Let's see how the targets are distributed."},{"metadata":{"trusted":true,"_uuid":"af5a54796ad2566eb9d13b82da7e139b7a70cf6c"},"cell_type":"code","source":"labels_num = [value.split() for value in train['Target']]\nlabels_num_flat = list(map(int, [item for sublist in labels_num for item in sublist]))\nlabels = [\"\" for _ in range(len(labels_num_flat))]\nfor i in range(len(labels_num_flat)):\n    labels[i] = subcell_locs[labels_num_flat[i]]\n\nfig, ax = plt.subplots(figsize=(15, 5))\npd.Series(labels).value_counts().plot('bar', fontsize=14)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2d73cfab246649aa37f09b750003122ee35bbee7"},"cell_type":"markdown","source":"According to [Chen *et al*. 2007](https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btm206), if images are segmented into single cell regions, additional features that are not appropriate for whole fields can be calculated after *seeded watershed segmentation*. Nucleus images provide a means to identify each cell, so image segmentation may start by identification of nuclei in images.  The function `cv2.connectedComponents` provides a simple and effective means to label nuclei in images. Conversely, as shown on the following notebook cell, identification of whole cells using `cv2.connectedComponents` is not as efficient, due to the less homogeneous signal in the yellow channel of the image."},{"metadata":{"trusted":true,"_uuid":"e92502b354d1a92a01f2afda8fd1fdb647f41be4"},"cell_type":"code","source":"#apply threshold on the nucleus image\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,"_uuid":"332fde916dc7f91baeae8875985a35f972ba1e1f"},"cell_type":"code","source":"#apply threshold on the endoplasmic reticulum image\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":{"_uuid":"2e2c581d22e9e9d0d47857dab3b7f05ad770ff8b"},"cell_type":"markdown","source":"Let's try different simple thresholding methods. Description of threshold types can be found [here](https://docs.opencv.org/3.4/d7/d4d/tutorial_py_thresholding.html) and [here](https://docs.opencv.org/3.4/d7/d1b/group__imgproc__misc.html#gaa9e58d2860d4afa658ef70a9b1115576)."},{"metadata":{"trusted":true,"_uuid":"00ea8d2bb3097e312345e03e5e48e42689868cba"},"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\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":{"_uuid":"4e9aa538dc6c5a6aacd9535226bfca2adbdefe70"},"cell_type":"markdown","source":"*To zero* simple thresholding is not adapted at all for identifying cell boundaries based on the yellow channel. Even after playing with the upper and lower parameter values, no satisfactory result is obtained. *Binary* and *truncate* methods work better. Let's see how *connectedComponents* work after both thresholding methods."},{"metadata":{"trusted":true,"_uuid":"79311045b4cb1c09db378ab0743c7758a051c72b"},"cell_type":"code","source":"fig, 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)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc8cb9295710299b8e605599d42fd0567e99e857"},"cell_type":"markdown","source":"At this point it's not clear if truncate thresholding is an improvement compared to binary thresholding. Some cells are fused to each other while they should not be.\n\nOn the other hand. Adaptive thresholding methods apply a different threshold on different parts of the image, let's see how well it does on our images. See [here](https://docs.opencv.org/3.4/d7/d1b/group__imgproc__misc.html#gaa42a3e6ef26247da787bf34030ed772c) for more explanations."},{"metadata":{"trusted":true,"_uuid":"694bf4d8a26991e12cc982c9bffa3be75f3c683d"},"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\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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"637e61a361d86422375239f4e8c78a420446adcf"},"cell_type":"code","source":"","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}