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

「コアレセントシミュレーターを実装する (1)」では突然変異がない場合の標準中立モデルに基づくコアレセントツリーの作成までを実装しました。普通の集団遺伝学の教科書でもそうであるように，次に考慮されるのは突然変異です。

ということで，今回は単一集団・集団サイズ一定・突然変異がある場合でのコアレセントシミュレーションが行えるようになるまでを実装します。 ms では ms nsam nreps -t theta に相当するものです。

突然変異モデル

集団遺伝学では様々な突然変異モデルを扱います。塩基配列レベルでは，無限サイトモデル (infinite-site model) と，マイクロサテライトで主に用いられるステップワイズモデル (stepwise model) がメジャーです。両方を考慮した実装も難しくはありませんが，今回は簡単な無限サイトモデルだけを考慮します。

遺伝子座

遺伝子座は，ここでは無限長の塩基配列です。突然変異を受けたサイトを内部に保持するコレクションと言い換えることができます。塩基サイトを 0 以上 1 未満の実数として，それをメンバーのコレクションに格納します。

public class Locus : IEnumerable<double>
{
   private readonly HashSet<double> mutations;
  
   public Locus()
   {
      this.mutations = new HashSet<double>();
   }
  
   public Locus(IEnumerable<double> mutations)
   {
      this.mutations = new HashSet<double>(mutations);
   }
  
   public void AddMutation(double site)
   {
      this.mutations.Add(site);
   }
  
   public bool IsMutant(double site)
   {
      return this.mutations.Contains(site);
   }
  
   public IEnumerator<double> GetEnumerator()
   {
      return this.mutations.GetEnumerator();
   }
   System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
   {
      return GetEnumerator();
   }
}

厳密には AddMutation メソッドで， Collection<T>.Contains を使って同じサイトに突然変異を発生させようとしているかをチェックする必要がありますが，ここでは無視しています。

遺伝子座はサンプル，つまりコアレセントツリーのノードが持っているはずです。ということで CoalescentTreeNode クラスのプロパティとして Locus を加えてやります。

private Locus locus;
public Locus Locus
{
   get
   {
      if (this.locus == null)
      {
         this.locus = new Locus();
      }
      return this.locus;
   }
   internal set { this.locus = value; }
}

遺伝

先祖で起きた突然変異は後代に遺伝します。したがってあるノードに対して突然変異を与えてやる場合は，子孫ノードにも同じサイトに変異を起こしてやって整合性を合わせてやる必要があります。

public void AddMutation(double site)
{
   foreach (var node in EnumerateAll(this))
   {
      node.Locus.AddMutation(site);
   }
}
 
internal static IEnumerable<CoalescentTreeNode> EnumerateAll(CoalescentTreeNode baseNode)
{
   if (baseNode == null)
   {
      throw new ArgumentNullException("baseNode");
   }
  
   yield return baseNode;
   if (baseNode.DescendantLeft != null)
   {
      foreach (var descendant in EnumerateAll(baseNode.DescendantLeft))
      {
         yield return descendant;
      }
   }
   if (baseNode.DescendantRight != null)
   {
      foreach (var descendant in EnumerateAll(baseNode.DescendantRight))
      {
         yield return descendant;
      }
   }
}

祖先で起きた突然変異は子孫に遺伝します。コアレセント理論では逆方向に考えることに注意すると，コアレセントイベント時に子孫同士の持つ突然変異を共通祖先に持たせてやる必要があります。

まずは Locus をマージさせるメソッドを定義します。ここでは演算子のオーバーロードとしています。

public static Locus operator |(Locus c1, Locus c2)
{
   if (c1 == null)
   {
      throw new ArgumentNullException("c1");
   }
   if (c2 == null)
   {
      throw new ArgumentNullException("c2");
   }
   return new Locus(c1.mutations.Union(c2.mutations));
}

コアレセントイベント時に共通祖先ノードを作る際に，これを利用します。 CoalescentTree.Coalesce メソッド中で次の修正を施します。

CoalescentTreeNode ancestor = new CoalescentTreeNode()
{
   DescendantLeft = left,
   DescendantRight = right,
   Locus = left.Locus | right.Locus,
};

乱数生成

シミュレーターに突然変異の機能を追加する前に，また 1 つ乱数生成の拡張メソッドを作ります。突然変異の個数はコアレセントイベントまでの待ち時間の長さを期待値とするようなポアソン分布にしたがいます。なのでポアソン分布にしたがう乱数を生成するメソッドを定義します。

public static int NextPoisson(this Random random, double mean)
{
   double lambda = Math.Exp(mean);
   int k = 0;
   if (!double.IsInfinity(lambda))
   {
      while ( (lambda *= random.NextDouble()) > 1.0)
      {
         k++;
      }
   }
   else
   {
      lambda = mean;
      while ( (lambda += Math.Log(random.NextDouble())) > 0.0)
      {
         k++;
      }
   }
   return k; 
}

期待値が大きいと計算機の都合で e の期待値乗が無限大に発散してしまうので，そのような場合は対数計算で処理しています^[A]。

シミュレーターの修正

シミュレーターに突然変異機能を実装する準備ができました。

public CoalescentTree PerformSimulation(IEnumerable<CoalescentTreeNode> samples, double theta = 0.0)
{
   if (theta < 0)
   {
      throw new ArgumentOutOfRangeException("theta");
   }
   CoalescentTree tree = new CoalescentTree(samples);
  
   double time = 0.0;
   while (tree.IsImcomplete)
   {
      IList<CoalescentTreeNode> nodes = tree.GetNodes();
      int count = nodes.Count;
      time = GenerateNextCoalescentEventTime(time, count);
      int leftIndex;
      int rightIndex;
      GetRandomTwoIndexes(count, out leftIndex, out rightIndex);
      CoalescentTreeNode leftNode = nodes[leftIndex];
      CoalescentTreeNode rightNode = nodes[rightIndex];
      tree.Coalesce(leftNode, rightNode, time);
      AddMutation(leftNode, theta);
      AddMutation(rightNode, theta);
   }
 
   return tree;
}
 
private void AddMutation(CoalescentTreeNode node, double theta)
{
   int count = random.NextPoisson(node.WaitingTime * theta);
   while (--count >= 0)
   {
      double site = random.NextDouble();
      node.AddMutation(site);
   }
}

theta は集団遺伝学の世界ではお馴染みの θ = 4Nu (4 × 集団サイズ × 突然変異率) です。前回の記事でも述べましたが， node.WaitingTime は 4N 世代単位時間ですので， node.WaitingTime * theta で 4N が消去されています。

使ってみる

前回と同様に使うことができます。変化したのは theta パラメーターが増えただけです。

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

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

次回予告

次回は集団のサイズ変化や分化といったデモグラフィーの影響を実装したいと思います。

脚注

余談ですが，私は期待値が大きい場合をテストしていなかったため無限ループにかなりハマりました。当たり前のことですが，テストはちゃんと定義域をある程度網羅するように作らないとダメですね。ちなみに無限ループに陥っているときは Visual Studio のデバッグ機能で現在のステートメントを表示することができます。このバグのおかげで知りました。 [↩]

捨てられたブログ