世の研究者たちは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イジリングツールはこれに落ち着きそうです。
