Hatena::ブログ(Diary)

Loud Minority このページをアンテナに追加 RSSフィード

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
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のインストール
    • (注1) which easy_install をして,"/usr/local/share/python/easy_install" である事を確認する."/usr/bin/easy_install" を指している場合には,shell profileを読みなおすか,shell を再起動
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
    • qtのインストールがスタートするので,非常〜に時間がかかります.諦めずに放置して寝ましょう.
brew install zmq
pip install pyzmq
pip install pygments
  • 起動!
ipython qtconsole --pylab=inline

こんな感じになります.

f:id:sesejun:20110828131405p:image

2011-08-20

[]RとPython(NumPy)の対応

> 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)]
>>> 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
  • 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] 
    • stringID一覧などで絞るとき用.もう少し簡単にできないものか.
>tf = [x in ["M","F"] for x in ary['gender']]
>ary[np.array(tf)]

2011-03-30

[]お茶大から東工大に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です.寄付金もお待ちしております:)

*1:お茶大は復職条件付きでの退職

*2:大学の教員はあまり上下関係は無く(特にお茶大の情報は上下関係がほぼ無く),私のような若造から定年退官される直前の先生までフラットに扱われますので,失礼を承知で同僚という言葉を使わせていただいております

*3:お茶大情報では,研究室配属は学生達が話しあって決めます.教官は受け入れ人数の上限を指定するだけです

2011-03-17 積ん読の解体

旅行地震の影響で,飛行機が飛ばずキャンセルになったので,積ん読していた本を順に読んでいました.感想文です.4冊とも面白く読めました.

闘う! ウイルスバスターズ 最先端医学からの挑戦 (河岡義裕,渡辺登喜子)

ウイルスで有名な河岡先生と渡辺先生の本.今,研究者に進むか,会社に務めるか迷っている人(特に生命系)が居たら,読んでみることをおすすめします.研究者が必ずしも小さい頃から「その道」を選んでいるわけではなく,環境によって,色々進路が分かれてきた様子が見て取れます.特にII部の3章のネイチャーに載るまでの紆余曲折や,II部4章の研究者による研究者の紹介(人生談)は楽しかったです(この章に限らず,渡辺先生の章は,旦那さんとの会話など生活感があって,和みます.女性で研究者を考えている人は参考になる所があるかもしれません).

「おわりに」の所で,河岡先生が実際に執筆した人の名前を出しているのは,多くの本がゴーストライター(と呼んで良いのか分かりませんが)を表に出さない中,河岡先生の真摯な態度が出ていて,好感がもてました.

こころ」は遺伝子でどこまで決まるのか―パーソナルゲノム時代の脳科学 (宮川 剛)

ゲノムと脳の関係は,私自身の興味でもあり,楽しく読めました.

また,著者の宮川先生研究で作成されている「網羅的行動テストバッテリー」は,私が大学院生だったか,助手時代だったか忘れましたが,参加させていただいた科学研究費の会議で,その話を拝見し「これはすごい!」と感動したのを覚えています.多くの発表が遺伝子ひとつに着目する中,「データを網羅的に取る」という意味重要性を本当に理解している研究者なのだと,その時思いました.

本書の中は,ゲノムと行動の関係について,入門からスタートし,体験談を含めつつ,分かっている事を順に説明しています.前半が入門,中盤は23andmeによる検査の話を挟みつつ,後半はSNPs(一塩基置換)と行動に関して,分かっている事を列挙している形です.これから,まだまだ分かることは増えていくでしょうが,現時点でのパーソナルゲノム時代に向けて分かっている事を整理してある本として,良書であると思います.

ビジネスのためのデザイン思考

この本と次の本はビジネス書です.いずれもタイトルは異なりますが,主軸としている内容は一緒で「ものつくりは成熟してきた.これから重要なのはそれらを活かすストーリーだ」というストーリーを主軸にする考え方で,ストーリーを考えることで,意味あるものつくりができるでしょう,という話です.研究者実験計画やデータ解析も全く同じなので,文言を読み替えながら読めば,十分参考になる点は多いと思います.

デザイン思考という言葉を初めて本書で知ったのですが,amazonとかを見ると,何冊かあるようで,知られている言葉なのですね.僕の感じとしては「ビジネスモデル」と呼ばれていたものが陳腐化してきたので,色をつけて「デザイン思考」と焼き直しているように感じます.本書でも後半は「ビジネスモデル」の説明が多々入っています.

