IA et bioinformatique : exploitons les réseaux convolutionnels (CNN)

Inter­es­sons-nous aujourd'hui aux séquences d’ADN. Nous uti­li­se­rons le data­set télé­char­geable ici : https://​www​.kaggle​.com/​d​a​t​a​s​e​t​s​/​n​a​g​e​s​h​s​i​n​g​h​/​d​n​a​-​s​e​q​u​e​n​c​e​-​d​a​t​a​set

L'ensemble des fichiers néces­saire à cet article sont dis­po­nibles ici.

Vous trou­ve­rez dans ce lot de don­nées un ensemble de séquences d’ADN issues de 3 espèces : l’homme, le chien et le chim­pan­zé. Cha­cune de ces séquences appar­tient à une des 7 familles de séquences pro­po­sées. Ces 7 classes cor­res­pondent res­pec­ti­ve­ment à :

• 0 : Récep­teurs cou­plés aux pro­téines G
• 1 : Recep­teurs-trans­mem­bra­naires à acti­vi­té tyro­sine kinase
• 2 : Tyro­sine phos­pha­tase
• 3 : Aspa­ra­gine syn­the­tase
• 4 : ADN mito­chon­drial
• 5 : Récep­teur vani­loïde
• 6 : Homéo-domaine TF1

L’idée sera donc ici de nous pla­cer dans un contexte de clas­si­fi­ca­tion à classes mul­tiples, qui dif­fère de la clas­si­fi­ca­tion binaire que nous avions réa­li­sée ici : https://​bioin​fo​-fr​.net/​b​i​o​i​n​f​o​r​m​a​t​i​q​u​e​-​e​t​-​i​a​-​u​n​-​p​r​e​m​i​e​r​-​pas.

Nous aurons recours à un enco­dage des séquences afin qu’elles puissent être ana­ly­sées par notre réseau. Je vous pro­pose aujourd’hui d’utiliser la théo­rie du chaos (ou jeux du chaos) que nous déve­lop­pe­rons un peu plus loin dans cet article. Pour le moment, étu­dions nos don­nées. 

Étude préalable des données

Nous nous attar­de­rons sur les séquences humaines qui sont au nombre de 4020, et que nous trai­te­rons selon une repré­sen­ta­tion en 5-mers (k=5), soit en matrices de 32x32 pixels. 

