塩基配列から塩基多様度などの統計量を求めるのは DnaSP がメジャーなプログラムですが，バージョン 5.10 では塩基多様度の分散を求めることができません。分散を求めるにはブートストラップを行えばよいのですが，ブートストラップ標本を生成するのはともかく，それをプログラムに渡していちいち塩基多様度を求めるのは非常に面倒だったので，自分で適当に実装してみました。

最初は実直に実装したため計算ごとに多型サイトを探していたのですが，サイトごとにペアワイズで比較していたために非常に計算コストが大きかったのです^[A]。サンプルが少量なら待てない時間ではないのですが，今後大量に処理したいことを考えると高速化の必要があったので，少し美しくなくなったのですが，多型サイトと配列長をキャッシュして計算してみました^[B]。すると劇的に速くなりました。よく考えれば当たり前のことです。配列長を L，多型サイト数を S とすれば，計算回数はほぼ S/L 倍になります。もし L = 500 bp で S = 10 だとしたら単純に考えて計算回数が 1/50 になるわけです。実際はブートストラップ標本の多型サイト数 S_B は S と等しくならないので再評価が必要になり，今の実装では S_B を 2 度評価しているため計算回数は 2S/L 程度です。さらに言えば多型サイトのキャッシュ作成が結構時間を食うため，全体の実行時間は 10 倍くらいにしかなりませんでした^[C]。

ともかくこれで大分実用的な速度になりました。今まではあまり計算速度とか気にしたことなかったのですが，やはり数値計算では効率を気にしないとだめですね。ちなみにマルチスレッド化は既に試していましたが^[D]， 4 コア CPU の Q6600 で 2.5 倍くらいでした^[C]。

[2010-01-21 追記] さらに見直して塩基多様度のペアワイズのイテレーション回数を減らすことでさらに高速化できました^[E]。

脚注

1. 単純なブートストラップ [↩]
2. 多型サイトをキャッシュしたブートストラップ [↩]
3. もちろんデータセットによって幅があります。 [↩] [↩]
4. マルチスレッドを用いたブートストラップ [↩]
5. 修正した塩基多様度計算方法 [↩]