はじめに

BRICSによって生成したフラグメントを組み合わせることで新しい分子構造を生成する。(BRICSについて)

前回は、二つのフラグメント(A,B)を結合することで、A-B型の分子を生成したが、今回は、三つのフラグメント（A,B,C）を結合させることで A-B-C型の分子を生成する。

水溶解度データ（金子研究室HPのデータセット）に含まれる1290分子の構造をフラグメント化し構造生成に利用する。

今回のコードは明治大学金子先生の"structure_generator_based_on_r_group"を参考に（コードを改変して）作成しました。(金子先生ありがとうございます)

前準備

モジュールのインストール

conda install -c conda-forge rdkit

コード

モジュールのインポート

import numpy as np
import itertools

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, BRICS, Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors

フラグメントの準備

# データの準備
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols_list = [mol for mol in suppl if mol is not None]

# １290分子を分解してフラグメントを準備する
fragment_set = set()
for mol in mols_list:
    fragment = BRICS.BRICSDecompose(mol)
    fragment_set = fragment_set | fragment

fragment_smiles = list(fragment_set)


# フラグメントのリストから、ダミーアトムの数が一つ、二つのものを抽出する。
frag_1dummy_s = [smiles for smiles in fragment_smiles if smiles.count('*')==1]
frag_2dummy_s = [smiles for smiles in fragment_smiles if smiles.count('*')==2]
frag_1dummy = np.array([Chem.MolFromSmiles(smiles) for smiles in frag_1dummy_s])
frag_2dummy = np.array([Chem.MolFromSmiles(smiles) for smiles in frag_2dummy_s])


# 複雑なフラグメント（分子量300以上）は削除する
descriptor_calc = MoleculeDescriptors.MolecularDescriptorCalculator(['MolWt'])
MolWt_list1 = np.array([descriptor_calc.CalcDescriptors(mol)[0] for mol in frag_1dummy])
frag_1dummy = frag_1dummy[np.where(MolWt_list1<=300)]

descriptor_calc = MoleculeDescriptors.MolecularDescriptorCalculator(['MolWt'])
MolWt_list2 = np.array([descriptor_calc.CalcDescriptors(mol)[0] for mol in frag_2dummy])
frag_2dummy = frag_2dummy[np.where(MolWt_list2<=300)]

print('number of fragment:',len(frag_1dummy))
# >>> number of fragment: 333
print('number of fragment:',len(frag_2dummy))
# >>> number of fragment: 130

ダミーアトム(*)が一つのフラグメント

img = Draw.MolsToGridImage(frag_1dummy[:3], molsPerRow=3,legends=frag_1dummy_s[:3])
img.save('frag_1dummy.png')
img

ダミーアトム(*)が二つのフラグメント

img = Draw.MolsToGridImage(frag_2dummy[:3], molsPerRow=3,legends=frag_2dummy_s[:3])
img.save('frag_2dummy.png')
img

構造生成する関数を作成する

フラグメントの一つを主骨格(main_mol,A)とみなしてそこにフラグメント(fragment1,2:B,C)を結合するテイ。

