Deep learning¶
Our toy example will be to classify species as endangered or not based genomic data. The rationale is that species with a small (effective population size) will have higher chances to be threatened. The amount of genomic variability (e.g. polymorphic sites and haplotype diversity) is taken as a proxy for the (effective) population size of each species. The genomic variation at marker loci is represented as an image.
Objective¶
We assume that we have collected some genomic data on many species which have already been categorised into 4 classes of conservation status: least concern, vulnerable, endangered, critically endagered. Our goal is to implement a classifier that given genomic data can predic whether that species is endangered or not.
Now the know how to process and manipulate images with python, we can actually do some science! As explained before our goal is to build a classifier to predict whether a certain species is endangered or not. We use genetic information as proxy for the ability of the species to react to novel conditions.
We assume we have 4 classes of conservation status (LC, VU, EN, CR) and 400 samples per class. Each data point is an image which represents the (biallelic) genetic variation (so it’s binary) across individuals (on the rows) over several genetic loci (on the columns). Images are double sorted to remove some noise associated to the order of samples and polymorphic sites.
In [1]:
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy
Training and testing data¶
The first think we do is to load images representing genomic data from various species. These images have already been converted into numpy arrays. (Note that we have already preprocessed our data, so pixel intensitities have been scaled between 0 and 1.)
In [ ]:
X_train = np.load(“Data/X_train.npy”)
print(X_train.shape)
If we want to plot one image, then we need to do a proper slicing.
In [ ]:
print(X_train[0,:,:,0].shape)
plt.imshow(X_train[0,:,:,0], cmap=plt.cm.binary)
which corresponds to the following matrix:
In [ ]:
X_train[0,:,:,0]
Apart from the data itself (called X), we also need labels (called y) associated to each data point. These labels indicate which class each data point belongs to. What is the rank and shape of this numpy array?
In [ ]:
y_train = np.load(“Data/y_train.npy”)
print(y_train.shape)
y_train
If you recall, we collected 400 samples per class so we should expect 1600 data points. We have 1200 entries in the training set. Why?
We split the data we have into training and testing. The learning will be done in the training set and the measurement of accuracy will happen on the testing set. What is the rank and shape of the these numpy objects?
In [ ]:
X_test = np.load(“Data/X_test.npy”)
print(X_test.shape)
y_test = np.load(“Data/y_test.npy”)
print(y_test.shape)
Finally, we want to predict the classification of an unknown species, the Marsican bear (labelled as ursus). Let’s load its image, check its rank and plot it.
In [ ]:
X_ursus = np.load(“Data/X_ursus.npy”)
print(X_ursus.shape)
plt.imshow(X_ursus[0,:,:,0], cmap=plt.cm.binary)
Nearest Neighbour Classifier¶
Now we wan to implement a NN classifier to assign a label to our image of the Marsican bear. The idea is to select the image(s) which is (are) closer to the one of interest. For instance, the simple elementwise difference between images can be calculated as:
In [ ]:
X_ursus[0,:,:,0] – X_train[0,:,:,0]
DIY (5) Implement a NN classifier using the L1 distance and predict the label for the Ursus image. A possible framework is given below but feel free to use your own creativity. (A possible solution is in the Solutions folder, learning.py file.)
In [ ]:
labels = [“LC”, “VU”, “EN”, “CR”] # these are the 4 possible classes
min_distance = 9999 # initialise a value for the distance
for i in range(0, X_train.shape[0]):
distance = …
if …
…
…
Now you can easily implement a NN classifer applied to the whole testing set. What’s the achieved accuracy? You can do it yourself, but the idea is that you have to execute the above operation for each sample from the testing set.
TensorFlow¶
We will be using TensorFlow to build and train our deep neural network. In particular, we will use Keras which is TensorFlow’s high-level API. Keras is used for fast prototyping and production since it is
• User friendly: a simple, consistent interface optimized for common use cases.
• Modular and composable: models are made by connecting building blocks together.
• Easy to extend: you can write custom building blocks, new layers, loss functions, etc. etc.
tf.keras is the implementation of Keras API specification. tf.keras can run any Keras-compatible code.
In [ ]:
import tensorflow as tf
Sequential model¶
Let’s build a simple model. In Keras, you assemble layers to build models. A model is a graph of layers. The most common type of model is a stack of layers: the tf.keras.Sequential model. For instance, we can build a simple, fully-connected network with the following code.
In [ ]:
model = tf.keras.Sequential()
# Adds a densely-connected layer with 16 units to the model:
model.add(tf.keras.layers.Dense(16, activation=’relu’))
# Add another:
model.add(tf.keras.layers.Dense(16, activation=’relu’))
# Add a softmax layer with 10 output units:
model.add(tf.keras.layers.Dense(5, activation=’softmax’))
Configure the layers¶
There are several tf.keras.layers available with some common constructor parameters:
• activation: set the activation function for the layer, specified by the name of a built-in function or as a callable object.
• kernel_initializer and bias_initializer: the initialization schemes that create the layer’s weights (kernel and bias).
• kernel_regularizer and bias_regularizer: the regularization schemes that apply the layer’s weights (kernel and bias). The following examples show how to create instances of tf.keras.layers.Dense layers using constructor arguments.
In [ ]:
# Create a sigmoid layer:
tf.keras.layers.Dense(16, activation=’sigmoid’)
# Or:
tf.keras.layers.Dense(16, activation=tf.sigmoid)
# A linear layer with L1 regularization of factor 0.01 applied to the kernel matrix:
tf.keras.layers.Dense(16, kernel_regularizer=tf.keras.regularizers.l1(0.01))
# A linear layer with L2 regularization of factor 0.01 applied to the bias vector:
tf.keras.layers.Dense(16, bias_regularizer=tf.keras.regularizers.l2(0.01))
# A linear layer with a kernel initialized to a random orthogonal matrix:
tf.keras.layers.Dense(16, kernel_initializer=’orthogonal’)
# A linear layer with a bias vector initialized to 2.0s:
tf.keras.layers.Dense(16, bias_initializer=tf.keras.initializers.constant(2.0))
Training¶
After the model is constructed, we can configure its learning process by calling the compile method.
In [ ]:
model = tf.keras.Sequential([
# Adds a densely-connected layer with 16 units to the model:
tf.keras.layers.Dense(16, activation=’relu’),
# Add another:
tf.keras.layers.Dense(16, activation=’relu’),
# Add a softmax layer with 5 output units:
tf.keras.layers.Dense(5, activation=’softmax’)])
model.compile(optimizer=tf.train.AdamOptimizer(0.001),
loss=’categorical_crossentropy’,
metrics=[‘accuracy’])
tf.keras.Model.compile takes three main arguments:
• optimizer: it specifies the training procedure, such as tf.train.AdamOptimizer, tf.train.RMSPropOptimizer, or tf.train.GradientDescentOptimizer.
• loss: function to minimize during optimization, such as mean square error (mse), categorical_crossentropy, and binary_crossentropy, specified by name or by passing a callable object from the tf.keras.losses module.
• metrics: to monitor training; string names or callables from the tf.keras.metrics module.
Let’s make a quick example using numpy to train and evaluate our model using the fit method.
In [ ]:
data = np.random.random((1000, 32))
labels = np.random.random((1000, 5))
model.fit(data, labels, epochs=10, batch_size=32)
tf.keras.Model.fit takes three main arguments:
• epochs: training is structured into epochs, iterations over the entire input data.
• batch_size: the model slices the data into smaller batches and iterates over these batches during training.
• validation_data: tuple of inputs and labels for validation.
When prototyping a model, you want to easily monitor its performance on some validation data. Passing the validation_data argument allows the model to display the loss and metrics in inference mode for the passed data, at the end of each epoch. Here’s an example using validation_data.
In [ ]:
data = np.random.random((1000, 32))
labels = np.random.random((1000, 5))
val_data = np.random.random((100, 32))
val_labels = np.random.random((100, 5))
model.fit(data, labels, epochs=10, batch_size=32,
validation_data=(val_data, val_labels))
For large datasets or multi-device training you can use the Datasets API.
Evaluate and predict¶
The tf.keras.Model.evaluate and tf.keras.Model.predict methods can use NumPy data and a tf.data.Dataset to evaluate the inference-mode loss and metrics.
In [ ]:
data = np.random.random((1000, 32))
labels = np.random.random((1000, 5))
model.evaluate(data, labels, batch_size=32)
We can predict the output for the data provided.
In [ ]:
result = model.predict(data, batch_size=32)
print(result.shape)
You can read the full documentation for building custom models and more advanced models. You may also want to check out TensorBoard and Low Level information which is helpful for debugging. A series of tutorials is also available.
DNN for predicting conservation status¶
Let’s build our CNN. First thing, we need to define the architecture. To do that, we need to import some modules from keras. Instead of using tensorflow, we can use the high-levekl API keras with a tensforlow backend.
In [16]:
import numpy as np
import matplotlib.pyplot as plt
import keras
import tensorflow as tf
X_train = np.load(“Data/X_train.npy”)
y_train = np.load(“Data/y_train.npy”)
X_test = np.load(“Data/X_test.npy”)
y_test = np.load(“Data/y_test.npy”)
X_ursus = np.load(“Data/X_ursus.npy”)
Architecture¶
Building the neural network requires configuring the layers of the model, then compiling the model. The sequential model type from keras is simply a linear stack of neural network layers. Each layer will be added (or stacked) to the initial model.
In [17]:
model = keras.Sequential()
The model needs also to know what input shape it should expect. The first layer (but only the first) needs to receive such information. Let’s add a convolutional layer. We can read how to do it from here. Also note that we can add activation and padding layers directly into this convolutional layer.
In [18]:
model.add(keras.layers.Conv2D(filters=16, kernel_size=(3,3), strides=(1,1), activation=”relu”, padding=”same”, input_shape=(64, 64, 1)))
After a convolutional layer it is often used a max-pooling layer. Read its definition here.
In [19]:
model.add(keras.layers.MaxPooling2D(pool_size=(2,2)))
To prevent overfitting, one trick is to use a Dropout layer.
In [20]:
model.add(keras.layers.Dropout(rate=0.5))
You can add several cycles of Conv-MaxPool-Dropout. If you then want to move towards the final fully connected neural network, you need to first flatten the network.
In [21]:
model.add(keras.layers.Flatten())
Finally we can add a fully connected network with a relu activation using a Dense layer, another dropout layer, and then the output layer with a softmax activation.
In [22]:
model.add(keras.layers.Dense(units=32, activation=”relu”))
# we need 4 units at the output since we have 4 classes
model.add(keras.layers.Dense(units=4, activation=”softmax”))
We can even print a summary of our architecture with the specification of the learnable parameters.
In [23]:
model.summary()
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
conv2d_2 (Conv2D) (None, 64, 64, 16) 160
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 32, 32, 16) 0
_________________________________________________________________
dropout_2 (Dropout) (None, 32, 32, 16) 0
_________________________________________________________________
flatten_2 (Flatten) (None, 16384) 0
_________________________________________________________________
dense_3 (Dense) (None, 32) 524320
_________________________________________________________________
dense_4 (Dense) (None, 4) 132
=================================================================
Total params: 524,612
Trainable params: 524,612
Non-trainable params: 0
_________________________________________________________________
Compiling¶
We need to “compile” the network and prepare it for the training. We do it by specifying the loss function and optimisation to use.
In [25]:
model.compile(optimizer=tf.train.AdamOptimizer(0.001),
loss=’categorical_crossentropy’,
metrics=[‘accuracy’])
Training¶
Let’s train our network. We pass the training data set and the network will optimise its parameters to minimise the loss function. We can allocate a portion of the training data as validation data. You can read more here.
In [26]:
hist = model.fit(X_train, y_train, batch_size=32, epochs=10, validation_split=0.20)
Train on 960 samples, validate on 240 samples
Epoch 1/10
960/960 [==============================] – 1s 858us/step – loss: 1.4510 – acc: 0.3219 – val_loss: 1.2742 – val_acc: 0.4833
Epoch 2/10
960/960 [==============================] – 0s 507us/step – loss: 1.1317 – acc: 0.5708 – val_loss: 1.0486 – val_acc: 0.6625
Epoch 3/10
960/960 [==============================] – 0s 512us/step – loss: 0.8701 – acc: 0.7448 – val_loss: 0.8598 – val_acc: 0.6750
Epoch 4/10
960/960 [==============================] – 0s 507us/step – loss: 0.6433 – acc: 0.8281 – val_loss: 0.7386 – val_acc: 0.6917
Epoch 5/10
960/960 [==============================] – 0s 510us/step – loss: 0.5189 – acc: 0.8510 – val_loss: 0.7071 – val_acc: 0.7375
Epoch 6/10
960/960 [==============================] – 0s 516us/step – loss: 0.4436 – acc: 0.8844 – val_loss: 0.5961 – val_acc: 0.7792
Epoch 7/10
960/960 [==============================] – 0s 518us/step – loss: 0.3714 – acc: 0.9240 – val_loss: 0.5457 – val_acc: 0.7917
Epoch 8/10
960/960 [==============================] – 0s 516us/step – loss: 0.3217 – acc: 0.9260 – val_loss: 0.5973 – val_acc: 0.7625
Epoch 9/10
960/960 [==============================] – 1s 529us/step – loss: 0.2730 – acc: 0.9469 – val_loss: 0.5566 – val_acc: 0.7667
Epoch 10/10
960/960 [==============================] – 1s 526us/step – loss: 0.2350 – acc: 0.9594 – val_loss: 0.5092 – val_acc: 0.7667
The validation accuracy is used to tune the hyper-parameters (e.g. learning rate, dropout rate). It’s convenient to plot the decay of loss and increase of accuracy for both the training and validation set.
In [27]:
from matplotlib import rcParams
train_loss = hist.history[‘loss’]
val_loss = hist.history[‘val_loss’]
train_acc = hist.history[‘acc’]
val_acc = hist.history[‘val_acc’]
xc = range(10)
x_axis = np.zeros(len(xc))
for x,i in enumerate (xc):
x_axis[i] = x + 1
rcParams[‘axes.titlepad’] = 20
plt.figure(1,figsize=(7,5),facecolor=’white’)
plt.plot(x_axis,train_loss)
plt.plot(x_axis,val_loss)
plt.xlabel(‘Epoch’, fontsize=12)
plt.ylabel(‘Loss’, fontsize=12)
plt.title(‘Training loss and validation loss’,fontsize=12)
plt.grid(True)
plt.legend([‘Training loss’,’Validation loss’],fontsize=12)
plt.style.use([‘classic’])
plt.figure(2,figsize=(7,5),facecolor=’white’)
plt.plot(x_axis,train_acc)
plt.plot(x_axis,val_acc)
plt.xlabel(‘Epoch’,fontsize=12)
plt.ylabel(‘Accuracy’,fontsize=12)
plt.title(‘Training accuracy and validation accuracy’,fontsize=12)
plt.grid(True)
plt.legend([‘Training accuracy’,’Validation accuracy’],fontsize=12,loc=4)
plt.style.use([‘classic’])


