読者です 読者をやめる 読者になる 読者になる

yokaのblog

琵琶湖で微生物の研究してます

学生最後の1年

 今日から新年度、博士課程も残り1年を切った。会社を退職して研究に戻ってきて2年、ここまでは研究の勘を取り戻すことに必死だったけど、だんだんこの生活が当たり前になってきて、良い意味でも、悪い意味でも、緊張がほぐれてきた気がする。

 研究に戻ってくる前から重々分かってはいたことだけど、研究者の将来はとても厳しい。しばらく研究していれば具体的なキャリアパスが見えてきて前向きになれるかとも思ったけど、むしろ世知辛い話を直接耳にする機会が増えて、ますます自分の将来を案ずるようになった。僕には、嫌なことが起こりそうになったら、最悪を想像することで予防線を張って、いざその事態が起こった時の精神的ダメージを軽減するような本能が備わっている。だから、研究に復帰するために無我夢中だったフェーズが終わった今、あまり将来に期待せず、いつ辞めることになっても良いように、淡々と研究に接するようになったと思う。

 もちろん研究は楽しいし、毎日本気でやっている。僕は今の仕事が自分に合っていて、前職を辞めて戻ってきて良かったと心の底から思っている。だけど、全身全霊で研究を愛することはしない。

「自分には研究しかない」と思うことをしなくなった

というほうが分かりやすいかもしれない。研究が続けられなくなったときの精神的ショックが怖いし、研究にしがみつくために不利な選択を強いられることが嫌だからだ。この点は、会社を辞めて「これからは命懸けで研究に身をささげるのだ!」と意気込んでいた2年前とは明らかに考え方が変わった。

 研究に対する情熱が無くなったということではない。むしろ、研究は当初の想定以上に色々なことが分かってきて、ますます面白くなってきている。変わったのは、本気で研究を楽しむ自分の上に、もう一段メタな自分を置くようになった、ということだと思う。

「本気」から「淡々と本気」になった

と言ってもいいかもしれない。冷めた目の自分が、没頭しすぎないように常に監視している。もっと言うと、

別に研究以外でも本気を出そうと思えば出せるし、何ならいつでも辞めてやるぞ

と思いながら研究している。というか、そう思うことで安心を得ようとしている。

 こういう考え方の人間は、業界には歓迎されないだろう。それは、何か業界全体にとって不利な出来事があった時に、皆で一丸になって抵抗しようとせず、手のひらを返して真っ先に逃げ出してしまうからだ。だけど、色々考えてみたけど、元来悲観的な性格の僕は、常に研究を見限る準備をした状態で研究を続けることでしか、この厳しい業界で心の安定を得る方法がないのだと思う。

 もちろん、続けるからには本気でやるし、行けるところまで行きたいと思う。こうやって偉そうなことを書いているけど、今の自分は業績的には研究業界から見たらゴミ以下で、いなくなったって全然困らない存在だ。学生最後の1年、もう少しで良い論文が出せそうだし、次の研究の面白いデータもたくさん出ているし、海外にも長期滞在する予定で、D論を書いて、将来を決めなければならない。長い一年になると思う。引き続き、淡々と本気で頑張りたい。

2月琵琶湖

今日は琵琶湖調査の日。昨日までの大雪で北湖は沿岸から真っ白。今年は去年よりも雪が多くて寒いように感じる。
f:id:yokazaki:20170213193207j:plain
全層循環が完了して、表層から深層まで酸素たっぷり、水温は表層から湖底まで7.9℃で均一な状態。透明度は7.9mと高くなっていて、ずっと優占していた緑藻は勢力が衰えてきた模様。濾過もいつもよりもスムーズだった。細菌数は3.8×106 cells/ml、ウイルス数は8.7×107 /ml。
f:id:yokazaki:20170213193533j:plain
これで今シーズン調査はおしまい。この手作り濾過システムで1年間、水を採りまくった。毎月10Lを1,2回採って、60L採った月もあるから、延べ200Lは濾過しているのではないかと思う。今後しばらくはこのサンプルから生み出された膨大なシーケンスデータの解析に注力する。お世話になった船ともしばらくお別れ。
f:id:yokazaki:20170213194029j:plain
船着場には春がきておりました。
f:id:yokazaki:20170213194117j:plain
やっぱり現場に出てこそ生物の研究者だと思うので、早いところ分析したデータから新しい仮説を立てて、次の調査にまた来たい。琵琶湖という恰好の研究フィールドにいつでもどこでも簡単にアクセスできるこの環境は本当に素晴らしいと思う。

