! Combined Linkage and Association analyses, parental genotypes available
! Fulker et al., 1999 ; Posthuma et al., 2004

#define n 3				! number of alleles is 3, coded 1, 2, 3
#define nvar 1				! one phenotype per individual 
#ngroups 5				! 5 Mx groups
#define ndef 2				! 2 covariates (definition variables - age and sex)

G1: calculation group between and within effects

Data Calc 

 Begin matrices;
  A Full 1 n  free			! additive allelic effects within
  C Full 1 n  free			! additive allelic effects between
  D Sdiag n n free 			! dominance deviations within
  F Sdiag n n free 			! dominance deviations between
  I Unit 1 n 				! one's
 end matrices;

Specify A 100 101 102
Specify C 200 201 202
Specify D 800 801 802
Specify F 900 901 902

 Begin algebra;
  K = (A'@I) + (A@I') ;			! Within effects, additive
  L = D + D' ;				! Within effects, dominance
  W = K+L ;				! Within effects - additive and dominance in one matrix

  M = (C'@I) + (C@I') ;			! Between effects, additive
  N = F + F' ;				! Between effects, dominance
  B = M+N ;				! Between effects - additive and dominance in one matrix
 End algebra ; 

st .2 all
end

G2: MZ TWINS
Data NInput=124
 Missing =-99.00              
 Rectangular File=example.dat
 Labels 
   famnr  zygos
   fata1 fata2 mota1 mota2		! fatherallele1 fatherallele2 motherallele1 mothereallele2
   pheno1 age1 sex1 tw1a1 tw1a2		! phenotypetwin1 agetwin1 sextwin1 twin1allele1 twin1allele2
   pheno2 age2 sex2 tw2a1 tw2a2		! phenotypetwin2 agetwin2 sextwin2 twin2allele1 twin2allele2

 !below probability ibd0(z0), ibd1(z1) and ibd2(z2) for positions 0 to 35

   z0_0 z1_0 z2_0 z0_1 z1_1 z2_1 z0_2 z1_2 z2_2 z0_3 z1_3 z2_3 z0_4 z1_4 z2_4 
   z0_5 z1_5 z2_5 z0_6 z1_6 z2_6 z0_7 z1_7 z2_7 z0_8 z1_8 z2_8 z0_9 z1_9 z2_9 
   z0_10 z1_10 z2_10 z0_11 z1_11 z2_11 z0_12 z1_12 z2_12 z0_13 z1_13 z2_13 z0_14 z1_14 z2_14
   z0_15 z1_15 z2_15 z0_16 z1_16 z2_16 z0_17 z1_17 z2_17 z0_18 z1_18 z2_18 z0_19 z1_19 z2_19
   z0_20 z1_20 z2_20 z0_21 z1_21 z2_21 z0_22 z1_22 z2_22 z0_23 z1_23 z2_23 z0_24 z1_24 z2_24
   z0_25 z1_25 z2_25 z0_26 z1_26 z2_26 z0_27 z1_27 z2_27 z0_28 z1_28 z2_28 z0_29 z1_29 z2_29
   z0_30 z1_30 z2_30 z0_31 z1_31 z2_31 z0_32 z1_32 z2_32 z0_33 z1_33 z2_33 z0_34 z1_34 z2_34
   z0_35 z1_35 z2_35

