Biopythonで塩基配列をアミノ酸配列にする

1
2
3
4

>>> from Bio.Seq import Seq
>>> mRNA = Seq("TATGAAAGT")
>>> mRNA.translate()
Seq('YES', ExtendedIUPACProtein())

Why Biopython?

xkcd: Python

http://xkcd.com/353/

Ref.



2013 Web Server Issue

http://nar.oxfordjournals.org/content/41/W1.cover-expansion

今年のWeb Server issue 2013が既に閲覧できるようになっていた．

Nucleic Acid Research誌は少し特殊な特別号を毎年出していて，1月にはDatabase issueが，7月にはWeb Server issueが発刊されるのが定番になっている．これらの特集号は，他の論文誌では論文になりにくいようなデータベース構築やウェブサービス系を論文として発表する場として，バイオインフォ系にとっては重要な投稿先の一つといえる．特に，ウェブサービスのアップデートなども同様の投稿ができるようになっており，新規性としては少し弱いウェブサービスの保守管理や機能追加に関しても評価されるのが大きな特徴だ（ただし2年間のインターバルが必要）．

今年のWeb Server issue：Table of Contents — 1 July 2013, 41 (W1)

今年は合計で95本のウェブサービスに関する論文が採択されたようだ．Editorialによると

For the 2013 issue, 293 summaries were submitted and 129, or 44%, were approved for manuscript submission. Of those approved, 95, or 73%, were ultimately accepted for publication.

http://nar.oxfordjournals.org/content/41/W1/W1.full

ということで，293本の投稿のうち最終的に95本が採択され，全体で見れば約32%の採択率となっている．


来年のWeb Server Issue 2014の締め切りは今年いっぱいまで(31 Dec. 2013)となっている．それまでに完全に動くウェブサービスと1ページのproposalを書いて投稿し，そこで通れば1ヶ月でmanuscriptを書いて7月に掲載という流れのようだ．



Tophatにおいて使用するスレッド数を上げすぎるとファイルディスクリプタの上限でエラーになる

久しぶりにTophatをいじったら一発目でエラーで落ちてしまい，出鼻をくじかれてしまった…．原因を調べてみると，どうやらTophatのコマンド実行時に指定する-pオプションの数字を調子に乗って上げすぎたのが良くなかったらしい．といっても計算機のリソースに問題があったわけではなく，ユーザが使えるシステムリソースの制限を超えてしまったために，Tophatが停止してしまったようだ．この問題はTophatの公式サイトのFAQでも取り上げられているが，このBlogでもエラーメッセージとともに原因と対策を書いておこうと思う．

Tophatの標準エラー出力のメッセージ

今回の場合は，Tophatのログの最後で以下のようなエラーメッセージが表示される．

1
2
3
4
5
6
7

[...]
[2013-06-23 00:22:33] Mapping right_kept_reads_seg3 to genome chlamydomonas_236 with Bowtie2 (3/4)
[2013-06-23 00:33:09] Mapping right_kept_reads_seg4 to genome chlamydomonas_236 with Bowtie2 (4/4)
[2013-06-23 00:44:57] Searching for junctions via segment mapping
        [FAILED]
Error: segment-based junction search failed with err =1
Error opening SAM file output/tmp/right_kept_reads_seg1.bam

これだけでは何のことだかよくわからないが，Tophatが出力するログのlogs/segment_juncs.logを見ると，このエラーに関してもう少し詳細な記述が出力されている．

1
2
3
4
5
6
7
8
9

    Loading scaffold_53...done
    Loading scaffold_54...done
    Loading ...done
>> Performing segment-search:
Loading left segment hits...
done.
Loading right segment hits...
open: Too many open files
Error opening SAM file output/tmp/right_kept_reads_seg2.bam

これを見ると，どうやらsegment_juncsの実行中に「open: Too many open files」が原因で実行が落ちたようだ．

原因と対策

