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ファイルの表示や分割

■stat：FASTAのステータスを確認

seqkit stat input.fa

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

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

cat input.fa | seqkit stat

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

seqkit head -n 5 input.fa

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

■split：FASTAを分割する

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のまま使えるのが便利。

■fx2tab：FASTA/Qをtsvに変換

seqkit fx2tab input.fa

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

(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をキーにしてハンドリングすることができる。

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

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

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

琵琶湖納め

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

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

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

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

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

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

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

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

12月びわこ

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

表層水温は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)でも数分で分析できた。
開発者ブログによるとバージョンアップを進めているらしいので、今後に期待です。

11月びわこ

今日は月例の琵琶湖調査。

気温14.7℃、水温15.9℃、透明度は5.5m。深層の細菌密度は1.0×10^6 cells/ml, ウイルス現存量は2.0×10^7 particles/ml。水温躍層は30m付近まで下がってきていて、昨年よりも確実に冷えるのが早い。去年が暖冬だったおかげで深層の水温が高いこともあって、今年は早めに全層循環が起こりそうだ。

　紅葉もちょうど見ごろを迎えていて、これも去年よりも1週間くらい早い。去年は温度の下がり方が鈍くて、紅葉が汚かったけど、今年は一気に冷えたおかげでとてもきれいだ。こちらは先週の信楽の様子↓

　今年は春が長かったし、梅雨もしっかりあって、8月はほとんど雨が降らず、その後しっかり秋雨が降り続けて、一気に気温が下がって、四季の移り変わりがはっきりしている気がする。冬も冷えるのかな。

系統地理ワークショップ＠京都

　京大のセミナーハウスで行われていた生物群横断系統地理ワークショップなるイベントに参加してきた。僕の研究では「なぜ同じ系統の細菌が世界中の別々の湖に共通して生息しているのか？」という疑問があって、他の生物の研究から何かヒントは得られないか、ということで今回情報収集を目的に参加してみた。

　自分の研究にダイレクトにつながるようなアイデアは得られなかったけど、普段とは全く異なる研究者コミュニティに入って、系統地理界隈のトレンド（何が分かっていて、何が分かっていないのか？）を大まかに把握できたのは良かった。

　もう少し突っ込んだ、内容とは別の感想では、直前に参加した微生物生態学会などと比べて、研究のモチベーションが純粋である、ということが印象深かった。講演やポスター発表でも何度も「面白い」という言葉が出てきたし、「その生き物が好きだから研究している」という人が多いように感じた。良くも悪くも「お金の臭い」が感じられない研究だと感じた。僕は基礎研究の将来には悲観的で、「面白いからやる研究」を「趣味の延長」みたいに捉えていて、税金を出すことを快く思わない人が社会にはたくさんいるし、これからも少子高齢化で財政が厳しくなっていく中で、そういう研究にお金を出す余裕はますます無くなっていくと思っている。だから、ワークショップの場でも議論があったけど、「なんの役に立つのか？」ということを急いで議論しないとマズいのではないかという不安を感じてしまった。新しいものを対象にすれば、面白い発見が見つかるのは当たり前であって、本当の戦いは「たくさんの面白い発見の中での優先順位をいかに上げるか」というところだ。面白いことももちろん大切なのだけど、残念ながら無い袖はふれない。「それは税金でやることではない」と言われたら、悔しいけど黙るしかない。もし「面白い」を前面に出すのであれば、クラウドファンディングみたいに「面白い」と感じてくれる人から直接研究費を貰うような方法を考えていくしかないと思う。

　じゃあ、どうやって「役に立つ」ことをアピールすればよいのだろうか。僕が今回この集会に素人なりに参加して感じたのは、「系統地理」って言うけど、色んな生き物、色んなレベルの変異、色んな地理スケールでやっている人が集まっているので意外とその定義は難しくて、結局のところ「野外での遺伝子の移動を追いかけるプロの集まり」がこの学問分野なのではないかということだ。だから、有害な遺伝子（耐性遺伝子とか）がどう広がっていくのか、とか、有用な遺伝子のソースが失われないようにするにはどうしたらいいか、みたいなところで、統一見解を形作っていければ、この分野のプレゼンスもより高まるのではないかと感じた。「そんなこととっくに考えてるよ」と言われそうですが。

　そしてやっぱり大きな生き物の研究と比較して感じるのは、「微生物って絶望的なくらい何もわかっていない、でもだからこそ、これからの伸びしろも大きそうだ」ということだ。そもそも、系統地理の研究で重要な単位である「種」の概念すら、細菌の世界では訳の分かっていない状況だ。大きな生き物の研究では当たり前のように出てくる「違う場所で採れるものは遺伝子レベルで見ると全然違う」みたいな話が、細菌でみつかるとトップジャーナルに乗るような大発見として扱われる*1。次世代シーケンサーのおかげで、今まさにボーナスステージで、そういう発見がものすごい勢いでなされている。僕も、大きな生き物の研究にヒントを得ながら「細菌の移動」の謎を解く発見をしていきたい。

*1:

Global biogeography of SAR11 marine bacteria. - PubMed - NCBI

Global biogeography of Prochlorococcus genome diversity in the surface ocean. - PubMed - NCBI

Global phylogeography of pelagic Polynucleobacter bacteria: restricted geographic distribution of subgroups, isolation by distance and influence of... - PubMed - NCBI

微生物生態学会@横須賀

　微生物生態学会＠横須賀に参加。この学会は、一番僕の研究分野にフィットした学会であると同時に、ものすごく守備範囲の広い学会でもある。感動するくらいマニアックなレベルで自分の研究を評価してくれる人もいる一方で、大多数の研究者は僕のポスターの前を素通りだ。

　微生物の研究は、「新手法⇒新発見」という図式が成り立ちやすい。だからどうしても記載的な研究が多い。（微生物以外の）生態学者がこの学会を「つまらない」と評するのを何度も聞いたことがあるけれど、この点が原因だと思う。だけど、だからこそ、少しでも生態や進化に絡めた議論に持ち込めればインパクトのある研究になる、まだまだ発展の余地がある面白い分野だと思う。僕の研究も今はとても記載的だ。だけどこれは基礎固めであって、今後、もっと一般的で面白い発見を量産できる系に発展させていくつもりだ。そうすれば、もう少しポスターの前で立ち止まってくれる人も増えるだろう。

　これまで参加してきた学会では人脈を広げることに積極的だったけど、今回はほとんど名刺交換をしなかった。それだけ顔見知りが増えたということもあるのだけど、それ以上に、自分の立ち位置が分かってきて、広げていくべき人脈の方向性や範囲が絞れてきたことが大きいように思う。新たな成果を出すことなしに、これ以上人脈だけ広げても意味がないし、迷惑なだけだ。すでに関連分野の研究者には自分の研究のことを一通り知ってもらえている。だからこの先は、成果で語るしかない。早く研究を先に進めて、記載的でない、この研究の本当の味を出したい。