Vous allez utiliser des outils informatiques qui vous permettront d’analyser des données brutes issues d’une analyse Shotgun Proteomics récemment publiée dans le journal Science sous le titre "Real-time visualization of drug resistance acquisition by horizontal gene transfer reveals a new role for AcrAB-TolC multidrug efflux pump".
Les données associées à cette publication sont publiques et accessibles sur la plateforme PRIDE. Le PDF de la publication est data/Nolivos_2019.pdf.
Vous forkerez le présent "repository" pour vous permettre de sauvegarder votre travail.
Vous le clonerez ensuite dans votre espace de travail.
Vous éditerez ce fichier README.md pour répondre aux questions dans les encarts prévus à cet effet et inserer les figures que vous aurez générées. Ce "repository" vous appartenant, vous pouvez créer tous les repertoires et fichiers necessaires à la conduite du TP.
Seul le langage Python (v3.X) est requis pour ce travail. Il vous est conseillé d'installer un environnement virtuel python pour installer les libraries requises independamment de votre systèmes d'exploitation.
- Le systême de gestion de paquets Conda est très pratique et disponible pour la plupat des systèmes d'exploitation. Une version légère suffisante pour nos besoin est téléchargeable ici
- Si vous disposez d'un interpreteur python 3.X installé sur votre systême virtualenv est désormais "built-in".
Depuis le repertoire de votre repository Git, installez le package scipy et lancez jupyter.
$PATH_TO_CONDA_DIR/bin/conda install -c conda-forge scipy notebook
$PATH_TO_CONDA_DIR/bin/jupyter notebookUne "appliance" IFB a été préparée avec les dépendances Python requises.
Elle est accessible ici.
Jupyter vous permettra d'ouvrir des terminaux SHELL et des notebook Python.
Le repertoire racine de Jupyter est /mnt/mydatalocal/
N/A
Jupyter est une environnement de type notebook permettant l'execution de code python dans des cellules avec une persitance des variables entre chaque evaluation de cellule. Jupyter fournit nativement le support de la librarie graphique matplotlib.
Dans l'inteface de jupyter, créez un nouveau fichier notebook (*.ipynb) localisé dans votre repertoire git. Dans la première cellule copiez le code suivant:
%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import normLa première ligne sert à activer le rendu graphique, pour tout le fichier notebook. Pour dessiner des graphiques, il vous est conseillé de suivre la méthode illustrée par le code suivant que vous pouvez executer dans une deuxième cellule notebook.
fig, ax = plt.subplots()
x = np.linspace(norm.ppf(0.01),
norm.ppf(0.99), 100)
ax.plot(x, norm.pdf(x),
'r-', lw=5, alpha=0.6)
ax.plot(x, np.full(len(x), 0.2),
'b-', lw=1)
fig.show()- Creation des objets
figetax - Ajout successif de graphiques sur la même figure par l'appel à des methodes de l'objet
ax - Affichage de la figure complète via
fig.show() - Evaluation de la cellule pour visualisation dans la cellule de résultats.
L'affichage dans la cellule de rendu du notebook devrait confirmer la bonne installation des dépendances.
La documentation matplotlib est bien faite, mais attention il vous est demandé, pour la construction des graphiques, d'ignorer les méthodes de l'objet plt (frequemment proposées sur le net) et d'utiliser uniquement les méthodes de l'objet Axes. plt est un objet global susceptible d'agir simultanement sur tous les graphiques d'un notebook. A ce stade, son utilisation est donc à éviter.
Un fichier data/TCL_wt1.tsv contient les données d'abondances differentielles mesurées sur une souche sauvage d'Escherichia coli entre deux conditions: avec Tetracycline et en milieu riche. Le contrôle est le milieu riche.
| Accession | Description | Gene Symbol | Corrected Abundance ratio (1.53) | Log2 Corrected Abundance Ratio | Abundance Ratio Adj. P-Value | -LOG10 Adj.P-val |
|---|---|---|---|---|---|---|
| Accesseur Uniprot | Texte libre | Texte libre |
Attention certaines valeurs numériques sont manquantes ou erronées, constatez par vous même en parcourant rapidement le fichier.
Les fiches de toutes les protéines de E.coli sont stockées dans un seul document XML data/uniprot-proteome_UP000000625.xml. Nous allons extraire de ce fichier les informations dont nous aurons besoin, à l'aide du module de la librarie standard XML.etree. Pour vous facilitez la tâche, les exemples suivants vous sont fournis. Prenez le temps de les executer et de les modifier dans un notebook. Vous vous familliariserez ainsi avec la structure du document XML que vous pouvez egalement inspecter dans un navigateur.
from xml.etree.ElementTree import parse, dump
# Parse the E.Coli proteome XML Document
tree = parse('data/uniprot-proteome_UP000000625.xml')
root = tree.getroot()
ns = '{http://uniprot.org/uniprot}' # MANDATORY PREFIX FOR ANY SEARCH within document
# Store all entries aka proteins in a list of xml nodes
proteins = root.findall(ns + 'entry')
# Display the xml subtree of the first protein
dump(proteins[0])# Find the xml subtree of a protein with accession "P31224"
for entry in proteins:
accessions = entry.findall(ns+"accession")
for acc in accessions:
if acc.text == "P31224":
dump(entry)
breakCherchez par exemple le sous arbre de la protéine nommée DACD_ECOLI
Les nombres d'occurences de tous les termes GO trouvés dans le protéome de E.coli sont stockés dans le fichier data/EColiK12_GOcounts.json. Ces statistiques ont été dérivées du fichier data/uniprot-proteome_UP000000625.xml.
- Manipuler les données experimentales tabulées
- Representer graphiquement les données d'abondance
- Construire la pvalue des fonctions biologiques (termes GO) portées par les protéines surabondantes
- Visualiser interactivement les pathways plus enrichis.
La lecture des données au format tabulé est l'occasion de se familliariser avec la librairie pandas.
La fonction read_csv accepte différents arguments de format de données très utiles.
df = pandas.read_csv()Quel est le type de l'objet df?
c'est une dataframe
Que permettent les méthodes suivantes?
permet d'avoir le nb de col lignes => attributs qui donne 1 tuple
les premières lignes
les dernières lignes
liste de headers
donne le type des données par col => object ou float
infos
pareil que dtypes mais donne que pour les col FLOAT
values = df[['Description', 'Gene Symbol']]Quel est le type de values ?
Verifiez si certaines méthodes de DataFrame lui sont applicables.
Ce type supporte l'accès par indice et les slice [a:b]
On peut accéder aux valeurs du DataFrame via des indices ou plages d'indice. La structure se comporte alors comme une matrice. La cellule en haut et à gauche est de coordonnées (0,0).
Il y a différentes manières de le faire, l'utilisation de .iloc[slice_ligne,slice_colonne] constitue une des solutions les plus simples. N'oublions pas que shape permet d'obtenir les dimensions (lignes et colonnes) du DataFrame.
df.iloc[ :5 , : ]df.iloc[ : , -1 ]df.iloc[ :5 , [0,2,3] ]Le type des valeurs d'une colonne peut être spécifiée:
- à la lecture
pandas.read_csv('data/TCL_wt1.tsv', sep="\t", dtype = {'Accession': str, 'Description': str, 'Gene Symbol': str,
'Corrected Abundance ratio (1.53)': np.float, 'Log2 Corrected Abundance Ratio': np.float,
'Abundance Ratio Adj. P-Value: (127. T3 Tc WT) / (126. T0 WT)': np.float, '-LOG10 Adj.P-val': np.float})- modifiée à la volée
df = df.astype({'Log2 Corrected Abundance Ratio': float, '-LOG10 Adj.P-val': float } )La méthode loc permet de selectionner toutes les lignes/colonnes respectant certaines contraintes
- Contraintes de valeurs continues
df.loc[(df['-LOG10 Adj.P-val'] < 0 ) & (df['Log2 Corrected Abundance Ratio'] > 0.0 ) ]- Contraintes de valeurs discrètes
df.loc[ df['Gene Symbol'].isin(['fadR', 'arcA'] ) ]1. Charger le contenu du fichier data/TCL_wt1.tsv dans un notebook en eliminant les lignes porteuses de valeurs numériques aberrantes
3. A partir de cette échantillon de ratio d'abondance, estimez la moyenne 
et l'ecart-type 
d'une loi normale.
4. Superposez la densité de probabilité de cette loi sur l'histogramme. Attention, la densité de probabilité devra être mis à l'echelle de l'histogramme (cf ci-dessous)
hist = ax.hist(_, bins=100) # draw histogram
x = np.linspace(min(_), max(_), 100) # generate PDF domain points
dx = hist[1][1] - hist[1][0] # Get single value bar height
scale = len(_)*dx # scale accordingly
ax.plot(x, norm.pdf(x, mu, sqrt(S_2))*scale) # compute theoritical PDF and draw it
A l'aide de la méthode scatter representer  = f(\text{Log}_2(\text{abundance ratio})))
Sont condidérées comme surabondantes les proteines remplissant ces deux critères:
Nous allons implementer une approche ORA (Over Representation Analysis) naive.
Quelles sont leurs identifiants UNIPROT ?
Les entry du fichier data/uniprot-proteome_UP000000625.xml présentent des balises de ce type:
<dbReference type="GO" id="GO:0005737">
<property type="term" value="C:cytoplasm"/>
<property type="evidence" value="ECO:0000501"/>
<property type="project" value="UniProtKB-SubCell"/>
</dbReference>A ce stade, on se contentera des identifiants GO (eg GO:0005737). Vous pouvez faire un set de cette liste d'identifiants GO pour en éliminer rapidement la redondance.
Nous evaluerons la significativité de la présence de tous les termes GO portés par les protéines surabondantes à l'aide d'un modèle hypergéometrique.
Si k protéines surabondantes porte un terme GO, la pvalue de ce terme sera équivalente à .
Completer le tableau ci-dessous avec les quantités vous semblant adéquates
| Symboles | Paramètres | Quantités Biologiques |
|---|---|---|
| k | nombre de succès observés | |
| K | nombre de succès possibles | |
| n | nombre d'observations | |
| N | nombre d'elements observables |
A l'aide du contenu de data/EColiK12_GOcounts.json parametrez la loi hypergeometrique et calculez la pvalue
de chaque terme GO portés par les protéines surabondantes. Vous reporterez ces données dans le tableau ci-dessous
| identifiant GO | définition | occurence | pvalue |
|---|---|---|---|
Quelle interpretation biologique faites-vous de cet enrichissement en termes GO spécifiques ?
CODE CORRECTION
%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
fig, ax = plt.subplots() # 1 seul dessin ; ax fait des dessins => ici plot choisi
x = np.linspace(norm.ppf(0.01), # norm = loi normale
norm.ppf(0.99), 100) # linspace = comme une array en python
ax.plot(x, norm.pdf(x),
'r-', lw=5, alpha=0.6) #plot fait du discret avec les données de x
ax.plot(x, np.full(len(x), 0.2),
'b-', lw=1) # ici in rajoute 1 dessin sur le premier dessin pour l'enrichir
fig.show()
import pandas
df = pandas.read_csv("./data/TCL_wt1.tsv", sep= "\t")
df
# à gauche c'est le num de ligne qui ne change jamais
#type(df) c'est 1 dataframe = 1 serie de vercteurs
#df.head() => head est une méthode !!!!
print(df.shape)
df= df.dropna() # permet de se débarasser de lignes bizarres
print(df.shape)
df = df.astype({
'Log2 Corrected Abundance Ratio': float,
'-LOG10 Adj.P-val': float } ) #permet de se débarasser des objets et d'avoir ici des floats
df.describe() # maintenant 2 cols en +
#values = df[['Description', 'Gene Symbol']]
#df.iloc[ :5 , : ] # accéder au 5 premières lignes de toutes les col
#df.iloc[ : , -1 ]
#df.iloc[ :5 , [0,2,3] ]
df.loc[(df['-LOG10 Adj.P-val'] > 0 ) & (df['Log2 Corrected Abundance Ratio'] > 0.0 ) ]
['-LOG10 Adj.P-val'].tolist() #toutes les lignes respectent les contraintes
fig, ax = plt.subplots()
field = 'Log2 Corrected Abundance Ratio'
hist= ax.hist(df[field].tolist(), bins=100, rwidth=0.4 )
ax.set_xlabel(field)
ax.set_ylabel('protein count')
fig.show()
#la position de la médiane => plus de wildtype rich que de tetracycline
#on va utiliser cet échantillon pour estimer la moyenne de la loi normale
from math import sqrt
_ = df[field].tolist()
n = len(_)
mu = np.mean(_)
s2 = np.std(_)*np.std(_)
sigma = sqrt((n/(n+1)*s2)) # correction
print(mu, sigma)
x= np.linspace(-4,1.5,100)
dx = hist[1][1] - hist[1][0]
bar_scale= dx * n # generate PDF domain points
ax.plot(x,norm.pdf(x,loc=mu,scale=sigma)* bar_scale)# scale (déviation standart)
# on rajoute ça ( loc=mu,scale=sigma) car mu et sigma ne sont pas les standard 0 et 1
fig.show()