コアレセントシミュレーターを実装する (1)

02 5月, 2011

コアレセントシミュレーションを行うアプリケーションで有名なものに ms[A] があります。シミュレーションを行って Tajima's D を計算するような単純なアプローチであれば ms (および同梱されている sample_stats) を使えば十分なのですが，モデルが複雑になってくると，いちいち ms が出力するテキストをパースして解析するといったアプローチはむしろ面倒になってきます。

そこで，コアレセントシミュレーターを C# で実装してみたいと思います。とりあえず最初は一番単純な，単一集団・集団サイズ一定・突然変異なしのモデルでの実装を目標にします。 ms では ms nsam nreps -T に相当するものです。

コアレセントツリー

コアレセント理論の根幹をなすのはコアレセントツリーです。コアレセントツリーは情報科学で言うところの二分木なので，実装は難しくありません。木そのものを表すクラスとノードを表すクラスを分けるのが自然です。

まずはノードクラスを実装します。ノードはサンプルを表します。葉ノード，つまり末端にあるノードは手元にあるサンプルです。また，内部ノードは，子ノード同士の最近の共通祖先 (MRCA) になります。

必要なものは，それぞれのノードの親ノード・子ノード・コアレセントイベントまでの待ち時間 (waiting time) ですので，それぞれ実装します。また，出力を考えてサンプル名を指定できるようにしています。基本的にノード操作を外部から行うことはないので，各プロパティのセッターは internal にして保護しています。

public class CoalescentTreeNode
{
   public string Name { get; set; }
 
   private CoalescentTree tree;
   public CoalescentTree Tree
   {
      get { return this.tree; }
      internal set { this.tree = value; }
   }
 
   private CoalescentTreeNode left;
   public CoalescentTreeNode DescendantLeft
   {
      get { return this.left; }
      internal set { this.left = value; }
   }
 
   private CoalescentTreeNode right;
   public CoalescentTreeNode DescendantRight
   {
      get { return this.right; }
      internal set { this.right = value; }
   }
 
   private CoalescentTreeNode ancestor;
   public CoalescentTreeNode Ancestor
   {
      get { return this.ancestor; }
      internal set { this.ancestor = value; }
   }
 
   private double time;
   public double WaitingTime
   {
      get { return this.time; }
      internal set { this.time = value; }
   }
 
   public bool IsLeaf
   {
      get
      {
         return DescendantLeft == null && DescendantRight == null;
      }
   }
 
   internal bool HasDescendant
   {
      get
      {
         return DescendantLeft != null && DescendantRight != null;
      }
   }
 
   public CoalescentTreeNode()
   {
   }
 
   internal double GetTotalTime()
   {
      double sum = 0.0;
      CoalescentTreeNode current = DescendantLeft;
      while (current != null)
      {
         sum += current.WaitingTime;
         current = current.DescendantLeft;
      }
      return sum;
   }
}

IsLeaf のようなプロパティは外部から呼ばれることが予想されますが， HasDescendant は内部的にしか使わないと考えられますので， internal にしてあります。 GetTotalTime は過去に遡る方向でのノードの分岐時間を取得するためのメソッドです。

次にツリーを実装します。

public class CoalescentTree
{
   private readonly List<CoalescentTreeNode> nodes;
 
   public CoalescentTreeNode Root
   {
      get { return IsIncomplete ? null : this.nodes[0]; }
   }
 
   public bool IsIncomplete
   {
      get { return this.nodes.Count != 1; }
   }
 
   public CoalescentTree()
   {
      this.nodes = new List<CoalescentTreeNode>();
   }
 
   internal CoalescentTree(IEnumerable<CoalescentTreeNode> nodes)
   {
      this.nodes = new List<CoalescentTreeNode>(nodes);
   }
 
   public CoalescentTreeNode Coalesce(CoalescentTreeNode left, CoalescentTreeNode right, double time)
   {
      if (!Contains(left) || left.Ancestor != null)
      {
         throw new ArgumentException("left");
      }
      if (!Contains(right) || right.Ancestor != null)
      {
         throw new ArgumentException("right");
      }
      double leftTime = left.GetTotalTime();
      if (leftTime > time)
      {
         throw new ArgumentException("left");
      }
      double rightTime = right.GetTotalTime();
      if (rightTime > time)
      {
         throw new ArgumentException("right");
      }
      
      left.WaitingTime = time - leftTime;
      right.WaitingTime = time - rightTime;
 
      CoalescentTreeNode ancestor = new CoalescentTreeNode()
      {
         DescendantLeft = left,
         DescendantRight = right,
      };
      left.Ancestor = ancestor;
      right.Ancestor = ancestor;
 
      Remove(left);
      Remove(right);
      Add(ancestor);
      return ancestor;
   }
 