def structure_generator(main_mol, fragment1, fragment2, r_position=1):
    
    """
    
    parameters
    ----------
    main_mol : RDkit Mol object
        主骨格（A）の構造（フラグメント）
    fragment1 : RDkit Mol object
        * の数が二つのフラグメント（Bにあたる）
    fragment2 : RDkit Mol object
        * の数が一つのフラグメント(Bにあたる）
    r_position : int (1 or 2) 
        fragment1 中の * のどちらに、main_mol, fragment1を結合させるか。
        
    Returns
    ----------
    generated_molecule : RDkit Mol object
        各フラグメントを結合させて生成したA-B-C型の分子
    
    
    """
    bond_list = [Chem.rdchem.BondType.UNSPECIFIED, Chem.rdchem.BondType.SINGLE, Chem.rdchem.BondType.DOUBLE,
             Chem.rdchem.BondType.TRIPLE, Chem.rdchem.BondType.QUADRUPLE, Chem.rdchem.BondType.QUINTUPLE,
             Chem.rdchem.BondType.HEXTUPLE, Chem.rdchem.BondType.ONEANDAHALF, Chem.rdchem.BondType.TWOANDAHALF,
             Chem.rdchem.BondType.THREEANDAHALF, Chem.rdchem.BondType.FOURANDAHALF, Chem.rdchem.BondType.FIVEANDAHALF,
             Chem.rdchem.BondType.AROMATIC, Chem.rdchem.BondType.IONIC, Chem.rdchem.BondType.HYDROGEN,
             Chem.rdchem.BondType.THREECENTER, Chem.rdchem.BondType.DATIVEONE, Chem.rdchem.BondType.DATIVE,
             Chem.rdchem.BondType.DATIVEL, Chem.rdchem.BondType.DATIVER, Chem.rdchem.BondType.OTHER,
             Chem.rdchem.BondType.ZERO]
    
    # make adjacency matrix and get atoms for main molecule
    level1_adjacency_matrix = Chem.rdmolops.GetAdjacencyMatrix(fragment1)

    for bond in fragment1.GetBonds():
        level1_adjacency_matrix[bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()] = bond_list.index(bond.GetBondType())
        level1_adjacency_matrix[bond.GetEndAtomIdx(), bond.GetBeginAtomIdx()] = bond_list.index(bond.GetBondType())
    level1_atoms = []
    for atom in fragment1.GetAtoms():
        level1_atoms.append(atom.GetSymbol())

    # * の位置を検索
    r_index_in_level1_molecule_old = [index for index, atom in enumerate(level1_atoms) if atom == '*']

    # matrix,atom の並び替え
    for index, r_index in enumerate(r_index_in_level1_molecule_old):
        modified_index = r_index - index
        atom = level1_atoms.pop(modified_index)
        level1_atoms.append(atom)
        tmp = level1_adjacency_matrix[:, modified_index:modified_index + 1].copy()
        level1_adjacency_matrix = np.delete(level1_adjacency_matrix, modified_index, 1)
        level1_adjacency_matrix = np.c_[level1_adjacency_matrix, tmp]
        tmp = level1_adjacency_matrix[modified_index:modified_index + 1, :].copy()
        level1_adjacency_matrix = np.delete(level1_adjacency_matrix, modified_index, 0)
        level1_adjacency_matrix = np.r_[level1_adjacency_matrix, tmp]

    r_index_in_level1_molecule_new = [index for index, atom in enumerate(level1_atoms) if atom == '*']

    # *がどの原子と結合しているか調べる。
    r_bonded_atom_index_in_level1_molecule = []
    for number in r_index_in_level1_molecule_new:
        r_bonded_atom_index_in_level1_molecule.append(np.where(level1_adjacency_matrix[number, :] != 0)[0][0])

    #結合のタイプを調べる
    r_bond_number_in_level1_molecule = level1_adjacency_matrix[
        r_index_in_level1_molecule_new, r_bonded_atom_index_in_level1_molecule]

    # *原子を削除
    level1_adjacency_matrix = np.delete(level1_adjacency_matrix, r_index_in_level1_molecule_new, 0)
    level1_adjacency_matrix = np.delete(level1_adjacency_matrix, r_index_in_level1_molecule_new, 1)

    for i in range(len(r_index_in_level1_molecule_new)):
        level1_atoms.remove('*')
    level1_size = level1_adjacency_matrix.shape[0]

    # level1フラグメントの ２つの* に残りのフラグメントを結合する
    generated_molecule_atoms = level1_atoms[:]
    generated_adjacency_matrix = level1_adjacency_matrix.copy()

    frag_permutations = list(itertools.permutations([main_mol, fragment2]))
    fragment_list = frag_permutations[r_position]

    for r_number_in_molecule, fragment_molecule in enumerate(fragment_list):  

        #ここから、主骨格の処理と同じ
        fragment_adjacency_matrix = Chem.rdmolops.GetAdjacencyMatrix(fragment_molecule)

        for bond in fragment_molecule.GetBonds():
            fragment_adjacency_matrix[bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()] = bond_list.index(
                bond.GetBondType())
            fragment_adjacency_matrix[bond.GetEndAtomIdx(), bond.GetBeginAtomIdx()] = bond_list.index(
                bond.GetBondType())
        fragment_atoms = []
        for atom in fragment_molecule.GetAtoms():
            fragment_atoms.append(atom.GetSymbol())

        # integrate adjacency matrix
        r_index_in_fragment_molecule = fragment_atoms.index('*')

        r_bonded_atom_index_in_fragment_molecule = \
            np.where(fragment_adjacency_matrix[r_index_in_fragment_molecule, :] != 0)[0][0]

        # 後で * を削除するのでそのための調整（たぶん）
        if r_bonded_atom_index_in_fragment_molecule > r_index_in_fragment_molecule:
            r_bonded_atom_index_in_fragment_molecule -= 1

        fragment_atoms.remove('*')
        fragment_adjacency_matrix = np.delete(fragment_adjacency_matrix, r_index_in_fragment_molecule, 0)
        fragment_adjacency_matrix = np.delete(fragment_adjacency_matrix, r_index_in_fragment_molecule, 1)

        #新たに生成する分子用のマトリックス作成
        main_size = generated_adjacency_matrix.shape[0]
        generated_adjacency_matrix = np.c_[generated_adjacency_matrix, np.zeros(
            [generated_adjacency_matrix.shape[0], fragment_adjacency_matrix.shape[0]], dtype='int32')]
        generated_adjacency_matrix = np.r_[generated_adjacency_matrix, np.zeros(
            [fragment_adjacency_matrix.shape[0], generated_adjacency_matrix.shape[1]], dtype='int32')]

        #マトリックスに結合のタイプを記述
        generated_adjacency_matrix[r_bonded_atom_index_in_level1_molecule[r_number_in_molecule], 
                                   r_bonded_atom_index_in_fragment_molecule + main_size] = \
            r_bond_number_in_level1_molecule[r_number_in_molecule]

        generated_adjacency_matrix[r_bonded_atom_index_in_fragment_molecule + main_size, 
                                   r_bonded_atom_index_in_level1_molecule[r_number_in_molecule]] = \
            r_bond_number_in_level1_molecule[r_number_in_molecule]

        generated_adjacency_matrix[main_size:, main_size:] = fragment_adjacency_matrix
        #フラグメントのマトリックスを入力(マトリックスの右下)

        # integrate atoms
        generated_molecule_atoms += fragment_atoms

    # generate structures 
    generated_molecule = Chem.RWMol()
    atom_index = []

    for atom_number in range(len(generated_molecule_atoms)):
        atom = Chem.Atom(generated_molecule_atoms[atom_number])
        molecular_index = generated_molecule.AddAtom(atom)
        atom_index.append(molecular_index)

    for index_x, row_vector in enumerate(generated_adjacency_matrix):    
        for index_y, bond in enumerate(row_vector):      
            if index_y <= index_x:
                continue
            if bond == 0:
                continue
            else:
                generated_molecule.AddBond(atom_index[index_x], atom_index[index_y], bond_list[bond])

    generated_molecule = generated_molecule.GetMol()

    return generated_molecule

テスト

main_mol = frag_1dummy[0]
fragment1 = frag_2dummy[1]
fragment2 = frag_1dummy[1]

mol = structure_generator(main_mol, fragment1, fragment2)

img = Draw.MolsToGridImage([main_mol, fragment1, fragment2, mol],
                           molsPerRow=4,
                           legends=['main_mol','fragment1','fragment2','generated_molecule'])
img.save('generated_molecule.png')
img

上手く分子を作れた。

参考

hkaneko1985/structure_generator_based_on_r_group