# new workflow compatible with ASREML-R 4 is described in the mtmm_workflow_as4.r script
# This script is optimized for data from the Arabidopsis 1001 Genomes project.

# EXAMPLE to run the data

source('scripts/prepare_mtmm.r')
load('data/MTMM_SAMPLE_DATA.Rdata')

# generate estimates of the variance components 
mtmm_estimates(Y,k=2,l=3,K,method='default',only.vca=FALSE) 
mydata<-'SD_SDV_mtmm_estimates.rda'
load(mydata)
# Now perform GWAS with this estimates 
results<-mtmm_part2(X,incl.singleGWAS=T)
plot_mtmm(name2='mtmm.pdf',incl.singleGWAS=T)

# all output data can also be found in the data folder

# old workflow with ASREML-R 3
# Load libraries and source needed functions
# The AsREML package needs a valid license that can be obtained at  http://www.vsni.co.uk/software/asreml

library(lattice)
library(asreml)

# msm and nadiv librarys are used to estimate SE of the correlation estimates, only used if run=FALSE 
#library(msm)
#library(nadiv)

source('mtmm_function.r')
source('emma.r')

# load your data (Phenotype(Y),Genotype(X) and Kinship(K))
# note you can calculate K using the emma package K<-emma.kinship(t(X)), make sure to set colnames(K)=rownames(K)=rownames(X)

# alternativley load the sample data
load('data/MTMM_SAMPLE_DATA.Rdata')

# different options include method(default or errorcorrelation, include.single.analysis, calculate.effect.size (if TRUE, #analysis is more time consuming) default for X is binary coding of 0 and 1, if your data are code  0,1 and 2 use #gen.data='heterozygot',  run=FALSE will not perform the GWAS, but only output the correlation estimates (fast)
mtmm(Y,X,K,method='default',include.single.analysis=T,calculate.effect.size=T,gen.data='binary',exclude=T,run=T)

# To only perform a Variance Coponent Analysis use the mtmm_estimate.r script with the flag only.vca=T set
VCA<-mtmm_estimates(Y,K=K,only.vca=T)

# the function outputs a list called results  ($phenotype ,$pvals, $statistics, $kinship)
output<-results$pvals

# manhattan plots
# default plots for include.single.analysis=T
par(mfrow=c(5,1),mar=c(3, 4, 1, 4))
plot_gwas(output,h=8)
plot_gwas(output,h=9)
plot_gwas(output,h=10)
plot_gwas(output,h=11)
plot_gwas(output,h=12)

#qq plots
par(mfrow=c(1,1),mar=c(3, 4, 1, 4))
qq_plot_all(output)