Public functions

build_model(model_equations::AbstractString,R=false; df::AbstractFloat=4.0, estimate_variance=true)

• Build a model from model equations with the residual variance R. In Bayesian analysis, R is the mean for the prior assigned for the residual variance with degree of freedom df, defaulting to 4.0. If R is not provided, a value is calculated from responses (phenotypes).
• By default, all variabels in modelequations are factors (categorical) and fixed. Set variables to be covariates (continuous) or random using functions `setcovariate()orset_random()`.
• The argument estimate_variance indicates whether to estimate the residual variance; estimate_variance=true is the default.

#single-trait
model_equations = "BW = intercept + age + sex"
R               = 6.72
models          = build_model(model_equations,R);

#multi-trait
model_equations = "BW = intercept + age + sex
                   CW = intercept + litter";
R               = [6.72   24.84
                   24.84  708.41]
models          = build_model(model_equations,R);

get_pedigree(pedfile::AbstractString;header=false,separator=',',missingstrings=["0"])

• Get pedigree informtion from a pedigree file with header (defaulting to false) , separator (defaulting to ,) and missing values (defaulting to ["0"])
• Pedigree file format:

a,0,0
c,a,b
d,a,c

get_genotypes(file::Union{AbstractString,Array{Float64,2},Array{Float32,2},Array{Int64,2}, Array{Int32,2}, Array{Any,2}, DataFrames.DataFrame}, G = false;
              ## method:
              method = "BayesC",Pi = 0.0,estimatePi = true, 
              ## variance:
              G_is_marker_variance = false, df = 4.0,
              estimate_variance=true, estimate_scale=false,
              constraint = false, #for multi-trait only, constraint=true means no genetic covariance among traits
              ## format:
              separator=',',header=true,
              ## quality control:
              quality_control=true, MAF = 0.01, missing_value = 9.0,
              ## others:
              center=true,starting_value=false)

• Get marker informtion from a genotype file/matrix. This file needs to be column-wise sorted by marker positions.
• If a text file is provided, the file format should be:

Animal,marker1,marker2,marker3,marker4,marker5
S1,1,0,1,1,1
D1,2,0,2,2,1
O1,1,2,0,1,0
O3,0,0,2,1,1

• If a DataFrame is provided, where n is the number of individuals and p is the number of markers,
  • This matrix needs to be column-wise sorted by marker positions.
  • The first column in the DataFrame should be individual IDs
  • The marker IDs can be provided as the header of the DataFrame. If omitted, marker IDs will be set to 1,2,3...
• If an nxp Matrix of genotypes (Array) is provided, where n is the number of individuals and p is the number of markers,
  • This matrix needs to be column-wise sorted by marker positions.
  • Individual IDs will be set to 1:n;
  • Marker IDs will be set to 1:p
• If quality_control=true, defaulting to true,
  • Missing genotypes should be denoted as 9, and will be replaced by column means. Users can also impute missing genotypes before the analysis.
  • Minor allele frequency MAF threshold, defaulting to 0.01, is uesd, and fixed loci are removed.
• G is the mean for the prior assigned for the genomic variance with degree of freedom df, defaulting to 4.0. If G is not provided, a value is calculated from responses (phenotypes).
• Available methods include "conventional (no markers)", "RR-BLUP", "BayesA", "BayesB", "BayesC", "Bayesian Lasso", and "GBLUP".
• In Bayesian variable selection methods, Pi for single-trait analyses is a number; Pi for multi-trait analyses is a dictionary such as Pi=Dict([1.0; 1.0]=>0.7,[1.0; 0.0]=>0.1,[0.0; 1.0]=>0.1,[0.0; 0.0]=>0.1), defaulting to all markers have effects (Pi = 0.0) in single-trait analysis and all markers have effects on all traits (Pi=Dict([1.0; 1.0]=>1.0,[0.0; 0.0]=>0.0)) in multi-trait analysis. Pi is estimated if estimatePi = true, , defaulting to false.
• Scale parameter for prior of marker effect variance is estimated if estimate_scale = true, defaulting to false.

set_covariate(model::MME,variables::AbstractString...)

• set variables as covariates; model is the output of function build_model().

#After running build_model, variabels age and year can be set to be covariates as
set_covariate(model,"age","year")
#or
set_covariate(model,"age year")

set_random(mme::MME,randomStr::AbstractString,ped::Pedigree, G;df=4)

• set variables as random polygenic effects with pedigree information ped. and variances G.
• G is the mean for the prior assigned for the variance with degree of freedom df, defaulting to 4.0. If G is not provided, a value is calculated from responses (phenotypes).

#single-trait (example 1)
model_equation  = "y = intercept + age + animal"
model           = build_model(model_equation,R)
ped             = get_pedigree(pedfile)
G               = 1.6
set_random(model,"animal", ped, G)

#single-trait (example 2)
model_equation  = "y = intercept + age + animal + animal*age"
model           = build_model(model_equation,R)
ped             = get_pedigree(pedfile)
G               = [1.6   0.2
                   0.2  1.0]
set_random(model,"animal animal*age", ped,G)

#multi-trait
model_equations = "BW = intercept + age + sex + animal
                   CW = intercept + age + sex + animal"
model           = build_model(model_equations,R);
ped             = get_pedigree(pedfile);
G               = [6.72   2.84
                   2.84  8.41]
set_random(model,"animal",ped,G)

set_random(mme::MME,randomStr::AbstractString,G;Vinv=0,names=[],df=4)

• set variables as random effects, defaulting to i.i.d effects, with variances G.
• G is the mean for the prior assigned for the variance with degree of freedom df, defaulting to 4.0. If G is not provided, a value is calculated from responses (phenotypes).
• the random effects are assumed to be i.i.d by default and it can be defined with any (inverse of) covariance structure Vinv with its index (row names) provided by names.

#single-trait (i.i.d randome effects)
model_equation  = "y = intercept + litter + sex"
model           = build_model(model_equation,R)
G               = 0.6
set_random(model,"litter",G)

#multi-trait (i.i.d randome effects)
model_equations = "BW = intercept + litter + sex
                   CW = intercept + litter + sex"
model           = build_model(model_equations,R);
G               = [3.72  1.84
                   1.84  3.41]
set_random(model,"litter",G)

#single-trait (randome effects with specific covariance structures)
model_equation  = "y = intercept + litter + sex"
model           = build_model(model_equation,R)
V               = [1.0  0.5 0.25
                   0.5  1.0 0.5
                   0.25 0.5 1.0]
G               = 0.6
set_random(model,"litter",G,Vinv=inv(V),names=[a1;a2;a3])

runMCMC(model::MME,df::DataFrame;
        #Data
        heterogeneous_residuals           = false,
        #MCMC
        chain_length::Integer             = 100,
        starting_value                    = false,
        burnin::Integer                   = 0,
        output_samples_frequency::Integer = chain_length>1000 ? div(chain_length,1000) : 1,
        update_priors_frequency::Integer  = 0,
        #Methods
        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                      = false,
        missing_phenotypes              = true,
        constraint                      = false,
        #Genomic Prediction
        outputEBV                       = true,
        output_heritability             = true,
        prediction_equation             = false,
        #MISC
        seed                            = false,
        printout_model_info             = true,
        printout_frequency              = chain_length+1,
        big_memory                      = false,
        double_precision                = false,
        ##MCMC samples (defaut to marker effects and hyperparametes (variance components))
        output_folder                     = "results",
        output_samples_for_all_parameters = false,
        ##for deprecated JWAS
        methods                         = "conventional (no markers)",
        Pi                              = 0.0,
        estimatePi                      = false)

Run MCMC for Bayesian Linear Mixed Models with or without estimation of variance components.

• Markov chain Monte Carlo
  • The first burnin iterations are discarded at the beginning of a MCMC chain of length chain_length.
  • Save MCMC samples every output_samples_frequency iterations, defaulting to chain_length/1000, to a folder output_folder, defaulting to results. MCMC samples for hyperparametes (variance componets) and marker effects are saved by default. MCMC samples for location parametes can be saved using function output_MCMC_samples(). Note that saving MCMC samples too frequently slows down the computation.
  • The starting_value can be provided as a vector for all location parameteres and marker effects, defaulting to 0.0s. The order of starting values for location parameters and marker effects should be the order of location parameters in the Mixed Model Equation for all traits (This can be obtained by getNames(model)) and then markers for all traits (all markers for trait 1 then all markers for trait 2...).
  • Miscellaneous Options
    • Priors are updated every update_priors_frequency iterations, defaulting to 0.