このエラーに関しては，実はTophatの公式サイトのFAQに「What should I do if I see a message like “Too many open files”?」という，まさに先ほどのエラーメッセージの内容そのままの項がある．

This usually happens when using “-p” option with a large value (many threads). TopHat may produce many intermediate files, the number of which is proportional to this value; sometimes the number of the files may go over the maximum number of files a process is allowed to open. The solution is to raise the limit to a higher number (e.g. 10000). For Mac, you can change this using a command, “sudo sysctl -w kern.maxfiles=10240”.

TopHat :: Center for Bioinformatics and Computational Biology

ざっと要約すると「Tophatは中間ファイルを大量に作るから時々許容数超えちゃうんだよね．扱えるファイル数の上限上げるか，Macなら次のコマンドで対処してね」ということになる．さすがにこれだけではよく分からないと思うので，SEQanswersの以下のスレッドも参考にしつつ，もう少し詳しく見ていこうと思う．

Tophat segment junction error 1, invalid BAM binary header - SEQanswers

limitとファイルディスクリプタ数の制限

さて，公式サイトのFAQにおいてファイル数の上限といった言葉が出てきたように，LinuxやMacではユーザごとに使用することのできる各種システムリソースに制限が設けられている．具体的にはユーザが使用できるCPU数やメモリの容量，プロセスの数などがそれに該当するのだが，その中に「ファイルディスクリプタ数」という項目があり，1プロセスが同時に開くことのできるファイル数の上限を定めている．これは主にひとつの計算機を複数人が使用するマルチユーザシステムにおいて，一人がリソースを独占するのを防いだりアプリケーションの暴走を止めたりするのに役立つのだが，今回のような大規模な計算を実行する際にはこれが邪魔になってしまう．

システムリソースの制限を確認したい場合には，「ulimit -a」または「limit」というコマンドを使う．例えば，私の環境では以下のようになる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# ソフトリミット
$ ulimit -a
-t: cpu time (seconds)              unlimited
-f: file size (blocks)              unlimited
-d: data seg size (kbytes)          unlimited
-s: stack size (kbytes)             8192
-c: core file size (blocks)         0
-m: resident set size (kbytes)      unlimited
-u: processes                       1024
-n: file descriptors                1024
-l: locked-in-memory size (kbytes)  64
-v: address space (kbytes)          unlimited
-x: file locks                      unlimited
-i: pending signals                 127413
-q: bytes in POSIX msg queues       819200
-e: max nice                        0
-r: max rt priority                 0

ここで表示されるのは厳密にはソフトリミットと呼ばれ，ユーザごとに設定されている制限である．一方でハードリミットと呼ばれる制限もあり，これは管理者(root)が定める制限となっている．ソフトリミットは，このハードリミットの範囲内でしか自由に制限値を変更することはできない．ハードリミットを確認したい場合には，-Hオプションを付ける．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# ハードリミット
$ ulimit -aH
-t: cpu time (seconds)              unlimited
-f: file size (blocks)              unlimited
-d: data seg size (kbytes)          unlimited
-s: stack size (kbytes)             unlimited
-c: core file size (blocks)         unlimited
-m: resident set size (kbytes)      unlimited
-u: processes                       127413
-n: file descriptors                4096
-l: locked-in-memory size (kbytes)  64
-v: address space (kbytes)          unlimited
-x: file locks                      unlimited
-i: pending signals                 127413
-q: bytes in POSIX msg queues       819200
-e: max nice                        0
-r: max rt priority                 0

さて，今回問題になっているファイルディスクリプタ数は「-n: file descriptors」で表示されている．上の例の場合，ソフトリミットでは1024，ハードリミットでは4096となっている．つまり，Tophatは1024個以上のファイルを1スレッドで開こうとしたために，このソフトリミットに引っかかってしまったようだ．

フィアルディスクリプタ数の制限への対策

この問題を回避するには主に2つの方法がある．

