Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JWAS update #112

Merged
merged 11 commits into from
Feb 22, 2022
Merged

JWAS update #112

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
17e37e0
allow continuous+ordinal
zhaotianjing Feb 12, 2022
eadf606
new interface for single and multiple censored traits
zhaotianjing Feb 15, 2022
b17740d
Merge branch 'master' into censored
zhaotianjing Feb 15, 2022
04c4621
allow mix traits
zhaotianjing Feb 17, 2022
82cdffb
Merge branch 'master' into censored
zhaotianjing Feb 17, 2022
872af1c
update
zhaotianjing Feb 17, 2022
0cb943f
clean code
zhaotianjing Feb 20, 2022
f32ef35
clean code again
zhaotianjing Feb 20, 2022
91da1f3
clean code again*2
zhaotianjing Feb 20, 2022
0ddbcab
binary trait
zhaotianjing Feb 21, 2022
bb2c281
print example for deprecated censored and categorical features
zhaotianjing Feb 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 24 additions & 13 deletions src/1.JWAS/src/JWAS.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,8 @@ include("markers/Pi.jl")
include("single_step/SSBR.jl")
include("single_step/SSGBLUP.jl")

#Categorical Traits
include("categorical_trait/categorical_trait.jl")

#Censored Traits
include("censored_trait/censored_trait.jl")
#Categorical and Censored Traits
include("categorical_and_censored_trait/categorical_and_censored_trait.jl")

#Structure Equation Models
include("structure_equation_model/SEM.jl")
Expand Down Expand Up @@ -158,8 +155,6 @@ function runMCMC(mme::MME,df;
single_step_analysis = false, #parameters for single-step analysis
pedigree = false, #parameters for single-step analysis
fitting_J_vector = true, #parameters for single-step analysis
categorical_trait = false,
censored_trait = false,
causal_structure = false,
mega_trait = mme.nonlinear_function == false ? false : true, #NNBayes -> default mega_trait=true
missing_phenotypes = true,
Expand All @@ -181,10 +176,13 @@ function runMCMC(mme::MME,df;
methods = "conventional (no markers)",
Pi = 0.0,
estimatePi = false,
estimateScale = false)

estimateScale = false,
categorical_trait = false, #this has been moved to build_model()
censored_trait = false) #this has been moved to build_model()

