Category Archives: Quimioinformàtica

Com buscar molècules similars

Moltes vegades tenim un conjunt de molècules i volem seleccionar, en un segon conjunt, aquelles parelles de molècules que més s’assemblen. Per exemple, tenim n molècules que hem vist que són actives per una determinada proteïna i volem buscar fàrmacs que s’hi assemblin (com més semblança, possiblement, millor semblança en activitat).

He fet un petit script que fa això fent servir el coeficient de Tanimoto (a partir dels fingerprints) per mesurar la semblança. El codi agafa les molècules query i les compara amb totes aquelles molècules de la base de dades on volem trobar molècules semblants (en l’exemple, els fàrmacs). Finalment, n’extreu aquelles parelles que estan dintre d’uns límits de semblança definits (éssent el coeficient 1 per dues molècules iguals i 0 per dues molècules completament diferents). Els resultats, les parelles de molècules, es desen en un sd file on s’ha afegint un camp relatiu a la semblança.

El codi està desenvolupat amb python fent ús de les llibreries de pybel. La crida del codi és ben senzilla, simplement cal indicar-li el fitxer amb les molècules query, el fitxer amb les molècules de referència i el límit de semblança (predefinit a 0.9).

1
python SimMolecules.py -i query.sdf -d db.sdf

A continuació teniu el codi:

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
"""
SimMolecules, a script to compare two molecular databases
Alfons Nonell-Canals - July 2010
"""
 
import fileinput
from pybel import *
import optparse
 
p = optparse.OptionParser()
p.add_option('--DBFile', '-d', default='0', help='Database of molecules used as a reference')
p.add_option('--sdf','-i',default='0', help='File with the query molecules')
p.add_option('--CutOff','-c',default='0.9', help = 'Similarity CutOff.')
options, arguments = p.parse_args()
 
#Read commandline options
db = options.DBFile
input= options.sdf
cutOff = float(options.CutOff)
 
RefDbfp = {}
RefDbmols = {}
 
for mol in readfile('sdf',db):
    fp = mol.calcfp()
    RefDbfp[mol.title] = fp
    RefDbmols[mol.title] = mol
 
for mol in readfile('sdf',input):
    fp = mol.calcfp()
    for RefMol in RefDbfp:
        tanimoto = fp|RefDbfp[RefMol]
        if tanimoto >= cutOff:
            print mol.title, RefMol, tanimoto
            mol.data['Similarity'] = tanimoto
            outRef = input.replace('sdf','Refs.sdf')
            out = open(outRef,'a')
            out.write(mol.write('sdf'))
            RefMolMol = RefDbmols[RefMol]
            out.write(RefMolMol.write('sdf'))
            out.close()

Pymol en mans d’Schrödinger

Segons llegeixo a Depth-First, després de la mort del desenvolupador de Pymol, el manteniment del programa passarà a mans de la coneguda empresa de programari quimioinformàtic Schrödinger.

On January 8, 2010 Schrödinger reached an agreement with the estate of the late Dr. Warren L. DeLano to acquire PyMOL. Schrödinger will take over continued development and maintenance, as well as support and sales of PyMOL, including all current subscriptions. Schrödinger will also continue to actively support the open-source community of PyMOL.

Prior to Warren’s tragic and unexpected passing, he had been working closely with Schrödinger to progressively integrate PyMOL with Schrödinger’s graphical interface, Maestro. With a great sense of humility, we will work hard to pick up as best we can where Warren left off and will strive to honor his memory by continuing in the spirit and tradition of PyMOL.

Afegir el títol de la molècula a l’sdf amb pybel

Pybel és un conjunt de classes escrites en Python per treballar amb l’OpenBabel. És una eina molt potent que et permet fer ús de totes les utilitats de l’OpenBabel des de Python. Jo l’empro molt i a diari!

Per, per exemple, llegir un smile i escriure’n un sdf només cal fer el següent:

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
#Importar el pybel
from pybel import *
 
#Ara llegirem l'smile
mol = readstring('smi', 'CCCC(Cl)CCF')
 