Select if zygos < 3 ; 			! MZ, zygos is coded as 1=MZM, 2=MZF, 3=DZM, 4=DZF, 5=DOS

 Select 	pheno1 pheno2
		fata1 fata2 
		mota1 mota2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2;

  Definition_variables 
		fata1 fata2 
		mota1 mota2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2;

 Begin Matrices;
  X Lower nvar nvar Free 	 	! Additive A
  Z Lower nvar nvar Free  		! Unshared E
  D Lower nvar nvar Free  		! Qtl Q
  
  G Full 1 nvar	Free			! grand means 
  F Full 1 ndef	Fix			! will contain age1 sex1   
  I Full 1 ndef 	Fix		! will contain age2 sex2 	
  R Full nvar ndef Free			! parameter estimates age sex
  
  B Computed n n  = B1 			! spurious and genuine genotypic effects
  W Computed n n  = W1			! genuine genotypic effects
   
  K Full 1 4 Fix			! tw1a1 tw1a2
  L Full 1 4 Fix          		! tw2a1 tw2a2  
  M Full 1 4 Fix			! fata1 mota1 
  N Full 1 4 Fix          		! fata1 mota2
  O Full 1 4 Fix          		! fata2 mota1
  P Full 1 4 Fix          		! fata2 mota2
  S Full 1 1 Fix			! to contain 4 
 End Matrices;

 Matrix S 4					
 Matrix K 1 1 1 1
 Matrix L 1 1 1 1
 Matrix M 1 1 1 1
 Matrix N 1 1 1 1
 Matrix O 1 1 1 1
 Matrix P 1 1 1 1
 	
 Specify K tw1a1 tw1a2 tw1a1 tw1a2
 Specify L tw2a1 tw2a2 tw2a1 tw2a2
 Specify M fata1 mota1 fata1 mota1
 Specify N fata1 mota2 fata1 mota2
 Specify O fata2 mota1 fata2 mota1
 Specify P fata2 mota2 fata2 mota2

 Specify F age1 sex1
 Specify I age2 sex2
 
 Begin Algebra;
  V = (\part(B,M) + \part(B,N) + \part(B,O) + \part(B,P)) % S ; !Between effects
  C = (\part(W,M) + \part(W,N) + \part(W,O) + \part(W,P)) % S ; !for Within effects
  A = X*X';
  E = Z*Z';
  Q = D*D';
 End Algebra;

 Means G+F*R'+V+(\part(W,K)-C)|G+I*R'+V+(\part(W,L)-C);
 Covariance 
	 A+E+Q| A+Q_
         A+Q | A+E+Q ;

st 2 X 2 1 1  
st 1 Z 2 1 1 
st .5 D 2 1 1
st 10 G  2 1 1  
st -.2 R 2 1 1 
st 1.5 R 2 1 2
 Option  nd=4       			! request 4 decimal places in output
 OPtion RS 	 			! request residuals 

End

G3: DZ TWINS
Data NInput=124
 Missing =-99.00              
 Rectangular File=example.dat
 Labels 
   famnr  zygos
   fata1 fata2 mota1 mota2		! fatherallele1 fatherallele2 motherallele1 mothereallele2
   pheno1 age1 sex1 tw1a1 tw1a2		! phenotypetwin1 agetwin1 sextwin1 twin1allele1 twin1allele2
   pheno2 age2 sex2 tw2a1 tw2a2		! phenotypetwin2 agetwin2 sextwin2 twin2allele1 twin2allele2

 !below probability ibd0(z0), ibd1(z1) and ibd2(z2) for positions 0 to 35

   z0_0 z1_0 z2_0 z0_1 z1_1 z2_1 z0_2 z1_2 z2_2 z0_3 z1_3 z2_3 z0_4 z1_4 z2_4 
   z0_5 z1_5 z2_5 z0_6 z1_6 z2_6 z0_7 z1_7 z2_7 z0_8 z1_8 z2_8 z0_9 z1_9 z2_9 
   z0_10 z1_10 z2_10 z0_11 z1_11 z2_11 z0_12 z1_12 z2_12 z0_13 z1_13 z2_13 z0_14 z1_14 z2_14
   z0_15 z1_15 z2_15 z0_16 z1_16 z2_16 z0_17 z1_17 z2_17 z0_18 z1_18 z2_18 z0_19 z1_19 z2_19
   z0_20 z1_20 z2_20 z0_21 z1_21 z2_21 z0_22 z1_22 z2_22 z0_23 z1_23 z2_23 z0_24 z1_24 z2_24
   z0_25 z1_25 z2_25 z0_26 z1_26 z2_26 z0_27 z1_27 z2_27 z0_28 z1_28 z2_28 z0_29 z1_29 z2_29
   z0_30 z1_30 z2_30 z0_31 z1_31 z2_31 z0_32 z1_32 z2_32 z0_33 z1_33 z2_33 z0_34 z1_34 z2_34
   z0_35 z1_35 z2_35

