(*****************************************************************************
Mathematica code for the paper titled "The Statistical Analysis of
  Spatially Clustered Genes under the Maximum Gap Criterion,"
  R. Hoberman, D. Sankoff and D. Durand. JCB 2005. This file contains the
  Mathematica code to compute Equations 1 through 8 in the above
  paper.  This code computes the probabilities of max-gap clusters for
  a marked-gene and a whole-genome comparison scenario.  No
  optimizations of running times have been made.  The code is simply a
  straightforward translation of the equations as given in the paper.
*****************************************************************************)


<<Binomial.m                    (* Load package to compuate FastBinomial *)
Off[General::spell];            (* Turn off spelling warnings *)
Off[General::spell1];

(***************************************************************************** 
                        Input Parameters to Functions
 *****************************************************************************
n: number of genes in the genome
m: number of marked genes (Marked Scenario)
   number of homologs shared between the two genomes (Whole Genome Comparison)
h: cluster size (number of genes in the cluster)
g: maximum allowed gap size between marked genes or homologs 
*)

(***************************************************************************** 
PosBin: n choose m
When m < n this function returns 0 instead of the standard negative value.
It requires FastBinomial (defined in Binomial.m) for efficiency in computing
large binomials. In addition, this function keeps a table of
previously computed values, so each binomial is only computed once.
If lack of memory becomes an issue then comment out "PosBin[n,m] ="
******************************************************************************)
PosBin[n_,m_] := PosBin[n,m] =
  Which[
    n < m, 0,
    n == m, 1,
    n > m, FastBinomial[n,m]
  ];

(***************************************************************************** 
Equation 1: d_e where e=2
Number of ways to generate a cluster of size m and length l.
******************************************************************************)   
d2[m_, g_, l_] := Sum[(-1)^i * PosBin[m-1, i] * PosBin[l-i*(g+1)-2,m-2], 
{i, 0, Floor[(l-m)/(g+1)]}];

(***************************************************************************** 
Equation 1: d_e where e=1
Number of ways to generate a cluster of size m and length <= l.
When l is greater than or equal to the maximum cluster length, a simpler
expression is available.
*****************************************************************************)   
d1[m_, g_, l_] :=
  If[l>=(m+(m-1)g),
    (g+1)^(m-1),
     Sum[(-1)^i * PosBin[m-1, i] * PosBin[l-i*(g+1)-1,m-1],
       {i, 0, Floor[(l-m)/(g+1)]}]
  ];

(***************************************************************************** 
Equation 1: d_e where e=0
Number of ways to fill a window of length l with m genes such that no
gap is larger than g.
******************************************************************************)   
d0[m_, g_, l_] := Sum[(-1)^i * PosBin[m-1, i] * PosBin[l-i*(g+1),m], 
{i, 0, Floor[(l-m)/(g+1)]}];

(***************************************************************************** 
Equation 2: Complete Cluster
Probability of observing all m marked genes such that the gap between adjacent
marked genes does not exceed g.
******************************************************************************)   
P[n_, m_, g_] := 
 If[(m-1)*g+m <= n+1,	(* maximum cluster length *)
   Return[N[(n - m + 1 - (m-1)g/2) * (g+1)^(m-1)/PosBin[n,m]]],
   Return[N[d0[m,g,n]/PosBin[n,m]]]
 ];


(***************************************************************************** 
Equation 3: Complete Cluster with Length Restriction
Probability of observing all m marked genes such that the gap between adjacent
marked genes does not exceed g, and the length of the cluster does not exceed r.
******************************************************************************)   
R[n_, m_, g_, r_] := N[(Max[0,n-r+1]*d1[m,g,r] + d0[m,g,r-1])/PosBin[n,m]];


(***************************************************************************** 
*    Equation 4: Incomplete Cluster of size at least h                       *
*    			                                                     *
*    The input arguments are maxN and maxM so that n and m can be used for   *
*    the matrix indices in building up the dynamic programming table.        *
*						                             * 
*    The probability of finding an incomplete cluster of size at least h     *
*    when there are m special genes in a genome of size n.  When h > m/2     *
*    the function simply calls Qhalf (Equation 4, defined below).	     *
*    Otherwise, it implements a dynamic program to enumerate all ways of     *
*    placing the m marked genes in a genome of size n in such a way that a   *
*    cluster of size h or larger is not created.  The probability of	     *
*    observing at least one cluster is then just 1 - the probability of	     *
*    observing no clusters.						     *
*    									     *
*    Notice that the dynamic program is not recursive, but implemented by    *
*    filling successive levels of the matrix.  To save memory, the code	     *
*    only remembers the values of the matrix for the current and previous    *
*    value of n (eta and prevEta).  Also, Mathematica's matrix indices start *
*    at one not zero so all the indices are one greater than the values	     *
*    shown in the paper. So eta[m,j,c] in the code is actually equivalent    *
*    to eta[n,m-1,j-1,c-1] as defined in the paper, and prevEta[m,j,c] is    *
*    equivalent to Eta[n-1,m-1,j-1,c-1]. 			             *
******************************************************************************)
Q[maxN_, maxM_, h_, g_] := Module[ {eta, prevEta, n,m,j, c, incJ},
  If[h > maxM/2,
    Return[Qhalf[maxN,maxM,h,g]];
  ];
  eta = Array[0,{maxM+1,g+2,h+1}];
  prevEta = Array[0,{maxM+1,g+2,h+1}];
  For[n=0, n <= maxN, n++,
    For[m=0, m <= maxM, m++,
      For[j=0, j <= g+1, j++,
        For[c=0, c <= h, c++,
          incJ = Min[j+1,g+1];
          eta[[m+1,j+1,c+1]] = Which[
            c == h, 0,
            n < m,  0,
            m == 0, 1,
            j <= g, prevEta[[m+1,incJ+1,c+1]] + prevEta[[m,1,c+2]],
            True,   prevEta[[m+1,incJ+1,c+1]] + prevEta[[m,1,2]]
          ];
        ];
      ];
    ];
    prevEta = eta;
  ];
  Return[N[1 - (eta[[maxM+1,g+2,1]]/PosBin[maxN,maxM])]];
];

