In [1]:

    

srand(314)
using(Distributions)
d = Binomial(2,0.5)

nObs     = 10
nMarkers = 100
X        = float(rand(d,(nObs,nMarkers)))
α        = randn(nMarkers)
a        = X*α
stdGen   = std(a)
a        = a/stdGen
y        = a + randn(nObs);


    

In [2]:

    

include("../0.deprecated_BR/BR.jl")


    

    
Out[2]:





BR

In [3]:

    

using QTL


    

In [4]:

    

input=InputParameters()
input.method       = "BayesC"
input.varGenotypic = 1.0
input.varResidual  = 1.0
input.probFixed    = 0.95
input.chainLength  = 50000
input.centering    = true
input.outFreq      = 5000
input.estimatePi   = true

MCMCinfo(input)


    

    




MCMC Information:
seed                        314
chainLength               50000
method                   BayesC
outFreq                    5000
probFixed                 0.950
varGenotypic              1.000
varResidual               1.000
estimateVariance           true
estimatePi                 true
estimateScale             false
dfEffectVar               4.000
nuRes                     4.000
nuGen                     4.000
centering                  true

In [5]:

    

beta2=BR.runBR(input,genotype=X,phenotype=y)


    

    




running BayesC with a MCMC of length 50000
Iteration 5000 with 27 loci included in the model, mean residual/marker effect  0.595/ 0.223 with mean π=  0.840.
Iteration 10000 with 8 loci included in the model, mean residual/marker effect  0.599/ 0.218 with mean π=  0.839.
Iteration 15000 with 9 loci included in the model, mean residual/marker effect  0.599/ 0.217 with mean π=  0.836.
Iteration 20000 with 42 loci included in the model, mean residual/marker effect  0.600/ 0.217 with mean π=  0.840.
Iteration 25000 with 0 loci included in the model, mean residual/marker effect  0.596/ 0.221 with mean π=  0.842.
Iteration 30000 with 1 loci included in the model, mean residual/marker effect  0.594/ 0.221 with mean π=  0.841.
Iteration 35000 with 12 loci included in the model, mean residual/marker effect  0.596/ 0.223 with mean π=  0.845.
Iteration 40000 with 1 loci included in the model, mean residual/marker effect  0.596/ 0.222 with mean π=  0.845.
Iteration 45000 with 27 loci included in the model, mean residual/marker effect  0.594/ 0.223 with mean π=  0.847.
Iteration 50000 with 9 loci included in the model, mean residual/marker effect  0.593/ 0.222 with mean π=  0.844.






    
Out[5]:





100-element Array{Float64,1}:
  0.000281512
 -0.0204212  
 -0.0178922  
 -0.0504015  
  0.0199567  
 -0.0270632  
  0.00833018 
  0.00846933 
  0.00570045 
  0.0185589  
 -0.0163994  
 -0.0310323  
  0.00593951 
  ⋮          
  0.0143341  
  0.0274117  
  0.0152736  
 -0.015015   
  0.0121579  
  0.0145877  
  0.0112964  
 -0.00902406 
  0.00162967 
  0.00613139 
 -0.00774564 
  0.0116662  

In [6]:

    

srand(314)
using(Distributions)
d = Binomial(2,0.5)

nObs     = 10
nMarkers = 100
X        = float(rand(d,(nObs,nMarkers)))
α        = randn(nMarkers)
a        = X*α
stdGen   = std(a)
a        = a/stdGen
y        = a + randn(nObs);


    

In [7]:

    

include("../0.deprecated_oldJWAS/oldJWAS.jl")
using oldJWAS

myOption=Dict()
myOption["run"]           = "BayesC"
myOption["seed"]          = 123
myOption["chainLength"]   = 50000
myOption["probFixed"]     = 0.95
myOption["estimatePi"]    = "yes"
myOption["estimateScale"] = "no"
myOption["varGenotypic"]  = 1
myOption["varResidual"]   = 1

output = runJWAS(myOption,X,y)


    

    




This is iteration 5000, number of loci 8, vara 1.1334569514909414, vare 0.13972895384576264






    




WARNING: imported binding for run overwritten in module oldJWAS






    




