In [1]:
import numpy as np
import matplotlib.pylab as plt

%matplotlib inline

Classification de gènes¶

On se propose d'éffectuer une classification de l'expression de gènes basée sur la distribution de leurs nucléotides. On possède un fichier data.dat dont voici le contenu.

In [2]:
cat "data/adn.dat"

0,gene1,GCACGGCATGCTGGCCTGCGCTAGCGAAGCCGACGCACGTCCTCGGCTCCAGCCGAACGTTTGGTCGAGGGTGGAGCACCGCCGGCACCTATTACTCGAC
0,gene2,GTACGCTCCGGTCGGGCAGGCCGTGGAGGGTCACATGCTCCCGCGACCAGGGCACAGGCCCCCATCGCACCCGCCATGACCACCCAGTCCGCCACGTCTC
0,gene3,TGACTTTCATTGTATGTCTTTTCGTTCGGTTCGTCGCGTGTTCTCCGTTGCGGTGGCCCCTCCACTTCTCGGTGGGTGCGGGGGGGGTTTTGTAGTTGGG
0,gene4,AAGATCGAACCCCCTCTGCCCCGCCGCCCACCAAACCCACCCACCGAGCCCCCCCGACGAGCGCAGCCCCACCCAGCCCGCAGCGACACACGGAGCACTG
0,gene5,GACACGCCCCAAGTTCCCCCCGCACCTGCTGGACAGCCTACTCTGCCCGTCGGACCCTCGCCGATTGACCTCCAGACACACCACCCCCGACAACTCCTCT
0,gene6,CCTGTCCCCGGCAGCACCCAGCGGCCCCGCCCCCCCGCGCCCGCACCGCCGCGCAACAGCTAACACGCACCCCAACGAGAAAAGGCCGCGCCCGCGGCCT
0,gene7,GGCCTTTCCCAGGCCAGGCCGGCCGCCCCACGGCCATCCTTCGGGGCGCAGATGAGCGCCCTCGCCGCGGTCCCTGATGACCGCTCGCCAGGAGAATCGC
0,gene8,CCGTCCCATGGCCCCCCTAACGCCGGCCAACCCGCATACTCCCGTTTCTACACCCCGCGCCTCGCGCACGCGCGCCATCACCCCGACCCCGGCCAACCCG
0,gene9,CCACCATTGATCGGGCTGAGCGCTTGTCGTCGTCGTCTAGTGTATGCCGGCCAGCTCACCCCTCCATGGTTCGGGTGTGGGCCGCCGCGTCGGCACGTCC
0,gene10,GCACTAAACAACACATGCCCGTGGCGGAGCCACACCCGCCCTCAACCATAGCGCCCGCGATCCAAACACGCGTCAACCATGGACCCCGACCCCCTCCACC
0,gene11,TCCTTACCTGGAACGCTATTACCTCATTCCAAGTCCGCACCTACCACCACCAAAACCACTCCACCTCCGGTGCGCATCCACCCCAACCACCATCCAAACT
0,gene12,CACCCACCGTCCAGCAGCGACCCCCATTGCCACCCACACCCTACCGGCCACCCCGTGCAACGACCTCCCACACACGAGACACCGCTGCAGACACGAGGCC
0,gene13,CAGCCCCGACGAAATATGCGCGCTCTGCACAGACACCTCCCCCTCTCTCTACCACAGCTATTCCCTACCACTGTACCCCCCAACCTGAAGGTCAGCGCTT
0,gene14,TACGCTATCCCTACTTTATGCTCTACATTATGTAACTTTTCTCTTTCTTTCACTGTCATTATCCCTCTAGTTCCCAGTCATTTACAACTTTCACTTTCTT
0,gene15,GTCCGCCGCCTTCTTGTCGGGTTTCTCTTGGTTTGGCGTCCGCTCGTCTGTATGCGCACCAGGGGTCTTTGTGGTGTGTGTCTTCATTGTCCGGTGGTGG
1,geneA,ACCAATTCCCTATTTAACCCTCTTCTACGGCGTACGATCCAATAGAGTAGGCGCCATGATGATCCACAAGACAAGGTAATAAACACACAGTTTACTCATG
1,geneB,TTAGAAAGATTGTATTCTCCTTCTTGGGAATATCGCGGTTTAATTGCATCGCATAGCTTCTAGCGGGAGGTGTTTTGTACGTATATCTCTTATTGTAATA
1,geneC,ACACAACCAAGCCACCGAAAGGCGAGAAGCTAGAAACCCGACCCGGATTCACCGGGAAAAGGAAACGCCACAACGATCACAACGGACATAGTACAATAAT
1,geneD,CAACCATCCATTATTTTACCAACTACCGTTCTACCATTTAGAATTCCTCGTACAATAGTATACTATTACCGACCACAGTGGTAACTTCATCGTACGTTCC
1,geneE,TACTTCTAATGGGGGAACGGAGCCAGGTCTTACTCCTTGTATCAGCGCAACGTATGCGATCGACGCGAAGTAAAACTCATATAACCAAGGTTAGAGGTAA
1,geneF,CCAGTTTGCGAGATAAGCAAATTCGACGAGGGTGGGTGGGTCTCGGTGAAATTGTCCGGTTAGGGGAGCGGACCGAGTGATAGCCCGAGTTGTGGGGGCC
1,geneG,TGTATGAGGTGGGTTGTTGGTTGAGGTGACGCTTGTGGTGGTATTTGGGACGGGGCTTGGTTTTTGGTGTGTTGTAAGGCTGAGTATGGTGGGGGGGAGT
1,geneH,GCGATCATATGTGCGATTGTGTTAAAGAGCGGAGGGGGGGTGGGCGTGTTGTTCTTGTTTATCTCTGTGGTTGGTGTTCAATGGAGTGGAGGAGGGTGGG
1,geneI,TTATCGTTAGTATCTTGGGAGGTTGTGTCTCCATAGCATGGACTGCAGTGATCATTAATTTCTTTGGTGCCTCTTTGTGGGTGTAGAGGGGAGGGCAGTC
1,geneJ,ACATACCATAACCATCGTACCGAACTTTTCTATAGTCTTAAAAAATAACAATTTACGTTCTCACAATATAGTAAGAATTCATCCTTAATCCTCATTTAAT
1,geneK,CTTCTAGCGAGAGGGCGAGGAGGAGGTGTAGGCAGGAGGCGGGAGGGCCTGTAGGGAGCGGCGGGGTAGGGGCTTGGTGCGCGGACAAGCCGGGGGGGTG
1,geneL,TGGCTGCTTTGGTCAGGGGCGAGGGCATCACGGGTGGCTAGAGTCTAGGTCTATATAGGGGTAGCAGCGGAGAAAGGCGTGGAAACGGCTAGGCGCCCGG
1,geneM,GGTGCGTTCTTAATTGGTCTATGGGTTTTTGGGATGTGTTGTCGTTTCCGCTGGGCTGGGCTATATGGAGCTCGAGGCTGTTTGGGGTTGTTGGCGATTG
1,geneN,GGCGGCAGACGGAGAACTCGTTATACCGGGTGCTAGGAGTGGAAGTATAGAGGGTATTTGTTCCTATATTGGTTTGTAAAGATTGTTGCCTAGACATGAA
1,geneO,GTTGTGGAAGGATCGTAATCTTGCTATTACATTAGGTAGGGTGTTATAGCTTTTTATGCTTTTTTTGTGACTGTGGATTTGACAGTATTGATTTACTGCA

