Récemment, j’ai eu à chercher une structure chimique donnée dans une liste de structures. En utilisant les librairies python de chimie informatique pybel et rdkit, je suis facilement arrivée à faire cette recherche, mais celle-ci prenait beaucup trop de temps à mon goût. En me demandant comment l’accélérer, je me suis souvenue de l’article de blog de Jean-Philippe intitulé « Faites travailler vos CPUs ! ». J’ai donc décidé de suivre ses instructions et de faire travailler mes CPUs!

But

Trouver une molécule (une structure chimique donnée) dans une liste de molécules.

Implémentation séquentielle

Notez que pour le présent projet, ma liste de molécules contient des objets pybel alors que ma requête est un objet rdkit. Sans cette contrainte, je n’aurais pas nécessairement utilisé les deux librairies: mon code aurait été écrit en choisissant l’une ou l’autre.


def exactSearch (q, start, end) : 
   ''' Chercher une molécule dans une liste de molécules'''
   ''' Retourne le nombre de fois que la molécule a été trouvée'''


    matches = []
    for mpybel in molecules[start, end]:  ## molécules pybel 
        fp1 = mpybel.calcfp()  ## pybel
        smi2 = Chem.MolToSmiles(q, isomericSmiles=True)  # molécule rdkit
        mpybel2 = pybel.readstring("smi", smi2 ) # convertit la molécule rdkit en molécule pybel
        if mpybel.formula == mpybel2.formula : 
            fp2 = mpybel2.calcfp()
            tanimoto = fp1 | fp2

            if tanimoto == 1:          
               matches.append(mpybel)               
    return len(matches)

n_matches = exactSearch (query, 0, len(molecules))

Ma fonction prend 29.8 secondes pour trouver les trois occurrences de ma requête parmi les 25791 molécules de la liste. Je suis patiente, mais pas assez pour attendre 30 secondes chaque fois que je dois chercher une molécule!

Implémentation avec multiprocessing

La librairie python multiprocessing permet d’utiliser plus d’un CPU pour faire le travail. Pour ce faire, je n’ai que deux modifications à faire à mon code : (1) je dois créer une fonction « wrapper » et (2) définir les tâches à exécuter en parallèle.


import multiprocessing
def exactSearch_wrapper (args) : 
    return exactSearch(*args)

def exactSearch (q, start, end) : 
   ... # same function as above

p = multiprocessing.Pool (processes=4)
tasks = [(query, i, i+1000]) for i in range(0, len(molecules), 1000)]
searches = p.map_async(exactSearch_wrapper, tasks)
p.close()
p.join()

Ce code trouve les trois occurrences de ma requête parmi ma liste de molécules. Toutefois, comme j’ai divisé la recherche en plusieurs recherches dans des listes de 1000 molécules et que j’utilise 4 CPUs, la recherche prend maintenant 7.8 secondes. J’ai pu diviser ma tâche principale en sous-tâches puisque les comparaisons de ma requête aux molécules de la liste sont des opérations indépendantes qui peuvent se faire en parallèle.

Dernières remarques

D’autres options m’auraient permis d’accélérer ma recherche. Cependant, mon but ici était de présenter un exemple simple (espérons!) d’utilisation de la librairie multiprocessing.

Parmi ces autres options, il y a l’utilisation d’un « profileur » de code (« code profiler ») pour identifier quelles opérations sont coûteuses en terme de temps et voir s’il est possible d’apporter des modifications dans le but d’optimiser le code.
Voici le résultat de profilage que j’ai obtenu pour ma fonction exactSearch :

Timer unit: 1e-06 s

Total time: 28.8646 s
File: 
Function: exactSearch at line 1

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def exactSearch0 (q, mols) : 
     2                                              
     3         1            1      1.0      0.0      matches = []
     4                                           
     5     25792        16267      0.6      0.1      for mpybel in mols:  ## pybel molecules
     6     25791     13097301    507.8     45.4          fp1 = mpybel.calcfp()  ## pybel
     7                                                   
     8                                               
     9     25791      1043002     40.4      3.6          smi2 = Chem.MolToSmiles(q, isomericSmiles=True)
    10     25791     14247872    552.4     49.4          mpybel2 = pybel.readstring("smi", smi2 )
    11                                           
    12     25791       459966     17.8      1.6          if mpybel.formula == mpybel2.formula : 
    13                                                       
    14                                           
    15         3          151     50.3      0.0              fp2 = mpybel2.calcfp()
    16         3           13      4.3      0.0              tanimoto = fp1 | fp2
    17                                           
    18         3            3      1.0      0.0              if tanimoto == 1: 
    19         3            2      0.7      0.0                  matches.append(mpybel)
    20                                                           
    21         1           67     67.0      0.0      print len(matches)   
    22         1            1      1.0      0.0      return matches

Je peux tout de suite voir que les opérations coûteuses sont la génération du fingerprint fp1 pour chacune des molécules de la liste ainsi que la conversion de la molécule de la requête. Ces deux opérations sont présentement effectuées 25791 fois et prennent environ 94% des 30 secondes nécessaires pour l’exécution. La seconde opération n’a clairement pas sa place dans la boucle puisque je n’ai besoin de la faire qu’une seule fois. Pour ce qui est de la génération de fp1, comme je fais une recherche exacte, je n’en ai besoin que si les formules chimiques des molécules comparées sont identiques. Si elles ne le sont pas, je sais que ce sont deux molécules différentes. Je n’ai donc pas besoin de calculer fp1 pour les 25791 molécules. L’optimisation du code de ma fonction, suite à cette inspection, réduit son temps d’exécution de 29.8 secondes à 0.137 secondes. Dans ce cas-ci, il n’y a aucun gain à utiliser plus d’un CPU: ça pourrait même prendre plus de temps!

0.137 secondes, c’est beaucoup plus raisonnable!!

Voici le code optimisé de ma fonction de recherche :


def exactSearch (q, start, end) : 
    matches = []
    smi2 = Chem.MolToSmiles(q, isomericSmiles=True)
    mpybel2 = pybel.readstring("smi", smi2 )
    for mpybel in molecules[start:end]: 
        
        if mpybel.formula == mpybel2.formula : 
            fp1 = mpybel.calcfp()  
            fp2 = mpybel2.calcfp()
            tanimoto = fp1 | fp2

            if tanimoto == 1: 
                matches.append(mpybel)
                
    return len(matches)