Introduction to JWAS
Contents
3. Introduction to JWAS#
JWAS is a package for Bayesian mixed effects model analyses
Can accommodate:
breeding value models (animal models)
maternal effects
multiple traits
SNP effects
Non-SNP random effects are assumed to be normallly distributed
SNP effects can be normal (RR-BLUP), t-distributed (BayesA), mixture (BayesB and BayesC), or Laplace (Bayesian LASSO)
Inference based on MCMC samples
3.1. Non-SNP part of the model: examples#
3.1.1. Single-trait models#
3.1.1.1. Data#
using CSV, DataFrames, Statistics
data1 = CSV.read("singleTraitEx.phen",DataFrame)
| Line | Age | Height | |
|---|---|---|---|
| Int64 | Int64 | Float64 | |
| 1 | 1 | 20 | 20.0 |
| 2 | 1 | 19 | 21.1 |
| 3 | 2 | 21 | 15.9 |
| 4 | 2 | 20 | 13.7 |
| 5 | 3 | 18 | 18.4 |
| 6 | 3 | 20 | 22.0 |
3.1.1.2. One-way model#
using JWAS, JWAS.Datasets
varRes = 10.0
model1 = build_model("Height = intercept + Line", varRes);
runMCMC(model1,data1,chain_length=500,output_samples_frequency=100);
running MCMC ...100%|███████████████████████████████████| Time: 0:00:00
The folder results is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
Height = intercept + Line
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Line factor fixed 3
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 100
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances: 10.000
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 4.000
The file results/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
model1.X
6×4 SparseArrays.SparseMatrixCSC{Float32, Int64} with 12 stored entries:
1.0 1.0 ⋅ ⋅
1.0 1.0 ⋅ ⋅
1.0 ⋅ 1.0 ⋅
1.0 ⋅ 1.0 ⋅
1.0 ⋅ ⋅ 1.0
1.0 ⋅ ⋅ 1.0
Matrix(model1.X)
6×4 Matrix{Float32}:
1.0 1.0 0.0 0.0
1.0 1.0 0.0 0.0
1.0 0.0 1.0 0.0
1.0 0.0 1.0 0.0
1.0 0.0 0.0 1.0
1.0 0.0 0.0 1.0
JWAS.getNames(model1)
4-element Vector{AbstractString}:
"Height:intercept:intercept"
"Height:Line:1"
"Height:Line:2"
"Height:Line:3"
3.1.1.3. One-way model without intercept#
varRes = 10.0
model1 = build_model("Height = Line", varRes);
runMCMC(model1,data1,chain_length=500,output_samples_frequency=10);
The folder results already exists.
The folder results1 is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
Height = Line
Model Information:
Term C/F F/R nLevels
Line factor fixed 3
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances: 10.000
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 4.000
The file results1/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
JWAS.getNames(model1)
3-element Vector{AbstractString}:
"Height:Line:1"
"Height:Line:2"
"Height:Line:3"
Matrix(model1.X)
6×3 Matrix{Float32}:
1.0 0.0 0.0
1.0 0.0 0.0
0.0 1.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
0.0 0.0 1.0
3.1.1.4. Model with covariate#
varRes = 10.0
model1 = build_model("Height = intercept + Line + Age", varRes)
set_covariate(model1,"Age");
runMCMC(model1,data1,chain_length=500,output_samples_frequency=10);
The folder results already exists.
The folder results1 already exists.
The folder results2 is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
Height = intercept + Line + Age
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Line factor fixed 3
Age covariate fixed 1
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances: 10.000
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 4.000
The file results2/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
3.1.1.5. This model has a common slope#
JWAS.getNames(model1)
5-element Vector{AbstractString}:
"Height:intercept:intercept"
"Height:Line:1"
"Height:Line:2"
"Height:Line:3"
"Height:Age:Age"
Matrix(model1.X)
6×5 Matrix{Float32}:
1.0 1.0 0.0 0.0 20.0
1.0 1.0 0.0 0.0 19.0
1.0 0.0 1.0 0.0 21.0
1.0 0.0 1.0 0.0 20.0
1.0 0.0 0.0 1.0 18.0
1.0 0.0 0.0 1.0 20.0
3.1.1.6. Model with Age by Line#
varRes = 10.0
model1 = build_model("Height = intercept + Line + Line*Age", varRes)
set_covariate(model1,"Age");
runMCMC(model1,data1,chain_length=500,output_samples_frequency=10);
The folder results already exists.
The folder results1 already exists.
The folder results2 already exists.
The folder results3 is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
Height = intercept + Line + Line*Age
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Line factor fixed 3
Line*Age interaction fixed 3
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances: 10.000
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 4.000
The file results3/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
3.1.1.7. This model has line-specific slopes#
JWAS.getNames(model1)
7-element Vector{AbstractString}:
"Height:intercept:intercept"
"Height:Line:1"
"Height:Line:2"
"Height:Line:3"
"Height:Line*Age:1 * Age"
"Height:Line*Age:2 * Age"
"Height:Line*Age:3 * Age"
Matrix(model1.X)
6×7 Matrix{Float32}:
1.0 1.0 0.0 0.0 20.0 0.0 0.0
1.0 1.0 0.0 0.0 19.0 0.0 0.0
1.0 0.0 1.0 0.0 0.0 21.0 0.0
1.0 0.0 1.0 0.0 0.0 20.0 0.0
1.0 0.0 0.0 1.0 0.0 0.0 18.0
1.0 0.0 0.0 1.0 0.0 0.0 20.0
3.1.1.8. Model with Age by Line and common intercept#
varRes = 10.0
model1 = build_model("Height = intercept + Line*Age", varRes)
set_covariate(model1,"Age");
runMCMC(model1,data1,chain_length=500,output_samples_frequency=10);
The folder results already exists.
The folder results1 already exists.
The folder results2 already exists.
The folder results3 already exists.
The folder results4 is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
Height = intercept + Line*Age
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Line*Age interaction fixed 3
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances: 10.000
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 4.000
The file results4/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
3.1.1.9. This model has common intercept and line-specific slopes#
JWAS.getNames(model1)
4-element Vector{AbstractString}:
"Height:intercept:intercept"
"Height:Line*Age:1 * Age"
"Height:Line*Age:2 * Age"
"Height:Line*Age:3 * Age"
Matrix(model1.X)
6×4 Matrix{Float32}:
1.0 20.0 0.0 0.0
1.0 19.0 0.0 0.0
1.0 0.0 21.0 0.0
1.0 0.0 20.0 0.0
1.0 0.0 0.0 18.0
1.0 0.0 0.0 20.0
3.1.2. Two-trait models#
3.1.2.1. Data#
data2 = CSV.read("twoTraitEx.phen",DataFrame)
| Line | Age | y1 | y2 | |
|---|---|---|---|---|
| Int64 | Int64 | Float64 | Float64 | |
| 1 | 1 | 20 | 20.0 | 6.2 |
| 2 | 1 | 19 | 21.1 | 5.9 |
| 3 | 2 | 21 | 15.9 | 10.0 |
| 4 | 2 | 20 | 13.7 | 8.2 |
| 5 | 3 | 18 | 18.4 | 9.6 |
| 6 | 3 | 20 | 22.0 | 11.0 |
cov(Matrix(data2[!,3:4]))
2×2 Matrix{Float64}:
10.2137 -0.805667
-0.805667 4.36967
3.1.2.2. Two-trait model Age by Line for trait 1#
modelEquation = "y1 = intercept + Line + Line*Age
y2 = intercept + Line"
varRes = [10.0 -1.0
-1.0 5.0]
model2 = build_model(modelEquation,varRes)
set_covariate(model2,"Age");
out = runMCMC(model2,data2,chain_length=500,output_samples_frequency=10);
running MCMC ...100%|███████████████████████████████████| Time: 0:00:00
The folder results already exists.
The folder results1 already exists.
The folder results2 already exists.
The folder results3 already exists.
The folder results4 already exists.
The folder results5 is created to save results.
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).Default or user-defined prediction equation are not available.
Phenotypes for 6 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
y1 = intercept + Line + Line*Age
y2 = intercept + Line
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Line factor fixed 3
Line*Age interaction fixed 3
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
residual variances:
10.0f0 -1.0f0
-1.0f0 5.0f0
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 6.000
The file results5/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
JWAS.getNames(model2)
11-element Vector{AbstractString}:
"y1:intercept:intercept"
"y1:Line:1"
"y1:Line:2"
"y1:Line:3"
"y1:Line*Age:1 * Age"
"y1:Line*Age:2 * Age"
"y1:Line*Age:3 * Age"
"y2:intercept:intercept"
"y2:Line:1"
"y2:Line:2"
"y2:Line:3"
Matrix(model2.X)
12×11 Matrix{Float32}:
1.0 1.0 0.0 0.0 20.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 0.0 0.0 19.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 0.0 1.0 0.0 0.0 21.0 0.0 0.0 0.0 0.0 0.0
1.0 0.0 1.0 0.0 0.0 20.0 0.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 1.0 0.0 0.0 18.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 1.0 0.0 0.0 20.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0
3.1.3. Breeding value models#
data3 = CSV.read("twoTraitMaternal.phen",DataFrame,missingstrings=["0"])
| Ind | Mat | y1 | y2 | |
|---|---|---|---|---|
| Int64 | Int64? | Float64 | Float64 | |
| 1 | 1 | missing | 10.0 | 11.0 |
| 2 | 2 | missing | 9.7 | 12.0 |
| 3 | 3 | 2 | 8.9 | 11.9 |
| 4 | 4 | 2 | 9.7 | 10.8 |
| 5 | 5 | 4 | 8.8 | 11.9 |
pedigree = get_pedigree("pedFile",separator=",",header=false);
The delimiter in pedFile is ','.
Pedigree information:
#individuals: 5
#sires: 2
#dams: 2
#founders: 2
varGen = [1.0 0.0 0.0
0.0 0.1 0.0
0.0 0.0 0.2];
varRes = [2.8 -0.74
-0.74 0.287]
modelEq3 = "y1 = intercept + Ind + Mat
y2 = intercept + Ind"
model3 = build_model(modelEq3,varRes)
set_random(model3,"Ind Mat",pedigree,varGen)
Mat is not found in model equation 2.
Can see below how “Ind” and “Mat” are ordered
model3.pedTrmVec
3-element Vector{Any}:
"y1:Ind"
"y2:Ind"
"y1:Mat"
out3 = runMCMC(model3,data3,chain_length=500,output_samples_frequency=10);
running MCMC ...100%|███████████████████████████████████| Time: 0:00:00
The folder results already exists.
The folder results1 already exists.
The folder results2 already exists.
The folder results3 already exists.
The folder results4 already exists.
The folder results5 already exists.
The folder results6 is created to save results.
Checking pedigree...
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any["y1:Ind", "y2:Ind", "y1:Mat"] (Note that genomic data is always included for now).Missing values are found in column Mat for some observations.Effects of this variable on such observations are considered as zeros. It will be used in the estimation of this effect. Users may impute missing values before the analysis
Missing values are found in column Mat for some observations.Effects of this variable on such observations are considered as zeros. It will be used in prediction. Users may impute missing values before the analysis.
Phenotypes for 5 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
A Linear Mixed Model was build using model equations:
y1 = intercept + Ind + Mat
y2 = intercept + Ind
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
Ind factor random 5
Mat factor random 5
MCMC Information:
chain_length 500
burnin 0
starting_value true
printout_frequency 501
output_samples_frequency 10
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
random effect variances (y1:Ind,y2:Ind,y1:Mat):
1.0f0 0.0f0 0.0f0
0.0f0 0.1f0 0.0f0
0.0f0 0.0f0 0.2f0
genetic variances (polygenic):
1.0f0 0.0f0 0.0f0
0.0f0 0.1f0 0.0f0
0.0f0 0.0f0 0.2f0
residual variances:
2.8f0 -0.74f0
-0.74f0 0.287f0
Genomic Information:
Degree of freedom for hyper-parameters:
residual variances: 6.000
polygenic effect variances: 7.000
The file results6/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The file results6/MCMC_samples_polygenic_effects_variance.txt is created to save MCMC samples for polygenic_effects_variance.
The file results6/MCMC_samples_y1.Ind_y2.Ind_y1.Mat_variances.txt is created to save MCMC samples for y1:Ind_y2:Ind_y1:Mat_variances.
The file results6/MCMC_samples_EBV_y1.txt is created to save MCMC samples for EBV_y1.
The file results6/MCMC_samples_EBV_y2.txt is created to save MCMC samples for EBV_y2.
The file results6/MCMC_samples_genetic_variance.txt is created to save MCMC samples for genetic_variance.
The file results6/MCMC_samples_heritability.txt is created to save MCMC samples for heritability.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
JWAS.getNames(model3)
17-element Vector{AbstractString}:
"y1:intercept:intercept"
"y1:Ind:1"
"y1:Ind:2"
"y1:Ind:4"
"y1:Ind:3"
"y1:Ind:5"
"y1:Mat:1"
"y1:Mat:2"
"y1:Mat:4"
"y1:Mat:3"
"y1:Mat:5"
"y2:intercept:intercept"
"y2:Ind:1"
"y2:Ind:2"
"y2:Ind:4"
"y2:Ind:3"
"y2:Ind:5"
3.2. Adding SNP effects to the model#
# 1. Load packages
using JWAS,DataFrames,CSV,Statistics,HTTP;
# 2. Read data
function getdata(file_name) # function to load data from github folder
http_obj = HTTP.get("https://raw.githubusercontent.com/zhaotianjing/bio_protocol/main/data/$file_name.txt")
data = CSV.read(http_obj.body, DataFrame,header=true,missingstrings=["."])
end
phenotypes = getdata("phenotypes")
genotypes = getdata("genotypes_5k")
pedigree = getdata("pedigree");
pedigree = get_pedigree(pedigree,separator=",",header=true)
genotypes = get_genotypes(genotypes,separator=',',method="BayesC")
# 3. Build Model Equations
model_equation ="t1 = intercept + genotypes
t2 = intercept + ID + genotypes
t3 = intercept + ID + genotypes";
model = build_model(model_equation);
# 4. Set Factors or Covariates (no covariates are fitted in the model specified in step 3
# example:
# set_covariate(model,"age");
# 5. Set Random or Fixed Effects
set_random(model,"ID",pedigree);
# other examples:
# set_random(model,"ID dam",pedigree);
# 6. Run Analysis
out=runMCMC(model,phenotypes,chain_length=1000);
# other examples:
# out=runMCMC(model, phenotypes, chain_length=50000, burnin=10000, output_samples_frequency= 100);
# out=runMCMC(model, phenotypes, chain_length=50000, single_step_analysis=true, pedigree=pedigree); for single-step Bayesian regression analysis
# 7. Check Results
# The estimates (i.e., posterior mean), standard errors, and MCMC samples for parameters of interest
# will be saved in the output folder after the analysis.
# example:
# -------------------------------------------------------------------------------------------------
# EBV_t1.txt | estimated breeding values for trait "t1"
# EBV_t2.txt | estimated breeding values for trait "t2"
# EBV_t3.txt | estimated breeding values for trait "t3"
# genetic_variance.txt | estimated genetic variance-covariance matrix across traits
# residual_variance.txt | estimated residual variance-covariance matrix across traits
# heritability.txt | estimated heritabilities for analyzed traits
# pi_genotypes.txt | estimated probability for marker inclusion patterns
# marker_effects_genotypes.txt | estimated marker effects for analyzed traits
# polygenic_effects_covariance_matrix.txt | estimated variance-covariance matrix for polygenic effects
# --------------------------------------------------------------------------------------------------
# Genomic prediction or GWAS analysis
# 1. Calculate accuracy for genomic prediction (use t3 as an example)
results = innerjoin(out["EBV_t3"], phenotypes, on = :ID)
ind_id = findall(x -> !ismissing(x), results[!,:t3])
accuruacy = cor(results[ind_id,:EBV],results[ind_id,:t3]) # pearson correlation between estimated breeding values and observed phentypes
# 2. GWAS
# (a) Compute the model frequency for each marker
#marker_effects_file = "results/MCMC_samples_marker_effects_genotypes_t1.txt" # MCMC samples of marker effects for trait t1
#GWAS(marker_effects_file,header=true)
# (b) Compute the posterior probability of association for genomic windows based on their explained proportion of total genetic variance
# example1:
# map_file = "map.csv"
# marker_effects_file = "results/MCMC_samples_marker_effects_genotypes_t1.txt"
# out=GWAS(model,map_file, marker_effects_file, header=true, window_size="1 Mb",threshold=0.01)
# example2:
# map_file="map.csv"
# marker_effects_file1="results/MCMC_samples_marker_effects_genotypes_t1.txt"
# marker_effects_file2="results/MCMC_samples_marker_effects_genotypes_t2.txt"
# marker_effects_file3="results/MCMC_samples_marker_effects_genotypes_t3.txt"
# out=GWAS(model,map_file,marker_effects_file1,marker_effects_file2,marker_effects_file3,header=true,window_size="1 Mb")
# (c) Compute the genetic correlation for each genomic window between two analyzed traits
# example:
# map_file="map.csv"
# marker_effects_file1 ="results/MCMC_samples_marker_effects_genotypes_y1.txt"
# marker_effects_file2 ="results/MCMC_samples_marker_effects_genotypes_y2.txt"
# out=GWAS(model,map_file,marker_effects_file1,marker_effects_file2,genetic_correlation=true,header=true,window_size="1 Mb")
running MCMC ...100%|███████████████████████████████████| Time: 0:02:02
Pedigree information:
#individuals: 6474
#sires: 1012
#dams: 3103
#founders: 1
The first column in the dataframe should be individual IDs.
The data type of markers should be Number.
Missing values (9.0) are replaced by column means.
217 loci which are fixed or have minor allele frequency < 0.01 are removed.
Genotype informatin:
#markers: 4783; #individuals: 3534
ID is not found in model equation 1.
The folder results already exists.
The folder results1 already exists.
The folder results2 already exists.
The folder results3 already exists.
The folder results4 already exists.
The folder results5 already exists.
The folder results6 already exists.
The folder results7 is created to save results.
Checking pedigree...
Checking genotypes...
Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Predicted values for individuals of interest will be obtained as the summation of Any["t2:ID", "t3:ID"] (Note that genomic data is always included for now).Phenotypes for 3459 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.
Prior information for genomic variance is not provided and is generated from the data.
Prior information for residual variance is not provided and is generated from the data.
Prior information for random effect variance is not provided and is generated from the data.
Pi (Π) is not provided.
Pi (Π) is generated assuming all markers have effects on all traits.
The prior for marker effects covariance matrix is calculated from genetic covariance matrix and Π.
The mean of the prior for the marker effects covariance matrix is:
0.000225 0.0 0.0
0.0 0.000195 0.0
0.0 0.0 0.000143
A Linear Mixed Model was build using model equations:
t1 = intercept + genotypes
t2 = intercept + ID + genotypes
t3 = intercept + ID + genotypes
Model Information:
Term C/F F/R nLevels
intercept factor fixed 1
ID factor random 6474
MCMC Information:
chain_length 1000
burnin 0
starting_value true
printout_frequency 1001
output_samples_frequency 1
constraint false
missing_phenotypes true
update_priors_frequency 0
seed false
Hyper-parameters Information:
random effect variances (t2:ID,t3:ID):
0.315 0.0
0.0 0.231
genetic variances (polygenic):
0.315 0.0
0.0 0.231
residual variances:
0.729 0.0 0.0
0.0 0.631 0.0
0.0 0.0 0.462
Genomic Information:
complete genomic data (i.e., non-single-step analysis)
Genomic Category genotypes
Method BayesC
genetic variances (genomic):
0.365 0.0 0.0
0.0 0.315 0.0
0.0 0.0 0.231
marker effect variances:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
Π: (Y(yes):included; N(no):excluded)
["t1", "t2", "t3"] probability
["N", "Y", "N"] 0.0
["Y", "N", "Y"] 0.0
["N", "N", "Y"] 0.0
["Y", "N", "N"] 0.0
["N", "N", "N"] 0.0
["Y", "Y", "Y"] 1.0
["N", "Y", "Y"] 0.0
["Y", "Y", "N"] 0.0
estimatePi true
estimateScale false
Degree of freedom for hyper-parameters:
residual variances: 7.000
polygenic effect variances: 6.000
marker effect variances: 4.000
The file results7/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.
The file results7/MCMC_samples_polygenic_effects_variance.txt is created to save MCMC samples for polygenic_effects_variance.
The file results7/MCMC_samples_marker_effects_genotypes_t1.txt is created to save MCMC samples for marker_effects_genotypes_t1.
The file results7/MCMC_samples_marker_effects_genotypes_t2.txt is created to save MCMC samples for marker_effects_genotypes_t2.
The file results7/MCMC_samples_marker_effects_genotypes_t3.txt is created to save MCMC samples for marker_effects_genotypes_t3.
The file results7/MCMC_samples_marker_effects_variances_genotypes.txt is created to save MCMC samples for marker_effects_variances_genotypes.
The file results7/MCMC_samples_pi_genotypes.txt is created to save MCMC samples for pi_genotypes.
The file results7/MCMC_samples_t2.ID_t3.ID_variances.txt is created to save MCMC samples for t2:ID_t3:ID_variances.
The file results7/MCMC_samples_EBV_t1.txt is created to save MCMC samples for EBV_t1.
The file results7/MCMC_samples_EBV_t2.txt is created to save MCMC samples for EBV_t2.
The file results7/MCMC_samples_EBV_t3.txt is created to save MCMC samples for EBV_t3.
The file results7/MCMC_samples_genetic_variance.txt is created to save MCMC samples for genetic_variance.
The file results7/MCMC_samples_heritability.txt is created to save MCMC samples for heritability.
The version of Julia and Platform in use:
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.2.0)
CPU: Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, cyclone)
The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.
0.8047499262574457