DeepInsightでちょっとよくわからないことのメモ

人工知能でゲノミクスをというプレスリリースでちょっとよくわからないことがあるのでメモ

最後の方の「図1 変数ベクトルxを変換Tで行列に変換する全体像と変換の具体的な手順」のところで1-aの具体的な手順としてtSNE/kPCAが提案されているが、これがよく理解できていない。

例えば化合物ライブラリの例だとそれぞれの化合物は2048次元の特徴ベクトル(フィンガープリント)を持つ。ただし二次元空間にマップされるのはそれぞれの化合物であって特徴(feature)ではない。

1-bで特徴がマップされるためには特徴自体が多次元ベクトルを持つ必要がある。同僚にN回測定のサンプルなんじゃないの?って言われたけど、それだったら平均とって終わりじゃない?ってことになった。

仮にGeneをマップしようとするばあいnサンプルを転置してベクトルにすればいいけどその場合は「訓練」「バリデーション」「テスト」にそれぞれtSNE画像ができてよくわからんことになる。

それから200x200の画像に変換するってあるんだけどデータの遺伝子が60483あるので、ピクセルに一つ一つに対応させても2万遺伝子くらいあふれるよなーと。黒く塗り潰されるか遺伝子の位置が重なって情報欠損すると思うんだけどそのあたりもよくわからん。

実装眺めるしかないなーとCode AvailabilityからURLたどって探したんだけど見つけることができなかった。

追記: コードがダウンロードできました 2019.08.09

http://www.alok-ai-lab.com/materials.php

のDeepInsight Package DeepInsight_Pkg.tar.gzだそうです。

実装はMatlabだったので手元で動かすことはできませんが、コードを読んでみました。

Cart2Pixel.mの3行目

% Q.data should be in no_of_genes x no_of_samples format

で、実際に50行目あたりで

Y=tsne(Q.data,'Algorithm','exact','Distance',Q.Dist);

となっているのでやっぱりtSNEのドットはサンプルか(転置した場合)遺伝子を表現していていて、「訓練」「バリデーション」「テスト」にそれぞれtSNE画像ができるような気がします。

Snorkeling at Ita Beach

Last weekend, I went snorkeling at the beach near my house for the first time this summer. One of the most famous diving and snorkeling spots is Ose-Zaki coast, which was already jammed with people as I had expected. I wanted to avoid heavy congestion, so I passed by Ose-Zaki coast and went to Ita coast, which is a smaller beach than Ose-Zaki.

The water and air temperatures were nice, but the clearness of the sea wasn't so nice.

It took about two hours on the way, but took about four hours the way back because there is only one road to the beach.

1565125372 1565125374

1565125376 1565125378

1565125381 1565125384

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

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

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

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

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

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

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

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

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

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

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

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

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

プログラミングまたは機械学習の高速道路

梅田望夫さんの「ウェブ進化論」という名著かつ古典があります。既に13年以上も前に出版された本ですが未読であれば読んでおくことをおすすめします。

このなかで「学習の高速道路」という話がでてきます。元ネタは羽生先生の「ITとネットの進化によって将棋の世界に起きた最大の変化は、将棋が強くなるための高速道路が一気に敷かれたということです。でも高速道路を走りぬけた先では大渋滞が起きている」っていう話で、特に学習が必要なものに関してはほぼ全てにおいて当てはまっているのかなーと思います。それは今でも変わらないかなーと。

あとはドラッカーだかだれだったか覚えてないんですが、「マネジメントがマネジメントするものは人ではなく環境である」みたいなことを言っていて、「そうよねー」って感銘をうけたので、私は基本的に学習の高速道路を引くことにこだわりがあります。優秀な人間は超スピードで駆け抜けて行くべきだし、そのためのインフラはマネジメントが用意すべきというのが私の持論であって、いつの日か離陸してほしい。

やっぱ空飛ぶ車は夢よねw

イントラGithubはプログラミング学習の高速道路としてはうまく機能しているし、Mattermostも若い人の高速道路としてうまく動いているかなーと。大渋滞のなかで頭一つ抜けるのは本人の頑張りなので、成果を出して是非その分野で存在感を出していただければなぁと思っています。

