Construcción de una librería de compuestos químicos comercialmente disponibles en MolPort en formato SDF y 3D destinada a simulaciones de molecular docking. A continuación se muestra scripts de Python para ejecutar en Windows, en mi caso por comodidad, desde Visial Studio Code y usando como kernel un entorno virtual creado con Miniconda3 que usa pandas, rdkit y otros módulos de Python.

Puedes empezar por solicitar un usuario a MolPort para que te deje acceso a su servidor FTP: altftp.molport.com de forma que puedas descargar los ficheros de texto con la información de los SMILES y Molport ID de las moleculas disponibles en cada momento a la venta, que lógicamente cambian, unas aparecen nuevas y otras dejan de estar disponibles. Para la construcción de esta versión de la base de datos que describo en este post, he descargado su base de datos completa, en forma de 102 ficheros comprimidos con tar.gz, cuya descompresión libera ficheros de texto con la información que ves en la figura inferior derecha.

Cada uno de los ficheros de texto descargado de MolPort tiene información en columnas seperadas por tabuladores. La primera tarea consiste en rescatar la información de las columnas: SMILES y MOLPORTID. Será muy util en este proceso largo y complejo guardar copias de los ficheros que generamos en cada paso. Algunos procesos tardan varias horas en ejecutarse completamente en una máquina con SO Windows 11 y procesador I7. Es probable que ejecutando en Linux se agilicen los tiempos de ejecuación pero será necesario adaptar los scripts de Python que están escritos para ejecutarse en Windows.

El siguiente script de Python, para ejecutar en Windows, busca en la carpeta D:\MolPort-database_v7\ todos los ficheros de texto con extensión *.txt y rescata de cada uno de ellos los valores de las columnas "SMILES" y "MOLPORTID". Genera de salida otro fichero de texto con el mismo nombre que el fichero de entrada más la cadena "-smiles". Ejemplo de nombre de fichero de entrada: fulldb_smiles-000-000-000--000-499-999.txt y de fichero de salida: fulldb_smiles-000-000-000--000-499-999-smiles.txt

Debes instalar pandas, si aún no lo tienes, usando desde la consola CMD (Administrador) la siguiente sintaxis:

pip install pandas

import os
import pandas as pd

# Ruta de la carpeta donde se encuentran los archivos
input_folder = r"D:\MolPort-database_v7"

