Univariate Linear Additive Genetic Model
(with Maternal Effects)

In [1]:
using DataFrames,JWAS,JWAS.Datasets
I. Univariate Linear Additive Genetic Model

In [2]:
pedfile   = Datasets.dataset("testMME","pedigree.txt")
phenofile = Datasets.dataset("testMME","simple.txt");

In [3]:
;cat $pedfile


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

In [4]:
d1 = readtable(phenofile)


Out[4]:
AnimalAgey
1S11-0.92
2D12-1.05
3O13-0.92
4O321.2

In [5]:
ped = get_pedigree(pedfile);


Finished!

In [6]:
varRes=1.0
model1 = build_model("y = intercept + Age + Animal",varRes)
set_covariate(model1,"Age")
G=2.5
set_random(model1,"Animal",ped,G)

In [7]:
out = solve(model1,d1,solver="GaussSeidel",printout_frequency=40)


40 3.3923200477199566e-5
Out[7]:
7x2 Array{Any,2}:
 "1:intercept : intercept"  -1.00435   
 "1:Age : Age"               0.204266  
 "1:Animal : S1"             0.116387  
 "1:Animal : D1"            -0.122438  
 "1:Animal : O1"            -0.294925  
 "1:Animal : O3"             0.996334  
 "1:Animal : O2"            -0.00302537
II. Univariate Linear Additive Genetic Model with Repeated Measures

In [8]:
phenofile = Datasets.dataset("testMME","repeated_measures.txt")
d2 = readtable(phenofile)


Out[8]:
AnimalAgey
1S11-0.92
2S12-1.35
3S13-0.33
4D11-0.3
5D12-1.05
6D130.56
7O11-0.09
8O120.44
9O13-0.92
10O31-0.2
11O321.2
12O330.25

In [9]:
varRes=1.0
model2 = build_model("y = intercept + Age + Animal + Animal*Age",varRes)
set_covariate(model2,"Age")
G = [1 0.1; 0.1 1.0]
set_random(model2,"Animal Animal*Age",ped,G)

In [10]:
out = solve(model2,d2,solver="Jacobi",printout_frequency=40)


40 0.004079837352814904
80 0.0008468686195409121
120 0.00026899969970286267
160 9.606419899377913e-5
200 3.649498312582733e-5
240 1.4825947812909268e-5
280 6.569648963398997e-6
320 3.2261059576412398e-6
360 1.7555444111655143e-6
400 1.0411088953585247e-6
Out[10]:
12x2 Array{Any,2}:
 "1:intercept : intercept"  -0.590464   
 "1:Age : Age"               0.0738261  
 "1:Animal : S1"            -0.0839778  
 "1:Animal : D1"             0.0822513  
 "1:Animal : O1"             0.191167   
 "1:Animal : O3"             0.204641   
 "1:Animal : O2"            -0.000908243
 "1:Animal*Age : S1"        -0.0943315  
 "1:Animal*Age : D1"         0.099499   
 "1:Animal*Age : O1"        -0.0351528  
 "1:Animal*Age : O3"         0.267083   
 "1:Animal*Age : O2"         0.00261206 

In [11]:
out = runMCMC(model2,d2,chain_length=1000);


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                             1001
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)...100%|██████| Time: 0:00:00
III. The Univariate Linear Additive Genetic Model with Maternal Effects

In [12]:
phenofile = Datasets.dataset("testMME","maternal_effects.txt")
d3 = readtable(phenofile)


Out[12]:
AnimalAgeyDam
1S11-0.920
2S12-1.350
3S13-0.330
4D11-0.30
5D12-1.050
6D130.560
7O11-0.09D1
8O120.44D1
9O13-0.92D1
10O31-0.2D1
11O321.2D1
12O330.25D1

In [13]:
;cat $pedfile


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

In [14]:
varRes = 1.0
model3 = build_model("y = intercept + Animal + Dam",varRes);

In [15]:
G = [1 0.1; 0.1 1.0]
set_random(model3,"Animal Dam",ped,G)

In [16]:
out = solve(model3,d3,solver="Gibbs",printout_frequency=1000)


at sample: 1000
at sample: 2000
at sample: 3000
at sample: 4000
at sample: 5000
Out[16]:
11x2 Array{Any,2}:
 "1:intercept : intercept"  -0.454248 
 "1:Animal : S1"            -0.282447 
 "1:Animal : D1"             0.184282 
 "1:Animal : O1"            -0.145197 
 "1:Animal : O3"             0.209786 
 "1:Animal : O2"            -0.0531358
 "1:Dam : S1"               -0.0683874
 "1:Dam : D1"                0.474598 
 "1:Dam : O1"                0.187702 
 "1:Dam : O3"                0.221185 
 "1:Dam : O2"                0.214125 

In [17]:
outputMCMCsamples(model3,"Animal");
out = runMCMC(model3,d3,chain_length=1000,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                             1001
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)...100%|██████| Time: 0:00:00

In [18]:
keys(out)


Out[18]:
Base.KeyIterator for a Dict{Any,Any} with 4 entries. Keys:
  "MCMC samples for residual variance"
  "MCMC samples for polygenic effects var-cov parameters"
  "MCMC samples for: 1:Animal"
  "Posterior mean of location parameters"

In [19]:
out["Posterior mean of location parameters"]


Out[19]:
11x2 Array{Any,2}:
 "1:intercept : intercept"  -0.422088  
 "1:Animal : S1"            -0.188057  
 "1:Animal : D1"             0.165045  
 "1:Animal : O1"            -0.0348238 
 "1:Animal : O3"             0.207546  
 "1:Animal : O2"            -0.0212569 
 "1:Dam : S1"                0.00715565
 "1:Dam : D1"                0.331039  
 "1:Dam : O1"                0.186682  
 "1:Dam : O3"                0.194227  
 "1:Dam : O2"                0.163045  

In [20]:
out["MCMC samples for: 1:Animal"]


Out[20]:
5x101 Array{Any,2}:
 "1:Animal : S1"  -0.411003   0.44366   …  -0.189172   -0.573982  -0.472356
 "1:Animal : D1"  -0.329615  -0.443609      0.396339   -0.626157   0.13115 
 "1:Animal : O1"  -0.552599  -0.313881      0.0682567  -0.562362  -0.257826
 "1:Animal : O3"  -0.363756   0.52053       0.534197   -0.662004   0.342766
 "1:Animal : O2"  -0.413678  -0.366833     -0.280669   -0.665699  -0.345913