Après vous avoir montré tout un tas de choses sur comment visualiser et extraire une carte, je vais vous montrer aujourd'hui comment une carte de contact chromosomique peut être observée sous la forme d'une structure 3D à partir de VMD. Pour cela nous allons partir de la méthode développée dans mon laboratoire : Shrec3D.
Cet article étant relativement technique, il se repose sur de nombreux concepts déjà expliqués dans d'autres articles du blog que je vous recommande fortement de lire : comprendre ce qu'on visualise, comment on traite les données en amont. On présentera aujourd'hui un tutoriel détaillant la méthodologie de réalisation d'un type de visualisation. Un autre article viendra ensuite prendre du recul et comparer ce type de résultat à d'autres méthodes.
De nombreux exemples de code seront glissés un peu partout dans l'article. Ces morceaux de code sont écrits en MATLAB mais restent facile à lire pour toutes personnes connaissant Python ou équivalent.
Explication (simplifiée) de l'algorithme
Transformation de la carte de contact en carte de distance
Pour cette partie, je vais directement partir d'un morceau de code de Shrec3D, dont l'implémentation la plus récente se trouve sur le dépôt github éponyme. J'utiliserai le chromosome 19 entier de la drosophile, sur les données de cellules embryonnaires issues de bonev à une résolution de 10kb. Une fois cette carte de contact chargée dans MATLAB la transformation peut commencer. À cette étape, nos données ressemblent à ça :
1 2 3 4 5 6 |
matrix=h5read('chr19.hdf5','/data'); %Lecture des données, que j'ai déjà filtré, le /data est la clef dans mon fichier hdf5 matnorm=SCN_sumV2(matrix); %Normalisation des données avant traitement figure,imagesc(log10(matnorm)) axis square colormap(jet) |
Ce que nous observons actuellement, c'est une carte de contact (rappel : chaque point lumineux dans l'image précédente représente la fréquence où deux régions ont été retrouvées en proximité spatiale/contact). Notre but c'est d'obtenir une représentation en 3D de notre chromosome, il nous faut donc la distance entre chaque région de notre carte. La première étape va donc être de réaliser une transformation de notre matrice de contact en matrice de distance. L'opération mathématique pour effectuer cette opération est de prendre l'inverse de chaque élément élevé à une puissance arbitraire alpha.
1 2 3 4 5 |
alpha=0.1 ; distmat = 1./matnorm.^alpha ; figure,imagesc(distmat) axis square |
À quoi sert l'alpha mit dans le code ?
Ce facteur est là pour moduler la proportion de contacts initiale utilisés pour réaliser la transformation. Cette valeur est généralement comprise entre 0,1 et 1. Plus alpha est grand, moins de valeurs seront utilisées (ne laissant que des informations de courtes proximités entre chaque région). Je recommande de mettre un alpha autour de 0,2. L'impact du alpha sur la structure 3D est trouvable ici.
Quand on regarde les deux images précédentes, on remarque tout de suite un problème. La matrice de distance est très… incomplète. La carte initiale contenait beaucoup de zéros (que je colorie en jaune sur l image suivante), ce sont des points dont la valeur après transformation en matrice de distance est "inf" (1/0^0.1 ça donne rien!). Pour tous ces points manquants, on va devoir inférer ce qui se passe. Pour cela on va faire de l'imputation en utilisant un algorithme du plus court chemin entre tous les points. L'idée derrière est assez simple, imaginez une carte de la France avec seulement la moitié des distances entre les villes dessus. Ici, on va simplement compléter l'information sur la carte en inférant les distances qui manquent à partir de celles qu'on a.
1 |
FastFloydedmat=FastFloyd(distmat); %Transformation en carte de distance et algorithme du plus court chemin (FastFloyd) |
Transformation de la carte de distance en coordonnée 3D
Maintenant, on a la distance entre tous les points de notre carte. Le problème ici c'est le nombre de dimension de cette matrice qui est, basiquement, son nombre de ligne/colonne. Pour avoir une structure 3D, on va donc devoir effectuer une réduction de dimension de la carte à l'aide de la fonction mdscale.
1 2 |
distmat = FastFloydedmat <span class="pl‑k">-</span> diag(diag(FastFloydedmat)); %Pour que mdscale marche, il faut que la diagonale de la matrice en entrée soit vide. XYZ = mdscale(distmat,3, 'criterion','strain'); |
À cette étape, le plus dur est fait ! Nous avons en effet transformé notre carte de contact en matrice de distance, puis en un système 3D ou XYZ nous donnant la position relative de chaque point dans un espace.
Le paramètre sur lesquels il est possible de jouer : le choix du critère de réduction de dimension.
Ce que fait cette fonction, c'est chercher à réduire le nombre de dimensions dans la matrice de distance au nombre indiqué, ici 3. Pour cela, cette fonction va chercher à minimiser l'erreur entre les vecteurs à 3 dimensions qu'elle aura trouvé et ceux à N dimensions donnés en entrée. Le critère par défaut est 'strain' et donne des résultats convaincants. 'Sammon' aura tendance à minimiser les écarts relatifs entre distances calculées dans les espaces ND et 3D. Donc dans le doute, si vous ne savez pas ce que vous manipulez avec cette fonction, utilisez le critère 'strain'.
Préparer à la visualisation :
Pour regarder le résultat voici deux option : restez en MATLAB comme vue à l'étape précédente ou allez dans VMD pour plus de pratique. Voyons comment faire pour sauvegarder le résultat en format pdb pour aller dans VMD :
1 2 3 4 5 6 7 8 9 10 |
[pts, Vb]= convhulln(XYZ); %astuce d’échelle scale=100/Vb^(1/3); XYZ=XYZ*scale ; data.X = XYZ(:,1); %on fait une structure de donnée que mat2pdb comprendra. data.Y = XYZ(:,2); data.Z = XYZ(:,3); data.outfile = 'chromosome19fast.pdb'; data.atomName = cell(1,length(XYZ)); data.atomName(1 :end) = {'CA'}; mat2pdb(data) |
Et là.. TADAAAA, on peut regarder ! Quel est l'intérêt ? Cette transformation est une bonne méthode de visualisation des données pour comprendre ce qui se passe, et comparer deux cartes de contact entre elles. Avoir un autre angle sur ces données, c'est toujours utile.
Dans cet article, j'ai beaucoup parlé de maths complexes sans les évoquer, en les montrant simplement comme des outils. Si vous souhaitez comprendre un peu plus le fond derrière, je ne peux que vous recommander de lire les articles associés évoqués plus haut (surtout les matériels supplémentaires en fait). Si l'article et le sujet vous intéressent, n'hésitez pas à me solliciter en commentaire, je me ferai une joie de partager d'autres articles sur ce thème !
Merci aux relecteurs et admin de la semaine : Yoann M., Guillaume Devailly, Julien Mozziconacci
Laisser un commentaire