{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Exploratory data analysis - Human Protein Atlas competition \nExploratory data analysis (EDA) for the Human Protein Atlas (HPA) competition.\n\nThis notebooks is based on:\n\n1. [Human Protein Atlas Cell Segmentation + [EDA] by Kuldeep Singh Chouhan](https://www.kaggle.com/kool777/human-protein-atlas-cell-segmentation-eda)\n\nAdditional Kaggle datasets for using HPA CellSegmentator without downloading files:\n1. [Nuclei and cell weights for HPA CellSegmentator by RDizzl3](https://www.kaggle.com/rdizzl3/hpacellsegmentatormodelweights).\n2. [Nuclei and cell for HPA CellSegmentator by daishuby](https://www.kaggle.com/daishu/hpacellsegmodel).  \n.","metadata":{}},{"cell_type":"code","source":"# libraries\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport glob\nimport warnings \nimport os,gc,cv2\nimport shutil\nfrom tqdm.notebook import tqdm\n\n%matplotlib inline\nwarnings.filterwarnings('ignore')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# HPA single-cell image segmentation \n!pip install https://github.com/CellProfiling/HPA-Cell-Segmentation/archive/master.zip\nimport hpacellseg.cellsegmentator as cellsegmentator\nfrom hpacellseg.utils import label_cell, label_nuclei","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions","metadata":{}},{"cell_type":"code","source":"# functions\n\n\n# read and visualize sample image\ndef read_sample_image(filename, path):    \n    '''\n    read individual images\n    of different filters (R, G, B, Y)\n    and stack them.\n    ---------------------------------\n    Arguments:\n    filename -- sample image path\n    \n    Returns:\n    stacked_images -- stacked (RGBY) image\n    '''\n    red = cv2.imread( path +'train/' + filename + \"_red.png\", cv2.IMREAD_UNCHANGED)\n    green = cv2.imread( path + 'train/' + filename + \"_green.png\", cv2.IMREAD_UNCHANGED)\n    blue = cv2.imread( path + 'train/' + filename + \"_blue.png\", cv2.IMREAD_UNCHANGED)\n    yellow = cv2.imread( path + 'train/' + filename + \"_yellow.png\", cv2.IMREAD_UNCHANGED)\n    stacked_images = np.transpose( np.array([red, green, blue, yellow] ), (1,2,0) )\n    return stacked_images\n\n\ndef plot_all(im, label):    \n    '''\n    plot all RGBY image,\n    Red, Green, Blue, Yellow, \n    filters images.\n    --------------------------\n    Argument:\n    im - image\n    ''' \n    plt.figure(figsize=(15, 15))\n    plt.subplot(1, 5, 1)\n    plt.imshow(im[:,:,:3])\n    plt.title('RGBY Image')\n    plt.axis('off')\n    plt.subplot(1, 5, 2)\n    plt.imshow(im[:,:,0], cmap='Reds')\n    plt.title('Microtubule channels')\n    plt.axis('off')\n    plt.subplot(1, 5, 3)\n    plt.imshow(im[:,:,1], cmap='Greens')\n    plt.title('Protein of Interest')\n    plt.axis('off')\n    plt.subplot(1, 5, 4)\n    plt.imshow(im[:,:,2], cmap='Blues')\n    plt.title('Nucleus')\n    plt.axis('off')\n    plt.subplot(1, 5, 5)\n    plt.imshow(im[:,:,3], cmap='Oranges')\n    plt.title('Endoplasmic Reticulum')\n    plt.axis('off')\n    plt.show()\n\n    \n# read and visualize sample image\ndef read_sample_image_seg(filename, path):    \n    '''\n    read individual images\n    of different filters (R, B, Y)\n    and stack them for segmentation.\n    ---------------------------------\n    Arguments:\n    filename -- sample image file path\n    \n    Returns:\n    stacked_images -- stacked (RBY) image path in lists.\n    '''\n    red = path + 'train/' + filename + \"_red.png\"\n    blue = path + 'train/' + filename + \"_blue.png\"\n    yellow = path + 'train/' + filename + \"_yellow.png\"\n    stacked_images = [[red], [yellow], [blue]]\n    return stacked_images, red, blue, yellow\n\n\n# segment cell \ndef segmentCell(image, segmentator):   \n    '''\n    segment cell and nuclei from\n    microtubules, endoplasmic reticulum,\n    and nuclei (R, B, Y) filters.\n    ------------------------------------\n    Argument:\n    image -- (R, B, Y) list of image arrays\n    segmentator -- CellSegmentator class object\n    \n    Returns:\n    cell_mask -- segmented cell mask\n    '''\n    nuc_segmentations = segmentator.pred_nuclei(image[2])\n    cell_segmentations = segmentator.pred_cells(image)\n    nuclei_mask, cell_mask = label_cell(nuc_segmentations[0], cell_segmentations[0])\n    gc.collect(); del nuc_segmentations; del cell_segmentations; del nuclei_mask\n    return cell_mask\n\n\n# plot segmented cells mask, image\ndef plot_cell_segments(mask, red, blue, yellow):    \n    '''\n    plot segmented cells\n    and images\n    ---------------------\n    Arguments:\n    mask -- cell mask\n    red -- red filter image path\n    blue -- blue filter image path\n    yellow -- yellow filter image path\n    '''\n    microtubule = plt.imread(r)    \n    endoplasmicrec = plt.imread(b)    \n    nuclei = plt.imread(y)\n    img = np.dstack((microtubule, endoplasmicrec, nuclei))\n    plt.figure(figsize=(15, 15))\n    plt.subplot(1, 3, 1)\n    plt.imshow(img)\n    plt.title('Image')\n    plt.axis('off')\n    plt.subplot(1, 3, 2)\n    plt.imshow(mask)\n    plt.title('Mask')\n    plt.axis('off')\n    plt.subplot(1, 3, 3)\n    plt.imshow(img)\n    plt.imshow(mask, alpha=0.6)\n    plt.title('Image + Mask')\n    plt.axis('off')\n    plt.show()\n\n    \n# plot single segmented cells mask, image\ndef plot_single_cell(mask, red, blue, yellow):\n    '''\n    plot single cell mask\n    and image\n    ---------------------\n    Arguments:\n    mask -- cell mask\n    red -- red filter image path\n    blue -- blue filter image path\n    yellow -- yellow filter image path\n    '''\n    microtubule = plt.imread(r)    \n    endoplasmicrec = plt.imread(b)    \n    nuclei = plt.imread(y)\n    img = np.dstack((microtubule, endoplasmicrec, nuclei))\n    \n    contours= cv2.findContours(mask.astype('uint8'),\n                               cv2.RETR_TREE, \n                               cv2.CHAIN_APPROX_SIMPLE)\n    areas = [cv2.contourArea(c) for c in contours[0]]\n    x = np.argsort(areas)\n    cnt = contours[0][x[-1]]\n    x,yc,w,h = cv2.boundingRect(cnt)  \n    plt.figure(figsize=(15, 15))\n    plt.subplot(1, 3, 1)\n    plt.imshow(img[yc:yc+h, x:x+w])\n    plt.title('Cell Image')\n    plt.axis('off')\n    plt.subplot(1, 3, 2)\n    plt.imshow(mask[yc:yc+h, x:x+w])\n    plt.title('Cell Mask')\n    plt.axis('off') \n    plt.subplot(1, 3, 3)\n    plt.imshow(img[yc:yc+h, x:x+w])\n    plt.imshow(mask[yc:yc+h, x:x+w], alpha=0.6)\n    plt.title('Cell Image + Mask')\n    plt.axis('off')\n    plt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. EDA for labels and images\n","metadata":{}},{"cell_type":"markdown","source":"### 1.1  EDA for labels\nCount label distribution across samples and number of labels per sample.","metadata":{}},{"cell_type":"code","source":"# data directory\nDIR = '../input/hpa-single-cell-image-classification/'\nos.listdir(DIR)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# csv files\ntrain_df = pd.read_csv( DIR + 'train.csv' )\nsample_submission = pd.read_csv( DIR + 'sample_submission.csv' )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot class counts\n\n# spliting label column\ntrain_df[\"Label\"] = train_df[\"Label\"].str.split(\"|\")\n\n# class labels\nclass_labels = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', '18']\n\n# binarizing each label/class\nfor label in tqdm(class_labels):\n    train_df[label] = train_df['Label'].map(lambda result: 1 if label in result else 0)\n\n# rename column\ntrain_df.columns = ['ID', 'Label', 'Nucleoplasm', 'Nuclear membrane', 'Nucleoli', 'Nucleoli fibrillar center',\n                    'Nuclear speckles', 'Nuclear bodies', 'Endoplasmic reticulum', 'Golgi apparatus', 'Intermediate filaments',\n                    'Actin filaments', 'Microtubules', 'Mitotic spindle', 'Centrosome', 'Plasma membrane', 'Mitochondria',\n                    'Aggresome', 'Cytosol', 'Vesicles and punctate cytosolic patterns', 'Negative']\n\n# class counts dataframe \nclass_counts = train_df.sum().drop(['ID', 'Label']).sort_values(ascending=False)\nprint(class_counts, '\\n')\n\n# generate plot\nplt.figure( figsize=(14,12) )\nax=sns.barplot(y=class_counts.index.values, x=class_counts.values, palette='tab10')\nplt.suptitle(\"Label Distribution\")\nplt.title(\"HPA single-cell classification Kaggle competition 2021\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* The most frequent label is Nucleoplasm, and there is clearly a class imbalance.\n* There are only 34 negative examples, which is the less frequent class.\n* It might be useful to try data augmentation methods in order to generate a more balanced set.","metadata":{}},{"cell_type":"code","source":"# plot labels per sample/image\nlabel_per_image = train_df.drop(['ID', 'Label'], axis=1).sum(axis=1)\nplt.figure(figsize=(16,10))\nax = sns.countplot(label_per_image, palette='Pastel1')\nfor p in ax.patches:\n    height = p.get_height()\n    ax.text( p.get_x()+p.get_width()/2.,\n             height + 3,\n             '{:1.2f}%'.format(height/len(label_per_image)*100),\n             ha=\"center\", fontsize=12)\n    plt.title(\"Label Per Sample/Image\", fontsize=16)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* A large percentage of images (48.19%) contain only one label.\n* The majority of images (88.88%) contain one or two labels.\n* This implies that only 11.12% of the images contains three or more labels.","metadata":{}},{"cell_type":"markdown","source":"### 1.2 EDA for images\nIn this subsection we want to understand how each sample image consists of four channels for each different cellular structures: Red - microtubules, Blue - nuclei, Yellow - endoplasmatic reticulum (ER), and Green - protein of interest. The classification task is to predict \"protein organelle localization labels for each cell in the image\" [see Kaggle notebook: 'Single-cell Patterns'](https://www.kaggle.com/lnhtrang/single-cell-patterns). This is why each image has multiple labels corresponding to these organelles. Only images with single labels (largest number in the training set) are considered. This means that each plot corresponds to one organelle label  and shows the stacked image (RGBY) along with each individual channel. For example the first image has the label 'Nucleoplasm' and means that the protein of interest is located only at this organelle according to the annotation.","metadata":{}},{"cell_type":"code","source":"# select images with single labels\ntrain = train_df.loc[train_df['Label'].apply(lambda x: len(x)==1)==True]\n\n# plot one example image per organelle\nfor label in train_df.drop(['ID', 'Label'], axis=1):\n    print(label)\n    im = read_sample_image(train[train[label]==1].sample(1).ID.to_string().split(' ')[4],  path=DIR)\n    plot_all(im, label)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1.3 Label-wise image cell segmentation \nEach image should be segmented into cells. The [HPA segmentation tool](https://github.com/CellProfiling/HPA-Cell-Segmentation) can be used to perform this segmentation task. Here we will show how to do it. Each plot corresponds to single label RBY image, the segmented mask, and the stacked image + mask. There is an issue when trying to download the HPA CellSegmentator model weights, reported [here](https://www.kaggle.com/c/hpa-single-cell-image-classification/discussion/230090) and [here](https://www.kaggle.com/c/hpa-single-cell-image-classification/discussion/231036). These weights are available as a [dataset](https://www.kaggle.com/rdizzl3/hpacellsegmentatormodelweights) (by RDizzl3) or in this [other dataset](https://www.kaggle.com/daishu/hpacellsegmodel) (by daishuby).  ","metadata":{}},{"cell_type":"code","source":"# Define CellSegmentator class\n# [source: https://github.com/CellProfiling/HPA-Cell-Segmentation]\n#-----------------------------------------------------------\n\n# [1] path to the nuclei model weights (from dataset):\nNUC_MODEL = '../input/hpacellsegmentatormodelweights/dpn_unet_nuclei_v1.pth'\n\n#-----------------------------------------------------------\n# [2] path to the cell model weights (from dataset):\nCELL_MODEL = '../input/hpacellsegmentatormodelweights/dpn_unet_cell_3ch_v1.pth'\n\n#-----------------------------------------------------------\n# [3] scale_factor: determines how much the images should be \n# scaled before being fed to the models. For HPA Cell images, \n# a value of 0.25 (default) is good.\n#-----------------------------------------------------------\n# [4] device: Inform Torch which device to put the model on. \n#Valid  values are ‘cpu’ or ‘cuda’ or pointed cuda device  \n# like 'cuda:0’. Defaults to cuda.\n#-----------------------------------------------------------\n# [5] padding: If True, add some padding before feeding the \n# images to the neural networks. This is not required but \n# can make segmentations, especially cell segmentations,\n# more accurate. Defaults to False. Note: If you have issues \n# running the segmentation due to image dimensions, setting \n# padding to True may help.\n#-----------------------------------------------------------\n# [6] multi_channel_model: If True, use the pretrained \n# three-channel version of the model. Having this set to \n# True gives you better cell segmentations but requires \n# you to give the model endoplasmic reticulum images as \n# part of the cell segmentation. Otherwise, the version \n# trained with only two channels, microtubules and nuclei, \n# will be used. Defaults to True\n#-----------------------------------------------------------\nsegmentator = cellsegmentator.CellSegmentator(\n    NUC_MODEL,\n    CELL_MODEL,\n    scale_factor=0.25,\n    device=\"cpu\",\n    padding=True,\n    multi_channel_model=True,\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# select images with single labels\nfor label in train_df.drop(['ID', 'Label'], axis=1):\n    print(label)\n    # plot one example image per organelle\n    im, r, b, y = read_sample_image_seg(train[train[label]==1].sample(1).ID.to_string().split(' ')[4], path=DIR)\n    mask = segmentCell(im, segmentator)\n    plot_cell_segments(mask, r, b, y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 1.4 Label-wise  single-cell segmentation","metadata":{}},{"cell_type":"code","source":"for label in train_df.drop(['ID', 'Label'], axis=1):\n    print(label)\n    im, r, b, y = read_sample_image_seg(train[train[label]==1].sample(1).ID.to_string().split(' ')[4],path=DIR)\n    mask = segmentCell(im, segmentator)\n    plot_single_cell(mask, r, b, y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}