今日は琵琶湖調査の日。昨日までの大雪で北湖は沿岸から真っ白。今年は去年よりも雪が多くて寒いように感じる。
全層循環が完了して、表層から深層まで酸素たっぷり、水温は表層から湖底まで7.9℃で均一な状態。透明度は7.9mと高くなっていて、ずっと優占していた緑藻は勢力が衰えてきた模様。濾過もいつもよりもスムーズだった。細菌数は3.8×106 cells/ml、ウイルス数は8.7×107 /ml。
これで今シーズン調査はおしまい。この手作り濾過システムで1年間、水を採りまくった。毎月10Lを1,2回採って、60L採った月もあるから、延べ200Lは濾過しているのではないかと思う。今後しばらくはこのサンプルから生み出された膨大なシーケンスデータの解析に注力する。お世話になった船ともしばらくお別れ。
船着場には春がきておりました。
やっぱり現場に出てこそ生物の研究者だと思うので、早いところ分析したデータから新しい仮説を立てて、次の調査にまた来たい。琵琶湖という恰好の研究フィールドにいつでもどこでも簡単にアクセスできるこの環境は本当に素晴らしいと思う。
1月琵琶湖 もうすぐ全循環
2017年の琵琶湖初め。先週末から続いていた冬の嵐で延期が続いていたけど、ようやく出ることができた。朝の気温は−3℃。船着場はスケートリンク状態で、ロープが凍った雪に埋まってしまっていた。
比良山も真っ白。12月に雪が少なくて苦戦していたスキー場も喜んでそう。
遠くの伊吹山もはっきりと。実は積雪世界一を記録した山だったりする。
この嵐でついに琵琶湖も全循環か、と期待したけど、惜しかった。船上で確認した溶存酸素のプロファイル(↓写真の赤い線)を見ても、湖底直上の数mだけ酸素が低くて、ぎりぎりまだ混ざりきっていないことが分かる。
表層水温が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"の項目をもとにポイントだけ訳すと
インストール
超簡単。ダウンロードページからバイナリファイルを手に入れて、パスの通っている場所に置くだけ。Linuxでしか動かしていないけど、Windows版もコマンドプロンプトでexeにアクセスするところまではできたので、同じように使えると思う。
使い方
基本形は
seqkit <コマンド> <オプション> <標準入力>
という形。これで標準出力がでるので、必要であれば"> output.fa"などを最後に付けて、アウトプットを保存する。入力にはFASTA, FASTQ, それを圧縮した.gzなどが入る。
(1)FASTAファイルの表示や分割
■stat
:FASTAのステータスを確認
seqkit stat input.faFASTAファイルのシーケンス数、長さの最小・最大・平均・合計を出してくれる。
ちなみに上記のとおり、パイプが使えるので、以下のように書いても同じ結果になる。複数コマンドを繋げたいときは便利(以下同様)
cat input.fa | seqkit stat
■head
:最初の数シーケンスだけ確認する
seqkit head -n 5 input.faシェルの
head
と似たようなコマンド。分析の進捗を確認するのに、最初の数シーケンスだけ出力する。-n
オプションで何シーケンス出すかを指定できる。
■split
:FASTAを分割する
seqkit split -s 100000 -2 input.fa100000シーケンスずつ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のまま使えるのが便利。
(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
で読み込んで一気に検索、という使い方もできる。
(4)その他
まだ使いこなせていないけど、各シーケンスを名前とは別に「ID」という概念で管理して、ソートしたりフィルタリングしたりすることもできる。例えば、シーケンス名が
>AAAA | BBBBB-CCCC #DDDD
みたいな構造になっているときに、AAAAとかDDDDをキーにしてハンドリングすることができる。
その他まだ使うチャンスがないのだけど、他にも以下のようなコマンドが用意されている。
コマンド | 機能 |
---|---|
replace | シーケンスやシーケンス名の置換。sed で出来ることの他に、リストを参照して複数の置換を一気に行ったり、通しで番号を振りなおしたりもできる |
subseq | 各シーケンスから特定の領域をのみを抜き出す。手入力の他、gtfやbedファイルでの指定も可能 |
sliding | 各シーケンスを指定した長さ毎に区切る。場所ごとにGC率などを出したいときに使える |
locate | パターンに合った配列の位置をテーブルにまとめて出力。プライマー領域・ORF・モチーフ等の抽出に使える |
rmdup | 重複するシーケンス(duplicate)を削除する |
common | 複数のFASTAファイルに共通して存在するシーケンスを抽出 |
これまでに試したどのツールよりも圧倒的に使いやすい。FASTAイジリングツールはこれに落ち着きそうです。
琵琶湖納め
今年最後の琵琶湖に行って来ました。今回はとにかく量が必要で、なんと60Lの採水。琵琶湖を60kgも軽くしてしまった。
気温7.7℃、水温11.6℃。透明度は5.5mとまだ低くて、緑藻のブルームが続いている。20µmのプランクトンネットをひくとあっという間に真緑に。陸の植物は枯葉になってしまっても、水の中はまだ青々としている↓
帰ってからはひたすら濾過。とにかく手を動かす。頭は使わない・・・
表層水を5.0μm のフィルターに通したところ。緑藻だけでなく、ミジンコもびっしり↓
ブルームの正体はMicrasterias。見た目はごついのだけど、横から見ると薄っぺらい。でも、顕微鏡で見るとみんな綺麗に向きがそろっていて、体が光を受けやすい向きになるように上手くできているのが分かる↓
こいつらどうやって増えるんだろう?と思っていたら、ちょうど分裂中のやつがいた。かわいい↓
ということで、今年の研究はおしまい。
来年もよろしくお願いします。年明けからかなり忙しそうな感じなので、年末年始はしっかりリフレッシュしてきます。
RNAmmerのインストール
FASTA形式のDNA配列からrRNA遺伝子を探してくるRNAmmer、元論文は2000回近く引用されていて、広く利用されているソフトだ。少量のシーケンス(1万配列、1000万塩基まで)であれば、ブラウザ上から直接配列を投げてサイト上で解析することもできて便利。
で、今回、大量のデータを分析しようと手元のLinuxにこれを入れようとして、色々と細かいところでハマったので備忘メモ。基本的には、開発者のブログに書いてある通りなのだけど、古いバージョンのソフトが必要だったり、手作業のステップもあったりして、結構ややこしい作業が必要。
- ダウンロードページの注意を読み、必要事項を入力して、tar.Zを一式をダウンロードして解凍
- RNAmmerを動かすのに必要なソフト、HMMERをインストールする。最新バージョンは3.1だけど、RNAmmerはv3以降では動かないので、ダウンロードページからv2.3.2を手に入れる。
- HMMERのマニュアルに従ってファイル一式を展開する。
- 必要なのは"hmmsearch"というバイナリファイルなので、これの場所を覚えておくか、任意の場所に移動させる。☆
- RNAmmerを解凍したフォルダに"core-rnammer"というスクリプトファイルがあるのを確認する。★
- RNAmmerを解凍したフォルダにある"rnammer"をテキストエディタで開き、
# the path of the program
の下にある、my $INSTALL_PATH =
以降のパスを、RNAmmerを解凍したフォルダ★に書き換える。- 同じく、
$HMMSEARCH_BINARY =
以降のパスをhmmsearchがある場所☆に書き換える。2か所あるが、Linuxを使っている場合はLinuxに関する部分だけ書き換えればOK。- rnammerを閉じ、保存。
- これで動くらしいのだけど、
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)でも数分で分析できた。
開発者ブログによるとバージョンアップを進めているらしいので、今後に期待です。
11月びわこ
今日は月例の琵琶湖調査。
気温14.7℃、水温15.9℃、透明度は5.5m。深層の細菌密度は1.0×10^6 cells/ml, ウイルス現存量は2.0×10^7 particles/ml。水温躍層は30m付近まで下がってきていて、昨年よりも確実に冷えるのが早い。去年が暖冬だったおかげで深層の水温が高いこともあって、今年は早めに全層循環が起こりそうだ。
紅葉もちょうど見ごろを迎えていて、これも去年よりも1週間くらい早い。去年は温度の下がり方が鈍くて、紅葉が汚かったけど、今年は一気に冷えたおかげでとてもきれいだ。こちらは先週の信楽の様子↓
今年は春が長かったし、梅雨もしっかりあって、8月はほとんど雨が降らず、その後しっかり秋雨が降り続けて、一気に気温が下がって、四季の移り変わりがはっきりしている気がする。冬も冷えるのかな。