Merge pull request #23 from reworkhow/julia1.0

update JWAS to Julia 0.7

.travis.yml

@@ -4,10 +4,14 @@ os:
   - osx
 julia:
   - 0.6
+  - 0.7
+  - 1.0
   - nightly
 matrix:
   allow_failures:
     - julia: nightly
+    - julia: 1.0
+    - julia: 0.6
 notifications:
   email: false
 after_success:

REQUIRE

@@ -1,4 +1,4 @@
-julia 0.6
+julia 0.7
 Distributions
 DataFrames
 ProgressMeter

src/1.JWAS/src/JWAS.jl

@@ -1,5 +1,8 @@
-using Distributions
+using Distributions,Printf
+using DelimitedFiles
 using DataFrames,CSV
+using SparseArrays
+using LinearAlgebra
 using ProgressMeter
 using .PedModule
 using .misc

src/1.JWAS/src/MCMC/DRY.jl

@@ -24,19 +24,19 @@ function pre_check(mme,df,sol)
     #mme.M.obsID is IDs for all genotyped animals in complete genomic data
     phenoID = map(String,df[:,1])
     if mme.M!=0 && !issubset(phenoID,mme.M.obsID)
-      warn("Phenotyped individuals are not a subset of either\n",
+      printstyled("Phenotyped individuals are not a subset of either\n",
       "genotyped individuals (complete genomic data,non-single-step) or\n",
       "individuals in pedigree (incomplete genomic data, single-step).\n",
-      "Only individuals with both information are used in the analysis.\n")
+      "Only individuals with both information are used in the analysis.\n",bold=false,color=:red)
       index = [phenoID[i] in mme.M.obsID for i=1:length(phenoID)]
       df    = df[index,:]
     end
 
     if mme.M!=0 && mme.output_ID!=0 && !issubset(mme.output_ID,mme.M.obsID)
-      warn("Testing individuals are not a subset of \n",
+      printstyled("Testing individuals are not a subset of \n",
       "genotyped individuals (complete genomic data,non-single-step) or\n",
       "individuals in pedigree (incomplete genomic data, single-step).\n",
-      "Only tesing individuals with both information are used in the analysis.\n")
+      "Only tesing individuals with both information are used in the analysis.\n",bold=false,color=:red)
       mme.output_ID = intersect(mme.output_ID,mme.M.obsID)
     end
 