Chaque ligne est composée d'une valeur booléenne indiquant l'expression du gène, du nom du gène et de sa séquence d’ADN. Chaque champ est séparé par une virgule.

Chargement des données¶

Commençons dans un permier temps par charger le fichier. On souhaite écrire loadData(filename: string) -> list qui prend en argument un nom de fichier puis charge les séquences de ce fichier sous forme d'une paire de liste. La première liste est une liste de valeur 0, 1 désignant l'expression du gène et la seconde est la liste des séquences. Deux méthodes des chaînes de caractères nous seront utiles ici:

  • strip() qui permet de supprimer les caractères blancs de fin et de début de ligne.
  • split(c) qui permet de couper une chaîne de caractères selon le caractère c et retourne un liste.

Le principe algorithmique est simple:

expr = Liste vide
genes = Liste vide
Pour chaque ligne $l$ du fichier:
    Enlever le caractère de fin de ligne
    champs = la ligne coupée par virgule
    Ajouter champs[0] dans la liste expr
    Ajouter champs[2] dans la liste genes

Ce qui donne en python:

In [3]:
def loadData(filename):
    with open(filename, "r") as f:
        exprs = []
        genes = []
        for line in f:
            line = line.strip()
            fields = line.split(",")
            exprs.append(int(fields[0]))
            genes.append(fields[2])
    return exprs, genes

Et lorsqu'on l'exécute sur le fichier de données:

In [4]:
expr_list, genes_list = loadData("data/adn.dat")
print("La liste d'expression:")
print(expr_list)
print("La liste de genes: ")
print(genes_list)