1．Tophatの-pオプションの値を小さくする

2．ソフトリミットのファイルディスクリプタの値を大きくする

まず1.では，Tophatが使用するスレッド数を少なくすることで，ファイルディスクリプタ数の上限に引っかからなくするというもの．一度-pオプションを無くして実行してみれば，おそらく今回のエラーには引っかからなくなるだろう．一度に開くファイル数が少なくなりほぼ確実に実行できるようにはなるが，並列処理数が減ってしまうのでTophatの実行時間は長くなってしまう．

そこで2.のようにの制限を無くして，-pオプションはそのままにファイルディスクリプタの上限を回避するという方法もある．実行コマンドや実行時間はそのままにエラーを回避することができる一方で，上限を上げたからといってもTophatがそれ以上の同時ファイルオープンをしてしまえば同様のエラーに引っかかってしまうほか，制限を上げたことにより計算機に負荷がかかる恐れもある．つまり，時と場合によっては成功するが確証は無いという感じだろうか．

ちなみに，私の場合はファイルディスクリプタ数を上げても以下のような別のエラーが出て実行できなかった．

1

Error: ReadStream::getRead() called with out-of-order id#!

ということで，結論としてはTophatの実行時間との兼ね合いを考えて，どちらかを選択したほうが良いだろう．素直に-pオプションの値を下げるほうが無難な気がする．

ちなみに，制限値を引き上げるには，ulimitで以下のように値を変更する．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

$ ulimit -n 2000
$ limit
-t: cpu time (seconds)              unlimited
-f: file size (blocks)              unlimited
-d: data seg size (kbytes)          unlimited
-s: stack size (kbytes)             8192
-c: core file size (blocks)         0
-m: resident set size (kbytes)      unlimited
-u: processes                       1024
-n: file descriptors                2000
-l: locked-in-memory size (kbytes)  unlimited
-v: address space (kbytes)          unlimited
-x: file locks                      unlimited
-i: pending signals                 16545839
-q: bytes in POSIX msg queues       819200
-e: max nice                        0
-r: max rt priority                 0

ulimitの後ろに該当するパラメータの値を指定することによって，上限を引き上げることができる．ただし，先ほど述べたようにハードリミットより上は指定することができないので注意が必要になる．

まとめ

Tophatのエラー「open: Too many open files」は，スレッドが一度に開くことのできるファイルディスクリプタ数がソフトリミットの上限に引っかかってしまったために起こる．-pオプションの値を下げて実行するか，ソフトリミットのファイルディスクリプタの上限を引き上げることによって回避することができる．まずは，-pオプションを指定せずに実行してみよう．

参考

実行環境

  • OS：RHEL 6.3
  • Tophat：v2.0.8b
  • Bowtie：version 2.1.0
  • Samtools：0.1.19.0


SRA Toolkitを使ってsraファイルからfastqファイルに変換する(更新版)

久しぶりにfastq-dumpを使ったら少し仕様が変わっていたのでメモ．だいぶ前に別の場所で書いた記事は，既に古くなってしまっている．何もせず記事を放置しておくのも申し訳ないし，どうやらfastq-dumpでググると上位にくるようなので，ここいらで更新しておかないと…．

(deprecated) Wolf Ears » SRA Toolkitを使ってsraファイルからfastqファイルに変換する

SRA Toolkitのインストール

sraからfastqに変換するプログラムは，SRA Toolkitの中にあるfastq-dumpを使う．まずはNCBIのサイトからOS環境に合わせてコンパイルされたSRA Toolkitをダウンロードしてインストールする．

Software : Software : Sequence Read Archive : NCBI/NLM/NIH

コンパイルされたバイナリとしては，Linux系列のCentOSやUbuntu，MacOSやMS Windowsに合わせたものが用意されている．私の環境はRed Hat Enterpriseなので，とりあえず「CentOS Linux 64 bit architecture」を選択する．

