Activity cliffs are pairs of structurally similar compounds with large differences in activity, like inhibitory potency. I searched them in ChEMBL database with pychembldb.

I haven't maintained it for a long time but It works well even in Python3.7 :-)

My strategy for finding them is filtering assay data by whether it has confidence_score 9 (Direct single protein target assigned) and has enough data points(>=10 and <100) and then generating SMILES and Activity file from each of assay data for processing Matched Molecular Pairs Analysis (MMPs) by mmpdb

from pychembldb import *

for assay in chembldb.query(Assay).filter(Assay.confidence_score > 8).all():
    smiles_file = "{}.smi".format(assay.chembl_id)
    act_file = "{}.tsv".format(assay.chembl_id)
    act_list = [a for a in assay.activities if a.standard_value is not None and  and a.standard_units == "nM"]
    if len(act_list) > 9 and len(act_list) < 100::
        with open(smiles_file, "w") as sf:
            with open(act_file, "w") as af:
                af.write("ID\tVal\n")
                for activity in assay.activities:
                    if activity.standard_value is not None and \
                    activity.compound.molecule.structure is not None:
                        sf.write("{} {}\n".format(
                            activity.compound.molecule.structure.canonical_smiles,
                            activity.compound.molecule.chembl_id))
                        af.write("{}\t{}\n".format(
                            activity.compound.molecule.chembl_id,
                            activity.standard_value
                        ))

It took a few hours, but finally, I got around 55000 assay data.

I used mmpdb for generating MMP, and prepared a script for batch processing.

from glob import glob
import subprocess

def make_mmpdb(chembl_id):
    smiles_file = "{}.smi".format(chembl_id)
    act_file = "{}.tsv".format(chembl_id)
    fragments_file = "{}.fragments".format(chembl_id)
    sqlite_file = "{}.mmpdb".format(chembl_id)

    subprocess.run([
        "mmpdb",
        "fragment",
        smiles_file,
        "-o",
        fragments_file])

    subprocess.run([
        "mmpdb",
        "index",
        fragments_file,
        "-o",
        sqlite_file,
        "--properties",
        act_file])

if __name__ == "__main__":
    for fn in glob("*.smi"):
        chembl_id = fn.split(".")[0]
        make_mmpdb(chembl_id)

Finally, I extracted the pairs whose activity difference is more than 2 by pIC50.

import os
from sqlalchemy import *
from sqlalchemy.orm import create_session, relationship
from sqlalchemy.ext.declarative import declarative_base
from math import log10
from glob import glob

def search_ac(mmpdb_name):
    uri = 'sqlite:///{}'.format(mmpdb_name)

    Base = declarative_base()
    engine = create_engine(uri)
    metadata = MetaData(bind=engine)

    class Dataset(Base):
        __table__ = Table('dataset', metadata, autoload=True)

    class Compound(Base):
        __table__ = Table('compound', metadata, autoload=True)

    class PropertyName(Base):
        __table__ = Table('property_name', metadata, autoload=True)

    class CompoundProperty(Base):
        __table__ = Table('compound_property', metadata, autoload=True)
        compound = relationship('Compound', backref='compound_properties')
        property_name = relationship('PropertyName', backref='compound_properties')

    class RuleEnvironment(Base):
        __table__ = Table('rule_environment', metadata, autoload=True)

    class EnvironmentFingerprint(Base):
        __table__ = Table('environment_fingerprint', metadata, autoload=True)

    class ConstantSmiles(Base):
        __table__ = Table('constant_smiles', metadata, autoload=True)

    class Pair(Base):
        __table__ = Table('pair', metadata, autoload=True)
        constant_smiles = relationship('ConstantSmiles', backref='pair')
        rule_environment = relationship('RuleEnvironment', backref='pair')

    class RuleEnvironmentStatics(Base):
        __table__ = Table('rule_environment_statistics', metadata, autoload=True)

    sa = create_session(bind=engine)
    pIC50_dict = {}
    public_id = {}

    for cp in sa.query(CompoundProperty).filter_by(property_name_id=0).all():
        #print(cp.value)
        pIC50_dict[cp.compound_id] = 9.0 - log10(cp.value)
        public_id[cp.compound_id] = cp.compound.public_id

    pairs = []
    for pair in sa.query(Pair).all():
        if (pair.compound1_id, pair.compound2_id) not in pairs:
            pairs.append((pair.compound1_id, pair.compound2_id))
            diff = abs(pIC50_dict[pair.compound1_id] - pIC50_dict[pair.compound2_id])
            if diff >= 2.0:
                print("{}: ({}, {}) {}".format(
                    mmpdb_name.split(".")[0],
                    public_id[pair.compound1_id],
                    public_id[pair.compound2_id],
                    diff))

