2011-03-17

積ん読の解体

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

2011-01-26

Rによるバイオインフォマティクスデータ解析（第２版）

book

著者の樋口さんより献本頂きました．ありがとうございます．
この本が合う人は明確です．

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

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

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

Rによるバイオインフォマティクスデータ解析第2版 −Bioconductorを用いたゲノムスケールのデータマイニング−

posted with amazlet at 11.01.26

樋口千洋
共立出版
売り上げランキング: 18450

Amazon.co.jp で詳細を見る

2010-11-10

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

hardware

材料
- 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は，引き出して，マウントする事が出来ます．

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

あとは、CentOSいれて、Software raidを組んだら終了！

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 (CPUは20%位、かなぁ)
- ソフトウエアRAID0の6TB。312 MB/s (CPUはDDが80%以上使う時間帯あり、他にディスク周りが20%前後)。

よかったところ
- 安い！安い！安い！
- 非常に静か（冷蔵庫より静か）
- 癖が無く、あっさりインストールが済む
- コンパクトなので、設置場所を選ばない
- 全面にUSBがある（結構重要）
気になるところ
- メモリが1GBは少ない
- CPUはやっぱり非力なのでsoftware raid+少しの作業で手一杯。でも用途を選べば十分。
  - 同時アクセス数の少ないWebサーバとか
- lustreのノードにしたくて作ってみました。実験結果はそのうち。
非常に素晴らしい参考資料：http://blog.keshi.org/hogememo/2010/09/24/hp-proliant-microserver-setup-1

2010-11-02

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

bioinformatics

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

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

server

速度が出なかったのは，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のインストール

bioinformatics ngs

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

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

export PATH=$HOME/bin:$PATH

- - ログアウトしてログインしなおすか、source ~/.bashrc (or source ~/.zshrc)をします。
Bowtie
- ダウンロードします。ファイルは、http://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12.5 にあります。CPU毎にわかれていますが、ここではコンパイルもする為、bowtie-0.12.5-src.zip　をダウンロードします。
- 展開してコンパイルします

% 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/

- いくつかのperl スクリプトが入っているので、もし使いたければコピーします。

% cp *.pl $HOME/bin

2010-07-09

膨大なFastaファイルを一定配列数で分割するスクリプト

bioinformatics

学生の練習問題的な感じですが，作ったのでメモを残しておきます．

以下の様にすると，input.fasta内の配列を1000本ずつに分けて，input.0001.fasta, input.0002.fasta,... というファイルに保存します．

ruby seq_split.rb input.fasta

オプションは

-o : 出力するディレクトリ．outdir下にファイルを出力します（予め作っておいてください）．デフォルトはカレントディレクトリ．
-d : ファイル名に振る番号の桁数を指定．デフォルトは4．
-s : 各ファイルに保存する配列の本数を指定．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