#Obrim el fitxer de sortida
out =Outputfile('sdf','NomDelFitxer.sdf')
#Hi escrivim i tanquem
 
out.write(mol)
out.close()

Això va molt bé però m’he trobat que escriu sdfs sense posar cap títol (nom) a les molècules. Pel que ho estic emprant ara però, necessito posar els InchiK com a títol de cada molècula. Hi he estat hores i hores, intentant entendre la classe “Outputfile”, però he desistit, no crec que es pugui fer, així, directament.

Com a alternativa he decidit emprar els mètodes natius de Python per escriure fitxers, barrejat amb el Pybel (l’opció molecula.write). Així doncs:

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#Importar el pybel
from pybel import *
 
#Ara llegirem l'smile
mol = readstring('smi', 'CCCC(Cl)CCF')
 
#Obrim el fitxer de sortida (tal com es fa amb Python). Li diem que és per escriure
out = open('NomDelFitxer.sdf','w')
 
#Hi escrivim el títol
out.write('TítolDeLaMolecula')
 
#I ara, la molècula
out.write(mol.write('sdf'))
 
#Finalment, tanquem
out.close()

De moment, és l’única solució que he trobat i va prou bé.

Es presenta la ChemPedia

chempediaAvui s’ha obert al públic la ChemPedia, un repositori obert, públic i gratuït de molècules. Quan algú hi afegeix una nova molècula, aquesta passa una revisió i se li assigna un codi únic. Posteriorment, els usuaris del portal li poden assignar noms.