1月琵琶湖 もうすぐ全循環

 2017年の琵琶湖初め。先週末から続いていた冬の嵐で延期が続いていたけど、ようやく出ることができた。朝の気温は−3℃。船着場はスケートリンク状態で、ロープが凍った雪に埋まってしまっていた。
f:id:yokazaki:20170118173837j:plain
比良山も真っ白。12月に雪が少なくて苦戦していたスキー場も喜んでそう。
f:id:yokazaki:20170118174129j:plain
遠くの伊吹山もはっきりと。実は積雪世界一を記録した山だったりする。
f:id:yokazaki:20170118174318j:plain
この嵐でついに琵琶湖も全循環か、と期待したけど、惜しかった。船上で確認した溶存酸素のプロファイル(↓写真の赤い線)を見ても、湖底直上の数mだけ酸素が低くて、ぎりぎりまだ混ざりきっていないことが分かる。
f:id:yokazaki:20170118180056j:plain
表層水温が9.3℃、湖底が8.9℃。今月中には完全に混ざりそうですね。透明度は7.0mでした。
 これまでの研究で、琵琶湖には真冬の全循環の期間限定で優占する細菌がいることが分かっている。Acidobacteria門という土壌ではよく見られる系統の細菌なのだけど、水中に出てくるのはあまり聞かなくて、生態は分かっていない。この細菌が一体何をしているのか、今後採るサンプルで明らかに出来ればと思っている。

seqkitを用いたFASTAのフィルタリング・ソート

世の研究者たちはFASTAファイルのハンドリングにはどのようなツールを使っているのだろうか。

  • 長さでソートしたい
  • 特定の文字列を持つものだけ抽出したい
  • 名前だけ抽出したい
  • 名前を置換したい/通しで付け直したい
  • でかすぎるファイルを分割したい
  • ランダムサンプリングしたい

といった欲望は、シーケンスを扱うすべての研究者たちが持っていそうなものだけど、調べてみても、意外と統一的な方法って無くて、各自が自分なりの方法で作業してるのが実情みたいだ。自分自身もbioperl, biopython, bioawkなど、色々なツールをかじってみたけれど、どれもとっつきにくくて、結局、一番使い慣れているRにパッケージ"seqinr"を入れて作業する方法に落ち着いていた。けどこの方法は、もともと文字列処理が得意でないRに無理矢理FASTAを読み込ませているところもあるので、次世代シーケンサーで出てくるバカでかいFASTA/Qをいじるのには手数もかかるし、重くて結構ストレスが溜まる。結局、もっと効率の良い方法はないか、都度ググって解決する毎日を送っていた。

 そんな中、seqkitというツールを見つけて、試しに使ってみたのだけど、これがシンプルで速くてマジおすすめなので、使い方の備忘のついでに紹介。
bioinf.shenwei.me

元文献

seqkitは以下の論文で2016年10月に発表されている。
journals.plos.org
何ができるのかは、この論文のTable1にまとめられているほか、公式マニュアルの各項目にサンプル付きで詳しく説明されている。

特長

公式ページの"Futures"の項目をもとにポイントだけ訳すと

  • Windows, Mac, Linuxに対応
  • 一つの軽いバイナリファイルだけ入れればOKコンパイルや付随ソフトのインストール必要なし
  • とても速い、マルチスレッドにも対応
  • 20のコマンドが使えて、マニュアルも充実
  • FASTAもFASTQも読める。gzipで圧縮されていてもそのまま読める
  • 標準入力/出力に対応していて、パイプでコマンドを繋げられる
  • 正規表現が使える
  • 各シーケンスに固有の「ID」を振って分析することもできる
  • ランダムサンプリングやシャッフルコマンドでは、seedを設定でき再現性確保

インストール

超簡単。ダウンロードページからバイナリファイルを手に入れて、パスの通っている場所に置くだけ。Linuxでしか動かしていないけど、Windows版もコマンドプロンプトでexeにアクセスするところまではできたので、同じように使えると思う。

