fantom_zona’s diary

Impact the world!!!

pubmedGPT動かしてみたメモ

これで動いた.
colabで動かそうと思ったけどモデルが重すぎてcolabでは動かなかったが,他の環境で動くことを確認した.
一応colab用のコードを貼っておく.

!wget -c https://repo.continuum.io/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 -q -y -c conda-forge rdkit python=3.8
import sys
sys.path.append('/usr/local/lib/python3.8/site-packages')

!conda create -n pubmedgpt python=3.8.12 pytorch=1.12.1 torchdata cudatoolkit=11.3 -c pytorch
!conda activate pubmedgpt
!pip install -r finetune/setup/requirements.txt

!git clone https://github.com/stanford-crfm/pubmedgpt
%cd pubmedgpt/finetune # sys.path.appendの方が良いかも

あとはgithubにあるものと同じようにして動かせばOK

import torch

from transformers import GPT2LMHeadModel, GPT2Tokenizer

device = torch.device("cuda")

tokenizer = GPT2Tokenizer.from_pretrained("stanford-crfm/pubmed_gpt_tokenizer")

model = GPT2LMHeadModel.from_pretrained("stanford-crfm/pubmedgpt").to(device)

input_ids = tokenizer.encode(
    "Photosynthesis is ", return_tensors="pt"
).to(device)

sample_output = model.generate(input_ids, do_sample=True, max_length=50, top_k=50)

print("Output:\n" + 100 * "-")
print(tokenizer.decode(sample_output[0], skip_special_tokens=True))

RNA描画ツールRNAcloudを作ってみた話

この記事はkyoto.bioinfoアドベントカレンダーの企画でお送りしております.
adventar.org

TL;DR

低分子集団をいい感じにプロットするmolcloudを改変して,RNAcloudという描画ツール作ってみた.
RNAcloudで君だけのRNAを描画しよう!!(molcloudのオプションとして使用できる)
github.com

動機・目的

発表スライドや資料に研究対象のいい感じの図を貼りたいことが往々にしてある.私の場合,RNAをやっているので,RNAの集合を描きたいのだが,手でちまちまやるのは結構めんどくさい,もういっそのこと自動的に作るツールを作りたいなと思っていた.単一RNAの可視化ツールは色々あるが,集合にするのは結構めんどくさそう...と思っていた時,Andrew White先生がmolcloudという描画ツールを作ったというツイートが流れてきた.


「おお,こんなのが欲しいな,というかもうこれじゃん!」と思い,molcloudを改変してRNAcloudを作ることにした.

インストール・使い方・出来上がり図

mocloudのオプションとしてインストールしてください.

pip install molcloud[all]
rnacloud fasta_file

fastaファイルは以下の形式で入力してください.

>sequence_id
AGUGUGUCGAUGCUAGUCGAUGC
...((((((((....))))))))

出来上がり図はこんな感じです.構造を取っているRNAが多いともっと見栄えがいいはず.

やったことの説明

そもそもmolcloudが何をしているのかと言う話に入らないといけない.中身を覗いてみると,どうやらnetworkxというpythonのライブラリを使って低分子をグラフに変換して描画をしている.networkxでの操作を手短に説明すると,グラフの集合から和を取って一つのグラフオブジェクトにして,ノードの初期位置ランダムでプロット,その後にノード間に斥力を入れて位置を調整するらしい.

なので,RNAをnetworkxのオブジェクトにするパートを書けば終わり.そういえば,PythonRNAの可視化ライブラリというと,Vienna大学TBIが出してるforgiが思いつく.こちらではどうなっているのだろうと思い見てみると,なんとnetworkxのオブジェクトを作成する関数が実装されているではないか!*1 これで殆ど話は終わりで,塩基対と共有結合を区別するようにしただけ.

感想・その他

正直,molcloudの発想が優れていて,それに相乗りしただけ.
RNAcloudに変更するにあたって殆ど書いてないが,それなりに満足のいくものができて嬉しい.
今回,この記事を書くにあたって色々とRNAcloudの細かいアップデートを行いたいなと思ったのだが,余裕がなくてできなかった.気になる点は色々とあるので,追々アップデートしていきたい.*2
大きな変更点で言うと,RNAの図をもっと良くしたい.R2Rという描画ソフトっぽくしたいのだが,molcloudの枠組みに乗って実現することはなかなか難しそう.
単純な発展版として糖鎖版のglycocloudも簡単にできそう.糖鎖はtokenの種類と対応する図のルール (Symbol Nomenclature for Glycans, SNFG) が決まっているのでこちらの方が一直線だと思う.
その他,改善点や何かこうした方が良い等と思われる点があったらご指摘して頂きたい.機会があったら是非使って欲しい.

