diff --git a/src/1.JWAS/src/build_MME.jl b/src/1.JWAS/src/build_MME.jl index 74acb7d0..0f8975a1 100644 --- a/src/1.JWAS/src/build_MME.jl +++ b/src/1.JWAS/src/build_MME.jl @@ -267,6 +267,13 @@ function getX(trm::ModelTerm,mme::MME) mme.mmePos += trm.nLevels end +function makeWeights(fw::T, useWeights::Bool, df::DataFrame, column="weights")::Vector{T} where T <: AbstractFloat + res = ones(typeof(fw),size(df,1)) + if useWeights + res = one(fw) ./ map((x)->(convert(typeof(fw),x)), df[:,column]) + end + return res +end """ Construct mixed model equations with @@ -279,13 +286,20 @@ function getMME(mme::MME, df::DataFrame) if mme.mmeLhs != false error("Please build your model again using the function build_model().") end + + use64 = mme.MCMCinfo == false || mme.MCMCinfo.double_precision + + useWeights = mme.MCMCinfo != false && mme.MCMCinfo.heterogeneous_residuals == true + #Heterogeneous residuals - if mme.MCMCinfo != false && mme.MCMCinfo.heterogeneous_residuals == true - invweights = 1 ./ convert(Array,df[!,Symbol("weights")]) + if use64 + invweights = makeWeights(Float64(1.0),useWeights, df, "weights") + mme.invweights =invweights else - invweights = ones(size(df,1)) + invweights = makeWeights(Float32(1.0),useWeights, df, "weights") + mme.invweights =invweights end - mme.invweights = (mme.MCMCinfo == false || mme.MCMCinfo.double_precision ? Float64.(invweights) : Float32.(invweights)) + #Make incidence matrices X for each term for trm in mme.modelTerms if trm.X == false @@ -321,7 +335,7 @@ function getMME(mme::MME, df::DataFrame) #rename genotype names mme.M[1].trait_names=mme.latent_traits #save omics data missing pattern - mme.missingPattern = .!ismissing.(convert(Matrix,df[!,mme.lhsVec])) + mme.missingPattern = .!ismissing.(Matrix(df[!,mme.lhsVec])) #replace missing data with values in yobs for i in mme.lhsVec #for each omics feature for j in 1:size(df,1) #for each observation @@ -339,13 +353,14 @@ function getMME(mme::MME, df::DataFrame) end ii = 1:length(y) jj = ones(length(y)) - vv = ((mme.MCMCinfo == false || mme.MCMCinfo.double_precision) ? Float64.(y) : Float32.(y)) + #vv = ((mme.MCMCinfo == false || mme.MCMCinfo.double_precision) ? Float64.(y) : Float32.(y)) + # This recoding should make it easy for julia to reason about the type of vv + vv::Union{Vector{Float64},Vector{Float32}} = use64 ? Float64.(y) : Float32.(y) ySparse = sparse(ii,jj,vv) #Make lhs and rhs for MME mme.X = X mme.ySparse = ySparse - if mme.nModels==1 #single-trait (lambda version) mme.mmeLhs = X'*Diagonal(mme.invweights)*X mme.mmeRhs = X'*Diagonal(mme.invweights)*ySparse @@ -363,7 +378,6 @@ function getMME(mme::MME, df::DataFrame) mme.mmeLhs = X'Ri*X mme.mmeRhs = X'Ri*ySparse end - #Random effects parts in MME if mme.nModels == 1 #random_term.GiNew*mme.R - random_term.GiOld*mme.ROld @@ -384,8 +398,8 @@ function getMME(mme::MME, df::DataFrame) #No phenotypic data for some levels of a factor in multi-trait analysis #e.g., y3:x3:f in https://github.com/reworkhow/JWAS.jl/blob/ #a6f4595796b70811c0b745af525e7e0a822bb954/src/5.Datasets/data/example/phenotypes.txt - for i in size(mme.mmeLhs,1) - if mme.mmeLhs[i,i] == 0.0 + for (i,d) in enumerate(Array(diag(mme.mmeLhs))) + if iszero(d) error("No phenotypic data for ",getNames(mme)[i]) end end diff --git a/src/1.JWAS/src/markers/readgenotypes.jl b/src/1.JWAS/src/markers/readgenotypes.jl index 1a4b9725..42f87298 100644 --- a/src/1.JWAS/src/markers/readgenotypes.jl +++ b/src/1.JWAS/src/markers/readgenotypes.jl @@ -103,13 +103,13 @@ function get_genotypes(file::Union{AbstractString,Array{Float64,2},Array{Float32 #set type for each column ncol= length(row1) etv = Array{DataType}(undef,ncol) - fill!(etv,Float64) + fill!(etv,Float32) etv[1]=String close(myfile) #read a large genotype file data = CSV.read(file,DataFrame,types=etv,delim = separator,header=false,skipto=(header==true ? 2 : 1)) - obsID = map(string,data[!,1]) - genotypes = map(Float32,Matrix(data[!,2:end])) + obsID = data[!,1] + genotypes = Matrix(data[!,2:end]) elseif typeof(file) == DataFrames.DataFrame #Datafarme println("The first column in the dataframe should be individual IDs.") println("The data type of markers should be Number.") @@ -119,7 +119,7 @@ function get_genotypes(file::Union{AbstractString,Array{Float64,2},Array{Float32 markerID = string.(1:(size(file,2)-1)) end obsID = map(string,file[!,1]) - genotypes = map(Float32,convert(Matrix,file[!,2:end])) + genotypes = map(Float32,Matrix(file[!,2:end])) elseif typeof(file) <: Union{Array{Float64,2},Array{Float32,2},Array{Any,2}} #Array if length(header) != (size(file,2)+1) header = ["id"; string.(1:size(file,2))] diff --git a/src/1.JWAS/src/residual.jl b/src/1.JWAS/src/residual.jl index 1b9e96d0..a1f79675 100644 --- a/src/1.JWAS/src/residual.jl +++ b/src/1.JWAS/src/residual.jl @@ -22,7 +22,7 @@ function mkRi(mme::MME,df::DataFrame,Rinv) nObs = size(tstMsng,1) ii = Array{Int64}(undef,nObs*ntrait^2) jj = Array{Int64}(undef,nObs*ntrait^2) - vv = Array{AbstractFloat}(undef,nObs*ntrait^2) + vv = Array{eltype(Rinv)}(undef,nObs*ntrait^2) pos = 1 for i=1:nObs sel = tstMsng[i,:]