! Script to estimate haplotype frequencies ! 3 diallelic loci #ngroups 5 #define nloci 4 #define twopowernloci 16 Title calculate predicted genotypic freqencies Calc Begin Matrices; I full 1 twopowernloci J full 1 twopowernloci K full 1 twopowernloci V full 1 twopowernloci U unit 1 twopowernloci End Matrices; Matrix I 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 Matrix J 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 Matrix J 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 Matrix K 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 Matrix V 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 Begin Algebra; L = I'*I; M = I'*(U-I)+(U-I)'*I; N = (U-I)'*(U-I); O = J'*J; P = J'*(U-J)+(U-J)'*J; Q = (U-J)'*(U-J); R = K'*K; S = K'*(U-K)+(U-K)'*K; T = (U-K)'*(U-K); W = V'*V; X = V'*(U-V)+(U-V)'*V; Y = (U-V)'*(U-V); End Algebra; End Group; Group 2; Estimate haplotype freqs #include chrnb.dat ! variables are allele present (1) or absent (0) at 3 loci Select if nic = 0 ; Select A1 A2 B1 B2 C1 C2 D1 D2 zero ; Definition A1 A2 B1 B2 C1 C2 D1 D2 ; Begin matrices = Group 1; A full 1 1 B full 1 1 C full 1 1 D full 1 1 E full 1 1 F full 1 1 G full 1 1 ! dummy covs (1/2pi) H full 1 twopowernloci free I full 1 1 J full 1 1 Z zero 1 1 ! dummy mean (0) U unit 1 1 ! dummy covs (1) End Matrices; Matrix G .1591550775 Specify A A1 Specify B A2 Specify C B1 Specify D B2 Specify E C1 Specify F C2 Specify I D1 Specify J D2 Start .0625 H 1 1 - H 1 twopowernloci Bound 0 1 H 1 1 - H 1 twopowernloci Begin Algebra; V = ((A.B)@L + ((U-A).B+A.(U-B)) @ M + (U-A).(U-B)@N ). ( (C.D)@O + ((U-C).D+C.(U-D)) @ P + (U-C).(U-D)@Q) . ( (E.F)@R + ((U-E).F+E.(U-F)) @ S + (U-E).(U-F)@T) . ( (I.J)@W + ((U-I).J+I.(U-J)) @ X + (U-I).(U-J)@Y) .(H'*H) ; End Algebra; Means Z; Covariance G; Weight \sum(V); End Group; Group 3; Estimate haplotype freqs second group #include chrnb.dat Select if nic = 1 ; Select A1 A2 B1 B2 C1 C2 D1 D2 zero ; Definition A1 A2 B1 B2 C1 C2 D1 D2 ; Begin matrices = Group 1; A full 1 1 B full 1 1 C full 1 1 D full 1 1 E full 1 1 F full 1 1 G full 1 1 ! dummy covs (1/2pi) H full 1 twopowernloci free I full 1 1 J full 1 1 Z zero 1 1 ! dummy mean (0) U unit 1 1 ! dummy covs (1) End Matrices; Matrix G .1591550775 Specify A A1 Specify B A2 Specify C B1 Specify D B2 Specify E C1 Specify F C2 Specify I D1 Specify J D2 Start .0625 H 1 1 - H 1 twopowernloci Bound 0 1 H 1 1 - H 1 twopowernloci Begin Algebra; V = ((A.B)@L + ((U-A).B+A.(U-B)) @ M + (U-A).(U-B)@N ). ( (C.D)@O + ((U-C).D+C.(U-D)) @ P + (U-C).(U-D)@Q) . ( (E.F)@R + ((U-E).F+E.(U-F)) @ S + (U-E).(U-F)@T) . ( (I.J)@W + ((U-I).J+I.(U-J)) @ X + (U-I).(U-J)@Y) .(H'*H) ; End Algebra; Means Z; Covariance G; Weight \sum(V); End Group; Group 4; Estimate haplotype freqs second group #include chrnb.dat Select if nic = 2 ; Select A1 A2 B1 B2 C1 C2 D1 D2 zero ; Definition A1 A2 B1 B2 C1 C2 D1 D2 ; Begin matrices = Group 1; A full 1 1 B full 1 1 C full 1 1 D full 1 1 E full 1 1 F full 1 1 G full 1 1 ! dummy covs (1/2pi) H full 1 twopowernloci free I full 1 1 J full 1 1 Z zero 1 1 ! dummy mean (0) U unit 1 1 ! dummy covs (1) End Matrices; Start .0625 H 1 1 - H 1 twopowernloci Bound 0 1 H 1 1 - H 1 twopowernloci Matrix G .1591550775 Specify A A1 Specify B A2 Specify C B1 Specify D B2 Specify E C1 Specify F C2 Specify I D1 Specify J D2 Begin Algebra; V = ((A.B)@L + ((U-A).B+A.(U-B)) @ M + (U-A).(U-B)@N ). ( (C.D)@O + ((U-C).D+C.(U-D)) @ P + (U-C).(U-D)@Q) . ( (E.F)@R + ((U-E).F+E.(U-F)) @ S + (U-E).(U-F)@T) . ( (I.J)@W + ((U-I).J+I.(U-J)) @ X + (U-I).(U-J)@Y) .(H'*H) ; End Algebra; Means Z; Covariance G; Weight \sum(V); End Group; Constrain haplotype freqs to sum to 1 Constraint Begin Matrices; H full 1 twopowernloci = H2 I full 1 twopowernloci = H3 J full 1 twopowernloci = H4 X unit 1 3 End Matrices; Constrain X=\sum(h)|\sum(i)|\sum(j); !Fix all Option append Option format=(16(F5.4,1x)) Option mxh=nic4.mat Option mxi=nic4.mat Option mxj=nic4.mat Option rs mu issat Option End Group; !LowND=HighND=Never Save nic4.mxs Specify 3 H 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Specify 4 H 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Start .0625 H 2 1 1 to H 2 1 twopowernloci Start .0625 H 3 1 1 to H 3 1 twopowernloci Start .0625 H 4 1 1 to H 4 1 twopowernloci Option df=-2 ! two fewer constraints End ! lowND=highND get nic4.mxs ! Specify 3 H 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Specify 4 H 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 Start .0625 H 2 1 1 to H 2 1 twopowernloci Start .0625 H 3 1 1 to H 3 1 twopowernloci Start .0625 H 4 1 1 to H 4 1 twopowernloci Option sat=2705.677,6678 Option df=-1 ! one less constraint End