Category Archives: Química

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()

Codi per posar un camp d’un sdf al títol de la molècula

Sovint, mentre manipulem fitxers sdf ens trobem que volem canviar el títol de la molècula. Moltes vegades, com a títol hi volem posar algun camp de les dades addicionals de la molècula en qüestió. M’he fet un petit script per fer aquesta tasca de manera automàtica i senzilla.

El programet llegeix l’sdf i, molècula a molècula, en llegeix el camp que volguem, l’assigna com a títol i desa la molècula en un nou fitxer. A més a més, agafa el títol original i el desa en un camp “old_title”.

?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
"""
Script to move a field to an id in a sdf
Alfons Nonell-Canals - April 2010
"""
 
import optparse
import sys
from pybel import *
 
p = optparse.OptionParser()
p.add_option('--sdf', '-i', help='SDF file to modify')
p.add_option('--field', '-f', help='Field to you want to set as molecule title')
options, arguments = p.parse_args()
 
if not options.sdf  or not options.field:
        p.error('An input sdf file and a field are required')
 
molfile = options.sdf
field = options.field
 
outFile = molfile.replace('.sdf','')+'.Mod.sdf'
out = Outputfile('sdf',outFile)
 
for mol in readfile("sdf", molfile):
        oldTitle = mol.title
        mol.data['old_title'] = oldTitle
        newtitle = mol.data[field]
        mol.OBMol.SetTitle(newtitle)
        out.write(mol)

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

A què m’he dedicat els darrers cinc anys?

Crec que és hora d’explicar amb calma i una mica detalladament, amb llenguatge planer, a què em dedico, què faig al laboratori/despatx. Primer però, per posar-nos en situació, començaré parlant de la tesi i miraré d’explicar a què m’he dedicat durant els darrers cinc anys.

El mes d’octubre passat vaig acabar una tesi en química computacional. Què vol dir això? doncs que vaig fer una tesi on aplicava mètodes de la química computacional a resoldre problemes químics o bioquímics reals. Com ho feia? introduint una sèrie de molècules (un fitxer amb coordenades espaials “xyz” de la molècula a estudiar) a l’ordinador que, mitjançant uns determinats programes, obtenia l’estructura molecular de menys energia.

Bé, i què és això de l’estructura molecular de menys energia? una molècula, com tot, tendeix a estar en una forma que impliqui la menor energia possible. Com a exemple, penseu en la cera d’una espelma. Quan tenim l’espelma en forma d’espelma, la cera és sòlida però, si li administrem energia, l’excitem, passarà a una forma de més energia, líquida que, quan es refredi, tornarà a l’estat de menys energia, sòlid (tot i que haurà patit una transformació i tindrà una forma diferent a la inicial). Els programes que he emprat per fer els meus “càlculs” busquen quina estructura molecular té la mínima energia.

I què n’he fet d’aquestes estructures de mínima energia? Estudiar reaccions. Una reacció química és un procés de tres etapes. Es comença amb uns reactius per separat (l’espelma i el foc -sense que hagin entrat en contacte-) i s’acaba amb uns productes (l’espelma ja apagada amb una forma diferent a la inicial). Pel mig hi ha hagut un estat excitat anomenat estat de trancisió (l’espelma cremant). Amb l’ordinador es poden estudiar aquests tres estats d’una reacció i analitzar-la així, amb molt detall.

Els mètodes computacionals, sovint, ajuden a entendre el perquè una reacció té lloc en unes circumstàncies i no en unes altres, a analitzar les diferents vies per les que pot anar una reacció, perquè el producte d’una reacció és diferent a l’esperat,…

Durant la tesi, he aplicat aquests mètodes a entendre diferents sistemes de química orgànica i bioquímica. El protagonista ha estat un aminoàcid anomenat selenocistgeïna, aquest és un aminoàcid peculiar, poc freqüent i present en proteïnes que catal·litzen (controlen) reaccions on hi ha intercanvi d’electrons. Doncs bé, la selenocisteïna és molt similar a un altre aminoàcid anomenat cisteïna, però el primer té seleni on el segon hi té sofre. Aquest apartat de la tesi compara (analitzant totes les etapes d’una reacció catal·litzada per un enzim que conté selenocisteïna) quins avantatges té el seleni en front del sofre. Els resultats que tinc de moment mostren que el sistema amb seleni perd els electrons més fàcilment que no pas el sistema amb sofre.

