<jupyter_start><jupyter_text>Machine Learning Engineer Nanodegree Capstone Project<NAME>  August 5th, 2019 I. Definition Project Overview[**CellSignal** - _Disentangling biological signal from experimental noise in cellular images_](https://www.rxrx.ai/) IntroductionMy final project is to participale in the NeurIPS competition on Kaggle called **CellSignal** - _Disentangling biological signal from experimental noise in cellular images_. More information about this competition is available here:  - Competition's website https://www.rxrx.ai  - Kaggle competition's link: https://www.kaggle.com/c/recursion-cellular-image-classification/overview Domain Background  Recursion Pharmaceuticals, creators of the industry’s largest dataset of biological images, generated entirely in-house, believes AI has the potential to dramatically improve and expedite the drug discovery process. More specifically, machine learning could help understand how drugs interact with human cells.This competition is designed to disentangle experimental noise from real biological signals. The goal is to classify images of cells under one of 1,108 different genetic perturbations (siRNA), and thus eliminate the noise introduced by technical execution and environmental variation between \[drug\] experiments.If successful, it could dramatically improve the industry’s ability to model cellular images according to their relevant biology. In turn, applying AI could greatly decrease the cost of treatments, and ensure these treatments get to patients faster. Problem StatementThis is a multiclass classification problem based on 6-channel fluorescent microscopy labeled images having a size 512x512 pixels.In plain english, based on 6-grayscale 512x512 images representing a single site of same-type cells, we must be able to detect the right genetic perturbation (siRNA) this site was altered with.Thus, the inputs will be 6x512x512 tensors and the output is a genetic perturbation (siRNA) represented by an integer ranging from 1 to 1,108.Below is an input example (which will be transformed in a tensor for training and predicting). Each channel represent different organelles of the cell - the nucleus, endoplasmic reticulum, actin cytoskeleton, nucleolus, mitochondria, and golgi apparatus.![cells in a well](https://raw.githubusercontent.com/michelml/ml-cellsignal/master/images/6channels.png)  And, next is an output example - an integer representing 1 of 1,108 genetic perturbations (siRNA) to classify. ```922``` Anticipated Solution for the problem at handThe solution for this problem will likely be resolved with the type of model architecture used in computer vision and image classification, e.g convolutional neural networks. This is a multiclass classification problem, but algorithms and model architectures we have seen in the dogs classification project https://github.com/MichelML/udacity-dog-project/, such as VGG-16 and ResNet-50, will be considered for this project. Metrics  We will use the cross entropy loss function and the accuracy metric to evaluate our model performance. [It is okay to stick with accuracy in our context](https://towardsdatascience.com/accuracy-paradox-897a69e2dd9b), as each class has a roughly equal number of data points.   Cross-entropy is the default loss function to use for multi-class classification problems ([source](https://machinelearningmastery.com/how-to-choose-loss-functions-when-training-deep-learning-neural-networks/)). Cross-entropy will calculate a score that summarizes the average difference between the actual and predicted probability distributions for all classes in the problem. The score is minimized and a perfect cross-entropy value is 0. II. Analysis<jupyter_code># all necessary imports for further sections
import numpy as np
import pandas as pd
import rxrx.io as rio
import geopandas as gpd
from shapely.geometry import Point
import matplotlib.pyplot as plt 

import torch
import torch.nn as nn
import torch.utils.data as D
from torch.optim.lr_scheduler import ExponentialLR
import torch.nn.functional as F
from torchvision import models, transforms

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Data Exploration    The data is available on the Kaggle's competition site https://www.kaggle.com/c/recursion-cellular-image-classification/data . For further information, see the [competition's website](https://rxrx.ai).![the dataset](https://raw.githubusercontent.com/michelml/ml-cellsignal/master/images/dataset_description.png)  One of the main challenges for applying AI to biological microscopy data is that even the most careful replicates of a process will not look identical. This dataset challenges you to develop a model for identifying replicates that is robust to experimental noise.The same siRNAs (effectively genetic perturbations) have been applied repeatedly to **four cell types** (HUVEC, HEPG2, RPE, U2OS) for a total of **51 experimental batches**. Each batch has **four plates**, each of which has **308 filled wells**. For each well, microscope images were taken at **two sites (per well)** and across **six imaging channels**.  Not every batch will necessarily have every well filled or every siRNA present, but those not filled were removed from the datasets.As explained earlier, a cell site is represented as 6 512x512 grayscale images. Each image (or channel) represents different organelles of the cell - the nucleus, endoplasmic reticulum, actin cytoskeleton, nucleolus, mitochondria, and golgi apparatus.  ```sh  identify images/example_input.png   --> example_input.png PNG 512x512 512x512+0+0 8-bit Gray 256c 58265B 0.000u 0:00.014```  Along with images and their labels, three other features are provided as metadata:  - `experiment`: the cell type and batch number  - `plate`: plate number within the experiment  - `well`: location on the plateThe outcome variable is siRNA, which is a multiclass variable which can be an integer ranging from 1 to 1,108, each integer representing a different genetic disruption.The data is already splitted into training and test sets. Further splitting will be performed for cross-validation purposes on the training set.     Controls  For each experiment, and thus each 1,108-siRNA set of wells, you have 30 positive controls and one or more negative control wells.  To understand the role of control wells, please refer to the beautiful explanation of DrHB in [this thread](https://www.kaggle.com/c/recursion-cellular-image-classification/discussion/101826latest-591591). Exploratory Visualization You can access further visualizations in the [full base-model notebook](https://github.com/MichelML/ml-cellsignal/blob/master/my_notebooks/basemodel.ipynb).  Note that the following visualization is greatly inspired by [Virginie Bonnefond's work](https://www.kaggle.com/pheaboo/a-journey-through-the-experiment-design).<jupyter_code># load metadata
md = rio.combine_metadata()
md.head()
# number of unique experiments
unique_experiments = md.experiment.unique()
len(unique_experiments)
# The same setup experiments have been tested on four cell types: HUVEC, RPE, HepG2, U2OS.
names, counts = np.unique([experiment.split('-')[0] for experiment in unique_experiments], return_counts=True)
for experiment_name, experiment_count in zip(names, counts):
    print('{} experiments focused on the cell type {}'.format(experiment_count, experiment_name))
one_experiment_md = md[(md.experiment=='RPE-01')]
one_experiment_md = one_experiment_md[(one_experiment_md.site == 1)]
one_experiment_md.groupby(['plate']).count()<jupyter_output><empty_output><jupyter_text>As discussed earlier, for each experiment, there are four plates in which 308 wells have been photographed (twice - two sites).<jupyter_code>def letter_to_int(letter):
    alphabet = list('abcdefghijklmnopqrstuvwxyz'.upper())
    return alphabet.index(letter) + 1

def well_to_point(well):
    letter = letter_to_int(well[0])
    number = int(well[1:])
    return Point(letter, number)

md['geometry'] = md.well.apply(lambda well: well_to_point(well))
md = gpd.GeoDataFrame(md)
def plot_well_type_positions_for_experiment(experiment_name):
    fig, axes = plt.subplots(1, 4, figsize=(20, 5))
    for plate, ax in zip([1,2,3,4], axes):
        one_plate_md = md[md.experiment==experiment_name]
        if plate==1:
            legend=True
        else:
            legend=False
        one_plate_md.plot(column='well_type',legend=legend, ax=ax);
        if plate==1:
            leg = ax.get_legend()
            leg.set_bbox_to_anchor((-0.3, 0., 0.2, 0.2))
    _ = fig.suptitle('Plate 1 to 4 - Experiment {}'.format(experiment_name))

plot_well_type_positions_for_experiment(md.experiment.unique()[0])
plot_well_type_positions_for_experiment(md.experiment.unique()[1])<jupyter_output><empty_output><jupyter_text>Each plate holds the same 30 control siRNA conditions (positive controlled - RED), 277 different non-control siRNA (treatment - TURQUOISE), and one untreated well (negative control - BLUE). The well types follow the same architecture from one plate to another and experiment to another.<jupyter_code>def plot_sirna_positions_for_experiment(experiment_name):
    fig, axes = plt.subplots(1, 4, figsize=(20, 5))
    for plate, ax in zip([1,2,3,4], axes):
        one_plate_md = md[md.experiment==experiment_name]
        one_plate_md = one_plate_md[one_plate_md.well_type=='positive_control'] 
        one_plate_md.plot(column='sirna', ax=ax, categorical=True, cmap='tab20c');
    _ = fig.suptitle('Plate 1 to 4 - Experiment {}'.format(experiment_name))
    
plot_sirna_positions_for_experiment(md.experiment.unique()[0])
plot_sirna_positions_for_experiment(md.experiment.unique()[-1])<jupyter_output><empty_output><jupyter_text>The location of each of the 1,108 non-control siRNA conditions is randomized in each experiment to prevent confounding effects of the location of a particular well.  Algorithms and Techniques  Our solution will use a convolutional neural network architecture, which is the defacto architecture for image classification ([source](https://medium.com/@ksusorokina/image-classification-with-convolutional-neural-networks-496815db12a8)).    To further justify our decision to go along with a CNN architecture, it is to be said ResNet-50 have achieved very good results in multiclass image classification problems in the recent past ([see source](https://arxiv.org/abs/1903.10035)), having reached 98.87% accuracy when classifying histopathology images.  We will thus use [transfer learning](https://en.wikipedia.org/wiki/Transfer_learning) with a pretrained ResNet-18 on ImageNet to see if we obtain interesting results, and will use complex versions of it (Resnet-50, Resnet-152 or Densenet(s)) if this is successful.   To use ResNet-18 properly, we will have to convert its first convolutional layer to support 6 channels instead of 3, and we will modify the final layer to support 1,108 output features (for the 1,108 possible siRNA classes). BenchmarkSince our solution uses a convolutional neural network architecture, we thus used a "vanilla" CNN model as our benchmark.  The base model (see [full base-model notebook](https://github.com/MichelML/ml-cellsignal/blob/master/my_notebooks/basemodel.ipynb)) was a very basic CNN model trained over 2 epochs on 6-channel 512x512 images. Here is the overall architecture:----------**Data**:6-channels images, both sites.Input resolution: 512x512**Base model**:Basic CNN model   ```  Net(  (pool): MaxPool2d(kernel_size=2, stride=2, padding=0, dilation=1, ceil_mode=False)  (conv1): Conv2d(6, 16, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)  (conv2): Conv2d(16, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)  (fc1): Linear(in_features=65536, out_features=2216, bias=False)  (fc2): Linear(in_features=2216, out_features=1108, bias=False))   ```   **Loss**:Cross Entropy Loss**Optimizer**:  Adam w/ lr=0.0003**Epochs**:2**Framework**:PyTorch----------This gave a cross-validation score of 0.0011 (.1% accuracy), and a test score of 0.002 (.2% accuracy). This score is a bit better than chance since we have 1108 classes. An accuracy reflecting chance would be 1/1108, which is equivalent to ~0.09% accuracy. III. MethodologyThe full notebook for the ResNet-50 implementation and training can be accessed [here](https://github.com/MichelML/ml-cellsignal/blob/master/my_notebooks/resnet50_firststep.ipynb). Data PreprocessingThe competition's authors already provided utils to load the 6-channel images for each cell site under the method `rio.load_site`. Here is an example where we convert a site into its 6-channel images.<jupyter_code>t = rio.load_site('train', 'RPE-05', 3, 'D19', 2, base_path="./data")
print(t.shape)
t_tensor_default = transforms.ToTensor()(t)
print(t_tensor_default.shape)<jupyter_output>(512, 512, 6)
torch.Size([6, 512, 512])<jupyter_text>This seems to work. Let's visualize individual channels.<jupyter_code>fig, axes = plt.subplots(2, 3, figsize=(24, 16))

for i, ax in enumerate(axes.flatten()):
  ax.axis('off')
  ax.set_title('channel {}'.format(i + 1))
  _ = ax.imshow(t[:, :, i], cmap='gray')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08765938069216758, 'fraction_numerical': 0.02663934426229508, 'mean_word_length': 5.332051898125901, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 18, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '565622', 'license': 'Apache-2.0', 'n_tokens_mistral': 3981, 'n_tokens_neox': 3540, 'n_words': 1426, 'path': 'my_notebooks/MachineLearningNanodegreeCapstoneReport.ipynb', 'repo_name': 'MichelML/ml-aging', 'size': '13183'}