La liste d'expression:
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
La liste de genes: 
['GCACGGCATGCTGGCCTGCGCTAGCGAAGCCGACGCACGTCCTCGGCTCCAGCCGAACGTTTGGTCGAGGGTGGAGCACCGCCGGCACCTATTACTCGAC', 'GTACGCTCCGGTCGGGCAGGCCGTGGAGGGTCACATGCTCCCGCGACCAGGGCACAGGCCCCCATCGCACCCGCCATGACCACCCAGTCCGCCACGTCTC', 'TGACTTTCATTGTATGTCTTTTCGTTCGGTTCGTCGCGTGTTCTCCGTTGCGGTGGCCCCTCCACTTCTCGGTGGGTGCGGGGGGGGTTTTGTAGTTGGG', 'AAGATCGAACCCCCTCTGCCCCGCCGCCCACCAAACCCACCCACCGAGCCCCCCCGACGAGCGCAGCCCCACCCAGCCCGCAGCGACACACGGAGCACTG', 'GACACGCCCCAAGTTCCCCCCGCACCTGCTGGACAGCCTACTCTGCCCGTCGGACCCTCGCCGATTGACCTCCAGACACACCACCCCCGACAACTCCTCT', 'CCTGTCCCCGGCAGCACCCAGCGGCCCCGCCCCCCCGCGCCCGCACCGCCGCGCAACAGCTAACACGCACCCCAACGAGAAAAGGCCGCGCCCGCGGCCT', 'GGCCTTTCCCAGGCCAGGCCGGCCGCCCCACGGCCATCCTTCGGGGCGCAGATGAGCGCCCTCGCCGCGGTCCCTGATGACCGCTCGCCAGGAGAATCGC', 'CCGTCCCATGGCCCCCCTAACGCCGGCCAACCCGCATACTCCCGTTTCTACACCCCGCGCCTCGCGCACGCGCGCCATCACCCCGACCCCGGCCAACCCG', 'CCACCATTGATCGGGCTGAGCGCTTGTCGTCGTCGTCTAGTGTATGCCGGCCAGCTCACCCCTCCATGGTTCGGGTGTGGGCCGCCGCGTCGGCACGTCC', 'GCACTAAACAACACATGCCCGTGGCGGAGCCACACCCGCCCTCAACCATAGCGCCCGCGATCCAAACACGCGTCAACCATGGACCCCGACCCCCTCCACC', 'TCCTTACCTGGAACGCTATTACCTCATTCCAAGTCCGCACCTACCACCACCAAAACCACTCCACCTCCGGTGCGCATCCACCCCAACCACCATCCAAACT', 'CACCCACCGTCCAGCAGCGACCCCCATTGCCACCCACACCCTACCGGCCACCCCGTGCAACGACCTCCCACACACGAGACACCGCTGCAGACACGAGGCC', 'CAGCCCCGACGAAATATGCGCGCTCTGCACAGACACCTCCCCCTCTCTCTACCACAGCTATTCCCTACCACTGTACCCCCCAACCTGAAGGTCAGCGCTT', 'TACGCTATCCCTACTTTATGCTCTACATTATGTAACTTTTCTCTTTCTTTCACTGTCATTATCCCTCTAGTTCCCAGTCATTTACAACTTTCACTTTCTT', 'GTCCGCCGCCTTCTTGTCGGGTTTCTCTTGGTTTGGCGTCCGCTCGTCTGTATGCGCACCAGGGGTCTTTGTGGTGTGTGTCTTCATTGTCCGGTGGTGG', 'ACCAATTCCCTATTTAACCCTCTTCTACGGCGTACGATCCAATAGAGTAGGCGCCATGATGATCCACAAGACAAGGTAATAAACACACAGTTTACTCATG', 'TTAGAAAGATTGTATTCTCCTTCTTGGGAATATCGCGGTTTAATTGCATCGCATAGCTTCTAGCGGGAGGTGTTTTGTACGTATATCTCTTATTGTAATA', 'ACACAACCAAGCCACCGAAAGGCGAGAAGCTAGAAACCCGACCCGGATTCACCGGGAAAAGGAAACGCCACAACGATCACAACGGACATAGTACAATAAT', 'CAACCATCCATTATTTTACCAACTACCGTTCTACCATTTAGAATTCCTCGTACAATAGTATACTATTACCGACCACAGTGGTAACTTCATCGTACGTTCC', 'TACTTCTAATGGGGGAACGGAGCCAGGTCTTACTCCTTGTATCAGCGCAACGTATGCGATCGACGCGAAGTAAAACTCATATAACCAAGGTTAGAGGTAA', 'CCAGTTTGCGAGATAAGCAAATTCGACGAGGGTGGGTGGGTCTCGGTGAAATTGTCCGGTTAGGGGAGCGGACCGAGTGATAGCCCGAGTTGTGGGGGCC', 'TGTATGAGGTGGGTTGTTGGTTGAGGTGACGCTTGTGGTGGTATTTGGGACGGGGCTTGGTTTTTGGTGTGTTGTAAGGCTGAGTATGGTGGGGGGGAGT', 'GCGATCATATGTGCGATTGTGTTAAAGAGCGGAGGGGGGGTGGGCGTGTTGTTCTTGTTTATCTCTGTGGTTGGTGTTCAATGGAGTGGAGGAGGGTGGG', 'TTATCGTTAGTATCTTGGGAGGTTGTGTCTCCATAGCATGGACTGCAGTGATCATTAATTTCTTTGGTGCCTCTTTGTGGGTGTAGAGGGGAGGGCAGTC', 'ACATACCATAACCATCGTACCGAACTTTTCTATAGTCTTAAAAAATAACAATTTACGTTCTCACAATATAGTAAGAATTCATCCTTAATCCTCATTTAAT', 'CTTCTAGCGAGAGGGCGAGGAGGAGGTGTAGGCAGGAGGCGGGAGGGCCTGTAGGGAGCGGCGGGGTAGGGGCTTGGTGCGCGGACAAGCCGGGGGGGTG', 'TGGCTGCTTTGGTCAGGGGCGAGGGCATCACGGGTGGCTAGAGTCTAGGTCTATATAGGGGTAGCAGCGGAGAAAGGCGTGGAAACGGCTAGGCGCCCGG', 'GGTGCGTTCTTAATTGGTCTATGGGTTTTTGGGATGTGTTGTCGTTTCCGCTGGGCTGGGCTATATGGAGCTCGAGGCTGTTTGGGGTTGTTGGCGATTG', 'GGCGGCAGACGGAGAACTCGTTATACCGGGTGCTAGGAGTGGAAGTATAGAGGGTATTTGTTCCTATATTGGTTTGTAAAGATTGTTGCCTAGACATGAA', 'GTTGTGGAAGGATCGTAATCTTGCTATTACATTAGGTAGGGTGTTATAGCTTTTTATGCTTTTTTTGTGACTGTGGATTTGACAGTATTGATTTACTGCA']

