2011-08-28
■homebrewでnumpy/matplotlabのインストール
以下を参照．パーミッションエラーが起きるけど，無視して進んでOK．パーミッションを正しく(?)設定していると，brew linkをしてくださいと言われるので，従う．
http://www.thisisthegreenroom.com/2011/installing-python-numpy-scipy-matplotlib-and-ipython-on-lion/
- Homebrewのインストール
/usr/bin/ruby -e "$(curl -fsSL https://raw.github.com/gist/323731)"
- Pathの設定：
export PATH=/usr/local/bin:$PATH
- Pythonのインストール
brew install git brew install readline sqlite gdbm pkg-config brew install python --framework --universal
利用されるpythonのバージョンを変更する
cd /System/Library/Frameworks/Python.framework/Versions sudo rm Current ln -s /usr/local/Cellar/python/2.7.2/Frameworks/Python.framework/Versions/Current
- pip&numpyのインストール
easy_install pip
- (注2) which pip をして，/usr/local/share/python/pip を指している事を確認する
pip install numpy
- scipyのインストール
brew install gfortran pip install -e git+https://github.com/scipy/scipy#egg=scipy-dev
- matplotlibのインストール
pip install -e git+https://github.com/matplotlib/matplotlib#egg=matplotlib-dev
- ipythonのインストール
pip install ipython brew install pyqt
brew install zmq pip install pyzmq pip install pygments
- 起動！
ipython qtconsole --pylab=inline
こんな感じになります．
エキスパートPythonプログラミング
2011-08-20
■[python]RとPython(NumPy)の対応
- 対応表
- NumPy for R users: http://mathesaurus.sourceforge.net/r-numpy.html
- Rのdata.frameに相当する事を実現するには
- 作成：NumPyのarrayで，各列にclassを指定できる (structured arrays)
> import numpy as np > C0=[1,2,3] > C1=['M','F','M'] > C2=[3.3,4.4,5.5] > x=np.array(zip(C0,C1,C2),dtype=[('id', '<i4'), ('gender', '|S1'), ('value', '<f8')]) > x[1] (2, 'F', 4.4) > print(x[x['gender']=='M']) [(1, 'M', 3.3) (3, 'M', 5.5)]
- Viewを切る
>>> print(x[[1,2]]) [(2, 'F', 4.4) (3, 'M', 5.5)] >>> print(x[['gender','value']]) [('M', 3.3) ('F', 4.4) ('M', 5.5)] >>> print(x[[1,2]][['gender','value']]) [('F', 4.4) ('M', 5.5)] > print(x[[1,2],['gender','value']]) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: invalid literal for long() with base 10: 'gender' > print(x[[1,2],[0,2]]) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: too many indices for array
- 最後の2つが動いてくれないのがいまいち．
- Join:xとyを特定のキーを基に結合
>y=np.array(zip([1,2,4],[6.6,7.7,8.8]),dtype=[('id','i4'),('height','f8')])
>from numpy.lib import recfunctions as rfn
>rfn.join_by('id',x,y) # inner join相当 = xとyの両方にあるもののみ返す
>rfn.join_by('id',x,y, jointype="outer") # outer join相当 = xとyのいずれかにあるものを全て返す．埋まらない値は"-"が入る
- 列結合：cbind相当(c_)....どうやるんだろう．vstackとかhstack動かないし，c_も動かない
- ファイルの読み込み：
dtype=[('id','i4'), ('height', 'f8')]
ary=np.genfromtxt(fil, dtype=dtype, skip_header=1)
- csvなら，delimiter="," をgenfromtxtのオプションで付ける
- header行がなければ，skip_header=1は要らない
>print "\t".join(ary.dtype.names) >for i in range(ary.shape[0]): > print "\t".join([str(x) for x in ary[i]])
- もう少し複雑にviewを切る
- 値での絞り込み
>ary[ary['id']==1]
>tf = [x in ["M","F"] for x in ary['gender']] >ary[np.array(tf)]
2011-03-30
■[research]お茶大から東工大に3年間レンタル移籍します
エープリルフールネタではないので，少しフライングで公開します．
3月31日で丸5年間お世話になったお茶の水女子大学を退職し，4月1日より東京工業大学 大学院情報理工学研究科 計算工学専攻の准教授として赴任する事になりました．ただし，通常の異動ではなく，お茶の水女子大学と東京工業大学との間で結ばれた教員の人材交流協定に則っての「転籍出向*1」となります．期間は3年間．2011年4月から2014年3月までです．新しい環境で，気持ちも新たに研究及び教育をしていきたいと思いますので，今まで同様，そして，今まで以上にお付き合いくださいますようよろしくお願い致します．
- 修士，博士の学生の受け入れ(東工大の学生として)も行いますので，瀬々研で研究したい学生さんはコンタクトをください．場所は大岡山キャンパスです．
- 各所先生方には，住所や電話番号など決まってから連絡しようと思っておりますので，先にブログ上での公開をお許しください．
さて，今回の異動は大学同士の人材交流協定に則ったものとはいえ，様々な人のご好意があって初めて実現したものでした．
まず，研究室の学生．本件を告げたときに多くの学生は狐につままれた様な顔をしていました．教員及び研究者は異動が多いとは言え，そのような事実に対して学生に免疫があるとは思われず，自分の研究室が突然無くなるという事態を把握するのに時間がかかった事と思います．また，本件を告げる事ができたのが卒論，修論の追い込みの時期，就活スタートの時期でもあり，精神的にも大変辛い決断を強いてしまいました．私の異動により，学生に精神的な負荷がかかることは十分予想してましたが，それでも異動を決断できたのは，研究室の学生が他の研究室でも十分に研究していけるほどシッカリとしたメンバーだと確信できていたからです．研究室のメンバーに恵まれていたことのには感謝してもしきれません．卒業して社会に飛び立つメンバーを除き，大学院に在籍するメンバーは，ある者は私について東工大で研究をし，ある者は他の研究室へと移って行きます．一人ひとりが最高の大学院生活を送れることを願っていますし，今まで以上に出来ることはサポートを続けたいと思います．
次に，お茶大の同僚の皆様*2．学科運営，副専攻運営等，本当にぎりぎりの人数で（あるいは，現時点でも人数が足りない状況で）運営している状況で，私ひとりが抜けることで，私の担当していた仕事を他の先生に任せる形となり，様々な先生方にご迷惑がかかる事になりました．どの先生も快く仕事を引き受けた上で「頑張ってきてください」という言葉をかけてくださいました．その言語を受けるたびにうれしさと共に非常に申し訳ない気持ちが溢れました．同僚の先生方は本当に皆様良い人ばかりで，本当に本当に良い人ばかりで，積極的に持分以上の仕事をしてくれます．これ以上同僚に恵まれる環境は無いのだろうな，と思います．今の教員メンバーなら，お茶大（特にお茶大情報）は，今後も学生が甘えながらも伸び伸びと勉強や研究できる環境が引き継がれていくと確信しています．
三番目に，東工大の先生方，事務の方．お茶大でPIをする事も，私の能力を遥かに超えたものだと今でも思っていますが，それ以上に東工大でのPIも能力を超えたもので，登用してもらって有り難い気持ちだけでなく，怖さもあります．そのような私が，更に出向という，ある種の外様の者に対して通常の教員と同様に扱えるよう手続きを行っていただいています．おかげで本業である研究教育以外の細々した不安少なくスタートが切れそうです．受け入れ態勢に関わってくださった，全ての先生方，事務の方には赴任前から大変感謝しております．
最後に，お茶大の3年生．私の研究室を志望してくれていた学生が何名かいた事も聞いていましたが，突然異動に伴い研究室が志望できなくなる事態になっても，それほどパニックを起こさず，研究室の配属を決めてくれました*3．研究室を選ぶ事は人によっては「人生をかけた」重要な事項と考えるかもしれませんが，お茶大情報の先生方は誰をとっても面倒見よく，学生を大切にしてくれますので，多少希望と異なる研究室へ行っても良い研究ができると思います．各自，目の前の研究に邁進してくれることを期待しています．研究会や国際会議で会いましょう．
以上の様々な問題が起きることを覚悟した上で，やはり異動をしようと思ったのは，環境を変えることでもう一段"上"の研究を目指そうとおもったからです．お茶大にいて「変わらない」なら今まで通りの研究成果を残していけるとおもいます．5年間で今までの学生が積みあげてくれた知識や設備の蓄積もあります．これらを壊すことになるのは怖いですが，変わらないで後悔するより，変わって後悔したいので，新しい環境で，新しい仲間と研究し，今までより質高く，今までより幅広い研究にチャレンジしてみようと思いました．
非常に運のよいことに，Human Frontier Science ProgramのYoung Investigator Grantも採択されました．こちらの研究期間も3年間と東工大に在籍する予定期間と同じです．幸先よい出だしとなりました．東工大の在籍中は，こちらのグラントの採択内容をはじめとして，今まで通り次世代シークエンサやマイクロアレイを用いた遺伝子発現の解析，ネットワーク解析，生命情報に関するデータマイニング手法の開発などを行っていきます．今までの共同研究も引き続き行っていきますので，よろしくおねがいいたします．
おまけ：
- 実は場所の関係で，物理的な引越しは3月末ではなく4月に入って暫くしてからになりそうです．というわけで，しばらくはお茶大と東工大の両方に出入りしています．
- 企業の方々との共同研究も今まで通りwelcomeです．寄付金もお待ちしております:)
2011-03-17 積ん読の解体
■闘う！ ウイルス・バスターズ 最先端医学からの挑戦 (河岡義裕，渡辺登喜子)
ウイルスで有名な河岡先生と渡辺先生の本．今，研究者に進むか，会社に務めるか迷っている人（特に生命系）が居たら，読んでみることをおすすめします．研究者が必ずしも小さい頃から「その道」を選んでいるわけではなく，環境によって，色々進路が分かれてきた様子が見て取れます．特にII部の3章のネイチャーに載るまでの紆余曲折や，II部4章の研究者による研究者の紹介（人生談）は楽しかったです（この章に限らず，渡辺先生の章は，旦那さんとの会話など生活感があって，和みます．女性で研究者を考えている人は参考になる所があるかもしれません）．
「おわりに」の所で，河岡先生が実際に執筆した人の名前を出しているのは，多くの本がゴーストライター（と呼んで良いのか分かりませんが）を表に出さない中，河岡先生の真摯な態度が出ていて，好感がもてました．
■「こころ」は遺伝子でどこまで決まるのか―パーソナルゲノム時代の脳科学 (宮川 剛)
ゲノムと脳の関係は，私自身の興味でもあり，楽しく読めました．
また，著者の宮川先生の研究で作成されている「網羅的行動テストバッテリー」は，私が大学院生だったか，助手時代だったか忘れましたが，参加させていただいた科学研究費の会議で，その話を拝見し「これはすごい！」と感動したのを覚えています．多くの発表が遺伝子ひとつに着目する中，「データを網羅的に取る」という意味の重要性を本当に理解している研究者なのだと，その時思いました．
本書の中は，ゲノムと行動の関係について，入門からスタートし，体験談を含めつつ，分かっている事を順に説明しています．前半が入門，中盤は23andmeによる検査の話を挟みつつ，後半はSNPs(一塩基置換)と行動に関して，分かっている事を列挙している形です．これから，まだまだ分かることは増えていくでしょうが，現時点でのパーソナルゲノム時代に向けて分かっている事を整理してある本として，良書であると思います．
■ビジネスのためのデザイン思考
この本と次の本はビジネス書です．いずれもタイトルは異なりますが，主軸としている内容は一緒で「ものつくりは成熟してきた．これから重要なのはそれらを活かすストーリーだ」というストーリーを主軸にする考え方で，ストーリーを考えることで，意味あるものつくりができるでしょう，という話です．研究者の実験計画やデータ解析も全く同じなので，文言を読み替えながら読めば，十分参考になる点は多いと思います．
デザイン思考という言葉を初めて本書で知ったのですが，amazonとかを見ると，何冊かあるようで，知られている言葉なのですね．僕の感じとしては「ビジネスモデル」と呼ばれていたものが陳腐化してきたので，色をつけて「デザイン思考」と焼き直しているように感じます．本書でも後半は「ビジネスモデル」の説明が多々入っています．
本書の利点は，図示が全般的に的確で非常に読みやすいです．雰囲気をつかむには，非常に良い本でした．それに対して，具体例はあまりありません．後半が具体例になっているはずなのですが，私には後出しジャンケンな雰囲気が感じられました．
本書の主題の中に「不確実な状況の中で，ストーリーを作成し，ソレを遂行するためのパーツは，柔軟に組み替える」という事があるとおもうのですが，最も難しいのはストーリーを設定することだと思います．既存の「ものつくり」は，今あるものを発展させる（例えば，CPUという既存の物の速度を速くするとか）事に主眼がおかれ，このような状況では物事を演繹的に考えていけば，次に行うことがわかったわけですが，現代では演繹的に考えても次に行うことがわからず，「何をすべきか」を考えた上で「アブダクション」をしないといけない状況になっている．やはり何をすべきかをどのように見つけるかは難しいなぁと思います．この点は，次に示す本の方が，よりつっこんでかいてあります．
本書は，この点の指摘してある方法や内容は非常に読みやすいので，どのように現代になって問題が変わっているのかを知る人には，読みやすく，良い本だと思いました．
■イシューからはじめよ―知的生産の「シンプルな本質」
本書はビジネス書ではあるのですが，著者が脳神経科学で博士号を有しており，本書の内容も随所に研究の話が出てきます．なので，研究者にもおすすめです．
特に，プロジェクトを進めたり，論文を書いたりする上で，ストーリーを作って，それに必要なパーツを考えるために軸を整理する方法や，定量分析には「比較」「構成」「変化」の切り口があるなど，参考になるknow-howがあるとおもいます．研究をする，あるいは，プロジェクト運営をする上で，いずれも必要なスキルなのに，大学などでは教えられることが無いものですので，このような本は参考になります．
本書にあったほうがよかったなぁと思うものに，索引があります．本書の中には有用な図や単語が多いと思うのですが，索引が無いのは後から辞書的に使うときに使いにくいです．その分目次がしっかりしているとも言えるのですが，索引が欲しかったです．
また，全般的に，十分わかりやすいのですが，前述したビジネスのためのデザイン思考の後に読んでしまうと，全般的に無機的で冷たい感じを受けました．両方の本の良いところが合わさると，最高の本なのになぁと思います．
2011-01-26
■[book]Rによるバイオインフォマティクスデータ解析（第２版）
著者の樋口さんより献本頂きました．ありがとうございます．
この本が合う人は明確です．
- Rで何が出来るか眺めたい．
- Rを使い始めて，使い方を検索エンジンで検索したが，全然見つからない（入れるべき単語が分からない）
- RJpWikiとか眺めたけど，専門用語っぽいものがわからないので，俯瞰したい．
- 理論はどうでもよいので，Rで出来そうな事を一通り見てみたい
というタイプの人です．これ一冊でRが1から10まで分かる事はありませんし，データ解析のイロハが分かる事もありません．ましてや，バイオインフォマティクスが何かにも答えてはくれません．例題も時々DNAの塩基配列やBioconductor(Rを用いたバイオデータ解析用の一連のコマンド群)が出てきますが，例題として出ている感じです．
なので，はじめから最後まで通して読むのではなく，辞書的に使う本です．特に図を眺めて，自分の目的にあった検索用語を調べるための簡易辞書として，使えるのではないでしょうか．
2010-11-10
■[hardware]8万円以内で8TBのサーバを作る
- 材料
- HP ProLiant MicroServer (\35,700) http://h50146.www5.hp.com/products/servers/proliant/micro/
- 2TBのHDD * 4本 (4万円前後): 今回は，日立製 0S02602 Deskstar
- バルク品を買うときには，HDDを留めるネジが別途必要（HPの本体に付いてこない）なので注意
- 設置
- 箱明けて，本体取り出します．正面の扉を開ける鍵は，背面に着いてます。写真左の銀色は3.5インチHDD
- 開けて見ます．HDDは，引き出して，マウントする事が出来ます．
- petitベンチマーク（11/21追記）
- よかったところ
- 気になるところ
- 非常に素晴らしい参考資料：http://blog.keshi.org/hogememo/2010/09/24/hp-proliant-microserver-setup-1
2010-11-02
■[bioinformatics]FASTAファイルをGREPするスクリプト
FASTAファイルをアノテーションや配列でgrepする事は度々あるのですが，シェルのgrepだと合致した行だけが出てきて，そのアノテーションを含む配列が何であったか，あるいは，その配列を含む配列の注釈は何かを調べるには一手間必要です．以下のスクリプトは，FASTAを検索して，該当の（行ではなく）エントリを返すスクリプトになります．
- 最後のエントリが表示されないバグ，修正しました．@hinaichigoさん，ありがとうございます！！！
- 使い方
引数: FASTAフォーマットファイル
オプション(特に指定しなければ，すべての電車がおくれます):
-a: 検索範囲をannotationだけに絞ります -s: 検索範囲を配列だけに絞ります．
例：
% ruby grep_fasta.rb -a ATNEK5 TAIR9_cds_20090619 >AT3G20860.1 | Symbols: ATNEK5 | ATNEK5 (NIMA-RELATED KINASE5); ATP binding / kinase/ protein kinase/ protein serine/threonine kinase/ protein tyrosine kinase | chr3:7306147-7308434 FORWARD ATGGCAAACAAGATAAGTGAAACTGCTTCCAAGATGGACGATTATGAGGTCGTGGAACAG ATCGGACGAGGCGCTTTTGGTTCTGCTTTTCTTGTAATTCACAAGTCCGAGAGACGAAAG TATGTGGTTAAGAAGATTCGTTTGGCTAAACAAACAGAGAGATGCAAGCTTGCTGCAATC
#!/usr/bin/env ruby require 'optparse' is_annot_search = true is_seq_search = true opt = OptionParser.new opt.on("-a"){|v| is_seq_search = false} opt.on("-s"){|v| is_annot_search = true} opt.parse!(ARGV) exp = ARGV.shift file = ARGV.shift fp = open(file, "r") name = "" seq = "" cat_seq = "" prev_name = "" while line = fp.gets if line =~ />(.+)$/ || fp.eof if fp.eof seq += line cat_seq += line.chomp end name = $1 if (is_annot_search == true && prev_name =~ /#{exp}/) || (is_seq_search == true && cat_seq =~ /#{exp}/) puts ">#{prev_name}" puts seq end prev_name = name seq = "" cat_seq = "" else seq += line cat_seq += line.chomp end end fp.close()
2010-09-02
■[server]並列分散ファイルシステムLustreを試してみた
- 速度が出なかったのは，100Mのハブに何故か繋がってたからでした．ということで速度書き直しました．
1台では容量が足りないという場合や，速度が足りないという場合用の並列ファイルシステム．最近クラスタやスパコンで使われるようになったlustreをインストールしてみました．既に2.0が出ていますが，ここでは多少ノウハウが溜まっている少し古い1.8の方をインストールしました．
- 準備
- インストール
- CentOS 5.5 (x86_64)を，サーバ設定でインストール
- NISを設定します．「CentOS サーバ設定のメモ」さんのページに良くまとまっているので，この通りに行っています．http://tmcosmos.org/linux/centos/co5_server.html
- Lustreのページ http://wiki.lustre.org/index.php/Main_Page から，Downloadに進んで，RedHat Enterprise の x86_64用に進み，以下のものを持ってきます．
- 全てのホストで，持ってきたrpmをインストールします(下のコマンド一覧参照)．
- 全てのホストの/etc/modprobe.confに"options lnet networks=tcp0(eth0)"を追加します
- modprobe lnet をします(再起動してもOK)
- MGSサーバを作ります．mkfs.lustreコマンドがあるので，それを利用します（下のコマンド一覧参照）．fsname は後の設定で使うので，分かりやすいものを．
- OSSサーバを作ります．mkfs.lustreで作成します．fsnameは参照するMGSサーバのものを作成．
- 必要なら，/etc/fstab を書き換えて自動マウントするようにします．OSSが起動するとき，MGSが無いとfailメッセージを出して終了してしまうので注意．
- クライアントからマウント．普通のファイルシステムと同じ感じで．（下記のコマンド参照）
まとめてコマンド
all% rpm -ivh kernel-2.6.18-194.3.1.el5_lustre.1.8.4.x86_64.rpm lustre-modules-1.8.4-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm lustre-ldiskfs-3.1.3-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm
all% rpm -ivh lustre-1.8.4-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm
all% rpm -Uvh e2fsprogs-1.41.10.sun2-0redhat.rhel5.x86_64.rpm
mgs% /usr/sbin/mkfs.lustre --fsname=lustre --mgs --mdt /dev/sdb1
oss% /usr/sbin/mkfs.lustre --ost --fsname=lustre --mgsnode="oss1" /dev/sda1
client% mount -t lustre oss1@tcp0:/lustre /mnt/
確認．動作を試してみます．今回は，サーバを2台，クライアント1台構成で，MGSに160G, OSS3.2TB(以上2つが同じサーバ"oss1"上)，OSS 2TB("oss2"上)，あとクライアントという3台があります．いずれもOSSは/exportにマウントしました．
- dfしてみる(該当箇所のみ表示)
oss1% df -h
Filesystem Size Used Avail Use% Mounted on
/dev/sdb1 163G 460M 154G 1% /mnt/mgs
/dev/sdb2 3.2T 471M 3.0T 1% /export
oss2% df -h
Filesystem Size Used Avail Use% Mounted on
/dev/sda1 2.0T 471M 1.9T 1% /export
client% df -h
Filesystem Size Used Avail Use% Mounted on
192.168.xxx.xxx@tcp0:/lustre
5.1T 941M 4.9T 1% /mnt
確かに，ほぼOSSノード2個がくっついた容量を示している．
- lustre提供のコマンドでdfしてみる．それぞれのMGSとOSSの容量が分かります．
client% lfs df -h
UUID bytes Used Available Use% Mounted on
UUID bytes Used Available Use% Mounted on
lustre-MDT0000_UUID 163.0G 459.9M 153.2G 0% /mnt[MDT:0]
lustre-OST0000_UUID 3.1T 470.8M 3.0T 0% /mnt[OST:0]
lustre-OST0001_UUID 2.0T 470.2M 1.9T 0% /mnt[OST:1]
filesystem summary: 5.1T 941.0M 4.8T 0% /mnt
- ddで書き込んでみる．
client% dd if=/dev/zero of=/mnt/zero.dat bs=500M count=2
2+0 records in
2+0 records out
NFSだと50MB/s弱なので，2倍以上出る！バックボーンは1Gbpsのetherなので、大体性能限界かな。
お遊びはこれからだ！
注意点
2010-07-28
■[bioinformatics][ngs] BowtieとBWAのインストール
次世代シークエンサのマッピング関係のソフトは成熟度がまだであるせいか、パイプラインの一部にしか成らないのを自覚しているせいか、インストーラが付いていません。なので、普通にコンパイルして、PATHの通っているディレクトリに置きましょう。
- Bowtie, BWA共通
export PATH=$HOME/bin:$PATH
- Bowtie
% unzip bowtie-0.12.5-src.zip % cd bowtie-0.12.5 % make
- 必要なファイルを移動します
% cp bowtie bowtie-build bowtie-inspect $HOME/bin/
- 移動しなくても、./bowtie ... とかすると動きます。また、マニュアルに書かれているサンプルデータは、bowtie-0.12.5のディレクトリのindexesやreadsに入っています。
- BWA
- ダウンロードします。ファイルはhttp://sourceforge.net/projects/bio-bwa/files/ にあります。こちらもCPU毎にわかれていますが、ここではコンパイルもしますので、bwa-0.5.8a.tar.bz2 をダウンロードします。
- 展開してコンパイルします。
% bzip -dc bwa-0.5.8a.tar.bz2 | tar xvf - % cd bwa-0.5.8a % make % cp bwa $HOME/bin/
% cp *.pl $HOME/bin
2010-07-09
■[bioinformatics]膨大なFastaファイルを一定配列数で分割するスクリプト
学生の練習問題的な感じですが，作ったのでメモを残しておきます．
以下の様にすると，input.fasta内の配列を1000本ずつに分けて，input.0001.fasta, input.0002.fasta,... というファイルに保存します．
ruby seq_split.rb input.fasta
オプションは
-o <outdir>: 出力するディレクトリ．outdir下にファイルを出力します（予め作っておいてください）．デフォルトはカレントディレクトリ．
-d <num>: ファイル名に振る番号の桁数を指定．デフォルトは4．
-s <num>: 各ファイルに保存する配列の本数を指定．2000とすると2000本ずつ分けます．デフォルトは1000
なので，
ruby seq_split.rb -o outdir -d 3 -s 2000 input.fasta
の様に使えます．
#!/usr/bin/env ruby require 'optparse' step = 1000 output_dir = "." digit = 4 opt = OptionParser.new opt.on("-d", "--digit [VALUE]", Integer){|v| digit = v.to_i} opt.on("-o", "--output_dir [VALUE]", String){|v| output_dir = v} opt.on("-s", "--step [VALUE]", Integer){|v| step = v.to_i} opt.parse!(ARGV) filename = ARGV[0] prefix = "" suffix = "" if filename =~ /^(.+)\.([^\.]+)$/ prefix = $1 suffix = $2 end fp = open(filename, "r") count = 0 split_num = 0 cur_filename = "" cur_fp = nil while line = fp.gets if line[0,1] == ">" if count % step == 0 STDERR.puts cur_filename cur_fp.close unless cur_fp.nil? split_num += 1 cur_filename = prefix + "." + "%0#{digit}d" % split_num + "." + suffix cur_fp = open(output_dir + "/" + cur_filename, "w") count = 0 end count += 1 end cur_fp.print line end STDERR.puts cur_filename cur_fp.close
