はじめに

BRICSによって生成したフラグメントを組み合わせることで新しい分子構造を生成する。(BRICSについて) 以前は、BRICSBuildによってランダムに構造を生成したが、今回は、二つのフラグメント（A, B）を結合させることで新たな分子（A-B）を生成する。

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

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

前準備

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

conda install -c conda-forge rdkit

コード

モジュールのインポート

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

フラグメントの準備

二つのフラグメントをつなぐので、ダミーアトム（結合可能部位）を一つ持つ分子をピックアップする。
複雑な構造のフラグメント（分子量300以上）は削除する。

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_1dummy = np.array([Chem.MolFromSmiles(smiles) for smiles in frag_1dummy_s])


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

frag_1dummy = frag_1dummy[np.where(MolWt_list<=300)]

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

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

フラグメント二つから分子一つを作る

#フラグメントを選ぶ（フラグメントの一つを主骨格(main_mol)とみなしてそこにフラグメント(fragment)を結合するテイ。
main_mol = frag_1dummy[0]
fragment = frag_1dummy[1]

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]

# main_mol の操作 =======================================================
# main_mol のadjacency matrixの作成と原子の情報を取得する
main_adjacency_matrix = Chem.rdmolops.GetAdjacencyMatrix(main_mol)

for bond in main_mol.GetBonds():
    main_adjacency_matrix[bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()] = bond_list.index(bond.GetBondType())
    main_adjacency_matrix[bond.GetEndAtomIdx(), bond.GetBeginAtomIdx()] = bond_list.index(bond.GetBondType())
    
main_atoms = []
for atom in main_mol.GetAtoms():
    main_atoms.append(atom.GetSymbol())

# ダミーアトム（*）の位置を調べる
r_index_in_main_molecule_old = [index for index, atom in enumerate(main_atoms) if atom == '*']

# matrix,atom を並び替えて新しく分子を生成する。
for index, r_index in enumerate(r_index_in_main_molecule_old):
    modified_index = r_index - index
    atom = main_atoms.pop(modified_index)
    main_atoms.append(atom)
    tmp = main_adjacency_matrix[:, modified_index:modified_index + 1].copy()
    main_adjacency_matrix = np.delete(main_adjacency_matrix, modified_index, 1)
    main_adjacency_matrix = np.c_[main_adjacency_matrix, tmp]
    tmp = main_adjacency_matrix[modified_index:modified_index + 1, :].copy()
    main_adjacency_matrix = np.delete(main_adjacency_matrix, modified_index, 0)
    main_adjacency_matrix = np.r_[main_adjacency_matrix, tmp]

#新しく生成した分子でのダミーアトム(*)の位置を取得する
r_index_in_main_molecule_new = [index for index, atom in enumerate(main_atoms) if atom == '*']

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

#結合のタイプを調べる
r_bond_number_in_main_molecule = main_adjacency_matrix[
    r_index_in_main_molecule_new, r_bonded_atom_index_in_main_molecule]

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

for i in range(len(r_index_in_main_molecule_new)):
    main_atoms.remove('*')
main_size = main_adjacency_matrix.shape[0]
# Frag_1 の操作ここまで ========================================================

#Frag_2の操作ここまで===========================================================
r_number_in_molecule = 0

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

for bond in fragment.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.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)
#Frag_2の操作ここまで===========================================================

#Frag_1とFrag_2をつなげる=========================================
generated_molecule_atoms = main_atoms[:]
generated_adjacency_matrix = main_adjacency_matrix.copy()

#新たに生成する分子用のマトリックス作成
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_main_molecule[r_number_in_molecule], 
                           r_bonded_atom_index_in_fragment_molecule + main_size] = \
    r_bond_number_in_main_molecule[r_number_in_molecule]

generated_adjacency_matrix[r_bonded_atom_index_in_fragment_molecule + main_size, 
                           r_bonded_atom_index_in_main_molecule[r_number_in_molecule]] = \
    r_bond_number_in_main_molecule[r_number_in_molecule]

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

# 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()
smiles = Chem.MolToSmiles(generated_molecule)