Evaluation¶
Now we need to test our network. In other words we evaluate it using the testing data set.
In [28]:
score = model.evaluate(X_test, y_test)
print(score) # will return the test loss and test accuracy
400/400 [==============================] – 0s 167us/step
[0.4578259766101837, 0.805]
We can even plot a confusion matrix.
In [29]:
from sklearn.metrics import confusion_matrix
# if you don’t have sklearn install it with: conda install scikit-learn
Y_pred = model.predict(X_test, batch_size=None, verbose=1)
y_pred = np.argmax(Y_pred, axis=1)
classes = [“LC”, “VU”, “EN”, “CR”] # these are the 4 possible classes
cm = confusion_matrix(np.argmax(y_test,axis=1), y_pred)
np.set_printoptions(precision=2)
fig = plt.figure(facecolor=’white’)
title = ‘Normalized confusion matrix’
cmap = plt.cm.Blues
cm = cm.astype(‘float’) / cm.sum(axis=1)[:, np.newaxis]
plt.imshow(cm, interpolation=’nearest’, cmap=cmap)
plt.title(title)
plt.colorbar()
#plt.colorbar(shrink=0.7) # alternative
tick_marks = np.arange(len(classes))
plt.xticks(tick_marks, classes, rotation=90, fontsize=8)
#plt.xticks(tick_marks, rotation=45, fontsize=6) # alternative
plt.yticks(tick_marks, classes, fontsize=8)
plt.tight_layout()
plt.ylabel(‘True label’)
plt.xlabel(‘Predicted label’)
400/400 [==============================] – 0s 283us/step
Out[29]:
Text(0.5, 25.374999999999986, ‘Predicted label’)

