

BeginPackage["Mathgenome`",
	     {"Graphics`Graphics`",
	     "Statistics`ContinuousDistributions`",
	     "Statistics`DiscreteDistributions`",
		"Statistics`HypothesisTests`",
	      "Statistics`NormalDistribution`",
	      "Statistics`DataManipulation`",
	      "Statistics`DescriptiveStatistics`",
		 "DiscreteMath`Combinatorica`",
	    "Statistics`Common`DistributionsCommon`"}]


HFun::usage=
	"HFun[r] gives genetic distance in morgans for a given recombination
	 fraction r, based on Haldane mapping function"

IHFun::usage=
	"IHFun[d] gives a recombination fraction for a given genetic distance
	 in morgans d, based on Haldane mapping function"

KFun::usage=
	"KFun[r] gives genetic distance in morgans for a given recombination
	 fraction r, based on Kosambi mapping function"

IKFun::usage=
	"IKFun[d] gives a recombination fraction for a given genetic distance
	 in morgans d, based on Kosambi mapping function"

CFun::usage=
	"CFun[r] gives genetic distance in morgans for a given recombination
	 fraction r, based on Carter and Falconer mapping function"

ICFun::usage=
	"ICFun[d] gives a recombination fraction for a given genetic distance
	 in morgans d, based on Carter and Falconer mapping function"

PMFun::usage=
	"PMFun[r] gives genetic distance in morgans for a given recombination
	 fraction r, based on Pascoe and Morton mapping function"

IPMFun::usage=
	"IPMFun[d] gives a recombination fraction for a given genetic distance
	 in morgans d, based on Pascoe and Morton mapping function"


Ssd::usage=
	"Ssd[selection pressure as a fraction] gives the standardized
	 selection differential associated with the given selection 
	 pressure" 

Abo::usage=
	"Abo[# of A, # of B, # of AB, # of O] gives maximum likelihood
	 estimates of the allele frequencies for the ABO system f(A),
	 f(B), f(O) and the support"

HwAbo::usage=
	"HwAbo[# of A, # of B, # of AB, # of O] Neyman-Pearson test for 
	 H-W in ABO system. Gives:
	 {f(A), f(B), f(O), support, {chi square, p}}"

Comp::usage=
	"Comp[{#A1, #B1, #AB1, #O1}, {#A2, #B2, #AB2, #O2}] compares the
	 allele frequencies for ABO system between two populations, by the
	 Neyman-Pearson test. Gives:
	 {chi square, p}"

Iden::usage=
	"Iden[{p1, p2, ..., pn}, {q1, q2, ..., qn}] calculates the genetic
	 identity between two populations with allele frequencies pi and
	 qi respectively, for one locus"

GenDis::usage=
	"GenDis[{p1, p2, ..., pn}, {q1, q2, ..., qn}] calculates the genetic
	 distance between two populations with allele frequencies pi and
	 qi respectively, for one locus"

Nsim::usage=
	"Nsim[mean, standard deviation, n] gives a list of n random numbers
	 from a Normal distribution with the especified mean and standard
	 deviation"

ChiTest::usage=
	"ChiTest[o_,e_,df_,yates_] performs a goodnes of fittness test with
	 Chi square between a list of observed data o, a list of expected
	 data e and with df degrees of freedom. The last argument is either
	 n for no Yates correction, or y for Yates correction. Gives:
	 {Chi square, p}"

Hw::usage=
	"Hw[p] gives Hardy-Weinberg genotype frequencies por the case of
	one locus and two alleles. Gives: {p^2, 2pq, q^2}"

Hwmean::usage=
	"Hwmean[p, genotypic value of AA, genot. val Aa, genot. val aa]
	gives the population mean under Hardy-Weinberg equilibrium"

Avef::usage=
	"Avef[p, genotypic value of AA, genot. val Aa, genot. val aa]
	gives {average effect of A, average effect of a} under 
	Hardy-Weinberg equilibrium"
	 
Eftab::usage=
	"Eftab[p, genotypic value of AA, genot. val Aa, genot. val aa]
	gives a matrix whose rows for {AA, Aa, aa} are:
	Hardy-Weinberg frequencies
	Genotypic values
	Kempthorne-coded genot. values
	Breeding values
	Dominance effects"

