{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Downloading HPA public data"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\n# Install W&B\n!pip install wandb --upgrade\n# Install HPA Cell Segmentation tool.\n!pip install https://github.com/CellProfiling/HPA-Cell-Segmentation/archive/master.zip","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# General imports.\nimport io\nimport re\nimport os\nimport cv2\nimport glob\nimport gzip\nimport imageio\nimport pathlib\nimport requests\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom PIL import Image\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nfrom skimage.transform import resize\nfrom sklearn.model_selection import train_test_split\n\n%matplotlib inline\n\n# HPA Segmentation tool related imports\nimport hpacellseg.cellsegmentator as cellsegmentator\nfrom hpacellseg.utils import label_cell, label_nuclei\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n%matplotlib inline\n\n\ndef tif_gzip_to_png(tif_path):\n    '''Function to convert .tif.gz to .png and put it in the same folder\n    Eg. for working in local work station\n    '''\n    png_path = pathlib.Path(tif_path.replace('.tif.gz','.png'))\n    tf = gzip.open(tif_path).read()\n    img = imageio.imread(tf, 'tiff')\n    imageio.imwrite(png_path, img)\n    \ndef download_and_convert_tifgzip_to_png(url, target_path):    \n    '''Function to convert .tif.gz to .png and put it in the same folder\n    Eg. in Kaggle notebook\n    '''\n    r = requests.get(url)\n    f = io.BytesIO(r.content)\n    tf = gzip.open(f).read()\n    img = imageio.imread(tf, 'tiff')\n    imageio.imwrite(target_path, img)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# W&B import and login\nimport wandb\nfrom kaggle_secrets import UserSecretsClient\n\nuser_secrets = UserSecretsClient()\nwandb_api = user_secrets.get_secret(\"wandb_api\")\n\nwandb.login(key=wandb_api)\n\nLOG_WANDB = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# All label names in the public HPA and their corresponding index. \nall_locations = dict({\n    \"Nucleoplasm\": 0,\n    \"Nuclear membrane\": 1,\n    \"Nucleoli\": 2,\n    \"Nucleoli fibrillar center\": 3,\n    \"Nuclear speckles\": 4,\n    \"Nuclear bodies\": 5,\n    \"Endoplasmic reticulum\": 6,\n    \"Golgi apparatus\": 7,\n    \"Intermediate filaments\": 8,\n    \"Actin filaments\": 9,\n    \"Focal adhesion sites\": 9,\n    \"Microtubules\": 10,\n    \"Mitotic spindle\": 11,\n    \"Centrosome\": 12,\n    \"Centriolar satellite\": 12,\n    \"Plasma membrane\": 13,\n    \"Cell Junctions\": 13,\n    \"Mitochondria\": 14,\n    \"Aggresome\": 15,\n    \"Cytosol\": 16,\n    \"Vesicles\": 17,\n    \"Peroxisomes\": 17,\n    \"Endosomes\": 17,\n    \"Lysosomes\": 17,\n    \"Lipid droplets\": 17,\n    \"Cytoplasmic bodies\": 17,\n    \"No staining\": 18\n})\n\n\ndef add_label_idx(df, all_locations):\n    '''Function to convert label name to index\n    '''\n    df[\"Label_idx\"] = None\n    for i, row in df.iterrows():\n        labels = row.Label.split(',')\n        idx = []\n        for l in labels:\n            if l in all_locations.keys():\n                idx.append(str(all_locations[l]))\n        if len(idx)>0:\n            df.loc[i,\"Label_idx\"] = \"|\".join(idx)\n            \n        print(df.loc[i,\"Label\"], df.loc[i,\"Label_idx\"])\n    return df\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"public_hpa_df = pd.read_csv('../input/publichpa-withcellline/kaggle_2021.tsv')\n# Remove all images overlapping with Training set\npublic_hpa_df = public_hpa_df[public_hpa_df.in_trainset == False]\n\n# Remove all images with only labels that are not in this competition\npublic_hpa_df = public_hpa_df[~public_hpa_df.Label_idx.isna()]\n\ncolors = ['blue', 'red', 'green', 'yellow']\ncelllines = ['A-431', 'A549', 'EFO-21', 'HAP1', 'HEK 293', 'HUVEC TERT2', 'HaCaT', 'HeLa', 'PC-3', 'RH-30', 'RPTEC TERT1', 'SH-SY5Y', 'SK-MEL-30', 'SiHa', 'U-2 OS', 'U-251 MG', 'hTCEpi']\npublic_hpa_df_17 = public_hpa_df[public_hpa_df.Cellline.isin(celllines)]\nlen(public_hpa_df), len(public_hpa_df_17)\n\npublic_hpa_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if LOG_WANDB:\n    public_hpa_df.to_csv('public_hpa.csv', index=False)\n\n    # log as artifact\n    run = wandb.init(project='hpa', job_type='hpa_public_dataset_creation')\n    artifact = wandb.Artifact('hpa_public_data', type='dataset')\n    artifact.add_file('public_hpa.csv')\n    run.log_artifact(artifact)\n    run.join()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_labels = public_hpa_df.Label_idx.unique().tolist()\nall_labels = '|'.join(all_labels)\nall_labels = all_labels.split('|')\nall_labels = list(set(all_labels))\nnum_unique_labels = len(all_labels)\nall_labels = sorted(all_labels, key=int)\nall_labels = ' '.join(all_labels)\nprint(f'{num_unique_labels} unique labels, values: {all_labels}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labels = [str(i) for i in range(19)]\n\nunique_counts = {}\nfor lbl in labels:\n    unique_counts[lbl] = len(public_hpa_df[public_hpa_df.Label_idx == lbl])\n\nfull_counts = {}\nfor lbl in labels:\n    count = 0\n    for row_label in public_hpa_df['Label_idx']:\n        if lbl in row_label.split('|'): count += 1\n    full_counts[lbl] = count\n    \ncounts = list(zip(full_counts.keys(), full_counts.values(), unique_counts.values()))\ncounts = np.array(sorted(counts, key=lambda x:-x[1]))\ncounts = pd.DataFrame(counts, columns=['label', 'full_count', 'unique_count'])\n\nsns.set(style=\"whitegrid\")\nf, ax = plt.subplots(figsize=(16, 12))\n\nsns.set_color_codes(\"pastel\")\nsns.barplot(x=\"full_count\", y=\"label\", data=counts, order=counts.label.values,\n            label=\"full count\", color=\"b\", orient = 'h')\n\n# Plot the crashes where alcohol was involved\nsns.set_color_codes(\"muted\")\nsns.barplot(x=\"unique_count\", y=\"label\", data=counts, order=counts.label.values,\n            label=\"unique count\", color=\"b\", orient = 'h')\n\n# Add a legend and informative axis label\nax.legend(ncol=2, loc=\"lower right\", frameon=True)\nax.set(ylabel=\"\",\n       xlabel=\"Counts\")\nsns.despine(left=True, bottom=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Get Single Labeled Images"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Add a column - num_classes\npublic_hpa_df['num_classes'] = public_hpa_df['Label_idx'].apply(lambda r: len(r.split('|')))\nprint(f'Total number of images: {len(public_hpa_df)}')\npublic_hpa_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# interested in single labels only\nprint(public_hpa_df['num_classes'].value_counts())\npublic_hpa_df['num_classes'].value_counts().plot.bar(title='Examples with multiple labels', xlabel='number of labels per example', ylabel='# train examples')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_one_label = public_hpa_df.loc[public_hpa_df['num_classes'] == 1]\nprint(f'Number of images with one image-level labels: {len(df_one_label)}')\ndf_one_label.Label_idx = df_one_label.Label_idx.astype('int64')\ndf_one_label.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Ref: https://www.kaggle.com/divyanshuusingh/eda-image-segmentation\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}\n\nlabels, counts = np.unique(df_one_label.Label_idx.values, return_counts=True)\nprint(f'The unique labels are: {labels} and there values are: {counts}')\n\nplt.figure(figsize=(15,5))\nplt.bar(labels, counts)\n\nfor index, value in enumerate(counts):\n    plt.text(index-0.25, value, str(value), fontdict=dict(fontsize=10))\n\nplt.xticks(np.arange(len(labels)), labels=label_names.values(), rotation=85)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"! Warning: These are raw, full size images. Each channel is approximately 8MB. With ~~82495~~ 77668 images * 4 channels, everything amounts to around ~~2.6TB~~ 2.4TB of data.\n\nWe have 17 cell lines in the training set and test set. So if you want to download extra public data for training and downloading is too slow because of big size, you probably want to consider:\n- Downloading just 17 cell lines (67k images * 4 channels, instead of 77.6k images)\n- Sampling according to label (eg. You have a lot of Nucleoplasm and Cytosol in the training set already so maybe you just need more rare labels)\n- Using jpeg images. These were purposely created for visualization on HPA website. They are much smaller in size (you will loose some info compared to tif file, but maybe that's enough for your model. Your call). To download jpeg, simply change `.tif.gz` to `.jpg` in the url. For example: https://images.proteinatlas.org/10005/921_B9_1_blue.jpg\t"},{"metadata":{"trusted":true},"cell_type":"code","source":"DESIRED_LABEL = 1\nLABEL_NAME = '_'.join(label_names[DESIRED_LABEL].split(' ')).lower()\n\ntmp_df = df_one_label.loc[df_one_label.Label_idx==DESIRED_LABEL]\ntmp_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if LOG_WANDB:\n    tmp_df.to_csv(f'pubic_{LABEL_NAME}.csv')\n    \n    # log as artifact\n    run = wandb.init(entity='ayush-thakur', project='hpa', job_type=f'public_{LABEL_NAME}')\n    artifact_csv = run.use_artifact('ayush-thakur/hpa/hpa_public_data:v0', type='dataset')\n\n    artifact = wandb.Artifact(f'{LABEL_NAME}', type='dataset')\n    artifact.add_file(f'pubic_{LABEL_NAME}.csv')\n    run.log_artifact(artifact)\n    run.join()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"SAVE_DIR = f'/kaggle/tmp/hpa_public_cell_level_{LABEL_NAME}_128x128/'\n\nos.makedirs(SAVE_DIR+'rgb', exist_ok=True)\nos.makedirs(SAVE_DIR+'protein', exist_ok=True)\n\n!ls {SAVE_DIR}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"NUC_MODEL = \"./nuclei-model.pth\"\nCELL_MODEL = \"./cell-model.pth\"\n\nsegmentator = cellsegmentator.CellSegmentator(\n    NUC_MODEL,\n    CELL_MODEL,\n    scale_factor=0.25,\n    device=\"cuda\",\n    padding=True,\n    multi_channel_model=True,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"IMAGE_HEIGHT = 128\nIMAGE_WIDTH = 128\n\n# check if the contour for nuclie is touching the image boundary. \ndef is_border_nuclei(contour_points):\n    unique_points = np.unique(contour_points)\n    # basically if any point is 0 that means its touching the edge of the image.\n    if 0 in unique_points:\n        return True\n    return False\n\ndef extract_cell_level_images(paths, image_id):\n    '''\n    Input\n    \n    paths - absolute path of colored images. \n    '''    \n    # LOAD IMAGES\n    for path in paths:\n        if 'red' in path:\n            mt = path\n        elif 'blue' in path:\n            nu = path\n        elif 'yellow' in path:\n            er = path\n        elif 'green' in path:\n            pr = path\n    \n    # Get all channel images.\n    microtubule = np.array(Image.open(mt))\n    endoplasmicrec = np.array(Image.open(er))\n    nuclei = np.array(Image.open(nu))\n    protein = np.array(Image.open(pr))\n    protein = ((protein/np.max(protein))*255).astype('uint8')\n\n    # Stack the channels to form image. Using Red, Yellow and Blue.\n    image = np.dstack((microtubule, endoplasmicrec, nuclei))\n    image = ((image/np.max(image))*255).astype('uint8')\n\n    # PERFORM SEGMENTATION\n    \n    # For nuclei segmentation only blue channel is required.\n    nuc_segmentation = segmentator.pred_nuclei([nu])\n    # For full cells all the three reference(except green) channels are required.\n    cell_segmentation = segmentator.pred_cells([[mt], [er], [nu]])\n    # get cell mask\n    nuclei_mask, cell_mask = label_cell(nuc_segmentation[0], cell_segmentation[0])\n\n    # GET IMDIVIDUAL CELLS\n    \n    # Count the number of cells.\n    cells = np.unique(cell_mask)\n    nuclei = np.unique(nuclei_mask)\n\n    cell_count = 0\n    for cell_index in cells[1:]:\n        # Get cell and nucleus mask for one cell.\n        single_cell_mask = np.where(cell_mask==cell_index, 1,0).astype('uint8')\n        if cell_index in nuclei:\n            nucleus_mask = np.where(nuclei_mask==cell_index, 1,0).astype('uint8')\n        else:\n            continue\n        \n        # get contour for cell and nucleus\n        cell_cnts, _ = cv2.findContours(single_cell_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n        nucleus_cnts, _ = cv2.findContours(nucleus_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n\n        # Check if the nucleus is touching the boundary of image.\n        if not is_border_nuclei(nucleus_cnts[0]): # If not touching the boundary \n            # Mask the cell to be cropped.\n            cell = cv2.bitwise_and(image, image, mask=single_cell_mask)\n            # Mask the protein to be cropped.\n            cell_protein = cv2.bitwise_and(protein, protein, mask=nucleus_mask)\n            \n            # Get bounding box covering the contour. \n            x,y,w,h = cv2.boundingRect(cell_cnts[0])\n            # Crop cell\n            cell = cell[y:y+h, x:x+w]\n            # Resize cell\n            cell = cv2.resize(cell, (IMAGE_HEIGHT, IMAGE_WIDTH), cv2.INTER_AREA)\n            \n            # Crop protein\n            cell_protein = cell_protein[y:y+h, x:x+w]\n            # Resize protein\n            cell_protein = cv2.resize(cell_protein, (IMAGE_HEIGHT, IMAGE_WIDTH), cv2.INTER_AREA)\n            \n            # Save images\n            cv2.imwrite(SAVE_DIR+'rgb/'+image_id+f'_{cell_count}.png', cell)\n            cv2.imwrite(SAVE_DIR+'protein/'+image_id+f'_{cell_count}.png', cell_protein)\n            \n            cell_count+=1\n            \n    print('.', end='')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"save_dir = os.path.join(os.getcwd(),'publichpa1')\nif not os.path.exists(save_dir):\n    os.makedirs(save_dir)\n    \nfor i, row in tqdm(tmp_df.iterrows()):\n    try:\n        # Download Image\n        img = row.Image\n        color_image_paths = []\n        for color in colors:\n            img_url = f'{img}_{color}.tif.gz'\n            save_path = os.path.join(save_dir,  f'{os.path.basename(img)}_{color}.png')\n            download_and_convert_tifgzip_to_png(img_url, save_path)            \n            color_image_paths.append(save_path)\n    except:\n        print(f'failed to download: {img}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Number of images: {len(os.listdir(save_dir))//4}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"images = os.listdir(save_dir)\nfor i, row in tqdm(tmp_df.iterrows()):\n    img = row.Image.split('/')[-1]\n    channel_imgs = []\n    for image in images:\n        if img in image:\n            channel_imgs.append(save_dir+'/'+image)\n    \n    extract_cell_level_images(channel_imgs, img)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Copy Kaggle API token to ~/.kaggle\n! mkdir -p /root/.kaggle/\n! cp ../input/apitoken/kaggle.json /root/.kaggle/kaggle.json\n# Initialize dataset creation\n! kaggle datasets init -p {SAVE_DIR}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls {SAVE_DIR}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%bash\necho \"{\n  \\\"title\\\": \\\"HPA: Single Label Public Nuclear Membrane\\\",\n  \\\"id\\\": \\\"ayuraj/SingleLabelPublicNuclearMembrane\\\",\n  \\\"licenses\\\": [\n    {\n      \\\"name\\\": \\\"CC0-1.0\\\"\n    }\n  ]\n}\" > /kaggle/tmp/hpa_public_cell_level_nuclear_membrane_128x128/dataset-metadata.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !kaggle datasets create -p {SAVE_DIR} -u --dir-mode tar\n!kaggle datasets version -p {SAVE_DIR} -m \"add complete images\"  --dir-mode tar","execution_count":null,"outputs":[]},{"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}