Greatな岡山の学会に行ってGreatな焼き鳥を味わった

great岡山の学会に参加してきた。

夜は大井町の鳥たかというお店に通った。うますぎ天国。

1568212594 1568212601

晩杯屋はうるさくなければいいけど、酔っぱらい学生とか来たら出る感じ。

1568212599

お昼はラーメン。

1568212597

CWLに入信したの図

1568212604

たまごかけ大学行くのはYAPC以来やなぁとか言っていたけどそのあとIBISで参加していたことを後に思い出した。

ま、どうでもいい。

ワークフロー型プログラミングツールで能力ブーストは良いのか悪いのか?

プログラミングまたは機械学習の高速道路は学習の高速道路に関しての話ですが、ライブラリとかツールキットでブースト問題ってのもあります。

大体どの言語でもライブラリが充実していて、エコシステムを形成していて生産性が上がるんですが行き過ぎたプログラミング環境はどうなんかなーと思っていました。

つまりノンプログラマーのためのプログラミング環境(Pipeline Pilot/KNIME)とかRobot Process Automation(RPA)ってやつです。

ある程度の処理がコンポーネントという単位にまとめられていて、ユーザーはブロックを組み合わせる感覚でプログラミングができて、しかも処理が視覚的にわかりやすいというメリットがあります。さらに学習の初期コストが低いのでプログラミングの素養がなくてもすぐにある程度のことはできるようになります。

が、ちょっと複雑なことをしようとするととたんに複雑怪奇なワークフローになるのと、やりたいことをするためのコンポーネントがないと詰みます。

  • 前者は 結局コード書いてないだけでプログラミングはしているので、解決パターンを知らないと回り道することになる。デザインパターンを知らないプログラマーが陥るのと似た感じ
  • 後者はTensorflowのチュートリアル楽勝だけど実際自分の問題解こうとするとまったくわからんってことになるのと一緒

というわけで、ワークフローツールはどうなんかなーって思っていたのだけど、逆にノンプログラマーが利用して効率的に仕事をするのにはとてもよい選択肢だなと改めて思いました。

じゃぁなんワークフローツールをあんなに嫌っていたんだ?

よくよく考えてたら、ツールに技術がロックされることに嫌悪感を抱いていたのでした。Pipeline Pilot年間使用料がそれなりに大きいのと、すごく使いやすくてこれでしか仕事ができないケモインフォマティストが増えてしまったという製薬業界の語られない闇があります。維持費だけでも馬鹿にならないし、転職先がPP使ってないとキャリアが生きないというねw それってキャリアブルスキルなのか?っていう。一方でPPができたから良いとこに転職できたという事例もないことはないので、それがいいことなのか悪いかは一概にいえないけどね。

で、RPAツールはまさにこのサービスロックを狙っている気がします。そもそもクラス3ちらみせしてRPAのクラス1のツールに年間何百万も払わせるなんてあれなんじゃないのかな。

話がそれたけど、ノンプラグラマーなのにバイオインフォにアサインされた人間のみで構成されているうちのチームはとりあえずKNIMEブーストかけておけばいいだろうっってことでそのための学習道路の整備を急遽していますが、これはこれで色々楽しいですね。既存のKNIMEノード組み合わせてもできない処理は僕がNode作ってやればいいだけだからね。

尚、昨日始めてKNIMEデビューしました。どんだけワークフローツール嫌っていたんだ?って話w

Finding activity cliffs in ChEMBL

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.

QSP

医薬品研究開発における 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がおすすめされるようです。

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

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


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

Mentality rather than modality

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

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

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

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

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

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


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

創薬プロジェクトのITSが欲しい

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

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

ITS

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

ADMET Prediction with PotentialNet

うちの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に取り入れられるのでしょうか?

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

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

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


AI in Medicinal Chemistry終わりました

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

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

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

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

1553687389

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

1553687384 1553687387

美味しかった。

ADと特徴量

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

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

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

ファーマコフォアについてのよもやま話

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

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

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

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

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

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

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