Microsoft Biology Foundation (MBF) という Microsoft Biology Initiative 製の .NET バイオインフォマティクスライブラリがあります。これを使って塩基多様度を計算してみます。
塩基配列のデータはウェブ上に公開されている資料「適応進化遺伝学」の 34 ページものを利用します。具体的な塩基多様度計算例も載っています。
string fasta = @">Individual_1 ggtactccgTttgctcagcaTaaacttgccccaAtggactggaCctgatgggagaaactCgacataatTagatctacgaGtcatcagctc >Individual_2 ggtactccgTttgctcagcaTaaacttgccccaAtggactggaCctgatgggagaaactCgacataatTagatctacgaGtcatcagctc >Individual_3 ggtactccgTttgctcagcaTaaacttgccccaAtggactggaCctgatgggagaaactCgacataatTagatctacgaGtcatcagctc >Individual_4 ggtactccgTttgctcagcaCaaacttgccccaTtggactggaActgatgggagaaactCgacataatCagatctacgaGtcatcagctc >Individual_5 ggtactccgTttgctcagcaCaaacttgccccaTtggactggaActgatgggagaaactCgacataatCagatctacgaGtcatcagctc >Individual_6 ggtactccgCttgctcagcaCaaacttgccccaCtggactggaCctgatgggagaaactCgacataatCagatctacgaGtcatcagctc >Individual_7 ggtactccgCttgctcagcaCaaacttgccccaCtggactggaCctgatgggagaaactCgacataatCagatctacgaAtcatcagctc >Individual_8 ggtactccgCttgctcagcaCaaacttgccccaCtggactggaCctgatgggagaaactTgacataatCagatctacgaGtcatcagctc"; FastaParser parser = new FastaParser(); IList<ISequence> sequences = parser.Parse(new StringReader(fasta), true); int count = sequences.Count; int length = sequences[0].Count; // すべての配列の長さが等しく,位置もそろっていると仮定 double k = ParallelEnumerable.Range(0, length).Select(nucleotideIndex => { var nucleotideCounts = (from n in Enumerable.Range(0, sequences.Count).Select(sequenceIndex => sequences[sequenceIndex][nucleotideIndex]) group n by n.Symbol into g select g.Count()).ToArray(); if (nucleotideCounts.Length <= 1) { return 0; } int sum = nucleotideCounts.Sum(); return nucleotideCounts.Sum(c => c * (sum - c)); }).Sum() / (double)(count * (count - 1)); double pi = k / length; Console.WriteLine("π: {0:e1}", pi); // 3.7e-002
配列のペアワイズでの相違塩基数を数えるのではなく,あるインデックスにおける塩基を種類ごとにカウントして並列化しています (23, 33 行目)。これは無限サイトモデル (infinite-site model) を意識すれば自然な実装だと思います。