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.

Full-Sib Mating

Initialization

ClearAll["Global`*"];

Number of loci

n=3

General definitions

<<"definitions.ma"

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

mat[i_,j_,k_,l_,x_,y_,u_,v_]:=
(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]])

matgen[i_,j_,k_,l_,x_,y_]:=
(G[i,j,k,l][[x]][[y]]+delta[x,y]*G[i,j,k,l][[y]][[x]])

Index for genotypes

S1=
Table[{i,i},{i,1,ngam/2}];
S2=Flatten[
Table[{i,j},
{i,1,ngam/2-1},{j,i+1,ngam-i}],1];
S3=Table[{ngam/2+1-i,ngam/2+i},{i,1,ngam/2}];
S4=Flatten[
Table[{i,j},
{j,ngam/2+2,ngam},{i,ngam+2-j,j-1}],1];
S5=Table[{i,i},{i,ngam/2+1,ngam}];
L1=Length[S1];L2=Length[S2];L3=Length[S3];

genlist=Join[S1,S2,S3,S4,S5]

ngen=Length[genlist]

Clear[index];
index=Table["+",{ngam},{ngam}];
Do[index[[genlist[[i,1]],genlist[[i,2]]]]=i,{i,1,ngen}];

TableForm[index,
TableHeadings->{Table[gamname[i],{i,ngam}],
Table[gamname[i],{i,ngam}]},
TableAlignments->Center]

Index for couples with symmetries :

ulist=Flatten[
Table[Join[genlist[[i]],genlist[[j]]],
{i,1,ngam(ngam+2)/4},{j,i,Max[ngam(ngam+2)/4,ngen+1-i]}],
1];

Table[{genlist[[i]],Table[genlist[[j]],{j,i,Max[ngam(ngam+2)/4,ngen+1-i]}]},
{i,1,ngam(ngam+2)/4}]//ColumnForm;

nsym=Length[ulist]

cp[i_]:=ulist[[i]]

gen[i_]:=genlist[[i]]

Symmetries male/female

sym1[cpI_]:=Flatten[Reverse[Partition[cpI,2]]]

Symmetries on genotypes

sym2[cpI_]:=ngam+1-Reverse[cpI]

a[cpI_,cpX_]:=Apply[mat,Flatten[{cpI,cpX}]]

b[cpI_,genX_]:=Apply[matgen,Flatten[{cpI,genX}]]

Recursion on couples frequencies (with symetries)

h[t_,x_]:=h[t,x]=
Simplify[
Sum[
(a[cp[i],cp[x]] + 
delta[ cp[i] , sym1[ cp[i] ] ] a[ sym1[ cp[i] ] , cp[x] ] +
delta[ sym1[ cp[i] ] , sym2[ cp[i] ] ] * delta[ cp[i] , sym2[ cp[i] ] ] *
( a[ sym2[ cp[i] ] , cp[x] ] +
delta[ cp[i] , sym1[ cp[i] ] ] *
a[ sym1[ sym2[ cp[i] ] ] , cp[x] ] ) ) *
h[t-1, i],
{i,1,nsym}]]

cprec[t_,x_]:=cprec[t,x]=
Sum[
(a[cp[i],cp[x]] + 
delta[ cp[i] , sym1[ cp[i] ] ] a[ sym1[ cp[i] ] , cp[x] ] +
delta[ sym1[ cp[i] ] , sym2[ cp[i] ] ] * delta[ cp[i] , sym2[ cp[i] ] ] *
( a[ sym2[ cp[i] ] , cp[x] ] +
delta[ cp[i] , sym1[ cp[i] ] ] *
a[ sym1[ sym2[ cp[i] ] ] , cp[x] ] ) ) *
z[t-1, i],
{i,1,nsym}]

Format[z[t_,x_]]:=Apply[Subscripted["h"[t]],ulist[[x]]]

Do[h[1,i]=0,{i,nsym}];
h[1,nsym]=1;

htot[t_]:=Table[h[t,i],{i,nsym}]

htot[1]

Recursion on genotype frequencies (with symmetries)

f[t_,x_]:=Simplify[
Sum[
(b[cp[i],gen[x]] + 
delta[ cp[i] , sym1[ cp[i] ] ] b[ sym1[ cp[i] ] , gen[x] ] +
delta[ sym1[ cp[i] ] , sym2[ cp[i] ] ] * delta[ cp[i] , sym2[ cp[i] ] ] *
( b[ sym2[ cp[i] ] , gen[x] ] +
delta[ cp[i] , sym1[ cp[i] ] ] *
b[ sym1[ sym2[ cp[i] ] ] , gen[x] ] ) ) *
h[t-1, i],
{i,1,nsym}]]

freq[t_]:=Table[f[t,x],{x,ngen}]

freq[2]