Select if zygos > 2 ; 			! DZ, zygos is coded as 1=MZM, 2=MZF, 3=DZM, 4=DZF, 5=DOS

 Select 	pheno1 pheno2
		fata1 fata2 
		mota1 mota2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2
		z1_20 z2_20;		! this is for position 20 only, but can also be used in batch mode

  Definition_variables 
		fata1 fata2 
		mota1 mota2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2
		z1_20 z2_20;


 Begin Matrices;
  X Lower nvar nvar = X2  		! Additive A
  Z Lower nvar nvar = Z2  		! Unshared E
  D Lower nvar nvar = D2  		! Qtl Q
  
  G Full 1 nvar	    = G2		! grand means 
  F Full 1 ndef				! will contain age1 sex1   
  I Full 1 ndef 			!  will contain age2 sex2 	
  R Full nvar ndef  = R2		! parameter estimates age sex 	
  
  B Computed n n    = B1 		! spurious and genuine genotypic effects
  W Computed n n    = W1		! genuine genotypic effects
   
  K Full 1 4 Fix			! tw1a1 tw1a2
  L Full 1 4 Fix          		! tw2a1 tw2a2  
  M Full 1 4 Fix			! fata1 mota1 
  N Full 1 4 Fix          		! fata1 mota2
  O Full 1 4 Fix          		! fata2 mota1
  P Full 1 4 Fix          		! fata2 mota2
  S Full 1 1 Fix			! to contain 4 

  T Full 1 1 fix			! p_ibd1
  U Full 1 1 fix 			! p_ibd2	
  H Full 1 1 fix			! .5	

 End Matrices;

 Matrix S 4					
 Matrix H 0.5
 Matrix K 1 1 1 1
 Matrix L 1 1 1 1
 Matrix M 1 1 1 1
 Matrix N 1 1 1 1
 Matrix O 1 1 1 1
 Matrix P 1 1 1 1
 	
 Specify K tw1a1 tw1a2 tw1a1 tw1a2
 Specify L tw2a1 tw2a2 tw2a1 tw2a2
 Specify M fata1 mota1 fata1 mota1
 Specify N fata1 mota2 fata1 mota2
 Specify O fata2 mota1 fata2 mota1
 Specify P fata2 mota2 fata2 mota2

 Specify F age1 sex1
 Specify I age2 sex2

 Specify T z1_20
 Specify U z2_20

 
 Begin Algebra;
  V = (\part(B,M) + \part(B,N) + \part(B,O) + \part(B,P)) % S ; !Between effects
  C = (\part(W,M) + \part(W,N) + \part(W,O) + \part(W,P)) % S ; !for Within effects
  A = X*X';
  E = Z*Z';
  Q = D*D';
  J = (H*T) + U ;
 End Algebra;

 Means G+F*R'+V+(\part(W,K)-C)|G+I*R'+V+(\part(W,L)-C);

 Covariance 
	 A+E+Q| H@A+J@Q_
         H@A+J@Q | A+E+Q ;
Options rs
 End

Constrain sum allelic effects = 0 
Constraint ni=1
Begin Matrices; 
  A full 1 n = A1	
  O zero 1 1  
End Matrices;
Begin algebra;
  B = \sum(A) ;
End Algebra;
Constraint O = B ;
end

Constrain sum allelic effects = 0 
Constraint ni=1
Begin Matrices;
 C full 1 n = C1	
 O zero 1 1   
End Matrices;
Begin algebra;
 B = \sum(C) ;
End Algebra;
Constraint O = B ;
option multiple issat 	! this is saturated model for submodel comparison
end


!test for linkage in presence of association 
Drop D 2 1 1 
end

!Test for population stratification: between effects = within effects
Specify 1 A 100 101 102
Specify 1 C 100 101 202			!note: last parameter need not be equal because of constraint group
Specify 1 D 800 801 802
Specify 1 F 800 801 802
end

!Test for presence of dominance 
Specify 1 A 100 101 102
Specify 1 C 100 101 202
Specify 1 D 800 801 802
Specify 1 F 800 801 802
Drop @0 800 801 802
end

!Test for presence of full association
Specify 1 A 100 101 102
Specify 1 C 100 101 202
Specify 1 D 800 801 802
Specify 1 F 800 801 802
Drop @0 800 801 802 100 101	!note: parameters 102 and 202 need not be dropped because of constraint group
end

! test for linkage in absence of association, compare with -2ll from one model above
Free D 2 1 1 
end