本書の利点は,図示が全般的に的確で非常に読みやすいです.雰囲気をつかむには,非常に良い本でした.それに対して,具体例はあまりありません.後半が具体例になっているはずなのですが,私には後出しジャンケンな雰囲気が感じられました.

本書の主題の中に「不確実な状況の中で,ストーリーを作成し,ソレを遂行するためのパーツは,柔軟に組み替える」という事があるとおもうのですが,最も難しいのはストーリーを設定することだと思います.既存の「ものつくり」は,今あるものを発展させる(例えば,CPUという既存の物の速度を速くするとか)事に主眼がおかれ,このような状況では物事を演繹的に考えていけば,次に行うことがわかったわけですが,現代では演繹的に考えても次に行うことがわからず,「何をすべきか」を考えた上で「アブダクション」をしないといけない状況になっている.やはり何をすべきかをどのように見つけるかは難しいなぁと思います.この点は,次に示す本の方が,よりつっこんでかいてあります.

本書は,この点の指摘してある方法や内容は非常に読みやすいので,どのように現代になって問題が変わっているのかを知る人には,読みやすく,良い本だと思いました.

イシューからはじめよ―知的生産の「シンプル本質

本書はビジネス書ではあるのですが,著者が脳神経科学博士号を有しており,本書の内容も随所に研究の話が出てきます.なので,研究者にもおすすめです.

特に,プロジェクトを進めたり,論文を書いたりする上で,ストーリーを作って,それに必要なパーツを考えるために軸を整理する方法や,定量分析には「比較」「構成」「変化」の切り口があるなど,参考になるknow-howがあるとおもいます.研究をする,あるいは,プロジェクト運営をする上で,いずれも必要なスキルなのに,大学などでは教えられることが無いものですので,このような本は参考になります.

本書にあったほうがよかったなぁと思うものに,索引があります.本書の中には有用な図や単語が多いと思うのですが,索引が無いのは後から辞書的に使うときに使いにくいです.その分目次がしっかりしているとも言えるのですが,索引が欲しかったです.

また,全般的に,十分わかりやすいのですが,前述したビジネスのためのデザイン思考の後に読んでしまうと,全般的に無機的で冷たい感じを受けました.両方の本の良いところが合わさると,最高の本なのになぁと思います.

2011-01-26

[]Rによるバイオインフォマティクスデータ解析(第2版)

著者の樋口さんより献本頂きました.ありがとうございます

この本が合う人は明確です.

  • Rで何が出来るか眺めたい.
  • Rを使い始めて,使い方を検索エンジン検索したが,全然見つからない(入れるべき単語が分からない)
  • RJpWikiとか眺めたけど,専門用語っぽいものがわからないので,俯瞰したい.
  • 理論はどうでもよいので,Rで出来そうな事を一通り見てみたい

というタイプの人です.これ一冊でRが1から10まで分かる事はありませんし,データ解析のイロハが分かる事もありません.ましてや,バイオインフォマティクスが何かにも答えてはくれません.例題も時々DNA塩基配列Bioconductor(Rを用いたバイオデータ解析用の一連のコマンド群)が出てきますが,例題として出ている感じです.

なので,はじめから最後まで通して読むのではなく,辞書的に使う本です.特に図を眺めて,自分目的にあった検索用語を調べるための簡易辞書として,使えるのではないでしょうか.

2010-11-10

[]8万円以内で8TBのサーバを作る

f:id:sesejun:20101110120531j:image

  • 開けて見ます.HDDは,引き出して,マウントする事が出来ます.

f:id:sesejun:20101110120719j:image

    • HDD以外(メモリ等)を増設するには,本体下部のマザーボードを引き出します.トルクスネジですが,工具は蓋に着いてます.
    • ケーブルが押し込めてあるので,結構取り出す&しまうのに苦労します.

f:id:sesejun:20101110122132j:image

  • petitベンチマーク(11/21追記)
    • 運用では、デフォルトで付いてくる160GBのディスクを / に、それ以外に入れた2TB (計6TB)を Software RAID0 で運用はじめました。また、メモリは増設して2GBにしてあります。
    • dd if=/dev/zero of=test bs=1G count=5 でテスト
    • RAIDしてない160GB のディスク。88.3 MB/s (CPU20%位、かなぁ)
    • ソフトウエアRAID0の6TB。312 MB/s (CPUはDDが80%以上使う時間帯あり、他にディスク周りが20%前後)。