# Iterar sobre todos los archivos en la carpeta
for file_name in os.listdir(input_folder):
    if file_name.endswith(".txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        # Leer el archivo de entrada
        try:
            data = pd.read_csv(input_file_path, sep="\t")  # Ajustar el separador si no es tabulación
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue
        
        # Verificar si las columnas "SMILES" y "MOLPORTID" existen
        if "SMILES" in data.columns and "MOLPORTID" in data.columns:
            output_data = data[["SMILES", "MOLPORTID"]]
            output_file_name = f"{os.path.splitext(file_name)[0]}-smiles.txt"
            output_file_path = os.path.join(input_folder, output_file_name)
            
            # Guardar el archivo de salida
            try:
                output_data.to_csv(output_file_path, sep="\t", index=False)
                print(f"Archivo generado: {output_file_path}")
            except Exception as e:
                print(f"Error al escribir el archivo {output_file_name}: {e}")
        else:
            print(f"El archivo {file_name} no contiene las columnas necesarias.")

Este ejemplo te muestra la información que has rescatado y te permite adevertir que MolPort te deja información que será necesario eliminar pensado en que vas a usar la base de datos para simulaciones de docking y necesitas ficheros SDF en los que cada molécula sea única, luego tendremos que eliminar iones, sales, etc. Los SMILES separados por "." (punto) indican que hay más de una molécula. Mi criterio es quedarme con la cadena de texto más grande de cada línea que tenga "." (punto) y eliminar todas las demás.

Necesito eliminar de cada fichero *-smiles.txt, generado con el script anterior, en la columna de SMILES aquellas cadenas de texto anteriores o posteriores a un punto "." La columna de SMILES tiene información que permite generar la estructura de moléculas en 3D y algunos valores de esta columna tienen más de una molécula, separadas por punto "."; por tanto, el script debe buscar en la ruta D:\MolPort-database_v7, en cada fichero cuyo nombre termine en "-smiles.txt" cadenas de texto que contengan punto "." , elimiminar todas las más pequeñas antes o después del "." y quedarse con la cadena más grande. De salida el script debe generar otro fichero de texto cuyo nombre sea el del fichero de entrada terminado en _sin-puntos.txt.

Ejemplo de nombre de fichero de entrada: fulldb_smiles-000-000-000--000-499-999-smiles.txt y de fichero de salida: fulldb_smiles-000-000-000--000-499-999_sin-puntos.txt. Vease la diferencia...

import os
import pandas as pd

# Ruta de la carpeta donde se encuentran los archivos
input_folder = r"D:\MolPort-database_v7"

# Iterar sobre todos los archivos en la carpeta
for file_name in os.listdir(input_folder):
    if file_name.endswith("-smiles.txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue
        
        # Verificar si la columna "SMILES" existe
        if "SMILES" in data.columns:
            smiles_column = data["SMILES"]
            processed_smiles = []
            removed_smiles = []

            for entry in smiles_column:
                if "." in entry:
                    # Dividir las cadenas por el punto y seleccionar la más larga
                    parts = entry.split(".")
                    largest_part = max(parts, key=len)
                    removed_parts = [part for part in parts if part != largest_part]
                    processed_smiles.append(largest_part)
                    removed_smiles.append(".".join(removed_parts))
                else:
                    processed_smiles.append(entry)
                    removed_smiles.append("")

            # Crear los DataFrames de salida
            data["SMILES"] = processed_smiles
            removed_data = pd.DataFrame({"SMILES Eliminado": removed_smiles})

            # Generar nombres de archivo de salida
            output_processed_file = os.path.join(input_folder, file_name.replace("-smiles.txt", "_sin-puntos.txt"))
            output_removed_file = os.path.join(input_folder, file_name.replace("-smiles.txt", "_eliminado.txt"))

            # Guardar los archivos de salida
            try:
                data.to_csv(output_processed_file, sep="\t", index=False)
                removed_data.to_csv(output_removed_file, sep="\t", index=False)
                print(f"Archivos generados: {output_processed_file}, {output_removed_file}")
            except Exception as e:
                print(f"Error al escribir archivos de salida para {file_name}: {e}")
        else:
            print(f"El archivo {file_name} no contiene la columna 'SMILES'.")

Necesitamos filtar esta libreria "monstruosa" y reducir el número de moléculas disponibles para simulaciones de molecular docking y para ello vamos a usar varios desciptores de esas moléculas que nos permitan descartar aquellas moléculas que no cumplan los rangos de valores de esos descriptores. Podemos tener en cuenta las reglas de Lipinski.

El siguiente script de Python, para ejecutar en Windows, busca en la ruta D:\MolPort-database_v7 todos los ficheros de texto cuyo nombre termine en _sin-puntos.txt y en los que en la primera columna encuentre SMILES y en la segunda MOLPORTID, para calcular el peso molecular de la molécula a la que da lugar cada smiles, los valores de logP, topological polar surface area: TPSA y número de enlaces rotables: RotatableBonds. Una vez calculados estos parámetros los colocará en columnas sucesivas, generando un nuevo fichero de salida cuyo nombre sea el del fichero de entrada añadiendo _MW_logP_TPSA_RotatableBonds.txt.

import os
import pandas as pd
from rdkit import Chem
from rdkit.Chem import Descriptors

# Ruta de la carpeta donde se encuentran los archivos
input_folder = r"D:\MolPort-database_v7"

# Iterar sobre todos los archivos en la carpeta
for file_name in os.listdir(input_folder):
    if file_name.endswith("_sin-puntos.txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue
        
        # Verificar si las columnas "SMILES" y "MOLPORTID" existen
        if "SMILES" in data.columns and "MOLPORTID" in data.columns:
            molecular_weights = []
            logP_values = []
            tpsa_values = []
            # num_atoms_values = []
            rotatable_bonds_values = []

            for smiles in data["SMILES"]:
                try:
                    # Convertir SMILES a una molécula RDKit
                    mol = Chem.MolFromSmiles(smiles)
                    if mol:
                        # Calcular descriptores
                        mw = Descriptors.MolWt(mol)
                        logP = Descriptors.MolLogP(mol)
                        tpsa = Descriptors.TPSA(mol)
                        rotatable_bonds = Descriptors.NumRotatableBonds(mol)
                    else:
                        mw = logP = tpsa = rotatable_bonds = None
                except Exception as e:
                    print(f"Error al procesar SMILES {smiles} en archivo {file_name}: {e}")
                    mw = logP = tpsa = rotatable_bonds = None

                molecular_weights.append(mw)
                logP_values.append(logP)
                tpsa_values.append(tpsa)
                rotatable_bonds_values.append(rotatable_bonds)

            # Agregar las columnas calculadas al DataFrame
            data["MW"] = molecular_weights
            data["logP"] = logP_values
            data["TPSA"] = tpsa_values
            data["RotatableBonds"] = rotatable_bonds_values

            # Generar nombre de archivo de salida
            output_file_path = os.path.join(input_folder, file_name.replace("_sin-puntos.txt", "_MW_logP_TPSA_RotatableBonds.txt"))
            
            # Guardar el archivo de salida
            try:
                data.to_csv(output_file_path, sep="\t", index=False)
                print(f"Archivo generado: {output_file_path}")
            except Exception as e:
                print(f"Error al escribir el archivo de salida para {file_name}: {e}")
        else:
            print(f"El archivo {file_name} no contiene las columnas necesarias 'SMILES' y 'MOLPORTID'.")

Cuando ya disponga de los ficheros de texto con los cuatro parámetros calculados: MW, logP, TPSA y RotatableBonds usaré el siguiente script de Python para filtrar compuestos con MW entre 300-750, logP menor o igual a 3, TPSA menor o igual a 140 y RotatableBonds menor o igual a 7.

El siguiente script de Python para ejecutar en Windows, lee todos los ficheros de texto en la ruta D:\MolPort-database_v7, cuyo nombre termine en "_MW_logP_TPSA_RotatableBonds.txt", que contienen seis columnas y en la primera columna tienen el valor de SMILES, en la tercera columna el peso molecular (MW), en la cuarta columna tiene el valor de logP, en la quinta columna tiene el valor de TPSA y detecta las moléculas con peso molecular mayor o igual a 300 y menor de 750, que además tengan un valor de logP menor o igual a 3, un valor de TPSA menor de 140 y 7 o menos enlaces rotables; generando un fichero de salida cuyo nombre sea el nombre el fichero de entrada sustituyendo "_MW_logP_TPSA_RotatableBonds.txt" por "_300-750_menor3_menor140_RotBonds7.txt"

import os
import pandas as pd

# Ruta de la carpeta donde se encuentran los archivos
input_folder = r"D:\MolPort-database_v7"

# Iterar sobre todos los archivos en la carpeta
for file_name in os.listdir(input_folder):
    if file_name.endswith("_MW_logP_TPSA_RotatableBonds.txt"):
        input_file_path = os.path.join(input_folder, file_name)

        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue

        # Verificar si las columnas necesarias existen
        required_columns = ["SMILES", "MW", "logP", "TPSA", "RotatableBonds"]
        if all(column in data.columns for column in required_columns):
            try:
                # Aplicar filtros
                filtered_data = data[
                    (data["MW"] >= 300) & (data["MW"] < 750) &
                    (data["logP"] <= 3) &
                    (data["TPSA"] < 140) &
                    (data["RotatableBonds"] <= 7)
                ]

                # Verificar si hay datos filtrados
                if not filtered_data.empty:
                    # Generar nombre de archivo de salida
                    output_file_path = os.path.join(
                        input_folder, 
                        file_name.replace("_MW_logP_TPSA_RotatableBonds.txt", "_300-750_menor3_menor140_RotBonds7.txt")
                    )

                    # Guardar el archivo de salida
                    filtered_data.to_csv(output_file_path, sep="\t", index=False)
                    print(f"Archivo generado: {output_file_path}")
                else:
                    print(f"No se encontraron moléculas que cumplan los criterios en {file_name}.")

            except Exception as e:
                print(f"Error al filtrar o guardar datos para {file_name}: {e}")
        else:
            print(f"El archivo {file_name} no contiene las columnas necesarias {required_columns}.")

Las moléculas con carbonos asimétricos (4 sustituyentes diferentes o sea, con centros quirales) son un "dolor" en química médica... voy a separar las moléculas que llevo filtradas hasta este punto del proceso en "quirales" y "no quirales", de las primeras mi base de datos sólo contiene un estereosiómero de todos los posibles.

El siguiente script de Python para ejecutar en Windows, lee todos los ficheros de texto en la ruta D:\MolPort-database_v7 cuyo nombre termine en "_300-750_menor3_menor140_RotBonds7.txt" que contienen seis columnas y en la primera columna tienen el valor de SMILES, en la tercera columna el peso molecular (MW), en la cuarta columna tiene el valor de logP, en la quinta columna tiene el valor de TPSA y detecta las moléculas no quirales; generando un fichero de salida cuyo nombre sea el nombre el fichero de entrada en el que se sustituye "_300-750_menor3_menor140_RotBonds7.txt" por "_quiral.txt".

import os
import pandas as pd
from rdkit import Chem

# Ruta de entrada
directory = r"D:\\MolPort-database_v7"

# Función para detectar carbonos quirales
def detect_chiral(smiles):
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol:
            chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)
            return len(chiral_centers) > 0
        return False
    except Exception as e:
        print(f"Error procesando SMILES: {smiles}, {e}")
        return False

# Procesar los archivos en la carpeta
for filename in os.listdir(directory):
    if filename.endswith("_300-750_menor3_menor140_RotBonds7.txt"):
        filepath = os.path.join(directory, filename)

        # Leer el archivo en un DataFrame
        try:
            data = pd.read_csv(filepath, sep="\t", header=0)

            # Comprobar que tiene las columnas esperadas
            if "SMILES" in data.columns and "MOLPORTID" in data.columns:

                # Filtrar moléculas quirales
                chiral_molecules = data[data["SMILES"].apply(detect_chiral)]

                # Generar el nombre del archivo de salida
                output_filepath = os.path.join(directory, filename.replace("_300-750_menor3_menor140_RotBonds7.txt", "_quiral.txt"))

                # Guardar el resultado en un nuevo archivo
                if not chiral_molecules.empty:
                    chiral_molecules.to_csv(output_filepath, sep="\t", index=False)
                    print(f"Archivo generado: {output_filepath}")
                else:
                    print(f"No se encontraron moléculas quirales en: {filename}")
            else:
                print(f"Archivo {filename} no contiene las columnas esperadas.")

        except Exception as e:
            print(f"Error procesando el archivo {filename}: {e}")

El script equivalente al anterior para detectar las moléculas no quirales.

import os
import pandas as pd
from rdkit import Chem

# Ruta de la carpeta donde se encuentran los archivos
input_folder = r"D:\MolPort-database_v7"

# Iterar sobre todos los archivos en la carpeta
for file_name in os.listdir(input_folder):
    if file_name.endswith("_300-750_menor3_menor140_RotBonds7.txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue

        # Verificar si la columna "SMILES" existe
        if "SMILES" in data.columns:
            non_chiral_rows = []

            for _, row in data.iterrows():
                smiles = row["SMILES"]

                try:
                    # Convertir SMILES a una molécula RDKit
                    mol = Chem.MolFromSmiles(smiles)
                    if mol:
                        # Detectar si no tiene carbonos quirales
                        if not Chem.FindMolChiralCenters(mol, includeUnassigned=True):
                            non_chiral_rows.append(row)
                except Exception as e:
                    print(f"Error al procesar SMILES {smiles} en archivo {file_name}: {e}")

            # Crear un DataFrame con las filas no quirales
            non_chiral_data = pd.DataFrame(non_chiral_rows)

            # Generar el nombre del archivo de salida
            output_file_path = os.path.join(input_folder, file_name.replace("_300-750_menor3_menor140_RotBonds7.txt", "_NO-quirales.txt"))
            
            # Guardar el archivo de salida
            try:
                non_chiral_data.to_csv(output_file_path, sep="\t", index=False)
                print(f"Archivo generado: {output_file_path}")
            except Exception as e:
                print(f"Error al escribir el archivo de salida para {file_name}: {e}")
        else:
            print(f"El archivo {file_name} no contiene la columna necesaria 'SMILES'.")

Para hacer simulaciones de acoplamiento molecular (docking) usando el software de YASARA necesitamos ficheros SDF que no contengan más de 3000 moléculas (pueden usarse ficheros con mayor número de moléculas pero la probabilidad de que el proceso se rompa es alta).

El siguiente script de Python lee en la ruta D:\MolPort-database_v7\07__filter_300-750_menor3_menor140_RotBonds7__kirales todos los ficheros cuyo nombre termine en _quiral.txt los divide en tantos ficheros de 3000 moléculas como sea necesario y los guarda, en la ruta D:\MolPort-database_v7\09_MP7_quiral_3000\, con las columnas SMILES y MOLPORTID.

import os
import pandas as pd

# Ruta de entrada y salida
input_folder = r"D:\MolPort-database_v7\07__filter_300-750_menor3_menor140_RotBonds7__kirales"
output_folder = r"D:\MolPort-database_v7\09_MP7_quiral_3000"

# Crear la carpeta de salida si no existe
os.makedirs(output_folder, exist_ok=True)

# Límite de moléculas por archivo
LIMIT = 3000

# Contador para los archivos de salida
file_counter = 1

# Buffer para acumular datos
accumulated_data = []

# Formatear el nombre del archivo de salida
def generate_output_filename(counter):
    return os.path.join(output_folder, f"MP7-Q-{counter:04d}.txt")

# Función para guardar el buffer en un archivo
def save_to_file(data, counter):
    output_file = generate_output_filename(counter)
    with open(output_file, "w") as f:
        for row in data:
            f.write("\t".join(row) + "\n")
    print(f"Archivo generado: {output_file}")

# Procesar los archivos de entrada
for file_name in os.listdir(input_folder):
    if file_name.endswith("_quiral.txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue

        # Verificar si las columnas SMILES y MOLPORTID existen
        if "SMILES" in data.columns and "MOLPORTID" in data.columns:
            for _, row in data.iterrows():
                accumulated_data.append((row["SMILES"], row["MOLPORTID"]))

                # Si alcanzamos el límite, guardar y resetear el buffer
                if len(accumulated_data) >= LIMIT:
                    save_to_file(accumulated_data, file_counter)
                    file_counter += 1
                    accumulated_data = []
        else:
            print(f"El archivo {file_name} no contiene las columnas necesarias 'SMILES' y 'MOLPORTID'.")

# Guardar los datos restantes si hay alguno
if accumulated_data:
    save_to_file(accumulated_data, file_counter)

El siguiente script de Python lee en la ruta D:\MolPort-database_v7\08__filter_300-750_menor3_menor140_RotBonds7__NO-kirales todos los ficheros cuyo nombre termine en _NO-quirales.txt los divide en tantos ficheros de 3000 moléculas como sea necesario y los guarda, en la ruta D:\MolPort-database_v7\10_MP7_no-quiral_3000\, con las columnas SMILES y MOLPORTID.

import os
import pandas as pd

# Ruta de entrada y salida
input_folder = r"D:\MolPort-database_v7\08__filter_300-750_menor3_menor140_RotBonds7__NO-kirales"
output_folder = r"D:\MolPort-database_v7\10_MP7_no-quiral_3000"

# Crear la carpeta de salida si no existe
os.makedirs(output_folder, exist_ok=True)

# Límite de moléculas por archivo
LIMIT = 3000

# Contador para los archivos de salida
file_counter = 1

# Buffer para acumular datos
accumulated_data = []

# Formatear el nombre del archivo de salida
def generate_output_filename(counter):
    return os.path.join(output_folder, f"MP7-noQ-{counter:04d}.txt")

# Función para guardar el buffer en un archivo
def save_to_file(data, counter):
    output_file = generate_output_filename(counter)
    with open(output_file, "w") as f:
        for row in data:
            f.write("\t".join(row) + "\n")
    print(f"Archivo generado: {output_file}")

# Procesar los archivos de entrada
for file_name in os.listdir(input_folder):
    if file_name.endswith("_NO-quirales.txt"):
        input_file_path = os.path.join(input_folder, file_name)
        
        try:
            # Leer el archivo de entrada
            data = pd.read_csv(input_file_path, sep="\t")
        except Exception as e:
            print(f"Error al leer el archivo {file_name}: {e}")
            continue

        # Verificar si las columnas SMILES y MOLPORTID existen
        if "SMILES" in data.columns and "MOLPORTID" in data.columns:
            for _, row in data.iterrows():
                accumulated_data.append((row["SMILES"], row["MOLPORTID"]))

                # Si alcanzamos el límite, guardar y resetear el buffer
                if len(accumulated_data) >= LIMIT:
                    save_to_file(accumulated_data, file_counter)
                    file_counter += 1
                    accumulated_data = []
        else:
            print(f"El archivo {file_name} no contiene las columnas necesarias 'SMILES' y 'MOLPORTID'.")

# Guardar los datos restantes si hay alguno
if accumulated_data:
    save_to_file(accumulated_data, file_counter)

Ahora viene cuando observamos los números de vértigo: en la sección quirales tenemos 1678 ficheros de 3000 moléculas cada uno y 1 fichero con 2418 moléculas, todos ellos suman 5.0336.418 moléculas. En la sección no quirales tenemos 1633 ficheros de 3000 moléculas cada uno y 1 fichero con 1062 moléculas, todos ellos suman 4.900.062 moléculas.

Nuestra última tarea consiste en convertir esas 5.0336.418 + 4.900.062 moléculas en estructuras 3D protonadas, por ejemplo a pH 7.4 usando un script de YASARA.

libname='MP7-Q-0001'
Clear
for line in file (libname).txt
  # Extraer SMILES y nombre
  smiles,name=split line
  # Construir SMILES
  obj = BuildSMILES '(smiles)'
  # Minimizar con NOVA para cerrar los enlaces largos que BuildSMILES dejó
  #  (La optimización QM no sabe dónde deben estar los enlaces)
  Cell Auto,Extension=10
  ph=7.4
  ForceField NOVA,SetPar=yes
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  # Optimizar con MOPAC para afinar la geometría
  OptimizeObj (obj)
  # Establecer el nombre del compuesto
  CompoundMol Obj (obj),'(name)'
  RemoveObj (obj)
# Guardar un único archivo SDF
AddObj all
SaveSDF !SimCell,(libname)

Pero YASARA es un poco lento ... y me parece mejor idea usar RDkit para conseguir convertir los ficheros de texto con SMILES y MOLPORTID en ficheros SDF con moléculas en 3D y utilizar MMFF94 como método de optimización geométrica. Para eso sirve el siguiente script de Python de ejecución en Windows y con los módulos pandas, numpy, RDkit y tqdm. El script descarta en la salida moléculas que contienen átomos no soportados antes de optimizar la geometría. Sólo las moléculas válidas se procesarán normalmente con hidrógenos añadidos y geometrías optimizadas.

import os
from rdkit import Chem
from rdkit.Chem import AllChem
from tqdm import tqdm

# Ruta de entrada y salida
ruta_entrada = r"D:\aaa"
ruta_salida = r"D:\aaa\salida_sdf"

# Crear la carpeta de salida si no existe
os.makedirs(ruta_salida, exist_ok=True)

# Función para optimizar geometría con MMFF94
def optimizar_con_mmff94(molecula):
    try:
        # Generar propiedades del campo de fuerza MMFF
        props = AllChem.MMFFGetMoleculeProperties(molecula)
        if props is None:
            return None  # Si no se pueden generar propiedades, devolver None
        
        # Crear el campo de fuerza y optimizar
        ff = AllChem.MMFFGetMoleculeForceField(molecula, props)
        if ff is not None:
            ff.Minimize()
            return molecula
        return None
    except Exception as e:
        print(f"Error durante la optimización con MMFF94: {e}")
        return None

# Iterar sobre los archivos de la ruta de entrada
for archivo in os.listdir(ruta_entrada):
    if archivo.endswith(".txt"):
        ruta_archivo = os.path.join(ruta_entrada, archivo)
        
        with open(ruta_archivo, 'r') as f:
            lineas = f.readlines()
        
        # Crear una lista para almacenar moléculas optimizadas
        moleculas_optimizadas = []
        
        print(f"Procesando archivo: {archivo}")
        
        # Crear barra de progreso
        for linea in tqdm(lineas, desc=f"Procesando moléculas en {archivo}", unit="mol"):
            try:
                smiles, nombre = linea.strip().split('\t')
                molecula = Chem.MolFromSmiles(smiles)
                if molecula is not None:
                    # Añadir hidrógenos
                    molecula = Chem.AddHs(molecula)
                    
                    # Generar coordenadas 3D iniciales
                    if AllChem.EmbedMolecule(molecula, randomSeed=42) == 0:
                        # Optimizar con MMFF94
                        molecula_optimizada = optimizar_con_mmff94(molecula)
                        if molecula_optimizada is not None:
                            molecula_optimizada.SetProp("_Name", nombre)
                            moleculas_optimizadas.append(molecula_optimizada)
                else:
                    print(f"No se pudo generar la molécula a partir de SMILES: {smiles}")
            except Exception as e:
                print(f"Error procesando la molécula {nombre} ({smiles}): {e}")
                continue
        
        # Escribir las moléculas optimizadas en formato SDF
        if moleculas_optimizadas:
            nombre_salida = os.path.splitext(archivo)[0] + ".sdf"
            ruta_salida_archivo = os.path.join(ruta_salida, nombre_salida)
            with Chem.SDWriter(ruta_salida_archivo) as escritor:
                for mol in moleculas_optimizadas:
                    escritor.write(mol)
        
        print(f"Procesado: {archivo}, moléculas optimizadas guardadas en {ruta_salida_archivo}")

El siguiente script de Python, para ejecutar en Linux, usa una lista contenida en la ruta /home/jant.encinar/py-linux/36_smiles-3D-sdf/lista_01.txt, con el nombre de los ficheros de texto en columnas separadas por tabulador, que contienen en la primera columna SMILES y en la segunda MOLPORTID, localizados en la misma ruta, para realizar la operación de conversión de SMILES a ficheros SDF con moléculas en 3D y utilizar MMFF94 como método de optimización geométrica.

# -*- coding: iso-8859-1 -*-
import os
from rdkit import Chem
from rdkit.Chem import AllChem
from tqdm import tqdm

# Rutas de entrada y salida
ruta_lista_archivos = "/home/jant.encinar/py-linux/36_smiles-3D-sdf/lista_01.txt"
ruta_salida = "/home/jant.encinar/py-linux/36_smiles-3D-sdf/salida_sdf"

# Crear la carpeta de salida si no existe
os.makedirs(ruta_salida, exist_ok=True)

# Función para optimizar geometría con MMFF94
def optimizar_con_mmff94(molecula):
    try:
        # Generar propiedades del campo de fuerza MMFF
        props = AllChem.MMFFGetMoleculeProperties(molecula)
        if props is None:
            return None  # Si no se pueden generar propiedades, devolver None
        
        # Crear el campo de fuerza y optimizar
        ff = AllChem.MMFFGetMoleculeForceField(molecula, props)
        if ff is not None:
            ff.Minimize()
            return molecula
        return None
    except Exception as e:
        print(f"Error durante la optimización con MMFF94: {e}")
        return None

# Leer la lista de archivos
try:
    with open(ruta_lista_archivos, 'r') as f:
        lista_archivos = [linea.strip() for linea in f.readlines() if linea.strip()]
except FileNotFoundError:
    print(f"No se encontró el archivo de lista: {ruta_lista_archivos}")
    exit(1)

# Procesar cada archivo en la lista
for archivo in lista_archivos:
    ruta_archivo = f"/home/jant.encinar/py-linux/36_smiles-3D-sdf/{archivo}.txt"

    if not os.path.exists(ruta_archivo):
        print(f"Archivo no encontrado: {ruta_archivo}")
        continue

    with open(ruta_archivo, 'r') as f:
        lineas = f.readlines()

    # Crear una lista para almacenar moléculas optimizadas
    moleculas_optimizadas = []

    print(f"Procesando archivo: {archivo}.txt")

    # Crear barra de progreso
    for linea in tqdm(lineas, desc=f"Procesando moléculas en {archivo}.txt", unit="mol"):
        try:
            smiles, nombre = linea.strip().split('\t')
            molecula = Chem.MolFromSmiles(smiles)
            if molecula is not None:
                # Añadir hidrógenos
                molecula = Chem.AddHs(molecula)

                # Generar coordenadas 3D iniciales
                if AllChem.EmbedMolecule(molecula, randomSeed=42) == 0:
                    # Optimizar con MMFF94
                    molecula_optimizada = optimizar_con_mmff94(molecula)
                    if molecula_optimizada is not None:
                        molecula_optimizada.SetProp("_Name", nombre)
                        moleculas_optimizadas.append(molecula_optimizada)
            else:
                print(f"No se pudo generar la molécula a partir de SMILES: {smiles}")
        except Exception as e:
            print(f"Error procesando la molécula {nombre} ({smiles}): {e}")
            continue

    # Escribir las moléculas optimizadas en formato SDF
    if moleculas_optimizadas:
        nombre_salida = f"{archivo}.sdf"
        ruta_salida_archivo = os.path.join(ruta_salida, nombre_salida)
        with Chem.SDWriter(ruta_salida_archivo) as escritor:
            for mol in moleculas_optimizadas:
                escritor.write(mol)

        print(f"Procesado: {archivo}.txt, moléculas optimizadas guardadas en {ruta_salida_archivo}")

Antes de usar el script anterior, activa el entorno virtual (en mi caso se llama jant-molport) y ejecuta el script verifica.py para asegurarte de que tienen instalado pandas, numpy, tqdm y rdkit.

(base) [jant.encinar@castleblack 36_smiles-3D-sdf]$ conda install -c conda-forge rdkit
(base) [jant.encinar@castleblack 36_smiles-3D-sdf]$ conda update -n base -c conda-forge conda
(base) [jant.encinar@castleblack 36_smiles-3D-sdf]$ conda install -c conda-forge tqdm
(base) [jant.encinar@castleblack 36_smiles-3D-sdf]$ conda install pandas numpy



(base) [jant.encinar@castleblack 36_smiles-3D-sdf]$ conda activate jant-molport
(jant-molport) [jant.encinar@castleblack 36_smiles-3D-sdf]$ python3 verifica.py
Version instalada de pandas:
2.2.3
Version instalada de numpy:
2.2.0
Version instalada de tqdm:
4.67.1
Version instalada de RDkit:
2024.09.3
(jant-molport) [jant.encinar@castleblack 36_smiles-3D-sdf]$

Script verifica.py justo debajo...

import pandas as pd
import numpy as np
import tqdm

print("Version instalada de pandas")
print(pd.__version__)
print("Version instalada de numpy")
print(np.__version__)
print("Version instalada de tqdm")
print(tqdm.__version__)
from rdkit import rdBase
print("Version instalada de RDkit")
print(rdBase.rdkitVersion)

Con el siguiente script de Python para ejecutar en Windows vamos a contar el número de moléculas de nuestra libreria, tanto para cada fichero *.sdf como el total.

import os
from tqdm import tqdm  # Para la barra de progreso

# Ruta donde se encuentran los archivos *.sdf
ruta_sdf = r"D:\MolPort-database_v7\11_MP7_quiral_sdf"

# Archivo de salida
archivo_salida = r"D:\MolPort-database_v7\11_MP7_quiral_sdf\_lista-sdf-quiral.txt"

def contar_moleculas_en_sdf(ruta_archivo):
    """Cuenta el número de moléculas en un archivo .sdf usando $$$$ como delimitador."""
    try:
        with open(ruta_archivo, 'r') as archivo:
            contenido = archivo.read()
        # Contar ocurrencias del delimitador $$$$
        return contenido.count("$$$$")
    except Exception as e:
        print(f"Error leyendo {ruta_archivo}: {e}")
        return 0

def procesar_archivos_sdf(ruta, salida):
    """Procesa todos los archivos .sdf en la ruta dada y genera un reporte."""
    try:
        # Lista de archivos en la ruta
        archivos_sdf = [f for f in os.listdir(ruta) if f.endswith('.sdf')]
        total_archivos = len(archivos_sdf)

        # Verificar si hay archivos para procesar
        if total_archivos == 0:
            print("No se encontraron archivos .sdf en la ruta especificada.")
            return

        # Inicializar contador de moléculas totales
        total_moleculas = 0

        # Abrir archivo de salida para escribir resultados
        with open(salida, 'w') as salida_file:
            # Escribir encabezado
            salida_file.write("Fichero\tNumero_de_moleculas\n")
            
            # Procesar cada archivo .sdf con barra de progreso
            for archivo in tqdm(archivos_sdf, desc="Procesando archivos", unit="archivo"):
                ruta_completa = os.path.join(ruta, archivo)
                num_moleculas = contar_moleculas_en_sdf(ruta_completa)
                total_moleculas += num_moleculas
                salida_file.write(f"{archivo}\t{num_moleculas}\n")
            
            # Escribir total de moléculas al final del archivo
            salida_file.write(f"\nTOTAL\t{total_moleculas}\n")

        # Mostrar total en consola
        print(f"\nProcesamiento completado. Total de moléculas: {total_moleculas}")
        print(f"Reporte generado en: {salida}")
    
    except Exception as e:
        print(f"Error procesando archivos en {ruta}: {e}")

# Ejecutar la función principal
if __name__ == "__main__":
    procesar_archivos_sdf(ruta_sdf, archivo_salida)

Puedes descargar una copia de la base de datos MolPort v7 (5,73 GB).