• Methods
  • Single step analysis is allowed if single_step_analysis = true and pedigree is provided.
  • Miscellaneous Options
    • Missing phenotypes are allowed in multi-trait analysis with missing_phenotypes=true, defaulting to true.
    • Catogorical Traits are allowed if categorical_trait=true, defaulting to false. Phenotypes should be coded as 1,2,3...
    • Censored traits are allowed if the upper bounds are provided in censored_trait as an array, and lower bounds are provided as phenotypes.
    • If constraint=true, defaulting to false, constrain residual covariances between traits to be zeros.
    • If causal_structure is provided, e.g., causal_structure = [0.0 0.0 0.0;1.0 0.0 0.0;1.0 0.0 0.0] for trait 1 -> trait 2 and trait 1 -> trait 3 (column index affacts row index, and a lower triangular matrix is required), phenotypic causal networks will be incorporated using structure equation models.
• Genomic Prediction
  • Predicted values for individuals of interest can be obtained based on a user-defined prediction equation prediction_equation, e.g., "y1:animal + y1:age".
  For now, genomic data is always included. Genetic values including effects defined with genotype and pedigree information are returned if prediction_equation= false, defaulting to false.
  • Individual estimted genetic values and prediction error variances (PEVs) are returned if outputEBV=true, defaulting to true. Heritability and genetic
  variances are returned if output_heritability=true, defaulting to true. Note that estimation of heritability is computaionally intensive.
• Miscellaneous Options
  • Print out the model information in REPL if printout_model_info=true; print out the monte carlo mean in REPL with printout_frequency, defaulting to false.
  • If seed, defaulting to false, is provided, a reproducible sequence of numbers will be generated for random number generation.
  • If big_memory=true, defaulting to false, a machine with lots of memory is assumed which may speed up the analysis.

outputEBV(model,IDs::Array)

Output estimated breeding values and prediction error variances for IDs.

outputMCMCsamples(mme::MME,trmStrs::AbstractString...)

Get MCMC samples for specific location parameters.

solve(mme::MME,df::DataFrame;solver="default",printout_frequency=100,tolerance = 0.000001,maxiter = 5000)

• Solve the mixed model equations (no marker information) without estimating variance components.

Available solvers include default, Jacobi, Gauss-Seidel, and Gibbs sampler.

GWAS(marker_effects_file,model;header=true,window_size=100,threshold=0.001)

run genomic window-based GWAS without a map file

• MCMC samples of marker effects are stored in markereffectsfile with delimiter ','.
• window_size is either a constant (identical number of markers in each window) or an array of number of markers in each window
• model is either the model::MME used in analysis or the genotypic covariate matrix M::Array
• File format:

GWAS(marker_effects_file;header=true)

Compute the model frequency for each marker (the probability the marker is included in the model) using samples of marker effects stored in markereffectsfile.

GWAS(model,map_file,marker_effects_file...;
     window_size = "1 Mb",sliding_window = false,
     GWAS = true, threshold = 0.001,
     genetic_correlation = false,
     header = true)

run genomic window-based GWAS

• MCMC samples of marker effects are stored in markereffectsfile with delimiter ','.
• model is either the model::MME used in analysis or the genotype cavariate matrix M::Array
• map_file has the (sorted) marker position information with delimiter ','. If the map file is not provided, i.e., map_file=false, a fake map file will be generated with window_size markers in each 1 Mb window, and each 1 Mb window will be tested.
• If two markereffectsfile are provided, and genetic_correlation = true, genomic correlation for each window is calculated.
• Statistics are computed for nonoverlapping windows of size window_size by default. If sliding_window = true, those for overlapping sliding windows are calculated.
• map file format:

markerID,chromosome,position
m1,1,16977
m2,1,434311
m3,1,1025513
m4,2,70350
m5,2,101135