Workflow+for+TIFF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import skimage.io\n",
    "import numpy\n",
    "import tensorflow as tf\n",
    "import pandas\n",
    "import pickle\n",
    "import os\n",
    "import os.path\n",
    "import math\n",
    "from itertools import groupby"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This function is to make sure all images stay at the same size:\n",
    "def __pad_or_crop(image, image_size):\n",
    "    bigger = max(image.shape[0], image.shape[1], image_size)\n",
    "\n",
    "    pad_x = float(bigger - image.shape[0])\n",
    "    pad_y = float(bigger - image.shape[1])\n",
    "\n",
    "    pad_width_x = (int(math.floor(pad_x / 2)), int(math.ceil(pad_x / 2)))\n",
    "    pad_width_y = (int(math.floor(pad_y / 2)), int(math.ceil(pad_y / 2)))\n",
    "    sample = image[int(image.shape[0]/2)-4:int(image.shape[0]/2)+4, :8]\n",
    "\n",
    "    std = numpy.std(sample)\n",
    "\n",
    "    mean = numpy.mean(sample)\n",
    "\n",
    "    def normal(vector, pad_width, iaxis, kwargs):\n",
    "        vector[:pad_width[0]] = numpy.random.normal(mean, std, vector[:pad_width[0]].shape)\n",
    "        vector[-pad_width[1]:] = numpy.random.normal(mean, std, vector[-pad_width[1]:].shape)\n",
    "        return vector\n",
    "\n",
    "    if (image_size > image.shape[0]) & (image_size > image.shape[1]):\n",
    "        return numpy.pad(image, (pad_width_x, pad_width_y), normal)\n",
    "    else:\n",
    "        if bigger > image.shape[1]:\n",
    "            temp_image = numpy.pad(image, (pad_width_y), normal)\n",
    "        else:\n",
    "            if bigger > image.shape[0]:\n",
    "                temp_image = numpy.pad(image, (pad_width_x), normal)\n",
    "            else:\n",
    "                temp_image = image\n",
    "        return temp_image[int((temp_image.shape[0] - image_size)/2):int((temp_image.shape[0] + image_size)/2),int((temp_image.shape[1] - image_size)/2):int((temp_image.shape[1] + image_size)/2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------- Labels of classes --------#\n",
    "classes = ['G1','S_phase','G2','Prophase','Metaphase','Anaphase','Telophase']\n",
    "\n",
    "#------- Specify here which channels to use-------#\n",
    "wanted_channel = ['Ch3','Ch4','Ch6']\n",
    "\n",
    "#------- Specify here the size of each image element -------#\n",
    "image_size = 40\n",
    "\n",
    "#------- Split data into training and validating set -------#\n",
    "split = 0.8\n",
    "\n",
    "#------- Parent input folder -------#\n",
    "path = '/home/minh-doan/Cell_cycle/Step2_input_single_tifs'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "G1_files = [os.path.join(path+'/G1/' ,file) for file in os.listdir(path+'/G1/')]\n",
    "S_files = [os.path.join(path+'/S_phase/' ,file) for file in os.listdir(path+'/S_phase/')]\n",
    "G2_files = [os.path.join(path+'/G2/' ,file) for file in os.listdir(path+'/G2/')]\n",
    "Prophase_files = [os.path.join(path+'/Prophase/' ,file) for file in os.listdir(path+'/Prophase/')]\n",
    "Metaphase_files = [os.path.join(path+'/Metaphase/' ,file) for file in os.listdir(path+'/Metaphase/')]\n",
    "Anaphase_files = [os.path.join(path+'/Anaphase/' ,file) for file in os.listdir(path+'/Anaphase/')]\n",
    "Telophase_files = [os.path.join(path+'/Telophase/' ,file) for file in os.listdir(path+'/Telophase/')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filelist = G1_files + S_files + G2_files + Prophase_files + Anaphase_files + Metaphase_files + Telophase_files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nested_filelist = [G1_files,S_files,G2_files,Prophase_files,Anaphase_files,Metaphase_files,Telophase_files]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Only care about those files that has the wanted channels\n",
    "wanted_filelist = [x for x in filelist if numpy.any([z in x for z in wanted_channel])]\n",
    "\n",
    "# Sort them into a good order. Channels might be specified near a group of letter. In my case it's '.ome'  \n",
    "wanted_filelist_channel_order = [list(j) for i, j in groupby(sorted(wanted_filelist), key=lambda x: x[x.index('.ome') - 1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you have a bunch of small TIFF files, each TIFF contains one object, use the following script to make tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Group up ALL the images into 1 single vector :\n",
    "images = [skimage.io.imread(wanted_filelist_channel_order[k][0]) for k in range(len(wanted_filelist_channel_order)) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Crop images and put each CROPPED image into its slot in the tensor T,[channel, image_index, width, height] with corresponding labels:\n",
    "\n",
    "multichannel_tensors = []\n",
    "labels = []\n",
    "\n",
    "for i in range(len(nested_filelist)): # each list in the nest is 1 class, and it contains multiple channels\n",
    "\n",
    "    single_channel_tensors = []\n",
    "    for j in range(len(wanted_channel)):\n",
    "\n",
    "        images_this_class_this_channel = [__pad_or_crop(image,image_size) for image in images[j::len(wanted_channel)][:int(len(nested_filelist[i])/len(wanted_channel))] ]       \n",
    "        \n",
    "        single_channel_tensors.append(numpy.array(images_this_class_this_channel)) # figure \n",
    "    \n",
    "    print('Transform nested single_channel_tensor of class ', i,' to a multi_channel_tensor')\n",
    "    multichannel_tensor = numpy.array(single_channel_tensors)\n",
    "    multichannel_tensors.append(multichannel_tensor) # pay attention to the plural \"tensorsssss\"\n",
    "    \n",
    "    label = numpy.zeros((multichannel_tensor.shape[1],len(nested_filelist)))\n",
    "    label[:multichannel_tensor.shape[0],i]=1\n",
    "    \n",
    "    labels.append(label)\n",
    "\n",
    "# Final tensor and label:  \n",
    "T = numpy.concatenate((multichannel_tensors), axis = 1)\n",
    "L = numpy.concatenate((labels))\n",
    "print('T.shape: ' + str(T.shape))\n",
    "print('L.shape: ' + str(L.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identified cellular objects could also be exported as one object per TIFF image. In this case, it is often inconvenient to store and transfer such multiple TIFF files. If user already stitched those small TIFF image together into a form of montage (grid), or user has a gallery of identified objects, please use this block of codes instead, which will transform a montage into corresponding tensors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Transform a montage of stitched images (grid) into a tensor (single channel)\n",
    "def im2tensor(image, cellWidth, cellHeight):\n",
    "    imgHeight,imgWidth = image.shape\n",
    "    rows = int(imgHeight/cellHeight)\n",
    "    cols = int(imgWidth/cellWidth)\n",
    "    tensor = image.reshape(rows,cellHeight,cols,cellWidth)\n",
    "    tensor = tensor.swapaxes(1,2)\n",
    "    tensor = tensor.reshape(-1, cellHeight,cellWidth)\n",
    "    return tensor[:,:,:,numpy.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make tensor and corresponding labels\n",
    "multichannel_tensors = []\n",
    "labels = []\n",
    "single_channel_tensors = []\n",
    "\n",
    "for i, j in groupby(sorted(wanted_filelist), key=lambda x: x[x.index('_CP') - 1]): # Channels specified at the position \"1 character from the left of _CP\" of filename. Only true if Channel < 10 sadly...   \n",
    "    list_temp = list(j)\n",
    "    print('channel: ' + i)\n",
    "       \n",
    "    cropped_images_this_channel = []\n",
    "    for k in range(len(list_temp)):\n",
    "\n",
    "        img = skimage.io.imread(list_temp[k])\n",
    "        images = im2tensor(img,image_size,image_size)\n",
    "        \n",
    "        # Filter empty tiles:\n",
    "        non_empty_tile = [image for image in images if numpy.std(image[:,int(image_size/2),:]) > 0]\n",
    "    \n",
    "        cropped_images_this_channel.append(numpy.array(non_empty_tile))       \n",
    "    \n",
    "    single_channel_tensors.append(numpy.concatenate((cropped_images_this_channel)))\n",
    "    \n",
    "    if len(single_channel_tensors) == len(wanted_channel): # it means we complete one folder      \n",
    "        print('Transform nested single_channel_tensor to multi_channel_tensor')\n",
    "        multichannel_tensor = numpy.concatenate((single_channel_tensors), axis = 3)\n",
    "        single_channel_tensors = []\n",
    "        multichannel_tensors.append(multichannel_tensor)\n",
    "        print('multichannel_tensor: ' + str(multichannel_tensor.shape))\n",
    "#         print('where am I? ' + list_temp[0])\n",
    "        \n",
    "        label = numpy.zeros((multichannel_tensor.shape[0],len(classes)))\n",
    "        for m in range(len(classes)):\n",
    "            if classes[m] in list_temp[k]:\n",
    "                label[:multichannel_tensor.shape[0],m]=1\n",
    "                labels.append(label)        \n",
    "\n",
    "# Final tensor and label:                  \n",
    "T = numpy.concatenate((multichannel_tensors))\n",
    "L = numpy.concatenate((labels))\n",
    "print('T.shape: ' + str(T.shape))\n",
    "print('L.shape: ' + str(L.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# filename = 'T' + '.sav'\n",
    "# pickle.dump(T, open(filename, 'wb'))\n",
    "# print(filename + \" saved !\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# filename = 'L' + '.sav'\n",
    "# pickle.dump(L, open(filename, 'wb'))\n",
    "# print(filename + \" saved !\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count how many cells in each classes, this block can be used to calculate class_weights"
   ]
  }, 
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Check for imbalance:\n",
    "\n",
    "from collections import Counter\n",
    "tuple_L = [tuple(element) for element in L]\n",
    "\n",
    "freq = Counter(tuple_L)\n",
    "\n",
    "import collections\n",
    "\n",
    "size = []\n",
    "for l in list(set(tuple_L)):\n",
    "    print( '%s : %d' % (l, freq[l]))\n",
    "    size.append(freq[l])"
   ]
  },  
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Split data set\n",
    "\n",
    "index = numpy.arange(L.shape[0])\n",
    "numpy.random.shuffle(index)\n",
    "\n",
    "Training_index = index[:int(L.shape[0]*split)]\n",
    "Training_set_temp = T[Training_index,:,:,:]\n",
    "Training_labels_set_temp = L[Training_index,:]\n",
    "Training_set_temp_index = numpy.arange(Training_set_temp.shape[0])\n",
    "\n",
    "Testing_index = index[:-int(L.shape[0]*split)]\n",
    "Testing_set = T[Testing_index,:,:,:]\n",
    "Testing_labels_set = L[Testing_index,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Placeholders for tensorflow\n",
    "\n",
    "x = tf.placeholder(tf.float32, shape=[None, image_size,image_size,len(wanted_channel)]) # Last number is number of color channels\n",
    "y_ = tf.placeholder(tf.float32, shape=[None, len(classes)]) # Last number is number of classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BASIC DEEP LEARNING\n",
    "\n",
    "This block of code define the network. This is the standard LeNet model described in the Google Tensorflow tutorials. However, there is one important difference between this model and the stock model provided by Google: Here we compute a final fully connected layer 'h_fc2' separate from the softmax layer. This is what we will use to get the high-level representations the network is learning and visualize it, as learned from [this tutorial](https://medium.com/@awjuliani/visualizing-deep-learning-with-t-sne-tutorial-and-video-e7c59ee4080c#.k75ll068e) "
   ]
  },  
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Kernels, weights, biases, settings for Deep learning\n",
    "\n",
    "def weight_variable(shape):\n",
    "    initial = tf.truncated_normal(shape, stddev=0.1)\n",
    "    return tf.Variable(initial)\n",
    "\n",
    "def bias_variable(shape):\n",
    "    initial = tf.constant(0.1, shape=shape)\n",
    "    return tf.Variable(initial)\n",
    "\n",
    "def conv2d(x, W):\n",
    "    return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME')\n",
    "\n",
    "def max_pool_2x2(x):\n",
    "    return tf.nn.max_pool(x, ksize=[1, 2, 2, 1],\n",
    "                        strides=[1, 2, 2, 1], padding='SAME')\n",
    "\n",
    "# shape of [5, 5, len(wanted_channel), 32]. the third is the number of input channels (colors)\n",
    "W_conv1 = weight_variable([5, 5, len(wanted_channel), 32])\n",
    "b_conv1 = bias_variable([32])\n",
    "\n",
    "# First convolutional layer\n",
    "h_conv1 = tf.nn.relu(conv2d(x, W_conv1) + b_conv1)\n",
    "h_pool1 = max_pool_2x2(h_conv1)\n",
    "\n",
    "# Second convolutional layer\n",
    "W_conv2 = weight_variable([5, 5, 32, 64])\n",
    "b_conv2 = bias_variable([64])\n",
    "\n",
    "h_conv2 = tf.nn.relu(conv2d(h_pool1, W_conv2) + b_conv2)\n",
    "h_pool2 = max_pool_2x2(h_conv2)\n",
    "\n",
    "# Densely Connected Layer\n",
    "W_fc1 = weight_variable([int(image_size/4)*int(image_size/4)*64, 1024]) # Max_pool is 2x2, do it twice, so have to divide by 4\n",
    "b_fc1 = bias_variable([1024])\n",
    "h_pool2_flat = tf.reshape(h_pool2, [-1, int(image_size/4)*int(image_size/4)*64])\n",
    "h_fc1 = tf.nn.relu(tf.matmul(h_pool2_flat, W_fc1) + b_fc1)\n",
    "\n",
    "# Dropout\n",
    "keep_prob = tf.placeholder(tf.float32)\n",
    "h_fc1_drop = tf.nn.dropout(h_fc1, keep_prob)\n",
    "\n",
    "W_fc2 = weight_variable([1024, len(classes)]) # Last number is number of classes\n",
    "b_fc2 = bias_variable([len(classes)]) # Number is number of classes\n",
    "\n",
    "h_fc2 = tf.matmul(h_fc1_drop, W_fc2) + b_fc2\n",
    "y_conv=tf.nn.softmax(h_fc2)\n",
    "\n",
    "cross_entropy = -tf.reduce_sum(y_*tf.log(y_conv))\n",
    "\n",
    "# y_conv = tf.matmul(h_fc1_drop, W_fc2) + b_fc2\n",
    "\n",
    "# cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=y_conv, labels=y_))\n",
    "\n",
    "correct_prediction = tf.equal(tf.argmax(y_conv,1), tf.argmax(y_,1))\n",
    "accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "epoch = int(Training_labels_set_temp.shape[0]/1000)*1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Configure which GPU/CPU device to use tensorflow\n",
    "configuration = tf.ConfigProto()\n",
    "configuration.gpu_options.allow_growth = True\n",
    "configuration.gpu_options.visible_device_list = \"2\"\n",
    "sess = tf.InteractiveSession(config = configuration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_step = tf.train.AdamOptimizer(1e-4).minimize(cross_entropy)\n",
    "\n",
    "# sess = tf.InteractiveSession()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "finalRepresentations = []\n",
    "\n",
    "for i in range(epoch):\n",
    "\n",
    "    numpy.random.shuffle(Training_set_temp_index)\n",
    "    batchIndex = Training_set_temp_index[0:30]\n",
    "\n",
    "    Training_set = Training_set_temp[batchIndex,:,:,:]\n",
    "    Training_labels_set = Training_labels_set_temp[batchIndex,:]\n",
    "\n",
    "    if i % 1000 == 0:\n",
    "        train_accuracy = accuracy.eval(feed_dict={x:Training_set, y_: Training_labels_set, keep_prob: 1.0})\n",
    "        print(\"step %d, training accuracy %g\"%(i, train_accuracy))\n",   
    "        finalRepresentations.append(h_fc2.eval(session=sess, feed_dict={x:Testing_set, keep_prob:1.0}))\n",
    "        # current representation of the network. Represent how the network \"see\" the data.\n",    
    "    train_step.run(feed_dict={x: Training_set, y_: Training_labels_set, keep_prob: 0.5})\n",
    "\n",
    "print(\"test accuracy %g\"%accuracy.eval(session=sess,feed_dict={x: Testing_set, y_: Testing_labels_set, keep_prob: 1.0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SAVE the wanted stage of Deep Learning filter\n",
    "\n",
    "filename = 'Last_DL_filter' + '.sav'\n",
    "pickle.dump(h_fc2.eval(session=sess, feed_dict={x:Testing_set, keep_prob:1.0}), open(filename, 'wb'))\n",
    "print(filename + \" saved !\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Save metadata for Tensorboard\n",
    "\n",
    "metadata = os.path.join('metadata.tsv')\n",
    "\n",
    "images_tb = tf.Variable(finalRepresentations[-1], name='images') # SPECIFY HERE WHICH DEEP LEARNING FILTER YOU WOULD LIKE TO USE\n",
    "            \n",
    "def save_metadata(file):\n",
    "    with open(file, 'w') as f:\n",
    "        for i in range(Testing_labels_set.shape[0]):\n",
    "            c = numpy.nonzero(Testing_labels_set[::1])[1:][0][i]\n",
    "            f.write('{}\\n'.format(c))\n",
    "            \n",
    "save_metadata('metadata.tsv')\n",
    "\n",
    "\n",
    "with tf.Session() as sess:\n",
    "    saver = tf.train.Saver([images_tb])\n",
    "\n",
    "    sess.run(images_tb.initializer)\n",
    "    saver.save(sess, os.path.join('images_tb.ckpt'))\n",
    "\n",
    "    config = projector.ProjectorConfig()\n",
    "    # One can add multiple embeddings.\n",
    "    embedding = config.embeddings.add()\n",
    "    embedding.tensor_name = images_tb.name\n",
    "    # Link this tensor to its metadata file (e.g. labels).\n",
    "    embedding.metadata_path = metadata\n",
    "    # Saves a config file that TensorBoard will read during startup.\n",
    "    projector.visualize_embeddings(tf.summary.FileWriter('.'), config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Close Tensorflow session, save memory !\n",
    "sess.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check carefully this output : \"projector_config.pbtxt\", you have to specify absolute path \"/path/to/logdir/metadata.tsv\", CANNOT be relative path \"./metadata.tsv\", nor \"~/metadata.tsv\". Next, use Terminal for this command: tensorboard --logdir=\"/path/to/logdir\". Finally, open web-browser, connect to http://localhost:6006 to visualize the results. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}