## Detailed description of the data set is available in the online GEMMA user manual
## Each of the following steps may take over one minute to run






## To calculate a centered relatedness matrix:
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -a mouse_hs1940.anno.txt -gk -o mouse_hs1940

# The estimated relatedness matrix should look like this:
0.3350590  -0.0227226  0.0103535 ...
-0.0227226  0.3035960 -0.0253762 ...
0.0103535  -0.0253762  0.3536100 ...
....................................







## To perform association tests with a linear mixed model:
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 1 -a mouse_hs1940.anno.txt -k ./output/mouse_hs1940.cXX.txt -lmm -o mouse_hs1940_CD8_lmm

# The result for top 5 SNPs should look like this:
chr	rs	ps	n_miss	beta	se	l_remle	p_wald
1	rs3683945	3197400	0	-7.788676e-02	6.193497e-02	4.317971e+00	2.087605e-01
1	rs3707673	3407393	0	-6.654296e-02	6.210229e-02	4.316122e+00	2.841257e-01
1	rs6269442	3492195	0	-5.344258e-02	5.377462e-02	4.323589e+00	3.204786e-01
1	rs6336442	3580634	0	-6.770167e-02	6.209261e-02	4.315691e+00	2.757527e-01
1	rs13475700	4098402	0	-5.659112e-02	7.175371e-02	4.340122e+00	4.304285e-01

# If you do a manhattan plot in R, you will see a strong signal in chr17

# The log file also contains pve estimates and its standard error
pve estimate in the null model = 0.608891
se(pve) in the null model = 0.0327712











## To fit BSLMM in the training set:

## To fit a quantitative trait
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 2 -a mouse_hs1940.anno.txt -bslmm -o mouse_hs1940_CD8_bslmm -w 1000 -s 10000 -seed 1

# the following three files may be of most importance:
# the *.hyp.txt contains a column for pve and pge
# the *.param.txt contains estimates for betas, gammas and alphas
# the *.bv.txt contains breeding value estimates

## To fit a binary trait using a linear model
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 4 -a mouse_hs1940.anno.txt -bslmm -o mouse_hs1940_CD8_bslmm_cc1 -w 1000 -s 10000 -seed 1

## To fit a binary trait using a probit model instead
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 4 -a mouse_hs1940.anno.txt -bslmm 3 -o mouse_hs1940_CD8_bslmm_cc3 -w 1000 -s 10000 -seed 1

# The pve estimates in the log file are based on the standard linear model (i.e. on the observed scale), and so you will need to properly transform it back to the liability scale

## To generate relatedness matrix based on the training data
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 2 -a mouse_hs1940.anno.txt -gk 1 -o mouse_hs1940_CD8_train

# This matrix will only be required if you want to do prediction based on estimated breeding values
# Prediction can also be done without using the breeding values but instead using the alphas. This later approach does not appear to lose much accuracy in many examples we have encountered, although this may not be the case in your data.










## To obtain predicted values for the test set using estimates from BSLMM
## To do prediction in the test set for quantitative traits
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 2 -epm ./output/mouse_hs1940_CD8_bslmm.param.txt -emu ./output/mouse_hs1940_CD8_bslmm.log.txt -ebv ./output/mouse_hs1940_CD8_bslmm.bv.txt -k ./output/mouse_hs1940_CD8_train.cXX.txt -predict -o mouse_hs1940_CD8_prdt_k

## or use the alphas instead of breeding values
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 2 -epm ./output/mouse_hs1940_CD8_bslmm.param.txt -emu ./output/mouse_hs1940_CD8_bslmm.log.txt -predict -o mouse_hs1940_CD8_prdt

# The results will be inside ./output/*.prdt.txt
# If you load both results in R and check the mean squared error or correlation, you will find that both ways give very similar results. Both the correlation and the mean squared error should be around 0.65

## Now, do prediction in the test set for the binary traits
## If the traits were fitted using the linear model, then:
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 4 -epm ./output/mouse_hs1940_CD8_bslmm_cc1.param.txt -emu ./output/mouse_hs1940_CD8_bslmm_cc1.log.txt -predict -o mouse_hs1940_CD8_prdt_cc1

## If the traits were fitted using the probit model, then use predict option 2:
../bin/gemma -g mouse_hs1940.geno.txt -p mouse_hs1940_CD8.pheno.txt -n 4 -epm ./output/mouse_hs1940_CD8_bslmm_cc3.param.txt -emu ./output/mouse_hs1940_CD8_bslmm_cc3.log.txt -predict 2 -o mouse_hs1940_CD8_prdt_cc3

# You will find that fitting the binary traits using either the linear version or the probit version of BSLMM gives similar results. The brier scores should be around 0.19 and the area under the curve (AUC) should be around 0.78.