if __name__ == "__main__":
    for fn in glob("*.mmpdb"):
        search_ac(fn)

After this calculation I realized that I needed to eliminate cell based assays that caused AC noises.

医薬品研究開発における Modeling and Simulation（M & S）手法の紹介という大変に面白い総説を教えてもらって読んでいる。

これは、事前にTrkAがpainのターゲットであることを知っているからこういう結論に出来たのか、そういう事前知識なしに自然に導かれたのか興味あるところではある。モデルからターゲット候補が出てくるなんて素敵。

QSP モデルによって，今までのモデル解析ではできなかった予測を行った事例を以下に 2 例紹介する．1 つ目は，MID3 白書に取り上げられている神経成長因子（NGF）パスウェイに対する QSP 解析 により，疼痛治療の新規の創薬ターゲットを検討した事例（Example 3）である．この事例では，すで に報告されていた NGF パスウェイを ODE で記述される数理モデル（Model 1）に変換し，さらにその パスウェイモデルを関連因子の神経細胞内外での挙動を反映させたモデルに連結して生理学的に拡張し たモデル（Model 2）を作成して NGF パスウェイの下流での創薬ターゲット候補の探索を行っている． 効果測定の指標として NGF 刺激の結果として観察され，疼痛発現につながると考えられる Diphosphorylated extracellular signal-regulated kinase（dppERK）の蓄積量を用い，モデルの妥当性検討は すでに報告されている NGF 阻害抗体で得られている結果と比較することにより行った．こうして構築 されたモデル内に含まれる因子を，in silico で仮想的に様々な阻害強度で阻害する感度分析により， Tropomyosin receptor kinase A（TrkA）が薬物治療のターゲット候補となる可能性が予測された．

尚、MoAと書くと現状はBabymetalがおすすめされるようです。

「LEGEND - S - BAPTISM XX - 」 (LIVE AT HIROSHIMA GREEN ARENA) [Blu-ray]

トイズファクトリー / 6471円 ( 2018-08-01 )

• (MID3) Good Practice - European Medicines Agency
• Good Practices in Model-Informed Drug Discovery and Development: Practice, Application, and Documentation.
• An introduction to PK-Sim®: the open source platform for PBPK modeling

PK-SimってOSSになったの？

最近、バイオインフォマティクスのチームに兼務になりました。このバイオインフォマティクスチームに求められてることが、いわゆる普通のターゲットファインディング的なバイオインフォだけじゃなくてPKPD(今だとQSPって言うの？)でのインフォマティクス側からの貢献を求められたりするので10連休は色々勉強するよいチャンスだったりします。しかし、PKモデリングは基本的なところはわかるけど、PDは全然だから着いていくのつらいっす。

それから、組織内でモダリティのブームが起きているみたいなので、そっちのほうも色々キャッチアップしなきゃいけない。

で、自分の解釈でいうところのモダリティというのは、「従来のkey&lock型の低分子創薬を超えて、MOAをうまく利用して生体に干渉するような創薬をやりましょう」ってことだと思うんですよね。つまりライフハック。

そういう意味ではこれからの分子設計は低分子だけじゃなく、中分子や核酸やdegradationでも貢献する必要があるだろうし、創薬ターゲットのMOAに応じて柔軟に対応出来るような組織になっているべきでしょうね。と考えるとMOAを理解できることが必須になるんじゃないですかね。

というわけで、分子生物学も少し復習しないとなぁという気分になったのでそれっぽい本を読んでいました。遺伝子ドライブ面白いですね。

ゲノム編集の基本原理と応用: ZFN,TALEN,CRISPR-Cas9
山本 卓
裳華房 / 2808円 ( 2018-06-06 )

コラムにCRISPRdirectが紹介されていて、「あー作者の人知っとるわ」ってなった。

創薬プロジェクト用のIssue Tracking System(ITS)を探しているのだけど良いものがみつからない。

