次回のMishima.sykは5/28の予定です

昨晩の(遅すぎる)新年会できまりました。

内容はscikit-learnのハンズオンがメインですので、予測モデルとか作ってみたいヒトは参加するといいと思います。

尚、静岡のPythonの機運が高まってきたのでShizuoka.pyも近々予定していますので皆さんネタのほうよろしくお願いします。

さて、昨日行ったお店は三島のBeerバーでBeer Windでした。

居心地は良くて楽しかったです。

ゆずホワイト、美味い

1457774309

ミルクスタウト、超美味い。香りが良すぎ

1457774310

プリムス、無難

1457774312

ジェットスターインペリアル、そこそこ美味い

1457774313

今日はちょっと車をディーラーに見せにいったついでに、昔住んでた家めぐりをしてきましたw 徒歩圏内を3箇所くらい移り住んだので(ペット可マンションに移りたくて数百メートル移動したこともあるw)

娘とワンコを連れてよく遊びに行った公園。何もないけどw

1457774314

帰りに弥栄に寄って手火山式の塩ラーメンを食べた。

1457774316

最近色々と思うことのある日々を過ごしている。

AJACS三島のNGS解析のハンズオンに参加してきた

色々勉強になり充実した二日間でした。皆様お疲れ様でした。

こういう立場のヒトもいるってことで参加に至った背景(+α)をちょっと書いておきますね。

僕自身は完全にドライな立場のバイオ(ケモ)インフォマティシャンで、手火山式やでってくらいカッチカチの乾燥タイプです。そして、 うちの職場にはウェットなNGSの人材(普通のマイクロアレイとかDNAチップを扱う人材も)はいません(重要)

じゃぁ、なんでNGSのデータ解析すんの?っていうことになるんですが、僕は企画調査的な立場で公共のデータをゴニョゴニョしたりする必要がたまにあるんだけど、最近NGSのデータって増えてきているわけです。公共データの使い方を知らなければ「分からない☆」と無知で通せるので幸せなんですが、流石にそれは自分のキャリアを考える上でよろしくないわけです。で、年末にちょっとそういった仕事が入ってきたので、いいタイミングだなぁと年始にDRY本を読んでいたんですがところどころ実験を理解していないとわからないようなところがあって消化不良だったところにハンズオンが丁度開催されたので参加したわけです。

ProductName 次世代シークエンサーDRY解析教本 (細胞工学別冊)

学研メディカル秀潤社 / 5832円 ( 2015-10-15 )


具体的にはQCのところだったんですが、実際に実験やっている人だったら機械や手法にどういう癖があって、だからこういうあたりに注意してQCするとかそういうことを意識しながらトリミングしたりするんでしょうが、そこら辺の間隔が分からなかったので参加して非常に勉強になりました。それからウェットの参加者の割合が比較高かったせいか講師の方が「実験やっている方ならわかると思いますが」っていうあたりは押さえておくべきことがらなんだろうなとメモった。

それから、基本的に発現解析だけ押さえておけばいいかなと思っていたんだが、やはりChIP-seqも勉強しとかなあかんかなと思った。推薦されていたのはこの本だったかな?

あと DAVIDを知ることができたのが良かった。あれは便利そう。他にはMetascapeも良いらしいく、調査が捗る感ある。

最後にハンズオンやっていてちょっと気になったあたりをメモっておきます。

Pipe(|)

隣のヒトにPipeを使ったコマンドの意味を聞かれたのだけど、unix初心者にパイプって馴染み薄いよなと思った。level2で

find * | grep gtf

みたいな記述がいくつかみられるんだけど

find . -name '*.gtf'

のほうがわかりやすいかも

R関係

  • p.117のextdataPathをpasteする必要はないと思う
  • p.118はsend <- csDendro(my.genes)とするべき

ダウンロード関係

p.105でebiにwgetするという記述があって、年始に写経した時にこの部分で結局3日取られてデータ取得もDRY本のネックやなと感じたけど、これはよく考えたら(考えなくてもw)DDBJからダウンロードするべきで、DRA Searchにアクセッション番号を入れて検索すればいいとのこと。

ただ、これはちょっとステップを踏むのでスクリプトを書いてみた。

wget `python draget.py ERR266338`

こんな感じでデータが取れる

お食事