Comme pour toute approche d’apprentissage nous divi­sons notre jeu de don­nées en 3 sets : entraî­ne­ment, test et vali­da­tion selon le décou­page 80%, 10% et 10% (pour plus d'information sur le décou­page des don­nées c'est ici : https://fr.wikipedia.org/wiki/Jeux_d%27entrainement,_de_validation_et_de_test). Nous obtien­drons après redi­men­sion­ne­ment des ten­seurs de taille 3216x32x32x1 pour le jeu d’entraînement, 402x32x32x1 pour les jeux de test et de vali­da­tion avec une répar­ti­tion stra­ti­fiée et res­pec­tueuse des pro­por­tions des dif­fé­rentes classes. 

En nous inté­res­sant à la répar­ti­tion des classes de séquences à l’échelle du jeu de don­nées glo­bal, on remarque une sur­re­pré­sen­ta­tion des séquences de classe 0 et une sous repré­sen­ta­tion des classes 3 et 6.

Le code néces­saire au char­ge­ment des don­nées et à l'étude préa­lable des don­nées est le sui­vant :

from collections import defaultdict
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from sklearn.model_selection import train_test_split

import tensorflow as tf
import keras
from tensorflow.keras.layers import Conv2D, Dense, Dropout, BatchNormalization, MaxPool2D
from tensorflow.keras.models import load_model
from keras.utils import to_categorical

# SVM part

from sklearn.datasets import make_classification
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score

from sklearn.metrics import ConfusionMatrixDisplay

"""# Mount drive"""

from google.colab import drive
drive.mount('/content/drive')

"""# Load data"""

working_folder = "/content/drive/MyDrive/Sequences/"

human_sequences = "/content/drive/MyDrive/Sequences/human.txt"
dog_sequences = "/content/drive/MyDrive/Sequences/dog.txt"
chimpanze_sequences = "/content/drive/MyDrive/Sequences/chimpanzee.txt"

df_human = pd.read_csv(human_sequences, sep = "\t")
df_dog = pd.read_csv(dog_sequences, sep = "\t")
df_chimpanze = pd.read_csv(chimpanze_sequences, sep = "\t")


"""# Exploratory Data Analysis"""

labels = [1, 2, 3, 4, 5, 6]
plt.figure(figsize = (15, 10))
plt.subplot(131)
plt.pie(df_human["class"].value_counts(sort = False).sort_index() , autopct='%1.1f%%')
plt.title("Proportion of sequences classes in human")
plt.legend(title="Class",
           labels=labels,
          loc="lower right",)

plt.subplot(132)
labels_dog = df_dog["class"].unique()
plt.pie(df_dog["class"].value_counts(sort = False).sort_index() , autopct='%1.1f%%')
plt.title("Proportion of sequences classes in dog")
plt.legend(title="Class",
           labels=labels,
          loc="lower right",)

plt.subplot(133)
labels_chimpanze = df_chimpanze["class"].unique()
plt.pie(df_chimpanze["class"].value_counts(sort = False).sort_index() , autopct='%1.1f%%')
plt.title("Proportion of sequences classes in chimpanze")
plt.legend(title="Class",
           labels=labels,
          loc="lower right",)

plt.savefig(working_folder + "EDA.png")
plt.show()

La théorie du chaos

Très étu­diée en phy­sique depuis un paquet d’années, la théo­rie du chaos est étroi­te­ment liée à la créa­tion de frac­tales (https://​fr​.wiki​pe​dia​.org/​w​i​k​i​/​F​r​a​c​t​ale). 

Dans notre cas nous allons créer des repré­sen­ta­tions matri­cielles car­rées (gra­phiques) des séquences. Ces repré­sen­ta­tions à l'allure de frac­tales sont éga­le­ment étroi­te­ment liées à la com­po­si­tion en k-mers de nos séquences. Voyons cela plus en détail. 

Pour cela nous allons défi­nir un car­ré dont cha­cun des som­mets cor­res­pond à l’un des 4 fan­tas­tiques nucléo­tides qui font la richesse de notre métier, j’ai nom­mé A, T, C et G. Le centre du car­ré sera notre ori­gine de repère soit (0, 0), ain­si nous dis­po­se­rons les 4 nucléo­tides men­tion­nées aux coor­don­nées A (-1, 1), G (1, 1), T (1, -1) et enfin C (-1, -1). 

L’idée est de se dépla­cer au sein de ce car­ré, à mesure que nous par­cou­rons la séquence. Oui mais com­ment ? Et bien c’est très simple. En par­tant du centre du car­ré, nous allons nous dépla­cer jusqu’au point situé à mi-par­cours sur le seg­ment qui relie l’origine de notre repère (O(0, 0)) au pre­mier nucléo­tide de la séquence. Pre­nons l’exemple où T serait le pre­mier nucléo­tide. Alors nous nous dépla­ce­rons jusqu’au point médian entre (0, 0) et (1, -1), soit (0.5, -0.5). Appe­lons ce point Q puis diri­geons nous dans la direc­tion du nucléo­tide sui­vant. Ce sera un G, alors nous sui­vrons le seg­ment reliant Q (0.5, -0.5)  à G(1, 1), pour atteindre le point médian du seg­ment, soit le point R (0.75, 0.25). Le nucléo­tide d’après est un A. Nous par­cou­rons alors le seg­ment RA en direc­tion de A (-1, 1) jusqu’au point médian S(-0.125, 0.625), et ain­si de suite…

Ain­si nous pla­ce­rons un ensemble de points (autant que de nucléo­tides dans la séquence) dans le car­ré de départ. La repré­sen­ta­tion obte­nue sera alors rela­tive à l'abondance de l’ensemble des k-mers ( à k préa­la­ble­ment défi­ni).

Les repré­sen­ta­tions numé­riques des espaces n’étant pas conti­nues mais dis­crètes, il convient de défi­nir notre espace numé­rique, autre­ment dit la taille de la matrice de sor­tie. Pour cela, il convient de défi­nir la taille k  des oli­go­mères que nous consi­dé­re­rons. Ain­si la matrice (ou espace) de sor­tie aura une car­di­na­li­té (ou aire) de 4^k élé­ments, soit un car­ré avec un côté de lon­gueur (sqrt(4^k)). 

Le par­cours que nous avons vu pré­cé­dem­ment per­met ain­si de loca­li­ser les k-mers au sein de cette matrice. En effet, si nous repre­nons le che­min préa­la­ble­ment sui­vi nous nous diri­geons d’abord vers le cadran infé­rieur droit du car­ré pour le T, puis dans le cadran supé­rieur droit du cadran pré­cé­dent pour le G, puis dans le cadran supé­rieur gauche de ce der­nier cadran et ain­si de suite. À chaque nucléo­tide par­cou­ru le cadran cou­rant est alors re-divi­sé en 4 en conser­vant les orien­ta­tions rela­tives aux som­mets des nucléo­tides. 

On obser­ve­ra alors que l’ensemble des séquences, et a for­tio­ri des k-mers com­men­çant par TG se situe­ront dans le cadran supé­rieur droit du cadran infé­rieur droit. Selon cette découpe suc­ces­sive en 4 cadrans, il est alors pos­sible de pla­cer l’ensemble des k-mers dans espace préa­la­ble­ment défi­ni (ou élé­ments de la matrice i.e. pixels de l’image). 

Pour un espace repré­sen­ta­tif des 2-mers et des 3-mers, nous aurions les décou­pages sui­vants : 

Très bien, mais quelles valeurs affec­ter à ces cadrans ou élé­ments de la matrice ? Et bien main­te­nant que l’ensemble des espaces dédiées aux k-mers sont défi­nis, nous affec­te­rons aux élé­ments cor­res­pon­dants dans la matrice, la pro­ba­bi­li­té d’apparition de ces k-mers dans la matrice. Pour plus d’informations sur les k-mers et le cal­cul des pro­ba­bi­li­té d’occurence, je vous ren­voie ici (https://fr.wikipedia.org/wiki/K-m%C3%A8re).
━━

Ain­si, dans le cadre d’une repré­sen­ta­tion en 3-mers, nous obtien­drons une matrice de dimen­sion 8x8 dont les élé­ments auront une valeur com­prise entre 0 et 1. Pour notre tuto­riel ce sera des 5-mers avec une repré­sen­ta­tion 32x32.

Pour en savoir plus sur les repré­sen­ta­tions de la théo­rie du chaos c’est ici : https://​pub​med​.ncbi​.nlm​.nih​.gov/​2​3​3​6​3​93/.

Le code néces­saire à la pro­duc­tion de telles repré­sen­ta­tion est le sui­vant : 

"""## Sanitize data
### Human
"""
## Delete sequences that contains N
nb_sequences_to_drop = 0
rows_indexes_to_drop = list()
for idx, seq in enumerate(df_human['sequence']):
    if 'N' in seq:
      nb_sequences_to_drop +=1
      # display(df.loc[df['sequence'] == seq])
      rows_indexes_to_drop.append(idx)
print(f"Number of sequence to be dropped : {nb_sequences_to_drop}")
print(f"Rows to be dropped : {rows_indexes_to_drop}")
print(f"Shape of the data before dropping sequences containing Ns : {df_human.shape}")
df_human.drop(rows_indexes_to_drop, inplace = True)
print(f"Shape of the data after dropping sequences containing Ns : {df_human.shape}")

"""### Dog"""

nb_sequences_to_drop = 0
rows_indexes_to_drop = list()
for idx, seq in enumerate(df_dog['sequence']):
    if 'N' in seq:
      nb_sequences_to_drop +=1
      rows_indexes_to_drop.append(idx)
print(f"Number of sequence to be dropped : {nb_sequences_to_drop}")
print(f"Rows to be dropped : {rows_indexes_to_drop}")
print(f"Shape of the data before dropping sequences containing Ns : {df_dog.shape}")
df_dog = df_dog.drop(rows_indexes_to_drop, inplace = True)

"""### Chimpanze"""

nb_sequences_to_drop = 0
rows_indexes_to_drop = list()
for idx, seq in enumerate(df_chimpanze['sequence']):
    if 'N' in seq:
      nb_sequences_to_drop +=1
      rows_indexes_to_drop.append(idx)
print(f"Number of sequence to be dropped : {nb_sequences_to_drop}")
print(f"Rows to be dropped : {rows_indexes_to_drop}")
print(f"Shape of the data before dropping sequences containing Ns : {df_human.shape}")
df_chimpanze.drop(rows_indexes_to_drop, inplace = True)

"""## CGR transformation"""

def cgr_encoding(sequence, k):
  """
  sequence : Sequence to encode as CGR
  k = k-mers size
  """

  kmer_count = defaultdict(int)
  for i in range(len(sequence)-(k-1)):
      kmer_count[str(sequence[i:i+k])] +=1
  for key in kmer_count.keys():
      if "N" in key:
          del kmer_count[key]
  probabilities = defaultdict(float)
  N = len(sequence)
  for key, value in kmer_count.items():
      probabilities[key] = float(value) / (N - k + 1)
  array_size = int(np.sqrt(4**k))
  chaos = []
  for i in range(array_size):
      chaos.append([0]*array_size)
  maxx = array_size
  maxy = array_size
  posx = 1
  posy = 1
  for key, value in probabilities.items():
      for char in key:
          if char == "T":
              posx += maxx / 2
          elif char == "C":
              posy += maxy / 2
          elif char == "G":
              posx += maxx / 2
              posy += maxy / 2
          maxx = maxx / 2
          maxy /= 2
      chaos[int(posy)-1][int(posx)-1] = value
      maxx = array_size
      maxy = array_size
      posx = 1
      posy = 1
  return tf.constant(chaos)

Voi­ci ce à quoi res­semblent nos séquences (non, non il ne s’agit pas de la TV qui ne capte pas).

Le code néces­saire à la pro­duc­tion de ces images est le sui­vant : 

"""## Plot data"""
# Randomly selected one example per class
random_idx_class_0 = np.random.choice(df_human[df_human["class"] == 0].index)
random_idx_class_1 = np.random.choice(df_human[df_human["class"] == 1].index)
random_idx_class_2 = np.random.choice(df_human[df_human["class"] == 2].index)
random_idx_class_3 = np.random.choice(df_human[df_human["class"] == 3].index)
random_idx_class_4 = np.random.choice(df_human[df_human["class"] == 4].index)
random_idx_class_5 = np.random.choice(df_human[df_human["class"] == 5].index)
# Plot some data
plt.figure(figsize = (10, 10))
plt.subplot(231)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_0], k = 5), cmap = "gray")
plt.title("CGR of class 0")
plt.subplot(232)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_1], k = 5), cmap = "gray")
plt.title("CGR of class 1")
plt.subplot(233)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_2], k = 5), cmap = "gray")
plt.title("CGR of class 2")
plt.subplot(234)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_3], k = 5), cmap = "gray")
plt.title("CGR of class 3")
plt.subplot(235)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_4], k = 5), cmap = "gray")
plt.title("CGR of class 4")
plt.subplot(236)
plt.imshow(cgr_encoding(sequence = df_human["sequence"][random_idx_class_5], k = 5), cmap = "gray")
plt.title("CGR of class 5")
plt.subplots_adjust(wspace=0.12, hspace=0)
plt.show()

