Spectres RMN avec Python

Il est difficile de trouver des données de spectres RMN. Contrairement à des bases de données de spectres UV-visible ou infrarouge, certains sites affichent un spectre uniquement sous forme de fichiers images, dont la mauvaise résolution limite l’insertion dans un document pédagogique.

J’ai trouvé deux sites de bases de données pour spectres RMN :

• https://np-mrd.org/classyfication
• https://sdbs.db.aist.go.jp/sdbs/cgi-bin/direct_frame_top.cgi

Le script Python proposé ici permet d’insérer les données numériques disponibles sur ces sites, afin de tracer des spectres personnalisables.

Récupération de données chimiques de spectre RMN

Commençons par le site https://sdbs.db.aist.go.jp/sdbs/cgi-bin/direct_frame_top.cgi. Cherchons l’acide ascorbique.

On trouve plusieurs résultats de recherche. La formule brute peut servir de confirmation du produit recherché. En cliquant ensuite sur la ligne du produit souhaité, on visualise le spectre, ainsi que la formule détaillée de la molécule, avec l’assignation des protons aux signaux.

L’image n’est pas personnalisable (pas de zoom possible sur les axes). On a parfois du mal à distinguer la multiplicité des signaux. Pour aller plus loin, il faut cliquer sur le bouton gris “peak data” sous la formule :

On a les colonnes de donnés “ppm” (abscisses des pics), Hz (fréquence) et “Int.” (intensité).

Pour réaliser un spectre personnalisé, on copie-colle ces données dans un fichier de texte brut, avec un logiciel très simple, comme il y en a sur tous les systèmes d’exploitation, comme Gedit (linux), Bloc-Notes (Windows) ou TextEdit (MacOS) :

Pour que le script Python retrouve facilement les données de ce fichier, il faut d’une part faire quelques choix définitifs pour les noms des colonnes (autant garder ceux utilisés par le site internet), et d’autre part supprimer tout caractère inutile ou qui pourrait empêcher la lecture des données. Ici, par exemple, il faut supprimer tous les espaces doubles séparant les colonnes, et les remplacer par un caractère unique (virgule, ou tabulation). Le logiciel de texte brut a un outil “rechercher/remplacer” assez efficace.

Enregistrer ensuite le fichier en lui donnant le nom de la molécule, et avec l’extension .csv

Création du graphique avec le script Python

À présent, le script Python va faire la suite, et mettre en forme les données dans un spectre RMN personnalisable. Ce script utilise les modules suivants : numpy (gestion de listes et tableaux de données), matplotlib (création de graphiques) et pandas (gestion de tableaux de données).

Voici les grandes lignes du programme :

• importer les données trouvées sur le site web, à partir du fichier .csv créé précédemment.
• plutôt que de créer un graphique avec des pics fins verticaux selon les coordonnées du fichier, l’idée est de tracer pour chaque pic une courbe de Gauss très fine, centrée sur l’abscisse de chaque ligne du fichier. La courbe de Gauss est du type I*exp(-a*(x-x0)^2) où x0 est l’abscisse d’un pic et I son ordonnée. Le paramètre “a” détermine la largeur de la courbe de Gauss. Avec a = 5000 ou 10000, les courbes de Gauss simulent assez bien des pics de signaux RMN.

Voici le résultat obtenu : le spectre comporte bien l’essentiel des données spectroscopiques, mais il est personnalisable (zoom, grille secondaire…).

Les données du fichier sont importées grâce au module pandas dans un DataFrame appelé “spectre”. C’est un tableau comportant un pic par ligne, avec en colonnes les abscisses (ppm) et les intensités (Int.).

Le script crée une liste numpy “X” de valeurs rapprochées entre le minimum et le maximum défini par l’utilisateur. Puis une liste numpy “Y” comportant uniquement des zéros, et de même longueur que X. Le module numpy permet très facilement ces créations de listes. La liste Y sert de ligne de base pour la courbe définitive.

Ensuite, une boucle for parcourt les ligne du tableau “spectre” ; pour chaque ligne (pic), on calcule les ordonnées Y de la courbe de Gauss, et on les ajoute à la liste Y précédente. La ligne de base est donc augmentée avec chaque pic. À l’issue de la boucle for, la liste Y comporte les données de tous les pics.

Il ne reste plus qu’à tracer Y en fonction de X dans le graphique créé avec matplotlib.

fichiers du spectre RMN de l’acide ascorbique :

Autre “mode” du script : création de multiplets

Prenons un autre exemple sur le site https://np-mrd.org/classyfication. Le choix s’est porté sur la lysine, un acide aminé.

On trouve la L-lysine dans la base de données. Il faut cliquer sur le nom. On arrive à une page présentant la molécule, avec notamment une formule topologique. En haut, un menu sur fond vert permet de cliquer sur spectra.

La page propose alors un ou plusieurs spectres (RMN C, RMN H, simulés ou expérimentaux, parfois avec un choix de conditions expérimentales). En cliquant sur l’un des spectres proposés, on arrive sur la page suivante :

Comme on peut le voir, le spectre est peu présentable, même après avoir zoomé sur l’axe horizontal. Cependant, le tableau fournit toutes les données nécessaires pour tracer un spectre personnalisé :

• l’abscisse du centre du signal multiplet
• le type de multiplet (s :singulet, d : doublet, etc.)
• le nombre d’atomes H du signal.

On crée alors un fichier texte brut. Nommons ce fichier du nom de la molécule (sans accents).

Dans une colonne “ppm” on copie l’abscisse du centre du multiplet. Dans la colonne “mult” : une lettre pour le multiplet correspondant, enfin dans la colonne “Int.” l’intensité du signal en nombre d’atomes H.

Le script Python comporte un autre mode de création de courbe : si les données comportent une colonne “mult”, alors la lettre ‘s’, ‘t’, ‘q’ pour singulet, triplet, quadruplet est utilisée dans une fonction multiplet qui crée automatiquement les pics de la forme voulue, et renvoie les valeurs de Y correspondantes. Pour la lysine, avec les données du fichier précédent, cela donne :

un triplet à 3.74 ppm, un autre à 3.02 ppm, et trois singulets à 1.89, 1.71 et 1.46 ppm.

