fantom_zona’s diary

Impact the world!!!

Cross validation + OptunaでMLPのモデル選択してみた話

やったこと

とあるテーブルコンペに出た時に, sklearnのMLPRegressorを使ってみました.
完全に素人がやってみた内容なので, 是非ともアドバイスください.

コード内容

optunaは最小化のためのパラメータ調整ライブラリらしいです.
sklearnのKFoldと合わせて使ってみました.

import pandas as pd
import numpy as np 
from sklearn.model_selection import KFold
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import r2_score


import lightgbm as lgb
import optuna

以下で, objectという関数を定義します. objectはパラメータを含む関数で, 関数内で調整するパラメータを明示的に記述する必要があります. optunaはobjectの出力を最小化します. 今回は決定係数でモデルの評価をしてみました. "X_"がinput, "y_"というのがoutputです.

kf = KFold(n_splits = 5, shuffle = True, random_state = 0)

def objective(trial):
    '''
    trial:set of hyperparameter    
    '''
    # hyper params
    max_layer_size = trial.suggest_int('hidden_layer_sizes', 50, 300)
    activation_func = trial.suggest_categorical('activation', ['relu', 'tanh'])
    alpha = trial.suggest_loguniform('alpha', 1e-5, 1e-3)
 
    # model
    model = MLPRegressor(hidden_layer_sizes=(max_layer_size, ),
                                 activation = activation_func,
                                  alpha=alpha)
    
    MLP_scores = []
    for train_index, test_index in kf.split(X_):
        X_train, X_test = X_[train_index], X_[test_index]
        y_train, y_test = y_[train_index], y_[test_index]

        model = MLPRegressor(hidden_layer_sizes=(max_layer_size, ),
                                 activation = activation_func,
                                  alpha=alpha)
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        MLP_scores.append(r2_score(y_test, y_pred))
        print("max_layer_size :", max_layer_size)
        print("activation_func :", activation_func)
        print("alpha :", alpha)
        print("result :", r2_score(y_test, y_pred))

 
    score =-np.mean(MLP_scores)  #optunaは最小化なので符号を反転させる.
    return score

関数を定義後, optunaを動かします.

study = optuna.create_study()
study.optimize(func=objective, # 実行する関数
               n_trials=20, 
               timeout=None, #default: None
               n_jobs=-1 #multiprocess
              )

PythonでRNAlibを動かすときのメモ

RNAlibのdocumentなりhelpがあまりにも欠如しているのでちまちま調べていたことを軽くまとめる。
githubのissueを見ると割と勉強になることを学んだ。

温度の設定

import RNA

print(RNA.__version__) #2.4.11

seq = RNA.random_string(20, "AUGC")
md = RNA.md()
print(md.temperature) #デフォルトで37.0 
fc_37 = RNA.fold_compound(seq, md)

md.temperature = 20
fc_20 = RNA.fold_compound(seq, md)

md.temperature = 4
fc_4 = RNA.fold_compound(seq, md)

print(fc_37.mfe())
print(fc_20.mfe())
print(fc_4.mfe())

制約条件下のfolding

制約条件の下でのfoldingの記載も見つけた。
constraint folding in python · Issue #31 · ViennaRNA/ViennaRNA · GitHub

constraint = ".......................xxxxxxx" #制約条件

sequence = RNA.random_string(len(constraint), "AUGC") #適当に生成

# create fold_compound object
fc = RNA.fold_compound(sequence)
# add hard constraints, specified by dot-bracket like string
fc.hc_add_from_db(constraint)
# compute partition function
(prob, energy) = fc.pf()

fc.hc_init() #初期化

duplexfoldの使い方について

fold_compoundを使うのが新しいRNAlibの使い方らしいが、duplexfoldみたいな古いもの未だに残っているらしい。
&を使ったstructureの出力が謎だったが、これを読んではっきり理解することが出来た。
Matching duplexfold structure to sequence · Issue #38 · ViennaRNA/ViennaRNA · GitHub

d = RNA.duplexfold('CGT', 'AAG')
print(d.structure)
d.energy
print(d.i, d.j)

#(.&.)
#2 2

(.&.) の &の前後の番号を示しているらしい。なるほど。流石にdocumentに書け。


suboptの使い方

suboptの引数の数字deltaはx0.01 kcal/molの間にあるsuboptの構造を出力してくれる

        a = RNA.fold_compound(sequence)
        solution = a.subopt(500)
        print "%d suboptimals" % len(solution)
        for s in solution:
            print "%s %6.2f" % (s.structure,s.energy)

deepchem.models.sequentialではグラフ構造を読み取れなかった話

注意

この記事を読んでも特に得るものはありません。ただdeepchem2.1.0ではsequentialにグラフ畳み込みが出来なかったというだけの話です。

概要

deepchem 2.1.0
python2.7

model=deepchem.models.Sequential()
"""
model.fit()
→ model._create_graph(self, feature_shape, label_shape)
  shapeは勝手に推論してくれる.
→ model._layer_listをfor文で回して各layer.in_layerに1個前のlayerとして定義していく. 最終的なlayerをoutputして保存.
→ TensorGraph.built()で各layerがbuild()される.
→以下略
"""

model._create_graphの中身

  def _create_graph(self, feature_shape, label_shape):
    """This is called to create the full TensorGraph from the added layers."""
    if self.built:
      return  # The graph has already been created.
    # Add in features
    features = Feature(shape=feature_shape)
    # Add in labels
    labels = Label(shape=label_shape)

    # Add in all layers
    prev_layer = features
    if len(self._layer_list) == 0:
      raise ValueError("No layers have been added to model.")
    for ind, layer in enumerate(self._layer_list):
      if len(layer.in_layers) > 1:
        raise ValueError("Cannot specify more than one "
                         "in_layer for Sequential.")
      layer.in_layers += [prev_layer]
      prev_layer = layer
    # The last layer is the output of the model
    self.outputs.append(prev_layer)

    if self._loss_function == "binary_crossentropy":
      smce = SoftMaxCrossEntropy(in_layers=[labels, prev_layer])
      self.set_loss(ReduceMean(in_layers=[smce]))
    elif self._loss_function == "mse":
      mse = ReduceSquareDifference(in_layers=[prev_layer, labels])
      self.set_loss(mse)
    else:
      # TODO(rbharath): Add in support for additional
      # losses.
      raise ValueError("Unsupported loss.")

    self.build()

debugの試み(まだ何も気づいていない時)

model=Sequential(loss="binary_crossentropy")
model.add(GraphConv(out_channel=2))
model.add(GraphGather(32, activation_fn="sigmoid"))
model.fit(train, nb_epoch=1)

以下がエラー

TypeError                                 Traceback (most recent call last)
<ipython-input-151-eccebc8e26f4> in <module>()
      2 model.add(GraphConv(out_channel=2))
      3 model.add(GraphGather(32, activation_fn="sigmoid"))
----> 4 model.fit(train, nb_epoch=1)
      5 

~/deepchem/models/tensorgraph/sequential.pyc in fit(self, dataset, **kwargs)
     78     """
     79     X_shape, y_shape, _, _ = dataset.get_shape()
---> 80     self._create_graph((None,) + X_shape[1:], (None,) + y_shape[1:])
     81     super(Sequential, self).fit(dataset, **kwargs)
     82 

~/deepchem/models/tensorgraph/sequential.pyc in _create_graph(self, feature_shape, label_shape)
    114       raise ValueError("Unsupported loss.")
    115 
--> 116     self.build()
    117 
    118   def make_estimator(self,

~/deepchem/models/tensorgraph/tensor_graph.pyc in build(self)
    665       for layer in self.topsort():
    666         with tf.name_scope(layer.name):
--> 667           layer.create_tensor(training=self._training_placeholder)
    668           self.rnn_initial_states += layer.rnn_initial_states
    669           self.rnn_final_states += layer.rnn_final_states

~/deepchem/models/tensorgraph/layers.pyc in create_tensor(self, in_layers, set_tensors, **kwargs)
   2586         b_list = self.variables[self.num_deg:]
   2587     else:
-> 2588       W_list, b_list = self._create_variables(in_channels)
   2589 
   2590     # Extract atom_features

~/deepchem/models/tensorgraph/layers.pyc in _create_variables(self, in_channels)
   2564         initializations.glorot_uniform(
   2565             [in_channels, self.out_channel], name='kernel')
-> 2566         for k in range(self.num_deg)
   2567     ]
   2568     b_list = [

~/deepchem/models/tensorgraph/initializations.pyc in glorot_uniform(shape, name)
     59 def glorot_uniform(shape, name=None):
     60   fan_in, fan_out = get_fans(shape)
---> 61   s = np.sqrt(6. / (fan_in + fan_out))
     62   return uniform(shape, s, name=name)
     63 

TypeError: unsupported operand type(s) for +: 'NoneType' and 'int'

どうやら_create_graphが出来てないみたいなので少し確認.
layerのリスト表示.

model._layer_list
[<deepchem.models.tensorgraph.layers.GraphConv at 0x1a395577d0>,
 <deepchem.models.tensorgraph.layers.GraphGather at 0x1a32d2f8d0>]

さらに、in_layer, outputsの表示. エラー的にはどこかのin_layerの定義がおかしそう?

print(model._layer_list[0].in_layers)
print(model._layer_list[1].in_layers)
print(model.outputs)
[<deepchem.models.tensorgraph.layers.Feature object at 0x1a3f89f450>]
[<deepchem.models.tensorgraph.layers.GraphConv object at 0x1a3c72f3d0>]
[<deepchem.models.tensorgraph.layers.GraphGather at 0x1a3c72f990>]

ソースコードを確認してみるとどうやらGraphConvのテンソル構築がうまくいっていないようです.

gc= model._layer_list[0]
inputs = gc._get_input_tensors(gc.in_layers[0]()) #[<tf.Tensor 'Placeholder_20:0' shape=(?,) dtype=float32>]
in_channels = inputs[0].get_shape()                       #TensorShape([Dimension(None)])

どうやらGraphConvに入れるサンプルの次元が無いよ、ということのようです。
そういえば、Sequentialのatom_featureの引数が無い......?

どうやらSequentialにはグラフ組めない......

ということで色々調べていたら、deepchem1.4まではSequentialGraphといってkearsっぽく書けたみたいなんですが、今となっては普通のnerural networkしか無理みたいですね. 以前触っていた時はdeepchem 1.4とかだったので気づいたら削除されていたようです. 悲しい. SequentialGraph絶対便利なのになんで消してしまったんだろうか。。。

deepchemのインストールについてのメモ

deepchemで遊んでいるうちに環境を壊してしまい、改めてインストールすることになりました。なのですが、これまでpython3系で動いていたはずなのにうまく行きません。deepchemのWebページでは3系でインストールできるみたいな雰囲気をかましているので、メモとして記録します。

2018年の4月くらいからこの問題は生じているみたいなのですが、ずっと解決していないようです。どうやっても3.5で動かすことができなかったので泣く泣く2.7で導入しました。

化合物-タンパク質相互作用の予測器を動かしてみた

はじめに

この論文を読んでみました。githubにコードが落ちていたのでcolaboratory上で動かしてみました。

colabでの環境構築

rdkitを入れなくてはいけないのでやむなくminicondaを入れました。
conda-forgeから入れるとこのバグは生じませんでした。

!wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
!chmod +x Miniconda3-latest-Linux-x86_64.sh
!bash ./Miniconda3-latest-Linux-x86_64.sh -b -f -p /usr/local

conda install -c conda-forge rdkit pytorch scikit-learn

実行手順

あとはgithubに書いてある通りです。
wgetでフォルダを入れた後にCPI_prediction/codeまで降ります。ちなみに、colabはディレクトリの移動は%cdじゃなくてはいけないようですので注意してください。なぜかは知りません。

!mkdir cpi
%cd cpi
!wget -c https://github.com/masashitsubaki/CPI_prediction/archive/master.zip
!unzip master.zip
%cd CPI_prediction-master/code
!bash preprocess_data.sh
!bash run_training.sh

結果

githubに書いてある通りです。
参考までに、epoch50でAUC=0.9596、Precision=0.9240、Recalll=0.8827です。
100 epoch回し終わるまでに1時間くらいだったかな?

感想

すごい。chemblなどの非常に大きなデータセットに対しても有効なのかどうか気になるところです。元データの論文を次は読んでみたいと思います。

びまん性正中グリオーマに対するCAR-T療法

この記事は、今年読んだ一番好きな論文2018の12月15日の記事です。
adventar.org

この一年は色々面白い話を見聞きしたのですが、いざ一番面白いのを選んで紹介するとなると、色々それぞれの良さがあるわけで選定するのって難しいですね*1。そういうわけで、色々悩んだ挙句、Potent antitumor efficacy of anti-GD2 CAR T cells in H3-K27M+ diffuse midline gliomasを紹介しようと思います。Bench-to-Bedside感あるやつで僕はCAR-T好きなので選定しました。

15秒で終わる紹介

・H3-K27M変異*2のびまん性正中グリオーマではGD2が高発現していることを発見した。
・GD2に対するCAR-T細胞を作成した。
・同所性に腫瘍細胞を移植されたマウスに抗GD2 CAR-T細胞を投与したところ、腫瘍細胞が除去された。
・一部マウスでは抗腫瘍効果による炎症反応で水頭症*3を発症し、致死的であった。

導入

まず、「キメラ抗原受容体T細胞(chimeric antigen receptor T cell, CAR-T細胞)」というワードと「びまん性正中グリオーマ(diffusive midline glioma, DMG)」というワードから説明したいと思います。

近年、B細胞性の急性リンパ球性白血病(B-ALL)に対する抗CD19 CAR-T細胞を始めとして、クローナル?*4な病変に対するChimeric antigen receptor T cell(CAR-T細胞)の有効性が注目されてきています。少し前にノバルティスからキムリアが出されて話題になっていたことを記憶している人もいるでしょう。CAR-T細胞のコンセプトは、細胞外領域に単鎖抗体、細胞内ドメインにCD3ζをくっ付けた受容体をT細胞に発現させることで、任意の細胞外抗原依存的にT細胞を活性化させるというものです。したがって、腫瘍細胞特異的な膜タンパクに対する人工受容体を作成・導入すれば、腫瘍細胞に対する炎症反応が惹起され、腫瘍細胞は除去されるだろうというのがCAR-T療法です*5

f:id:fantom_zona:20181217083117p:plain
CAR-T細胞のコンセプト概略図
https://www.amed.go.jp/news/release_20171107.html より拝借)

分子標的薬など治療法の発展によって様々な腫瘍の生存期間が延長している一方、DMGという腫瘍ではそうではありません。こちらのページ*6を参照して、論文紹介に必要な項目だけ抜き出すと、DMGは、
・橋を中心とした脳幹部にびまん性に発生する小児腫瘍
・現時点の有効な治療法は、緩和的な放射線照射のみ
・予後は極めて不良
という腫瘍です。DMGは橋に発生することが多く、この場合、びまん性内在性橋グリオーマ(DIPG)という名前が付いています。橋には生命維持に重要な中枢が多く存在するために、切除することはできません。DMGは小児脳腫瘍の10%を占めているにも関わらず、DMGに対する治療法はほとんど進歩していないため、長期生存が見込める様な革新的な治療法が長らく求められています。

以上を踏まえて、DMGに対するCAR-T細胞を作成すれば新たな治療戦略になるのでは?という疑問に対する研究を今回紹介します。

結果

びまん性内在性橋グリオーマではGD2が高発現している

著者たちは、まず、DIPGの検体を用いて抗原を探しました。抗体のライブラリを用いて抗原になりそうな表面抗原を同定したところ、GD2が複数の検体に共通して最も高発現でした。GD2はdisialogangliosideというガングリオシドの仲間だそうです

f:id:fantom_zona:20181215064946p:plain:h600
GD2はDIPGのマーカー候補である
また、免疫染色でもH3-K27M陽性の細胞にGD2が発現していること、H3-K27M陰性DIPG(VUMC-DIPG10)や高分化型グリオーマ(SU-pcGBM2)ではGD2の発現が弱いことが確認できました。
f:id:fantom_zona:20181215065533p:plain:w300f:id:fantom_zona:20181215065537p:plain:w300
免疫染色・フローサイトメトリによるGD2発現の確認

抗GD2 CAR-T細胞は、DIPGに対してGD2依存的に細胞障害性を示す

GD2が良い候補であったため、次にこれに対するCAR-T細胞がワークするのかを検討しました。作成したCAR-T細胞は、vitroでGD2陽性DIPGを除去することが分かりました。また、GD2依存的な活性化であることを、GD2陽性/陰性含む複数の検体やGD2合成酵素をノックアウトした細胞を作成して確認しています。人工受容体のネガティブコントロールとして抗CD19に対するCARを使用しています。

f:id:fantom_zona:20181215070145p:plain:w500
f:id:fantom_zona:20181215070158p:plain:w700
in vitroで、抗GD2 CAR-T細胞はGD2特異的に細胞障害を引き起こす
GD2は当然正常組織にも発現しているため、正常組織に対しても細胞障害性を引き起こすように思われます。今回の人工受容体と同様のデザインによる抗GD2 CAR-T療法は他の疾患で先行研究があり、ヒトでの臨床研究が行われていますが、大きな神経毒性を持つことは報告されておらず、ヒトにおける安全性はある程度担保されていると考えられるようです。

CAR-T細胞は、DIPG同所性移植マウスで抗腫瘍効果を発揮する

次にin vivoでCAR-T細胞の効果を検証するため、ルシフェラーゼを発現させたDIPGを橋に移植したモデルマウスを作成しました(ちなみに、移植してもきちんとびまん性に浸潤する)。移植後7-8週間の後にCAR-T細胞を静注、静注後40日まで(40 DPT, 40 days post-treatment)追跡を行いました。すると、ルシフェラーゼの発光が時間とともに減少していき、実際にDPT50ではH3-K27M陽性細胞が組織はほとんど認められなくなりました。

f:id:fantom_zona:20181215072510p:plain:w500f:id:fantom_zona:20181215072513p:plain:w500
抗GD2 CAR-T細胞投与マウスではルシフェラーゼの発光が減少する

実際の組織はこちら。(腫瘍はGFPを発現している)

f:id:fantom_zona:20181215080534p:plain:w600
組織像でも、抗GD2 CAR-T細胞投与マウスでは腫瘍が除去されている

CAR-T細胞は、DIPG同所性移植マウスの生存期間を延長する

最も重要な問いは、「実際にCAR-T細胞はマウスを救うのか?」というものですが、そのためにはモデルマウスが死んでくれないと困るわけです。ですが、先ほどのモデルマウスではDPT50を生きながらえてしまうようなので、著者たちは、腫瘍の増生が激しいタイプの患者由来細胞(SU-DIPG13P)を移植して、ひと月ほどで死ぬモデルマウスを作成し、これを用いて生存期間の解析を行いました。すると、CAR-T細胞によって腫瘍が除去され、生存率の改善を認めることができました。

f:id:fantom_zona:20181215075216p:plain:w700
生存解析のための短命モデルマウスの作成

しかし、CAR-T細胞によって治療したにも関わらず、死んでしまったマウスも存在します。このマウスの組織をみたところ、水頭症を示唆するような脳室の拡大が認められていました。移植部位以外での非特異的な脳実質の炎症などは認められませんでした。

f:id:fantom_zona:20181215075552p:plain:w400
抗GD2 CAR-T療法による水頭症の発症

CAR-T細胞は、他のタイプのDMG同所性移植マウスの生存期間も延長するか?

これまでは、DMGの中でも橋に発生するDIPGという腫瘍について検討してきましたが、同様の知見がDIPG以外でも成立するのかを最後に検討しました。
以下駆け足になりますが、視床と脊髄に発生したDMGでもGD2が発現していることを確認し、それぞれ同所性に移植しその後にCAR-Tを投与して反応をみました。

f:id:fantom_zona:20181215081937p:plain:w700
その他DMGでのGD2発現確認
f:id:fantom_zona:20181215082119p:plain:w600
視床DMGに対するCAR-T療法は、腫瘍を除去するものの生存率を悪くする
脊髄の場合はDIPGのように上手く行ったのですが(略)、視床の場合、きちんと腫瘍が除去されているにも関わらず、マウスの生存率は下がってしまいました(図のk)。これは、炎症による腫瘍の腫脹によって*7、脳が押しやられ逸脱する現象(テント切痕ヘルニア)によって死亡するためのようです(図のj)。腫瘍の炎症による腫脹はDIPGの場合でも認められ、それは水頭症として現れましたが、視床は特にテント切痕ヘルニアが起こりやすい位置に存在するため、腫瘍よりも脳ヘルニアで死亡する割合が増えてしまうようです。モデルは免疫不全マウスに腫瘍を移植して作成したため、実際に移植した際の炎症反応よりも抑えられている可能性があると著者らは論文内で考察しています。

まとめ

GD2に対するCAR-T療法は、DMGを除去することが分かりました。しかしながら、橋では水頭症視床ではテント切痕ヘルニアによって死亡例が出てきてしまいました。

感想とか、自由記載欄的な

DIPGというとmiserableな予後の腫瘍という印象がすごく強かったんですが、論文を読む限り割と希望が持てそうな感じだったので感動しました。サイエンスとしての面白さという点ではもっと面白いことしてるのは沢山あったんですが、実際に医療にインパクトを及ぼしそうな印象を受けたので今回紹介しました。
水頭症・脳浮腫対策とか取ったらどれだけ生き残れたんですかね。考察では、シャント作ったり減圧すると良いんじゃないかと書いてありましたが、気になるところなので早く続報が出てきてほしいです。

*1:あと、他の日本語メディアで紹介されてるのを避けてくと割と減ってしまう。仕方ないけど。

*2:この変異が多く8割ほどで見られる

*3:脳に髄液が貯まる病気。炎症による腫瘍周囲組織の腫脹で第四脳室あたりが閉塞し、髄液の循環が阻害されたため?

*4:自己免疫性B細胞など(https://www.ncbi.nlm.nih.gov/pubmed/27365313)

*5:所謂、がん免疫療法です。

*6:めちゃくちゃ詳しい。読むべし。

*7:免疫チェックポイント阻害剤使用時でも、"pseudoprogression"と言われて似たような現象が起きることが知られているらしい。

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