{"cells":[{"metadata":{"trusted":true,"_uuid":"eac732664992ec7aa7f1f5e653e6615e89a45315"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport missingno as msn \nimport seaborn as sns\nsns.set(style=\"whitegrid\")\nimport matplotlib.pylab as plt\n%pylab inline\n\nPATH = '../input/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b8c61051b2371a2b6f654eb5be2da4db1d985a17"},"cell_type":"code","source":"label = {\n\"0\" : \"Nucleoplasm\", \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\nreversed_label = dict()\nfor i in label.keys():\n    reversed_label[label[i]] = i\n    \nlabeled_columns = [i for i in reversed_label.keys()]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5d0a2d524984de1daf8d6e4532cb4c8c0a6c17f0"},"cell_type":"code","source":"df_train_labels = pd.read_csv(PATH+'train.csv', sep=',')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2aec9b05f31f2657a588f297455e6fd3ffd976df"},"cell_type":"code","source":"display(df_train_labels.describe(include='all'))\n\ndisplay(df_train_labels.sample(5))\n\nfor i in range(25):\n    loca = [t for t in df_train_labels[\"Target\"].value_counts()[:25].index]\n    loca_name = []\n    for n in loca:\n        temp = [label[j] for j in n.split(' ')]\n        loca_name.append(temp)\n    counts = df_train_labels[\"Target\"].value_counts()[i]  \n        \n    print(\"{} occurs {} times in data.\".format((\", \").join(loca_name[i]), counts))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e4b25f5c55ed6c468e4036ea0b9d795e8619dbbb"},"cell_type":"code","source":"for j in label.values():\n    df_train_labels[j] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0971c8fa4bc8047f8a78affea1e083ca836e23d8"},"cell_type":"code","source":"def fill_rows(row):\n    #function which fills dataframe based on target label\n    for i in row[\"Target\"].split(\" \"):\n        name = label[i]\n        row.loc[name] = 1\n    return row\n        \ndf_train_labels = df_train_labels.apply(fill_rows, axis=1)\n\ndf_train_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"647ce334a63531df371a8d1e87cf959e30bd2911"},"cell_type":"code","source":"df_train_labels[labeled_columns].sum(0).sort_values(ascending=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"646a7bd19eaf2864320f264f44429ea22fe69caf"},"cell_type":"markdown","source":"### Nucleoplasm occurs most frequently in the data. When comparing the single locations with the combined locations one can see that there must be a frequent co-labeling of the nucleoplasm with other organelles. Those might be other nuclear sites like nucleoli or nuclear speckles but also others for example the most frequent co-labeling occurs with cytosol. "},{"metadata":{"trusted":true,"_uuid":"446c7217c91c65187889dfa1aa75f9df2dacb0dd"},"cell_type":"code","source":"plt.figure(figsize=(7,7))\nax = plt.subplot()\nplt.title(\"Fig1. # annotations for image\")\nplt.xlabel(\"# annot\")\nax.spines[\"top\"].set_visible(False)\nax.spines[\"right\"].set_visible(False)\nplt.grid(False)\nsns.countplot(df_train_labels.sum(1))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"45a2ee0de2cb5e7a623359720e6f1c849835483c"},"cell_type":"markdown","source":"### Huge differences in the label. Have to come up with appropiate sampeling strategy. Going to try a combination of undersampling high count labels and artificial increasing amoung of samples for low count labels (i.e. image augmentation)"},{"metadata":{"trusted":true,"_uuid":"5a0e2562c4f9da8f8d7e913634ee483f8ec16d47"},"cell_type":"code","source":"plt.figure(figsize=(10,10))\nplt.title(\"Fig2. Correlation of label\")\nsns.heatmap(df_train_labels[labeled_columns].corr(), cmap=\"PiYG\", linewidths=.05,\n           linecolor='b',square=True)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"dabfd1a713b67404952b3014b8f304c23e171d33"},"cell_type":"code","source":"plt.figure(figsize=(25,15))\nfor i, loca in enumerate(labeled_columns):\n    if i==3:\n        plt.title(\"How likely is the co-labeling for a certain compartment\")\n    ax = plt.subplot(5,6,i+1)\n    ax.spines[\"top\"].set_visible(False)\n    ax.spines[\"right\"].set_visible(False)\n    plt.grid(False)\n    only = df_train_labels[(df_train_labels[loca] == 1) & \n                   (df_train_labels.sum(1) == 1)].shape[0]\n    co_labeled = df_train_labels[(df_train_labels[loca] == 1) & \n                   (df_train_labels.sum(1) > 1)].shape[0]\n\n    plt.title(loca)\n    plt.bar(x=[0,1], height=[only,co_labeled], tick_label=[\"only\", \"colabeled\"], color=[\"purple\", \"green\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"330d227d26fe0f3b66254d54348ff0d004258e7b"},"cell_type":"code","source":"plt.figure(figsize=(7,7))\nplt.title(\"Fig3. Barplot of label counts\")\nsns.barplot(x=df_train_labels[labeled_columns].sum(0).sort_values(ascending=False).index,\n           y=df_train_labels[labeled_columns].sum(0).sort_values(ascending=False))\nplt.xticks(rotation=90);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e64f0b52212b2aa50fd99de1eae9684b963ddb77"},"cell_type":"markdown","source":"# EDA summary\n<p>We are dealing with a training dataset which contains 31072 unique entries. Those entries are labeled with an id which relates to the respective images and a 'Target' column. 'Target' column shows the subcellular localization of an protein of interest (POI). A POI can be located at more than one subcellular compartment. However, localization at more than two compartments becomes increasingly unlikely (Fig1). The investigation of correaltion of certain subcellular compartments shows that endosomes and lysosomes are highly correlated (Fig2). This is not further suprising since they look very similiar. Also cytokinetic bridge and mictrotubules/microtubule ends show a positive correlation. We will look at this in a little bit more detail when investigating the actual images. The frequency of the labels is highly variable. Nucleoplasm is by far the most common label, while Peroxisomes, Endosomes, Lysosomes, Microtubule ends and Rods&Rings are scarce (Fig3). In order to not over- or underrespresent those labels an appropiate sampling strategy is mandatory.</p>"},{"metadata":{"trusted":true,"_uuid":"c657f4f7d335e0c7e84d68db0764336780542376"},"cell_type":"code","source":"import os, sys\nimport cv2\nimport gc\nimport random","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6881c9ed04621d89968f2dcdbbecccc5a4b21d85"},"cell_type":"code","source":"def show_random_img(noi = 1):\n    #Function  that shows a given number of random imgs\n    \n    colors = [\"_blue.png\", \"_green.png\", \"_red.png\", \"_yellow.png\"]\n    cmaps =[\"Blues\", \"Greens\", \"Oranges\", \"Reds\"]\n    \n    rnd_imgs = []\n    for i in range(noi):\n        rnd = random.randint(0,len(os.listdir(PATH+\"train\")))\n        rnd_imgs.append(os.listdir(PATH+\"train\")[rnd].split(\"_\")[0])\n    \n    for rnd_img in rnd_imgs:\n        targets = (df_train_labels[\"Target\"][df_train_labels[\"Id\"].str.contains(rnd_img)])\n        for i in targets.iteritems():\n            l = [label[j] for j in i[1].split(\" \")]\n\n        plt.figure(figsize=(20,10))\n        for j,color in enumerate(colors):\n            plt.subplot(1,4,j+1)\n            if j == 0:\n                plt.title(\"Nucleus\")\n            if j == 1:\n                plt.title(l)\n            if j == 2:\n                plt.title(\"Microtubules\")\n            if j == 3:\n                plt.title(\"ER\")\n            plt.grid(False)\n            img = cv2.imread(PATH+\"train/\"+rnd_img+color, 0)\n            plt.imshow(img, cmap=cmaps[j])\n","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"c5ac4f89657f35a667d45b83146e7a8cfcc00aa3"},"cell_type":"code","source":"show_random_img(4)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e8870a3148381eefec340b6a89302dc793d8be41"},"cell_type":"markdown","source":"## Looking at a couple of images shows that images were aquired using various magnifications. Also density of the cells and the overall cellshape varies. Different cell types can vastly differ (see HumanProteinAtlas ENSG00000167552-TUBA1A) in the content of tubulins (which are the major constituent of microtubule (e.g. uniprot/Q71U36)). These differences could possibly be problematic when using those as references in order to determine the subcellular localization of the POI. Might be worth to cluster images based on the auxiliary channels prior to training.  \n## Lets further take a quick look at an overlay of the channels to appreciate nature's beauty :)"},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"44543ee927ef3a5d977aab18a55e52542e7d431d"},"cell_type":"code","source":"rnd = random.randint(0,len(os.listdir(PATH+\"train\")))\nrnd_img = os.listdir(PATH+\"train\")[rnd].split(\"_\")[0]\n\nnuc = cv2.imread(PATH+\"train/\"+rnd_img+\"_blue.png\", 0)\nnuc = cv2.cvtColor(nuc, cv2.COLOR_GRAY2BGR)\npoi = cv2.imread(PATH+\"train/\"+rnd_img+\"_green.png\", 0)\npoi = cv2.cvtColor(poi, cv2.COLOR_GRAY2BGR)\nmt = cv2.imread(PATH+\"train/\"+rnd_img+\"_yellow.png\", 0)\nmt = cv2.cvtColor(mt, cv2.COLOR_GRAY2BGR)\ner = cv2.imread(PATH+\"train/\"+rnd_img+\"_red.png\", 0)\ner = cv2.cvtColor(er, cv2.COLOR_GRAY2BGR)\n\nnuc[:,:,:2] = 0\npoi[:,:,0] = 0\nmt[:,:,2] = 0\ner[:,:,1:] = 0\n\nimg3 = (0.25*(nuc/255) + 0.25*(poi/255) + 0.25*(mt/255) + 0.25 *(er/255))\nplt.figure(figsize=(8,8))\nplt.grid(False)\nplt.imshow(3*img3)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"f65a42e7dd5f86a11843903d4127ab45f4e3d2e5"},"cell_type":"code","source":"plt.figure(figsize=(10,10))\nplt.title(\"Nucleus and microtubule\")\nplt.grid(False)\nnuc_mt = (0.5*(nuc/255) + 0.5*(mt/255))\nplt.imshow(3*nuc_mt)\n\nplt.figure(figsize=(10,10))\nplt.title(\"Nucleus and ER\")\nplt.grid(False)\nnuc_er = (0.5*(nuc/255) + 0.5*(er/255))\nplt.imshow(3*nuc_er)\n\nplt.figure(figsize=(10,10))\nplt.title(\"Nucleus, microtubule and ER\")\nplt.grid(False)\nnuc_er_mt = (0.33*(nuc/255) + 0.33*(er/255) + 0.33*(mt/255))\nplt.imshow(3*nuc_er_mt)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2f3b4dcbacab92e79870b4489ca098a64e83e60f"},"cell_type":"markdown","source":"# Cell segmentation\n\n### Segmentation of single cells could increase predictive power of ML algorithms since the global distribution of proteins can vary. Looking at various images we saw that the amount and size of cells and therefore the POI shows a great deal of variance. The relative position of the POI in it's respective cell is of much more importance. Let us try some segmentation methods. Combining the nuclear, microtubule and ER staining gives us a comprehensive picture of the cell. So we will start by using a combination of those three channels for the segmentation. "},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"dca7a0a23763ba833d3460a2cdeac3b69a55dc96"},"cell_type":"code","source":"rnd = random.randint(0,len(os.listdir(PATH+\"train\")))\nrnd_img = os.listdir(PATH+\"train\")[rnd].split(\"_\")[0]\nnuc = cv2.imread(PATH+\"train/\"+rnd_img+\"_blue.png\", 0)\nmt = cv2.imread(PATH+\"train/\"+rnd_img+\"_yellow.png\", 0)\ner = cv2.imread(PATH+\"train/\"+rnd_img+\"_red.png\", 0)\ncomposit = cv2.add(nuc, mt, er)\ncomposit = cv2.resize(composit, (256,256))\n\nplt.figure()\nplt.grid(False)\nplt.title(\"original\")\nplt.imshow(composit, cmap=\"gray\")\n\n\n\n#test opencv threshold methods\nmethods = [\n(\"THRESH_BINARY\", cv2.THRESH_BINARY),\n(\"THRESH_BINARY_INV\", cv2.THRESH_BINARY_INV),\n(\"THRESH_TRUNC\", cv2.THRESH_TRUNC),\n(\"THRESH_TOZERO\", cv2.THRESH_TOZERO),\n(\"THRESH_TOZERO_INV\", cv2.THRESH_TOZERO_INV),\n(\"THRESH_OTSU\", cv2.THRESH_OTSU)]\n\nminval = 0\nmaxval = 255\n\nplt.figure(figsize=(10,10))\nfor i, (Name, Method) in enumerate(methods):\n    plt.subplot(3, 3, i+1)\n    plt.grid(False)\n    plt.title(Name)\n    (T, thresh) = cv2.threshold(composit, minval,maxval, type=Method)\n    plt.imshow(thresh, cmap=\"gray\")\n\nblur = cv2.GaussianBlur(composit,(5,5),0)\nplt.imshow(blur)\nt, thresh = cv2.threshold(blur, minval,maxval,cv2.THRESH_OTSU)\nplt.imshow(thresh)\nkernel = np.ones((5,5),np.uint8)\nclosing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel)\nplt.figure()\nplt.title(\"Thresholding with Closing\")\nplt.grid(False)\nplt.imshow(closing, cmap=\"gray\")\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8670ef7048f07940bc6a4bd623c7f7254d9cfea4"},"cell_type":"markdown","source":"# Using the complete composit image seemingly has the problem that the ER staining is leading to a very grainy thresholding of the image. Could try with stronger closing but probably could lead to problems when actually segmenting the cells. We will try another thing. Therefore we will threshold based on the nuclei, this should be much easier to do. Afterwards, we can measure the area of the nuclei and make a crude segmentation or at least image classification based on this. "},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"70830dbdbbccf6ca729beca82801f971375eeda4"},"cell_type":"code","source":"nuc = cv2.imread(PATH+\"train/\"+rnd_img+\"_blue.png\", 0)\nnuc = cv2.resize(nuc, (256,256))\nplt.figure(figsize=(15,5))\n\nplt.subplot(131)\nplt.grid(False)\nplt.title(\"Nuclear staining\")\nplt.imshow(nuc, cmap='gray')\n\nt, thresh = cv2.threshold(nuc, 0,255,cv2.THRESH_OTSU)\n\nplt.subplot(132)\nplt.grid(False)\nplt.title(\"OTSU thresholding\\nof nuclear staining\")\nplt.imshow(thresh)\n\nkernel = np.ones((4,4),np.uint8)\nclosing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel)\n\nplt.subplot(133)\nplt.grid(False)\nplt.title(\"OTSU thresholding\\nafter closing operation\")\nplt.imshow(closing)\n\nim, contours,hierarchy = cv2.findContours(closing, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n\ntest = cv2.cvtColor(nuc, cv2.COLOR_GRAY2BGR)\n\nt=cv2.drawContours(composit, contours, -1, (255,255,0), 1)\n\nplt.figure(figsize=(10,10))\nplt.grid(False)\nplt.title(\"Nuclear contours drawn in composit image\")\nplt.imshow(t)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ef25992381c916db1e90c71b15876a8e1d6a79a4"},"cell_type":"markdown","source":"### It is to note that depending on the closing kernel this would allow us to acutally not only identify the nucleus but also the nucleoli. It might be good to take a look at the protein of interest and if we could eliminate certain labels after a classification of images based on overlaps of POI and certain cellular compartments (e.g. if no overlap of POI and nucleus all the nuclear labels should be 0 in any case). "},{"metadata":{"trusted":true,"_uuid":"99039bd9146e756cc47f7ffa907049258ebe7be6"},"cell_type":"code","source":"rnd_imgs = []\nnoi = 50 \ntot_area = []\n\nfor i in range(noi):\n    rnd = random.randint(0,len(os.listdir(PATH+\"train\")))\n    rnd_imgs.append(os.listdir(PATH+\"train\")[rnd].split(\"_\")[0])\n    \nfor img_idx in rnd_imgs:\n    \n    nuc = cv2.imread(PATH+\"train/\"+img_idx+\"_blue.png\", 0)\n    nuc = cv2.resize(nuc, (256,256))\n    \n    t, thresh = cv2.threshold(nuc, 0,255,cv2.THRESH_OTSU)\n    \n    kernel = np.ones((4,4),np.uint8)\n    closing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel)\n    \n    im, contours,hierarchy = cv2.findContours(closing, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n\n    area = [cv2.contourArea(cnt) for cnt in contours]\n    tot_area.append(area)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2c057c15fe8779fcbc75b35f10b023b1007beffd"},"cell_type":"code","source":"check = [0,5,10,19,25,48]\nplt.figure()\nax = plt.subplot()\nmedian_areas = [np.median(area) for area in tot_area]\nplt.bar(range(len(median_areas)), (median_areas))\nfor i in check:\n    ax.patches[i].set_facecolor('r')\nax.spines[\"top\"].set_visible(False)\nax.spines[\"right\"].set_visible(False)\nplt.title(\"Median area of segmented nuclei of 50 random images\")\nplt.xlabel(\"# img\")\nplt.ylabel(\"Area (pxl)\")\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3d4c844d7dc886b0dd82dc32c6b35c14128447f9"},"cell_type":"markdown","source":"### Lets take a closer look look at some of the images in order to validate if this approach actually works"},{"metadata":{"trusted":true,"_uuid":"00d2222a073b8f9a10de8f92c5061cad87360dce"},"cell_type":"code","source":"check = [0,5,10,19,25,48]\n\nplt.figure(figsize=(10,10))\nfor i,idx in enumerate(check):\n    temp = cv2.imread(PATH+\"train/\"+rnd_imgs[idx]+\"_blue.png\", 0)\n    plt.subplot(2,3,i+1)\n    plt.grid(False)\n    plt.title(\"Image {}\".format(idx))\n    plt.imshow(temp,cmap='gray_r')","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}