{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Let's see what is going on here!"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.\nfrom pathlib import Path\nfrom PIL import Image\nfrom tqdm import tqdm_notebook\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nBASE_DIR = Path('../input/')\nCHANNELS = [1, 2, 3, 4, 5, 6]\nSITES = [1, 2]\nRGB_MAP = {\n    1: {\n        'rgb': np.array([19, 0, 249]),\n        'range': [0, 51]\n    },\n    2: {\n        'rgb': np.array([42, 255, 31]),\n        'range': [0, 107]\n    },\n    3: {\n        'rgb': np.array([255, 0, 25]),\n        'range': [0, 64]\n    },\n    4: {\n        'rgb': np.array([45, 255, 252]),\n        'range': [0, 191]\n    },\n    5: {\n        'rgb': np.array([250, 0, 253]),\n        'range': [0, 89]\n    },\n    6: {\n        'rgb': np.array([254, 255, 40]),\n        'range': [0, 191]\n    }\n}\n\ntrain_meta = pd.read_csv('../input/train.csv')\ntest_meta = pd.read_csv('../input/test.csv')\ntrain_meta['well_type'] = 'treatment'\ntest_meta['sirna'] = None\ntest_meta['well_type'] = 'treatment'\ntrain_controls = pd.read_csv('../input/train_controls.csv')\ntest_controls = pd.read_csv('../input/test_controls.csv')\ncommon_controls = [1108, 1109, 1115, 1116, 1117, 1121, 1123, 1124, 1125, 1126, 1128, 1129, 1131, 1134, 1135, 1136, 1137, 1138]\ntrain_meta = pd.concat([train_meta, train_controls])\ntest_meta = pd.concat([test_meta, test_controls])\n\nprint(train_meta.shape, test_meta.shape, train_controls.shape, test_controls.shape)\n\nNUM_WELL_ROWS = 24\nNUM_WELL_COLS = 16","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Define functions for visualization."},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_image_tensor(experiment, plate, well, site, channels=CHANNELS):\n    if os.path.exists(BASE_DIR / 'train' / experiment):\n        # It is from train.\n        d = BASE_DIR / 'train' / experiment / f'Plate{plate}'\n    else:\n        # It is from test.\n        d = BASE_DIR / 'test' / experiment / f'Plate{plate}'\n\n    img_channels = []\n    for channel in channels:\n        path = d / f'{well}_s{site}_w{channel}.png'\n        img_channels.append(np.array(Image.open(path)))\n    \n    return np.array(img_channels).transpose(1, 2, 0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Copied from\n# https://github.com/recursionpharma/rxrx1-utils/blob/master/rxrx/io.py\ndef convert_tensor_to_rgb(t, channels=CHANNELS, vmax=255, rgb_map=RGB_MAP):\n    \"\"\"\n    Converts and returns the image data as RGB image\n    Parameters\n    ----------\n    t : np.ndarray\n        original image data\n    channels : list of int\n        channels to include\n    vmax : int\n        the max value used for scaling\n    rgb_map : dict\n        the color mapping for each channel\n        See rxrx.io.RGB_MAP to see what the defaults are.\n    Returns\n    -------\n    np.ndarray the image data of the site as RGB channels\n    \"\"\"\n    colored_channels = []\n    for i, channel in enumerate(channels):\n        x = (t[:, :, i] / vmax) / \\\n            ((rgb_map[channel]['range'][1] - rgb_map[channel]['range'][0]) / 255) + \\\n            rgb_map[channel]['range'][0] / 255\n        x = np.where(x > 1., 1., x)\n        x_rgb = np.array(\n            np.outer(x, rgb_map[channel]['rgb']).reshape(512, 512, 3),\n            dtype=int)\n        colored_channels.append(x_rgb)\n    im = np.array(np.array(colored_channels).sum(axis=0), dtype=int)\n    im = np.where(im > 255, 255, im)\n    return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def visualize(experiment, plate, well, site, channels=CHANNELS, ax=None):\n    img_tensor = load_image_tensor(experiment, plate, well, site, channels)\n    img = convert_tensor_to_rgb(img_tensor, channels=channels)\n    \n    if ax is None:\n        fig = plt.figure()\n        ax = fig.add_subplot(111)\n    ax.imshow(img)\n    ax.set_axis_off()\n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"well_width = 24\n\nfig = plt.figure(figsize=(well_width*1, well_width*6))\nchannels = []\nfor i in range(1, 7):\n    ax = fig.add_subplot(1, 6, i)\n    channels.append(i)\n    ax.set_title('Channels: ' + ' '.join(map(str, channels)))\n    visualize('HUVEC-01', 1, 'B02', 1, channels=channels, ax=ax);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It works!"},{"metadata":{},"cell_type":"markdown","source":"# Visualization of pseudo-plate.\n\nHere we reconstruct and visualize plates in-silico. For that, we have to map a (0-based) coordinate of a well in a 384-well (16 x 24) plate to a well ID.\n\nAcross well IDs, there are 14 kinds of alphabets (B~O) and 22 kinds of numbers (02~23). With that information, we can make a lookup table!"},{"metadata":{"trusted":true},"cell_type":"code","source":"well2coord = dict()\nfor row, num in enumerate(range(1, 25)):  # We include outer rows and columns which were excluded due to the environmental effect here.\n    for col, alphabet in enumerate('ABCDEFGHIJKLMNOP'):\n        is_valid = (1 <= row <= 22) and (1 <= col <= 14)\n        well2coord[f'{alphabet}{num:02}'] = (row, col, is_valid)\n\n# Hold inverse mapping as well.\ncoord2well = {(r, c):w for w, (r, c, v) in well2coord.items()}\n\ndef is_valid_well(well):\n    return well2coord[well][2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_sirna(experiment, plate, well):\n    mask = (train_meta.experiment == experiment) & (train_meta.plate == plate) & (train_meta.well == well)\n    return train_meta[mask].sirna.values[0], train_meta[mask].well_type.values[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_pseudoplate(experiment, plate, site=1, channels=CHANNELS):\n    # Figure aesthetics.\n    well_size = 10\n    \n    # Plotting.\n    fig = plt.figure(figsize=(NUM_WELL_COLS * well_size, NUM_WELL_ROWS * well_size))\n\n    for row in tqdm_notebook(range(NUM_WELL_ROWS), desc='row'):\n        for col in range(NUM_WELL_COLS):\n            well = coord2well[(row, col)]\n            if not is_valid_well(well):\n                continue\n            ax = fig.add_subplot(NUM_WELL_ROWS, NUM_WELL_COLS, row * NUM_WELL_COLS + (col + 1))\n            sirna, well_type = get_sirna(experiment, plate, well)\n            ax.set_title(f'{sirna} ({well_type})')\n            ax.set_axis_off()\n            visualize(experiment, plate, well, site=site, channels=channels, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"show_pseudoplate('HUVEC-01', 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualizing all the controls.\n\nAre there any apparent commonalities of cells treated with control siRNAs?\n\nFirst, let's see how the controls are distributed throughout the plate."},{"metadata":{"trusted":true},"cell_type":"code","source":"def show_controls(experiment, plate, ax):\n    controls = train_meta[(train_meta.experiment == experiment) & (train_meta.plate == plate) & (train_meta.well_type.str.contains('control'))]\n\n    for i, row in controls.iterrows():\n        y, x, _ = well2coord[row.well]\n        color = 'red' if row.well_type == 'positive_control' else 'blue'\n        ax.scatter(x, NUM_WELL_ROWS - y, c=color)\n    ax.set_xlim([0, NUM_WELL_COLS])\n    ax.set_ylim([0, NUM_WELL_ROWS])\n    ax.set_aspect('equal')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 30))\n\nnum_experiments = 7\nfor experiment in range(1, num_experiments + 1):\n    exp = f'HEPG2-{experiment:02}'\n    for plate in range(1, 5):\n        ax = fig.add_subplot(num_experiments, 4, (experiment - 1) * 4 + plate)\n        ax.set_title(f'{exp} Plate{plate}')\n        show_controls(experiment=exp, plate=plate, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 60))\n\nnum_experiments = 16\nfor experiment in range(1, num_experiments + 1):\n    exp = f'HUVEC-{experiment:02}'\n    for plate in range(1, 5):\n        ax = fig.add_subplot(num_experiments, 4, (experiment - 1) * 4 + plate)\n        ax.set_title(f'{exp} Plate{plate}')\n        show_controls(experiment=exp, plate=plate, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 30))\n\nnum_experiments = 7\nfor experiment in range(1, num_experiments + 1):\n    exp = f'RPE-{experiment:02}'\n    for plate in range(1, 5):\n        ax = fig.add_subplot(num_experiments, 4, (experiment - 1) * 4 + plate)\n        ax.set_title(f'{exp} Plate{plate}')\n        show_controls(experiment=exp, plate=plate, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12, 13))\n\nnum_experiments = 3\nfor experiment in range(1, num_experiments + 1):\n    exp = f'U2OS-{experiment:02}'\n    for plate in range(1, 5):\n        ax = fig.add_subplot(num_experiments, 4, (experiment - 1) * 4 + plate)\n        ax.set_title(f'{exp} Plate{plate}')\n        show_controls(experiment=exp, plate=plate, ax=ax)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some plates have got additional negative controls (such as RPE-06 Plate 2), missing positive controls (such as HUVEC-13 Plate 4). But generally all the controls are seemingly at identical positions."},{"metadata":{},"cell_type":"markdown","source":"## Visualizing positive control siRNA-treated cell images.\n\nLet's compare images of cells that are treated with positive control siRNAs with nontreated cell lines (negative controls)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_well_with_sirna(experiment, plate, sirna):\n    mask = (train_meta.experiment == experiment) & (train_meta.plate == plate) & (train_meta.sirna == sirna)\n    return train_meta[mask].well.values[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def visualize_cells_with_sirna(experiment, sirna, channels=CHANNELS):\n    well_size = 10\n    title_fs = 20\n\n    fig = plt.figure(figsize=(4*well_size, 2*well_size))\n    experiment = experiment\n    negative_sirna = 1108\n    positive_sirna = sirna\n\n    for plate in range(1, 5):\n        ax1 = fig.add_subplot(2, 4, plate)\n        ax2 = fig.add_subplot(2, 4, 4 + plate)\n\n        neg_well = get_well_with_sirna(experiment, plate, negative_sirna)\n        pos_well = get_well_with_sirna(experiment, plate, positive_sirna)\n\n        visualize(experiment, plate, neg_well, site=1, channels=channels, ax=ax1)\n        visualize(experiment, plate, pos_well, site=1, channels=channels, ax=ax2)\n        ax1.set_title(f'{experiment} Plate{plate}, negative control', fontdict={'fontsize': title_fs})\n        ax2.set_title(f'{experiment} Plate{plate}, siRNA {positive_sirna}', fontdict={'fontsize': title_fs})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As an example, lets visualize cells treated with siRNA 1109 along with untreated cells."},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[4])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[6])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-01', 1109, channels=[1, 2, 3, 4, 5, 6])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-02', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HUVEC-03', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HEPG2-01', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HEPG2-02', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('HEPG2-03', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('RPE-01', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('RPE-02', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('RPE-03', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('U2OS-01', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('U2OS-02', 1109)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"visualize_cells_with_sirna('U2OS-03', 1109)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}