From e7b09bf6e3b483814730d88ab50ecf9184b9dd68 Mon Sep 17 00:00:00 2001 From: haocheng Date: Fri, 17 Aug 2018 21:57:35 -0700 Subject: [PATCH 1/8] fix warnings in julia0.7 --- src/1.JWAS/src/JWAS.jl | 3 +- src/1.JWAS/src/MCMC/MCMC.jl | 8 +-- src/1.JWAS/src/build_MME.jl | 4 +- src/1.JWAS/src/random_effects.jl | 6 +-- src/1.JWAS/src/types.jl | 14 ++--- src/2.PedModule/src/PedModule.jl | 14 ++--- src/3.SSBR/src/SSBR_types.jl | 20 +++---- src/4.misc/src/GWAS.jl | 2 +- src/4.misc/src/QTL_types.jl | 9 ++-- src/4.misc/src/files.jl | 91 +++++++++++++++++--------------- src/4.misc/src/fixed_effects.jl | 2 +- src/4.misc/src/genotypes.jl | 2 +- src/4.misc/src/qc/qc.jl | 1 - 13 files changed, 90 insertions(+), 86 deletions(-) diff --git a/src/1.JWAS/src/JWAS.jl b/src/1.JWAS/src/JWAS.jl index c85d8bb8..32761f7c 100644 --- a/src/1.JWAS/src/JWAS.jl +++ b/src/1.JWAS/src/JWAS.jl @@ -1,5 +1,6 @@ -using Distributions +using Distributions,Printf using DataFrames,CSV +using SparseArrays using ProgressMeter using .PedModule using .misc diff --git a/src/1.JWAS/src/MCMC/MCMC.jl b/src/1.JWAS/src/MCMC/MCMC.jl index f3a98420..dcdc4f8c 100644 --- a/src/1.JWAS/src/MCMC/MCMC.jl +++ b/src/1.JWAS/src/MCMC/MCMC.jl @@ -198,14 +198,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) diff --git a/src/1.JWAS/src/build_MME.jl b/src/1.JWAS/src/build_MME.jl index 7cc27dad..2595ac7c 100644 --- a/src/1.JWAS/src/build_MME.jl +++ b/src/1.JWAS/src/build_MME.jl @@ -309,8 +309,8 @@ function getinfo(model;data=false) end nLevels = i.nLevels - fixed = (term in random_effects)?"random":"fixed" - factor = (nLevels==1)?"covariate":"factor" + fixed = (term in random_effects) ? "random" : "fixed" + factor = (nLevels==1) ? "covariate" : "factor" if size(model.mmeRhs)==()&&data==false nLevels="NA" diff --git a/src/1.JWAS/src/random_effects.jl b/src/1.JWAS/src/random_effects.jl index 950541bc..de7e1a43 100644 --- a/src/1.JWAS/src/random_effects.jl +++ b/src/1.JWAS/src/random_effects.jl @@ -185,7 +185,7 @@ function addA(mme::MME) #two terms,e.g.,"animal" and "maternal" may not near in startPosj= pedTrmj.startPos endPosj = startPosj + pedTrmj.nLevels - 1 #lamda version (single trait) or not (multi-trait) - myaddA = (mme.nModels!=1)?(mme.Ai*mme.Gi[i,j]):(mme.Ai*(mme.GiNew[i,j]*mme.RNew - mme.GiOld[i,j]*mme.ROld)) + myaddA = (mme.nModels!=1) ? (mme.Ai*mme.Gi[i,j]) : (mme.Ai*(mme.GiNew[i,j]*mme.RNew - mme.GiOld[i,j]*mme.ROld)) mme.mmeLhs[startPosi:endPosi,startPosj:endPosj] = mme.mmeLhs[startPosi:endPosi,startPosj:endPosj] + myaddA end @@ -200,7 +200,7 @@ end function addLambdas(mme::MME) for random_term in mme.rndTrmVec term_array = random_term.term_array - Vi = (random_term.Vinv!=0)?random_term.Vinv:speye(mme.modelTermDict[term_array[1]].nLevels) + Vi = (random_term.Vinv!=0) ? random_term.Vinv : speye(mme.modelTermDict[term_array[1]].nLevels) for (i,termi) = enumerate(term_array) randTrmi = mme.modelTermDict[termi] startPosi = randTrmi.startPos @@ -209,7 +209,7 @@ function addLambdas(mme::MME) randTrmj = mme.modelTermDict[termj] startPosj = randTrmj.startPos endPosj = startPosj + randTrmj.nLevels - 1 - myaddLambda = (mme.nModels!=1)?(Vi*random_term.Gi[i,j]):(Vi*(random_term.GiNew[i,j]*mme.RNew - random_term.GiOld[i,j]*mme.ROld)) + myaddLambda = (mme.nModels!=1) ? (Vi*random_term.Gi[i,j]) : (Vi*(random_term.GiNew[i,j]*mme.RNew - random_term.GiOld[i,j]*mme.ROld)) mme.mmeLhs[startPosi:endPosi,startPosj:endPosj] = mme.mmeLhs[startPosi:endPosi,startPosj:endPosj] + myaddLambda end diff --git a/src/1.JWAS/src/types.jl b/src/1.JWAS/src/types.jl index dd81969f..52943fee 100644 --- a/src/1.JWAS/src/types.jl +++ b/src/1.JWAS/src/types.jl @@ -7,7 +7,7 @@ # terms: 1:A,1:B,2:A,2:B,2:A*B # ################################################################################ -type ModelTerm +mutable struct ModelTerm iModel::Int64 # 1st or 2nd model_equation # | trmStr | nFactors | factors | @@ -52,7 +52,7 @@ end #In JWAS, ONLY used when residual variance is constant #or missing phenotypes are not imputed at each step of MCMC (no marker effects). ################################################################################ -type ResVar +mutable struct ResVar R0::Array{Float64,2} RiDict::Dict{BitArray{1},Array{Float64,2}} end @@ -63,7 +63,7 @@ end #single-trait e.g. termarray: [ModelTerm(1:A)] #multi-trait e.g. termarray: [ModelTerm(1:A), ModelTerm(2:A)] ################################################################################ -type RandomEffect #Better to be a dict? key: term_array::Array{AbstractString,1}?? +mutable struct RandomEffect #Better to be a dict? key: term_array::Array{AbstractString,1}?? term_array::Array{AbstractString,1} Gi::Array{Float64,2} #covariance matrix (multi-trait) GiOld::Array{Float64,2} #specific for lambda version of MME (single-trait) @@ -75,12 +75,12 @@ type RandomEffect #Better to be a dict? key: term_array::Array{AbstractString, end #MCMC samplers for location parameters -type MCMCSamples +mutable struct MCMCSamples term::ModelTerm sampleArray::Array{Float64,2} end -type Genotypes +mutable struct Genotypes obsID::Array{AbstractString,1} #row ID of genotypes markerID nObs::Int64 @@ -94,7 +94,7 @@ type Genotypes Genotypes(a1,a2,a3,a4,a5,a6,a7,a8)=new(a1,a2,a3,a4,a5,a6,a7,a8,0.0,false) end -type DF +mutable struct DF residual::Float64 polygenic::Float64 marker::Float64 @@ -107,7 +107,7 @@ end # single-trait analyses: lambda version of MME # multi-trait analysis: formal version of MME ################################################################################ -type MME +mutable struct MME nModels::Int64 #number of model equations modelVec::Array{AbstractString,1} #["y1 = A + B","y2 = A + B + A*B"] modelTerms::Array{ModelTerm,1} #ModelTerms for "1:intercept","1:A","2:intercept","2:A","2:A*B"...; diff --git a/src/2.PedModule/src/PedModule.jl b/src/2.PedModule/src/PedModule.jl index 013bfa80..041fbd99 100644 --- a/src/2.PedModule/src/PedModule.jl +++ b/src/2.PedModule/src/PedModule.jl @@ -3,14 +3,14 @@ module PedModule using DataFrames,CSV using ProgressMeter -type PedNode +mutable struct PedNode seqID::Int64 sire::String dam::String f::Float64 end -type Pedigree +mutable struct Pedigree currentID::Int64 idMap::Dict{AbstractString,PedNode} aij::Dict{Int64, Float64} @@ -61,7 +61,7 @@ function calcAddRel!(ped::Pedigree,id1::AbstractString,id2::AbstractString) if id1=="0" || id2=="0" # zero return 0.0 end - old,yng = ped.idMap[id1].seqID0 && damPos>0 d = sqrt(4.0/(2 - ped.idMap[sire].f - ped.idMap[dam].f)) diff --git a/src/3.SSBR/src/SSBR_types.jl b/src/3.SSBR/src/SSBR_types.jl index 56067427..74442fa1 100644 --- a/src/3.SSBR/src/SSBR_types.jl +++ b/src/3.SSBR/src/SSBR_types.jl @@ -1,4 +1,4 @@ -type Numbers +mutable struct Numbers ped::Int64 #individulals in pedigree pedn::Int64 #non-genotyped individuals in pedigree pedg::Int64 #genotyped individuals in pedigree @@ -8,7 +8,7 @@ type Numbers markers::Int64 #number of markers end -type ZMats +mutable struct ZMats full::SparseMatrixCSC{Float64,Int64} #Z # Z= | Zn 0 |=|Z_n Z_g| n::SparseMatrixCSC{Float64,Int64} #Zn # | 0 Zg| g::SparseMatrixCSC{Float64,Int64} #Zg # @@ -16,44 +16,44 @@ type ZMats _g::SparseMatrixCSC{Float64,Int64} #Z_g # |0 | end -type AiMats +mutable struct AiMats full::SparseMatrixCSC{Float64,Int64} #Ai nn::SparseMatrixCSC{Float64,Int64} #Ai_nn ng::SparseMatrixCSC{Float64,Int64} #Ai_ng end -type YVecs +mutable struct YVecs full #::Array{Float64,1} #y n #::Array{Float64,1} #yn g #::Array{Float64,1} #yg #order of ids is same to order of y ids #::Array{String,1} #order of ids is nongeno then geno end -type MMats +mutable struct MMats full::Array{Float64,2} n::Array{Float64,2} #Mn g::Array{Float64,2} #Mg end -type JVecs +mutable struct JVecs full::Array{Float64,2} n::Array{Float64,2} #Jn g::Array{Float64,2} #Jg end -type XMats #sparse may be better +mutable struct XMats #sparse may be better full::Array{Float64,2} n::Array{Float64,2} #Xn g::Array{Float64,2} #Xg end -type WMats +mutable struct WMats full::Array{Float64,2} n::Array{Float64,2} #Wn g::Array{Float64,2} #Wg end -type HybridMatrices +mutable struct HybridMatrices Z::ZMats Ai::AiMats y::YVecs @@ -64,7 +64,7 @@ type HybridMatrices num::Numbers end -type Results +mutable struct Results mats::HybridMatrices posMeans::misc.Output ebv diff --git a/src/4.misc/src/GWAS.jl b/src/4.misc/src/GWAS.jl index 54273b27..d770eb7e 100644 --- a/src/4.misc/src/GWAS.jl +++ b/src/4.misc/src/GWAS.jl @@ -53,7 +53,7 @@ function GWAS(marker_file,mme;header=true,window_size=100,threshold=0.001) for win=1:nWindows wStart = wEnd + 1 wEnd += window_size[windowi] - wEnd = (wEnd > nMarkers) ? nMarkers:wEnd + wEnd = (wEnd > nMarkers) ? nMarkers : wEnd winVarProps[i,win] = var(X[:,wStart:wEnd]*α[wStart:wEnd])/genVar windowi +=1 diff --git a/src/4.misc/src/QTL_types.jl b/src/4.misc/src/QTL_types.jl index b7ed4282..aa42624d 100644 --- a/src/4.misc/src/QTL_types.jl +++ b/src/4.misc/src/QTL_types.jl @@ -1,4 +1,5 @@ -type InputParameters +using Printf +mutable struct InputParameters seed::Int64 # seed method::AbstractString # BaysABC chainLength::Int64 # number of iterations @@ -35,7 +36,7 @@ function MCMCinfo(input::InputParameters) @printf("%-20s %10s\n","centering",input.centering) end -type GibbsMats #in src/1.JWAS/markers/tools +mutable struct GibbsMats #in src/1.JWAS/markers/tools X::Array{Float64,2} nrows::Int64 ncols::Int64 @@ -49,7 +50,7 @@ type GibbsMats #in src/1.JWAS/markers/tools end end -type Current +mutable struct Current varGenotypic::Float64 #genotypic variance varResidual::Float64 #residual variance varEffect::Float64 # common marker variance @@ -102,7 +103,7 @@ type Current end end -type Output +mutable struct Output mean_fixed_effects::Array{Float64,1} meanMarkerEffects::Array{Float64,1} modelFreq::Array{Float64,1} diff --git a/src/4.misc/src/files.jl b/src/4.misc/src/files.jl index 7b77f598..d815d64d 100644 --- a/src/4.misc/src/files.jl +++ b/src/4.misc/src/files.jl @@ -19,47 +19,50 @@ end # colLabels = "", # rowLabels = "" -function latexTable{T}(A::Array{T,2}; - fileName::AbstractString = "", - colLabels = "", - rowLabels = "" - ) - if fileName == "" - outFile = STDOUT - else - outFile = open(fileName, "w") - end - rows, cols = size(A) - println(outFile,"\\begin{center}") - print(outFile,"\\begin{tabular}{",) - if rowLabels!="" - print(outFile,"l") - end - for j=1:cols - print(outFile,"r") - end - println(outFile,"}\\hline") - if colLabels!="" - nCol = length(colLabels) - for j = 1:(nCol-1) - print(outFile,colLabels[j]," & ") - end - print(outFile,colLabels[nCol]," \\\\ \\hline \n") - end - for i =1:rows - if rowLabels!="" - print(outFile,rowLabels[i]," & ") - end - for j = 1:(cols-1) - print(outFile,A[i,j]) - print(outFile," & ") - end - print(outFile,A[i,cols]," \\\\ \n") - end - println(outFile,"\\hline") - println(outFile,"\\end{tabular}") - println(outFile,"\\end{center}") - if fileName != "" - close(outFile) - end -end +# ┌ Warning: Deprecated syntax `parametric method syntax latexTable{T}(A::Array{T, 2},; fileName::AbstractString = "", colLabels = "", rowLabels = "")` around /home/jovyan/qtlcheng/Github/JWAS.jl/src/4.misc/src/files.jl:27. +# │ Use `latexTable(A::Array{T, 2},; fileName::AbstractString = "", colLabels = "", rowLabels = "") where T` instead. +# └ @ nothing /home/jovyan/qtlcheng/Github/JWAS.jl/src/4.misc/src/files.jl:27 +# function latexTable{T}(A::Array{T,2}; +# fileName::AbstractString = "", +# colLabels = "", +# rowLabels = "" +# ) +# if fileName == "" +# outFile = STDOUT +# else +# outFile = open(fileName, "w") +# end +# rows, cols = size(A) +# println(outFile,"\\begin{center}") +# print(outFile,"\\begin{tabular}{",) +# if rowLabels!="" +# print(outFile,"l") +# end +# for j=1:cols +# print(outFile,"r") +# end +# println(outFile,"}\\hline") +# if colLabels!="" +# nCol = length(colLabels) +# for j = 1:(nCol-1) +# print(outFile,colLabels[j]," & ") +# end +# print(outFile,colLabels[nCol]," \\\\ \\hline \n") +# end +# for i =1:rows +# if rowLabels!="" +# print(outFile,rowLabels[i]," & ") +# end +# for j = 1:(cols-1) +# print(outFile,A[i,j]) +# print(outFile," & ") +# end +# print(outFile,A[i,cols]," \\\\ \n") +# end +# println(outFile,"\\hline") +# println(outFile,"\\end{tabular}") +# println(outFile,"\\end{center}") +# if fileName != "" +# close(outFile) +# end +# end diff --git a/src/4.misc/src/fixed_effects.jl b/src/4.misc/src/fixed_effects.jl index d8c4375b..2e6eacf8 100644 --- a/src/4.misc/src/fixed_effects.jl +++ b/src/4.misc/src/fixed_effects.jl @@ -1,4 +1,4 @@ -type FixedMatrix +mutable struct FixedMatrix C::Array{Float64,2} variables::Array{Any,1} end diff --git a/src/4.misc/src/genotypes.jl b/src/4.misc/src/genotypes.jl index b2712cc1..2bb9bc1b 100644 --- a/src/4.misc/src/genotypes.jl +++ b/src/4.misc/src/genotypes.jl @@ -1,4 +1,4 @@ -type Genotypes +mutable struct Genotypes obsID::Array{AbstractString,1} #row ID of genotypes markerID nObs::Int64 diff --git a/src/4.misc/src/qc/qc.jl b/src/4.misc/src/qc/qc.jl index 5b6139d4..cac16f68 100644 --- a/src/4.misc/src/qc/qc.jl +++ b/src/4.misc/src/qc/qc.jl @@ -39,7 +39,6 @@ O1,1,2,0,1,0 O3,0,0,2,1,1 ``` """ - function QC(infile,outfile;separator=' ',header=true,missing=false,MAF=0.1) myfile = open(infile) From 521e677dbc1189684a1a3678d1c379f9b31c840d Mon Sep 17 00:00:00 2001 From: haocheng Date: Fri, 17 Aug 2018 23:27:02 -0700 Subject: [PATCH 2/8] fix warnings in julia0.7 --- src/1.JWAS/src/MCMC/DRY.jl | 2 +- src/1.JWAS/src/MCMC/MCMC_BayesC.jl | 8 ++++---- src/1.JWAS/src/MCMC/MCMC_GBLUP.jl | 2 +- src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl | 20 ++++++++++---------- src/1.JWAS/src/MCMC/outputMCMCsamples.jl | 2 +- src/1.JWAS/src/markers/tools4genotypes.jl | 2 +- src/1.JWAS/src/variance.jl | 2 +- src/3.SSBR/src/SSBR.jl | 2 ++ src/5.Datasets/src/Datasets.jl | 1 + 9 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/1.JWAS/src/MCMC/DRY.jl b/src/1.JWAS/src/MCMC/DRY.jl index 7b55a95f..37d7a61d 100644 --- a/src/1.JWAS/src/MCMC/DRY.jl +++ b/src/1.JWAS/src/MCMC/DRY.jl @@ -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 diff --git a/src/1.JWAS/src/MCMC/MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MCMC_BayesC.jl index ade2c2ff..2bee3527 100644 --- a/src/1.JWAS/src/MCMC/MCMC_BayesC.jl +++ b/src/1.JWAS/src/MCMC/MCMC_BayesC.jl @@ -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 @@ -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 diff --git a/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl b/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl index 00d77537..0c84bf2b 100644 --- a/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl +++ b/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl @@ -144,7 +144,7 @@ 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 diff --git a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl index e4ad2064..60b220c9 100644 --- a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl +++ b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl @@ -163,7 +163,7 @@ function MT_MCMC_BayesC(nIter,mme,df; if mme.M != 0 ycorr[:] = ycorr[:] - mme.X*sol iR0 = inv(mme.R) - iGM = (methods!="BayesB"?inv(mme.M.G):[inv(G) for G in arrayG]) + iGM = (methods!="BayesB" ? inv(mme.M.G) : [inv(G) for G in arrayG]) #WILL ADD BURIN INSIDE if methods == "BayesC" sampleMarkerEffectsBayesC!(mArray,mpm,wArray, @@ -201,7 +201,7 @@ function MT_MCMC_BayesC(nIter,mme,df; ######################################################################## # 2.1 Residual Covariance Matrix ######################################################################## - resVec = (mme.M==0?(mme.ySparse - mme.X*sol):ycorr) + resVec = (mme.M==0 ? (mme.ySparse - mme.X*sol) : ycorr) #here resVec is alias for ycor *** if missing_phenotypes==true @@ -257,7 +257,7 @@ function MT_MCMC_BayesC(nIter,mme,df; ycorr[:] = ycorr[:] + X*sol #same to ycorr[:]=resVec+X*sol end - mme.mmeRhs = (mme.M == 0?(X'Ri*mme.ySparse):(X'Ri*ycorr)) + mme.mmeRhs = (mme.M == 0 ? (X'Ri*mme.ySparse) : (X'Ri*ycorr)) ######################################################################## # 2.2 Genetic Covariance Matrix (Polygenic Effects) @@ -344,14 +344,14 @@ function MT_MCMC_BayesC(nIter,mme,df; if output_samples_frequency != 0 && (iter-burnin)%output_samples_frequency==0 && iter>burnin if mme.M != 0 if methods in ["BayesC","BayesCC"] - out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0?G0:false),BigPi,uArray,vec(mme.M.G),outfile) + out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0 ? G0 : false),BigPi,uArray,vec(mme.M.G),outfile) elseif methods == "RR-BLUP" - out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0?G0:false),BigPi,alphaArray,vec(mme.M.G),outfile) + out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0 ? G0 : false),BigPi,alphaArray,vec(mme.M.G),outfile) elseif methods == "BayesB" - out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0?G0:false),BigPi,uArray,false,outfile) + out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0 ? G0 : false),BigPi,uArray,false,outfile) end else - out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0?G0:false),false,false,false,outfile) + out_i=output_MCMC_samples(mme,out_i,sol,R0,(mme.pedTrmVec!=0 ? G0 : false),false,false,false,outfile) end end @@ -387,13 +387,13 @@ function MT_MCMC_BayesC(nIter,mme,df; end end if mme.M != 0 && methods == "RR-BLUP" - output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0?G0Mean:false),output_samples_frequency, + output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0 ? G0Mean : false),output_samples_frequency, meanAlphaArray,GMMean,estimatePi,false,output_file) elseif mme.M != 0 - output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0?G0Mean:false),output_samples_frequency, + output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0 ? G0Mean : false),output_samples_frequency, meanuArray,GMMean,estimatePi,BigPiMean,output_file) else - output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0?G0Mean:false),output_samples_frequency, + output=output_result(mme,solMean,R0Mean,(mme.pedTrmVec!=0 ? G0Mean : false),output_samples_frequency, false,false,false,false,output_file) end diff --git a/src/1.JWAS/src/MCMC/outputMCMCsamples.jl b/src/1.JWAS/src/MCMC/outputMCMCsamples.jl index 4cc24be9..442ff35f 100644 --- a/src/1.JWAS/src/MCMC/outputMCMCsamples.jl +++ b/src/1.JWAS/src/MCMC/outputMCMCsamples.jl @@ -126,7 +126,7 @@ function output_MCMC_samples(mme,out_i,sol,vRes,G0, writedlm(outfile[trmStri*"_variances"],vec(inv(effect.Gi))',',') end - writedlm(outfile["residual_variance"],issubtype(typeof(vRes),Number)?vRes:vec(vRes)',',') + writedlm(outfile["residual_variance"],issubtype(typeof(vRes),Number) ? vRes : vec(vRes)' ,',') #mme.samples4R[out_i,:]=vRes if mme.pedTrmVec != 0 diff --git a/src/1.JWAS/src/markers/tools4genotypes.jl b/src/1.JWAS/src/markers/tools4genotypes.jl index 0c2ef9dc..2c987515 100644 --- a/src/1.JWAS/src/markers/tools4genotypes.jl +++ b/src/1.JWAS/src/markers/tools4genotypes.jl @@ -37,7 +37,7 @@ function getXpRinvX(X) return XpRinvX end -type GibbsMats +mutable struct GibbsMats X::Array{Float64,2} nrows::Int64 ncols::Int64 diff --git a/src/1.JWAS/src/variance.jl b/src/1.JWAS/src/variance.jl index f1aef700..4162a283 100644 --- a/src/1.JWAS/src/variance.jl +++ b/src/1.JWAS/src/variance.jl @@ -22,7 +22,7 @@ end function sampleVCs(mme::MME,sol::Array{Float64,1}) for random_term in mme.rndTrmVec term_array = random_term.term_array - Vi = (random_term.Vinv!=0)?random_term.Vinv:speye(mme.modelTermDict[term_array[1]].nLevels) + Vi = (random_term.Vinv!=0) ? random_term.Vinv : speye(mme.modelTermDict[term_array[1]].nLevels) S = zeros(length(term_array),length(term_array)) for (i,termi) = enumerate(term_array) randTrmi = mme.modelTermDict[termi] diff --git a/src/3.SSBR/src/SSBR.jl b/src/3.SSBR/src/SSBR.jl index 6fffd8d9..4c93d61b 100644 --- a/src/3.SSBR/src/SSBR.jl +++ b/src/3.SSBR/src/SSBR.jl @@ -3,6 +3,8 @@ module SSBR using ..misc using ..PedModule using DataFrames +using SparseArrays +using Printf include("SSBR_types.jl") include("getMatrices.jl") diff --git a/src/5.Datasets/src/Datasets.jl b/src/5.Datasets/src/Datasets.jl index dc8d38ea..7de8b249 100644 --- a/src/5.Datasets/src/Datasets.jl +++ b/src/5.Datasets/src/Datasets.jl @@ -1,5 +1,6 @@ module Datasets +using Printf function dataset(dataset_name::AbstractString,file_name::AbstractString) basename = joinpath(dirname(@__FILE__), "..", "data", dataset_name) From ec594fc7944b9070549f65aeb78613dc17892950 Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 01:29:57 -0700 Subject: [PATCH 3/8] fix warnings in julia0.7 --- src/1.JWAS/src/JWAS.jl | 2 + src/1.JWAS/src/MCMC/DRY.jl | 20 +++++----- src/1.JWAS/src/MCMC/MCMC.jl | 47 ++++++++++++----------- src/1.JWAS/src/MCMC/MCMC_BayesC.jl | 12 +++--- src/1.JWAS/src/MCMC/MCMC_GBLUP.jl | 8 ++-- src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl | 27 +++++++------ src/1.JWAS/src/MCMC/outputMCMCsamples.jl | 12 +++--- src/1.JWAS/src/build_MME.jl | 16 ++++---- src/1.JWAS/src/markers/readgenotypes.jl | 6 +-- src/1.JWAS/src/markers/tools4genotypes.jl | 4 +- src/1.JWAS/src/printout.jl | 4 +- src/1.JWAS/src/random_effects.jl | 18 ++++----- src/1.JWAS/src/residual.jl | 6 +-- src/1.JWAS/src/variance.jl | 2 +- src/2.PedModule/src/PedModule.jl | 3 +- src/4.misc/src/GWAS.jl | 2 +- src/4.misc/src/stat.jl | 6 +-- 17 files changed, 102 insertions(+), 93 deletions(-) diff --git a/src/1.JWAS/src/JWAS.jl b/src/1.JWAS/src/JWAS.jl index 32761f7c..e385cc3d 100644 --- a/src/1.JWAS/src/JWAS.jl +++ b/src/1.JWAS/src/JWAS.jl @@ -1,6 +1,8 @@ using Distributions,Printf +using DelimitedFiles using DataFrames,CSV using SparseArrays +using LinearAlgebra using ProgressMeter using .PedModule using .misc diff --git a/src/1.JWAS/src/MCMC/DRY.jl b/src/1.JWAS/src/MCMC/DRY.jl index 37d7a61d..48ad9ebe 100644 --- a/src/1.JWAS/src/MCMC/DRY.jl +++ b/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 @@ -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 diff --git a/src/1.JWAS/src/MCMC/MCMC.jl b/src/1.JWAS/src/MCMC/MCMC.jl index dcdc4f8c..6b99551f 100644 --- a/src/1.JWAS/src/MCMC/MCMC.jl +++ b/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)"] @@ -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") diff --git a/src/1.JWAS/src/MCMC/MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MCMC_BayesC.jl index 2bee3527..ad91132c 100644 --- a/src/1.JWAS/src/MCMC/MCMC_BayesC.jl +++ b/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 @@ -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 diff --git a/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl b/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl index 0c84bf2b..5940a563 100644 --- a/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl +++ b/src/1.JWAS/src/MCMC/MCMC_GBLUP.jl @@ -151,12 +151,12 @@ function MCMC_GBLUP(nIter,mme,df; ######################################################################## 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 diff --git a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl index 60b220c9..3381903c 100644 --- a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl +++ b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl @@ -16,7 +16,7 @@ function MT_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 @@ -101,18 +101,18 @@ function MT_MCMC_BayesC(nIter,mme,df; ######################################################################## ##WORKING VECTORS ######################################################################## - ycorr = vec(full(mme.ySparse)) #ycorr for different traits - wArray = Array{Array{Float64,1}}(nTraits)#wArray is list reference of ycor + ycorr = vec(Matrix(mme.ySparse)) #ycorr for different traits + wArray = Array{Array{Float64,1}}(undef,nTraits)#wArray is list reference of ycor ######################################################################## #Arrays to save solutions for marker effects ######################################################################## #starting values for marker effects(zeros) and location parameters (sol) - alphaArray = Array{Array{Float64,1}}(nTraits) #BayesC,BayesC0 - meanAlphaArray = Array{Array{Float64,1}}(nTraits) #BayesC,BayesC0 - deltaArray = Array{Array{Float64,1}}(nTraits) #BayesC - meanDeltaArray = Array{Array{Float64,1}}(nTraits) #BayesC - uArray = Array{Array{Float64,1}}(nTraits) #BayesC - meanuArray = Array{Array{Float64,1}}(nTraits) #BayesC + alphaArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC,BayesC0 + meanAlphaArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC,BayesC0 + deltaArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC + meanDeltaArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC + uArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC + meanuArray = Array{Array{Float64,1}}(undef,nTraits) #BayesC for traiti = 1:nTraits startPosi = (traiti-1)*nObs + 1 @@ -153,7 +153,10 @@ function MT_MCMC_BayesC(nIter,mme,df; ######################################################################## # 1.1. Non-Marker Location Parameters ######################################################################## + println("HAHAHA ",iter) Gibbs(mme.mmeLhs,sol,mme.mmeRhs) + println("WAWAWA ",iter) + if iter > burnin solMean += (sol - solMean)/(iter-burnin) end @@ -357,15 +360,15 @@ function MT_MCMC_BayesC(nIter,mme,df; if iter%outFreq==0 && iter>burnin println("\nPosterior means at iteration: ",iter) - println("Residual covariance matrix: \n",round.(R0Mean,6)) + println("Residual covariance matrix: \n",round.(R0Mean,digits=6)) if mme.pedTrmVec !=0 - println("Polygenic effects covariance matrix \n",round.(G0Mean,6)) + println("Polygenic effects covariance matrix \n",round.(G0Mean,digits=6)) end if mme.M != 0 if methods != "BayesB" - println("Marker effects covariance matrix: \n",round.(GMMean,6)) + println("Marker effects covariance matrix: \n",round.(GMMean,digits=6)) end if methods in ["BayesC","BayesB"] && estimatePi == true println("π: \n",BigPiMean) diff --git a/src/1.JWAS/src/MCMC/outputMCMCsamples.jl b/src/1.JWAS/src/MCMC/outputMCMCsamples.jl index 442ff35f..bcc6d6dd 100644 --- a/src/1.JWAS/src/MCMC/outputMCMCsamples.jl +++ b/src/1.JWAS/src/MCMC/outputMCMCsamples.jl @@ -27,7 +27,7 @@ function outputSamplesFor(mme::MME,trmStr::AbstractString) for trmStr in res trm = mme.modelTermDict[trmStr] - samples = MCMCSamples(trm,Array{Float64}(1,1)) + samples = MCMCSamples(trm,Array{Float64}(undef,1,1)) push!(mme.outputSamplesVec,samples) end end @@ -64,9 +64,9 @@ function output_MCMC_samples_setup(mme,nIter,output_samples_frequency,file_name= for i in outvar file_i = file_name*"_"*i*".txt" if isfile(file_i) - warn("The file "*file_i*" already exists!!! It is overwritten by the new output.") + printstyled("The file "*file_i*" already exists!!! It is overwritten by the new output.\n",bold=false,color=:red) else - info("The file "*file_i*" is created to save MCMC samples for "*i*".") + printstyled("The file "*file_i*" is created to save MCMC samples for "*i*".\n",bold=false,color=:green) end outfile[i]=open(file_i,"w") end @@ -126,7 +126,7 @@ function output_MCMC_samples(mme,out_i,sol,vRes,G0, writedlm(outfile[trmStri*"_variances"],vec(inv(effect.Gi))',',') end - writedlm(outfile["residual_variance"],issubtype(typeof(vRes),Number) ? vRes : vec(vRes)' ,',') + writedlm(outfile["residual_variance"],(typeof(vRes) <: Number) ? vRes : vec(vRes)' ,',') #mme.samples4R[out_i,:]=vRes if mme.pedTrmVec != 0 @@ -145,7 +145,7 @@ function output_MCMC_samples(mme,out_i,sol,vRes,G0, writedlm(outfile["marker_effects_variances"],locusEffectVar',',') end writedlm(outfile["pi"],π,',') - if !issubtype(typeof(π),Number)#add a blank line + if !(typeof(π) <: Number) #add a blank line println(outfile["pi"]) end end @@ -168,7 +168,7 @@ end # and Array{SubString{String} for issue 8 function transubstrarr(vec) lvec=length(vec) - res =Array{String}(1,lvec) + res =Array{String}(undef,1,lvec) for i in 1:lvec res[1,i]=vec[i] end diff --git a/src/1.JWAS/src/build_MME.jl b/src/1.JWAS/src/build_MME.jl index 2595ac7c..8bc05eb1 100644 --- a/src/1.JWAS/src/build_MME.jl +++ b/src/1.JWAS/src/build_MME.jl @@ -3,7 +3,7 @@ function mkDict(a) aUnique = unique(a) d = Dict() - names = Array{Any}(size(aUnique,1)) + names = Array{Any}(undef,size(aUnique,1)) for (i,s) in enumerate(aUnique) names[i] = s d[s] = i @@ -39,7 +39,7 @@ function build_model(model_equations::AbstractString,R;df=4) To find an example, type ?build_model and press enter.\n") end - modelVec = [strip(i) for i in split(model_equations,[';','\n'],keep=false)] + modelVec = [strip(i) for i in split(model_equations,[';','\n'],keepempty=false)] nModels = size(modelVec,1) lhsVec = Symbol[] #:y, phenotypes modelTerms = ModelTerm[] #initialization outside for loop @@ -72,7 +72,7 @@ set_covariate(models,"age year") function set_covariate(mme::MME,covStr::AbstractString...) covVec=[] for i in covStr - covVec = [covVec;split(i," ",keep=false)] + covVec = [covVec;split(i," ",keepempty=false)] end mme.covVec = [mme.covVec;[Symbol(i) for i in covVec]] end @@ -84,8 +84,8 @@ end function getData(trm::ModelTerm,df::DataFrame,mme::MME) #ModelTerm("1:A*B") nObs = size(df,1) - trm.str = Array{AbstractString}(nObs) - trm.val = Array{Float64}(nObs) + trm.str = Array{AbstractString}(undef,nObs) + trm.val = Array{Float64}(undef,nObs) if trm.factors[1] == :intercept #for intercept str = fill("intercept",nObs) @@ -127,7 +127,7 @@ getFactor(str) = [strip(i) for i in split(str,"*")] function getX(trm::ModelTerm,mme::MME) #Row Index nObs = length(trm.str) - xi = (trm.iModel-1)*nObs + collect(1:nObs) + xi = (trm.iModel-1)*nObs .+ collect(1:nObs) #Value xv = trm.val #Column Index @@ -185,7 +185,7 @@ function getX(trm::ModelTerm,mme::MME) trm.nLevels = length(dict) dict["0"] = 1 #for missing data xj = round.(Int64,[dict[i] for i in trm.str]) #column index - xv[trm.str.=="0"]=0 #for missing data + xv[trm.str.=="0"] .= 0 #for missing data xi = [xi;1] # adding a zero to xj = [xj;trm.nLevels] # the last column in row 1 @@ -237,7 +237,7 @@ function getMME(mme::MME, df::DataFrame) y = [y; recode(df[mme.lhsVec[i]],missing=>0.0)] end ii = 1:length(y) - jj = ones(ii) + jj = ones(length(y)) vv = y ySparse = sparse(ii,jj,vv) diff --git a/src/1.JWAS/src/markers/readgenotypes.jl b/src/1.JWAS/src/markers/readgenotypes.jl index 8b704ecb..a7cfad80 100644 --- a/src/1.JWAS/src/markers/readgenotypes.jl +++ b/src/1.JWAS/src/markers/readgenotypes.jl @@ -3,7 +3,7 @@ function readgenotypes(file::AbstractString;separator=' ',header=false,center=tr myfile = open(file) #get number of columns - row1 = split(readline(myfile),[separator,'\n'],keep=false) + row1 = split(readline(myfile),[separator,'\n'],keepempty=false) if header==true markerID=row1[2:end] #skip header @@ -41,7 +41,7 @@ function readgenotypes(file::AbstractString;separator=' ',header=false,center=tr markerMeans = center!(copy(genotypes)) #get marker means end p = markerMeans/2.0 - sum2pq = (2*p*(1-p)')[1,1] + sum2pq = (2*p*(1 .- p)')[1,1] return Genotypes(obsID,markerID,nObs,nMarkers,p,sum2pq,center,genotypes) end @@ -60,7 +60,7 @@ function readgenotypes(M::Union{Array{Float64,2},DataFrames.DataFrame};header=fa markerMeans = center!(copy(genotypes)) #get marker means end p = markerMeans/2.0 - sum2pq = (2*p*(1-p)')[1,1] + sum2pq = (2*p*(1 .- p)')[1,1] return Genotypes(obsID,markerID,nObs,nMarkers,p,sum2pq,center,genotypes) end diff --git a/src/1.JWAS/src/markers/tools4genotypes.jl b/src/1.JWAS/src/markers/tools4genotypes.jl index 2c987515..3a9340d7 100644 --- a/src/1.JWAS/src/markers/tools4genotypes.jl +++ b/src/1.JWAS/src/markers/tools4genotypes.jl @@ -12,7 +12,7 @@ end function get_column_ref(X) ncol = size(X)[2] - xArray = Array{Array{Float64,1}}(ncol) + xArray = Array{Array{Float64,1}}(undef,ncol) for i=1:ncol xArray[i] = get_column(X,i) end @@ -21,7 +21,7 @@ end function center!(X) nrow,ncol = size(X) - colMeans = mean(X,1) + colMeans = mean(X,dims=1) BLAS.axpy!(-1,ones(nrow)*colMeans,X) return colMeans end diff --git a/src/1.JWAS/src/printout.jl b/src/1.JWAS/src/printout.jl index fd759422..bdcd810b 100644 --- a/src/1.JWAS/src/printout.jl +++ b/src/1.JWAS/src/printout.jl @@ -2,7 +2,7 @@ #Get names for variables in order ################################################################################ function getNames(mme::MME) - names = Array{AbstractString}(0) + names = Array{AbstractString}(undef,0) for trm in mme.modelTerms for name in trm.names push!(names,trm.trmStr*" : "*name) @@ -12,7 +12,7 @@ function getNames(mme::MME) end function getNames(trm::ModelTerm) - names = Array{AbstractString}(0) + names = Array{AbstractString}(undef,0) for name in trm.names push!(names,trm.trmStr*" : "*name) end diff --git a/src/1.JWAS/src/random_effects.jl b/src/1.JWAS/src/random_effects.jl index de7e1a43..ee2b5cbc 100644 --- a/src/1.JWAS/src/random_effects.jl +++ b/src/1.JWAS/src/random_effects.jl @@ -50,7 +50,7 @@ set_random(model,"Animal", ped,G) ``` """ function set_random(mme::MME,randomStr::AbstractString,ped::PedModule.Pedigree, G;df=4,output_samples=true) - pedTrmVec = split(randomStr," ",keep=false) # "animal animal*age" + pedTrmVec = split(randomStr," ",keepempty=false) # "animal animal*age" #add model equation number; "animal" => "1:animal" res = [] @@ -64,7 +64,7 @@ function set_random(mme::MME,randomStr::AbstractString,ped::PedModule.Pedigree, outputMCMCsamples(mme,trm) #output MCMC samples (used to calculate EBV,PEV) end else - info(trm," is not found in model equation ",string(m),".") + printstyled(trm," is not found in model equation ",string(m),".\n",bold=false,color=:red) end end end #"1:animal","1:animal*age" @@ -81,10 +81,10 @@ function set_random(mme::MME,randomStr::AbstractString,ped::PedModule.Pedigree, if mme.nModels!=1 #multi-trait mme.Gi = inv(G) else #single-trait - if issubtype(typeof(G),Number)==true #convert scalar G to 1x1 matrix + if (typeof(G)<:Number) ==true #convert scalar G to 1x1 matrix G=reshape([G],1,1) end - mme.GiOld = zeros(G) + mme.GiOld = zero(G) mme.GiNew = inv(G) end mme.df.polygenic=Float64(df) @@ -119,7 +119,7 @@ set_random(model,"litter",G) function set_random(mme::MME,randomStr::AbstractString,G;Vinv=0,names=[],df=4) G = map(Float64,G) df= Float64(df) - randTrmVec = split(randomStr," ",keep=false) # "herd" + randTrmVec = split(randomStr," ",keepempty=false) # "herd" for trm in randTrmVec res = [] for (m,model) = enumerate(mme.modelVec) @@ -133,13 +133,13 @@ function set_random(mme::MME,randomStr::AbstractString,G;Vinv=0,names=[],df=4) mme.modelTermDict[mtrm].names=names #********************* else - info(trm," is not found in model equation ",string(m),".") + printstyled(trm," is not found in model equation ",string(m),".\n",bold=false,color=:red) end end if length(res) != size(G,1) error("Dimensions must match. The covariance matrix (G) should be a ",length(res)," x ",length(res)," matrix.\n") end - if issubtype(typeof(G),Number)==true #convert scalar G to 1x1 matrix + if (typeof(G)<:Number) ==true #convert scalar G to 1x1 matrix #Need (here in julia0.7) G=reshape([G],1,1) end Gi = inv(G) @@ -147,7 +147,7 @@ function set_random(mme::MME,randomStr::AbstractString,G;Vinv=0,names=[],df=4) term_array = res df = df+length(term_array) scale = G*(df-length(term_array)-1) #G*(df-2)/df #from inv χ to inv-wishat - randomEffect = RandomEffect(term_array,Gi,zeros(Gi),Gi,df,scale,Vinv,names) + randomEffect = RandomEffect(term_array,Gi,zero(Gi),Gi,df,scale,Vinv,names) push!(mme.rndTrmVec,randomEffect) end nothing @@ -200,7 +200,7 @@ end function addLambdas(mme::MME) for random_term in mme.rndTrmVec term_array = random_term.term_array - Vi = (random_term.Vinv!=0) ? random_term.Vinv : speye(mme.modelTermDict[term_array[1]].nLevels) + Vi = (random_term.Vinv!=0) ? random_term.Vinv : SparseMatrixCSC{Float64}(I, mme.modelTermDict[term_array[1]].nLevels, mme.modelTermDict[term_array[1]].nLevels) for (i,termi) = enumerate(term_array) randTrmi = mme.modelTermDict[termi] startPosi = randTrmi.startPos diff --git a/src/1.JWAS/src/residual.jl b/src/1.JWAS/src/residual.jl index c9536bc1..387c9c49 100644 --- a/src/1.JWAS/src/residual.jl +++ b/src/1.JWAS/src/residual.jl @@ -24,9 +24,9 @@ function mkRi(mme::MME,df::DataFrame) mme.missingPattern = tstMsng n = size(tstMsng,2) nObs = size(tstMsng,1) - ii = Array{Int64}(nObs*n^2) - jj = Array{Int64}(nObs*n^2) - vv = Array{Float64}(nObs*n^2) + ii = Array{Int64}(undef,nObs*n^2) + jj = Array{Int64}(undef,nObs*n^2) + vv = Array{Float64}(undef,nObs*n^2) pos = 1 for i=1:size(tstMsng,1) sel = reshape(tstMsng[i,:],n) diff --git a/src/1.JWAS/src/variance.jl b/src/1.JWAS/src/variance.jl index 4162a283..4b7f62ae 100644 --- a/src/1.JWAS/src/variance.jl +++ b/src/1.JWAS/src/variance.jl @@ -22,7 +22,7 @@ end function sampleVCs(mme::MME,sol::Array{Float64,1}) for random_term in mme.rndTrmVec term_array = random_term.term_array - Vi = (random_term.Vinv!=0) ? random_term.Vinv : speye(mme.modelTermDict[term_array[1]].nLevels) + Vi = (random_term.Vinv!=0) ? random_term.Vinv : SparseMatrixCSC{Float64}(I, mme.modelTermDict[term_array[1]].nLevels, mme.modelTermDict[term_array[1]].nLevels) S = zeros(length(term_array),length(term_array)) for (i,termi) = enumerate(term_array) randTrmi = mme.modelTermDict[termi] diff --git a/src/2.PedModule/src/PedModule.jl b/src/2.PedModule/src/PedModule.jl index 041fbd99..2856def9 100644 --- a/src/2.PedModule/src/PedModule.jl +++ b/src/2.PedModule/src/PedModule.jl @@ -1,6 +1,7 @@ module PedModule using DataFrames,CSV +using SparseArrays using ProgressMeter mutable struct PedNode @@ -179,7 +180,7 @@ end function getIDs(ped::Pedigree) n = length(ped.idMap) - ids = Array{String}(n) + ids = Array{String}(undef,n) for i in ped.idMap ids[i[2].seqID] = i[1] end diff --git a/src/4.misc/src/GWAS.jl b/src/4.misc/src/GWAS.jl index d770eb7e..09974659 100644 --- a/src/4.misc/src/GWAS.jl +++ b/src/4.misc/src/GWAS.jl @@ -125,7 +125,7 @@ function GWAS(marker_effects_file,map_file,mme;header=false,window_size="1 Mb",t winVarProps = GWAS(marker_effects_file,mme,header=header,window_size=window_size,threshold=threshold) winVarProps[isnan.(winVarProps)]=0.0 #replace NaN caused by situations no markers are included in the model WPPA, prop_genvar = vec(mean(winVarProps .> threshold,1)), vec(mean(winVarProps,1)) - prop_genvar = round.(prop_genvar*100,2) + prop_genvar = round.(prop_genvar*100,digits=2) #bug in Julia, vcat too long #out = [["window";1:length(WPPA)] ["chr"; window_chr] ["start"; window_pos_start] ["end"; window_pos_end]["WPPA"; WPPA]] out1 = [["window";1:length(WPPA)] ["chr"; window_chr]] diff --git a/src/4.misc/src/stat.jl b/src/4.misc/src/stat.jl index 652a213a..633a4248 100644 --- a/src/4.misc/src/stat.jl +++ b/src/4.misc/src/stat.jl @@ -7,8 +7,8 @@ function report(X::Array{Array{Float64,2},1};index=false) Xmean = mean(X) Xvar = mean(X.^2)-Xmean.^2 println("Summary Stats:") - println("Mean:\n",round(Xmean,6)) - println("Variance:\n",round(Xvar,6)) + println("Mean:\n",round(Xmean,digits=6)) + println("Variance:\n",round(Xvar,digits=6)) else xArray = zeros(length(X)) i=1 @@ -30,7 +30,7 @@ end function report(X::Array{Array{Float64,1},1};index=false) if index==false println("Summary Stats:") - println("Mean:\n",round(mean(X),6)) + println("Mean:\n",round(mean(X),digits=6)) else xArray = zeros(length(X)) i=1 From 93975ea7244e48a54e9f19d1a1ba5233fdf54a2d Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 02:33:57 -0700 Subject: [PATCH 4/8] add error message for non positive definite covariance matrix --- src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl | 2 -- src/1.JWAS/src/build_MME.jl | 4 ++++ src/1.JWAS/src/random_effects.jl | 7 +++++++ 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl index 3381903c..015b3e02 100644 --- a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl +++ b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl @@ -153,9 +153,7 @@ function MT_MCMC_BayesC(nIter,mme,df; ######################################################################## # 1.1. Non-Marker Location Parameters ######################################################################## - println("HAHAHA ",iter) Gibbs(mme.mmeLhs,sol,mme.mmeRhs) - println("WAWAWA ",iter) if iter > burnin solMean += (sol - solMean)/(iter-burnin) diff --git a/src/1.JWAS/src/build_MME.jl b/src/1.JWAS/src/build_MME.jl index 8bc05eb1..1043d471 100644 --- a/src/1.JWAS/src/build_MME.jl +++ b/src/1.JWAS/src/build_MME.jl @@ -34,6 +34,10 @@ models = build_model(model_equations,R); ``` """ function build_model(model_equations::AbstractString,R;df=4) + if !isposdef(R) + @error "The covariance matrix is not positive definite." + end + if !(typeof(model_equations)<:AbstractString) || model_equations=="" error("Model equations are wrong.\n To find an example, type ?build_model and press enter.\n") diff --git a/src/1.JWAS/src/random_effects.jl b/src/1.JWAS/src/random_effects.jl index ee2b5cbc..235e20dc 100644 --- a/src/1.JWAS/src/random_effects.jl +++ b/src/1.JWAS/src/random_effects.jl @@ -50,6 +50,10 @@ set_random(model,"Animal", ped,G) ``` """ function set_random(mme::MME,randomStr::AbstractString,ped::PedModule.Pedigree, G;df=4,output_samples=true) + if !isposdef(G) + @error "The covariance matrix is not positive definite." + end + pedTrmVec = split(randomStr," ",keepempty=false) # "animal animal*age" #add model equation number; "animal" => "1:animal" @@ -117,6 +121,9 @@ set_random(model,"litter",G) ``` """ function set_random(mme::MME,randomStr::AbstractString,G;Vinv=0,names=[],df=4) + if !isposdef(G) + @error "The covariance matrix is not positive definite." + end G = map(Float64,G) df= Float64(df) randTrmVec = split(randomStr," ",keepempty=false) # "herd" From 85ba9418d29b627482ee19fb70c4d568bc6302ed Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 19:23:21 -0700 Subject: [PATCH 5/8] fix warnings in SSBR for julia0.7 --- src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl | 2 +- src/1.JWAS/src/Pi.jl | 2 +- src/1.JWAS/src/SSBR/SSBR.jl | 6 +++--- src/1.JWAS/src/markers/BayesianAlphabet/MTBayesC.jl | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl index 015b3e02..dd6257ef 100644 --- a/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl +++ b/src/1.JWAS/src/MCMC/MT_MCMC_BayesC.jl @@ -235,7 +235,7 @@ function MT_MCMC_BayesC(nIter,mme,df; if mme.M != 0 mme.R = R0 R0 = mme.R - Ri = kron(inv(R0),speye(nObs)) + Ri = kron(inv(R0),SparseMatrixCSC{Float64}(I, nObs, nObs)) RiNotUsing = mkRi(mme,df) #get small Ri (Resvar) used in imputation end diff --git a/src/1.JWAS/src/Pi.jl b/src/1.JWAS/src/Pi.jl index 988c8f0b..09283be7 100644 --- a/src/1.JWAS/src/Pi.jl +++ b/src/1.JWAS/src/Pi.jl @@ -9,7 +9,7 @@ function samplePi(deltaArray,BigPi,BigPiMean,iter) nLoci_array=zeros(2^nTraits) for i in keys(BigPi) #assume order of key won't change temp2 = broadcast(-,temp,i') - nLoci = sum(mean(abs.(temp2),2).==0.0) + nLoci = sum(mean(abs.(temp2),dims=2).==0.0) nLoci_array[iloci] = nLoci +1 iloci = iloci +1 end diff --git a/src/1.JWAS/src/SSBR/SSBR.jl b/src/1.JWAS/src/SSBR/SSBR.jl index b8108ff9..ecdd216a 100644 --- a/src/1.JWAS/src/SSBR/SSBR.jl +++ b/src/1.JWAS/src/SSBR/SSBR.jl @@ -48,8 +48,8 @@ end # Genotypes ############################################################################ function impute_genotypes(geno,ped,mme,Ai_nn,Ai_ng) - Mg = Array{Float64,2}(geno.nObs,geno.nMarkers) - MgID = Array{String,1}(geno.nObs) + Mg = Array{Float64,2}(undef,geno.nObs,geno.nMarkers) + MgID = Array{String,1}(undef,geno.nObs) num_pedn = size(Ai_nn,1) #reorder genotypes to get Mg with the same order as Ai_gg for i in 1:geno.nObs #Better to use Z matrix? @@ -65,7 +65,7 @@ function impute_genotypes(geno,ped,mme,Ai_nn,Ai_ng) mme.M.genotypes = Mfull mme.M.obsID = IDs mme.M.nObs = length(mme.M.obsID) - gc() + GC.gc() end ############################################################################ # Fixed effects (J) diff --git a/src/1.JWAS/src/markers/BayesianAlphabet/MTBayesC.jl b/src/1.JWAS/src/markers/BayesianAlphabet/MTBayesC.jl index b8d79ebc..5cea0ddb 100644 --- a/src/1.JWAS/src/markers/BayesianAlphabet/MTBayesC.jl +++ b/src/1.JWAS/src/markers/BayesianAlphabet/MTBayesC.jl @@ -30,7 +30,7 @@ function sampleMarkerEffectsBayesC!(xArray,xpx,wArray,alphaArray,meanAlphaArray, nok = deleteat!(collect(1:nTraits),k) Ginv12 = Ginv[k,nok] C11 = Ginv11+Rinv[k,k]*xpx[marker] - C12 = Ginv12+xpx[marker]*diagm(δ[nok])*Rinv[k,nok] + C12 = Ginv12+xpx[marker]*Matrix(Diagonal(δ[nok]))*Rinv[k,nok] #C12 = Ginv12+xpx[marker]*Rinv[k,nok].*δ[nok]' #δ[:,nok] : row vector, invLhs0 = 1/Ginv11 From 77c7ff140f8e46efc422779c0f5fc6f75a6299ab Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 22:54:46 -0700 Subject: [PATCH 6/8] update require --- REQUIRE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/REQUIRE b/REQUIRE index 3b465494..a683c4d2 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,4 @@ -julia 0.6 +julia 0.7 Distributions DataFrames ProgressMeter From 46f4d08eb3bf2249e171501552d15f916f40704a Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 23:01:13 -0700 Subject: [PATCH 7/8] update travis --- .travis.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.travis.yml b/.travis.yml index 15fe2cbe..420b007f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,10 +4,13 @@ os: - osx julia: - 0.6 + - 0.7 + - 1.0 - nightly matrix: allow_failures: - julia: nightly + - julia: 0.6 notifications: email: false after_success: From ad11154a6752a0e657fe95179cccb579deb1a316 Mon Sep 17 00:00:00 2001 From: haocheng Date: Sun, 19 Aug 2018 23:17:09 -0700 Subject: [PATCH 8/8] update travis --- .travis.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.travis.yml b/.travis.yml index 420b007f..5c5948d9 100644 --- a/.travis.yml +++ b/.travis.yml @@ -10,6 +10,7 @@ julia: matrix: allow_failures: - julia: nightly + - julia: 1.0 - julia: 0.6 notifications: email: false