Multiclass Bayesian analysis

Multi-class Bayesian Alphabet models are used when multiple genomic data are included in the model. For example, genome may be split into multiple regions and markers in each region will be fitted as one class with a class-specific prior.

Below we will split genome into 3 classes. This is used to demonstrate the multi-class analysis using JWAS. Users will assign markers to multiple groups given additional biological information.

Multi-class Bayesian Alphabet model

# Step 1: Load packages
using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets


# Step 2: Read data 
phenofile   = dataset("phenotypes.csv")
pedfile     = dataset("pedigree.csv")
genofile1   = dataset("genotypes_group1.csv")
genofile2   = dataset("genotypes_group2.csv")
genofile3   = dataset("genotypes_group3.csv")


phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])
pedigree   = get_pedigree(pedfile,separator=",",header=true);

#3 groups
geno1      = get_genotypes(genofile1,separator=',',method="BayesA");
geno2      = get_genotypes(genofile2,separator=',',method="BayesB");
geno3      = get_genotypes(genofile3,separator=',',method="BayesC");


# Step 3: Build Model Equations

model_equation  ="y1 = intercept + x1 + x2 + x2*x3 + ID + dam + geno1 + geno2 + geno3
                  y2 = intercept + x1 + x2 + ID + geno1 + geno2 + geno3
                  y3 = intercept + x1 + ID + geno1 + geno2 + geno3";
model = build_model(model_equation);

# 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);

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