Il appa­raît que les repré­sen­ta­tions des séquences d’ADN sont carac­té­ris­tiques des espèces, des types de séquences, etc. Voyons alors com­ment uti­li­ser ces repré­sen­ta­tions dans notre objec­tif de clas­si­fi­ca­tion. 

Pour ce faire, nous allons uti­li­ser un petit réseau à 2 blocs convo­lu­tifs sui­vis d’une couche dense et d’une seconde couche dense de clas­si­fi­ca­tion.

Afin d'entraîner notre réseau, il convient de pré­pa­rer les don­nées afin qu’elles puissent être four­nies au réseau. Il convient aus­si de divi­ser le jeu de don­née de départ en 3 sets : entrai­ne­ment, test et vali­da­tion.

"""## Learning part

### Generate data
"""

encoded_list = []

for seq in df_human["sequence"]:
  encoded_list.append(cgr_encoding(seq, k =5))

X, y  = np.array(encoded_list), to_categorical(df_human["class"])

"""### Data spliting"""

X = X.reshape(4020, 32, 32, 1)

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 42, stratify = y, shuffle = True, test_size = 0.2) # Split data into train and test with respect of class proportions

X_test, X_validation, y_test, y_validation = train_test_split(X_test, y_test, random_state = 42, stratify = y_test, shuffle = True, test_size = 0.5) # Split data into test and validation with respect of class proportions