   public void Add(CoalescentTreeNode node)
   {
      if (node == null)
      {
         throw new ArgumentNullException("node");
      }
      if (node.Tree != null)
      {
         throw new InvalidOperationException();
      }
      this.nodes.Add(node);
      node.Tree = this;
   }
 
   private bool Remove(CoalescentTreeNode node)
   {
      return this.nodes.Remove(node);
   }
 
   private bool Contains(CoalescentTreeNode node)
   {
      return this.nodes.Contains(node);
   }
 
   internal IList<CoalescentTreeNode> GetNodes()
   {
      return new ReadOnlyCollection<CoalescentTreeNode>(nodes);
   }
}

メンバー変数の nodes は，親が未決定のノードを持たせています。したがってツリーが完成すれば nodes の要素の数は 1 になります。通常の使用において，未完成のツリーは外部から呼ばれないはずなので， nodes の要素数をチェックする， IsIncomplete というプロパティを用意しています。

乱数生成

さて，シミュレーターを実装する前に，乱数生成の拡張メソッドを 1 つ用意しましょう。集団サイズが一定の場合，コアレセントイベントまでの待ち時間は指数分布することが知られています。そこで，指数分布にしたがう乱数を生成する拡張メソッドを定義します。

public static class RandomExtension
{
   public static double NextExponential(this Random random, double rate)
   {
      double r;
      do
      {
         r = random.NextDouble();
      } while (r == 0.0);
      return -Math.Log(r) / rate;
   }
}

ここで rate は指数分布のパラメーターで[B]，期待値は 1 / rate になります。

シミュレーター

いよいよ大詰めです。シミュレーターを実装しましょう。シミュレーションの基本的な流れは次のようになります。

  1. ランダムな 2 サンプルを選択する。
  2. サンプル数に比例した，コアレセントイベントまでの待ち時間を生成する。
  3. 選択した 2 つのサンプルを coalesce する。
  4. サンプルが 1 つになるまで繰り返す。
public class CoalescentSimulator
{
   private readonly Random random;
 
   public CoalescentSimulator()
   {
      this.random = new Random();
   }
 
   public CoalescentTree PerformSimulation(IEnumerable<CoalescentTreeNode> samples)
   {
      CoalescentTree tree = new CoalescentTree(samples);
 
      double time = 0.0;
      while (tree.IsIncomplete)
      {
         IList<CoalescentTreeNode> nodes = tree.GetNodes();
         int count = nodes.Count;
         time = GenerateNextCoalescentEventTime(time, count);
         int leftIndex;
         int rightIndex;
         GetRandomTwoIndexes(count, out leftIndex, out rightIndex);
         tree.Coalesce(nodes[leftIndex], nodes[rightIndex], time);
      }
 
      return tree;
   }
 
   private void GetRandomTwoIndexes(int count, out int index1, out int index2)
   {
      index1 = random.Next(count);
      index2 = random.Next(count - 1);
      if (index2 >= index1)
      {
         ++index2;
      }
   }
 
   private double GenerateNextCoalescentEventTime(double currentTime, int sampleCount)
   {
      double rate = (double)sampleCount * (sampleCount - 1);
      double waitingTime = random.NextExponential(rate);
      return currentTime + waitingTime;
   }
}

GenerateNextCoalescentEventTime メソッド内で指数分布のパラメーター rate に sampleCount * (sampleCount - 1) を与えているところに注意してください。コアレセント理論における時間は世代数単位で表されます。サイズ N の二倍体集団において， n 個の系統がコアレセントイベントを経験する待ち時間の期待値は 4N/n(n-1) 世代なので， GenerateNextCoalescentEventTime メソッドで得られる時間は 4N 世代単位，つまり 4N 倍した値が真の値になります。

使ってみる

ということで，コアレセントシミュレーションにおける最も基本的な場合が実装できました。使用方法は次のようになります。

CoalescentSimulator simulator = new CoalescentSimulator();
var samples = Enumerable.Range(0, 6).Select(i => new CoalescentTreeNode() { Name = (i + 1).ToString() });
CoalescentTree tree = simulator.PerformSimulation(samples);

実際に動かした例を http://ideone.com/t1fkB に貼り付けましたので参考にしてみてください。ツリーの出力方法もそちらを参照。

次回予告

次回は突然変異を実装したいと思います。

更新履歴

  • [2011-05-04 06:50+0900] Coalesce メソッド内で double rightTime = left.GetTotalTime() となっていたのを double rightTime = right.GetTotalTime() に修正。

脚注

  1. Hudson RR (2002) Generating samples under a Wright-Fisher neutral model. Bioinformatics 18: 337–338. [↩]
  2. λ という文字で表現されることが多いです。 [↩]

© Copyright 2011 - 捨てられたブログ

Colorbox Theme byTheme Mansion, Proudly Powred by WordPress