In [1]:

    

using DataFrames,JWAS,JWAS.Datasets


    

In [2]:

    

phenofile = Datasets.dataset("testMME","simple.txt")
genofile  = Datasets.dataset("testMME","genotype.txt")
pedfile   = Datasets.dataset("testMME","pedigree.txt");


    

In [3]:

    

;cat $phenofile


    

    




Animal,Age,y
S1,1,-0.92
D1,2,-1.05
O1,3,-0.92
O3,2,1.2

In [4]:

    

;cat $genofile


    

    




Animal,x1,x2,x3,x4,x5
S1,1,0,1,1,1
D1,2,0,2,2,1
O1,1,2,0,1,0
O3,0,0,2,1,1

In [5]:

    

;cat $pedfile


    

    




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

In [6]:

    

phenotype = readtable(phenofile,separator = ',',header=true);


    

In [7]:

    

pedigree = get_pedigree(pedfile);


    

    




Finished!

In [8]:

    

residual_variance = 1.0
genetic_variance  = 2.5
genetic_variance_by_marker    = 1.5
genetic_variance_by_polygenic = 1.0;


    

In [9]:

    

model = build_model("y = intercept + Age + Animal",residual_variance)
set_covariate(model,"Age")
set_random(model,"Animal",pedigree,genetic_variance_by_polygenic)


    

In [10]:

    

add_markers(model,genofile,genetic_variance_by_marker,separator=',');


    

    




5 markers on 4 individuals were added.

In [11]:

    

output=runMCMC(model,phenotype,chain_length=5000,
methods="BayesC",Pi=0.9,estimatePi=true,output_samples_frequency=100);


    

    




Priors for marker effects covariance matrix were calculated from genetic covariance matrix and π.
Marker effects covariance matrix is 
6.575342.


MCMC Information:
methods                                      BayesC
chain_length                                   5000
estimatePi                                     true
constraint                                    false
missing_phenotypes                            false
starting_value                                false
output_samples_frequency                        100
printout_frequency                             5001
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 BayesC...100%|█████████████████████████| Time: 0:00:01

In [12]:

    

keys(output)


    

    
Out[12]:





Base.KeyIterator for a Dict{Any,Any} with 6 entries. Keys:
  "Posterior mean of marker effects"
  "MCMC samples for residual variance"
  "MCMC samples for: π"
  "MCMC samples for polygenic effects var-cov parameters"
  "Posterior mean of location parameters"
  "Posterior mean of Pi"

In [13]:

    

output["Posterior mean of location parameters"]


    

    
Out[13]:





7x2 Array{Any,2}:
 "1:intercept : intercept"  -1.54912   
 "1:Age : Age"               0.569232  
 "1:Animal : S1"            -0.00425909
 "1:Animal : D1"            -0.0447008 
 "1:Animal : O1"            -0.0674396 
 "1:Animal : O3"             0.083868  
 "1:Animal : O2"            -0.0286708 

In [14]:

    

output["Posterior mean of marker effects"]


    

    
Out[14]:





5x2 Array{Any,2}:
 "x1"  -0.496822
 "x2"  -0.318663
 "x3"   0.175454
 "x4"  -0.225796
 "x5"   0.221099

In [16]:

    

using Plots
histogram(output["MCMC samples for: π"])


    

    
Out[16]:




    

    

In [ ]: