【ケモインフォマティクス事例】MMPA(Matched Molecular Pair Analysis)で部分構造が類似する化合物を分類する

Matched Molecular Pair Analysis（MMPA）は、薬物設計、医薬品発見、化合物最適化の過程で使用されるコンピューターアシストツールです。このアプローチは、化合物ライブラリー内の化合物間の相互作用や性質の関連性を特定し、新しい医薬品の設計に役立てるために使用されます。

• MMPAの主要な特徴とプロセスは以下の通りです：
  1. 化合物の分解: MMPAの最初のステップは、薬物候補を構成する化合物を小さな部分構造（フラグメントまたはペア）に分解することです。これらの部分構造は、化学的な官能基、環状構造、置換基など、分子の特定の部分を表します。
  2. 部分構造のマッチング: これらの部分構造を異なる化合物間でマッチングし、共通の部分構造（"マッチド・モレキュラー・ペア"または"MMパーツ"と呼ばれる）を特定します。これにより、異なる化合物が共通の要素を持つ場合、その要素の影響を評価できるようになります。
  3. データの収集と分析: MMPAは、これらのMMパーツの化学的な特性や生物学的な活性、物性などのデータを収集し、異なるMMパーツの影響を比較・分析します。これにより、特定のMMパーツが薬物の効力、毒性、ADME（吸収、分布、代謝、排泄）特性にどのように影響を与えるかを理解することができます。
  4. 化合物設計への応用: MMPAの結果を活用して、薬物設計者は特定のMMパーツを修正することに焦点を当て、新しい化合物の設計や既存の化合物の最適化を行います。これにより、特定の化合物の特性を向上させるために必要な変更を予測し、効率的な医薬品開発をサポートします。

実装

まず、分類する化合物のリストをpubchemに保存されているsdfファイルから取得するための関数を準備します。

from bs4 import BeautifulSoup
import requests, traceback, os, re, shutil, pathlib, gzip, math
from tqdm import tqdm

def pubchemUrlList() -> list:
    """
    pubchemの特定のサイトからsdf.gzのurlを取得する。
    """
    url = "https://ftp.ncbi.nlm.nih.gov/pubchem/Compound/CURRENT-Full/SDF/"
    response = requests.get(url)
    soup = BeautifulSoup(response.text, "html.parser")

    urls = []
    for row in soup.find_all("a"):
        href = row.get("href")
        if href.endswith(".sdf.gz"):
            urls.append(url+href)
    return urls

def pubchemModifiedDict() -> dict:
    """
    pubchemの特定のサイトからファイル名とlast_modifiedを取得する。
    """
    url = "https://ftp.ncbi.nlm.nih.gov/pubchem/Compound/CURRENT-Full/SDF/"
    response = requests.get(url)
    soup = BeautifulSoup(response.text, "html.parser")
    
    # textからYYYY-MM-DD HH:MMの形式の文字列を抽出するlambda関数
    pattern = r'\d{4}-\d{2}-\d{2} \d{2}:\d{2}'
    get_date = lambda text: re.search(pattern, text).group() if re.search(pattern, text) else None

    rows = soup.find('pre').text.split('\n')
    return {r.split(' ')[0]:get_date(r) for r in rows}

def pubchemDownloadList() -> list:
    """
    pubchemの特定のサイトからファイル名、url、last_modified、db名を取得する。
    """
    db = "pubchem"
    urls = pubchemUrlList()
    modifids = pubchemModifiedDict()
    return [[os.path.basename(u), u, modifids[os.path.basename(u)], db] for u in urls]

def fetch(url: str, filename:str ='', directory:str ='') -> bool:
    """
    url先のファイルをダウンロードします。

    Parameters
    ----------
    url : str
    ダウンロードしたいファイルのURL

    filename : str, default ''
    保存するファイルの名称。指定していない場合、urlの末尾になる。

    directory : str, default ''
    保存先のディレクトリ。指定しない場合、プログラムファイルと同一の階層に保存。

    Returns
    -------
    bool
    fetchに成功するとTrue。失敗するとFalseを返す。
    """
    try:
        response = requests.head(url, allow_redirects=True)
        size = response.headers.get('content-length', -1)
        print('size', convertSize(int(size)))
        if filename == '':
          filename=os.path.basename(url)

        save_name = filename
        if str(directory) != '':
          save_name = os.path.join(str(directory), filename)
        if os.path.exists(save_name):
          print('already exists')
          return save_name
        r = requests.get(url, stream=True)
        with open(save_name, 'wb') as f:
          with tqdm(total=int(size), desc='download', unit='B', unit_scale=True, leave=False, ncols=100) as pbar:
            for chunk in r.iter_content(chunk_size=1024):
              if chunk:
                f.write(chunk)
                f.flush()
                pbar.update(len(chunk))
        return save_name
    except Exception as e:
        print(traceback.format_exc())
        return False

def splitMolBlock_exist(file):
    """
    .sdfまたは.sdf.gzファイルからCIDとSMILESを取得するジェネレータ。
    """
    gzip_file = True
    if file.endswith(".sdf.gz"):
        fi = gzip.open(file, 'r')
        gzip_file = True
    elif file.endswith(".sdf"):
        fi = open(file, 'r')
        gzip_file = False
    else:
        raise Exception("file extension is not sdf or sdf.gz")
    count = 0
    i = 1
    try:
        end_flag = True
        cid_flag = False
        smiles_flag = False
        for line in fi:
            # print(line)
            if gzip_file:
                line = line.decode('utf-8')

            if cid_flag:
                cid = int(line)
                cid_flag = False
            if smiles_flag:
                smiles = line.replace('\n', '')
                smiles_flag = False

            if 'PUBCHEM_COMPOUND_CID' in line:
                cid_flag = True

            if 'PUBCHEM_OPENEYE_CAN_SMILES' in line:
                smiles_flag = True

            if line[:4] == "$$$$":
                # if not cid in db_ids:
                yield cid, smiles
                count += 1
                i += 1

        if end_flag:
            print(f"block {count} is done")
            yield cid, smiles
            print("done")
    except Exception as e:
        print(traceback.format_exc())
        pass
    fi.close()

化合物を部分構造に分解するための関数を作成します。

from rdkit import Chem
from rdkit.Contrib.mmpa import rfrag
from rdkit.Chem import SaltRemover

def fragment_mol(smiles, cid):
    try:
        remover = SaltRemover.SaltRemover()
        mol = Chem.MolFromSmiles(smiles)
        mol = remover.StripMol(mol)
        smi = Chem.MolToSmiles(mol)
        out = rfrag.fragment_mol(smi, cid)
        return out
    except:
        return []

def scaffoldFromSmiles(smiles, cid):
    for line in fragment_mol(smiles, cid):
        broken = line.rstrip().split(",")
        if broken[2]: # single cut
            continue
        smiles = broken[-1].split(".")
        mols = [Chem.MolFromSmiles(smi) for smi in smiles]
        numAtoms = [mol.GetNumAtoms() for mol in mols]
        if len(numAtoms) < 2:
            continue
        if numAtoms[0] > 5 and numAtoms[1] < 12:
            cid = broken[1]
            scaffold = smiles[0]
            rgroup = smiles[1]
            yield (cid, scaffold, rgroup)
        if numAtoms[1] > 5 and numAtoms[0] < 12:
            cid = broken[1]
            scaffold = smiles[1]
            rgroup = smiles[0]
            yield (cid, scaffold, rgroup)

分解した部分構造をcsvに書き込む関数作成します。

• 部分構造が共通する化合物を同じcsvに書き込む

import csv
import uuid
import os

def mmpaWriter(cid, smiles, scaffold, rgroup, directory='mmpa_group'):
    if not os.path.exists(directory):
        os.makedirs(directory)
    id = str(uuid.uuid3(uuid.NAMESPACE_DNS, scaffold))
    path = f'{directory}/{id}.csv'
    if not os.path.exists(path):
        with open(path, 'w') as f:
            writer = csv.writer(f, lineterminator='\n')
            writer.writerow(['scaffold', scaffold])

    with open(path, 'a') as f:
        writer = csv.writer(f, lineterminator='\n')
        writer.writerow([cid, smiles, rgroup])

上記の関数を使用してダウンロードしたSDFファイル内の化合物を部分構造でグループ化する

urls = pubchemDownloadList()
url = urls[0]
file_path = fetch(url[1], url[0])

count = 0
for cid, smiles in splitMolBlock_exist(file_path):
    for cid, scaffold, rgroup in scaffoldFromSmiles(smiles, cid):
        mmpaWriter(cid, smiles, scaffold, rgroup)

    count += 1
    if count > 1000:
        break

上記を実行すると部分構造の数だけmmpa_groupディレクトリにcsvファイルが生成される。

csvファイルには部分構造が共通した化合物が記載される。

生成されるcsvファイルの例（mmpa_group/2ac702c0-c5fd-3ad5-8261-17edbc7da08a.csv）

scaffold,CC(=O)C(O)[*:1]
179,CC(C(=O)C)O,C[*:1]
337,CC(=O)C(C(CO)O)O,OCC(O)[*:1]
641,CC(=O)C(C(COP(=O)(O)O)O)O,O=P(O)(O)OCC(O)[*:1]
669,CC(=O)C(COP(=O)(O)O)O,O=P(O)(O)OC[*:1]

1行目：共通している部分構造CC(=O)C(O)[*:1]

2行目以降：部分構造が共通した化合物。1列目がCID、2列目がSMILESが記載される。

株式会社piponでは技術でお困りのことがある方はオンライン相談が可能です。
こちらから会社概要資料をDLできます！
お問い合わせ内容に「オンライン相談希望」とご記載ください。
https://share.hsforms.com/19XNce4U5TZuPebGH_BB9Igegfgt

株式会社piponでは定期的に技術勉強会を開催しています。
ChatGPT・AI・データサイエンスについてご興味がある方は是非、ご参加ください。
https://chatgptllm.connpass.com/

株式会社piponではChatGPT・AI・データサイエンスについて業界ごとの事例を紹介しています。ご興味ある方はこちらのオウンドメディアをご覧ください。
https://bigdata-tools.com/