Sans avoir trouvé comment exporter les données générée par MOE sous forme exploitable, ce script permet l'analyse complète des données obtenues après un copié - collé depuis le pdf des résultats vers une feuille excel enregistrée en csv. (Explications d'utilisation en bas de page).

CSV_full_analysis.py

#!/bin/python3
##Florent DUFOUR
###Novembre 2017


## -- IMPORTATIONS -- ##

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

input_raw = pd.read_csv('input.csv', header=2, sep=(";")) #Le dossier courant doit contenir 'input.csv'

## -- DEFINITION DES FONCTIONS POUR PLUS TARD -- ##

def cleaner (df):
    for k in df:
        if ((k != 'mol ') and (k!= 'S ')):
            del df[k]
    df = df.dropna(axis=0, how='any')
    df = df.drop(df[df['S '].map(lambda x: x.startswith('S'))].index)
    return df

def sorter (df):
    df['S '] = df['S '].astype('float')  #Convesion des str en float
    df = df.sort_values(by='S ')         #On peut maintenant trier les énergies
    return df

def enrichment(percent,name_list):       # D'après Lucas S.
    Size=(len(name_list)/100)*percent
    count, target_number, total_number = 0,0,0
    
    for j in name_list:
        if 'ZINC' not in j:
            target_number += 1
            total_number += 1
            count += 1
        else:
            total_number += 1
            count += 1
        if count>Size:
            break
    E=target_number/total_number
    print ('Enrichment for', percent,'% corresponds to',target_number,'target in',total_number,'total molecules')
    return E

## -- COEUR DU PROBLEME, calculer les valeurs pour tracer la courbe de ROC -- ##

output_cleaned = cleaner(input_raw)
output_sorted  = sorter(output_cleaned)

list_Mol = output_sorted['mol '].tolist() #Les molécules à tester sont dans une liste propre

name_List = [] #Initialisation axe des X
hauteur = [0]  #Initialisation axe des Y

for k in (list_Mol):         #Les poses à évaluer
    if (k not in name_List): #C'est une nouvelle molécule
        name_List += [k]     #On enregistre qu'on l'a déjà rencontré
        if ("CHEM" in k):    #C'est une nouvelle molécule qui est un positif
            hauteur += [hauteur[-1]+1] #On augmente d'un cran
        else:
            hauteur += [hauteur[-1]]   #On stagne

## -- COURBE DE ROC -- ##

X = np.linspace(0, 1, num=len(hauteur)) #Axe des X
Y = [x / 7 for x in hauteur]                       #Axe des Y pour les valeurs de ROC
Z = np.linspace(0, 1, num=len(hauteur))            #Axe des Y pour droite aléatoire

plt.plot(X,Y)       #Courbe de ROC
plt.plot(X,Z, '--') #Courbe de hasard

## -- STATISTIQUES -- ##

area_exp = np.trapz(Y,X,dx=1) #Aire sous la courbe expérimentale
area_ale = np.trapz(Z,X,dx=1) #Aire sous la courbe aléatoire
diff = (area_exp - area_ale)
Etot=7/107

## -- TOUS LES RESULTATS -- ##

output_sorted.to_csv("output.csv") #Le fichier CSV propre est ajouté dans le dossier de travail
print("\n ROC AUC=", area_exp,"\n","ROC=", diff,"\n", "Good protocole ? = ", diff>0, "\n") # ROC_AUC, ROC, Validation du protocole
print('EF5% =',enrichment(5,name_List)/Etot)  #EF5%
print('EF10%=',enrichment(10,name_List)/Etot) #EF10%
print('Nombre de molécules identifiées sur les 107:',len(Y)-1)

plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve obtained with our best protocole \n against the VS1 database')

plt.show() #Courbe de ROC

Comment utiliser ce script

• Ouvrir le fichier .mdbdans MOE
• L'imprimer en pdf avec pdf forge
  • Il est important de presser 'Print' et d'utiliser pdf forge sans quoi le copier-coller dans excel ne fonctionnera pas correctement
  • Le tableau de résultats dans le pdf ne doit pas déborder sur plusieurs pages en largeur
• Ouvrir le pdf depuis un mac et faire:
  • cmd-A pour sélectionner toutes les données
  • cmd-C pour copier toutes les données dans le presse papier
  • Dans une feuille excel faire cmd-V pour coller les données
  • Enregistrer la feuille en .csv sous input.csv avec ";" comme séparateur de valeur.
• Mettre ce script dans le même dossier
• Executer le fichier python (voir ci dessous depuis un mac avec un exemple de résultats)

Remerciements: Lucas S. pour l'enrichissement

Testé avec python 3.6.3 sous macOS et Ubuntu avec Jupyter Notebook et python3 en ligne de commande

Note de compatibilité

• Les dépendances suivantes doivent être installées sur la machine: pandas, numpy et matplotlib. Voir Ici pour les installer avec pip3
• Certaines versions de python ne supportent pas les affectations multiples. On a:

#La ligne 31:

count, target_number, total_number = 0,0,0

#Qui devient:

count = 0
target_number = 0 
total_number = 0

Si le CSV est généré autrement qu'avec Excel

Zamzar.com permet par exemple de traduire le pdf en .csv. Le script doit alors être ajusté:

• Notez que les entêtes des colonnes se terminent par un espace avec excel. On a en fait "mol " et "S " pour les 2 colonnes d'intérêt. il faut les enlever du script si on utilise Zamzar.
  • "mol " devient "mol" partout dans le script
  • "S " devient "S" partout dans le script
• La taille du header devra être ajusté (typiquement = 1)
• Le séparateur ne sera pas forcément ";" (typiquement = ",")

On obtient alors:

input_raw = pd.read_csv('input.csv', header=1, sep=(","))

Si vous n'arrivez pas à ouvrir le pdf généré par MOE

C'est que ce n'est en fait pas un .pdf mais un .ps. Il suffit d'ajouter l'extension qui convient et ça devrait rouler (mac & linux en tout cas)

Un petit supplément pour éviter de passer par un pdf (utilise les SMILES pour assigner les débuts de ID 'CHEMBL' ou 'ZINC' dans une colonne supplémentaire).

import pandas as pd

d = pd.read_csv("input_screening_results.csv", sep = ";") # Depends on the separator in the file
d1 = pd.read_csv("input_mol_screened.csv", sep = ";")

dico_mol = {d1['mol'][n]:d1['Activity'][n] for n in range(len(d1))}
print len(dico_mol.values())
dico_ID = {1:'CHEMBL', 0:'ZINC'}
SMILES_ID = [[d["mol"][n], dico_ID[dico_mol[d["mol"][n]]]] for n in range(len(d))]
    
file_ID = pd.DataFrame(SMILES_ID, columns = ['SMILES_2', 'mol'])
file_ID
#print d.columns
d = d.rename(columns ={'mol':'SMILES'})
#print d.columns
list(file_ID['mol'][:]).count('CHEMBL')
d2 = d.join(file_ID)
#d2
d2.to_csv('./filename.csv', sep = ",")

Avec du coup quelques modifications de 3 lignes du script de Florent:

input_raw = pd.read_csv('filename.csv', sep=(","))
list_Mol = input_raw['mol']
#output_sorted.to_csv("output_f.csv") #Le fichier CSV propre est ajouté dans le dossier de travail #****on inactive cette ligne