たまに僕のタイムラインに握りの盛り合わせが投稿されるw もろこしずしで地魚丼を食べた。寒すぎて握りの気分ではなかったので次回はおまかせにぎり8貫で攻めたい

1456660556

二日目の帰りにはつこでちょっと引っ掛けるかと思い寄ってみたが休みだった。facebookで告知して欲しいところ。

1456660558

その後やごみも振られて、蕎麦宗は休日だし、濃いものは嫌だとダラダラと駅に向かったら香香についてしまったw

ピータンと砂肝の唐揚げをつまんで帰った。

1456660559 1456660561

Bioinformatics版のProject Euler

RosalindというProject Eulerっぽいサイトがあったので少し遊んでみた。

楽しくbioinformaticsの問題を解きながらアルゴリズムを覚えていく感じなので、BioPythonを覚えるのにちょうどいいかなと思った。

GC Contentを計算するメソッドくらいあるだろうとは思うんだけど知らないので素朴に計算。リスト内包表記でmaxで取り出すのがいいかなと思ったけど、面倒くさいのでこれまたループで素朴に。

import sys
from Bio import SeqIO

def count_gc(seq):
    return float((seq_record.seq.count("G") + seq_record.seq.count("C"))) * 100 / len(seq_record.seq)

if __name__ == '__main__':
    fasta = sys.argv[1]
    highest_id = ""
    highest_content = 0
    for seq_record in SeqIO.parse(fasta, "fasta"):
        c = count_gc(seq_record)
        if c > highest_content:
            highest_id = seq_record.id
            highest_content = c

    print highest_id
    print "{0:.6f}".format(highest_content)

僕はBioPerlでバイオインフォマティクス関連のプログラミングをしてたので、BioPythonで配列解析をしたことは殆ど無いんですよね。

missing heritability

遺伝子と疾患(表現形)との強い相関と弱い相関が見られるが中程度の相関がみられないという謎(missing)のことらしい。

それだけ、環境要因が強いということなのだろうか?

The "missing heritability" problem[1][2][3][4][5] can be defined as the fact that individual genes cannot account for much of the heritability of diseases, behaviors, and other phenotypes. This is a problem that has significant implications for medicine, since a person's susceptibility to disease may depend more on "the combined effect of all the genes in the background than on the disease genes in the foreground".Missing heritability problem

高速文字列解析の世界を読み始めた

僕のフィールドには一応Bioinformatics,Chemoinformaticsも含まれているので、文字列だけでどこまでいけるのかは非常に興味がある。まぁ物理法則を無視できなくなってくると文字列処理ではどうしようもないんだけどね。

3章のBWTまで読んだけど、3-4のBWTの性質と復元のところがさらっと流されていてさっぱり分からなかったので、検索したらわかりやすい説明を見つけた。

が、全体的には丁寧に解説されていて、分かりやすいと思うし、読んでて楽しい。

個人的にはSMILESみたいな構造情報を文字列にしたものをBWTで扱えないかなぁと考えている(chemoinformaticsに応用してみたい)。構造変換ルールにもイディオムみたいなのあるしね。MMPを高速に検索できても嬉しいだろうし、なんか使い道がありそうな気がするんだけど。

ゲノム言語ATGC

プログラムとして実行できるfasta形式のプログラミング言語を作ってみた。いちおうチューリング完全(なはず)。

>HELLO_WORLD
ccggaccgcg gggcaccgcc ggcggaccgc cgccggaccg cccggcgacc gccgccccac
cgcccgccca ccgcggggga ccgccgcccc accgccgccg gaccgccgcc ggaccgccgg
cgcaccgcgg cgggaaccga cacatccata ccacagaacc caaaa

これはatgcというコマンドで解釈して実行します。

$ bin/atgc hello_world.fasta 
Hello world!

ゲノム的にGCリッチなほうがいいだろうということで0と1にg,cをそれぞれ割り当てて数字を表現するようにしてる。exitコマンドには終止コドンを割り当てようかとも思ったが、なんとなくaaaにしてみた(polyA)。

しかも(というか当たり前だけど)blastでホモロジーサーチがかけられるし、multifastaにしておけばソース管理もできるうえに、データベース化してインデックスはっておけば、NCBIのツール群でコマンド一発で取り出せる。

ただ今回作ったHELLO_WORLDの配列はblastnだといい感じにヒットしなくて悲しかったので、blastxかけたらブラックコットンウッドからなんかひっかかった。