Varc::usage=
	"Varc[p, genotypic value of AA, genot. val Aa, genot. val aa]
	gives {genotypic variance, additive variance, dominance variance}
	in a Hardy-Weinberg equilibrium population"

PathAnalysis::usage=
	"PathAnalysis[correlation matrix among basic variables, correlation
	vector between basic variables and response variable] gives
	{vector of path coefficients, error-associated path coefficient,
	table of direct and indirec effects}"


FreqAl::usage=
	"FreqAl[matrix of genotypic absolute frequencies] gives a vector of allele
	frequencies {A1, A2, ..., Ak} from a matrix of genotypic frequencies in a codominant, multiple-
	allele locus. The matrix is symetric of the form {{A1A1, A1A2,...,A1Ak},
	{A2A1, A2A2, ..., A2Ak}, ..., {AkA1, AkA2, ..., AkAk}} and contains the
	observed numbers of genotypes with no distinction between AiAj and AjAi"

HWMultAllele::usage=
	"HWMultAllele[vector of allelic frequencies] gives a matrix of Hardy-Weinberg
	genotypic frequencies from a vector of allele frequencies {A1, A2, ..., Ak}.
	The matrix is symetric of the form {{A1A1, A1A2,...,A1Ak},
	{A2A1, A2A2, ..., A2Ak}, ..., {AkA1, AkA2, ..., AkAk}} and contains the
	genotypic frequencies with no distinction between AiAj and AjAi"


ALC::usage=
	"ALC[model, means, standard errors] is a command for analysis of line crosses, 
	with a design matrix (model), a vector of generation means and their standard
	errors. It gives parameter estimates corresponding to the columns of the design
	matrix, the standard errors of the estimates, the predicted means, and the value 
	of the chi-square test statistic for the fitness of the model"

XAlFreq::usage=
	"XAlFreq[# dominant females, # recessive females, # dominant males, # recessive
	males] gives the array of maximum likelihood estimates {p, q} for the frequencies
	of the dominant (p) and recessive (q) alleles"

Gametes::usage=
	"Gametes[genotype] gives the array of possible gametes 
	{{Ai, Bi, ...},{Aj, Bj, ...}, ...} for the genotype {{Ai, Aj},{Bi, Bj}, ...}.
	Alleles must be coded as string variables"

Crossing::usage=
	"Crossing[genotype1, genotype2] works out the cross of genotype 1 and genotype 2 
	coded as :{{Ai, Aj},{Bi, Bj}, ...}. It gives the array of the progeny genotypes,
	their proportions and probabilities. Alleles must be coded as string variables"

SelecInd::usage=
	"SelecInd[p, id, g, a] calculates a classical selection index. The variable p is a matrix
	nxn of phenotypic covariane, id is an array with 0s and 1s, denoting with 1 those phenotypes
	to be included in the index, g is a matrix nxn of genotypic or additive covariance, a
	is a vector n of economical weights. It gives the coefficients of the index, general
	and individual genetic gains, and relative efficiency"

Cumul::usage=
	"Cumul[list, k] gives the sum of the first k elements of list"

TakeAllele::usage=
	"TakeAllele[alleles, freq] takes a random allele with the specified frequencies freq"


RGenotype::usage=
	"RGenotype[a, p] generates a random genotype from an H-W ppopulation with alleles a and
	allele frequencies p"   

PIC::usage=
	"PIC[p] calculates the Polymorphism Information Conent for an array p of allele frequencies in
	a population" 


HeterMLE::usage=
	"HeterMLE[p] gives the maximum likelihood estimate of heterozygosity for an array p of estimated 
	allele frequencies"



VarHeterMLE::usage=
	"VarHeterMLE[p, n] gives the sampling variance of the MLE estimate of heterozygosity for an array p
	alele frequency estimates and n sampled alleles. A codominant locus is assumed"

HeterUnbiased::usage=
	"HeterUnbiased[p, n] gives an unbiased estimate of heterozygosity for an array p of allele frequency
	estimates and n sampled alleles. A codominant locus is assumed"

VarHeterUnbiased::usage=
	"VarHeterUnbiased[p, n] gives the sampling variance of the unbiased estimate of heterozygosity for 
	an array p of allele frequency estimates and n sampled alleles. A codominant locus is assumed"