A banda del tema anterior, col·laborant amb dos grups experimentals (dels que treballen amb tubs d’assaig i provetes), la tesi presenta l’estudi de tres reaccions orgàniques: la síntesi de 4-helicenoquinones, la síntesi de 5-helicenoquinones i l’ús de sals de guanidini com a catal·litzadors. En aquests tres projectes s’ha estudiat el perquè de determinats resultats experimentals, clarificant els processos intermitjos de la reacció i justificant el perquè d’algunes de les característiques dels productes.

En la propera entrega parlaré del que faig actualment a la feina.

SMILES (II). Uns quants exemples

L’SMILE més senzill és el metà, representat per una simple C. El segueix l’età (CC) el propà (CCC) i, així, successivament. Molècules més complexes, com ara el l’acetat de metil (CC(=O)OC) ja introdueixen dobles enllaços (=) i cadenes laterals (es posen entre parèntesis.

Cada molècula té el seu SMILE únic i pot ser més o menys senzill representar-lo. Buscant alguns exemples per fer proves (i no trencar-me les banyes fent els SMILES) he trobat una web prou interessant ja que recull un fotimer d’exemples per diferents tipus de molècules. Ho podeu trobar aquí.

SMILES

SMILES, de l’anglès Simplified Molecular Input Line Entry Specification, és una notació per descriure, de manera única, l’estructura d’una molècula fent servir una simple cadena de caràcters alfanumèrics en format ASCII. Aquesta manera d’anomenar les molècules, va néixer als anys vuitanta de la mà de n’Arthur Weininger (qui va desenvolupar diferents programes i algoritmes per manejar-les).

Aquesta notació té un seguit de coses positives que la fan interessant. És un “llenguatge” fàcil d’entendre tant per ordinadors com per persones. És un nom “únic”, és a dir, és una manera universal i única d’anomenar una molècula. Requereix poc espai d’emmagatzematge en sistemes informàtics. Com a principal pega és que si bé és relativament senzill transformar un SMILE a coordenades en dues dimensions (2D), obtenir quelcom en tres dimensions ja és molt més complicat.

Fa poc he descobert aquesta nomenclatura i m’hi estic mirant d’acostumar i entendre ja que és la manera amb que remenaré les molècules durant el post-doc.

Aproximació de la càrrega al llarg d’una reacció

Aquests dies estic jugant a fer dinàmiques moleculars, té el seu què, a més a més, em permet seguir amb més detall com evoluciona una reacció tenint en compte el seu entorn. Una cosa que m’interessa saber, entre d’altres coses, és l’evolució de la càrrega de diferents àtoms al llarg de la reacció, partint de la base que només conec la càrrega inicial i la final.

Per fer aquesta aproximació (evolució de la càrrega) així com d’altres aproximacions de paràmetres que varien al llarg de la reacció, introdueixo el factor “lambda” (λ). Lambda va variant al llarg de la reacció, comença per un valor de zero (reactius) i va fins a 1 (productes). Així doncs, quan tenim un λ=0 estem a reactius, λ0.5, podem dir que és un punt intermedi de la reacció (no té per que ser l’estat de trancisió) i, amb λ=1 som a productes. Així doncs, la variació d’una propietat, per exemple, la càrrega (q) al llarg d¡una reacció es pot explicar com a:


q=qi(1-λ)+qiiλ

On qi és la càrrega inicial, qii la càrrega final i q la càrrega del moment. Si som a reactius, tenim que λ=0, per tant, qiiλ=0 i, per tant, q=qi(1-λ), que, com que λ és zero, queda com a q=qi(1-0) i, finalment, q=qi. Podeu fer la mateixa extrapolació per a λ=1. Així doncs, l’equació ens permet saber l’estat d’una propietat determinada en un moment concret de la reacció. En el meu cas, m’interessa saber la càrrega.

Per a fer un càlcul ràpid, es pot programar un script que no té massa complicació (aquí en teniu un exemple en php). Simplement cal modificar el codi per dir-li quina és la càrrega inicial, la final i a córrer!

Article interessant: Computer simulations on enzyme catalysis: methods, progress, and insights

Article, review1 molt interessant que fa una posada al dia dels mètodes teòrics existents per a estudiar una catàlisi enzimàtica. Està escrit pel Prof. Arieh Warshel, de la Universitat del Sud de Califòrnia, que, tot i centrar-se molt en els seus mètodes (ell s’ho pot permetre), fa un repàs de tot el que podem fer a l’hora d’estudiar com funciona un enzim.

Primer ens planteja el problema: què és una reacció enzimàtica? com la podem estudiar?. Seguidament, ens fa una introducció a com el podem abordar: els estudis amb models de centres actius són vàlids? cal treballar sempre amb tota la proteïna? com podem validar el nostre mètode?. I, finalment, la solució. Ens parla dels mètodes QM amb models de centres actius, després planteja treballar a nivell quàntic amb sistemes més grans, fa un repàs dels diferents mètodes QM/MM, una mica de base teòrica, com funcionen. Després se centra molt en els estudis EVB (Empirical Valence Bond), un dels seus mètodes estrella.

És un review amb les fórmules justes, ni moltes ni poques, no entra en profunditat en cap mètode però la informació que hi dóna és suficient. A mi m’ha agradat per que ofereix una visió de conjunt, no m’ha aborrit la teoria però tampoc l’he trobat fluix. En fi, crec que és una bona introducció (com a mínim pels que no som químics) a com podem estudiar un enzim. Cal dir que estic content ja que considera que els estudis quàntics amb petits models (el que faig jo) són vàlids només en el cas de les metal·lopreteïnes (amb les que treballo jo) i ho justifica dient que davant la presència d’un metall, l’efecte de la resta de la proteïna és molt més limitat.

1- Warshel A., Annu. Rev. Biophys. Biomol. Struct. 2003, 32:425-443

Relaxació en una dinàmica

Darrerament, a la recta final de la tesi, m’he dedicat a aprendre a fer dinàmiques moleculars[en]. És quelcom entretingut que té la seva gràcia. No és ni més difícil ni més senzill que la química teòrica de la quàntica pura i dura (que és al que m’he dedicat pràcticament durant aquests darrers quatre anys). Però potser sí que puc dir que és un xic més divertit. Et permet jugar molt més amb el sistema, treballar amb proteïnes senceres i, sense el detall que et proporciona la quàntica, jugar a fer reaccions, a escalfar la proteïna, a fer mutacions,… Concretament el que estic intentant fer és estudiar la reacció de la Format Deshidrogenasa H fent servir l’EVB (Emprirical Valence Bond), un mètode desenvolupat per en Warshel.

El pas previ a la dinàmica en sí, a la reacció, és la relaxació del sistema. Això consisteix, bàsicament i en grans trets, a agafar la proteïna, posar-la en aigua i anar augmentant gradualment la temperatura en periodes de temps determinats. Això em prepara la proteïna per a poder fer la “reacció” en unes condicions determinades, és a dir, si vull fer la “reacció” a 300 K (27ªC), primer he de preparar la proteïna, posar-la amb aigua i, començant des de 0, anar-li augmentant la temperatura fins als 300. Si comencés directament a 300 K,  l’estructura no s’aguantaria.

Quelcom “divertit” de les dinàmiques, són els vídeos que es poden fer. També es poden fer en quàntica, però tenen “menys coses”. Aquí hi ha una relaxació mooolt curteta, 0.5 fempto-segons. El víedo també és curtet, dotze segons. Més endavant, quan tingui ja una bona dinàmica, ja en faré un una mica més elaborat. Gaudiu-lo ràpid :-)

Enllaç al vídeo

Una explosió al despatx…

Aquesta nit hi ha hagut una explosió a la feina. Quan he arribat he mirat com anaven els càlculs que vaig enviar (posar a la cua de càlcul) ahir la nit. N’hi havia un, una optimització senzilla, que hauria d’haver acabat. En mirar com estava he vist que havia acabat malament. M’he extranyat, com he dit, no havia de tneir cap problema. He obert un programa gràfic per veure l’evolució de l’optimització de la geometria i… fixeu-vos en l’anell de la part superior dreta (una imidazola):

Per algun motiu, part de la gometria se n’ha anat a pendre pel sac. D’aquesta manera no m’ahvia passat mai! Però bé, per sort, això ha passat dins la cpu d’un ordinador i no ha tingut més conseqüències que la pèrdua d’algunes hores de càlcul.

Pel que fa al video, ho sé, pèssima qualitat. He fet un povray de cada cicle d’optimització i els he juntat en un video. Tot plegat amb una resolució que deixa molt que desitjar… és el primer que faig, a veure si més endavant milloro.