This is iteration 10000, number of loci 7, vara 0.8671078760995132, vare 0.1594504285534714
This is iteration 15000, number of loci 29, vara 0.6803892381588219, vare 0.2642661248658907
This is iteration 20000, number of loci 2, vara 0.14332177567986606, vare 0.4196453538422861
This is iteration 25000, number of loci 20, vara 0.860933035908188, vare 0.27306921591839145
This is iteration 30000, number of loci 0, vara 0.0, vare 1.6503014560881413
This is iteration 35000, number of loci 1, vara 0.014639877551054697, vare 0.885866535290019
This is iteration 40000, number of loci 23, vara 0.5378351107174109, vare 0.32855495722363925
This is iteration 45000, number of loci 20, vara 0.4568528726042872, vare 0.2414369534805469
This is iteration 50000, number of loci 12, vara 0.6276096551455, vare 1.2435607348104465






    
Out[7]:





Dict{Any,Any} with 8 entries:
  "Proportion of variance… => 50000x1 Array{Float64,2}:…
  "posterior sample of ge… => [0.453909,0.926868,0.795386,0.288816,0.0682186,0.…
  "posterior mean of fixe… => [1.6275769933768554]
  "posterior mean of mark… => [0.000171776,-0.0224289,-0.0170548,-0.0527296,0.0…
  "posterior sample of pi" => [0.97659,0.843546,0.886302,0.947962,0.941708,0.85…
  "posterior sample of re… => [0.494222,0.375909,0.506723,0.767744,0.555247,0.8…
  "posterior sample of sc… => [0.21097,0.21097,0.21097,0.21097,0.21097,0.21097,…
  "model frequency"        => [0.14994,0.16468,0.1671,0.20104,0.16138,0.17086,0…

In [8]:

    

beta1=output["posterior mean of marker effects"];
beta1


    

    
Out[8]:





100-element Array{Float64,1}:
  0.000171776
 -0.0224289  
 -0.0170548  
 -0.0527296  
  0.0211365  
 -0.0278139  
  0.00817019 
  0.00985067 
  0.00618653 
  0.020807   
 -0.0207543  
 -0.0336786  
  0.00608711 
  ⋮          
  0.0148481  
  0.0279333  
  0.0162836  
 -0.0156757  
  0.0145379  
  0.0144553  
  0.013036   
 -0.00927405 
  0.00410696 
  0.00379601 
 -0.00802192 
  0.0144613  

In [9]:

    

[beta1 beta2]


    

    
Out[9]:





100x2 Array{Float64,2}:
  0.000171776   0.000281512
 -0.0224289    -0.0204212  
 -0.0170548    -0.0178922  
 -0.0527296    -0.0504015  
  0.0211365     0.0199567  
 -0.0278139    -0.0270632  
  0.00817019    0.00833018 
  0.00985067    0.00846933 
  0.00618653    0.00570045 
  0.020807      0.0185589  
 -0.0207543    -0.0163994  
 -0.0336786    -0.0310323  
  0.00608711    0.00593951 
  ⋮                        
  0.0148481     0.0143341  
  0.0279333     0.0274117  
  0.0162836     0.0152736  
 -0.0156757    -0.015015   
  0.0145379     0.0121579  
  0.0144553     0.0145877  
  0.013036      0.0112964  
 -0.00927405   -0.00902406 
  0.00410696    0.00162967 
  0.00379601    0.00613139 
 -0.00802192   -0.00774564 
  0.0144613     0.0116662  

In [10]:

    

cor(X*beta1,X*beta2)


    

    
Out[10]:





0.9999913011034178

In [11]:

    

using DataFrames


    

In [12]:

    

#writedlm("geno.txt",X)


    

In [13]:

    

test=DataFrame(y=y)


    

    
Out[13]:




y
1
1.8031502046877939
2
3.376068808270318
3
0.05449720591496543
4
1.9547479445616438
5
0.510558668932955
6
1.552954739410723
7
1.618770528203041
8
2.3041115374993426
9
1.1268729585780282
10
1.9772296992742773

In [14]:

    

include("../../../MMEModule.jl/src/MMEModule.jl")


    

    
Out[14]:





MMEModule

In [15]:

    

using MMEModule


    

In [16]:

    

varRes=1.0
mme = MMEModule.build_model("y = intercept",varRes);


    

In [17]:

    

MMEModule.addMarkers(mme,"geno.txt",1.0,header=false,separator='\t');


    

    




The delimiters in file geno.txt is 	  .

In [18]:

    

output=MMEModule.runMCMC(mme,test,chain_length=100000,
                         printout_frequency=5000,methods="BayesC",Pi=0.95,
estimatePi=true)


    

    




at sample: 5000 with meanVare: 0.5780393961589684
at sample: 10000 with meanVare: 0.5836559699036918
at sample: 15000 with meanVare: 0.5845246683908016
at sample: 20000 with meanVare: 0.5872376800261305
at sample: 25000 with meanVare: 0.5861955212148773
at sample: 30000 with meanVare: 0.5890290536231467
at sample: 35000 with meanVare: 0.5883784057708568
at sample: 40000 with meanVare: 0.5890896188674953
at sample: 45000 with meanVare: 0.5885332488511633
at sample: 50000 with meanVare: 0.5872151165210724
at sample: 55000 with meanVare: 0.5864796006132
at sample: 60000 with meanVare: 0.5877527153281169
at sample: 65000 with meanVare: 0.5885307347575364
at sample: 70000 with meanVare: 0.5885050055684955
at sample: 75000 with meanVare: 0.5888472955190147
at sample: 80000 with meanVare: 0.5879852432499618
at sample: 85000 with meanVare: 0.587851340793821
at sample: 90000 with meanVare: 0.5881579539539179
at sample: 95000 with meanVare: 0.5889035567045622
at sample: 100000 with meanVare: 0.589698383903424






    
Out[18]:





Dict{Any,Any} with 4 entries:
  "Posterior Mean of Loca… => 1x2 Array{Any,2}:…
  "MCMC samples for resid… => [0.566783,0.348597,0.362075,0.330791,0.388863,0.1…
  "Posterior Mean of Mark… => [2.63001e-6,-0.0231469,-0.0178541,-0.0506663,0.01…
  "MCMC samples for: π"    => [0.902196,0.92549,0.885801,0.93874,0.924486,0.942…

In [19]:

    

keys(output)


    

    
Out[19]:





Base.KeyIterator for a Dict{Any,Any} with 4 entries. Keys:
  "Posterior Mean of Location Parameters"
  "MCMC samples for residual variance"
  "Posterior Mean of Marker Effects"
  "MCMC samples for: π"

In [20]:

    

beta3=output["Posterior Mean of Location Parameters"]


    

    
Out[20]:





1x2 Array{Any,2}:
 "intercept: intercept"  1.62761

In [21]:

    

beta3=output["Posterior Mean of Marker Effects"]


    

    
Out[21]:





100-element Array{Float64,1}:
  2.63001e-6
 -0.0231469 
 -0.0178541 
 -0.0506663 
  0.0196927 
 -0.0272377 
  0.00800267
  0.00818424
  0.00591928
  0.0204951 
 -0.0181813 
 -0.0335485 
  0.00627734
  ⋮         
  0.0149391 
  0.0276611 
  0.0162816 
 -0.0157178 
  0.0134405 
  0.0147612 
  0.0129465 
 -0.00772476
  0.00282223
  0.00439893
 -0.00824716
  0.0133199 

In [22]:

    

[beta1 beta2 beta3]


    

    
Out[22]:





100x3 Array{Float64,2}:
  0.000171776   0.000281512   2.63001e-6
 -0.0224289    -0.0204212    -0.0231469 
 -0.0170548    -0.0178922    -0.0178541 
 -0.0527296    -0.0504015    -0.0506663 
  0.0211365     0.0199567     0.0196927 
 -0.0278139    -0.0270632    -0.0272377 
  0.00817019    0.00833018    0.00800267
  0.00985067    0.00846933    0.00818424
  0.00618653    0.00570045    0.00591928
  0.020807      0.0185589     0.0204951 
 -0.0207543    -0.0163994    -0.0181813 
 -0.0336786    -0.0310323    -0.0335485 
  0.00608711    0.00593951    0.00627734
  ⋮                                     
  0.0148481     0.0143341     0.0149391 
  0.0279333     0.0274117     0.0276611 
  0.0162836     0.0152736     0.0162816 
 -0.0156757    -0.015015     -0.0157178 
  0.0145379     0.0121579     0.0134405 
  0.0144553     0.0145877     0.0147612 
  0.013036      0.0112964     0.0129465 
 -0.00927405   -0.00902406   -0.00772476
  0.00410696    0.00162967    0.00282223
  0.00379601    0.00613139    0.00439893
 -0.00802192   -0.00774564   -0.00824716
  0.0144613     0.0116662     0.0133199 

In [26]:

    

cor(X*beta2,X*beta1)


    

    
Out[26]:





0.9999913011034178

In [27]:

    

cor(X*beta2,X*beta3)


    

    
Out[27]:





0.9999922682855739

In [28]:

    

using oldJWAS

myOption=Dict()
myOption["run"]           = "BayesB"
myOption["seed"]          = 123
myOption["chainLength"]   = 50000
myOption["probFixed"]     = 0.95
myOption["estimatePi"]    = "no"
myOption["estimateScale"] = "no"
myOption["varGenotypic"]  = 1
myOption["varResidual"]   = 1

output = runJWAS(myOption,X,y)


    

    




This is iteration 5000, number of loci 3and vare 0.2903985218537527
This is iteration 10000, number of loci 6and vare 0.7742689898990693
This is iteration 15000, number of loci 5and vare 0.9599325469252361
This is iteration 20000, number of loci 4and vare 0.4232268341835519
This is iteration 25000, number of loci 5and vare 0.2588546300753118
This is iteration 30000, number of loci 7and vare 0.7562739905523324
This is iteration 35000, number of loci 3and vare 1.426257844833594
This is iteration 40000, number of loci 5and vare 0.3316108868243706
This is iteration 45000, number of loci 3and vare 0.9173205606441555
This is iteration 50000, number of loci 7and vare 0.46927460533298676






    
Out[28]:





Dict{Any,Any} with 8 entries:
  "Proportion of variance… => 10x1 Array{Float64,2}:…
  "posterior sample of ge… => [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
  "posterior mean of fixe… => [1.6274526471727675]
  "posterior mean of mark… => 100x1 Array{Float64,2}:…
  "posterior sample of pi" => [0.95,0.95,0.95,0.95,0.95,0.95,0.95,0.95,0.95,0.9…
  "posterior sample of re… => [0.912059,0.775106,1.17091,2.72984,1.79787,1.3286…
  "posterior sample of sc… => [0.21097,0.21097,0.21097,0.21097,0.21097,0.21097,…
  "model frequency"        => [0.03454,0.05202,0.04996,0.09394,0.04592,0.05684,…

In [29]:

    

output2=MMEModule.runMCMC(mme,test,chain_length=100000,
                         printout_frequency=5000,methods="BayesB",Pi=0.95)


    

    




at sample: 5000 with meanVare: 0.4850655956242786
at sample: 10000 with meanVare: 0.4857213455852994
at sample: 15000 with meanVare: 0.4844783222476542
at sample: 20000 with meanVare: 0.4866607089837925
at sample: 25000 with meanVare: 0.4851175084213007
at sample: 30000 with meanVare: 0.48643197635976904
at sample: 35000 with meanVare: 0.48801046708100176
at sample: 40000 with meanVare: 0.48804877581836204
at sample: 45000 with meanVare: 0.48916323177240534
at sample: 50000 with meanVare: 0.4895638076660683
at sample: 55000 with meanVare: 0.4885267165067474
at sample: 60000 with meanVare: 0.4886563804981134
at sample: 65000 with meanVare: 0.48706036687552823
at sample: 70000 with meanVare: 0.48749920171226885
at sample: 75000 with meanVare: 0.48813434839479536
at sample: 80000 with meanVare: 0.4875085203173253
at sample: 85000 with meanVare: 0.48747722202855476
at sample: 90000 with meanVare: 0.4873177672641523
at sample: 95000 with meanVare: 0.48688043240964146
at sample: 100000 with meanVare: 0.48657577582269623






    
Out[29]:





Dict{Any,Any} with 3 entries:
  "Posterior Mean of Loca… => 1x2 Array{Any,2}:…
  "MCMC samples for resid… => [0.704782,0.479053,0.234446,0.269246,0.166422,0.4…
  "Posterior Mean of Mark… => [-0.000503365,-0.0146095,-0.016149,-0.0541999,0.0…

In [30]:

    

beta1=output["posterior mean of marker effects"];
beta2=output2["Posterior Mean of Marker Effects"];


    

In [31]:

    

cor(beta1,beta2)


    

    
Out[31]:





1x1 Array{Float64,2}:
 0.996402

In [32]:

    

cor(X*beta1,X*beta2)


    

    
Out[32]:





1x1 Array{Float64,2}:
 0.999861