*1:この関数はforgi内部ではあまり使われていないので,本当に偶然か何かっぽい

*2:例えばmolcloudは画像のmaskingができるのだが,RNAcloudではできない.あとはSHAPE reactivityみたいな外部スコアを色で表現したい.

bioanalyzerのファイルパーサーを作った話

これはkyoto.bioinfoアドベントカレンダーの記事です.
adventar.org

TL;DR

bioanalyzerの結果を可視化するためのpythonスクリプト(PyBioAnalyzer)を書いた.
github.com

追記・編集履歴

pipで使えるようにしてみました.初めてsetup.py書きました.*1
それに伴って,一部編集しました.

目的・動機

用意した核酸がNGS等あとの実験に大きく影響を及ぼす場合,電気泳動装置を使ってサンプルがキチンとできているのかを確認することが多いです.電気泳動装置の中でそこそこ有名なものとしてagilent社のbioanalyzerというのがあり私はよく使っています*2マイクロチップ電気泳動というものらしいです.bioanalyzerの詳しい説明はなんとagilentの人が日本語の解説記事を過去に書いていたのを発見したので,そちらを参照ください.
www.jstage.jst.go.jp

さて,単一バンドの場合ならbioanalyzerで善悪がすぐにわかるのですが,長さの分布が複雑なライブラリの場合目視で確認するのは結構しんどいです.bioanalyzer付属のソフトウェアで可視化して検討することも可能ですが,付属のPCでしか扱えないし,なんだか使いづらくて何ができるのかよくわからないです.
ということで,pythonで好きにカスタムする方が良いかなと思い,PyBioAnalyzerというのを作ってみることにしました.

インストール

pip install pybioanalyzer

で入ります.やったね!

PyBioAnalyzerができること

pybioanalyzer --folder_name example/ --assay_type HS_DNA --min_lim 30 --max_lim 500 --disable_ladder

とすると,以下のようなplotが出てきます.

自分で色々カスタムしたい場合には,以下のような形で動かしてください.
また,一応ladderがまともに流れているかどうかを確認することもできます.

import matplotlib.pyplot as plt 
from pybioanalyzer import PyBioAnalyzer

pba = PyBioAnalyzer(folder_name, assay_type)

# 時間と長さの関係を確認
pba.linearity_check()
plt.show()

#結果のtable表示
pba.summary

# plot
pba.plot_samples(pba.BAfiles, [args.min_lim, args.max_lim], ladder = args.disable_ladder)

# 自分で追加
plt.~~~~

bioanalyzerの出力ファイルに関して

実のところ,ここまでで話は終わってしまっているのですが,流石にもうちょっと何か書こうと思い,bioanalyzerの出力形式でも書いてみることにします.
キャピラリー電気泳動なので,測定されているのは「どの時間に,どのくらいのシグナルが得られたか」という情報です.
時間を核酸の長さに変換するために,ラダーの情報が用いられます.
ラダーはbioanalyzerに付属のもので,ソフトウェア内部でラダーのピークの検出とピークと対応するべき長さのアノテーションが自動的に付けられるのだと思われます.
付属ソフトウェアから出力されるフォルダの構成は以下の通り.

.
├── 2100\ expert_High\ Sensitivity\ DNA\ Assay_DE13804795_2019-10-10_19-27-39_Ladder.csv
├── 2100\ expert_High\ Sensitivity\ DNA\ Assay_DE13804795_2019-10-10_19-27-39_Results.csv
├── 2100\ expert_High\ Sensitivity\ DNA\ Assay_DE13804795_2019-10-10_19-27-39_Sample1.csv
└── 2100\ expert_High\ Sensitivity\ DNA\ Assay_DE13804795_2019-10-10_19-27-39_Sample2.csv

Ladder, Sampleはサンプルごとのシグナル情報を含んでおり,Resultsには各サンプルのまとめが記載されています.内部ではこれらを読み込んでそれぞれ必要な情報を頑張って抜き出しています.
Resultsには,濃度などの情報が細かく記載されているのですが,今回はladderのSize-Timeの情報しか用いていません.
濃度情報を拾ってplotに反映させることも一瞬頭をよぎったのですが,実用的にはintensityの情報だけで十分なのでひとまずおいておきました.

ポエムなど

このスクリプト自体は以前書いたもので,今回の記事を書くにあたって見直したところ,今見ると色々と気になるところが出てきてすごく書き直したくなってしまいました*3

