複数の配列のFASTAを分割する - アメリエフの技術ブログ

1つのファイルに複数の配列が書いてある下記のようなFASTAファイルを考えます。

> seq1
ATTGGCTGCCGCGCGGGGCGGGGAGCGGGGTCGGCTCAGTGGCCCTGAGACCCTAGCTCT
GCTCTCGGTCCGCTCGCTGTCCGCTAGCCCGCTGCGATGTTGCGCGCTGCCGCCCGCTTC
GGGCCCCGCCTGGGCCGCCGCCTCTTGTCAGCCGCCGCCACCCAGGCCGTGCCTGCCCCC
> seq2
TTGGCGAAGCTCTTTCATTAACATCATCCCATTTAATCTTGGTTAAGCTGTAATGGATTG
GTTGAGTAATTTTTTCAAAACAAAGGATGAACAAGCTCAGAGAGCCTGAGAGACTTACTG
AAGATCACACAGTAAATGATAAACTATTTTACTTTATGTCTAAGGTCTTTCATAATATGA
> seq3
TGCTCAGGAAATAAGAGCTGCTATTATTACTACTCTTCATAGCATAATGTAAGGGTTTCC
TTCCTCTCATCCCAGCTCGGAAGCACCTTATCAAGGTGAACTGAGCCGCTTCTGTTTCTT
CTCTTCCCTCTAGCCACCACGCCCCTTCCTTTCCGCCATCCTGCCCCCAGTCTTAGCTGC

これを1配列ごとに、異なるファイルに分割したいことってよくありますよね。
こんな感じです。

seq1.fa

> seq1
ATTGGCTGCCGCGCGGGGCGGGGAGCGGGGTCGGCTCAGTGGCCCTGAGACCCTAGCTCT
GCTCTCGGTCCGCTCGCTGTCCGCTAGCCCGCTGCGATGTTGCGCGCTGCCGCCCGCTTC
GGGCCCCGCCTGGGCCGCCGCCTCTTGTCAGCCGCCGCCACCCAGGCCGTGCCTGCCCCC

seq2.fa

> seq2
TTGGCGAAGCTCTTTCATTAACATCATCCCATTTAATCTTGGTTAAGCTGTAATGGATTG
GTTGAGTAATTTTTTCAAAACAAAGGATGAACAAGCTCAGAGAGCCTGAGAGACTTACTG
AAGATCACACAGTAAATGATAAACTATTTTACTTTATGTCTAAGGTCTTTCATAATATGA

seq3.fa

> seq3
TGCTCAGGAAATAAGAGCTGCTATTATTACTACTCTTCATAGCATAATGTAAGGGTTTCC
TTCCTCTCATCCCAGCTCGGAAGCACCTTATCAAGGTGAACTGAGCCGCTTCTGTTTCTT
CTCTTCCCTCTAGCCACCACGCCCCTTCCTTTCCGCCATCCTGCCCCCAGTCTTAGCTGC

検索してみると、様々なツールで実現可能でしたので、それぞれの特徴を調べてみました。

今回は3つとりあげますが、ここに書いたツール以外にもたくさんあります。

BBtools
中身はjavaなので、java環境が必要です。
バイナリが配布されているので、ダウンロードして解凍してください。たくさんツールがあるのですが、そのなかの「partition.sh」を使います。
まず、 ways= で最終的に作成したいファイル数を指定する必要があります。grepコマンドなどで配列数を数えておく必要があります。-cオプションを使うと便利です（参考）。
また、今回GRCm38のゲノムを使ったのですが、デフォルトではメモリ不足で落ちてしまいました。ゲノムサイズが大きいときは-Xmxでjavaの最大ヒープサイズを適宜指定する必要があります。
出力ファイル名はout=で指定します。%に、通し番号が振られます。しかも0始まりです。元ファイルの配列名でファイルを作成したいときは少し不便。
ついでに、これによって困ることはないと思いますが、元FASTAファイルは50文字で折り返していたのですが、70文字折り返しになっていました。