(***************************************************************************** 
Equation 5: Incomplete Cluster when h > m/2
The probability of observing a cluster of at least size h, when h > m/2.  
If h <= m/2 this function will not return the correct probability.
Called by Q[n,m,h,g].
*****************************************************************************)
Qhalf[n_,m_,h_,g_] := Module[{num},
  If[h>m,Return[0]];
  num = 0;
  For[l=h,l <= Min[n,h+g(h-1)],l++,
    num += d2[h,g,l] * 
     ((Max[0,n-l-g] * PosBin[n-l-g-1,m-h]) + PosBin[n-l+1,m-h+1] - PosBin[n-l-g,m-h+1]);
  ];
  Return[N[num/PosBin[n,m]]];
];

(***************************************************************************** 
Whole Genome Comparison, complete cluster
Probability of all m homologs clustering together in both genomes.
*****************************************************************************)
P2[n_,m_,g_] := P[n,m,g]^2;

(*****************************************************************************
Equation 6: Whole genome comparison, upper bound on the probability of
observing an incomplete cluster of size exactly h.  To save some
computation, q1 and q2 (Equation 6) are computed simulataneously.
*****************************************************************************)
Pup[n_,m_,h_,g_] := Module[
{ans,t,L,q1,q2,l,b,q1L,q2L,c,d2val},
   ans = 0;
   For[t=0,t<=m-h,t++,
     L = Min[n,(h-1)g + h];
     q1 = 0; 
     q2 = 0;
     For[l=h,l<=L,l++,
       b = l+2(g+1);
       q1L = 0;  
       q2L = 0;
       If[n-b+1 > 0, 
         q1L = (n-b+1) * PosBin[b-h,t] * PosBin[n-b,  m-h-t];
         q2L = (n-b+1) * PosBin[n-b,t] * PosBin[n-h-t,m-h-t];
       ];
       For[j=0,j<=Min[g,n-l],j++,
        c = Min[n,l+j+g+1];
        q1L += If[n-l-j>g,2,1] * PosBin[c-h,t] * PosBin[n-c,  m-h-t];
        q2L += If[n-l-j>g,2,1] * PosBin[n-c,t] * PosBin[n-h-t,m-h-t];	
       ];
       d2val = d2[h,g,l];
       q1 += d2val * q1L;
       q2 += d2val * q2L;
     ];
     ans += Factorial[t] * Factorial[m-h-t] * q1 * q2;
   ];
   Return[N[ans * Factorial[h] / (Factorial[m] * (PosBin[n,m])^2)]];
];


(*****************************************************************************
Equation 8: Whole genome comparison, lower bound on the probability of
observing an incomplete cluster of size exactly h.  To save some
computation, q2,q3, and q4 (Equation 6) are computed together.
*****************************************************************************)
Plow[n_,m_,h_,g_] := Module[
{ans,t,L,q2,q3,q4,l,b,q2L,q3L,q4L,c,d2val},
   ans = 0;
   For[t=0,t<=m-h,t++,
     L = Min[n,(h-1)g + h];
     q2 = 0;
     q3 = 0;
     q4 = 0; 
     For[l=h,l<=L,l++,
       b = l+2(g+1);
       q2L = 0;
       q3L = 0;  
       q4L = 0;  
       If[n-b+1 > 0, 
         q2L = (n-b+1) * PosBin[n-b,t] * PosBin[n-h-t,m-h-t];
         q3L = (n-b+1) * PosBin[l-h,t] * PosBin[n-b,  m-h-t];
         q4L = (n-b+1) * PosBin[n-b,t] * PosBin[n-h-t-b+l, m-h-t];
       ];
       For[j=0,j<=Min[g,n-l],j++,
        c = Min[n,l+j+g+1];
        q2L += If[n-l-j>g,2,1] * PosBin[n-c,t] * PosBin[n-h-t,m-h-t];	
        q3L += If[n-l-j>g,2,1] * PosBin[l-h,t] * PosBin[n-c,  m-h-t];
        q4L += If[n-l-j>g,2,1] * PosBin[n-c,t] * PosBin[n-h-t-c+l, m-h-t];
       ];
       d2val = d2[h,g,l];
       q2 += d2val * q2L;
       q3 += d2val * q3L;
       q4 += d2val * q4L;
     ];
     ans += Factorial[t] * Factorial[m-h-t] * q3 * (2*q2 - q4);
   ];
   Return[N[ans * Factorial[h] / (Factorial[m] * (PosBin[n,m])^2)]];
];