まあこのツールはちょっとしたものなんですが,予想外に満足感が結構得られてしまって,こういう身の回りのちょいツールを作るのって時間対心理効果のコスパが良いなと思いました*4. これで気分が良くなって,ラボにQseq100というbioanalyzerに似た別のキャピラリー電気泳動が導入されたこともあり,こちらのファイルのパーサーも一応作成してみたりしましたが,色々なランモードが設定できるのとファイルの形式がややこしくて対応するのが大変で中途半端になっています*5. ファイル形式が扱いやすいと助かりますね.

もし使ってくれたら嬉しいです.また,今回できていない機能や改善点・問題点などあったらプルリクなりissuesに投稿するなりして頂けたら幸いです.*6

COI開示

agilent社, BiOptic社とのいかなる利益関係もありません.

*1:たのしい.

*2:tapestationというのも研究所にあるようです

*3:でも面倒だったのであまり困ってないので少し手直ししたくらいにとどまった

*4:自分へ.これは研究ではないので研究で満足しましょう.(つらい)

*5:それはそうと,Qseq100はbioanalyzerに対して短いDNAの波形が細かく見れてすごい

*6:今更ですが,これバイオインフォ?と不安になってきました.未定義語に漬け込んで好き放題言うマン

『タンパク質構造とトポロジー』読書記録①

最近データ解析に使えるかと思って共立出版から出ている『タンパク質構造とトポロジー』を読んでいるのでその超簡単なログ.
データ点の幾何的構造を抽出するためのパーシステントホモロジーという道具のための本.
数学に触れたのがかなり久しぶりのため,だいぶ忘れているが少しずつ思い出してきた気がする.
それと以前勉強していた時には何をしているのかよくわからないところもあったが,今思うとその時よりもやりたいことが明確に理解できているような気がする.その逆はあまりない気がするので,そういう認知バイアスがあるのだろうか?

第1章はまず単体複体を導入して最後にどんな感じのことがやりたいのかを感覚的に説明.
脈体定理を使って,図形とホモとピー同値な各種単体複体の構成方法を考える.脈体定理が非常に強力.証明大変そう.

第2章はホモロジー群を導入するための長めの背景知識の説明
今はR加群のところを読んだ.Z加群を導入するための代数学が簡単にまとまっていて復習にはちょうどよかった.
有限生成自由R加群のところまで行ってrankが定義できれば,あとのノリは線形代数と基本的に同じっぽい.(PIDになっている必要がある?)
また,ホモロジーのことを思い出すために『広がりゆくトポロジーの世界』の2章も軽く眺めてみた.
よくわからない単語が多かったが,幾何構造を調べるために部分構造同士の境界の関係を調べるというお気持ちがホモロジーに含まれているらしい.
なるほど?

初めての焼きなまし法やってみた

TL;DR

pythonのsimannealパッケージを使うと一瞬で実装できるます.
作ってくれた人ありがとうます.

背景

ラボメンに,遺伝子パーツの選択の相談をされた.
30種くらいの遺伝子パーツの中から,10個くらい良さげなペアを選んでほしいとのこと.
各々の遺伝子パーツは,入力分子とそれに応答する素子で構成される.応答素子は対応する入力分子以外にも応答する可能性があるため,この応答交差性を最小限にするようなパーツの組み合わせを選択したいという問題になる.
選択基準としては,遺伝子パーツ間の距離が定義されていて,選択したものの間の距離の最小値を最大にするというタスクになりそう.
ここcs.stackexchange.com
に同じような相談があったのでそれを参考にした.
DPでできるかなとも思ったが,よくわからなかったDPの更新式がめんどくさそうのでとりあえずDPはやめた.(できた人はおしえてください.)(1つメンバーを変えると,それ以外のメンバーの間の距離も変わってしまうため,DPには向かないタスクと思われた.)
全探索したらどうなるのか軽く概算したところ10^8 ~ 10^9くらいの計算量になるっぽい.
とりあえず近似解で良いかなと思い,焼きなまし法を選択することにした.

やったこと

pythonのsinannealパッケージがどうやら便利とのことなのでそれで実装した.確かにすごい便利.
バージョンとか書くべきなんだろうけど,とりあえずsimanneal-0.5.0使いました.
Annealerクラスを継承して,energyとmoveを書いてやればおしまい.
一部詳細を省いているが,以下のような感じになった.超便利.
やはり,subsetの組み合わせを計算しなければいけないので,結構時間が掛かる.
何回かtrialを重ねて最も良いsubsetを選択するような形でなんとか行けた.

from simanneal import Annealer
import random
from scipy.spatial import distance

