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.

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

"""
Script to compute the Lipibski descriptors
Alfons Nonell-Canals - October 2009
Chemogenomics Lab
 
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

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.

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.

Interessant utitlitat que fins ahir desconeixia totalment: els “sets”. Bàsicament, és un conjunt d’eines per a comparar, unir,… conjunts de dades. Amb les funcions de la classe podem, per exemple, obtenir de maneta senzilla la diferència entre dos conjunts de dades, obtenir les diferències, actualitzar-ne un a partir d’un segon.

Un exemple ben senzill és quan es té una llista A:

a = ['carboni', 'nitrogen', 'hidrogen']

I ho volem comparar amb una llista B:

b = ['oxigen', 'nitrogen', 'molibdè', 'calci', 'sofre']

Ara volem comparar les dues llistes, primer generem els dos sets

as = set(a)
bs = set(b)

Ara, per exemple, obtindrem una llista consens, és a dir, amb els elements comuns d’ambdues llistes:

as.intersection(b)

Això és un petit exemple, però, a partir d’aquí, es poden fer meravelles. Aquí hi ha un petit tutorial, molt didàctic.

Dies enrere ho vaig llegir al Nature News però fins avui no m’ho he mirat: Google flu trends. Ara Google també pot “predir” els brots de grip. Més que predir, l’eina permet avisar d’on i quan hi ha brots de grip i ho fa mitjançant les cerques de la gent.

Quan hi ha una “passa” de la grip, la gent té més tendència a buscar certes paraules, com ara, grip, antigripal, febre,… la sobreposició d’aquestes cerques amb la informació geogràfica del lloc des d’on es fan (quan es cerca alguna cosa a Google, aquest emmagatzema les paraules cercades i la IP de l’ordinador des d’on es fa la consulta) pot donar idea dels llocs on hi ha grip. Amb això, la gent de Google i Yahoo! han estudiat els brots de grip dels darrers anys (als Estats Units) i han comparat aquesta informació amb les dades que té el CDC d’Atlanta (organisme que controla l’evolució de les malalties). Els resultats han estat sorprenents ja que pràcticament coincideixen.

Els mètodes clàssics de control de brots de malalties es basen en la informació aportada pels diferents centres sanitaris als organismes de control. Això, requereix que passin uns dies abans no es pot declarar un brot en una zona concreta. L’eina de Google és molt més ràpida ja que, pràcticament, pot ser en temps real (la informació que donen avui és de fa tansols dos dies). Seguint amb les comparatives entre Google i CDC, el primer ha anunciat amb deu dies d’antelació, els casos de grip d’enguany.

Aquesta eina no pot substituir, en cap cas, els mètodes clàssics de control de malalties ja que són molt més fiables (recordo que Google flu trends funciona a partir de les cerques que fan els usuaris i usuàries del cercador i són estimacions) però sí que poden ajudar i anticipar la resposta en casos de pandèmies greus.

De moment només està disponible pels Estats Units.

Fins ahir tenia un problema prou complicat: estic treballant en un programa en Python que requereix múltiples i constants connexions a una base de dades mySQL. Quan treballava en local, no tenia cap problema però, quan executava el programa en màquines diferents apuntant cap a un servidor de bases de dades (com passaria, per exemple, en un cluster), de tant en tant i de manera aleatòria, el servidor mySQL es quedava “tonto” i no responia. Aquesta pèrdua de connexió feia que el programa “petés”.

Amb l’administrador de sistemes del grup hem treballat per mirar què passa amb el mySQL ja que sembla que el problema és allà, però no hem trobat solució. No hi ha ni problemes ni amb el límit màxim de connexions, ni amb la memòria, cpu,… en fi, per algun motiu que desconeixem, el servidor es satura i té petites penjades. Així dons, la solució passa per fer alguna cosa al programa perquè no peti quan la base de dades no respongui.

La solució està en les “excepcions” (de l’anglès exception, desconec si aquesta traducció literal és vàlida). És a dir, per defecte, quan python intenta accedir a la base de dades, si aquesta no respon, es crea una excepció que, si no se li diu res, provoca la sortida del programa. La idea està en modificar aquesta excepció perquè quan python vegi que la base de dades no respon, la solució no sigui sortir del programa sinó reintentar-ho. He jugat amb el try i except.

db = DB()
cursor = db.cursor()
sql = ’select count(*) from molecules’
number = cursor.execute(sql)

class DB:
def cursor(self):
try:
bd = MySQLdb.connect(host=HOSTNAME, user=”USERNAMEl”, passwd=”PASSWD”,db=”DB_NAME)
cursor = bd.cursor()
return cursor
except MySQLdb.OperationalError, message:
return self.cursor()

La primera part del codi inicia la classe DB i n’obté el “cursor”, operador amb el que realitzarem les crides sql. La part més interessant però és la classe DB, que conté la funció “cursor”. Aquesta funció connecta amb la base de dades i retorna el “cursor”, comanda dins l’ordre try. Però, el fet de jugar amb el try i except fa que el programa primer provi el que hi ha a try però, si falla i troba una excepció que coincideix amb la indicada a except, en l’exemple MySQLdb.OperationalError, en comptes d’aturar l’execució del programa (el que faria per defecte) fa el que se li diu, en aquest cas que torni a executar la funció.

Obprop és una aplicació inclosa al paquet Openbabel i que calcula alguns descriptors per una molècula donada. Cal passar-li la molècula en un fitxer sdf (hi pot haver diverses molècules dins el fitxer) i et retorna la seva fórmula, l’SMILE canònic, el pes molecular, el logP,…

Per fer-la servir simplement cal executar l’ordre:

obprop fitxer_amb_molècules.sdf

I, com a resultat, s’obté quelcom com ara:

name [Name]
formula [Formula]
mol_weight [Molecular Weight]
exact_mass [Isotopic Mass]
canonical_SMILES [String]
num_atoms [Number]
num_bonds [Number]
num_residues [Number]
sequence [Residue Sequence]
num_rings [Number of Rings (by SSSR)]
logP [Number (octanol-water partition)]
PSA [Number (topological polar surface area)]
MR [Number (molar refractivity)]
$$$

Eina útil, senzilla i ràpida per obtenir les quatre dades bàsiques d’una molècula.

Avui he descobert un mòdul de Python amb una funció molt i molt interessant re.split. Segurament és una xorrada, però és ben útil. Per què? per que separa elements d’una cadena fent ús de diferents patrons. La funció que jo emprava fins al moment, split, és útil quan per separar una cadena hi ha només un element comú.

Un exemple ben clar: una molècula. Tinc la cadena següent: C-C-N-C-O. Per separar els diferents àtoms és ben senzill, simplement cal fer un split indicant que el separador és el guionet -. La funció retorna un array.

molecule = ‘C-C-N-C-O’
atoms_array = molecule.split(’-',”)

Ara, imaginem una molècula on, a banda dels enllaços senzills (-) en tenim de dobles i triples (= i #, respectivament i en notació SMILE). La cadena exemple seria: C#C-N-C=O. Aquí, els marcadors que separen els diferents àtoms no segueixen un únic patró i poden ser -, = i #. La funció anterior no ens serveix. Aquí és on entra en joc la funció que he comentat:

import re
molecule = ‘C#C-N-C=O’
atoms_array = re.split(’[\-\=\#]‘,molecule)

El patró, aquí \-\=\#, pot ser divers: símbols, lletres, números, qualsevol dígit número enter,… En fi, per mi m’és molt còmode

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, 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.