#Neural Network
############################################################################
# Neural Network
############################################################################
is_nnbayes_partial = (mme.nonlinear_function != false && mme.is_fully_connected==false)
if mme.nonlinear_function != false #modify data to add phenotypes for hidden nodes
mme.yobs_name=Symbol(mme.lhsVec[1]) #e.g., lhsVec=[:y1,:y2,:y3], a number label has been added to original trait name in nnbayes_model_equation(),
Expand All @@ -207,7 +205,9 @@ function runMCMC(mme::MME,df;
mme.M[1].trait_names=mme.latent_traits
end
end
############################################################################
#for deprecated JWAS fucntions
############################################################################
if mme.M != 0
for Mi in mme.M
if Mi.name == false
Expand All @@ -219,6 +219,18 @@ function runMCMC(mme::MME,df;
end
end
end
if categorical_trait != false || censored_trait != false
print_single_categorical_censored_trait_example()
error("The arguments 'categorical_trait' and 'censored_trait' has been moved to build_model(). Please check our latest example.")
end
############################################################################
# censored traits
############################################################################
#add the column :traitname using trait's lower bound
censored_trait_index = findall(x -> x=="censored", mme.traits_type)
for i in censored_trait_index
df[!,mme.lhsVec[i]]= df[!,Symbol(mme.lhsVec[i],"_l")]
end
############################################################################
# Set a seed in the random number generator
############################################################################
Expand Down Expand Up @@ -251,7 +263,6 @@ function runMCMC(mme::MME,df;
printout_model_info,printout_frequency, single_step_analysis,
fitting_J_vector,missing_phenotypes,constraint,mega_trait,estimate_variance,
update_priors_frequency,outputEBV,output_heritability,prediction_equation,
categorical_trait,censored_trait,
seed,double_precision,output_folder)
############################################################################
# Check 1)Arguments; 2)Input Pedigree,Genotype,Phenotypes,
Expand Down Expand Up @@ -440,7 +451,7 @@ end
* (internal function) Print out MCMC information.
"""
function getMCMCinfo(mme)
is_nnbayes_partial = mme.nonlinear_function != false && mme.is_fully_connected==false
is_nnbayes_partial = mme.nonlinear_function != false && mme.is_fully_connected==false
if mme.MCMCinfo == false
printstyled("MCMC information is not available\n\n",bold=true)
return
Expand Down Expand Up @@ -474,7 +485,7 @@ function getMCMCinfo(mme)
println()
end
if mme.nModels == 1
@printf("%-30s %20.3f\n","residual variances:", (MCMCinfo.categorical_trait ? 1.0 : mme.R))
@printf("%-30s %20.3f\n","residual variances:", mme.R)
else
@printf("%-30s\n","residual variances:")
Base.print_matrix(stdout,round.(mme.R,digits=3))
Expand Down
56 changes: 22 additions & 34 deletions src/1.JWAS/src/MCMC/MCMC_BayesianAlphabet.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,9 @@ function MCMC_BayesianAlphabet(mme,df)
estimate_variance = mme.MCMCinfo.estimate_variance
invweights = mme.invweights
update_priors_frequency = mme.MCMCinfo.update_priors_frequency
categorical_trait = mme.MCMCinfo.categorical_trait
censored_trait = mme.MCMCinfo.censored_trait
has_categorical_trait = "categorical" ∈ mme.traits_type
has_censored_trait = "censored" ∈ mme.traits_type
categorical_trait_index = findall(x -> x=="categorical", mme.traits_type)
missing_phenotypes = mme.MCMCinfo.missing_phenotypes
constraint = mme.MCMCinfo.constraint
causal_structure = mme.causal_structure
Expand All @@ -23,11 +24,8 @@ function MCMC_BayesianAlphabet(mme,df)
############################################################################
# Categorical Traits (starting values for maker effects defaulting to 0s)
############################################################################
if categorical_trait == true
category_obs,threshold = categorical_trait_setup!(mme)
end
if censored_trait != false
lower_bound,upper_bound = censored_trait_setup!(mme)
if has_categorical_trait || has_censored_trait
category_obs,thresholds,lower_bound,upper_bound = categorical_censored_traits_setup!(mme,df)
end
############################################################################
# Working Variables
Expand All @@ -39,12 +37,9 @@ function MCMC_BayesianAlphabet(mme,df)
#mme.sol (its starting values were set in runMCMC)
mme.solMean, mme.solMean2 = zero(mme.sol),zero(mme.sol)
#residual variance
if categorical_trait == true
mme.R = mme.meanVare = mme.meanVare2 = 1.0
else
mme.meanVare = zero(mme.R)
mme.meanVare2 = zero(mme.R)
end
mme.meanVare = zero(mme.R)
mme.meanVare2 = zero(mme.R)

#polygenic effects (pedigree), e.g, Animal+ Maternal
if mme.pedTrmVec != 0
mme.G0Mean,mme.G0Mean2 = zero(mme.Gi),zero(mme.Gi)
Expand Down Expand Up @@ -156,25 +151,18 @@ function MCMC_BayesianAlphabet(mme,df)
if nonlinear_function != false
mme.weights_NN = vcat(mean(mme.ySparse),zeros(mme.nModels))
end

@showprogress "running MCMC ..." for iter=1:chain_length
########################################################################
# 0. Categorical traits (liabilities)
# 0. Categorical and censored traits
########################################################################
if categorical_trait == true
ycorr = categorical_trait_sample_liabilities(mme,ycorr,category_obs,threshold)
writedlm(outfile["threshold"],threshold',',')
end
if censored_trait != false
ycorr = is_multi_trait ? MT_censored_trait_sample_liabilities(mme,ycorr,lower_bound,upper_bound) : censored_trait_sample_liabilities(mme,ycorr,lower_bound,upper_bound)
if has_categorical_trait || has_censored_trait
sample_liabilities!(mme,ycorr,lower_bound,upper_bound) #update mme.ySparse, ycorr
writedlm(outfile["liabilities"],mme.ySparse',',')
end
#update wArray for multi-trait
if is_multi_trait
for traiti = 1:mme.nModels
startPosi = (traiti-1)*length(mme.obsID) + 1
ptr = pointer(ycorr,startPosi)
wArray[traiti] = unsafe_wrap(Array,ptr,length(mme.obsID))
if has_categorical_trait
#sample threshold for categorical traits
thresholds=categorical_trait_sample_threshold(mme, thresholds, category_obs) #update thresholds
writedlm(outfile["threshold"],vcat(thresholds...),',')
update_bounds_from_threshold!(lower_bound,upper_bound,category_obs,thresholds,categorical_trait_index) # update lower_bound, upper_bound
end
end
########################################################################
Expand Down Expand Up @@ -317,8 +305,8 @@ function MCMC_BayesianAlphabet(mme,df)
mme.df.residual, mme.scaleR,
invweights,constraint)
Ri = kron(inv(mme.R),spdiagm(0=>invweights))
else
if categorical_trait == false # fixed mme.R for categorical_trait
else #single trait
if !has_categorical_trait # fixed mme.R=1 for single categorical trait
mme.ROld = mme.R
mme.R = sample_variance(ycorr,length(ycorr), mme.df.residual, mme.scaleR, invweights)
end
Expand Down Expand Up @@ -388,11 +376,11 @@ function MCMC_BayesianAlphabet(mme,df)
if causal_structure != false
close(causal_structure_outfile)
end
if categorical_trait == true
close(outfile["threshold"])
end
if censored_trait != false
if has_categorical_trait || has_censored_trait
close(outfile["liabilities"])
if has_categorical_trait
close(outfile["threshold"])
end
end
end
if methods == "GBLUP"
Expand Down
14 changes: 12 additions & 2 deletions src/1.JWAS/src/build_MME.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,9 @@ models = build_model(model_equations,R);
```
"""
function build_model(model_equations::AbstractString, R = false; df = 4.0,
num_hidden_nodes = false, nonlinear_function = false, latent_traits=false,
user_σ2_yobs = false, user_σ2_weightsNN = false) #nonlinear_function(x1,x2) = x1+x2
num_hidden_nodes = false, nonlinear_function = false, latent_traits=false, #nonlinear_function(x1,x2) = x1+x2
user_σ2_yobs = false, user_σ2_weightsNN = false,
censored_trait = false, categorical_trait = false)

if R != false && !isposdef(map(AbstractFloat,R))
error("The covariance matrix is not positive definite.")
Expand Down Expand Up @@ -137,6 +138,15 @@ function build_model(model_equations::AbstractString, R = false; df = 4.0,
end
end

#setup traits_type (by default is "continuous")
for t in 1:mme.nModels
if string(mme.lhsVec[t]) ∈ censored_trait
mme.traits_type[t]="censored"
elseif string(mme.lhsVec[t]) ∈ categorical_trait
mme.traits_type[t]="categorical"
end
end

return mme
end

Expand Down
Loading