NMarker::usage=
	"NMarker[P, d, L] calculates the minimum number of marker loci for an arbitrary locus being within d
	cM from the nearest marker in a a genome of L cM, with a probability P" 

R::usage=
	"R[n1, n2, n3, n4] estimates the relative risk (relative incidence) from a sample, where n1 is the number 
	of patients with a given characteristic (marker), n2 is the number of control individuals with the 
	characteristic, n3 is the number of patients without the characteristic, and n4 is the number of control
	individuals without the characteristic"

RO::usage=
	"RO[n1, n2, n3, n4] estimates the odds ratio or approximate relative risk, where n1 is the number 
	of patients with a given characteristic (marker), n2 is the number of control individuals with the 
	characteristic, n3 is the number of patients without the characteristic, and n4 is the number of control
	individuals without the characteristic"

CorrAssoc::usage=
	"CorrAssoc[n1, n2, n3, n4] estimates a correlation coefficiente for two qualitative traits, where n1 is the number 
	of patients with a given characteristic (marker), n2 is the number of control individuals with the 
	characteristic, n3 is the number of patients without the characteristic, and n4 is the number of control
	individuals without the characteristic"


GamDesT::usage=
	"GamDesT[r, D, t] predicts gametic disequilibrium after t random mating generations in two loci with a recombination
	fraction r, and with an initial gametic (linkage) disequilibrium D"

GamDes::usage=
	"GamDes[AB, Ab, aB, ab] calculates the gametic (linkage) disequilibrium parameter D for two loci with haplotypic
	frequencies AB, Ab, aB and ab"

GamDesSt::usage=
	"GamDesSt[AB, Ab, aB, ab] calculates the normalized gametic (linkage) disequilibrium parameter D' for two loci
	with haplotypic frequencies AB, Ab, aB and ab"

GamDesRoS::usage=
	"GamDesRoS[AB, Ab, aB, ab] calculates the gametic (linkage) disequilibrium parameter rho square for two loci
	with haplotypic frequencies AB, Ab, aB and ab"

GamDesP::usage=
	"GamDesP[AB, Ab, aB, ab, n] calculates the empirical significance level for a chi square test of
	gametic (linkage) disequilibrium, for two loci with haplotypic frequencies AB, Ab, aB and ab, with a sample
	size n"

NClon::usage=
	"NClon[g, v, p] calculates the number of clones of a genome of size g, with a vector that can carry
	inserts of size v, to build a library with a probability p of containing a given point of the target genome"

EDist::usage=
	"EDist[vector] calculates the euclidean distance of the vector from the origin"

EuDist::usage=
	"EuDist[vectora, vectorb] calculates the euclidean distance between vectora and vectorb"

RandomAlleleArray::usage=
	"RandomAlleleArray[n] gives an array of n random allele frequencies with normal distribution and mean
	1/n and variance (1/n)(1-1/n). It may take a while for allele numbers greater than 15"

TRandomAlleleArray::usage=
	"TRandomAlleleArray[n] gives an array of a random number of allele frequencies with normal distribution and mean
	1/n and variance (1/n)(1-1/n). The number has a uniform distribution between 1 and n. It may take a while for
	allele numbers greater than 15"

EPIC::usage=
	"EPIC[p] calculates the Entropy-based Polymorphism Information Content for an array p of allele frequencies"

EntropyD::usage=
	"EntropyD[p]calculates the multivariate Shannon Entropy for a random discrete variable with a vector of
	probabilities p"








Begin["`Private`"]




Ssd[p_]:=(nst=NormalDistribution[0,1];
PDF[nst,z/.FindRoot[(1-CDF[nst,z])==p,{z,0}]]/p)

Abo[a_,b_,h_,o_]:=(
n=a+b+h+o;
ia[fa_,fo_]:=(2a(fa^2/(fa^2+2fa fo))+
a(2fa fo/(fa^2+2fa fo))+h)/(2n);
ib[fb_,fo_]:=(2b(fb^2/(fb^2+2fb fo))+
b(2fb fo/(fb^2+2fb fo))+h)/(2n);
sup[x_,y_]:=a Log[x^2+2x(1-x-y)]+b Log[y^2+2y(1-x-y)]+
h Log[2 x y]+o Log[(1-x-y)^2];
x=.3;y=.2;z=.5;
cia=ia[x,z];
cib=ib[y,z];
While[({x,y}-{cia,cib}).({x,y}-{cia,cib})>10^-20,
x=cia;y=cib;
cia=ia[x,1-x-y];
cib=ib[y,1-x-y]];
{x,y,1-x-y,sup[x,y]}
)

