fantom_zona’s diary

Impact the world!!!

deepchemを動かしてみました?②

久しぶりの更新です。
deepchemのチュートリアルを動かそうとしたら、タンパク質の構造を表示してくれませんでした。ここで少しつまずいてしまったので、記事にまとめようと思います。deepchemと言うかdeepchem以前なのですが、、、。

はじめに

このチュートリアルをやろうと思いました。足りないパッケージは適宜追加しながらやれば問題ないだろうと思っていたのですが、途中、

ngltraj = visualize_complex(complex_mdtraj)
ngltraj

と言う部分で、

Failed to display Jupyter Widget of type NGLWidget.

みたいな文字が表示され、複合体の構造が表示されませんでした。「あれ?NGLviewは途中で入れたはずなのにな?」と思い色々調べたことをまとめます。

結論

NGLViewを入れる際に、!pipで適当に付け足すだけではダメだったのです。

!pip install ipywidgets
!jupyter nbextension enable --py widgetsnbextension --sys-prefix

!pip install nglview
!jupyter-nbextension enable nglview --py --sys-prefix

これをコピペしましょう。終わり。
試しに、

import nglview as nv
view=nv.demo()
view

としてタンパク質が出てくればオッケーです。クリクリして遊びましょう。f:id:fantom_zona:20181104043347p:plain

医薬品開発の流れについて

基礎研究によって発見された薬剤候補は、非臨床試験で有効性・安全性が確認された後に臨床試験に移る。臨床試験は1,2,3相に分かれ、人における有効性・安全性を試験する。臨床試験を経て、国による承認を得たものが新薬として発売される。発売後も医薬品が適正に使用されるように医薬情報は収集される(育薬研究)。

臨床試験

第Ⅰ相試験

原則として少人数の健康な成人を対象として、治験薬の安全性とADME(吸収、分布、代謝、排泄)を調べる。抗がん剤などでは患者を対象とし、予備的な有効性についても検証する。

第Ⅱ相試験

 比較的少人数の患者を対象として、いくつかの使用方法(投与量・投与間隔。投与期間)を検討し、使用方法を決定する。

第Ⅲ相試験

 多数の患者を対象として、既存薬やプラセボとの比較を行う。

育薬研究

5toos

・Too few
治験は1000程度の限られた患者数で得られた情報

・Too simple
年齢や併用薬、合併症に制限が加えられた患者群の成績

・Too brief
投与期間も長期でない

・Too median-age
高齢者や小児などへの適用例は少ない

・Too narrow
腎機能・肝機能障害を合併する患者や妊婦は除外され、評価は専門家によって行われる

骨リモデリング・骨代謝における力の役割について

骨は、周囲の力学的刺激に誘導されるリモデリングによって、構造を変化させ適応的に身体を支える構造を維持している。例えば、宇宙飛行士や寝たきりの患者では、骨への力学的負荷が掛からないため、骨量の減少が認められる。このように、力学的刺激は骨構造の維持に重要な因子であることが分かっている。

骨リモデリングは、破骨細胞による骨吸収と骨芽細胞による骨形成のバランスによって制御されている。骨系細胞の一つである骨細胞は、骨への力学的負荷に応答しリモデリングを担う細胞群の活動を調節するメカノセンサーとしての役割を果たすと考えられている。骨細胞が力学的負荷を感知すると、NOやプロスタグランジンを産生することで、骨形成の方向に促進させることが知られている。骨細胞の力学刺激の形態としては、間質液流や静水圧上昇などが提案されているものの、未だそのメカニズムに関しては不明な点が多い。

 

**参考文献

CiNii Dissertations - 骨細胞の力学刺激感知およびカルシウム応答伝播のバイオメカニクス

Bone remodelling at a glance. - PubMed - NCBI

アンブロキソールの効果

今回はこの論文

アンブロキソールと言えば痰切りとして誰もが知るところだが、具体的な機序についてはそういえばよく知らない。医学生程度ではこういう機序の分からないものは結構あり、実際教科書に載っているような生理学や分子生物学に基づかない機序であることが多いようだ。こういうのは幾らでもあり、アンブロキソールを始め、カルボシステインだとかマグミットだとかミヤBMだとかかなり頻繁に使用される薬剤にも存在すると思われる。だから、アンブロキソールなんてものがautophagy-lysosome系に作用するなんてのは僕にとってかなり驚きである。

Gaucher病の原因変異として知られていたGBA1変異が、GWASの結果から、実はParkinson病の相当なリスクファクターであることは界隈では有名だが、まさかアンブロキソールがGCaseのシャペロンだなんてことは全く知らなかった。論文の内容としては、アンブロキソールがautophagyを抑制しexocytosisを促進したり、ミトコンドリアの生合成を促進するという内容。α-synuclein代謝まで言及されており、リン酸化α-synucleinが減少することからParkinson病を抑制する方向に行くのではないかと示唆されている。(僕はそれは怪しいと思うが)

discussionで機序の仮説として書かれているのは、アンブロキソールがlysosomeに移行しexocytosisやCa2+放出を促進し、autophagyの抑制やTFEBの活性化が起きるということだそうだ。(仮説を示すデータはこの論文では示されていない。)ちなみに、実際にアンブロキソールを内服したら神経細胞に影響を与えるか?という考察もされており、内服だけではこの実験の濃度には達しないだろうと書かれている。

GBAやアンブロキソールの先行研究を知らないのでなんとも言えないとこだが、非常に面白い内容だと思った。アンブロキソールがいつ発売されたのか知らないが、おそらく新しい薬ではないだろうと思う。こういう昔ながらの薬の機序は意外と面白いような気がしていて、マクロライドの抗炎症作用だとかメトホルミンのGLP-1分泌促進作用だとか意外な関係が調べてて面白い。余談だが、この前薬屋さんの説明で聞いた話では、メトホルミンにGLP-1促進作用があるので、DPP4との組み合わせが非常に良いらしい。あと、ビグアナイドを内服してる際は休薬しないと造影剤を基本的には使ってはいけないが、実のところ、メトホルミンではそんなに乳酸アシドーシスが起きないのではないかとかなんとか。データがきちんとあるのかは知らないが。でも乳酸アシドーシス怖いもんね。

THのIPFに対する効果の論文

これ。

 

肺の遺伝子発現データベースからIPFにDIO2が高発現していることを発見し、患者検体も用いて検証。DIO2はT4→T3への変換を担うことから、THがIPFになんらかの修飾を加えているだろうということでブレオマイシンやTGF-β1の肺線維症モデルマウスに対してTHを試したところ、実はTHが線維化を抑制していたという。しかも、全身投与じゃなくて吸入だけでもニンテタニブやピルフェニドンに劣らない効果を示してる。電顕で見てみたところ、ミトコンドリアが変性しアポトーシスしていたとこのことで、PINK1やPCG-1αの関与を疑い、それを立証している。

ストーリーとしてはTHは変性したミトコンドリアのクリアランスと正常ミトコンドリアの増殖を促すというところ。

 

驚きの結果。組織像が綺麗だ。血ガス所見がどうなのかが気になるところ。マウスの血ガスをしてないなんてことはないだろうから、血ガスは良いデータが得られなかったのかな?こういう意外な関係の発見、大好きなのでもっともっと探索して行きたい。

deepchem動かしてみました①

deepchemというライブラリを勉強してみようと思います。チュートリアルからまずやっていこうと思います。日本語の記事が少ないので、自分用のメモとして記録に残そうと思います。deepchemの溶解度予測のチュートリアルを一個一個読んでいこうと思います。

データの読み込み

%load_ext autoreload
%autoreload 2
%pdb off

しょっぱなのautoreloadがまず不明でしたが、このページに解説がありました。チュートリアルの方では%load_ext autoreloadはありませんでしたが、これを加えるとipython notebookではエラーが出ないようです。
%pdb offの方は不明です。pdbはデバッガらしいですが......誰かわかったら教えてください。

from deepchem.utils.save import load_from_disk

dataset_file= "./deepchem/datasets/delaney-processed.csv"
dataset = load_from_disk(dataset_file)
#ソースコードを見るとload_from_diskはcsvじゃなくても読み込める関数
dataset.head()

f:id:fantom_zona:20180322070209p:plain
読み込んだcsvの形式はこういう形のようです。今回行いたいのは、smiles形式から溶解度の予測をすることなので、それ以外のものは全て使いません。

描画

import tempfile
from rdkit import Chem
from rdkit.Chem import Draw
from itertools import islice
from IPython.display import Image, HTML, display

def display_images(filenames):
    """Helper to pretty-print images."""
    imagesList=''.join(
        ["<img style='width: 140px; margin: 0px; float: left; border: 1px solid black;' src='%s' />"
         % str(s) for s in sorted(filenames)])
    display(HTML(imagesList))
#空白にpngをjoinで足していってHTML表示する
    
def mols_to_pngs(mols, basename="test"):
    """Helper to write RDKit mols to png files."""
    filenames = []
    for i, mol in enumerate(mols):
        filename = "%s%d.png" % (basename, i)
        Draw.MolToFile(mol, filename)
        filenames.append(filename)
    return filenames
#Draw.MolToFile: mol object->pngにする
#ファイル名のリストにする

num_to_display = 8
molecules = []
for _, data in islice(dataset.iterrows(), num_to_display):
    molecules.append(Chem.MolFromSmiles(data["smiles"]))
display_images(mols_to_pngs(molecules))
#MolFromSmiles: smiles->mol object
#display_imagesはファイル名をしているすると表示してくれる関数

RDKitではmol objectというのを生成して、そこから描画ができるようです。
上の流れでは、
Chem.MolFromSmiles: smiles→mol object
Draw.MolToFile : mol object→png
として画像を出力し、display_images(filename)で画像を表示しています。
f:id:fantom_zona:20180324080224p:plain

特徴量の取得

このチュートリアルでは、ECFP4という手法で化合物の特徴量を取得しています。
化合物の局所的な構造を元にした特徴量をfingerprintというらしく、ECFPは主だった特徴量として有名なようです。

import deepchem as dc

featurizer = dc.feat.CircularFingerprint(size=1024)

loader = dc.data.CSVLoader(
      tasks=["measured log solubility in mols per litre"], smiles_field="smiles",
      featurizer=featurizer)
dataset = loader.featurize(dataset_file)

dc.data.CSVLoader()は、tasksとsmikes_filedを指定するとcsv形式データを取り込むclassですが、methodとしてfeaturizeが用意されています。featurizeは、事前に定義したfeaturizerを利用してsmiles形式から勝手に特徴量を取得してくれるため便利です。
CSVloaderで出力されるオブジェクトは、deepchemではDatasetオブジェクトとして扱われています。Datasetオブジェクトは、(X, y, w, ids) という値を持っています。Xは特徴量、yは出力、wは重み、idsはID(今回の場合はsmiles形式を返す)です。deepchemのDatasetに関するチュートリアルもあるので今後それも記事にできたらと思います。

前処理

データをtrain/splitに分け、前処理は正規化を行います。

splitter = dc.splits.ScaffoldSplitter(dataset_file)
train_dataset, valid_dataset, test_dataset = splitter.train_valid_test_split(dataset)

transformers = [
    dc.trans.NormalizationTransformer(transform_y=True, dataset=train_dataset)]

dc.splits.ScaffoldSplitterでは、分子構造の偏り?が生まれないようにうまくsplitしてくれているらしいです。他にもrandom splitなども選択できるようです。

学習部分

deepchemはsklearnをラップしています。sklearnの異なるモデルに対して、SklearnModelに突っ込めば、fit(dataset object)で学習してくれます。学習したモデルはsave()で保存してくれます。save()では、sklearnのjoblib.dumpが背後で動いてくれます。

from sklearn.ensemble import RandomForestRegressor

sklearn_model = RandomForestRegressor(n_estimators=100)
model = dc.models.SklearnModel(sklearn_model)
model.fit(train_dataset)

評価部分

from deepchem.utils.evaluate import Evaluator

metric = dc.metrics.Metric(dc.metrics.r2_score)
#r2_scoreはdc.metrics内でimportしてあるので呼び出せる

evaluator = Evaluator(model, valid_dataset, transformers)
#class Evaluator(model, dataset, transformers, verbose=False)

r2score = evaluator.compute_model_performance([metric])
print(r2score)

Evaluator.compute_model_performance()で、metricの計算をしてくれます。
適当に決めた木の数(=100)では、決定係数は0.2くらいでした。

チューニング

木の数と特徴量の数を適当に振り、最適なものを選び出します。deepchem.hyper.grid_search moduleにあるHyperparamOptでグリッドサーチしてくれます。

def rf_model_builder(model_params, model_dir):
    sklearn_model = RandomForestRegressor(**model_params)
    return dc.models.SklearnModel(sklearn_model, model_dir)
params_dict = {
    "n_estimators": [10, 20,40, 80, 160, 320],
    "max_features": ["auto", "sqrt", "log2", None],
}

metric = dc.metrics.Metric(dc.metrics.r2_score)
optimizer = dc.hyper.HyperparamOpt(rf_model_builder)
#ただのclass
#methodでhyperparam_searchが使用できる
#dictionaryでパラメータを指定すれば検索してくれる
best_rf, best_rf_hyperparams, all_rf_results = optimizer.hyperparam_search(
    params_dict, train_dataset, valid_dataset, transformers,
    metric=metric)

木は10個、特徴量の数はルートが最適でした。R2の値は0.3くらいでした。
結果を出力してみるとそこそこ学習しているかなといったところでしょうか。ただし、大小関係のみが保たれている傾向があり、実際の値の予測はだめそうです。
deepchemのチュートリアルではランダムフォレストに加えて、ニューラルネットも使っていましたが、deepじゃなさそうなのにdeepって書いてあり、しかも使っている関数がdeepchem 2.0.0では削除されているようなので割愛しました。

task = "measured log solubility in mols per litre"
predicted_test = best_rf.predict(test_dataset)
true_test = test_dataset.y
plt.scatter(predicted_test, true_test)
plt.xlabel('Predicted log-solubility in mols/liter')
plt.ylabel('True log-solubility in mols/liter')
plt.title(r'RF- predicted vs. true log-solubilities')
plt.show()

f:id:fantom_zona:20180324080257p:plain

感想

以上、csvから予測までのチュートリアルを動かしてみたという話でした。
少しだけdeepchemが使えるようになった気がしてきました。今回はグラフ構造の学習をしなかったので、グラフ構造の学習とかもやっていきたいと思いました。