class FindSubsetProblem(Annealer):
    def __init__(self, data, state):
        super(FindSubsetProblem, self).__init__(state)
        self.state = state
        self.data = data
        self.state_data = hogehoge # dataからstateに属するsubsetの選択. 

    def pair_energy(self, g1, g2):
        return # subsetを選択したときの, g1-g2距離の定義.
    
    def energy(self):
        min_dist = float("inf")
        for i in range(len(self.state)):
            for j in range(len(self.state)):
                if i > j:
                    d = self.pair_energy(self.state[i], self.state[j])
                    if d < min_dist:
                        min_dist = d
        return min_dist

    def move(self):
        E_pre = self.energy()
        rm_i = random.randint(0, len(self.state)-1)
        add_gene = random.choice(list(set(self.data.index) - set(self.state)))
        self.state[rm_i] = add_gene
        self.state_data = hogehoge # initでも書いているが, ここで改めて記述しないと更新されない仕様になっている. 
        E_post = self.energy()
        return E_pre - E_post

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('--input', '-i')    # 必須の引数を追加
    parser.add_argument('--n_members', '-n', default = 13, type = int)    # 必須の引数を追加
    parser.add_argument('--steps', type = int)    # 必須の引数を追加
    parser.add_argument('--random_seed', default = 42, type = int)    # 必須の引数を追加
    args = parser.parse_args()

    random.seed(args.random_seed)

    fsp = FindSubsetProblem(df, init_state)
    fsp.steps = args.steps
    fsp.copy_strategy = "slice"
    state, e = fsp.anneal()
    print("Selected members", state)
    f = FindSubsetProblem(df, state)
    print("Minimum distance", f.energy())

結語

pythonのsimannealパッケージ,神.

condaでRの環境構築+jupyter notebookでRを起動

目的

タイトルの通り. devtools周りで途中で躓いたのでメモしておく.

TL;DR

conda create -n r_env r-essentials r-base
conda activate r_env
conda install 
r

以下, R.

install.packages(c('repr', 'IRdisplay', 'evaluate', 'crayon', 'pbdZMQ', 'uuid', 'digest'))

conda経由でdevtoolsをinstall.

conda install -c conda-forge r-devtools

再びRでIRkernelをinstall.

devtools::install_github('IRkernel/IRkernel')
IRkernel::installspec() # install for the current user
IRkernel::installspec(user = FALSE) # install system-wide

これでjupyter notebookからRが使用可能.

0. 環境

OS : macOS High Sierra 10.13.6
conda version : 4.8.3
R version : 3.6.1 (2019-07-05)

1. condaでRの環境を構築

この記事を参考にcondaでRの環境構築をした. しかし, Rをinstall/uninstallを繰り返すと以下のようなエラーが出てきてしまい, 面倒だった. 出てきた際には, 環境丸ごと削除・作成を繰り返した...(疲弊)

dyld: Library not loaded: @rpath/libreadline.6.2.dylib 
    Referenced from: /Users/path/to/envs/r_env/lib/R/lib/libR.dylib 
    Reason: image not found 
Abort trap: 6

2. R kernelをjupyterに導入

googleで検索するといくつも記事(例えばこの記事)がヒットしたが, devtoolsの導入時にエラーが出てしまい, とても困ってしまった.

checking whether the C compiler works... no
configure: error: C compiler cannot create executables
ERROR: configuration failed for package ‘git2r’
ERROR: dependencies ‘usethis’, ‘git2r’ are not available for package ‘devtools’

どうやら, conda環境のinstall.packages経由ではdevtoolsを入れられないらしい....
色々模索したが, なぜCのコンパイラでエラーが出るのか正直よくわからなかった.
解決策と思われるものもいくつか見つけたが, 環境全体が汚れてしまうんじゃないかと心配だったので他の方法を調べた.
そしたら, conda経由でdevtoolsが入れられるようだったのでこちらを採用した.
condaなら環境をよしなにしてくれるだろう......(祈り)

conda install -c conda-forge r-devtools

Rで以下を実行.

devtools::install_github('IRkernel/IRkernel')
IRkernel::installspec() # install for the current user
IRkernel::installspec(user = FALSE) # install system-wide

これで, jupyter notebookからRが使用できた.

構造アラートに関するメモ

目的

化合物を選択する際に避けたい「忌避構造」を定義するが「構造アラート」である。
これまでに色々な構造アラートが考案されていているが、体系的にまとまった資料が欲しいなと思ったのでそのメモ。
構造アラート含めた化合物選択に関する総論的資料としてはこれが良さげだった。

ChEMBLの資料

"ChEMBL & Structural Alerts"とググると出てくるpdf。有難や...。
特に最後の方の追加資料がかなりありがたい。

"Alarms about structural alerts"

タイトルでググると論文と発表スライドが両方出てくる。ありがとう。