Multivariate Linear Additive Genetic Model
(with Maternal Effects)



In [1]:

    
using DataFrames,JWAS,JWAS.Datasets

I. Multivariate Linear Additive Genetic Model



In [2]:

    
phenofile = Datasets.dataset("testMT","phenotype.txt")
pedfile   = Datasets.dataset("testMT","pedigree.txt");

phenotypes



In [3]:

    
;cat $phenofile









    



Animal,BW,CW,age,sex
S1,100.0,10.0,8,M
D1,50.0,12.9,7,F
O1,150.0,13.0,3,M
O3,40.0,5.0,4,F

pedigree



In [4]:

    
;cat $pedfile









    



S1 0 0
D1 0 0
O1 S1 D1
O2 S1 D1
O3 S1 D1



In [5]:

    
data=readtable(phenofile)









    Out[5]:




Animal BW CW age sex
1 S1 100.0 10.0 8 M
2 D1 50.0 12.9 7 F
3 O1 150.0 13.0 3 M
4 O3 40.0 5.0 4 F

Genetic covariance matrix and residual covariance matrix



In [22]:

    
R      = [10.0 2.0
           2.0 1.0]
G      = [20.0 1.0
           1.0 2.0];



In [23]:

    
model_equations = "BW = intercept + age + sex + Animal;
                   CW = intercept + age + sex + Animal";



In [24]:

    
model    = build_model(model_equations,R);



In [25]:

    
set_covariate(model,"age");

set variables as random variables

random variables whose covariance matrices are numerator relationship matrix



In [26]:

    
ped = get_pedigree(pedfile);









    



Finished!



In [27]:

    
set_random(model,"Animal", ped,G)



In [28]:

    
out = runMCMC(model,data,chain_length=1000,printout_frequency=500,output_samples_frequency=10);









    



MCMC Information:
methods                        conventional (no markers)
chain_length                                   1000
estimatePi                                    false
constraint                                    false
missing_phenotypes                            false
starting_value                                false
output_samples_frequency                         10
printout_frequency                              500
update_priors_frequency                           0

Degree of freedom for hyper-parameters:
residual variances:                           4.000
iid random effect variances:                  4.000
polygenic effect variances:                   4.000
marker effect variances:                      4.000




Posterior means at iteration: 500
Residual covariance matrix: 
[39.861169 7.645041
 7.645041 2.173226]
Polygenic effects covariance matrix 
[101.435623 21.418891
 21.418891 6.633986]

running MCMC for conventional (no markers)... 53%|███   |  ETA: 0:00:00
Posterior means at iteration: 1000
Residual covariance matrix: 
[63.809784 14.233363
 14.233363 3.877617]
Polygenic effects covariance matrix 
[65.196935 12.086433
 12.086433 4.264667]

running MCMC for conventional (no markers)...100%|██████| Time: 0:00:00



In [29]:

    
keys(out)









    Out[29]:





Base.KeyIterator for a Dict{Any,Any} with 5 entries. Keys:
  "Posterior mean of polygenic effects covariance matrix"
  "Posterior mean of residual covariance matrix"
  "MCMC samples for residual covariance matrix"
  "MCMC samples for polygenic effects covariance matrix"
  "Posterior mean of location parameters"



In [30]:

    
out["MCMC samples for residual covariance matrix"]









    Out[30]:





4x100 Array{Float64,2}:
 362.562    27.2403   8.10654  13.218    …  29.9768   36.3408   48.5839 
  26.3573    5.41863  2.5419    3.22084      6.84176   8.60629  13.4308 
  26.3573    5.41863  2.5419    3.22084      6.84176   8.60629  13.4308 
   2.72721   1.18553  1.05415   1.09727      1.76212   2.62938   4.00713

II. The Multivariate Linear Additive Genetic Model with Maternal Effects



In [31]:

    
phenofile = Datasets.dataset("testMT","maternal.txt")
data=readtable(phenofile)









    Out[31]:




Animal BW CW age sex dam
1 S1 100.0 10.0 8 M 0
2 D1 50.0 12.9 7 F 0
3 O1 150.0 13.0 3 M D1
4 O3 40.0 5.0 4 F D1



In [32]:

    
model_equations = "BW = intercept + age + sex + Animal+ dam;
                   CW = intercept + age + sex + Animal";



In [33]:

    
model = build_model(model_equations,R);
set_covariate(model,"age");



In [34]:

    
# order is BW:Animal, BW:Dam, CW: Animal
G0 = [5   1    0.1
      1   1    0.01
      0.1 0.01 0.5] 
set_random(model,"Animal dam", ped,G0)



In [35]:

    
out=runMCMC(model,data,chain_length=1000,printout_frequency=500);









    



MCMC Information:
methods                        conventional (no markers)
chain_length                                   1000
estimatePi                                    false
constraint                                    false
missing_phenotypes                            false
starting_value                                false
output_samples_frequency                          0
printout_frequency                              500
update_priors_frequency                           0

Degree of freedom for hyper-parameters:
residual variances:                           4.000
iid random effect variances:                  4.000
polygenic effect variances:                   4.000
marker effect variances:                      4.000



running MCMC for conventional (no markers)... 39%|██    |  ETA: 0:00:00
Posterior means at iteration: 500
Residual covariance matrix: 
[117.436423 27.621048
 27.621048 7.033014]
Polygenic effects covariance matrix 
[4.674556 0.852744 0.097393
 0.852744 0.829121 -0.042028
 0.097393 -0.042028 0.440499]

running MCMC for conventional (no markers)... 81%|█████ |  ETA: 0:00:00
Posterior means at iteration: 1000
Residual covariance matrix: 
[127.724727 30.475367
 30.475367 7.767695]
Polygenic effects covariance matrix 
[4.788765 0.966935 0.175261
 0.966935 0.870485 0.024887
 0.175261 0.024887 0.495291]

running MCMC for conventional (no markers)...100%|██████| Time: 0:00:00



In [36]:

    
keys(out)









    Out[36]:





Base.KeyIterator for a Dict{Any,Any} with 3 entries. Keys:
  "Posterior mean of polygenic effects covariance matrix"
  "Posterior mean of residual covariance matrix"
  "Posterior mean of location parameters"