使い方

基本形は

seqkit <コマンド> <オプション> <標準入力> 

という形。これで標準出力がでるので、必要であれば"> output.fa"などを最後に付けて、アウトプットを保存する。入力にはFASTA, FASTQ, それを圧縮した.gzなどが入る。

(1)FASTAファイルの表示や分割

statFASTAのステータスを確認

seqkit stat input.fa

FASTAファイルのシーケンス数、長さの最小・最大・平均・合計を出してくれる。

ちなみに上記のとおり、パイプが使えるので、以下のように書いても同じ結果になる。複数コマンドを繋げたいときは便利(以下同様)

cat input.fa | seqkit stat

head:最初の数シーケンスだけ確認する

seqkit head -n 5 input.fa

シェルのheadと似たようなコマンド。分析の進捗を確認するのに、最初の数シーケンスだけ出力する。-nオプションで何シーケンス出すかを指定できる。

splitFASTAを分割する

seqkit split -s 100000 -2 input.fa

100000シーケンスずつFASTAを切って、カレントディレクトリに出力する。-sでファイル当たりのシーケンス数を指定。-2はtwo-pass modeといって、メモリの消費を抑えることができるらしい。なので、基本は指定しておいた方が良い(ただしFASTAでしか使えないとのこと)。-Oでアウトプットフォルダの指定もできる。

なお、シーケンス数でなく、ファイルの個数で分けたい場合は、-pを使う。例えば、5つに分割したければ、

seqkit split -p 5 -2 input.fa

その他、multipulexバーコードのような、指定した位置の塩基(アミノ酸)配列ごとにファイルを分割したり(-r)、各シーケンスの名前(ID)ごとに分割したり(-i)もできる。

sample:ランダムサンプリング

seqkit sample -n 100 -s 1234 -2 input.fa

-nでシーケンス数を指定する。この場合、インプットからランダムに100シーケンスを選んで出力する。-2は先に書いたものと一緒。-sはランダムサンプリングのseed。defaultは-s 11なので、何も指定しないと何度やっても同じランダムサンプリングの結果が出てくる。毎回結果を変えたいときは毎回違う値をセットする。

なお、シーケンス数ではなく、比率でサンプリング数を指定することもできる。例えば、全シーケンスの10%を出力したければ、

seqkit sample -p 0.1 -2 input.fa

ちなみに、FASTA内のシーケンスをランダムに並び替えるshuffleというコマンドもある。使い方は公式マニュアルを参照

fq2fa:FASTQをFASTAに変換

seqkit fq2fa input.fq.gz > output.fa.gz

そのまんま。.gzのまま使えるのが便利。

fx2tabFASTA/Qをtsvに変換

seqkit fx2tab input.fa

-nで名前だけの出力(シーケンス情報は載せない)にできて、-l-gなどのオプションを使って、長さやGC率などの情報とセットになったテーブルを出すことも可能。シーケンス名を一覧化してエクセルなどで触りたいときに使える。この逆の操作でtab2fxというのもあって、FASTAを長さでフィルタリングするときに使っている(以下参照)。

(2)FASTAの中身をいじる

seq:シーケンスの各種変換

このコマンドはたくさんオプションがあるので、詳細は公式マニュアルを参照。よく使いそうなやつだけリストにしておく。

オプション 機能
-w 改行までの文字数を数値で指定。0で改行無しのFASTAにできる
-p 相補配列を出力
-r 逆配列を出力
-pr 逆相補を出力
-l 小文字に変換
-u 大文字に変換
-n シーケンス名(">"以降)のみ抽出
-s シーケンス本体(">"の無い行)のみ出力
-g ギャップを削除して出力

なお、-wについてはglobal flagといって、seqkitの全コマンドに使えるオプションで、かつdefaultが60に設定されているので、パイプで処理を繋げるときは最後のコマンドに付けないと出力に反映されないので注意。
 ちなみにglobal flagは他にもいくつかあって(詳細は公式マニュアル参照)、そのうち-jではスレッド数を指定できる。defaultは2なんだけど、すでに十分に速いので、値を大きくしても(12とかにしても)、体感速度はそこまで変わらなかった。