1
2

$ wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.3.2-5/sratoolkit.2.3.2-5-centos_linux64.tar.gz
$ tar zxvf sratoolkit.2.3.2-5-centos_linux64.tar.gz

今回使用するfastq-dumpは展開したディレクトリのbin以下に入っている．以下のように実行して動作するか確認する．

1
2
3
4
5
6
7
8
9

$ sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump

Usage:
  sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump [options] <path [path...]>
  sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump [options] [ -A ] <accession>

Use option --help for more information

sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump : 2.3.2

fastq-dumpの使い方

基本的な使い方としては，以下のようにfastq-dumpにsraファイルを指定すればよい…

1

$ fastq-dump ./DRR002191.sra

…のだが，ここでひとつ注意が必要になる．上のコマンドを指定すると，データがシングルエンドでもペアエンドでも同様に一つのファイルとして出力されてしまう．実際に変換されたfastqファイルの中身を見てみると，上のDRR002191.sraは本当は90bpのペアエンドなのだが，以下のように結合された配列として出力される．

1
2
3
4
5

$ head -n 4 DRR002191.fastq
@DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=180
NCAATGGCAATAGCAATGCATTGAAATGAAAAGGCATTTACCAGGAGCAGGAAAGCCAGAAAGAGGAGCAGTGNNCNNGGAGTTGCNNNNTGCTTGGCTTCATCTTTTGCAATTCTGCATCTTTTGCATTTCCCTTCTCGCTCTGCTGCTCGCTCGCCTTTCTTNNNCGNCCTTTTCCCC
+DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=180
#1=DDFFFHFFHHIIJJJIJJIJJIJBCHHIHGCC>FGIGIHDH?HJIEGEGI@FEHIJIIE@HI9=EH=CDD#################B@CFFDDDDFHHHGIJJIJJGGIGJIIHGHIGIJIHIJIEIIHIJJIGGHHEGI@DGCGEFHHIGIEFEFFBCB################

ということで，きちんとペアエンドを2つのfastqファイルに変換するには，–split-filesというオプションを使用する．

1

$ fastq-dump --split-files ./DRR002191.sra

これで，きちんと2つのfastqファイルが出力される．個別のfastqファイルを見ると，やはり先程の変換はペアエンドの配列を結合して出力されていることがわかる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

$ ls *.fastq
DRR002191_1.fastq  DRR002191_2.fastq

$ head -n 4 DRR002191_1.fastq DRR002191_2.fastq
==> DRR002191_1.fastq <==
@DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=90
NCAATGGCAATAGCAATGCATTGAAATGAAAAGGCATTTACCAGGAGCAGGAAAGCCAGAAAGAGGAGCAGTGNNCNNGGAGTTGCNNNN
+DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=90
#1=DDFFFHFFHHIIJJJIJJIJJIJBCHHIHGCC>FGIGIHDH?HJIEGEGI@FEHIJIIE@HI9=EH=CDD#################

==> DRR002191_2.fastq <==
@DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=90
TGCTTGGCTTCATCTTTTGCAATTCTGCATCTTTTGCATTTCCCTTCTCGCTCTGCTGCTCGCTCGCCTTTCTTNNNCGNCCTTTTCCCC
+DRR002191.1 FCD0BMAACXX:5:1101:1133:1889# length=90
B@CFFDDDDFHHHGIJJIJJGGIGJIIHGHIGIJIHIJIEIIHIJJIGGHHEGI@DGCGEFHHIGIEFEFFBCB################

その他の注意点

fastq-dumpに指定するsraファイルは，きちんとパスを指定しないといけないようだ．以下のようなエラーが出る場合は，絶対パスを指定するか，カレントディレクトリにsraファイルがある場合にはファイル名の前に「./」を付け加える．

パスが解決できなくて実行できない例

1
2
3