Calcul des histogrammes¶

Une fois le fichier chargé, nous allons calculer l'histogramme de chaque séquence (i.e. le nombre de "A", "T", "C", "G") et nous allons les stocker dans un tableau numpy pour pouvoir les manipuler plus facilement par la suite. On écrit la fonction histGenes(expr: List, genes: List) -> ndarray qui prend la liste d'expression des gènes et la liste des séquences et renvoie une matrice $n \times 5$. La première colonne correspond à l'expression du gène (0 ou 1), les suivantes correspondent au nombre de "A", "T", "C", "G".

On utilise ici la méthode count(c) des chaînes de caractères qui permet de compter le nombre d'occurence du caractère c.

Le principe algorithmique, on suppose que $n$ est le nombre de gènes.

  1. Créer une matrice vide $M$ de taille $n \times 5$
  2. Pour chaque gène n°$i$, d'expression $e$ et de séquence $s$:
    • La $i$ ème ligne de $M$ se compose de [$e$, #A dans $s$, #T dans $s$, #C dans $s$, #G dans $s$]
In [5]:
def histGenes(expr, genes):
    M = []
    for e,g in zip(expr, genes):
        na, nt, nc, ng = g.count("A"), g.count("T"), g.count("C"), g.count("G")
        M.append([e, na, nt, nc, ng])
    return np.array(M)

En exécutant cette fonction sur nos deux listes:

In [6]:
H = histGenes(expr_list, genes_list)
H
Out[6]:
array([[ 0, 17, 16, 35, 32],
       [ 0, 16, 12, 43, 29],
       [ 0,  5, 38, 22, 35],
       [ 0, 23,  4, 53, 20],
       [ 0, 19, 15, 49, 17],
       [ 0, 18,  4, 53, 25],
       [ 0, 13, 13, 42, 32],
       [ 0, 15, 11, 55, 19],
       [ 0, 10, 23, 35, 32],
       [ 0, 26,  9, 47, 18],
       [ 0, 27, 18, 46,  9],
       [ 0, 24,  7, 51, 18],
       [ 0, 22, 19, 44, 15],
       [ 0, 19, 46, 29,  6],
       [ 0,  4, 37, 25, 34],
       [ 1, 33, 24, 27, 16],
       [ 1, 23, 40, 15, 22],
       [ 1, 42,  8, 29, 21],
       [ 1, 29, 32, 29, 10],
       [ 1, 31, 23, 21, 25],
       [ 1, 20, 21, 18, 41],
       [ 1, 11, 36,  5, 48],
       [ 1, 15, 32,  9, 44],
       [ 1, 17, 36, 15, 32],
       [ 1, 38, 34, 22,  6],
       [ 1, 16, 12, 16, 56],
       [ 1, 20, 18, 19, 43],
       [ 1,  9, 39, 13, 39],
       [ 1, 26, 29, 13, 32],
       [ 1, 21, 44, 10, 25]])

Normalisation des histogrammes¶

Les histogrammes ne sont pas normalisés (la somme ne fait pas 1). Nous allons les normaliser afin que ce soit de vraies distributions. La fonction normHist(H: ndarray) -> ndarray renvoie une copie de H normalisé.

On utilisera les fonctions suivantes:

  • copy(A) qui permet de créer une copie d'un tableau.
  • sum(A) qui fait la somme d'un tableau.

Note: ce sont aussi des méthodes. On peut écrire sum(A) ou A.sum().

Pour normaliser un histogramme, il suffit de sommer chaque /bin/ de l'histogramme puis de diviser chaque /bin/ par cette somme.

Le slicing est trés utile ici pour récupérer une ligne entière. Par exemple, pour récupérer la 10ème ligne:

In [7]:
H[9,:]
Out[7]:
array([ 0, 26,  9, 47, 18])

Comme on s'intéresse uniquement à l'histogramme, la première colonne est inutile. On peut donc slicer à partir du 1er élément, sommer, puis normaliser:

In [8]:
print("10e histogramme:", H[9,1:])
print("10e histogramme (somme):", H[9,1:].sum())
print("10e histogramme (normalisé):", H[9,1:] / H[9,1:].sum())

10e histogramme: [26  9 47 18]
10e histogramme (somme): 100
10e histogramme (normalisé): [0.26 0.09 0.47 0.18]

Faisons une boucle sur les lignes et appliquons ce que l'on vient de faire pour chaque histogramme.

In [9]:
def normHist(A):
    B = A.astype(float).copy()
    n,m = A.shape
    for i in range(n):
        B[i,1:] = A[i,1:] / A[i,1:].sum()
    return B

En l'appliquant sur notre matrice d'histogramme:

In [10]:
h = normHist(H)
h
Out[10]:
array([[0.  , 0.17, 0.16, 0.35, 0.32],
       [0.  , 0.16, 0.12, 0.43, 0.29],
       [0.  , 0.05, 0.38, 0.22, 0.35],
       [0.  , 0.23, 0.04, 0.53, 0.2 ],
       [0.  , 0.19, 0.15, 0.49, 0.17],
       [0.  , 0.18, 0.04, 0.53, 0.25],
       [0.  , 0.13, 0.13, 0.42, 0.32],
       [0.  , 0.15, 0.11, 0.55, 0.19],
       [0.  , 0.1 , 0.23, 0.35, 0.32],
       [0.  , 0.26, 0.09, 0.47, 0.18],
       [0.  , 0.27, 0.18, 0.46, 0.09],
       [0.  , 0.24, 0.07, 0.51, 0.18],
       [0.  , 0.22, 0.19, 0.44, 0.15],
       [0.  , 0.19, 0.46, 0.29, 0.06],
       [0.  , 0.04, 0.37, 0.25, 0.34],
       [1.  , 0.33, 0.24, 0.27, 0.16],
       [1.  , 0.23, 0.4 , 0.15, 0.22],
       [1.  , 0.42, 0.08, 0.29, 0.21],
       [1.  , 0.29, 0.32, 0.29, 0.1 ],
       [1.  , 0.31, 0.23, 0.21, 0.25],
       [1.  , 0.2 , 0.21, 0.18, 0.41],
       [1.  , 0.11, 0.36, 0.05, 0.48],
       [1.  , 0.15, 0.32, 0.09, 0.44],
       [1.  , 0.17, 0.36, 0.15, 0.32],
       [1.  , 0.38, 0.34, 0.22, 0.06],
       [1.  , 0.16, 0.12, 0.16, 0.56],
       [1.  , 0.2 , 0.18, 0.19, 0.43],
       [1.  , 0.09, 0.39, 0.13, 0.39],
       [1.  , 0.26, 0.29, 0.13, 0.32],
       [1.  , 0.21, 0.44, 0.1 , 0.25]])

Affichage des hisogrammes¶

Nous allons maintenant en profiter pour afficher la distribution moyenne et l'écart type pour les gènes exprimés puis pour les gènes non-exprimés.

D'abord filtrons les données pour récupérer uniquement les gènes non-exprimés (ceux dont la 1ère colonne est à 0). Pour récupérer cette sous-matrice on peut alors écrire:

In [11]:
h0 = h[h[:,0] == 0]
h0
Out[11]:
array([[0.  , 0.17, 0.16, 0.35, 0.32],
       [0.  , 0.16, 0.12, 0.43, 0.29],
       [0.  , 0.05, 0.38, 0.22, 0.35],
       [0.  , 0.23, 0.04, 0.53, 0.2 ],
       [0.  , 0.19, 0.15, 0.49, 0.17],
       [0.  , 0.18, 0.04, 0.53, 0.25],
       [0.  , 0.13, 0.13, 0.42, 0.32],
       [0.  , 0.15, 0.11, 0.55, 0.19],
       [0.  , 0.1 , 0.23, 0.35, 0.32],
       [0.  , 0.26, 0.09, 0.47, 0.18],
       [0.  , 0.27, 0.18, 0.46, 0.09],
       [0.  , 0.24, 0.07, 0.51, 0.18],
       [0.  , 0.22, 0.19, 0.44, 0.15],
       [0.  , 0.19, 0.46, 0.29, 0.06],
       [0.  , 0.04, 0.37, 0.25, 0.34]])

Calculons ensuite la distribution moyenne. La méthode mean permet de calculer la moyenne du tableau entier, ce n'est pas ce que nous voulons. Cette méthode possède un argument axis qui permet d'effectuer l'opération dans une direction. On veut éffectuer la moyenne par colonne (axe 0), on écrira donc:

In [12]:
vmean = h0[:,1:].mean(axis=0)
vmean
Out[12]:
array([0.172     , 0.18133333, 0.41933333, 0.22733333])

De même pour l'écart type:

In [13]:
vstd = h0[:,1:].std(axis=0)
vstd
Out[13]:
array([0.06744875, 0.12349719, 0.10142101, 0.08970074])

On finit par afficher l'histogramme moyen avec bar. La signature est bar(left, height, width=0.8) où left est un tableau des coordonnées $x$ du début des barres (on a 4 colonnes donc width = [0,1,2,3]), height est la hauteur des barres et width la largeur. Les options yerr et ecolor permettent de mettre des barres d'erreurs et de leur donner une couleur. Enfin xticks(postion, labels) permet de placer des labels sur l'axe $x$ aux positions voulues.

In [14]:
plt.bar(range(4), vmean, yerr = vstd, ecolor="r")
plt.title(u"Distribution Marginale des nucléotides pour les gènes non exprimés")
plt.xticks((0, 1, 2, 3)  , ("A","T","C","G"));
No description has been provided for this image

En faisant de même, pour les gènes exprimés et en s'inspirant-vous du code http://matplotlib.org/examples/api/barchart_demo.html, on peut créer un joli graphique avec les histogrammes côte-à-côte.

In [15]:
h0, h1 = h[h[:,0] == 0], h[h[:,0] == 1]
vmean0 = h0[:,1:5].mean(axis=0)
vmean1 = h1[:,1:5].mean(axis=0)
vstd0 = h0[:,1:5].std(axis=0)
vstd1 = h1[:,1:5].std(axis=0)
plt.bar(np.arange(4), vmean0, width=0.4, yerr=vstd0, ecolor="r", label="Gènes non-exprimés")
plt.bar(np.arange(4) + 0.4, vmean1, width=0.4, color='g', yerr=vstd1, ecolor="r", label="Gènes exprimés")
plt.legend()
plt.xticks(np.arange(4) + 0.2, ("A","T","C","G"));
No description has been provided for this image

Il y a un gap important dans les distibutions du C, il doit sûrement jouer un rôle dans la discrimination.

Affichage des gènes en fonction de la distribution des nucléotides¶

On souhaite écrire la fonction plotGenes(H: numpy array, nuca: string, nucb: string) où

  • $nuc_a$ et $nuc_b$ sont deux nucléotides
  • $H$ est la matrice d'histogrammes des gènes.

On note $n_{nuc_a}(g)$ qui est le nombre de nucléotides $nuc_a$ dans le gène $g$. Cette fonction crée un graphique où:

  • pour chaque gène $g$ qui s'exprime, un point bleu est affiché à la coordonnée $(n_{nuc_a}(g), n_{nuc_b}(g))$
  • pour chaque gène $h$ qui ne s'exprime pas, un carré rouge est affiché à la coordonnée $(n_{nuc_a}(h), n_{nuc_b}(h))$

On utilise pour cela les fonctions $scatter$ qui permettent d'afficher des points, ainsi:

  • scatter(x, y, c="b", marker="o") permet d'afficher une séquence de points de coordonnées $x$ et $y$ sous forme de points bleus
  • scatter(x, y, c="r", marker="s") permet d'afficher une séquence de points de coordonnées $x$ et $y$ sous forme de carrés rouges.
In [16]:
def plotGenes(H, a, b):
    idx = { "A" : 1, "T" : 2, "C" : 3, "G" : 4}
    i,j = idx[a], idx[b]
    matA = H[H[:,0] == 0]
    matB = H[H[:,0] == 1]
    plt.scatter(matA[:,i], matA[:,j], c="b", marker="o")
    plt.scatter(matB[:,i], matB[:,j], c="r", marker="s")
    plt.xlabel("Pourcentage de " + a)
    plt.ylabel("Pourcentage de " + b)

On affiche ensuite pour chaque couple de nucléotides $(a,b)$, le graphique qui correspond à plotGenes(H, a, b). On utilisera la fonction subplot qui permet de créer des sous graphiques.

In [17]:
plt.figure(figsize=(15,5))
i=1
for a in "ACGT":
    for b in "ACGT":
        if b >= a:
            plt.subplot(2,5,i)
            plotGenes(H, a, b)
            i+=1
No description has been provided for this image

Lors de l'affichage des histogrammes, on avait déjà perçu le fait que le nombre de Cytosine pourait être discriminant, on le voit de nouveau sur le graphique 5 de la 1ere ligne. Au deça de 30% environ de Cytosines dans la séquence, le gène est non exprimé, en dessous le gène est exprimé. Cette manière de discriminer n'est cependant pas parfaite, et crée un taux d'erreurs non nul, puisqu'on voit clairement que ces deux ensembles ne sont pas linéairement séparable sur le graphique 5.

Le seul graphique où les deux ensembles sont linéairement séparables est le second graphique de la première. Approximativement, on peut tracer une droite d'équation $p_C = \frac{2}{3}.p_A + \frac{35}{3}$, où $p_C$ est la proportion de Cytosines et $p_A$ la proportion d'Adénine. En d'autres termes:

  • si $p_C > \frac{2}{3}.p_A + \frac{35}{3}$, le gène est non exprimé.
  • si $p_C < \frac{2}{3}.p_A + \frac{35}{3}$, le gène est exprimé.

On peut tracer cette droite.

In [18]:
plotGenes(H, "A", "C")
x=np.array([0, 45])
y=2/3 * x + 35/3
plt.plot(x, y)
plt.xlim((0,45))
Out[18]:
(0.0, 45.0)
No description has been provided for this image

Vers un classifieur¶

A partir de l'observation précédente, il est maintenant trés facile d'en déduire un classifieur qui à partir de la séquence d'un gène, va dire s'il sera exprimé ou pas. Pour cela, on a juste à calculer la proportion de Cytosines $pC$ et d'Adénines $pA$ puis:

  • si $p_C > \frac{2}{3}.p_A + \frac{35}{3}$, le gène est non exprimé.
  • si $p_C < \frac{2}{3}.p_A + \frac{35}{3}$, le gène est exprimé.
In [19]:
def discriminate(seq):
    n = float(len(seq))
    na = seq.count("A") / n * 100
    nc = seq.count("C") / n * 100
    a=2.0/3
    b=35.0/3
    if nc > (a * na + b):
        return "Non Exprimé"
    else:
        return "Exprimé"

Testons maintenant sur nos gènes de départ: (on rapelle que zip permet de parcourir deux séquences en même temps)

In [20]:
for expr, gene in zip(expr_list, genes_list):
    print(f"Le gene {gene} \ndont l'expression est {expr} a été détecté : {discriminate(gene)}")

Le gene GCACGGCATGCTGGCCTGCGCTAGCGAAGCCGACGCACGTCCTCGGCTCCAGCCGAACGTTTGGTCGAGGGTGGAGCACCGCCGGCACCTATTACTCGAC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene GTACGCTCCGGTCGGGCAGGCCGTGGAGGGTCACATGCTCCCGCGACCAGGGCACAGGCCCCCATCGCACCCGCCATGACCACCCAGTCCGCCACGTCTC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene TGACTTTCATTGTATGTCTTTTCGTTCGGTTCGTCGCGTGTTCTCCGTTGCGGTGGCCCCTCCACTTCTCGGTGGGTGCGGGGGGGGTTTTGTAGTTGGG 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene AAGATCGAACCCCCTCTGCCCCGCCGCCCACCAAACCCACCCACCGAGCCCCCCCGACGAGCGCAGCCCCACCCAGCCCGCAGCGACACACGGAGCACTG 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene GACACGCCCCAAGTTCCCCCCGCACCTGCTGGACAGCCTACTCTGCCCGTCGGACCCTCGCCGATTGACCTCCAGACACACCACCCCCGACAACTCCTCT 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene CCTGTCCCCGGCAGCACCCAGCGGCCCCGCCCCCCCGCGCCCGCACCGCCGCGCAACAGCTAACACGCACCCCAACGAGAAAAGGCCGCGCCCGCGGCCT 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene GGCCTTTCCCAGGCCAGGCCGGCCGCCCCACGGCCATCCTTCGGGGCGCAGATGAGCGCCCTCGCCGCGGTCCCTGATGACCGCTCGCCAGGAGAATCGC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene CCGTCCCATGGCCCCCCTAACGCCGGCCAACCCGCATACTCCCGTTTCTACACCCCGCGCCTCGCGCACGCGCGCCATCACCCCGACCCCGGCCAACCCG 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene CCACCATTGATCGGGCTGAGCGCTTGTCGTCGTCGTCTAGTGTATGCCGGCCAGCTCACCCCTCCATGGTTCGGGTGTGGGCCGCCGCGTCGGCACGTCC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene GCACTAAACAACACATGCCCGTGGCGGAGCCACACCCGCCCTCAACCATAGCGCCCGCGATCCAAACACGCGTCAACCATGGACCCCGACCCCCTCCACC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene TCCTTACCTGGAACGCTATTACCTCATTCCAAGTCCGCACCTACCACCACCAAAACCACTCCACCTCCGGTGCGCATCCACCCCAACCACCATCCAAACT 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene CACCCACCGTCCAGCAGCGACCCCCATTGCCACCCACACCCTACCGGCCACCCCGTGCAACGACCTCCCACACACGAGACACCGCTGCAGACACGAGGCC 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene CAGCCCCGACGAAATATGCGCGCTCTGCACAGACACCTCCCCCTCTCTCTACCACAGCTATTCCCTACCACTGTACCCCCCAACCTGAAGGTCAGCGCTT 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene TACGCTATCCCTACTTTATGCTCTACATTATGTAACTTTTCTCTTTCTTTCACTGTCATTATCCCTCTAGTTCCCAGTCATTTACAACTTTCACTTTCTT 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene GTCCGCCGCCTTCTTGTCGGGTTTCTCTTGGTTTGGCGTCCGCTCGTCTGTATGCGCACCAGGGGTCTTTGTGGTGTGTGTCTTCATTGTCCGGTGGTGG 
dont l'expression est 0 a été détecté : Non Exprimé
Le gene ACCAATTCCCTATTTAACCCTCTTCTACGGCGTACGATCCAATAGAGTAGGCGCCATGATGATCCACAAGACAAGGTAATAAACACACAGTTTACTCATG 
dont l'expression est 1 a été détecté : Exprimé
Le gene TTAGAAAGATTGTATTCTCCTTCTTGGGAATATCGCGGTTTAATTGCATCGCATAGCTTCTAGCGGGAGGTGTTTTGTACGTATATCTCTTATTGTAATA 
dont l'expression est 1 a été détecté : Exprimé
Le gene ACACAACCAAGCCACCGAAAGGCGAGAAGCTAGAAACCCGACCCGGATTCACCGGGAAAAGGAAACGCCACAACGATCACAACGGACATAGTACAATAAT 
dont l'expression est 1 a été détecté : Exprimé
Le gene CAACCATCCATTATTTTACCAACTACCGTTCTACCATTTAGAATTCCTCGTACAATAGTATACTATTACCGACCACAGTGGTAACTTCATCGTACGTTCC 
dont l'expression est 1 a été détecté : Exprimé
Le gene TACTTCTAATGGGGGAACGGAGCCAGGTCTTACTCCTTGTATCAGCGCAACGTATGCGATCGACGCGAAGTAAAACTCATATAACCAAGGTTAGAGGTAA 
dont l'expression est 1 a été détecté : Exprimé
Le gene CCAGTTTGCGAGATAAGCAAATTCGACGAGGGTGGGTGGGTCTCGGTGAAATTGTCCGGTTAGGGGAGCGGACCGAGTGATAGCCCGAGTTGTGGGGGCC 
dont l'expression est 1 a été détecté : Exprimé
Le gene TGTATGAGGTGGGTTGTTGGTTGAGGTGACGCTTGTGGTGGTATTTGGGACGGGGCTTGGTTTTTGGTGTGTTGTAAGGCTGAGTATGGTGGGGGGGAGT 
dont l'expression est 1 a été détecté : Exprimé
Le gene GCGATCATATGTGCGATTGTGTTAAAGAGCGGAGGGGGGGTGGGCGTGTTGTTCTTGTTTATCTCTGTGGTTGGTGTTCAATGGAGTGGAGGAGGGTGGG 
dont l'expression est 1 a été détecté : Exprimé
Le gene TTATCGTTAGTATCTTGGGAGGTTGTGTCTCCATAGCATGGACTGCAGTGATCATTAATTTCTTTGGTGCCTCTTTGTGGGTGTAGAGGGGAGGGCAGTC 
dont l'expression est 1 a été détecté : Exprimé
Le gene ACATACCATAACCATCGTACCGAACTTTTCTATAGTCTTAAAAAATAACAATTTACGTTCTCACAATATAGTAAGAATTCATCCTTAATCCTCATTTAAT 
dont l'expression est 1 a été détecté : Exprimé
Le gene CTTCTAGCGAGAGGGCGAGGAGGAGGTGTAGGCAGGAGGCGGGAGGGCCTGTAGGGAGCGGCGGGGTAGGGGCTTGGTGCGCGGACAAGCCGGGGGGGTG 
dont l'expression est 1 a été détecté : Exprimé
Le gene TGGCTGCTTTGGTCAGGGGCGAGGGCATCACGGGTGGCTAGAGTCTAGGTCTATATAGGGGTAGCAGCGGAGAAAGGCGTGGAAACGGCTAGGCGCCCGG 
dont l'expression est 1 a été détecté : Exprimé
Le gene GGTGCGTTCTTAATTGGTCTATGGGTTTTTGGGATGTGTTGTCGTTTCCGCTGGGCTGGGCTATATGGAGCTCGAGGCTGTTTGGGGTTGTTGGCGATTG 
dont l'expression est 1 a été détecté : Exprimé
Le gene GGCGGCAGACGGAGAACTCGTTATACCGGGTGCTAGGAGTGGAAGTATAGAGGGTATTTGTTCCTATATTGGTTTGTAAAGATTGTTGCCTAGACATGAA 
dont l'expression est 1 a été détecté : Exprimé
Le gene GTTGTGGAAGGATCGTAATCTTGCTATTACATTAGGTAGGGTGTTATAGCTTTTTATGCTTTTTTTGTGACTGTGGATTTGACAGTATTGATTTACTGCA 
dont l'expression est 1 a été détecté : Exprimé

Auteurs : cette feuille a été rédigée par Edwin Carlinet et mise en Python moderne par Olivier Ricou.

In [ ]: