fantom_zona’s diary

合成生物学(と最近では神経変性疾患)

アンブロキソールの効果

今回はこの論文

アンブロキソールと言えば痰切りとして誰もが知るところだが、具体的な機序についてはそういえばよく知らない。医学生程度ではこういう機序の分からないものは結構あり、実際教科書に載っているような生理学や分子生物学に基づかない機序であることが多いようだ。こういうのは幾らでもあり、アンブロキソールを始め、カルボシステインだとかマグミットだとかミヤ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との組み合わせが非常に良いらしい。あと、ビグアナイドを内服してる際は48時間休薬しないと造影剤を基本的には使ってはいけないが、実のところ、メトホルミンではそんなに乳酸アシドーシスが起きないのではないかとかなんとか。データがきちんとあるのかは知らないが。でも乳酸アシドーシス怖いもんね。

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が使えるようになった気がしてきました。今回はグラフ構造の学習をしなかったので、グラフ構造の学習とかもやっていきたいと思いました。

KaggleのTitanicにトライしてみた話

右も左も分からない人間が、kaggleに挑戦してみた話です。

とりあえず登録してこのページに行ってtrain.csv, test.csvをダウンロード。
train.csvで予測モデルを構築して、test.csvを予想して提出という形。
とりあえず色々データを覗いてみる。

import pandas as pd 
import matplotlib.pyplot as plt
%matplotlib inline

train=pd.read_csv("./Titanic/train.csv")
train.head(20)

こんな感じに先頭の20個を表示してみると、Cabinに欠損値が多いことが分かる。Name、Sex、Ticket、Cabin、Embarkedは数値じゃないので、捨てるか数値に変換したい。
f:id:fantom_zona:20180321042244p:plain

ヒストグラムにしてもう少し覗いてみる前に、SexとEmbarkedを数値に変換して、Name、Ticket、Cabinのデータは捨てることにする。ヒストグラムにする際には、欠損値はdropna()で落としておく。

train.replace({"male": 0, "female": 1}, inplace=True);
train.replace({"Q": 0, "S": 1, "C": 2}, inplace=True);
train.drop(["Name", "Ticket", "Cabin"], axis=1, inplace=True);

dead_data=train[train.Survived==0]
survived_data=train[train.Survived==1]

plt.hist([dead_data.Age.dropna(), survived_data.Age.dropna()], label=["dead"," survived"])
plt.legend()
plt.titile("age")
plt.show()

plt.hist([dead_data.Fare.dropna(), survived_data.Fare.dropna()], label=["dead"," survived"])
plt.legend()
plt.titile("fare")
plt.show()


plt.hist([dead_data.Embarked.dropna(), survived_data.Embarked.dropna()], label=["dead","survived"])
plt.legend()
plt.titile("embarked")
plt.show()

f:id:fantom_zona:20180321043953p:plainf:id:fantom_zona:20180321044002p:plainf:id:fantom_zona:20180321044011p:plain

欠損値の補完は適当に中央値にしておく。

train.Age.fillna(test.Age.median(), inplace=True);
train.Fare.fillna(test.Fare.median(), inplace=True);
train.Embarked.fillna(1, inplace=True);

今回行う学習は2値分類なので、有名どころのSVMとRandomforestの比較をしてみる。
データをtrain/testに分けたいときは、sklearnのtran_test_splitが便利。
分けたあとは、ラベルと入力を分けて学習させて評価すれば終わり。

import numpy as np
from sklearn.cross_validation import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn import svm
from sklearn.metrics import classification_report, accuracy_score

train1, test1=train_test_split(train, test_size=0.2, random_state=None)
train_feature=train1.drop("Survived", axis=1)
train_label=train1.Survived
test_feature=test1.drop("Survived", axis=1)
test_label=test1.Survived

clf = svm.SVC()
clf.fit(train_feature, train_label)
svm_pred=clf.predict(test_feature)
print(classification_report(test_label, svm_pred))
print(accuracy_score(test_label, svm_pred))

rf=RandomForestClassifier()
rf.fit(train_feature, train_label)
rf_pred=rf.predict(test_feature)
print(classification_report(test_label, rf_pred))
print(accuracy_score(test_label, rf_pred))

こちらの手元では、SVMが0.6、RFが0.8弱だったので、RFを採用して提出することにした。

#testデータの前処理
test.replace({"male": 0, "female": 1}, inplace=True);
test.replace({"Q": 0, "S": 1, "C": 2}, inplace=True);
test.drop(["Name", "Ticket", "Cabin"], axis=1, inplace=True);
test.Age.fillna(test.Age.median(), inplace=True);
test.Fare.fillna(test.Fare.median(), inplace=True);
test.Embarked.fillna(1, inplace=True);

rf=RandomForestClassifier()
rf.fit(train.drop("Survived", axis=1), train.Survived)
rf_pred_all=rf.predict(test)

submission = pd.DataFrame()
submission['PassengerId'] = test.PassengerId
submission['Survived'] = rf_pred_all
submission.to_csv('submission.csv', index=False)

が、、、ここまで来てKaggleのページからどうやってsubmitすればいいのか分からなくなった。。。
調べたところ、最近になってKaggleはAPIで提出する形になったようなので仕方なくpip install kaggle。kaggleのページに行って、"My account">"create new api token"に進むと、kaggle.jsonのダウンロードが開始される。homeディレクトリに.kaggleを作成してkaggle.jsonをそこに入れたのち、

!kaggle competitions submit -c titanic -f submission.csv -m "RF all data"

で提出完了!データを落としてくるところもAPIでいけるらしい。ほー。

感想とか。
Embarkを0,1,2にしちゃダメそうな気がしたけど、まあとりあえずやっちゃった感ある。
SVMの精度がめちゃめちゃ悪かったのがなんでなのかよく分からない。ここでは結構精度が出てるので、どこかが変なのかもしれない。
前処理はもっと工夫するべきなんだろうけど、初めてのKaggleということで大目に見てもらいたい。

コネクトームをシーケンサーに掛けたい

今年読んだ一番好きな論文2017の12月21日の記事です。

初めまして。kira_kira_worldです。合成生物学と基礎医学が好きな学部生です。この企画は以前から見ていて、ずっと見る専だったんですが、今年は景品欲しさに登録してみました。

 

紹介論文

紹介する論文は、"Using high-throughput barcode sequencing to efficiently map connectomes"です。バーコード技術(後述)を使ってコネクトームを決定する方法を考えましたっていう論文です。

 

前置き

 Wikipediaの"コネクトーム"というページを引用しますと、

「コネクトーム(connectome)とは、生物の神経系内の各要素(ニューロンニューロン群、領野など)の間の詳細な接続状態を表した地図、つまり神経回路の地図のこと。」

とのことです。今回の論文の中では、コネクトームとは神経細胞の全接続のことを指しますが、コネクトームの全体像を描くのはかなり難しいようで、線虫(めっちゃ小さい生き物)とかでしかコネクトームは分かっていないようです。現状、コネクトームを網羅的に決定するには顕微鏡を用いるしかありませんが、脳全体を顕微鏡で少しずつ見ていくことはほぼ不可能なので、コネクトーム決定の困難さはこのことからもわかります。

 そこで、いい感じにハイスループットにポンポンポン!ってなんねぇかな〜〜ってことで、著者たちがいい感じの方法を思いついたみたいです。ワクワクしますね、一体どんな方法なんでしょう。

 

基本的なアイデア

 アイデアとしては、シナプス前膜と後膜それぞれを異なる配列の短いRNA(バーコード)で以って標識してやり、なんとかして一つの配列に繋げてシーケンサーに突っ込めば、各シナプス結合がタンデム配列として認識できるだろうというものです。

 具体的には、NeurexinとNeuroliginというシナプス前後に局在して手をつなぐタンパク質を利用します。RNA結合モチーフを細胞質側に融合し、細胞外にはリンカー結合部位を融合します。すると、バーコード配列はシナプスに局在し、ビオチンの付いたリンカーを加えてやればタンパク質は架橋されます。その後、細胞を溶解させ、ビオチンで免疫沈降すれば、RNAのペアが取得できます。名付けてSYNseqと言うらしいです。なんて知的でアクロバティックかつ繊細なアイデアなんでしょう!

f:id:fantom_zona:20171221210451p:plain

イデアの概略図

 

実験結果

 新しくてかなりアクロバットな手法を開発しようとしているため、細かいvalidationが多かったみたいです。論文に記述しきれていない実験が大量にある雰囲気が読んでてバシバシ伝わってきます。実験内容について簡単に説明すると、著者たちは、

 

RNARNA結合モチーフが結合して細胞膜に局在する(HEK細胞&ニューロンにおいて)

シナプスで両方のタンパクは近傍に存在している

✅ 架橋し、免疫沈降で落ちてくる

 

という確認をしたようです。(詳細は省略)

 続いて、RNA同士を結合させるという部分ですが......すでにここまでの実験でかなり大変だったみたいなのに、ここが最大の壁だったようです。著者たちは、色んな方法を試してみて(splinted ligation, ribozyme)、エマルジョン overlap&RT-PCRという方法に行き着きました。

 エマルジョンPCRとは、油中に非常に小さな水滴を作り、その水滴を小さな反応炉として利用することで、沢山の種類の核酸を別個に増幅させる方法です。通常、次世代シークエンサーのサンプル調整で使用されます(例えばこのページを参照)。

 

f:id:fantom_zona:20171221221044p:plain

 

 ということで、全体のフローがようやく完成しました。とりあえずHEK細胞でワークするか試した結果が下の図です。preまたはpostのタンパク質を発現させた細胞を混ぜたあとにシークエンスしてあります。縦軸・横軸がそれぞれHEK細胞にトランスフェクションしたバーコード配列(即ち個々の細胞)、各点がそれらの接続を意味します。

f:id:fantom_zona:20171221223740p:plain

 おお、なんかそれっぽそう!

 

 では、ニューロンでは?というと、HEK細胞での実験と似たような手法を用いて実験を行いましたが、感度が非常に悪くなってしまい、RNAをほとんど検出できなかったようです。図すらありませんでした。これはどうやらRNAの結合効率が悪かったせいのようです。著者らの実験はここで終了しています。*1

(他にも色々書かれているけど省略)

 

結論

 今回の結果ではコネクトームは読めない。

 

感想とか

 手法の正しさをどうやって証明するのか、preとpostのタンパク質は本当にニューロンで正しい局在をするのか (細胞膜全体に局在していそうな図がある)、1細胞内でcross-linkしてしまうのではないか、などと色々と疑問が多い論文でした。

 結果は全然ダメでしたが、発想においては今年読んだ論文で一番面白かったです。今回の記事ではかなり省略しましたが、色々な分子生物学的なギミックが本論文では使われており、そういうツールを使って生物学に切り込んでいくのは楽しそうだなと思いました(結果が伴えば)。僕にとっては技術的に難しい部分が多く、非常に勉強になりました。気になった方は読んでみてください!(そして間違いがあったら指摘してください。)

*1:結局、RNA結合率の良い手法を開発できなかったものと思われます。

UCSD Roy lab その3

これ

顕微鏡の自動サンプリングシステムを使って、aSynによるvesicleダイナミクスのunbiasedな解析をした論文。系が面白い。もう写真撮る過程のbiasは無いものと信じるしかないからな。「PDと小胞輸送の関連はもう前から言われていたけど、aSynはpreにいるからER→Golgiの話をするのはおかしい」とのこと。当たり前だけど、文献のアブストの内容くらいしか覚えてないようだとここら辺はあまり意識しないかもな。まあでも最近はdendriteに蛋白合成系のオルガネラは揃ってるみたいな話もあったりするから、もしかしたらpreにも系が全部揃ってる可能性はある気がする。オートファジーとかなら、どこのオルガネラがオートファジーを受けるのかってことになってくるのかな。

やっぱりこの論文の1stの方、この分野続けて欲しかったなぁ。