$ partition.sh in=genome.fa out=chr%.fa ways=22 -Xmx50G

Kent's src
ソースはここから。OSに合わせたバイナリ版があります。
今回使用したいツールは faSplit 。自分の環境に合わせたバイナリをダウンロードした後 chmod 755 faSplit で実行権限を与えてください。
ファイル名と出力ディレクトリ名を指定して実行します。また、FASTAの分割の方法について 'about' 'byname' 'base' 'gap' 'sequence' 'size' のいずれかが選択できます。今回は配列ごとに分けたいので'byname'を選択しますが、「100個のファイルに分けたい」「2000bpごとに分割したい」などの方法も選べるわけです。
'byname'を選択しているときは、配列名.faが作成されます。

$ faSplit byname genome.fa ./

GenomeTools
ソースコードが配布されているので、
- ダウンロードする
- 解凍する
- makeでコンパイルする
  のいつもの3ステップでOKです。binディレクトリにあるgtファイルが実行ファイルです。たぶんGenomeToolsの略でしょう。たくさんの機能があるので-helpでusageを確認して、使いたい機能を探します。今回は'splitfasta'機能を使います。
  これも分割したい数や塩基配列長など、分割の方法を選べます。今回のように、配列ごとに分割する場合は-splitdescオプションを使用し、出力ディレクトリを指定します。これも配列名.faが作成されます。
  なお、元FASTAは50文字折り返しだったのですが、出力ファイルは改行無し（1行1配列）になっていました。これはちょっと使いにくいことがありそう。

$ gt splitfasta -splitdesc  ./ genome.fa

比較
この3ツール、「何か便利なのないかな～」という動機で調べただけなので1回実行しただけでまじめな検証はしていませんが、GRCm38（約2.7Gbp）を分割した時にかかる時間と最大メモリ使用量を計測してみたところ、以下のとおりになりました（計測方法はこちら）。
実際にはFASTAファイルのサイズのほか、配列の数など、様々な要因が影響するかもしれません。

ツール名	実行時間（s）	最大メモリ使用量（KB）
BBtools	24.75	6,624,768
Genome Tools	171.39	985,804
Kent's src	25.67	263,824

あ、ちなみにseqkitという、FASTAやFASTQの操作を行える、機能豊富なツールも人気がありますが、小さいFASTAファイルは問題なく分割できたのですが、マウスの全ゲノムFASTAを配列ごとに分割しようとしたら、メモリ不足でエラーが出てしまいました（メモリ使用量を指定するオプションもありませんでした）。
塩基数ごとなど、別の分割方法なら問題なく使えるかも？

では最後の選択肢。

自分で書く

前の記事で、「バイオインフォマティシャンは、つい自分でファイル変換スクリプトを書こうと思ってしまう……」的なことを言ったばかりですが……FASTAについては事情があります。

FASTAファイルは割と無法地帯。
あるいは独特のルールが乱立している状態なので、世の中のFASTAは配列名などに思いもよらない文字が混ざっていたり、望んでいた配列名ではなかったりします。
上記のどのツールでも、自動分割後のファイル名は、だいたい配列名をもとに新しいファイル名を決めるので、変な文字が入っているときは新しいファイルを作成する際に不都合があったり、変な名前のファイルが生成されて後で一生懸命リネームしなければならなくなったりします（経験談）。
そんなトラブルに臨機応変に対応することができるのが、自分で書いたプログラムの強み！

webで検索すると、FASTAを自動で分割してくれるperlやpythonスクリプトを公開している人がたくさんいて、みんな自分で書いてるんだな…と実感しました。いずれもそれほど複雑ではありませんでした。
BioPerlやBioPythonを使えば簡単に書けそうですし、正直それらを使うまでもないほど簡単。プログラミングが得意な人なら、これくらいのシンプルな機能は自分で書く、という人も多いようです。

今回は「マルチFASTAを配列ごとに分割する」方法について紹介しましたが、FASTAを分割するときはほかにも100kbpごとの分割、など様々な目的があると思います。