塩基多様度の計算 - MBF Cooking


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) を意識すれば自然な実装だと思います。