Prediction¶
Finally we used the final optimised network to predict the label for unknown entries. More info here.
In [30]:
y_ursus = model.predict(X_ursus)
print(y_ursus)
[[3.44e-05 8.72e-03 5.63e-01 4.28e-01]]
The nice thing is that these are proper posterior probabilities. Therefore we can even calculate Bayes factors and do model testing. Imagine that we want to test whether this species is not threatened (so not “least concern”). We can calculate this Baye factor as follows:
In [31]:
(1 – y_ursus[0,0])/(3/4) / (y_ursus[0,0])/(1/4)
Out[31]:
155033.13120789762
BYON (Build Your Own Network)¶
Can you do better than this network? Can you build a network which leads to a higher testing accuracy than this one? Try to build your own architecture and compiler. Some tips:
• you can use multiple cycles of Conv+MaxPool+Dropout layers
• you can play with the dropout rates to prevent overfitting
• you can tune the learning rate
• you can add more/less units (or neurons) in the dense layers
• you can add more/less filters or change their size (but be careful with it!)
• try a leakyReLu or other activation functions
• change how you initialise the weights
• …
Be aware that the accuracy won’t improve by default if you make the network deeper: you have more parameters to optimise with the same data set!
Ideally, you may want to try different configurations and retain the one with the highest validation accuracy. Then this network will be passed to the testing set. For instance, assume that you have one hyper-parameter and you want to estimate it. How would you build such pipeline?
The best way of learning is by doing. Read the keras manual https://keras.io/ to learn how to implement your ideas. Pick a partner or more and start a team.
It’s a competition.
Good luck!
In [ ]:
# …