$ fastq-dump DRR002191.sra
2013-06-20T02:11:42 fastq-dump.2.3.2 err: name not found while resolving tree within virtual file system module - failed to open 'DRR002191.sra'
Written 0 spots total

絶対パスを指定するか./を付け加える

1
2

$ fastq-dump /path/to/DRR002191.sra
$ fastq-dump ./DRR002191.sra

また，他のオプションなどの情報はfastq-dumpのhelpから見ることができる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14

$ fastq-dump --help

Usage:
  /home/is/yuki-ok/biolocal/src/sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump [options] <path [path...]>
  /home/is/yuki-ok/biolocal/src/sratoolkit.2.3.2-5-centos_linux64/bin/fastq-dump [options] [ -A ] <accession>

INPUT
  -A|--accession <accession>       Replaces accession derived from <path> in
                                   filename(s) and deflines (only for single
                                   table dump)
  --table <table-name>             Table name within cSRA object, default is
                                   "SEQUENCE"

[...]

環境

  • NCBI SRA Toolkit : May 9 2013, version 2.3.2-5 release
  • fastq-dump : 2.3.2

参考



MacPortsで使用していない過去のバージョンのパッケージを消す＆ちょっと可視化

MacBook Airの128GBしかない非力なSSDもそろそろ満杯になってきた．最近ではスワップが溜まっては再起動するということを繰り返していたので，色々と容量の食っているファイルをちまちま消していた．そういった地道な努力をしつつ無駄なキャッシュや不必要なファイルなどをちゃんと調べてみると，どうやらMacPortsが10GB近く占めていることがわかり，しかもバージョンアップによって使用されなくなったgccやllvm，boostなどの古いバージョンのせいだとわかったので，思い切って消してみることにした．

古いバージョンのパッケージを消す

MacPortsでは，バージョンが上がってアップデートする際に，インストールされているバージョンのパッケージをInactiveとし最新バージョンをActiveにするという方法で，過去のバージョンのパッケージが残される仕組みになっている．ということで，Inactiveなパッケージを全て消去するために，以下のコマンドを使用する．

1

$ sudo port -u uninstall

これでInactiveなパッケージを全て消去することができる．なお，古いバージョンとして残っていたパッケージを消すということは気軽にバージョンを変更したりできなくなるということなので，もしバージョンアップによる不具合に対応する必要が出てくると後々面倒になることを注意しておきたい．

結果

/opt/localが10GBから5GBになり，かなりのダイエットに成功した！容量の少ないMacBook Airでも，これでもうちょっと耐えられそう．

DaisyDiskによる可視化

さて，ここからは完全に蛇足なのだが，最近使用しているのDaisyDiskという可視化ソフトウェアがなかなか気に入っているので，今回のファイル消去の前後で/opt/localがどれだけ変化するかを確かめてみた．DaisyDiskは本来はインタラクティブに操作することができて，気になる箇所をクリックしたりして拡大することができるのだが，今回は前後の画像比較だけとなっている．使い心地が気になる人は以下のリンクから公式ページに飛んで紹介動画を見ていただきたい．

DaisyDisk - Analyze disk usage and free up disk space on Mac

Before

After

まず全体の容量が10GBから5GBに変わっており，それぞれの面積は単純に割合を示していることに注意．Inactiveなパッケージを消去する前は少数の比較的大きな領域がいくつも見られる．具体的には，左側の外縁から2つ目のピンクの領域がgccやllvm，boost，emacsなどのパッケージのディレクトリを表しており，一番外縁にある灰色がそれぞれのパッケージのファイルになっている．一つのパッケージに複数の灰色で区切られた細かいディレクトリが含まれているのは，それだけアップデートして古いバージョンのパッケージが溜まっていったからだ．そしてInactiveなパッケージを消去した後は，複数の小さな領域が組み合わさっていることがわかり，それぞれのパッケージに含まれるファイルの容量が小さくなったことを示している．

参考