fantom_zona’s diary

Impact the world!!!

DGX Sparkを買ってTailscaleでSSH接続できるようにした

DGX Sparkを購入した。NVIDIAのパーソナルAIスーパーコンピュータで、Grace Blackwellアーキテクチャを搭載したデスクトップサイズのマシンだ。

Tailscaleとは

TailscaleはWireGuardベースのメッシュVPNサービス。従来のVPNと違って中央サーバーを経由せず、端末同士が直接通信する。主な特徴:

  • ゼロコンフィグ: インストールしてログインするだけでVPNネットワークが構築される
  • NAT越え: ルーターのポート開放やDDNSの設定が不要。自宅のNAT配下にあるマシンにも外出先からアクセスできる
  • MagicDNS: Tailscaleネットワーク内のマシンにホスト名でアクセスできる（IPアドレスを覚える必要がない）
  • Tailscale SSH: SSHの鍵管理をTailscaleに委譲でき、鍵の配置やauthorized_keysの管理が不要になる

個人利用なら無料プランで十分使える。

DGX SparkにTailscaleをインストール

DGX Spark上でTailscaleをインストールしてSSHを有効化する:

curl -fsSL https://tailscale.com/install.sh | sh
sudo tailscale up --ssh

--sshフラグを付けることで、Tailscale SSHが有効になる。これにより、OpenSSHの設定や鍵の配布なしで、Tailscaleネットワーク内の他のマシンからSSH接続できるようになる。

接続元のマシンにもTailscaleをインストールしておく。

SSH configの設定

毎回ホスト名を打つのは面倒なので、~/.ssh/configに設定を追加しておく:

Host dgx
    HostName dgx-spark
    User <ユーザー名>

HostNameにはTailscaleのMagicDNSで割り当てられたホスト名を指定する。確認するには:

tailscale status

これで簡潔にSSH接続できる:

ssh dgx

VS CodeのRemote SSH拡張からも同じ設定で接続できるので、リモート開発環境としても快適に使える。

『物理学者、ウォール街を往く』読了

10年以上前に一度bookoffで購入してみたものの、金融商品に関する内容が全然入ってこずなかなか読み進められないうちに途中で読むのをやめてしまった本。
なんとなくまた読み始めてみたら、今度は比較的スルスルと読み進めることができた。以前とは文章の読解能力が上がっているせいもあるかもしれない。

今アカデミアに身を置いている者としては前半までのアカデミア時代の話は、読んでいて非常にわかるものだった。
大学院時代に平たく言えばbrilliant jerk的な教授がよくない雰囲気を作り出していたり、論文が出なかったりポストがなかったりする話はいつの時代も一緒なのだなぁと思ったりした。
特に前半に関しては、その教授を半ば崇拝している大学院生がいたり、そういった雰囲気を中の人間が否定して寄せ付けない様な雰囲気を出すのは難しいことなどは、一度見聞きしないとわからない様な気がする。
後半は自分には全くわからない世界で、ゴールドマンサックスでのオプションモデルの開発について述べられている。
ここで、ゴールドマンサックスでの理論研究だけじゃなく、トレーダーとのコミュニケーションを通じて理論と実践の架け橋になろうとする姿勢は学ぶ所が多かった。
コンピューターオタクやPh.D.持ちは侮蔑の対象だったという記述もあって、アカデミアにいるとあまり感じないがこういう世界の方が広いのかもなぁと思ったり。。。
ビジネスの世界に身を置きながらも一貫して著者の学問的態度が見られ、それによる素晴らしい成果が強調されている。（競争相手に対する自らの優位性をはっきり書いていたりもする点はやや気になる。私にそれを否定する根拠はないのだが。。。）おそらく、ポスドクをやっていた時代が悪かっただけで、かなり優秀なのだろう。
著者はどうやら40そこそこで銀行に就職したらしく、そういう意味でも何か励みになる本だった。ポスドク以上のアカデミアの人に特におすすめできる。

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のオブジェクトにするパートを書けば終わり．そういえば，PythonでRNAの可視化ライブラリというと，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:今更ですが，これバイオインフォ？と不安になってきました．未定義語に漬け込んで好き放題言うマン

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

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

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

第２章はホモロジー群を導入するための長めの背景知識の説明
今は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パッケージ，神．