HwAbo[a_,b_,h_,o_]:=(
arr=Abo[a,b,h,o];
somega = Select[{a,b,h,o}, # > 0 & ] . 
Log[Select[{a,b,h,o}, # > 0 & ]/n];
chi=2(somega-sup[x,y])//N;
{arr,{chi,1-CDF[ChiSquareDistribution[1],chi]}}
)

Comp[ar1_,ar2_]:=(
tst=2(Apply[Abo,ar1][[4]]+Apply[Abo,ar2][[4]]-
Apply[Abo,ar1+ar2][[4]]);
{tst,1-CDF[ChiSquareDistribution[2],tst]})

Iden[a_,b_]:=
Sum[a[[i]]b[[i]],{i,Length[a]}]/
Sqrt[Sum[a[[i]]^2,{i,Length[a]}]*
Sum[b[[i]]^2,{i,Length[a]}]]

GenDis[a_,b_]:=-Log[Iden[a,b]]

Nsim[mu_,sigma_,n_]:=
Table[Random[NormalDistribution[mu,sigma]],{n}]

ChiTest[o_,e_,df_,yates_]:=
(If[yates=="y"|| yates=="Y",y=1/2,y=0];
chis=Apply[Plus,(Abs[o-e]-y)^2/e];
{chis,1-CDF[ChiSquareDistribution[df],chis]})

Hw[p_]:={p^2,2p(1-p), (1-p)^2}

Hwmean[p_,d_,h_,r_]:=Hw[p].{d,h,r}

TK[p_,d_,h_,r_]:={d,h,r}-Hwmean[p,d,h,r]

Avef[p_,d_,h_,r_]:=
{p TK[p,d,h,r][[1]]+(1-p) TK[p,d,h,r][[2]],
p TK[p,d,h,r][[2]]+(1-p) TK[p,d,h,r][[3]]}

Eftab[p_,d_,h_,r_]:=(
rv={2 Avef[p,d,h,r][[1]],
Avef[p,d,h,r][[1]]+Avef[p,d,h,r][[2]],
2 Avef[p,d,h,r][[2]]};
{Hw[p],{d,h,r},TK[p,d,h,r],rv,TK[p,d,h,r]-rv})

Varc[p_,d_,h_,r_]:=
{Hw[p].(TK[p,d,h,r]^2),Hw[p].(Eftab[p,d,h,r][[4]]^2),
Hw[p].(Eftab[p,d,h,r][[5]]^2)}

PathAnalysis[a_,r_]:=(b=Inverse[a].r;
{b,Sqrt[1-b.r],a.DiagonalMatrix[b]})

convert[q_]:=Transpose[Table[Table[If[q[[i,j]]==".",".",ToExpression[q
[[i,j]]]],
{i,1,Dimensions[q][[1]]}],{j,1,Dimensions[q][[2]]}]]

ReadConv[file_]:=convert[ReadList[file,Word,RecordLists->True]]

HFun[y_]:=-(1/2)Log[1-2y]

IHFun[m_]:=(1-Exp[-2m])/2

KFun[y_]:=ArcTanh[2y]/2

IKFun[m_]:=Tanh[2m]/2

CFun[y_]:=(1/4)(ArcTanh[2y]+ArcTan[2y])

ficf[y_,m_]:=CFun[y]-m

ICFun[m_]:=
(a=0;b=.499999999999999999999999999999;
While[Abs[a-b]>.000001,
If[ficf[(a+b)/2,m]>0,b=(a+b)/2,a=(a+b)/2]];
(a+b)/2)//N

PMFun[y_]:= (-Log[(1-2y)^2/(1+2y+4y^2)]/12+Sqrt[3]ArcTan[(1+4y)
/Sqrt[3]]/6-0.15115)//N

fipm[y_,m_]:=PMFun[y]-m

IPMFun[m_]:=
(a=0;b=.499999999999999999999999999999;
While[Abs[a-b]>.000001,
If[fipm[(a+b)/2,m]>0,b=(a+b)/2,a=(a+b)/2]];
(a+b)/2)//N

FreqAl[s_] := 
  Table[(s[[i, i]] + 
          Sum[s[[i, j]], {j, 1, Length[s]}])/(2 Sum[
            s[[k, l]], {k, 1, Length[s]}, {l, k, Length[s]}]),
    {i, 1, Length[s]}]

HWMultAllele[s_] := 2Outer[Times, s, s] - DiagonalMatrix[s s]


ALC[md_, means_, se_] := (
    var = DiagonalMatrix[se^2];
    scovar = Inverse[Transpose[md].Inverse[var].md];
    estimates = scovar.Transpose[md].Inverse[var].means;
    sterrors = Table[scovar[[i, i]], {i, 1, Length[scovar]}];
    meanest = md.estimates;
    chicua = Apply[Plus, (means - meanest)^2/se^2];
    {{"Parameter estimates->"}, estimates, {"Standard errors->"}, sterrors, 
    {"Mean estimates->"},meanest,{"Chi Square->"},chicua})


XAlFreq[DF_, dF_, DM_, 
    dM_] := (pxal = (4*dF + 4*DF + 2*dM + 3*DM - 
            Sqrt[16*dF^2 + 16*dF*DF + 16*dF*dM + 8*DF*dM + 4*dM^2 + 8*dF*DM + 
                4*dM*DM + DM^2])/(2.*(2*dF + 2*DF + dM + DM)); 
    qxal = 1 - pxal; {pxal, qxal})

Gametes[genotype_] := (Gam[lst__] := 
      Union[Flatten[Outer[List, lst], Length[{lst}] - 1]];
    Apply[Gam, genotype])

Crossing[a_, 
    b_] := (Cross1[x_, y_] := 
      Flatten[Outer[List, Gametes[a], Gametes[b], 1], 1];
    rawlst = Cross1[a, b];
    corazon = Table[Transpose[rawlst[[i]]], {i, 1, Length[rawlst]}]; 
    corazon = 
      Reverse[Table[
          Reverse[Sort[corazon[[i, j]]]], {i, 1, Length[rawlst]}, {j, 1, 
            Length[a]}]];
    ncor = 
      Reverse[Union[
          Table[Reverse[Sort[corazon[[i, j]]]], {i, 1, Length[rawlst]}, {j, 1,
               Length[a]}]]];
    myfuncount[x_] := Count[corazon, x];
    {ncor, Map[myfuncount, ncor], Map[myfuncount, ncor]/Length[rawlst]}
    )

SelecInd[p_, id_, g_, a_] := (x = Position[id, 
    1]; newp = Transpose[p[[Flatten[x]]]][[Flatten[x]]]; newg = g[[Flatten[
      x]]]; thisb = Inverse[newp].newg.a;
    thisgan = "k"
        Sqrt[thisb.newp.thisb]; indgan = (thisb.newg)/Sqrt[thisb.newp.thisb];
    erel = Sqrt[thisb.newp.thisb]/Sqrt[a.g.a];
    Print["b   ", thisb]; Print["Expected gain    ", thisgan];
    Print["Expected ind. gains     ", "k"indgan];
    Print["Relative efficiency     ", erel])

Cumul[p_, k_] := Apply[Plus, Table[p[[i]], {i, 1, k}]]

TakeAllele[a_, p_] := (rndm = Random[]; 
    landscape = Prepend[Table[Cumul[p, i], {i, 1, Length[p]}], 0]; 
    Do[If[(rndm >= landscape[[i]] && rndm < landscape[[i + 1]]), 
        monallele = a[[i]]], {i, 1, Length[a]}]; monallele)

RGenotype[a_, p_] := {TakeAllele[a, p], TakeAllele[a, p]}

PIC[p_] := 1 - Sum[p[[
  i]]^2, {i, 1, Length[p]}] - Sum[(1/2)*(2p[[i]]p[[j]])^2, {i, 1, Length[p] - 
    1}, {j, i + 1, Length[p]}]

HeterMLE[p_] := 1 - Sum[p[[i]]^2, {i, 1, Length[p]}]


VarHeterMLE[p_, n_] := (2(n - 
    1)/n^3)((3 - 2n)(Sum[p[[i]]^2, {i, 
          1, Length[p]}])^2 + 
              2(n - 2)Sum[p[[i]]^3, {i, 1, Length[
                p]}] + Sum[p[[i]]^2, {i, 1, Length[p]}])


HeterUnbiased[p_, n_] := HeterMLE[p]n/(n - 1)

VarHeterUnbiased[p_, n_] := VarHeterMLE[p, n](n^2)/(n - 1)^2

NMarker[p_, d_, l_] := Log[1 - p]/Log[1 - 2d/l]

R[n1_, n2_, n3_, n4_] := n1(n2 + n4)/(n2(n1 + n3))

RO[n1_, n2_, n3_, n4_] := n1 n4/(n2 n3)

CorrAssoc[n1_, n2_, n3_, n4_] := ((n1 + n2 + n3 + n4) *n1 - (n1 + 
  n2)(n1 + n3))/Sqrt[(n1 + n2)(n3 + n4)(n1 + n3)(n2 + n4)]

GamDesT[r_, de_, t_] := (1 - r)^t de

GamDes[AB_, Ab_, aB_, ab_] := AB ab - Ab aB

GamDesSt[AB_, Ab_, aB_, ab_] := (pofA = AB + Ab; pofB = AB + aB;
   If[GamDes[AB, Ab, aB, ab] < 0, GamDes[AB, Ab, aB, ab]/Min[pofA pofB, (1 - pofA)(1 - pofB)], 
GamDes[AB, Ab, aB, ab]/Min[pofA(1 -pofB), (1 - pofA)pofB]])

GamDesHalf[r_] := IntegerPart[-Log[2]/Log[1 - r]] + 1

GamDesRoS[AB_, Ab_, aB_, ab_] := (pofA = AB + Ab; pofB = AB + aB; 
  GamDes[AB, Ab, aB, ab]^2/(pofA (1 - pofA)pofB(1 - pofB)))

GamDesP[AB_, Ab_, aB_, ab_, n_] := 1 - CDF[ChiSquareDistribution[1],
   GamDesRoS[AB, Ab, aB, ab]n] // N

NClon[g_, v_, p_] := Log[1 - p]/Log[1 - (v/g)]

EDist[x_] := Sqrt[Apply[Plus, x^2]]

EuDist[x_, y_] := EDist[x - y]

NRandAlleles[x_, y_] := Random[Integer, {x, y}] (* a random number of alleles or anything between x and y *)

NorFreq[n_] := Random[NormalDistribution[1/n, Sqrt[(1/n)(1 - 1/n)]]] (*a random normal number with mean 1/n
and variance (1/n)(1-1/n) *)

t[x_] := (x < 0 || x > 1) (* gives True if x is not between 0 and 1 *)

CNorFreq[n_] := (Label[uuu]; myranumb = NorFreq[n]; If[t[myranumb], Goto[uuu], myranumb]) (* a nonnegative random normal 
number with mean 1/nand variance (1/n)(1-1/n) *)

RandomAlleleArray[n_] := (Label[uuu]; firstable = 
  Table[CNorFreq[n], {n - 
    1}]; Clear[myranumb]; 
      firstable = Append[firstable, 1 - Apply[Plus, firstable]]; If[
    MemberQ[firstable, x_ /; x < 0], Goto[uuu], firstable])

TRandomAlleleArray[n_] := RandomAlleleArray[NRandAlleles[1, n]]

EPIC[p_] := 2 - Sum[p[[i]]^2p[[
    j]](3p[[j]] + 4p[[i]]), {i, 
          1, Length[p]}, {j, 1, Length[p]}] + 5Sum[p[[i]]^4, {i, 1, Length[
        p]}]

MyLog2[x_]:=If[x\[Equal]0,0,x Log[2,x]]

EntropyD[p_]:=-Apply[Plus,Map[MyLog2,p]]


Show[Graphics[{{GrayLevel[0], 
        Text[StyleForm["MathGenome", FontSize -> 45, 
            FontFamily -> "Bookman"], {1, 0}]}}, 
    Background -> 
      RGBColor[0, .6, .2]]]; Print[""]; Print["This Mathematica application \
is being developed by"];
Print["Humberto Reyes-Valdés"];
Print["All rights reserved"]

End[]

EndPackage[]