L’eina crec que està molt bé i pot ser molt útil per a gestionar els sinònims (el paracetamol, per exemple, es coneix amb els noms de paracetamol i acetaminophen). Els usuaris, a banda d’assignar noms, poden votar pels noms ja assignats. Amb tot, els usuaris van obtenint una puntuació (en funció de si s’assignen noms a molècules seves, si assignen noms a d’altres molècules. Com més puntuació té l’usuari, més confiança se li dóna.

Com a punt fort hi destacaria el tema dels sinònims i la facilitat de ser-ne usuari (funciona amb openID, usuari Google, Yahoo!, Wordpres,…). Com a punts negatius però, el fet de crear una enèssima classificació (que els gestors del portal li assignen un rang d’única) i no fer ús d’eines existents i universals com els Inchis. Per altra banda, per afegir una nova molècula cal dibuixar-la en un “applet”, no permetent afegir ni smiles ni “pujar” sdf (cosa que simplificaria enormement la tasca).

La idea és bona però crec que caldrà veure com evoluciona.

Script per calcular els descriptors de lipinski

Els “descriptors de Lipinski”, de la norma del cinc, inclouen del pes molecular, el coeficient octanol/aigua (cLogP), els punts donadors d’enllaç d’hidrogen (HBD) i els punts acceptors d’enllaç d’hidrogen (HBA).

Aquests dies he hagut de treballar amb diferents molècules d’una base de dades gegant i em calia calcular diferents descriptors. Per als de Lipinski, he fet aquest petit script (en python). Requereix pybel i es basa en un que hi ha publicat a l’article que presentapybel.

?View Code PYTHON
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
"""
Script to compute the Lipibski descriptors
Alfons Nonell-Canals - October 2009
 
USAGE:
python Lipinski.py -i sdfFile
 
"""
 
import optparse
from pybel import *
from openbabel import * 
 
p = optparse.OptionParser()
p.add_option('--infile','-i',default='default')
options, arguments = p.parse_args()
#files
infile = options.infile
 
#function to compute the Lipinski descriptors
HBD = Smarts("[#7,#8;!H0]")
HBA = Smarts("[#7,#8]")
 
def lipinski(mol):
   desc = []
   desc.append(mol.molwt)
   desc.append(mol.calcdesc(['LogP']) ['LogP'])
   desc.append(len(HBA.findall(mol)))
   desc.append(len(HBD.findall(mol)))
   return desc
 
for mol in readfile("sdf",infile):
    molName = mol.title
    desc = lipinski(mol)
    print molName,desc

Congressos per compartir coneixement

La setmana passada vaig estar de congrés. Vaig assistir al “27th Noordwijkerhout-Camerino-Cyprus Symposium, Trends in drug research”. Era el meu primer congrés com a post-doc i, per tant, la primera vegada que assistia a un congrés de química mèdica. El congrés en general em va agradar, vaig aprendre (com dic, era el primer congrés del tema al que anava) però em va decebre el poc intercanvi de coneixement que hi va haver.

M’explico. Fins ara, els congressos als que he anat, a gent, en fer una xerrada exposava un problema, explicava amb força detall la metodologia emprada per abordar-lo i, finalment, mosatrava uns resultats i conclusions. En el congrés de la setmana passada, en general, no va ser així. En la majoria de xerrades se’ns exposava un problema, una caixa negra el resolia i al final, se’ns mostraven uns resultats i unes conclusions. I, a més a més, sovint, sense massa detalls per qüestions de “propietat intel·lectual”.

El per què? doncs, perquè la majoria de xerrades la feien laboratoris farmacèutics. I és aquí on rau el problema. Entenc que en aquests tiups de congressos hi ha d’haver presència de la indústria, em sembla bé i tot! també puc arribar a entendre que, quan un laboratori fa una xerrada, no expliqui tots els detalls. Però, el que ja no entenc és com pot ser que en un congrés, la gran majoria de xerrades les facin laboratoris farmacèutics i que, per tant, el nivell d’intercanvi de coneixement sigui tant baix.

Això és el que m’ha falalt d’aquest congrés. Com en molts congressos, hi havia expositors de diferents laboratoris i empreses, cosa que em sembla bé ja que són patrocinadors. Però, d’aquests expositors, tots, van xerrar (això ja em sona a “si pagues per un expositor, podràs xerrar”), però, a més a més, hi havia patrocinadors dels que… tots, també van xerrar.

Tinc la sensació que, qui ha xerrat ha estat perquè ha pagat. La gran majoria de xerrades eren d’empreses i, per tant, explicàven ben poc. I, de les poques universitats que han xerrat, quasi totes han mostrat resultats de col·laboracions fetes amb els laboratoris i on molts detalls no eren explicats per “no violar la propietat intel·lectual”.

Quan pagues per anar a un congrés pots entendre que hi haurà uns patrocinadors, que algunes xerrades les faran ells,… però, a canvi, esperes aprendre, veure coneixement. Aquests dies, he après, no ho nego (com he dit, era el meu primer congrés del tema), però m’he quedat una mica amb la sensació que hem pagat un congrés per anar a una fira de mostres on, tot de comercials encorbatats (sí, la majoria homes i ben “entrajats”), ens expliquin les meravelles que fan les seves empreses, però que ens les expliquin agafant com a exemples antics projectes solucionats per caixes negres.

Calcular el coeficient de Tanimoto amb Pybel

Necessito comparar diverses molècules. Encara no hi entenc massa, però, crec que una de les millors maneres és fer servir el coeficient de Tanimoto [en]. Això, el càlcul d’aquest coeficient l’he de fer de forma reiterada, necessito anar calculant el coeficient per cada molècula que generi i comparar-la així amb una molècula final. A més a més, vull integrar-ho dins un programa que tinc fet amb Python.

La millor manera, doncs, és fent ús de Pybel, un mòdul Python per fer servir eines de l’Openbabel des de codi Python. La manera de fer-ho és ben senzilla. Des de Python:

>>> import pybel #importarem el mòdul
>>> smiles = [‘CCCC’, ‘CCCN’] #Array de molècules
>>> mols = [pybel.readstring(“smi”, x) for x in smiles] #Llegim les molècules
>>> fps = [x.calcfp() for x in mols] #calculem els fingerprints (descriptors)
>>> print fps[0] | fps[1] #finalment, calculem el coeficient
0.333333333333

I… llestos! suposo que hi deuen haver maneres més acurades i que, com tota cosa estadística, com més descriptors es tinguin per cada molècula, millor, però, per a començar, no va gens malament i el mètode és molt fàcil d’implementar.