Multilocus genotype frequencies

Published results using this notebook should cite:
Hospital, F., Dillmann, C. and Melchinger, A.E. (1996) A general algorithm to compute multilocus genotype frequencies under various mating systems. Comput. Appl. Biosci, 12(6):455-462.

General Definitions

ngam=2^n;
ngen=ngam(ngam+1)/2;

Imn=IdentityMatrix[ngam];

gam[i_]:=gam[i]=Table[
(Mod[(i-1),2^(n-l+1)]-Mod[(i-1),2^(n-l)])/2^(n-l),{l,1,n}];

Phi[i_]:=Phi[i]=1+gam[i];

zyg[i_,j_]:=zyg[i,j]={gam[i],gam[j]};

gamprod[i_,j_]:=gamprod[i,j]=
Table[zyg[i,j][[ Phi[k][[l]] ]] [[l]],{k,1,ngam},{l,1,n}];

rho[k_,l_]:=rho[k,l]=
Abs[gam[k][[l]]-gam[k][[l+1]]];

g[x_]:=g[x]=1+Sum[2^k x[[n-k]],{k,0,n-1}];

delta[i_,j_]:=If[i==j,0,1]

R[k_]:=R[k]=
1/2*Product[rho[k,l]*r[l]+(1-rho[k,l])*(1-r[l]),{l,1,n-1}];

P[i_,j_]:=P[i,j]=
Sum[R[k]*Imn[[g[gamprod[i,j][[k]]]]],{k,1,ngam}];

gamname[i_]:=StringJoin[Map[ToString,gam[i]]]

If[n>2,
Format[r[x_]]:=Subscripted[r[x]],
Format[r[x_]]:=r];
Format[s[x_]]:=Subscripted[s[x]];
Format[u[x_,y_]]:=Subscripted[u[x,y]];
Format[v[x_]]:=Subscripted[v[x]]

If[n>2,
Format[r[x_],TeXForm]:=Subscripted[r[x]],
Format[r[x_],TeXForm]:=r];
Format[s[x_],TeXForm]:=Subscripted[s[x]];
Format[u[x_,y_],TeXForm]:=Subscripted[u[x,y]];
Format[v[x_],TeXForm]:=Subscripted[v[x]]

rulers={
r[x_]+r[y_]-2r[x_]r[y_]->1-s[x]s[y]-r[x]r[y],
-(r[x_]+r[y_]-2r[x_]r[y_])->-(1-s[x]s[y]-r[x]r[y]),
r[x_]+r[y_]-r[x_]r[y_]->1-s[x]s[y],
-r[x_]-r[y_]+r[x_]r[y_]->s[x]s[y]-1,
1 - r[x_]->s[x],
-1+r[x_]->-s[x]}

Full Sib mating

Pf[i_,j_]:=Pf[i,j]=P[i,j]/.r[x_]->rf[x];
Pm[i_,j_]:=Pm[i,j]=P[i,j]/.r[x_]->rm[x]

Applications

Below are simple implementations of the formulae in Hospital et al. (1996).
More sophisticated implementations taking symmetries into account can be
found in files: 'random','backcross','selfing', and 'fullsib'.

Random mating

f[t_,x_,y_]:=f[t,x,y]=
2^delta[x,y]*
Sum[P[i,j][[x]]*f[t-1,i,j],{i,1,ngam},{j,i,ngam}]*
Sum[P[i,j][[y]]*f[t-1,i,j],{i,1,ngam},{j,i,ngam}]

Backcrossing

B=Union[Map[g,gamprod[b1,b2]]];
nb=Length[B];

f[t_,x_,y_]:=f[t,x,y]=
Sum[
(P[i,B[[j]]][[x]]*P[b1,b2][[y]]+
delta[x,y]*P[i,B[[j]]][[y]]*P[b1,b2][[x]])*
f[t-1,Min[i,B[[j]]],Max[i,B[[j]]]],{i,1,ngam},{j,1,nb}]

Selfing

f[t_,x_,y_]:=f[t,x,y]=
2^delta[x,y]*
Sum[P[i,j][[x]]*P[i,j][[y]]*f[t-1,i,j],{i,1,ngam},{j,i,ngam}]

Full-Sibs mating

G[i_,j_,k_,l_]:=G[i,j,k,l]=
Outer[Times,Pf[i,j],Pm[k,l]]

h[t_,x_,y_,u_,v_]:=h[t,x,y,u,v]=
Sum[
(G[i,j,k,l][[x]][[y]]+delta[x,y]*G[i,j,k,l][[y]][[x]])*
(G[i,j,k,l][[u]][[v]]+delta[u,v]*G[i,j,k,l][[v]][[u]])*
h[t-1,i,j,k,l],
{i,1,ngam},{j,i,ngam},{k,1,ngam},{l,k,ngam}];


f[t_,x_,y_]:=f[t,x,y]=
Sum[
(G[i,j,k,l][[x]][[y]]+delta[x,y]*G[i,j,k,l][[y]][[x]])*
h[t-1,i,j,k,l],
{i,1,ngam},{j,i,ngam},{k,1,ngam},{l,k,ngam}];