(*^ ::[ frontEndVersion = "Microsoft Windows Mathematica Notebook Front End Version 2.2"; microsoftWindowsStandardFontEncoding; fontset = title, "Times New Roman", 24, L0, center, nohscroll, bold; fontset = subtitle, "Times New Roman", 18, L0, center, nohscroll, bold; fontset = subsubtitle, "Times New Roman", 14, L0, center, nohscroll, bold; fontset = section, "Times New Roman", 14, L0, bold, grayBox; fontset = subsection, "Times New Roman", 12, L0, bold, blackBox; fontset = subsubsection, "Times New Roman", 10, L0, bold, whiteBox; fontset = text, "Times New Roman", 12, L0; fontset = smalltext, "Times New Roman", 10, L0; fontset = input, "Courier New", 12, L0, nowordwrap, bold; fontset = output, "Courier New", 12, L0, nowordwrap; fontset = message, "Courier New", 10, L0, nowordwrap, R65535; fontset = print, "Courier New", 10, L0, nowordwrap; fontset = info, "Courier New", 10, L0, nowordwrap; fontset = postscript, "Courier New", 8, L0, nowordwrap; fontset = name, "Times New Roman", 10, L0, nohscroll, italic, B65535; fontset = header, "Times New Roman", 10, L0, right, nohscroll; fontset = footer, "Times New Roman", 10, L0, right, nohscroll; fontset = help, "Times New Roman", 10, L0, nohscroll; fontset = clipboard, "Times New Roman", 12, L0, nohscroll; fontset = completions, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = graphics, "Courier New", 10, L0, nowordwrap, nohscroll; fontset = special1, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = special2, "Times New Roman", 12, L0, center, nowordwrap, nohscroll; fontset = special3, "Times New Roman", 12, L0, right, nowordwrap, nohscroll; fontset = special4, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = special5, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = leftheader, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = leftfooter, "Times New Roman", 12, L0, nowordwrap, nohscroll; fontset = reserved1, "Courier New", 10, L0, nowordwrap, nohscroll;] :[font = input; initialization; nowordwrap; ] *) chiasma:=(dial:=(question= InputString["Do you want to repeat the calculations with different guesses(y/n)?"]; If[question=="y" || question=="Y", (nguess=Table[Input[{"Your guess for segment", S[[i]]}],{i,Length[S]}]; mapa2[S,m,n,nguess ]),(Print[MatrixForm[{{"Option","Type"}, {"1D likelihood plots","plots Shift-Return"}, {"1D color likelihood plots","plotsc Shift-Return"}, {"3D likelihood plots","plot3d Shift-Return"}, {"3D color likelihood plots","plot3dc Shift-Return"}, {"Critical point test","test Shift-Return"}, {"Plots of mapping functions","mfun Shift-Return"}}]])]); mapa2[S_,m_,n_,guess_]:=( ta[x_]:=Complement[S,x]; (*complement of set x*) tab=Map[ta,m,{2}]; (*table of comb. of seg. with no chiasma*) l[x_]:=1-x; tab=Map[l,tab,{2}]; (*same table as 1-A...*) a=Table[Apply[Times, Join[tab[[i,j]],m[[i,j]]]],{i,1,Length[m]}, {j,1,Length[m[[i]]]}]; (*table of expected freq. of comb.*) b=Table[Apply[Plus,a[[i]]],{i,1,Length[m]}]; (*table of exp. val. of config.*) c=Table[a[[i,j]]/b[[i]],{i,1,Length[m]}, {j,1,Length[m[[i]]]}]; (*conditional expected freq. of combinations*) c=Table[c[[i]]n[[i]],{i,1,Length[n]}]; (*conditional expected numbers of combinations*) d=Table[ Table[If[MemberQ[m[[i,j]],S[[k]]],1,0], (*checks if a given segment is in a given combination*) {i,1,Length[m]} ,{j,1,Length[m[[i]]]}],{k,1,Length[S]}]; (*gives a list of lists like m. It has 1 if chiasma is present at the given segment, 0 if not*) f=Table[Flatten[c] Flatten[d[[i]]],{i,1,Length[S]}]; (*list of lists. every list contains expected # if chiasma is present in segment, 0 if not. the sum of each primary list is the count for the segment*) g=guess; LL1=N[n.Log[b]/.Table[S[[j]]->g[[j]],{j,1,Length[S]}],8]; gan=N[n.Log[b]/.Table[S[[j]]->g[[j]],{j,1,Length[S]}],10]; liston={gan}; prueba=1; contador=0; Print["Iterations and log-likelihoods"]; Print[""]; Print["Iteration ",0," Support ",gan]; While[prueba>10^(-25) && contador<300,contador=contador+1; ga=gan; g=Table[Apply[Plus,f[[i]]/.Table[S[[j]]->g[[j]] ,{j,1,Length[S]}]],{i,1,Length[S]}]/ Apply[Plus,n]; gan=N[n.Log[b]/.Table[S[[j]]->g[[j]],{j,1,Length[S]}],16]; Print["Iteration ",contador," Support ",gan]; liston=Append[liston,gan]; prueba=gan-ga]; (*iteration, first term is parameter, second term is likelihood*) LL2=N[n.Log[b]/.Table[S[[j]]->g[[j]],{j,1,Length[S]}],16]; expected=Apply[Plus,n] b/.Table[S[[i]]->g[[i]], {i,1,Length[S]}]; (*expected #s of configurations*) hal=-(1/2)Log[1-g]100//N; (*map distances, haldane*) kos=ArcTanh[g]50//N; (*map distances, kosambi*) carfal=25(ArcTanh[g]+ArcTan[g])//N; (*map distances, carter and falconer*) d1=Table[D[expected. Log[b],S[[i]]],{i,1,Length[S]}]; (*derivating the log-lik, evaluated in exp. #s*) d2=Table[D[d1,S[[i]]],{i,1,Length[S]}]; (*Hessian matrix*) var=-Inverse[d2/.Table[S[[j]]->g[[j]],{j,1,Length[S]}]]; (*inv. of expected Hessian*) stdev=Sqrt[Table[var[[i,i]],{i,1,Length[S]}]]; (*square roots of diag. elements*) chis=Apply[Plus,((n-expected)^2)/expected]; (*chi square value*) df=Length[n]-Length[S]-1; lf=Select[n,#>0&]. (*selecting nonzero elements..*) Log[(Select[n,#>0&]/Apply[Plus,n])]; (*to calculate log-lik full model*) np=2(lf-n.Log[b]/.Table[S[[j]]->g[[j]],{j,1,Length[S]}])//N; Print[contador," iterations"]; Print["*************"]; Print["Seed values"]; Print[MatrixForm[{Prepend[S,"Segments"], Prepend[guess,"Init.values"]}]]; Print[""]; Print["Log-likelihood = ",LL1]; Print["*************"];Print[""]; Print["Maximum likelihood estimates"]; Print[""]; Print[MatrixForm[{Prepend[S,"Segment"], Prepend[g//N,"P(chiasma)"], Prepend[stdev//N,"S.D."], Prepend[hal//N,"Haldane (cM)"], Prepend[kos//N,"Kosambi (cM)"], Prepend[carfal//N,"Cart.-Falc. (cM)"]}]]; Print["Log-likelihood = ",LL2]; Print["*************"];Print[""]; Print["Observed and expected frequencies of meiotic configurations"]; Print[""]; Print[MatrixForm[Transpose[{Prepend[n,"Observed"], Prepend[expected//N,"Expected"]}]]]; Print["*************"];Print[""]; Print["Chi Squared ",df," df"," = ",chis//N]; Print["Neyman Pearson ",df," df"," = ",np//N]; Print["*************"]; dial ); pregunta=InputString["Do you have your data in a file (y/n)"]; If[pregunta=="n" || pregunta=="N", (S=Input["Set of segments"]; size=Input["Number of classes"]; m=Table[Input[{"Array of chiasmate sets for class",i}],{i,size}]; n=Table[Input[{"Nuber of observations for class",i}],{i,size}]; guess=Table[Input[{"Your guess for segment",S[[i]]} ],{i,Length[S]}]; cuestion=InputString["Do you want to save your data(y/n)?"]; If[cuestion=="y" || cuestion=="Y", Save[InputString["Name of your file"],S,m,n,guess]]), (archivo=InputString["File name"];Get[archivo])]; mapa2[S,m,n,guess]) (* :[font = input; initialization; nowordwrap; ] *) llk:=n.Log[b] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) lik[k_]:=n.Log[b]/.Table[S[[j]]-> ReplacePart[g,x,k][[j]],{j,1,Length[S]}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) ploli[k_]:=Plot[lik[k],{x,.001,.999}, PlotPoints->200,PlotLabel->FontForm[S[[k]], {"Times-Bold",13}], AxesLabel->{"P(chiasma)","Log-likelihood"}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plolic[k_]:=Plot[lik[k],{x,.001,.999}, PlotPoints->200,PlotLabel->FontForm[S[[k]], {"Times-Bold",13}], AxesLabel->{"P(chiasma)","Log-likelihood"}, PlotStyle->Hue[.65]] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plots:=Do[ploli[k],{k,1,Length[S]}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plotsc:=Do[plolic[k],{k,1,Length[S]}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) test:=(likelihood=n.Log[b]; Print["Support function"]; Print[""]; Print[likelihood]; fd=Table[D[likelihood,S[[i]]],{i,1,Length[S]}]; Print[""]; Print["First derivatives"]; Print[""]; Print[MatrixForm[fd/.Table[S[[j]]->g[[j]],{j,1,Length[S]}]]]; hessian=Table[ D[fd,S[[i]]],{i,1,Length[S]}]; hessian2=hessian/.Table[S[[j]]->g[[j]],{j,1,Length[S]}]; Print[""]; Print["Hessian Matrix"]; Print[""]; Print[MatrixForm[hessian2//N]]; Print[""]; eig=Eigenvalues[N[hessian2]]; Print["Eigenvalues"]; Print[""]; Print[MatrixForm[eig]]; Print[""]; Print["Note: When all your estimates are different form zero and one,"]; Print["a nonegative eigenvalue indicates that the critical point is not a maximum"] ) (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) lk[k_,l_]:=( lis1=ReplacePart[g,x,k]; lis2=ReplacePart[lis1,y,l]; n.Log[b]/.Table[S[[j]]->lis2[[j]],{j,1,Length[S]}]) (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plk[k_,l_]:=Plot3D[lk[k,l],{x,0.000001,.999999}, {y,0.000001,.99999}, AxesLabel->{S[[k]],S[[l]],""}, PlotLabel->"Log-likelihood", DefaultFont->{"Times",10}, ColorFunction->GrayLevel] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plkc[k_,l_]:=Plot3D[lk[k,l],{x,0.000001,.999999}, {y,0.000001,.99999}, AxesLabel->{S[[k]],S[[l]],""}, PlotLabel->"Log-likelihood", ColorFunction->Hue,DefaultFont->{"Times",10}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plot3d:= Do[Do[plk[i,j],{i,j+1,Length[S]}],{j,1,Length[S]-1}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) plot3dc:= Do[Do[plkc[i,j],{i,j+1,Length[S]}],{j,1,Length[S]-1}] (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) mfun:=( H[y_]:=-(1/2)Log[1-y]; K[y_]:=ArcTanh[y]/2; F[y_]:=(1/4)(ArcTanh[y]+ArcTan[y]); ploty=Plot[{100H[y],100K[y],100F[y]}, {y,0,.99}, PlotPoints->500,PlotRange->All, PlotStyle->{{Hue[.5]},{Hue[.2]},{Hue[.38]}}, Background->GrayLevel[0], AxesLabel->{"P(chiasmate)","Map Distance"}, DisplayFunction->Identity]; haldy:=Graphics[{{Hue[.5],Text [FontForm["Haldane",{"Courier-Bold",8}], {.4,220}]}}]; kosy:=Graphics[{{Hue[.2],Text [FontForm["Kosambi",{"Courier-Bold",8}], {.4,160}]}}]; cafy:=Graphics[{{Hue[.38],Text [FontForm["Cart-Falc",{"Courier-Bold",8}], {.4,100}]}}]; Show[ploty ,haldy,kosy,cafy,DisplayFunction->$DisplayFunction]) (* :[font = input; initialization; nowordwrap; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 65535; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; bold; fontName = "Courier New"; fontSize = 12; ] *) (Show[Graphics[{{Hue[.1],Text[FontForm["WELCOME!",{"Courier-Bold",20}], {1,0}]}},Background->Hue[.8,.8,.9]] ];Print[""];Print["This Mathematica application was developed by"]; Print["M. Humberto Reyes-Valdés and David M. Stelly"]; Print["Texas A&M University, Department of Soil & Crop Sciences"]; Print[""]; Print["To start the program, write: chiasma Shift-Return"]) (* ^*)