Workflow

A step by step workflow for how to run JWAS is shown in this section. The workflow below is used to demonstrate a three-trait Bayesian linear mixed model fitting fixed effects (x1, x3), random effects (x2), direct genetic effects (ID), maternal genetic effects (dam) and genomic information.

• Workflow
• Available Models
• Get Data Ready
• Step 1: Load Packages
• Step 2: Read Data
• Step 3: Build Model Equations
• Step 4: Set Factors or Covariate
• Step 5: Set Random or Fixed Effects
• Step 6: Run Bayesian Analysis
• check results

Available Models

Given the data and model equations, several different types of models are available in JWAS as shown below. In the table below, "X" denotes the type of available data, and "Y<=A" denotes that Y individuals is a subset of A individuals.

Get Data Ready

By default, input data files are comma-separated values (CSV) files, where each line of the file consists of one or more fields, separated by commas. Other field separators such as space (' ') or tab ('\t') can be used if you supply the keyword argument, e.g, CSV.read(...,delim='\t') or add_genotypes(...,separator='\t')

Click on the buttons inside the tabbed menu to see the data:

Step 1: Load Packages

using JWAS,CSV,DataFrames

The JWAS package is loaded, as well as the CSV and DataFrame packages for reading text files.

Step 2: Read Data

Read Phenotypic Data

phenotypes = CSV.read("phenotypes.txt",DataFrame,delim = ',',header=true,,missingstrings=["NA"])
first(phenotypes)

output:

6×8 DataFrames.DataFrame
│ Row │ ID │ y1    │ y2   │ y3    │ x1  │ x2 │ x3 │ dam       │
├─────┼────┼───────┼──────┼───────┼─────┼────┼────┼───────────┤
│ 1   │ a1 │ -0.06 │ 3.58 │ -1.18 │ 0.9 │ 2  │ m  │ missing   │
│ 2   │ a2 │ -0.6  │ 4.9  │ 0.88  │ 0.3 │ 1  │ f  │ missing   │
│ 3   │ a3 │ -2.07 │ 3.19 │ 0.73  │ 0.7 │ 2  │ f  │ missing   │
│ 4   │ a4 │ -2.63 │ 6.97 │ -0.83 │ 0.6 │ 1  │ m  │ a2        │
│ 5   │ a5 │ 2.31  │ 3.5  │ -1.52 │ 0.4 │ 2  │ m  │ a2        │
│ 6   │ a6 │ 0.93  │ 4.87 │ -0.01 │ 5.0 │ 2  │ f  │ a3        │

The phenotypic data is read on line 1. On line 2, the first several rows of the phenotypic data are shown.

Read Pedigree Data

pedigree   = get_pedigree("pedigree.txt",separator=",",header=true)

• link to documentation for get_pedigree

The pedigree data is read on line 1.

Read Genomic Data

genotypes  = get_genotypes("genotypes.txt",G,method="BayesC",separator=",",header=true) #G is optional

• link to documentation for get_genotypes

On line 1, the genomic information is read on line 1 with the genotype file. and variance G (a 3x3 matrix). In Bayesian analysis, the 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).

Step 3: Build Model Equations

model_equation = "y1 = intercept + x1 + x3 + ID + dam + genotypes
                  y2 = intercept + x1 + x2 + x3 + ID + genotypes
                  y3 = intercept + x1 + x1*x3 + x2 + ID + genotypes"
model=build_model(model_equation,R) #R is optional

• link to documentation for build_model

The model equation for a 3-trait analysis is defined on the first 3 lines.

• The effects fitted in the model for trait y1 are the intercept, x1, x3, direct genetic effects (ID), maternal genetic effects (dam), and molecular marker effects (genotypes).
• The effects fitted in the model for trait y2 are the intercept, x1, x2, x3, direct genetic effects (ID), and molecular marker effects (genotypes).
• The effects fitted in the model for trait y3 are the intercept, x1, the interaction between x1 and x3, x2, direct genetic effects (ID), and , and molecular marker effects (genotypes).

On the last line, the model is built given the model equation and residual variance R (a 3x3 matrix). 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 effects are treated as fixed and classed as factors (categorical variables) rather than covariates (quantitative variables).

Step 4: Set Factors or Covariate

set_covariate(model,"x1")

• link to documentation for set_covariate

On line 1, the effect x1 is defined to be a covariate (a quantitative variable) rather than a factor (a categorical variable).

Step 5: Set Random or Fixed Effects

set_random(model,"x2",G1) #G1 is optional
set_random(model,"ID dam",pedigree,G2) #G2 is optional

• link to documentation for set_random

On line 1, the x2 class effect is defined as random with variance G1(a 2x2 matrix). On line 2, direct genetic effects and maternal genetic effects are fitted as ID and dam with G2 (a 4x4 matrix) and the inverse of the numerator relationship matrix defined from pedigree. In Bayesian analysis, G1 and G2 are the means for the priors assigned for the variances with degree of freedom df, defaulting to 4.0. If G1 or G2 is not provided, a value is calculated from responses (phenotypes).

Step 6: Run Bayesian Analysis

outputMCMCsamples(model,"dam")
out=runMCMC(model,phenotypes)

• link to documentation for outputMCMCsamples
• link to documentation for runMCMC

On line 1, MCMC samples from runMCMC for x2 is saved to a file, where each row represents one sample from the MCMC. On line 2, a multi-trait BayesC analysis is performed with model and phenotypes as had been defined in step 1-6. MCMC samples for marker effects, location parameters specified on line 1, and all variance components from this analysis are saved every output_samples_frequency iterations to files.

Several steps above can be skipped if no related information is available, e.g., step 4 is skipped if all effects are classed as factors. Several detailed examples are available in the examples section. Here is the link to documentation for all Public functions.

check results

Posterior means (estimate) and standard deviations (SD) of location parameters, most variance components, and marker effects are saved as the variable out and in text files. They can be listed and obtained as

keys(out)

# output:
#
# Base.KeyIterator for a Dict{Any,Any} with 7 entries. Keys:
#   "polygenic effects covariance matrix"
#   "Model frequency"
#   "residual covariance matrix"
#   "marker effects"
#   "marker effects variance"
#   "location parameters"
#   "Pi"

out["residual variance"]

# output:
#
#Covariance	Estimate	SD
#y1_y1	1.65265	0.29405
#y1_y2	-0.0290279	0.02347
#y1_y3	-0.252009	0.048289
#y2_y1	-0.0290279	0.02347
#y2_y2	0.977405	0.009732
#y2_y3	0.0451994	0.095828
#y3_y1	-0.252009	0.048289
#y3_y2	0.0451994	0.095828
#y3_y3	0.363878	0.049278

In addition, MCMC samples from posterior distributions of marker effects, all variance components, and location parameters specified in step 7, are saved to text files in your working directory. Users can compute the posterior distributions of parameters of interest using these MCMC samples files. A list of output files are shown below.

Output files:

Below is a list of files containing estimates and standard deviations for variables of interest.

Below is a list of files containing MCMC samples for variables of interest.