>gb|ABK94795.1|  unknown [Populus trichocarpa]
Length=229

 Score = 33.5 bits (75),  Expect = 5.9
 Identities = 14/21 (66%), Positives = 15/21 (71%), Gaps = 0/21 (0%)
 Frame = +1

Query  100  GPPPDRRRTAAGTDTSIPQNP  162
            GPPPDRRRT  GT  S P +P
Sbjct  209  GPPPDRRRTRQGTTKSEPASP  22

VMとかは特にいじってないのでEsotericの本を参照のこと。

ProductName Rubyで作る奇妙なプログラミング言語 ~Esoteric Language~
原 悠
毎日コミュニケーションズ / ?円 ( 2008-12-20 )


VMを使った中間言語方式の強力さを理解した。

追記 2012.07.31

ソースをGitHubに移した

$ ./bin/atgc examples/shizuoka.fasta

                                                  **
                                                 **
                                              ****
                                              ******
                                              ******
                                            ********
                                              ****                  ******
                                            ********              ******      ************
                                          ********              **************************
                                     ****************            **********************
                              **********************        **********************
                          **************************      ************************
                        **************************          **********************
                        **********************************************************
                    ******************************************   ******************
                  ****************************************                ************
                  **************************************                ****************
                ************************************                    ****************
              **************************************               ********************
           **************************************                ************************
        ************************************                      ************************
        **************************************                    **********************
  ******************************************                    ********************
************************************************                  ********************
****      ********************************                        ********************
****     **************** **************                              ****************
**************************************                              **************
                 ****************************                                ******
                                          **********

気になる創薬関連書籍

2012年2月発刊の創薬関連書籍を見てたら面白そうなのを二冊見つけた。

バイオマーカーのマイニングの本は興味があるなぁ。

FBDDの本は事例が多めなんだけど、目新しさは感じられない。会社に買ってもらおうっと。

PyMOL v1.5をosx10.6.8に入れた

リリースしている安定版はインストールできなかったので、svnのtrunk(r3983)を入れた。それからbrewを使っているのでsetup.pyをちょっと変える必要があった。

/opt/localはmacports用の設定だと思うので/usr/X11に変更しないとGL/gl.hがないとかそんなエラーを吐くはず。

156         EXT = "/usr/X11"
157         inc_dirs=["ov/src",
158                   "layer0","layer1","layer2",
159                   "layer3","layer4","layer5",
160                   EXT+"/include",
161                   EXT+"/include/GL",
162                   EXT+"/include/freetype2",
163                   "modules/cealign/src",
164                   "modules/cealign/src/tnt",
165                   "generated/include",
166                   "generated/src",
167                   ]

それからPmwはソースをダウンロードしてきて入れた。なくても動くと思うが、import errorが出るので気持ち悪かったから後から入れておいた。

PythonでFasta fileを扱う

pyfastaってのがあった。

Biopythonでいいんじゃないかと思ったが、pyfastaのほうがメモリ効率がいいのかな。ゲノムを扱うような気がする。

あとでちゃんと見てみる。

PDBデータに化合物情報を追記して一つのファイルにしたい

mol形式(sdf形式)のデータだと化合物の区切りが$$$$なので、化合物を追加したい場合は何も考えずにファイルに追記するだけでいいのでよいですね。

PDB形式のデータにsdf形式の化合物情報をマージしたいんだけど、いい方法ないかなぁと調べてみたところmol2でOKだった。

両方mol2形式にして

cat compounds.mol2 >> protein_data.mol2

ってやればマージできる。

どういう用途を想定しているかっていうと

ある適当な部分構造(substructure)を持っている化合物の複合体結晶構造に、同じsubstructureをもつ別の化合物のコンフォマーを発生しつつ複合体のsubstructureの座標でalignする

つまり

  • obconformerでコンフォマーを発生
  • 複合体結晶構造のリガンドの部分構造を使ってobfitでコンフォマーをアライン
  • 一つのファイルにまとめてドッキングモデル完成

みたいなことをやりたかったわけです。こういうのはファイルが2つに分かれてるとユーザーのヒトとか使いにくいしどういう計算したのかわからなくなっちゃうからね。

ProductName Bioinformatics Programming Using Python
Mitchell L. Model
Oreilly & Associates Inc / 5119円 ( 2009-12-23 )