Introduction

Le cours d'aujourd'hui reprend, avec des détails supplémentaires, un cours d'« étude de cas » que j'ai donné à des M1 d'ingénierie mathématique (où les codes sont en R au lieu de Python).

Nous allons nous occuper aujourd'hui avec de vraies données recueillies dans une étude de la dynamique des ions calcium $Ca^{2+}$ dans des neurones de l'hypothalamus de souris. Ces neurones libèrent des hormones peptidiques dérivées de la Pro-opiomélanocortine (POMC) dont certains sont fortement impliqués dans la régulation de l'appétit. Comme les problèmes d'obésité deviennent préoccupant dans nos sociétés, les physiologistes s'intéressent de plus en plus à ces problèmes -- comprenez qu'ils peuvent obtenir de l'argent (essentiellement de sources publics) pour faire « tourner » leurs labos quand ils travaillent sur ces problèmes.

Mes collègues du laboratoire de Peter Kloppenburg à l'Université de Cologne (Allemagne) travaillent sur ces questions en collaboration avec les généticiens du laboratoire de Jens Brüning (même université). Ils étudient l'effet du régime alimentaire imposé aux souris (basses calories, niveau normal, hautes calories) sur l'activité des neurones synthétisant la POMC. Plus particulièrement ils étudient l'effet du régime sur la dynamique des ions calcium -- que nous savons par ailleurs avoir une influence prépondérante sur l'activité de ces neurones. Les données que nous allons utiliser ici ont été enregistrées par Andreas Pippow lors de sa thèse dans le labo de Peter Kloppenburg. Nous allons utilser deux jeux de données :

Si vous souhaitez en savoir plus sur les raisons qui poussent les physiologistes à étudier la dynamique des ions calcium en général, consulter mon cours sur le sujet.

Chargement des données « POMC » dans Python

Nous procédons de la « manière habituelle » après avoir démarré IPython :

In [1]:
%pylab
from __future__ import print_function, division, unicode_literals, absolute_import
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

