Joint Analysis of Continuous, Categorical, and Censored Traits

Multi-trait analysis (continuous & censored & categorical traits)

For single/multiple continuous trait analysis, please see our example "Multiple Trait Analysis".

For single/multiple censored trait analysis, please see our example "Censored Trait Analysis".

For single/multiple categorical trait analysis, please see our example "Categorical Trait Analysis".

For joint analysis of continuous, censored and categorical traits, the categorical traits and censored traits should be identified in the argument categorical_trait and censored_trait in build_model().

Below is an example for three-trait analysis, with

• y1: censored trait
• y2: categorical trait
• y3: continuous trait

Data simulation

using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets,Random
Random.seed!(1)

phenofile  = dataset("phenotypes.csv")
phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])

#y1 - censored trait
phenotypes[!, :lower_bound1] .= round.(phenotypes[!, :y1],digits=2);
phenotypes[!, :upper_bound1] .= round.(phenotypes[!, :y1],digits=2);
phenotypes[0.0 .< phenotypes[!, :y1] .< 1.0, :upper_bound1] .= phenotypes[0.0 .< phenotypes[!, :y1] .< 1.0, :y1] .+ 0.5;
phenotypes[0.0 .< phenotypes[!, :y1] .< 1.0, :lower_bound1] .= phenotypes[0.0 .< phenotypes[!, :y1] .< 1.0, :y1] .- 0.5;
phenotypes[phenotypes[!, :y1] .> 2.0,:upper_bound1] .= Inf;
phenotypes[phenotypes[!, :y1] .> 2.0,:lower_bound1] .= 2.0;
phenotypes[phenotypes[!, :y1] .< -3.0,:upper_bound1].= -3.0;
phenotypes[phenotypes[!, :y1] .< -3.0,:lower_bound1].= -Inf;
phenotypes=phenotypes[:, Not(:y1)]  #drop :y1 column

#y2 -  categorical trait
category1_sel = phenotypes[!,:y2] .< -1.7
category2_sel = -1.7 .< phenotypes[!,:y2] .< 0.0
category3_sel = phenotypes[!,:y2] .> 0.0
phenotypes[category1_sel, :y2] .= 1
phenotypes[category2_sel, :y2] .= 2
phenotypes[category3_sel, :y2] .= 3

JWAS analysis:

# Step 0: Load packages
using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets,Random
Random.seed!(1)

# Step 1: for censored traits, rename the lower bound as "traitname_l", and upper bound as "traitname_u"
rename!(phenotypes, :lower_bound1 => :y1_l, :upper_bound1 => :y1_u)

# Step 2: Read data
pedfile    = dataset("pedigree.csv")
pedigree   = get_pedigree(pedfile,separator=",",header=true);
genofile   = dataset("genotypes.csv")
genotypes  = get_genotypes(genofile,separator=',',method="BayesC");


# Step 3: Build Model Equations
model_equation  ="y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes
                  y2 = intercept + x1 + x2 + ID + genotypes
                  y3 = intercept + x1 + ID + genotypes"
model = build_model(model_equation,
                    censored_trait=["y1"],
                    categorical_trait=["y2"])

# Step 4: Set Factors or Covariates
set_covariate(model,"x1");

# Step 5: Set Random or Fixed Effects
set_random(model,"x2");
set_random(model,"ID dam",pedigree);

# Step 6: Run Analysis
out=runMCMC(model,phenotypes,chain_length=5000);

# Step 7: Check Accuracy
results    = innerjoin(out["EBV_y1"], phenotypes, on = :ID)
accuruacy  = cor(results[!,:EBV],results[!,:bv1]) 

results    = innerjoin(out["EBV_y2"], phenotypes, on = :ID)
accuruacy  = cor(results[!,:EBV],results[!,:bv2]) 

results    = innerjoin(out["EBV_y3"], phenotypes, on = :ID)
accuruacy  = cor(results[!,:EBV],results[!,:bv3])