で、ないのなら作ってしまえということでとりあえずE-R図でも書こうかとnode.jsインストールしてwwwsqldesignerでやってみた

テーブルスキーマ出来たんであとはFlaskでサクサクつくればよろしい。

うちのSlackオルタネティブのAI創薬絡みのチャンネルではArxivの機械学習新着が流れてくるようになっていて、たまに大きな桃がどんぶらこどんぶらこと流れてきては、中身をぶった切ってチャンネル参加者の皆でキャッキャウフフしてるわけです。

で、今回流れてきたのがStep Change Improvement in ADMET Prediction with PotentialNet Deep FeaturizationというPande先生のとこで開発したPotentialNetっていう手法をメルクのデータで検証したら調子良かった（良すぎ！）っていう論文です。

なんかfuとかmicrosomeのクリアランスとかADMEの予測がすごい改善していて期待度高いです。PotentialNetってのはGCNN(Neural FP)のとこにLeRUじゃなくてGRUつかうことで特徴量抽出を工夫しているようです。具体的になにがよろしくなってどういう特徴量が抽出されているのかイメージつかめないんですが、何が良くなるんでしょうね？

それから実装はPyTorchらしいけどGithubには上がってないようなので自分で再実装する必要ありそうだし、最初にAtomにどういう情報持たせるのかに工夫の余地はありそう。

PotentialNetはいずれdeepchemに取り入れられるのでしょうか？

• PotentialNet for Molecular Property Prediction
• PotentialNet for Molecular Property Prediction(ACS)
• GATED GRAPH SEQUENCE NEURAL NETWORKS
• Deep learning for molecules, introduction to chainer chemistry
• A PyTorch Implementation of GGNN

深層学習も積みっぱなしだから読まないと、、、

深層学習 (アスキードワンゴ)

ドワンゴ / ?円 ( 2018-08-27 )

参加者の皆様お疲れ様でした。楽しんでいただけて且つ何か得るものがあったのなら企画の甲斐があったなと思います。

それにしてもauto-DMTAやばかったですね。我々も頑張らないといけないなーと思いました。

さて、Thierryさんとはどういう経由で呼んだの？とかなり聞かれたのですが、実は去年のICCSでポスター出してたのでネホリンハホリンしたら次の日のexcursionで話しかけてきたので、「来年AI in Medicinal Chemistry企画してるんだけど興味ある？」って聞いたら「めっちゃある、呼ばれたら絶対行くで！」って、まぁそんな感じで決まりました。ICCSは面白い割に日本からの参加者はほぼいないので機会があれば参加すると良いと思います（3年ごとに開催）。

尚、ランチははん亭にしました。

演者は懇親会は質問攻めにあってあまり食べられないので、皆で帰りに「にし乃」に寄りました。

美味しかった。

ちょっと質問されたので、py4chemoinformaticsの説明追加してみたんだけど自分の中で問題が整理されて良かった。

結局acivity cliffとかmagic methylってのはApplicability domainの問題ではなくて、特徴量の抽出のほうの問題なんですよね。グラフ類似性が実際の三次元の構造の類似性とは微妙に違うし、結合モードの類似性とかとも異なるのに、その類似性が活性の類似性と相関するという仮説そのものが正しくない可能性があるということを認識した上でモデルを作らんといかんよなと再認識させられた。

尚、CoMFAはぶつかって活性消失したという事実がないとモデリングできないし、ファーマコフォアの排除体積も同様。ドッキングのスコア関数はそれっぽい特徴表現だけど、粗すぎてネガティブスクリーニングにしか使えんけど、ファーマコフォア表現ってのが特徴量としては一番しっくりくるかなと言ったところ。

ファーマコフォアはもともとは、同じターゲットで複数の異なる骨格の薬剤とか開発化合物を3次元的に重ね合わせて共通の特徴を抽出するLBDD的な手法により推定されていてCatalystを使って計算するのが一般的だったように思います。

で、ファーマコフォアを作ろうとする話で触れられているように、ファーマコフォアをSBDDの側から解釈するとターゲット蛋白質でのポケット内で、強く相互作用している残基とリガンドの官能基の相互作用を見ていることになるわけです。