Paramètres du réseau

Comme men­tion­né pré­cé­dem­ment , nous alors uti­li­ser un petit réseau de type convo­lu­tion­nel. Contrai­re­ment au réseau uti­li­sé la fois der­nière, celui ci sera basé sur des bloc convo­lu­tion­nels 2D, néna­moins le fonc­tion­ne­ment est simi­laire. Pour quelques rap­pels sur la convo­lu­tion c’est ici : https://medium.com/@CharlesCrouspeyre/comment-les-r%C3%A9seaux-de-neurones-%C3%A0-convolution-fonctionnent-b288519dbcf8.

Notre réseau com­pren­dra donc deux blocs convo­lu­tion­nels (une couche de convo­lu­tion sui­vie d’une couche de poo­ling elle même sui­vie par une couche de dro­pout), les­quels seront sui­vis par une pre­mière couche dense conte­nant et apla­tis­sant les filtres pré­cé­dem­ment obte­nus par convo­lu­tion, et une couche dense à 7 neu­rones pour la clas­si­fi­ca­tion. 

Contrai­re­ment à la fois pré­cé­dente, la couche de clas­si­fi­ca­tion n’utilisera pas la fonc­tion sig­moïde comme acti­va­tion, mais la fonc­tion soft­max afin de rendre une pro­ba­bi­li­té d'appartenance à une classe don­née plu­tôt qu’une réponse de type “oui/​non”.

Pour ce qui concerne les para­mètres d’optimisation, nous uti­li­se­rons une des­cente de gra­dient (pour l’ajustement des poids du réseau) de type ADAM avec un taux d’apprentissage de 0.001. La fonc­tion de perte sera la “cate­go­ri­cal cross entro­py” afin de tenir compte de l’erreur dans un contexte de clas­si­fi­ca­tion à plu­sieurs classes.

Le code néces­saire à la construc­tion du réseau et à son entrai­ne­ment est le sui­vant : 

model = tf.keras.Sequential()

model.add(Conv2D(filters = 32, kernel_size = 3, activation = 'relu', padding = "valid", input_shape = (32, 32, 1), name = "Conv_1")) # Input shape : Batch_size, width, height, channels
model.add(Dropout(0.3, name = "Dropout_1"))
model.add(MaxPool2D(pool_size=4, strides=None, padding='valid', name = "MaxPool_1"))
model.add(Conv2D(filters = 64, kernel_size = 2, activation = 'relu', padding = "valid", name = "Conv_2")) # Input shape : Batch_size, width, height, channels
model.add(Dropout(0.3, name = "Dropout_2"))
model.add(tf.keras.layers.Flatten(name = "Flatten"))
model.add(tf.keras.layers.Dense(units = 7, activation = 'softmax', name = "Output"))

model.summary()

# Compile model
model.compile(loss='categorical_crossentropy',
              optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
              metrics=['accuracy'])

# Define early stopping to avoid overfitting

early_stop = keras.callbacks.EarlyStopping(monitor = 'val_accuracy', min_delta = 0.0005, patience=8,
                                           restore_best_weights=True )

# Learning from data
history = model.fit(X_train, y_train, batch_size = 8, validation_data= (X_test, y_test),epochs=200)

Analyse des résultats

"""## Plot results"""