J’ai essayé de faire en sorte que la hauteur des pics correspondent au nombre de H du signal (à une constante multiplicative près). J’ai pris en compte le fait que dans une gaussienne bornée en y à une valeur I, l’aire sous la courbe vaut I (la constante multiplicative est sigma, écart-type, fois racine de 2 pi. Pour les multiplets, j’ai tenu compte des termes du triangle de Pascal pour calculer la hauteur max du pic le plus faible. J’espère ne pas m’être trompé !

Une dernière fonctionnalité du script : une fonction peut être appelée, qui crée un DataFrame contenant les valeurs X et Y de la courbe tracée, et les exporte dans un fichier .csv. Cela est utile pour réutiliser ces données dans un document LaTeX.

Le script à télécharger :

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Créé le Wed May 31 22:07:04 2023

@auteur: David ALBERTO
(www.astrolabe-science.fr)
Ce script trace des spectres RMN, soit d'après des valeurs expérimentales
trouvées sur internet, soit des spectres simulés.
"""
import matplotlib.pyplot as plt # graphiques
import numpy as np # listes
import pandas as pd # tableaux de données et import/export


# plt.rcParams["font.family"] = "Roboto" #ou autre police installée
plt.rcParams["font.size"] = 6

molecule = 'L-lysine'
xmin, xmax = 0.5, 4 # à ajuster aux données

NomFichieBrut = molecule + '.csv'
NomFichierExport = molecule + '_continu.csv'

prop_ligne = dict(lw=0.9, c= 'C0')# mise en forme de la ligne de la courbe

spectre = pd.read_csv(molecule+'.csv', sep=',')# import des données

ppm = spectre['ppm']
if 'mult' in spectre.columns:
    mult = spectre['mult']
H = spectre['Int.']

delta = 0.06 # écart entre les pics, pour les multiplets modélisés
a = 10000 # inverse de la largeur des pics (+ grand, + fin)

def multiplet(m, ppm, H):
    """
    m : str (multiplicité)
    H : (entier) : intensité du signal (nombre de protons)
    ppm (flottant): abscisse du signal
    renvoie une liste de valeurs Y en ordonnées
    I est l'ordonnée du pic le + faible du multiplet
    """
    if m=='d': # doublet
        I = H / 2
        Y = I * np.exp( -a*(X-ppm-delta)**2) + I * np.exp( -a*(X-ppm+delta)**2)
        return Y
    elif m=='t': # triplet
        I = H / 4
        Y = I * (np.exp( -a*(X-ppm-delta)**2) + 2 * np.exp( -a*(X-ppm)**2) + np.exp( -a*(X-ppm+delta)**2))
        return Y
    elif m=='q': # quadruplet
        I = H / 8
        Y = I * (np.exp(-a*(X-ppm-1.5*delta)**2) + 2 * np.exp(-a*(X-ppm-0.5*delta)**2) + 2 * np.exp( -a*(X-ppm+0.5*delta)**2) +  np.exp(-a*(X-ppm+1.5*delta)**2))
        return Y
    elif m=='qt': # quintuplet
        I = H / 16
        Y = I * (np.exp(-a*(X-ppm-2*delta)**2) + 4 * np.exp(-a*(X-ppm-delta)**2) + 6 * np.exp( -a*(X-ppm)**2) +  4 * np.exp(-a*(X-ppm+delta)**2) + np.exp(-a*(X-ppm+2*delta)**2))
        return Y
    elif m=='sx':# sextuplet
        I = H / 32
        Y = I * (np.exp(-a*(X-ppm-2.5*delta)**2) + 5 * np.exp(-a*(X-ppm-1.5*delta)**2) + 10 * np.exp( -a*(X-ppm-0.5*delta)**2) +  10 * np.exp(-a*(X-ppm+0.5*delta)**2) + 5 * np.exp(-a*(X-ppm+1.5*delta)**2) + np.exp(-a*(X-ppm+2.5*delta)**2))
        return Y
    else:
        Y = H*np.exp( -a*(X-ppm)**2)
        return Y

def pic_isole(ppm, H):
    """
    H : (entier) : intensité du signal (nombre de protons)
    ppm (flottant): abscisse du signal
    renvoie une liste de valeurs Y en ordonnées
    I est l'ordonnée du pic le + faible du multiplet
    """
    Y = H * np.exp( -a*(X-ppm)**2)
    return Y


def export_fichier(X,Y):
    """
    crée un DataFrame puis un fichier csv pour exporter les valeurs calculées.
    """
    df = pd.DataFrame({'ppm': X,
                        'I': Y})
    df.to_csv(NomFichierExport, index=False, sep='\t')

fig, ax = plt.subplots()

ax.invert_xaxis()
# ax.set_xlim(xmax,xmin)
ax.set_xlabel(r'$\delta$ (ppm)')
ax.set_yticks([])

X = np.arange(xmin,xmax,0.001)
Y = np.zeros_like(X)

for N in range(len(spectre)):
    x = ppm[N]
    intensite = H[N]
    if 'mult' in spectre.columns:
        m = mult[N]
        Y = Y + multiplet(m, x, intensite)
    else:
        Y = Y + pic_isole(x, intensite)

ax.plot(X,Y, **prop_ligne)
ax.text(0.01, 0.99, molecule, transform=ax.transAxes,va='top',c='firebrick')

export_fichier(X, Y)

fig.savefig('spectreRMN' + molecule + '.png',dpi=200)
fig.savefig('spectreRMN' + molecule + '.pdf')