@@ -103,7 +103,7 @@ function output_result(mme,solMean,meanVare,G0Mean,output_samples_frequency,
         end
     end
 
-    output["Posterior mean of marker effects"] = (mme.nModels==1)?markerout[1]:markerout
+    output["Posterior mean of marker effects"] = (mme.nModels==1) ? markerout[1] : markerout
     output["Posterior mean of marker effects variance"] = meanVara
     if estimatePi == true
         output["Posterior mean of Pi"] = mean_pi
@@ -143,13 +143,13 @@ end
 ################################################################################
 function reformat2DataFrame(res::Array)
     ##SLOW, 130s
-    #out_names=[strip(i) for i in split(res[1,1],':',keep=false)]
+    #out_names=[strip(i) for i in split(res[1,1],':',keepempty=false)]
     #for rowi in 2:size(res,1)
-    #    out_names=[out_names [strip(i) for i in split(res[rowi,1],':',keep=false)]]#hcat two vectors
+    #    out_names=[out_names [strip(i) for i in split(res[rowi,1],':',keepempty=false)]]#hcat two vectors
     #end
-    out_names = Array{String}(size(res,1),3)
+    out_names = Array{String}(undef,size(res,1),3)
     for rowi in 1:size(res,1)
-        out_names[rowi,:]=[strip(i) for i in split(res[rowi,1],':',keep=false)]
+        out_names[rowi,:]=[strip(i) for i in split(res[rowi,1],':',keepempty=false)]
     end
 
     if size(out_names,2)==1 #convert vector to matrix
@@ -171,7 +171,9 @@ function genetic2marker(M::Genotypes,Pi::Dict)
   denom   = zeros(nTraits,nTraits)
   for i in 1:nTraits
     for j in i:nTraits
-      pi_selected = filter((k,v)->k[i]==1.0 && k[j]==1.0,Pi)
+      #pi_selected = filter((k,v)->k[i]==1.0 && k[j]==1.0,Pi)
+      pi_selected = filter(d->d.first[i]==1.0 && d.first[j]==1.0,Pi)
+
       denom[i,j] = M.sum2pq*sum(values(pi_selected))
       denom[j,i] = denom[i,j]
     end

src/1.JWAS/src/MCMC/MCMC.jl

@@ -90,13 +90,14 @@ function runMCMC(mme,df;
     if mme.M!=0
         #set up Pi
         if mme.nModels !=1 && Pi==0.0
-            info("Pi (Π) is not provided.","\n")
-            info("Pi (Π) is generated assuming all markers have effects on all traits.","\n")
-            mykey=Array{Float64}(0)
+            printstyled("Pi (Π) is not provided.\n",bold=false,color=:green)
+            printstyled("Pi (Π) is generated assuming all markers have effects on all traits.\n",bold=false,color=:green)
+            mykey=Array{Float64}(undef,0)
             ntraits=mme.nModels
             Pi=Dict{Array{Float64,1},Float64}()
-            for i in [ bin(n,ntraits) for n in 0:2^ntraits-1 ]
-              Pi[float(split(i,""))]=0.0
+            #for i in [ bin(n,ntraits) for n in 0:2^ntraits-1 ] `bin(n, pad)` is deprecated, use `string(n, base=2, pad=pad)
+            for i in [ string(n,base=2,pad=ntraits) for n in 0:2^ntraits-1 ]
+              Pi[parse.(Float64, split(i,""))]=0.0
             end
             Pi[ones(ntraits)]=1.0
         end
@@ -106,22 +107,22 @@ function runMCMC(mme,df;
             println()
             if mme.nModels != 1
               if !isposdef(mme.M.G) #also work for scalar
-                error("Marker effects covariance matrix is not postive definite! Please modify the argument: Pi.")
+                @error "Marker effects covariance matrix is not postive definite! Please modify the argument: Pi."
               end
-              println("The prior for marker effects covariance matrix is calculated from ")
-              println("genetic covariance matrix and Π. The prior for the marker effects ")
-              println("covariance matrix is: \n")
-              Base.print_matrix(STDOUT,round.(mme.M.G,6))
+              print("The prior for marker effects covariance matrix is calculated from ")
+              print("genetic covariance matrix and Π. The mean of the prior for the marker effects ")
+              print("covariance matrix is: \n")
+              Base.print_matrix(stdout,round.(mme.M.G,digits=6))
             else
               if !isposdef(mme.M.G) #positive scalar (>0)
-                error("Marker effects variance is negative!")
+                @error "Marker effects variance is negative!"
               end
-              println("The prior for marker effects variance is calculated from ")
-              println("the genetic variance and π. The prior for the marker effects variance ")
-              println("is: ",round.(mme.M.G,6))
+              print("The prior for marker effects variance is calculated from ")
+              print("the genetic variance and π. The mean of the prior for the marker effects variance ")
+              print("is: ",round.(mme.M.G,digits=6))
             end
         elseif methods=="GBLUP" && mme.M.G_is_marker_variance==true
-            error("Please provide genetic variance for GBLUP analysis")
+            @error "Please provide genetic variance for GBLUP analysis"
         end
         println("\n\n")
     end
@@ -160,7 +161,7 @@ function runMCMC(mme,df;
             error("No options!!!")
         end
     elseif mme.nModels > 1
-        if Pi != 0.0 && round(sum(values(Pi)),2)!=1.0
+        if Pi != 0.0 && round(sum(values(Pi)),digits=2)!=1.0
           error("Summation of probabilities of Pi is not equal to one.")
         end
         if methods in ["BayesC","BayesCC","BayesB","RR-BLUP","conventional (no markers)"]
@@ -198,14 +199,14 @@ function MCMCinfo(methods,Pi,chain_length,burnin,starting_value,printout_frequen
     @printf("%-30s %20s\n","chain_length",chain_length)
     @printf("%-30s %20s\n","burnin",burnin)
     if !(methods in ["conventional (no markers)", "GBLUP"])
-      @printf("%-30s %20s\n","estimatePi",estimatePi?"true":"false")
+      @printf("%-30s %20s\n","estimatePi",estimatePi ? "true" : "false")
     end
-    @printf("%-30s %20s\n","starting_value",starting_value?"true":"false")
+    @printf("%-30s %20s\n","starting_value",starting_value ? "true" : "false")
     @printf("%-30s %20d\n","printout_frequency",printout_frequency)
     @printf("%-30s %20d\n","output_samples_frequency",output_samples_frequency)
 
-    @printf("%-30s %20s\n","constraint",constraint?"true":"false")
-    @printf("%-30s %20s\n","missing_phenotypes",missing_phenotypes?"true":"false")
+    @printf("%-30s %20s\n","constraint",constraint ? "true" : "false")
+    @printf("%-30s %20s\n","missing_phenotypes",missing_phenotypes ? "true" : "false")
     @printf("%-30s %20d\n","update_priors_frequency",update_priors_frequency)
 
 
@@ -214,11 +215,11 @@ function MCMCinfo(methods,Pi,chain_length,burnin,starting_value,printout_frequen
     if mme.nModels==1
         for i in mme.rndTrmVec
             thisterm=split(i.term_array[1],':')[end]
-            @printf("%-30s %20s\n","random effect variances ("*thisterm*"):",round.(inv(i.GiNew),3))
+            @printf("%-30s %20s\n","random effect variances ("*thisterm*"):",round.(inv(i.GiNew),digits=3))
         end
         @printf("%-30s %20.3f\n","residual variances:",mme.RNew)
         if mme.pedTrmVec!=0
-            @printf("%-30s\n %50s\n","genetic variances (polygenic):",round.(inv(mme.GiNew),3))
+            @printf("%-30s\n %50s\n","genetic variances (polygenic):",round.(inv(mme.GiNew),digits=3))
         end
         if !(methods in ["conventional (no markers)", "GBLUP"])
             if mme.M == 0
@@ -232,23 +233,23 @@ function MCMCinfo(methods,Pi,chain_length,burnin,starting_value,printout_frequen
         for i in mme.rndTrmVec
             thisterm=split(i.term_array[1],':')[end]
             @printf("%-30s\n","random effect variances ("*thisterm*"):")
-            Base.print_matrix(STDOUT,round.(inv(i.GiNew),3))
+            Base.print_matrix(stdout,round.(inv(i.GiNew),digits=3))
             println()
         end
         @printf("%-30s\n","residual variances:")
-        Base.print_matrix(STDOUT,round.(mme.R,3))
+        Base.print_matrix(stdout,round.(mme.R,digits=3))
         println()
         if mme.pedTrmVec!=0
             @printf("%-30s\n","genetic variances (polygenic):")
-            Base.print_matrix(STDOUT,round.(inv(mme.Gi),3))
+            Base.print_matrix(stdout,round.(inv(mme.Gi),digits=3))
             println()
         end
         if !(methods in ["conventional (no markers)", "GBLUP"])
             @printf("%-30s\n","genetic variances (genomic):")
-            Base.print_matrix(STDOUT,round.(mme.M.G,3))
+            Base.print_matrix(stdout,round.(mme.M.G,digits=3))
             println()
             @printf("%-30s\n","marker effect variances:")
-            Base.print_matrix(STDOUT,round.(mme.M.G,3))
+            Base.print_matrix(stdout,round.(mme.M.G,digits=3))
             println()
             println("Π:")
             @printf("%-20s %12s\n","combinations","probability")

src/1.JWAS/src/MCMC/MCMC_BayesC.jl

@@ -13,7 +13,7 @@ function MCMC_BayesC(nIter,mme,df;
     # Pre-Check
     ############################################################################
     #starting values for location parameters(no marker) are sol
-    solMean     = zeros(sol)
+    solMean     = zero(sol)
 
     if methods == "conventional (no markers)"
         if mme.M!=0
@@ -83,7 +83,7 @@ function MCMC_BayesC(nIter,mme,df;
     #  WORKING VECTORS (ycor, saving values)
     ############################################################################
     #adjust y for starting values
-    ycorr       = vec(full(mme.ySparse)-mme.X*sol)
+    ycorr       = vec(Matrix(mme.ySparse)-mme.X*sol)
 
     ############################################################################
     #  SET UP OUTPUT MCMC samples
@@ -178,9 +178,9 @@ function MCMC_BayesC(nIter,mme,df;
         ########################################################################
         if output_samples_frequency != 0 && (iter-burnin)%output_samples_frequency==0 && iter>burnin
             if mme.M != 0
-                out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0?G0:false),π,α,vEff,outfile)
+                out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0 ? G0 : false),π,α,vEff,outfile)
             else
-                out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0?G0:false),false,false,false,outfile)
+                out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0 ? G0 : false),false,false,false,outfile)
             end
         end
 
@@ -189,16 +189,16 @@ function MCMC_BayesC(nIter,mme,df;
         ########################################################################
         if iter%outFreq==0 && iter>burnin
             println("\nPosterior means at iteration: ",iter)
-            println("Residual variance: ",round(meanVare,6))
+            println("Residual variance: ",round(meanVare,digits=6))
             if mme.pedTrmVec !=0
-                println("Polygenic effects covariance matrix \n",round.(G0Mean,3))
+                println("Polygenic effects covariance matrix \n",round.(G0Mean,digits=3))
             end
             if mme.M != 0
                 if methods!="BayesB"
-                    println("Marker effects variance: ",round(meanVara,6))
+                    println("Marker effects variance: ",round(meanVara,digits=6))
                 end
                 if estimatePi == true
-                    println("π: ", round(mean_pi,3))
+                    println("π: ", round(mean_pi,digits=3))
                 end
             end
         end
@@ -213,10 +213,10 @@ function MCMC_BayesC(nIter,mme,df;
       end
     end
     if mme.M != 0
-        output=output_result(mme,solMean,meanVare,(mme.pedTrmVec!=0?G0Mean:false),output_samples_frequency,
+        output=output_result(mme,solMean,meanVare,(mme.pedTrmVec!=0 ? G0Mean : false),output_samples_frequency,
                              meanAlpha,meanVara,estimatePi,mean_pi,output_file)
     else
-        output=output_result(mme,solMean,meanVare,(mme.pedTrmVec!=0?G0Mean:false),output_samples_frequency,
+        output=output_result(mme,solMean,meanVare,(mme.pedTrmVec!=0 ? G0Mean : false),output_samples_frequency,
                              false,false,false,false,output_file)
     end
     return output

src/1.JWAS/src/MCMC/MCMC_GBLUP.jl

@@ -144,19 +144,19 @@ function MCMC_GBLUP(nIter,mme,df;
         # 3.1 Save MCMC samples
         ########################################################################
         if output_samples_frequency != 0 && iter%output_samples_frequency==0 && iter>burnin
-            out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0?G0:false),0)
+            out_i=output_MCMC_samples(mme,out_i,sol,vRes,(mme.pedTrmVec!=0 ? G0 : false),0)
         end
         ########################################################################
         # 3.2 Printout
         ########################################################################
         if iter%outFreq==0
             println("\nPosterior means at iteration: ",iter)
-            println("Residual variance: ",round(meanVare,6))
+            println("Residual variance: ",round(meanVare,digits=6))
             if mme.pedTrmVec !=0
-              println("Polygenic effects covariance matrix \n",round(G0Mean,3))
+              println("Polygenic effects covariance matrix \n",round(G0Mean,digits=3))
             end
-            println("Genetic variance (G matrix): ",round(meanVara,6))
-            println("Genetic variance (GenSel): ",round(meanVarg,6))
+            println("Genetic variance (G matrix): ",round(meanVara,digits=6))
+            println("Genetic variance (GenSel): ",round(meanVarg,digits=6))
         end
     end