plt.figure(figsize = (10,4))
plt.subplot(121)
plt.plot(history.epoch, history.history["loss"], 'g', label='Training loss')
plt.plot(history.epoch, history.history["val_loss"], 'b', label='Validation loss')
plt.title('Training loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()

plt.subplot(122)
plt.plot(history.epoch, history.history["accuracy"], 'r', label='Training accuracy')
plt.plot(history.epoch, history.history["val_accuracy"], 'orange', label='Validation accuracy')
plt.title('Training accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.savefig(working_folder + "base_modelpng")
plt.show()

"""### Save model"""

working_folder = "/content/drive/MyDrive/Sequences/"
model.save(working_folder)

"""### Model evaluation"""

loss, acc = model.evaluate(X_validation, y_validation, verbose=2)
print("Basic model, accuracy: {:5.2f}%".format(100 * acc))
"""## Plot results"""

plt.figure(figsize = (10,4))
plt.subplot(121)
plt.plot(history.epoch, history.history["loss"], 'g', label='Training loss')
plt.plot(history.epoch, history.history["val_loss"], 'b', label='Validation loss')
plt.title('Training loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()

plt.subplot(122)
plt.plot(history.epoch, history.history["accuracy"], 'r', label='Training accuracy')
plt.plot(history.epoch, history.history["val_accuracy"], 'orange', label='Validation accuracy')
plt.title('Training accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.savefig(working_folder + "base_modelpng")
plt.show()

"""### Save model"""

working_folder = "/content/drive/MyDrive/Sequences/"
model.save(working_folder)

"""### Model evaluation"""

loss, acc = model.evaluate(X_validation, y_validation, verbose=2)
print("Basic model, accuracy: {:5.2f}%".format(100 * acc))

Notre entraî­ne­ment est ter­mi­né, que pou­vons-nous en déduire ?

La qua­li­té glo­bale de l'entraînement est plu­tôt satis­fai­sante. Nous attei­gnons un début de pla­teau au niveau de la jus­tesse, avec 0.9 au bout de 200 epochs (phases d’apprentissages) en entraî­ne­ment. La jus­tesse en vali­da­tion est moindre, mais conve­nable. Notre réseau semble plu­tôt bien avoir appris. En éva­luant le modèle avec le set de don­née de vali­da­tion, nous obte­nons une per­for­mance glo­bale de 74% ce qui est rela­ti­ve­ment satis­fai­sant, mais pour­rait-on faire mieux ? Voyons si le trans­fert lear­ning (ou appren­tis­sage trans­fé­ré) pour­rait nous aider.

Tranfert d'apprentissage

Le trans­fert lear­ning repré­sente l'utilisation d’un appren­tis­sage préa­la­ble­ment réa­li­sé afin de répondre à notre ques­tion. Pour autant nous pour­rions être les pre­miers à nous poser cette ques­tion de clas­si­fi­ca­tion de séquences, et encore plus si l’on consi­dère notre enco­dage selon la théo­rie du chaos.

En pra­tique, le trans­fert lear­ning est plu­tôt mis en place afin de répondre à des contraintes de temps et/​ou de per­for­mances. Il s’agit de réuti­li­ser des modèles sem­blables au nôtre et ayant déjà été entraî­nés et don­nant des résul­tats satis­fai­sants. Vous me direz que tout cela est bien joli, mais com­ment trou­ver un modèle stric­te­ment iden­tique au mien ? Et bien il n’est pas néces­saire que ce modèle soit stric­te­ment iden­tique. Il faut sim­ple­ment que les dimen­sions d’entrée du réseau soient com­pa­tibles avec nos don­nées, ou alors, adap­ter nos don­nées au réseau pré-entraî­né. 

Se pose alors la ques­tion de la ques­tion auquel répond notre réseau ver­sus celle à laquelle répond le réseau pré-entrai­né. Pour répondre à ces inter­ro­ga­tions, je vous pro­pose de plon­ger direc­te­ment dans notre exemple.

Nous allons uti­li­ser un réseau convo­lu­tion­nel déjà entraî­né sur le set de don­nées MNIST et l’adapter à nos besoins. MNIST ? Vous l’avez sans doute déjà croi­sé il s’agit de ce célèbre jeu de don­nées conte­nant des images de chiffres écrits à la main, célèbre dans tous les tuto­riaux d’apprentissage pro­fond. En prin­cipe, ces images sont de taille 28x28x1. Je me suis per­mis d’adapter ce réseau afin qu’il accepte nos repré­sen­ta­tions frac­tales de séquences de dimen­sions 32x32x1. Le réseau que nous allons uti­li­ser est une adap­ta­tion du réseau LeNet5 pro­po­sé par Yann Lecun (plus d'informations).

Cepen­dant, une dif­fé­rence majeure se pré­sente. En effet, la couche de clas­si­fi­ca­tion de ce réseau com­porte 10 neu­rones afin de rendre compte de l’appartenance de l’une des images de chiffres à une classe entre 0 et 9. Pas de panique nous allons adap­ter cela à nos besoins.
Voyons d’abord dans les grandes lignes ce que nous allons faire. Tout d’abord il va s’agir de rechar­ger le modèle déjà entraî­né. (pour plus d’informations sur le modèle et sa concep­tion, je vous ren­voie vers ce note­book : https://​github​.com/​b​i​o​i​n​f​o​-​f​r​/​d​a​t​a​_​i​a​/​b​l​o​b​/​m​a​i​n​/​a​r​t​i​c​l​e​_​i​a​_​2​/​L​e​N​e​t​5​.​i​p​ynb).

Nous sou­hai­tons conser­ver l’ensemble des couches du réseau, mais aus­si pro­fi­ter des poids déjà pré-éta­blis par l'entraînement pré­cé­dent. Pour autant le réseau n’est pas prêt à l’emploi pour notre pro­blé­ma­tique. Il va s’agir de “fine-tuner” ce modèle afin de tirer par­tie des des­crip­teurs (filtres) éta­blis lors de l'entraînement sur la base de don­nées MNIST, en n'entraînant (ajus­tant les poids) que d’un sous-ensemble de couches de notre réseau. Dans notre exemple, le LeNet étant un modèle plu­tôt petit, nous ne conser­ve­rons que les filtres de la pre­mière couche. Il va alors fal­loir “geler” les poids (coef­fi­cients des filtres convo­lu­tifs) de la pre­mière couche et lais­ser le réseau n’ajuster que les autres poids.

Par ailleurs, nous avons évo­qué que le LeNet5 per­met une clas­si­fi­ca­tion dans un contexte de 10 classes. Nous allons alors ôter la der­nière couche du réseau dont le rôle est de réa­li­ser cette clas­si­fi­ca­tion, afin de la rem­pla­cer avec une couche dense de clas­si­fi­ca­tion à 7 classes adap­tée à notre pro­blème, qui sera à entraî­ner. Nous allons gar­der les mêmes para­mètres d’optimisation que ceux de notre pre­mier réseau :  ADAM avec un taux d’apprentissage de 0.001, fonc­tion de perte de type “cate­go­ri­cal cross entro­py”.

Le code néces­saire à la réa­li­sa­tion de ces opé­ra­tions est le sui­vant :

"""## Load LeNet5"""

lenet_folder = "/content/drive/MyDrive/LeNet/"

lenet_model = load_model(lenet_folder)

"""### Check LeNet model"""

lenet_model.summary()

"""### Freeze bottom layers for transfer learning"""

# Freeze the first layers
for idx, layer in enumerate(lenet_model.layers):
  if idx <= 1:
    layer.trainable = False

"""### Pop last dense layer and replace it for relevent goal"""

# Remove last layer
lenet_model.pop()

# Replace it with a more appropriate
lenet_model.add(tf.keras.layers.Dense(units = 7, activation = 'softmax', name = "Output"))

"""### Check modified model"""

lenet_model.summary()

"""### Compile model and train"""

lenet_model.compile(loss='categorical_crossentropy',
              optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),
              metrics=['accuracy'])

# Define early stopping to avoid overfitting
early_stop = keras.callbacks.EarlyStopping(monitor = 'val_accuracy', min_delta = 0.0005, patience=8,
                                           restore_best_weights=True )

history_lenet = lenet_model.fit(X_train, y_train, batch_size = 8, validation_data= (X_test, y_test),
                        epochs=200)

transfer_learning_directory = "/content/drive/MyDrive/Sequences/transfert_learning/"


lenet_model.save(transfer_learning_directory)

Analyse des résultats

Voyons nos résul­tats :

"""# Plot final results"""

plt.figure(figsize = (14,4))
plt.subplot(121)
plt.plot(history.epoch, history.history["loss"], 'g', linestyle = 'dotted', label='Training - LeNet')
plt.plot(history_lenet.epoch, history_lenet.history["loss"], 'g', label='Training - LeNet')
plt.plot(history.epoch, history.history["val_loss"], 'b', label='Validation loss')
plt.plot(history_lenet.epoch, history_lenet.history["val_loss"], 'b' , linestyle = 'dotted', label='Validation  - LeNet')
plt.title('Training loss')
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()

plt.subplot(122)
plt.plot(history.epoch, history.history["accuracy"], 'r', label='Training accuracy')
plt.plot(history_lenet.epoch, history_lenet.history["accuracy"], 'r',linestyle = 'dotted', label='Training - LeNet')

plt.plot(history.epoch, history.history["val_accuracy"], 'orange', label='Validation accuracy')
plt.plot(history_lenet.epoch, history_lenet.history["val_accuracy"], 'orange', linestyle = 'dotted', label='Validation - LeNet')

plt.title('Training accuracy')
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()
plt.savefig(transfer_learning_directory + "dna_lenet_model.png")
plt.show()

"""### Evaluate model"""

loss, acc = lenet_model.evaluate(X_validation, y_validation, verbose=2)
print("DNA LeNet model, accuracy: {:5.2f}%".format(100 * acc))

On note ici deux dif­fé­rences avec le modèle pré­cé­dent. La pre­mière étant que la jus­tesse à l’issue de nos 200 epochs est légè­re­ment supé­rieure en uti­li­sant notre DNA-LeNet. La jus­tesse en vali­da­tion est sem­blable. Cepen­dant la valeur de la fonc­tion de perte sur la fin de l'entraînement semble s’affoler pour ce même réseau.

Par ailleurs, on observe un petit pla­teau au début de l'entraînement de nos deux modèles. Tou­te­fois, le pla­teau dans le cadre de notre DNA-LeNet est plus court, le réseau se met plus rapi­de­ment à apprendre, mais apprend plus len­te­ment que notre pre­mière ver­sion de réseau. 

En éva­luant ce nou­veau modèle avec le set de don­nées de vali­da­tion, nous obte­nons une per­for­mance glo­bale de 72%. Il semble que chan­ger de modèle ne nous ait que peu aidé.

Peut être que cela est impu­table à notre couche de clas­si­fi­ca­tion qui peine à déduire la bonne classe de l’ensemble des des­crip­teurs four­nis par les couches convo­lu­tives ? Voyons si nous pou­vons clas­ser nos séquences autre­ment que via cette couche.

Couplage apprentissage profond et apprentissage machine

À la dif­fé­rence de l’apprentissage machine, l’apprentissage pro­fond pré­sente l’avantage de nous affran­chir par­tiel­le­ment de l’extraction de des­crip­teurs de nos don­nées. En effet, dans le cadre d’un réseau convo­lu­tif, ces der­niers sont appris. Il s’agit des filtres de convo­lu­tion dont les coef­fi­cients sont ajus­tés tout au long de l'entraînement.

Pour­rait-on alors ne gar­der que ces des­crip­teurs, et les four­nir à un algo­rithme de clas­si­fi­ca­tion plus conven­tion­nel ? C’est ce que nous allons voir. 

Nous uti­li­se­rons la SVM (Machine à Sup­port de Vec­teur). Il s’agit d’un algo­rithme de clas­si­fi­ca­tion binaire dont le prin­cipe géné­ral est assez simple. Il va s’agir de pla­cer nos don­nées dans un espace vec­to­riel de grande dimen­sion afin d'établir l’équation d’un hyper­plan, lequel repré­sen­te­ra la fron­tière sépa­rant au mieux les dif­fé­rentes classes liées au dif­fé­rents points dans cet espace. Autre­ment dit dans notre cas, sépa­rant l’ensemble des repré­sen­ta­tions des dif­fé­rentes classes de séquence, en maxi­mi­sant la marge entre ces repré­sen­ta­tions à la marge. 

Puisqu’il s’agit d’un hyper­plan sépa­ra­teur, nous ne pour­rons que sépa­rer deux classes à la fois. Néan­moins les SVMs offrent une stra­té­gie appe­lée “one vs one” qui per­met de construire un sys­tème de clas­si­fi­ca­tion à plu­sieurs classes. La sépa­ra­tion de classe va donc s’effectuer deux à deux dans notre espace cher­chant à d’abord sépa­rer la classe 0 de la classe 1, puis la classe 0 de la classe 2, puis la classe 0 de la classe 3 et ain­si de suite. Nous allons alors entraî­ner autant de clas­si­fieurs binaires qu’il existe de com­bi­nai­sons de classes. 

Une fois ces clas­si­fieurs entraî­nés, l’étape de clas­si­fi­ca­tion se fera ain­si. Une des repré­sen­ta­tions sera pla­cée dans l’espace de clas­si­fi­ca­tion. Une classe sera attri­buée selon cha­cun des hyper­plans cal­cu­lés. Pour une repré­sen­ta­tion don­née, nous aurons donc un ensemble de classes plau­sibles. La classe pré­dite sera alors déter­mi­née à la majo­ri­té. En d’autres termes, si la classe 0 est pré­dite 9 fois, la classe 1 elle 3 fois et la classe 5 à son tour 5 fois, alors la classe 0 sera affec­tée à notre repré­sen­ta­tion. 

Trêve de théo­rie, pas­sons à la pra­tique. 

En pre­mier lieu il s’agit de trans­for­mer notre réseau convo­lu­tion­nel ini­tial en extrac­teur de des­crip­teurs. Pour cela nous allons reti­rer la couche de clas­si­fi­ca­tion. Ain­si la sor­tie du réseau sera la sor­tie de la couche “Flat­ten” soit un vec­teur de dimen­sion 2304 par séquence.

"""## Change classification layer for SVM

### Adapt intial CNN for SVM classification
"""

# Chop off head of CNN

modified_model = model
modified_model.pop()
modified_model.summary()

L’étape sui­vante consiste en le ré-enco­dage des don­nées de sor­ties. On va alors four­nir l’ensemble de nos repré­sen­ta­tions CGR à notre réseau afin d’obtenir les vec­teurs cor­res­pon­dants, et ce pour les sets d'entraînement, de test et de vali­da­tion : 

"""### Prepare data for svm"""
# Go back on labels from one hot encoding
X_train_svm  = modified_model.predict(X_train)
y_train_svm = np.argmax(y_train, axis=1)

X_test_svm = modified_model.predict(X_test)
y_test_svm = np.argmax(y_test, axis=1)

X_validation_svm = modified_model.predict(X_validation)
y_validation_svm = np.argmax(y_validation, axis=1)

"""### Checking dimensions"""

print(f"Shape of x_train_svm : {X_train_svm.shape}")
print(f"Shape of y_train_svm : {y_train_svm.shape}")

print(f"Shape of x_test_svm : {X_test_svm.shape}")
print(f"Shape of y_test_svm : {y_test_svm.shape}")

print(f"Shape of x_validation_svm : {X_validation_svm.shape}")
print(f"Shape of y_validation_svm : {y_validation_svm.shape}")

Enfin nous allons défi­nir notre clas­si­fieur, soit une SVM en mode “one vs one” et entraî­ner ce modèle : 

"""### Instanciate SVM classifier and train it"""

SVM_classifier = SVC(decision_function_shape = 'ovo', class_weight = "balanced", C = 125.0)

Notre modèle étant à pré­sent entraî­né en uti­li­sant les des­crip­teurs extraits par le le réseau convo­lu­tion­nel ini­tial modi­fié, exa­mi­nons les résul­tats : 

SVM_classifier.fit(X_train_svm, y_train_svm,)

Analyse des résultats

"""### Evaluate classification"""

# Make predictions on test set

y_pred_svm = SVM_classifier.predict(X_test_svm)

accuracy = accuracy_score(y_test_svm, y_pred_svm)

print(f"Acuracy using CNN coupled to SVM : {accuracy}")

"""## Confusion matrix"""

class_names = np.arange(7)

# Plot non-normalized confusion matrix
titles_options = [
    ("Confusion matrix, without normalization", None),
    ("Normalized confusion matrix", "true"),
]
for title, normalize in titles_options:
    disp = ConfusionMatrixDisplay.from_estimator(
        SVM_classifier,
        X_test_svm,
        y_test_svm,
        display_labels=class_names,
        cmap=plt.cm.Blues,
        normalize=normalize,
    )
    disp.ax_.set_title(title)

    print(title)
    print(disp.confusion_matrix)

plt.savefig(working_folder + "svm_confusion_matrix.png")
plt.show()

Comme toute pro­cé­dure d’apprentissage machine, il aura fal­lu para­mé­trer les hyper-para­mètres du modèle. Ici, je ne me suis inté­res­sé qu’au coef­fi­cient “C” qui est le coef­fi­cient de marge d’un clas­si­fieur type SVM. 

Les résul­tats sont ici très inté­res­sants. Nous obte­nons une jus­tesse de 83.8% contre envi­ron 70-75% avec nos deux pré­cé­dentes démarches. 

En obser­vant la matrice de confu­sion, on peut s'apercevoir que les classes 0, 1 et 6 sont bien sépa­rées des autres, tan­dis que les classes 2 et 3 semblent être plus dif­fi­ci­le­ment dis­tin­guables. Néan­moins notre sys­tème de clas­si­fi­ca­tion cou­plant CNN et SVM semble effi­cace. 

Conclusion

Au tra­vers de ce tuto­riel, nous aurons par­cou­ru la chaîne stan­dard d’une (mini) ana­lyse basée sur l’apprentissage auto­ma­tique. 

La pre­mière étape aura consis­té en l’import des don­nées, une ana­lyse explo­ra­toire des don­nées et leur enco­dage selon la théo­rie des jeux du chaos. Dans une deuxième étape nous aurons entraî­né un petit réseau de convo­lu­tion avec nos séquences trans­for­mées en images. 

Dans une troi­sième étape, nous aurons ten­té d’améliorer les per­for­mances de ce modèle, en uti­li­sant un autre modèle du même aca­bit entraî­né lui, sur des images de chiffres, afin d’aborder la notion de trans­fert d’apprentissage. 

Il sem­ble­rait que les étapes de clas­si­fi­ca­tion de ces deux réseaux pré­cé­dents aient du mal à rem­plir leur rôle. Cela peut éven­tuel­le­ment être lié à l'entraînement pour lequel 200 epochs ne seraient pas suf­fi­sants. En effet, on peut obser­ver sur les courbes d’apprentissage que le pla­teau n’est pas tout à fait atteint. 

Enfin dans une der­nière étape, nous aurons cher­ché à amé­lio­rer notre clas­si­fieur, en uti­li­sant notre réseau convo­lu­tion­nel ini­tial en tant qu’extracteur de des­crip­teurs et en réa­li­sant la clas­si­fi­ca­tion avec une machine à vec­teur de sup­port, abor­dant ain­si le cou­plage appren­tis­sage pro­fond et appren­tis­sage machine. 

L'ensemble des don­nées, des scripts, note­books et sau­ve­gardes de modèles sont dis­po­nibles sur ce repo git.

N.B. : Si vous repro­dui­sez ce tuto­riel il se peut que vos résul­tats dif­fèrent de ceux pré­sen­tés ici en rai­son de l'aspect sto­chas­tique des pro­cé­dures d'apprentissage et découpe du jeu de don­née ini­tial.

Un grand mer­ci à nos relec­teurs pour cet article : Léo­pold, Samuel, Aze­rin & ZaZo0o  !