10年以上前に書いたhcaslという小ツールについて「これなあに?」と問い合わせが来た。備忘録として裏事情を含めて書き残しておこうと思う。
hcaslには元ネタとなる自作ツールがある。元ネタの自作ツールは、かつて「教えて! goo」が生きていた頃にプログラミングのカテゴリに登録された質問を見て個人的に実装してみた代物だ。
2000年代後半のことだが、当時の私はプログラミングを始めてから日が浅く、しかもC言語ぐらいしか使える言語が無かった。そこで一時期、修行の一環として、「教えて! goo」などの質問掲示板に投稿されたプログラミング系の質問を読んで解答となるプログラムをC言語で書く、といったことを繰り返していた。
パズル的な課題や、現在で言えば競プロで出てくるような問題が、自分自身の好みから外れていたこともあるが……修行において主眼に置いていた項目が「実用に耐えうる堅牢なコードを書く」「分かりやすく明快なコードを書く」「自分以外のプログラマがメンテしやすいコードを書く」であったこともあり、どちらかといえば「何かしら実用的な小ツールを作る」系の課題を求めていた。で、当時の質問掲示板には時々その手の質問が投稿されていたのだ。
Wayback Machineに当時の質問が残っていた:
若干エスパーしつつ内容をまとめると、DNAをACGTの4文字で表現したテキストファイル――恐らくFASTAフォーマットからヘッダ行を取り除いたようなものだと思われる――を入力としてk-merを生成して、その出現回数をカウントしたい、ということだろう。質問ではabcdの4文字で構成されたテキストファイルとされているが、実際にはDNAを構成する代表的な塩基を示すACGTの4文字で構成されていたはずだ(何しろ質問者は生化学系の分野の人らしいのだから)。部分配列の長さは8文字、つまり8-merの頻度解析である。
この課題を解くにあたり、最初に思いついたアイデアは「CLIアプリを2つ作ってパイプで連結する」だった。概ねこんな感じだ:
- CLIアプリ
8merを実装する。8merは標準入力からテキストファイルを読み込み、8-merを「1行に1個」の形式で標準出力に書き出す。 - CLIアプリ
8freqを実装する。8freqは標準入力から8merの出力を受け取り、8-merの出現回数をカウントして出力する。 - 上記のプログラムを
./8mer <input.txt | ./8freqのようにパイプで組み合わせて使用する。
ここで言うCLIアプリ8merを見直して、任意の文字数の部分文字列を出力可能にしたのがhcaslだ。
当時も今も、2つに分割するアイデアは「性能面のボトルネックが問題とならなければ」悪くないアイデアだと考えている。
8merは単純に8-merを生成するだけだが、言い換えれば組み合わせるツール次第で「出現回数のカウント」以外の作業にも流用できる。当初は8-mer固定で流用性は低かったが、再実装したhcaslのように任意のk-merを生成できるのならば、流用可能な範囲は広がるはずだ。
8merもhcaslも、生成されたk-merの出力順序は入力データに依存する。もしも「入力シーケンスの中のどのあたりに出現したか?」という情報が必要なケースが存在するならば、そういった類の作業にも使えるはずだ。不要ならばsortなどの他のツールで整列させてしまえばよい。
また8freqについては、例えば出現回数が0の8-merを表示しなくても構わないのなら、Unix環境では./8mer <input.txt | sort | uniq -cと処理してしまえばよい。つまり8freqは不要となる。
出現回数が0の8-merも含めて結果を出力したい場合も、一時ファイルを使用して構わないならば、一般的なPOSIXコマンドの組み合わせで実現できる。8freqのような専用ツールを書く必要はない。
例えば以下のような感じに8-merの全パターンを網羅したソート済みファイル8-mer-all.txtを用意しておいて:
AAAAAAAA AAAAAAAC AAAAAAAG AAAAAAAT AAAAAACA AAAAAACC
(中略)
TTTTTTGG TTTTTTGT TTTTTTTA TTTTTTTC TTTTTTTG TTTTTTTT
先述のsort | uniq -cとの組み合わせで得られた8-merの出現回数の結果とjoinしてしまえばよい:
# 分かりやすく分割して書くとこんな感じ。 # # 8-mer-all.txtの各行の形式 : `<8-mer>` # 8-freq.txtの各行の形式 : `<出現回数> <8-mer>` # # join(1) のオプション: # * `-1 1 -2 2` より、8-mer-all.txtの1列目と8-freq.txtの2列目をキーとして結合する。 # * `-o 0,2.1` より、出力の形式は `<キーとする列> <8-freq.txtの1列目>`、 # すなわち `<8-mer> <8-freq.txtの<出現回数>>` とする。 # ちなみに `-o 1.1,2.1` と指定してもよい。 # * `-a 1` より、キーが一致しなかった場合には、 # 該当する8-mer-all.txtの列を出力させる。 # これにより全ての `<8-mer>` が必ず出力されることになる。 # * `-e 0` より、出力時に「空のフィールド」があったら `0` を出力させる。 # `-a 1` との組み合わせより、キーが一致しなかった場合に # `<8-freq.txtの<出現回数>>` の代わりに `0` が出力される。 # すなわち `<8-mer> <0>` となる。 # ./8mer <input.txt | sort | uniq -c >8-freq.txt join -a 1 -e 0 -o 0,2.1 -1 1 -2 2 8-mer-all.txt 8-freq.txt # ワンライナー化したバージョン: ./8mer <input.txt | sort | uniq -c | join -a 1 -e 0 -o 0,2.1 -1 1 -2 2 8-mer-all.txt -
この件から得られる教訓は何だろうか? それは「足りない部分だけプログラムを書く」ということだ。
とかく我らサンデープログラマはせっかちですぐに手を動かしたくなるものだが、一呼吸いれてから解くべき課題を注意深く観察してみると、本当にコードを書かなくてはならない部分は少なく、残りの大半は既存のツールの組み合わせで済んでしまう、といったケースは意外と多い。
今回の件で言えば、8merにあたる機能は用意する必要があったが、8freqについては実装せずに済ます余地があった。実装する範囲が狭まれば作業が楽になるし、バグを仕込んでしまう可能性も少なくなるはずだ。
その上で、「足りない部分だけプログラムを書く」にあたり、他のプログラムと容易に組み合わせられるように設計する必要がある。そうしなければ、他のプログラムと連携できない。
8merもhcaslもテキストフィルタだ。Unix環境においては様々なコマンドと組み合わせて扱うことができるだろう。
職業プログラマの場合は諸事情により「他の出来合いのプログラムと組み合せる」ことが厳禁なことも多いが、個人的な課題を解決するためのプログラミングにおいては、その辺の自由度は高いはずだ。
……まあ、今となっては生化学も含めてどこもかしこもPythonを使っているだろうから、hcaslなんてお呼びじゃないだろう。
蛇足:ちなみに8merやhcaslを実装した際の裏テーマは「入力ファイルに書かれているテキストが無茶苦茶長くても異常終了しない」だった。元の質問に入力ファイルの最大サイズの指定が無かったので、ヒトゲノムを基準として約32億文字のテキストが突っ込まれてもハングアップしない(ただし処理完了まで無茶苦茶時間がかかるけど……)ように実装した。ここまでデータ量が多い場合には、sort | uniq -cと組み合せると大変な事態に陥りそうなので、出現回数をカウントするために8freqのような専用ツールを実装する意味があるはずだ。