やごみとリパブリュー

やごみで反省会と次回の日程決め。次回は9/21となりました

エンガワ

1561883380 1561883371

磯辺揚げとタコ唐揚げ

1561883390 1561883375

定番のアジ刺

1561883382

週末に会場を予約した帰りにリパブリューに寄った。

1561883386 1561883388

そして3Lのビールを飲んだ。

1561883377 1561883373

台北で買った魯肉飯ミックスが残っていたので豚ひき肉で魯肉飯を作ってみた。 これはなかなかいけた。

1561883384

つみき

久しぶりに出張に行ってきたました。

昼は御茶ノ水のエチオピアに行くつもりだったんだけど、混んでたので隣のしゃりんという店に変えた。 シャリンってなんやろか?と思ったらどうも六厘舎のセカンドブランドっぽい感じでした。

あんま並ばなかったのでこっちでもいいかもでも小くらいでちょうどいいかな。並は多かった

1561882911

帰りに神田のつみきに寄ったら、店が一杯だったので、外のテーブルで飲んだ。ここはコスパ良いよね。

1561882908 1561882913

このタンおろしが美味しい。

1561882918

次の日はリパブリューでデトックス

1561882916 1561882923

四神湯と魯肉飯

四神湯と魯肉飯を買ったので早速作ってみました。

四神湯は豚モツだけを用意すればいいので湯でモツを買ってきて白い脂肪を丁寧にとって臭みを抜きました。 右の液体は日本酒に胡椒を漬けたものですが、本場の謎液はもっとディープな奥深さがありましたね。

1561882009 1561882006

四神湯は美味しかったけど、本場には敵わんのでまた飲みに行かなければと思いました。

1561882001

魯肉飯は豚バラを細かく切って炒めてから、砂糖と醤油と水を入れて煮込みつつ、ミックススパイスを投入して煮込むんだけど、豚バラだとざっくりしていて歯ごたえが出るので本場感がなかった。

1561882011 1561881998

これはこれでうまいんだけど、豚バラの脂肪分多めの部位をミンチにして使うほうがよいんかな?

1561882003

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.

ミームを信じるか否か

ちょっと最近思うことがあったので記しておく。

私に大きな影響を与えたサイトの一つとしてjoyrideがある。まぁその当時はしょっちゅう見ていた。ちなみに更新日を見てもらうとわかるけど20年以上前のサイトだ。そのうちhigh gearになってなんか心境変わったのかなーと思っていたのだけど、byebyeのページをちょっと待っていれば理由がわかります。

この出来事は自分に結構大きな影響を与えたと思う。心残りはちゃんとゲストブックに記入しなかったことかな。

「どこでも誰とでも働ける」を読んだ

台湾に行く飛行機の中で読んだのですが、これは大変良かったです。いまならPrimeReadingで無料で読めるのでサクッとチェックしておくとよいかと思います。

自分がもつ知識はできる限りオープンにしたほうが得をする

これは昔から実践していて、効果を実感している。

  • 壁をつくって知識を隠すことのメリットはどんどんなくなっています。
  • 自分の持つ知識をオープンにすると「旗を立てる」という効果もあります。

プロとしてやっていくならアカウンタビリティが必須になる

自分の仕事を自分の言葉でクライアントに説明できないといけない

  • アカウンタブルかどうかを決めるのは自分ではなく、あくまで受け手

18 会議術

会議はなかなか難しい

  • 会議の意図をメモする、共有する

会社を辞めるつもりはなくても転職活動は毎年する

いつでもやめる覚悟を持って、会社を良くしようとするから貢献できるし評価されるというのはそのとおりだと思う。最初から会社にしがみつくのが目的の人がよい貢献をするというのはちょっと考えられないかな。

  • 会社と個人はあくまでフラット
  • 会社にしがみついている人には、本気で会社を帰ることはできません

サイモン・シネックが出てたw

38 ハッカーのように課題を発見し解決を楽しむ思考法

  • 何が世の中の問題となっているのか、解決すべき課題が何かを発見する能力が重要