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

はじめに

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

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

[過去の検討例]

  1. RDKitの類似度マップを用いて原子ごとの寄与を可視化する
  2. NNの予測根拠可視化をライブラリ化する
  3. Visualize atom weight of AttentiveFP #DGL #RDKit #Chemoinformatics
  4. RDKitでフィンガープリントの可視化
  5. Using the new fingerprint bit rendering code
  6. Visualize important features of machine leaning #RDKit
  7. #6 RDKitのDrawMorganBitsはFCFPやunhashed FPでもできるのか確認した

 可視化の方向性は大きく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を用いた寄与度可視化 (比較用)

f:id:rkakamilan:20191214035624p:plain
logp_similaritymap

 今回の手法を適用した可視化です。ヒートマップは左軸に原子IDが振ってあるのですが、構造式中のどれに対応するか分かりづらいので原子IDを埋め込んだ図も合わせて出力するようにしました。

f:id:rkakamilan:20191214035812p:plain f:id:rkakamilan:20191214035832p:plain f:id:rkakamilan:20191214035849p:plain

 ...原子IDを埋め込んでもまだ見づらいですね。このあたりは改善の余地がありそうです。それでも構造式だけよりは寄与度の大小などが見やすくなったのではないでしょうか。

 logP以外にもTPSAやGasteiger Chargeでも同様に描画してみました。コードは以下にアップしています。

nbviewer.jupyter.org

 今回はレポジトリーの手法を原子ごとの重みが計算できる簡単な手法に適用してみました。次は機械学習モデルの結果解釈に本手法を使ってみようと思います。

(2019/12/14 23:00, 誤字修正)

(2019/12/17 07:30, コードをnbviewerに変更)