少し拡張した化合物の可視化を試してみた

はじめに

　創薬 (dry) Advent Calendar 2019 第13日目の記事です。

　最近では機械学習の結果の解釈性に注目が集まっています。その中でケムインフォマティクスの世界でもモデルの予測解釈性を上げるために、判断の根拠となる部分構造を可視化する手法が検討されています。

[過去の検討例]

　可視化の方向性は大きく2つあります。(A) 1つの化合物に対して重要度の高い原子/判断の根拠となった原子をハイライトする (記事1-3); (B) 重要なfingerprintのbitを抜き出して描画する (記事4-7)。各化合物ごとにどの部分構造が寄与しているかみたい時には(A)の手法が良いのですが、ハイライトした構造式だけではやや定性的過ぎて個人的にはもう少し情報が欲しいなぁと思ってました。

f:id:rkakamilan:20191213165052j:plain — 記事2より転載

　そんな中、構造式だけではなく、原子ごとの重みを別途可視化した図を使っている論文/レポが最近出ていたのでその手法を使って簡単な可視化を試してみました。

[1905.13686] Explainability Techniques for Graph Convolutional Networks

github.com

　元論文はGCNで予測したモデルの解釈性に関する手法を提案したものですが、詳細は本記事の主題ではないので元論文をご参照頂ければと思います。

　個人的に気に入ったのは原子/結合ごとの重みを可視化したヒートマップ+棒グラフで、構造式中で原子をハイライトしただけでは分かりづらい色合いの微妙な違いも分かりやすくなっています。

f:id:rkakamilan:20191213170540p:plain — 原子ごとの重みを構造式に埋め込んだもの

f:id:rkakamilan:20191213170543p:plain — 原子ごとのヒートマップ/棒グラフ

f:id:rkakamilan:20191213170555p:plain — 結合ごとのヒートマップ/棒グラフ

脂溶性の寄与度可視化への適用

　RDkitのドキュメントにも載っているlogPを事例として試しました。元論文とは異なり結合の重みは算出されてないので原子の重みだけを使って少し体裁を変えています。

　まず、各原子の重みはドキュメント同様に以下で計算できます。

def calc_crippen_weight(mol):
    contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
    weights = [x for x,y in contribs]
    return weights

　ヒートマップに用いる重みのデータは、以下のようにデータフレームとして準備しました。(numpy.arrayのままでも良かった)

symbols = [f'{mol.GetAtomWithIdx(i).GetSymbol()}_{i}' for i in range(mol.GetNumAtoms())] 
df = pd.DataFrame(columns=atoms)
contribs = weight_fn(mol) # calc_crippen_weight
num_atoms = mol.GetNumAtoms()
arr = np.zeros((num_atoms, len(atoms)))
for i in range(mol.GetNumAtoms()):
    _a = mol.GetAtomWithIdx(i).GetSymbol()
    arr[i,atoms.index(_a)] = contribs[i]
df = pd.DataFrame(arr, index=symbols, columns=atoms)

　構造式は Chem.Draw.rdMolDraw2D.MolDraw2DCairo で描画します。描画の際の色は原子ごとの重みから計算します。今回は chainer-chemistryで使われているものを拝借しました。

# from chainer-chemistryの visualizer_utils.pyより
def red_blue_cmap(x):
    """Red to Blue color map
    Args:
        x (float): value between -1 ~ 1, represents normalized saliency score
    Returns (tuple): tuple of 3 float values representing R, G, B.
    """
    if x > 0:
        # Red for positive value
        # x=0 -> 1, 1, 1 (white)
        # x=1 -> 1, 0, 0 (red)
        return 1.0, 1.0 - x, 1.0 - x
    else:
        # Blue for negative value
        x *= -1
        return 1.0 - x, 1.0 - x, 1.0

def is_visible(begin, end):
    if begin <= 0 or end <= 0:
        return 0
    elif begin >= 1 or end >= 1:
        return 1
    else:
        return (begin + end) * 0.5
    
def color_bond(bond, saliency, color_fn):
    begin = saliency[bond.GetBeginAtomIdx()]
    end = saliency[bond.GetEndAtomIdx()]
    return color_fn(is_visible(begin, end))

# 重みを絶対値最大の値でスケーリング
weights, vmax = SimilarityMaps.GetStandardizedWeights(contribs)
vmin = -vmax
# ハイライトする色を算出
atom_colors = {i: red_blue_cmap(e) for i, e in enumerate(weights)}
bondlist = [bond.GetIdx() for bond in mol.GetBonds()]
bond_colors = {i: color_bond(bond, weights, red_blue_cmap) for i, bond in enumerate(mol.GetBonds())}

　算出したものを用い、描画します。

drawer = rdMolDraw2D.MolDraw2DCairo(molSize[0],molSize[1])
drawer.drawOptions().useBWAtomPalette()
drawer.drawOptions().padding = .2
drawer.DrawMolecule(
    mc,
    highlightAtoms=[i for i in range(len(atom_colors))], # 上で計算した原子の色
    highlightAtomColors=atom_colors, 
    highlightBonds=[i for i in range(len(bond_colors))], 
    highlightBondColors=bond_colors,
    highlightAtomRadii={i: .5 for i in range(len(atom_colors))}
)
drawer.FinishDrawing()

　これで材料は揃ったので後はレポジトリーにあるコードを適宜変更して描画していくだけです。

　RDKitのドキュメントにあるSimilarityMapsを用いた寄与度可視化 (比較用)