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 :

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 : 

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 : 

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.

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 : 

Analyse des résultats

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 :

Analyse des résultats

Voyons nos résul­tats :

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.

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 : 

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

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 : 

Analyse des résultats

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  !