sort:シーケンスのソート

seqkit sort -lr2 input.fa

これでシーケンスを長い順にソートできる。-lは長さでソートするためのオプション。そのままだと、昇順(短い順)で出てきてしまうので、-rを付けることで、逆向きにソートするよう指定すれば長い順になる。-2は上記と同じくtwo-pass modeの指定。
 あまり使わなさそうだけど、シーケンスの名前(-n)や、シーケンスそのもの(-s)の順に基づいてソートすることもできる。

(3)特定条件を満たすFASTAのフィルタリング(抽出)

grep:一致する文字パターンを持つシーケンスを抽出

seqkit grep -nrp mojiretsu input.fa

一番感動したのがこれ。これで"mojiretsu"を名前に持つシーケンスだけをフィルタリングできる。この1行で完結できる作業のためにこれまでどれだけの時間を費やしたことか。しかも速い!
-nでシーケンス名(">"以降の列)を対象にして、-pで検索対象の文字列パターンを指定。さらに-rを入れることで、正規表現が使えるようになる。

シーケンスを対象にしたフィルタリングも可能で、例えば

seqkit grep -sirp ^aggtctaa input.fa

で「AGGTCTAA」で始まるシーケンスだけを抽出できる。-sがシーケンスを検索対象にするオプション、-iで大文字小文字の違いを無視するオプション。

さらに縮重塩基や、検索対象の塩基の位置も指定できる。例えば、

seqkit grep -sidp GGGWGG -R 20:100 input.fa

で、「GGGWGG」が20~100塩基目の間に出現する配列だけを抽出できる。-dが縮重塩基を入力にするためのオプション、-Rがコロン区切りで検索対象塩基の位置の最初と最後を指定するオプション。ちなみに、-r-dを一緒に使うと怒られるんだけど、軽く使ってみた感じでは、-dを入れた状態で基本的な正規表現(^とか)は使えそう。完全に上位互換かどうかは調べられていないけど。
 あとはシェルのgrepと同じで、「一致しないもの」を抽出したいときは-vが使える。そのほか、検索したい複数のパターンをリスト化しておいて-fで読み込んで一気に検索、という使い方もできる。

■シーケンスを長さで抽出

2017/03/17追記
以前のバージョンでは長さでの抽出には一手間必要だったのだけど、v0.5.0からこれも一発でできるようになった。便利すぎる!

seqkit seq -m 1000 -M 10000 input.fa > output.fa

-mで最小の長さ、-Mで最大の長さを指定するだけ。上記の例では1000塩基(アミノ酸)以上、10000塩基(アミノ酸)以下のシーケンスのみが抽出できる。

(4)その他

まだ使いこなせていないけど、各シーケンスを名前とは別に「ID」という概念で管理して、ソートしたりフィルタリングしたりすることもできる。例えば、シーケンス名が

>AAAA | BBBBB-CCCC #DDDD

みたいな構造になっているときに、AAAAとかDDDDをキーにしてハンドリングすることができる。

その他まだ使うチャンスがないのだけど、他にも以下のようなコマンドが用意されている。

コマンド 機能
subseq 各シーケンスから特定の領域をのみを抜き出す。手入力の他、gtfやbedファイルでの指定も可能
sliding 各シーケンスを指定した長さ毎に区切る。場所ごとにGC率などを出したいときに使える
locate パターンに合った配列の位置をテーブルにまとめて出力。プライマー領域・ORF・モチーフ等の抽出に使える
rmdup 重複するシーケンス(duplicate)を削除する
common 複数のFASTAファイルに共通して存在するシーケンスを抽出
replace シーケンスやシーケンス名の置換。sedで出来ることの他に、通しで番号を振りなおしたりもできる


これまでに試したどのツールよりも圧倒的に使いやすい。FASTAイジリングツールはこれに落ち着きそうです。

琵琶湖納め

 今年最後の琵琶湖に行って来ました。今回はとにかく量が必要で、なんと60Lの採水。琵琶湖を60kgも軽くしてしまった。

f:id:yokazaki:20161226231925j:plain

気温7.7℃、水温11.6℃。透明度は5.5mとまだ低くて、緑藻のブルームが続いている。20µmのプランクトンネットをひくとあっという間に真緑に。陸の植物は枯葉になってしまっても、水の中はまだ青々としている↓