というわけで、ターゲットのアポ体の構造がわかっていたり、新規なスキャフォールドが欲しかったり（知財の関係で）する場合にターゲットのポケットの中で水素結合をする可能性のある残基や主査に水素結合ポイントを定義して三点ファーマコフォアでスクリーニングするというようなこともよくやっていましたが、これは結構ヒット率低かったです。

そもそも、ファーマコフォアポイントじゃないところにファーマコフォアを仮定する時点で外れ確定なのでそのあたりが難しい要因でした。

• リガンド重ね合わせのアプローチの場合、似た構造を重ね合わせに使ってしまい、アーティファクトが出やすい（ベンゼン環起因のアロマティックなポイントとか）、そもそも異なるポケットに結合しているとか別の残基と相互作用している化合物が混じっているとかで共通項が取れない場合もある。
• 結晶構造から推定する場合は、水素結合はいいとして、アロマティックなフィーチャーは設定しにくい。Catalystが不安定なコンフォマーもファーマコアマッピングするので、それヒットにしていいの？と疑問を呈したい場合も多々あった。

なので、もし複合体結晶構造がある場合は一度FMOをかけて、きちんと相互作用を確認してからファーマコフォアぽいイントを設定するなり、ドッキングシミュレーションをしたほうが成功確率が上がります。

こんな感じのFMO計算用のAMIがあるとクラスターもってなくてもサクッと計算できるから良いかもしれませんね。てか、今度作ろっと。

PsikitでEnergy Decomposition Analysis(EDA)をやりたくてKitaura-Morokuma Analysisとか実装されてないかなーと調べてたんですが、SAPTってのを使えばよろしいらしいところまでは到達した。

でもって、そのあたり使ってる論文がJ.C.I.Mあたりに投稿されてないかなーとさらに調べた結果、Tuning Stacking Interactions between Asp–Arg Salt Bridges and Heterocyclic Drug Fragmentsという論文をみつけて読んだらやばかったという話。

Asp–Arg塩橋とヘテロ環のスタッキングってよくみられるし強い相互作用っぽいんだけど、よくわかってないよね？だから量子化学計算で調べるわっていう内容で、単環、二環、三環のヘテロ環（窒素シャッフリングが多い）の63モデルでローカルミニマム、グローバルミニマムを探索して、どういう位置関係になっているかとかEDAしてみてどういう成分が効いているのかを調べていた。

得られた結論としては塩橋だけどヘテロ環はカチオニックな残基(LYS,ARG)との相互作用のほうが安定に形成しそうだってことと、ElectroStaticな項が相互作用に支配的で、DispersionInteractionがまぁまぁ効いてそうってこと。それから双極子モーメントはほとんど寄与していないってことから無指向性の相互作用ってこと。

論文では塩橋をモデル化しているけど、内容をよく考えるてみるとこれはおそらくカチオニックな残基単体でも成り立つだろうし、そっちのほうが強いんだろうなぁというところまで理解できたので非常に良かった。うまくやればプロアクティブにヘテロ環導入して活性向上できるかなと思う。一方でなんでカチオニックな残基との相互作用が強いの？ってことに関してはおそらく窒素シャッフリングの環だからだろうなぁと思った。

じゃぁアニオニックな残基(ASP,GLU)と強く相互作用するようなヘテロ環、アリルは何よ？っていう疑問が湧くと思うんだけど、そっちに関しては既に理解しているのでそのあたりまとめてどこかで発表しようかなと思っている。ちょっとコントラバーシャルかなーと思わないこともないけど、量子化学系の人たちがもっとSBDDに参入するきっかけになってくれると嬉しいなーと。

ドッキングは反応物と生成物が変化しない化学反応であるという主張を広めたい

— kzfm (@fmkz___) 2019年3月1日

Considering the conformational effects of the compound is important in Structure Based Drug Design, this paper discussed about it, in terms Protein-Ligand binding using torsional scan of each ligands(PDB:2JH0,2JH5,2JH6). They calculated torsional energy, and explained the relation between inhibitory activity and torsional energy.

Torsional scanning is the task of the quantum chemistry rather than that of chemoinformatics. But I wanted to conduct quantum chemical calculation as an extention of chemoinformatics way, so I implemented it in Psikit

• Torsion scan example

日本語訳

jupyterでRDKitからのB3LYP/6-31G*でのDFT計算さっくりうごくの素晴らしい。