Si nous souhaitons avoir les figures dans notre notebook IPython nous utilisons aussi (avec une commande additionnelle pour améliorer la résolution des graphiques que j'ai trouvée sur un excellent site) :

In [2]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'

Nous allons maintenant voir successivement comment charger les données dans les trois formats disponibles.

Lecture des données au format text

In [3]:
from urllib2 import urlopen
u = urlopen("http://xtof.disque.math.cnrs.fr/data/POMC.txt")

Les 6 premières lignes du fichier contiennent une description des données :

In [4]:
u.readline()
Out[4]:
'POMC data set recorded by Andreas Pippow (Kloppenburg Laboratory Cologne University, http://cecad.uni-koeln.de/Prof-Peter-Kloppenburg.82.0.html)\n'
In [5]:
u.readline()
Out[5]:
'168 measurements performed with a CCD camera recording Fura-2 fluorescence (excitation wavelength: 340 nm).\n'
In [6]:
u.readline()
Out[6]:
'The size of the CCD chip is 60 x 80 pixels. A stimulation (depolarization induced calcium entry) comes at time 527.\n'
In [7]:
u.readline()
Out[7]:
'Details about this data set can be found in: Joucla et al (2013) Estimating background-subtracted fluorescence transients in calcium imaging experiments: A quantitative approach. Cell Calcium. In Press; in the CRAN package FluOMatic (http://cran.at.r-project.org/web/packages/FluOMatic/index.html).\n'
In [8]:
u.readline()
Out[8]:
'The first 168 numbers are the time of measurement (in s). The next numbers are the ADU readings frame by frame the pixel matrix in row major mode; that is, the first 60 numbers are the sixty readings of the first coumn at the first time, etc.\n'
In [9]:
u.readline()
Out[9]:
'To read these data into R (assuming this file seats in your working directory) do: {toto <- scan("POMC.txt",skip=6);time <- toto[1:168];stack <- array(toto[-(1:168)],dim=c(60,80,168));rm(toto)}\n'

Les 168 lignes suivantes contiennent le temps (en secondes) auquel chacune des 168 expositions a été effectuée. Nous stockons ces temps dans un vecteur :

In [10]:
temps_POMC = np.zeros((168,))
for i in range(168):
    temps_POMC[i] = float(u.readline())

Le capteur CCD est une matrice de 60 x 80 pixels. Les 60 x 80 x 168 lignes suivantes contiennent les lectures en sortie de capteur de pour chaque pixel lors de chacune des 168 expositions. Nous les stockons dans un array à trois dimensions :

In [11]:
pile_POMC = np.zeros((60,80,168),dtype=int)
for k in range(168):
    for j in range(80):
        for i in range(60):
            pile_POMC[i,j,k] = int(u.readline())

Une fois la lecture finie, nous fermons le « fichier » :

In [12]:
u.close()

À ce stade, c'est une bonne idée de regarder les données pour s'assurer qu'aucune grosse bêtise n'a été commise lors de leur chargement. Une façon rapide de faire est de calculer la moyenne (au cours du temps) du signal pour chaque pixel et de regarder ce que ça donne.

La moyenne par pixel s'obtient simplement avec la méthode mean :

In [13]:
pile_POMC_moyenne = pile_POMC.mean(2)

On visualise alors avec la fonction imshow de pyplot :

In [14]:
imshow(transpose(pile_POMC_moyenne),origin='lower')
set_cmap('gray')
colorbar()
Out[14]:
<matplotlib.colorbar.Colorbar instance at 0x7f176fbe7320>

Lecture des données au format npz

La version npz des données créée avec la fonction savez_compressed de numpy doit d'abord être téléchargée, ce qui se fait simplement avec :

In [15]:
import urllib
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/Data_POMC.npz","Data_POMC.npz")
Out[15]:
(u'Data_POMC.npz', <httplib.HTTPMessage instance at 0x7f3d933a6368>)

Le fichier Data_POMC.npz contient deux arrays, le premier avec le temps des expositions, le second avec la pile d'images. On charge les données de la façon suivante :

In [15]:
POMC = load("Data_POMC.npz")
temps_POMC_npz = POMC['arr_0']
pile_POMC_npz = POMC['arr_1']
POMC.close()

Ici une rapide manipulation nous permet de vérifier que nous avons chargé les mêmes données :

In [16]:
pile_POMC_npz_moyenne = pile_POMC_npz.mean(2)
imshow(transpose(pile_POMC_npz_moyenne),origin='lower')
colorbar()
Out[16]:
<matplotlib.colorbar.Colorbar instance at 0x7f176d820680>

Lecture des données au format hdf5

Une troisième version des données est disponible en format HDF5 ; un format de stockage hiérarchique -- l'organisation interne du fichier est similaire à celle des fichiers et répertoires sur votre ordinateur -- qui permet d'annoter les jeux de données -- il est ainsi possible de sauvegarder les informations relatives à l'acquisition des données avec les données elles-mêmes. Un module Python, le module h5py permet de manipuler ces fichiers sans peine.

Nous procédons comme dans le cas précédent en téléchargeant les données :

In [18]:
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/Data_POMC.hdf5","Data_POMC.hdf5")
Out[18]:
(u'Data_POMC.hdf5', <httplib.HTTPMessage instance at 0x7f3d912db4d0>)

Nous continuons en chargeant le module puis nous ouvrons le fichier en mode lecture :

In [17]:
import h5py
pomc = h5py.File("Data_POMC.hdf5","r")

Pour obtenir la liste des jeux de données dataset et des groups (équivalents aux répertoires d'un système de fichier), on emploie la fonction list :

In [18]:
list(pomc)
Out[18]:
[u'stack', u'time']

Nous retrouvons le vecteur des temps, time, et la pile d'images, stack. Pour extraire ces deux jeux de données nous employons :

In [19]:
temps_POMC_hdf5 = pomc['time'][...]
pile_POMC_hdf5 = pomc['stack'][...]

Une fois les données extraites, nous fermons le fichier :

In [20]:
pomc.close()

Comme précédemment, nous vérifions que l'importation s'est « bien » passée avec :

In [21]:
pile_POMC_hdf5_moyenne = pile_POMC_hdf5.mean(2)
imshow(transpose(pile_POMC_hdf5_moyenne),origin='lower')
colorbar()
Out[21]:
<matplotlib.colorbar.Colorbar instance at 0x7f1765db2f80>

Représentation des données

Le plus souvent, les physiologistes vont montrer des piles d'images sous forme de films. C'est possible de le faire dans Python et si vous souhaitez le faire vous trouverez sur le oueb les indications nécessaires. Je préfère néanmoins une représentation qui montre sous forme matricielle, les éléments de la matrice correspondant aux pixels du capteur CCD, le décours temporel de la fluorescence. Pour obtenir ce type de représentation, nous allons définir une fonction :

In [22]:
def plotSignal(stack,lw=1):
    n_x, n_y, n_t = stack.shape
    amp_min = np.min(stack)
    amp_max = np.max(stack)
    amp_diff = np.ptp(stack)
    x_domain = np.arange(n_t)/n_t
    y_domain = (0,n_y)
    for r_idx in range(n_x):
        for c_idx in range(n_y):
            y_min = n_x - r_idx - 1
            sig = stack[r_idx,c_idx,:]
            Y = (sig-amp_min)/amp_diff + y_min
            X = x_domain + c_idx
            plt.plot(X,Y,lw=lw,color='black')
    plt.axis('off')

Nous appelons notre fonction sur notre pile en regardant ce qui se passe du coté du soma :

In [23]:
plotSignal(pile_POMC[20:41,30:51,:],lw=0.5)