f:id:yokazaki:20161226231946j:plain

帰ってからはひたすら濾過。とにかく手を動かす。頭は使わない・・・

表層水を5.0μm のフィルターに通したところ。緑藻だけでなく、ミジンコもびっしり↓

f:id:yokazaki:20161226231948j:plain

ブルームの正体はMicrasterias。見た目はごついのだけど、横から見ると薄っぺらい。でも、顕微鏡で見るとみんな綺麗に向きがそろっていて、体が光を受けやすい向きになるように上手くできているのが分かる↓

f:id:yokazaki:20161226231952j:plain

こいつらどうやって増えるんだろう?と思っていたら、ちょうど分裂中のやつがいた。かわいい↓

f:id:yokazaki:20161226231950j:plain

ということで、今年の研究はおしまい。

来年もよろしくお願いします。年明けからかなり忙しそうな感じなので、年末年始はしっかりリフレッシュしてきます。

12月びわこ

今日は月例の琵琶湖の日。車も地面も凍り付く朝。比良山も白くなっていた。

f:id:yokazaki:20161212214526j:plain

表層水温は13.1℃、透明度は5mで緑藻が多い様子。前日の嵐のおかげか、表層から水深40mまで綺麗に水が混ざっている状態だった。全層循環までもう少し。しかし湖上は寒い・・・

RNAmmerのインストール

 FASTA形式のDNA配列からrRNA遺伝子を探してくるRNAmmer、元論文は2000回近く引用されていて、広く利用されているソフトだ。少量のシーケンス(1万配列、1000万塩基まで)であれば、ブラウザ上から直接配列を投げてサイト上で解析することもできて便利。

RNAmmer 1.2 Server

で、今回、大量のデータを分析しようと手元のLinuxにこれを入れようとして、色々と細かいところでハマったので備忘メモ。基本的には、開発者のブログに書いてある通りなのだけど、古いバージョンのソフトが必要だったり、手作業のステップもあったりして、結構ややこしい作業が必要。

  1. ダウンロードページの注意を読み、必要事項を入力して、tar.Zを一式をダウンロードして解凍
  2. RNAmmerを動かすのに必要なソフト、HMMERをインストールする。最新バージョンは3.1だけど、RNAmmerはv3以降では動かないので、ダウンロードページからv2.3.2を手に入れる。
  3. HMMERのマニュアルに従ってファイル一式を展開する。
  4. 必要なのは"hmmsearch"というバイナリファイルなので、これの場所を覚えておくか、任意の場所に移動させる。☆
  5. RNAmmerを解凍したフォルダに"core-rnammer"というスクリプトファイルがあるのを確認する。★
  6. RNAmmerを解凍したフォルダにある"rnammer"をテキストエディタで開き、# the path of the programの下にある、my $INSTALL_PATH =以降のパスを、RNAmmerを解凍したフォルダ★に書き換える。
  7. 同じく、 $HMMSEARCH_BINARY = 以降のパスをhmmsearchがある場所☆に書き換える。2か所あるが、Linuxを使っている場合はLinuxに関する部分だけ書き換えればOK。
  8. rnammerを閉じ、保存。
  9. これで動くらしいのだけど、
FATAL: POSIX threads support is not compiled into HMMER; --cpu doesn't have any effect

というエラーが出ることがあって(自分は出た)、その場合は"core-rnammer"をテキストエディタで開き、スクリプト内にある--cpu 1という記述(2か所)を削除し、保存。

と、これでようやく使えるようになる。

 アウトプットは様々な形式で出せるけど、見つかったrRNAの場所を記載したGFFと、その配列を抽出したFASTAを出したい場合は、作業フォルダに移動して、

perl RNAmmer解凍フォルダ/rnammer -S bac -m lsu,ssu,tsu -gff output.gff -f output.fasta input.fasta

とすればOK。-Sがどのドメインの生物を検索するか、-mがどのサブユニットのrRNAを検索するかというオプション。マルチスレッドには非対応だけど、そこそこ大きなFASTAファイル(10 MB)でも数分で分析できた。
開発者ブログによるとバージョンアップを進めているらしいので、今後に期待です。