2010-11-02

[]FASTAファイルGREPするスクリプト

FASTAファイルをアノテーション配列でgrepする事は度々あるのですが,シェルのgrepだと合致した行だけが出てきて,そのアノテーションを含む配列が何であったか,あるいは,その配列を含む配列の注釈は何かを調べるには一手間必要です.以下のスクリプトは,FASTAを検索して,該当の(行ではなく)エントリを返すスクリプトになります.

  • 使い方

引数: 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

[]並列分散ファイルシステムLustreを試してみた

  • 速度が出なかったのは,100Mのハブに何故か繋がってたからでした.ということで速度書き直しました.

1台では容量が足りないという場合や,速度が足りないという場合用の並列ファイルシステム最近クラスタスパコンで使われるようになったlustreをインストールしてみました.既に2.0が出ていますが,ここでは多少ノウハウが溜まっている少し古い1.8の方をインストールしました.

  • 準備
    • マシンを3台以上用意します.
      • 1台はMGS(/MDT)と呼ばれるメタデータを入れる所.ファイルシステム全体の1-2%位を用意.ディスクの追加はできるが,MGSの追加は出来ない(or 厄介?)ようなので,少し多めに2%を取りましょう.
      • 1台はOSSと呼ばれる実際にデータ入れる所.MGSと兼ねることができます.
      • 1台はクライアント(MGSやOSSと兼ねる事は推奨しません)
      • 今回はノード間はGigabit etherで繋がっています.
  • インストール
    1. CentOS 5.5 (x86_64)を,サーバ設定でインストール
    2. NISを設定します.「CentOS サーバ設定のメモ」さんのページに良くまとまっているので,この通りに行っています.http://tmcosmos.org/linux/centos/co5_server.html
    3. Lustreのページ http://wiki.lustre.org/index.php/Main_Page から,Downloadに進んで,RedHat Enterprise の x86_64用に進み,以下のものを持ってきます.
      1. kernel-2.6.18-194.3.1.el5_lustre.1.8.4.x86_64.rpm
      2. lustre-modules-1.8.4-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm
      3. lustre-ldiskfs-3.1.3-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm
      4. lustre-1.8.4-2.6.18_194.3.1.el5_lustre.1.8.4.x86_64.rpm
      5. e2fsprogs-1.41.10.sun2-0redhat.rhel5.x86_64.rpm
    4. 全てのホストで,持ってきたrpmをインストールします(下のコマンド一覧参照).
    5. 全てのホストの/etc/modprobe.confに"options lnet networks=tcp0(eth0)"を追加します
    6. modprobe lnet をします(再起動してもOK)
    7. MGSサーバを作ります.mkfs.lustreコマンドがあるので,それを利用します(下のコマンド一覧参照).fsname は後の設定で使うので,分かりやすいものを.
    8. OSSサーバを作ります.mkfs.lustreで作成します.fsnameは参照するMGSサーバのものを作成.
    9. 必要なら,/etc/fstab を書き換えて自動マウントするようにします.OSSが起動するとき,MGSが無いとfailメッセージを出して終了してしまうので注意.
    10. クライアントからマウント.普通のファイルシステムと同じ感じで.(下記のコマンド参照)

まとめてコマンド

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

1048576000 bytes (1.0 GB) copied, 8.75246 seconds, 120 MB/s

NFSだと50MB/s弱なので,2倍以上出る!バックボーンは1Gbpsのetherなので、大体性能限界かな。

お遊びはこれからだ!

注意点

  • クライアントからマウントした状態で,OSSを切り離すとクライアントからファイルが見えなくなるだけでなく,コマンドが固まります.なので,切り離す時は注意.(OSSとクライアントを同一のPCで行うときも注意)

2010-07-28

[][] BowtieとBWAのインストール

次世代シークエンサのマッピング関係のソフト成熟度がまだであるせいか、パイプラインの一部にしか成らないのを自覚しているせいか、インストーラが付いていません。なので、普通にコンパイルして、PATHの通っているディレクトリに置きましょう。

  • Bowtie, BWA共通
    • パスを通します。ここではbash/zsh系を使っているとして、$HOME/bin にインストールすると仮定します。
      • .bashrc もしくは .profile (zshなら.zshrc)に以下を付け加えます。
export PATH=$HOME/bin:$PATH
% 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

[]膨大な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