ハヤブサなみの応答性でデータファイルの精密調査

 昨晩は、コマンドを打ってから答えが返ってくるのに15分ほどかかるのを数回繰り返して早朝になってしまった。寝不足で出勤し、営業部詰めの合間に、awkの一行野郎様でファイルの構造の調査を続ける。
 サンプルは1被験者あたり18か所くらいからとっているので、ほしい場所由来の配列は全体(7000万本)のうちの485万本くらいで、1被験者あたり2万本程度シーケンスされていることがうかがわれた。その一方、中にはデータが数塩基長のものも散見され、サイズでふるい分けすることも必要になりそうである(このあたり本家ではQIIMEを使っているのであろうか)。
 ほしい配列のサブセットを得ることができれば、データのサイズは1/14になって、コマンドに対する応答が14倍速くなるに違いない。サブセットをどうやって得るか。もちろん、テキストエディタに読み込んで…というのはもとから念頭にない。頼みのawkで、一行ずつデータを読んで、ほしいものなら書き出し、不要なら捨てるというのを繰り返していけばよさそうに思われるが、FASTAの2行/レコードの形式のデータをそのまま扱うのはかなり難しい。そこで、後で復元できるように目印をつけながら、データの1レコードを1行にするのであるが、まずは改行記号を全てはぎ取って一行40ギガ文字の状態にする。
$tr '\n' '#' < SRP002395-7514-cs-nbp-rc.fsa > result.txt
 次に、データの先頭にだけ改行を戻すとよいわけであるが、データを一行(\nが出てくるまで)バッファに読み込んで処理する設計のソフトウェアフィルタではうまく動かない。\nを探しているとファイルの終わり(40ギガ文字目)まで読み込もうとしてsegvるのである。getchar()系のものを使うか、awkでレコードの区切り文字を定義して読み込むという手もあるだろう。1行/レコードの状態にすれば、awkで、ほしいサンプリング部位のシーケンスを取りまとめて、リダイレクトで一つのファイルに書き出すのはわけもない。生成したファイルのサイズは予想通り2〜3 GB程度におさまっている。
 と書いたところで、最初からawkでレコード区切り文字を改行以外のキャラクタに替えたら、ほとんど労なくここまで来られたのではないかという疑問がわいてきたが、先を急ぐので未来の自分への宿題ということにしておく。
 

本ブログではamazon associate広告を利用しています。