generated_mol = Chem.MolFromSmiles(smiles)


#作った分子を描画する
img = Draw.MolsToGridImage([main_mol,fragment,generated_mol],
                           molsPerRow=3,
                           legends=['main_mol','fragment','generated_mol'])
img.save('generated_mol.png')
img

関数化してみた

def structure_generator(main_mol,fragment):

    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]

    # main_mol の操作 =======================================================
    # main_mol のadjacency matrixの作成と原子の情報を取得する
    main_adjacency_matrix = Chem.rdmolops.GetAdjacencyMatrix(main_mol)

    for bond in main_mol.GetBonds():
        main_adjacency_matrix[bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()] = bond_list.index(bond.GetBondType())
        main_adjacency_matrix[bond.GetEndAtomIdx(), bond.GetBeginAtomIdx()] = bond_list.index(bond.GetBondType())

    main_atoms = []
    for atom in main_mol.GetAtoms():
        main_atoms.append(atom.GetSymbol())

    # ダミーアトム（*）の位置を調べる
    r_index_in_main_molecule_old = [index for index, atom in enumerate(main_atoms) if atom == '*']

    # matrix,atom を並び替えて新しく分子を生成する。
    for index, r_index in enumerate(r_index_in_main_molecule_old):
        modified_index = r_index - index
        atom = main_atoms.pop(modified_index)
        main_atoms.append(atom)
        tmp = main_adjacency_matrix[:, modified_index:modified_index + 1].copy()
        main_adjacency_matrix = np.delete(main_adjacency_matrix, modified_index, 1)
        main_adjacency_matrix = np.c_[main_adjacency_matrix, tmp]
        tmp = main_adjacency_matrix[modified_index:modified_index + 1, :].copy()
        main_adjacency_matrix = np.delete(main_adjacency_matrix, modified_index, 0)
        main_adjacency_matrix = np.r_[main_adjacency_matrix, tmp]

    #新しく生成した分子でのダミーアトム(*)の位置を取得する
    r_index_in_main_molecule_new = [index for index, atom in enumerate(main_atoms) if atom == '*']

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

    #結合のタイプを調べる
    r_bond_number_in_main_molecule = main_adjacency_matrix[
        r_index_in_main_molecule_new, r_bonded_atom_index_in_main_molecule]

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

    for i in range(len(r_index_in_main_molecule_new)):
        main_atoms.remove('*')
    main_size = main_adjacency_matrix.shape[0]
    # Frag_1 の操作ここまで =======================================================

    #Frag_2の操作ここまで===========================================================
    r_number_in_molecule = 0

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

    for bond in fragment.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.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)
    #Frag_2の操作ここまで===========================================================

    #Frag_1とFrag_2をつなげる=========================================
    generated_molecule_atoms = main_atoms[:]
    generated_adjacency_matrix = main_adjacency_matrix.copy()

    #新たに生成する分子用のマトリックス作成
    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_main_molecule[r_number_in_molecule], 
                               r_bonded_atom_index_in_fragment_molecule + main_size] = \
        r_bond_number_in_main_molecule[r_number_in_molecule]

    generated_adjacency_matrix[r_bonded_atom_index_in_fragment_molecule + main_size, 
                               r_bonded_atom_index_in_main_molecule[r_number_in_molecule]] = \
        r_bond_number_in_main_molecule[r_number_in_molecule]

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

    # 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()
    smiles = Chem.MolToSmiles(generated_molecule)

    generated_mol = Chem.MolFromSmiles(smiles)

    return generated_mol

関数を使い、ランダムに構造を生成してみる

frag_1 = frag_1dummy[np.random.randint(len(frag_1dummy))]
frag_2 = frag_1dummy[np.random.randint(len(frag_1dummy))]

generated_mol = structure_generator(frag_1,frag_2)
#作った分子を描画する
img = Draw.MolsToGridImage([frag_1,frag_2,generated_mol],
                           molsPerRow=3,
                           legends=['main_mol','fragment','generated_mol'])
img.save('generated_mol_random.png')
img