Updated RRM_new branch merged with up-to-date RRM master branch

Project.toml

@@ -1,26 +1,29 @@
 name = "JWAS"
 uuid = "c9a035f4-d403-5e6b-8649-6be755bc4798"
-version = "0.13.0"
 authors = ["Hao Cheng <qtlcheng@ucdavis.edu>"]
+version = "0.14.4"
 
 [deps]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+CSV = "0.5, 0.6, 0.7, 0.8"
 DataFrames = "0.19, 0.20, 0.21, 0.22"
 Distributions = "0.21, 0.22, 0.23, 0.24"
+ForwardDiff = "0.10"
 ProgressMeter = "1.0, 1.1, 1.2"
-CSV = "0.5, 0.6, 0.7, 0.8"
 StatsBase = "0.30, 0.31, 0.32, 0.33"
 julia = "1"

README.md

@@ -1,6 +1,6 @@
 # JWAS.jl
 
-[![Build Status](https://travis-ci.org/reworkhow/JWAS.jl.svg?branch=master)](https://travis-ci.org/reworkhow/JWAS.jl)
+[![Build Status](https://travis-ci.com/reworkhow/JWAS.jl.svg?branch=master)](https://travis-ci.org/reworkhow/JWAS.jl)
 [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://reworkhow.github.io/JWAS.jl/latest)
 <!---[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://reworkhow.github.io/JWAS.jl/stable)--->
 
@@ -37,3 +37,21 @@ JWAS.jl
 1. Show this README file in REPL or notebook using `?JWAS`
 2. For help on a specific function above, type ? followed by its name, e.g. `?runMCMC` and press enter.
 3. Run `Pkg.add(PackageSpec(name="JWAS", rev="master"))` to get the newest unofficial JWAS. Run `Pkg.free("JWAS")` to go back to the official one.
+
+### Examples [available here](https://github.com/reworkhow/JWAS.jl/wiki)
+
+* Single Trait Analysis
+* Multiple Trait Analysis
+* Repeated Measures
+* Single Step Analysis
+* Categorical Trait Analysis
+* Censored Trait Analysis
+* Multi-class Bayesian Analysis
+* Neural Networks
+* Cross Validation
+* Genome Wide Association Study
+* Integrating Phenotypic Causal Networks in GWAS
+* single trait and multiple trait GBLUP by providing the relationship matrix directly
+* Description of Mixed Effects Model
+* Quality Control of Genotypes
+

docs/src/log.md

@@ -0,0 +1,17 @@
+#allow intermediate omics data in neural networks (8/31/2021)
+
+1. add an argument "latent_traits" in build_model() to allow users to provide
+   which columns in the phenotypic data will be used as the observed values for
+   intermediate traits. Add latent_traits as a member in the struct mme.
+
+2. Because multi-trait models with latent traits as responses will be used in neural
+   network (between input layer and hidden layer), we reassign column names of latent
+   traits to mme.lhsVec, which will be used to make the matrices in the multi-trait models.
+   In this case, we still need a variable to save the (empirical) phenotypes (i.e., output),
+   so mme.yobs is made to save it. Add yobs as a member in the struct mme.
+
+3. mme.ySparse is used to same values for latent traits. In intermediate omics data,
+   because some/many elements in mme.ySparse are observed, only missing values in mme.ySparse
+   are sampled.
+
+# add neural network

docs/src/manual/workflow.md

@@ -283,7 +283,7 @@ Several steps above can be skipped if no related information is available, e.g.,
 
 ## check results
 
-Posterior means and standard deviations of location parameters, most variance components, and marker effects are saved as the variable `out` and in text files.
+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
 ```julia
 keys(out)
@@ -316,9 +316,42 @@ out["residual variance"]
 
 ```
 
-MCMC samples for marker effects, location parameters specified in step 7, and all variance components are saved to text
-files in your working directory. They can be obtained as
-
-```julia
-res=readdlm("MCMC_samples_marker_effects_y1.txt",',',header=true)
-```
+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. 
+
+| file name      | description |
+| -----------    | ----------- |
+| EBV_y1.txt     | estimated breeding values for trait named "y1"       |
+| EBV_y2.txt     | estimated breeding values for trait named "y2"       |
+| EBV_y3.txt     | estimated breeding values for trait named "y3"       |
+|genetic_variance.txt                   | estimated genetic variance-covariance for all traits |
+|heritability.txt                       | estimated heritability                               |
+|location_parameters.txt                | estimated non-genetic effects                        |
+|pi_genotypes.txt                       | estimated pi                                         |
+|polygenic_effects_covariance_matrix.txt| estimated variance-covariance between polygenic effects (y1_ID, y2_ID, y3_ID, y1_dam)|
+|marker_effects_genotypes.txt           | estimated marker effects for all traits              |
+|residual_variance.txt                  | estimated residual variance-covariance for all traits| 
+
+Below is a list of files containing MCMC samples for variables of interest. 
+
+| file name      | description |
+| -----------    | ----------- |
+| MCMC_samples_EBV_y1.txt     | MCMC samples from the posterior distribution of breeding values for trait named "y1"       |
+| MCMC_samples_EBV_y2.txt     | MCMC samples from the posterior distribution of breeding values for trait named "y2"       |
+| MCMC_samples_EBV_y3.txt     | MCMC samples from the posterior distribution of breeding values for trait named "y3"       |
+|MCMC_samples_genetic_variance.txt                   | MCMC samples from the posterior distribution of genetic variance-covariance for all traits   |
+|MCMC_samples_heritability.txt                       | MCMC samples from the posterior distribution of heritability                                 | 
+|MCMC_samples_marker_effects_genotypes_y1            |MCMC samples from the posterior distribution of marker effects for trait named "y1"           |
+|MCMC_samples_marker_effects_genotypes_y2            |MCMC samples from the posterior distribution of marker effects for trait named "y2"           |
+|MCMC_samples_marker_effects_genotypes_y3            |MCMC samples from the posterior distribution of marker effects for trait named "y3"           |
+|MCMC_samples_marker_effects_variances_genotypes.txt | MCMC samples from the posterior distribution of marker effect variance for all traits        |
+|MCMC_samples_pi_genotypes.txt                       | MCMC samples from the posterior distribution of pi                                           |
+|MCMC_samples_polygenic_effects_variance.txt         |  MCMC samples from the posterior distribution of variance-covariance between y1_ID, y2_ID, y3_ID, y1_dam|
+|MCMC_samples_residual_variance.txt                  |  MCMC samples from the posterior distribution of residual variance-covariance for all traits |
+|MCMC_samples_y1.dam.txt                             | MCMC samples from the posterior distribution of dam effect for y1|
+|MCMC_samples_y1.ID_y2.ID_y3.ID_y1.dam_variances.txt | MCMC samples from the posterior distribution of variance-covariance between y1_ID, y2_ID, y3_ID, y1_dam|
+|MCMC_samples_y2.x2_y3.x2_variances.txt              | MCMC samples from the posterior distribution of variance-covariance between y2_x2 and y3_x2| 

src/1.JWAS/src/JWAS.jl

@@ -6,6 +6,7 @@ using SparseArrays
 using LinearAlgebra
 using ProgressMeter
 using .PedModule
+using ForwardDiff
 
 import StatsBase: describe #a new describe is exported
 
@@ -48,6 +49,7 @@ include("structure_equation_model/SEM.jl")
 #Latent Traits
 include("Nonlinear/nonlinear.jl")
 include("Nonlinear/bnn_hmc.jl")
+include("Nonlinear/nnbayes_check.jl")
 
 #input
 include("input_data_validation.jl")
@@ -139,10 +141,10 @@ export dataset
         * 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 2 -> trait 1 and trait 3 -> trait 1, phenotypic causal networks will be incorporated using structure equation models.
+          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 an user-defined prediction equation `prediction_equation`, e.g., "y1:animal + y1:geno + y1:age".
-    For now, genomic data is always included. Genetic values including effects defined with genotype and pedigre information are returned if `prediction_equation`= false, defaulting to `false`.
+    * 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
@@ -168,10 +170,10 @@ function runMCMC(mme::MME,df;
                 categorical_trait               = false,
                 censored_trait                  = false,
                 causal_structure                = false,
-                mega_trait                      = false,
-                RRM                             = false,
+                mega_trait                      = mme.nonlinear_function == false ? false : true, #NNBayes -> default mega_trait=true
                 missing_phenotypes              = true,
                 constraint                      = false,
+                RRM                             = false,
                 #Genomic Prediction
                 outputEBV                       = true,
                 output_heritability             = true,
@@ -191,9 +193,11 @@ function runMCMC(mme::MME,df;
                 estimatePi                      = false,
                 estimateScale                   = false)
 
-    #Nonlinear
-    if mme.latent_traits == true
-        yobs = df[!,Symbol(string(Symbol(mme.lhsVec[1]))[1:(end-1)])]
+
+    #Neural Network
+    is_nnbayes_partial = (mme.nonlinear_function != false && mme.is_fully_connected==false)
+    if mme.nonlinear_function != false #modify data to add phenotypes for hidden nodes
+        yobs = df[!,Symbol(string(Symbol(mme.lhsVec[1]))[1:(end-1)])]#a number label is added to original trait name in nnbayes_model_equation()
         for i in mme.lhsVec
             df[!,i]= yobs
         end
@@ -216,6 +220,13 @@ function runMCMC(mme::MME,df;
     if seed != false
         Random.seed!(seed)
     end
+    #when using multi-thread, make sure the results are reproducible for users
+    nThread = Threads.nthreads()
+    if nThread>1
+        Threads.@threads for i = 1:nThread
+              Random.seed!(seed+i)
+        end
+    end
     ############################################################################
     # Create an folder to save outputs
     ############################################################################
@@ -270,14 +281,11 @@ function runMCMC(mme::MME,df;
     mme.causal_structure = causal_structure
     if causal_structure != false
         #no missing phenotypes and residual covariance for identifiability
-        missing_phenotypes, constraint = false, true
+        mme.MCMCinfo.missing_phenotypes, mme.MCMCinfo.constraint = false, true
         if !istril(causal_structure)
             error("The causal structue needs to be a lower triangular matrix.")
         end
     end
-
-
-
     # Double Precision
     if double_precision == true
         if mme.M != 0
@@ -295,21 +303,19 @@ function runMCMC(mme::MME,df;
         end
         mme.Gi = map(Float64,mme.Gi)
     end
-    #mega_trait
-    if mme.MCMCinfo.mega_trait == true || mme.MCMCinfo.constraint == true
-        if mme.nModels == 1
-            error("more than 1 trait is required for MegaLMM analysis.")
-        end
-        mme.MCMCinfo.constraint = true
-        ##sample from scale-inv-⁠χ2, not InverseWishart
-        mme.df.residual  = mme.df.residual - mme.nModels
-        mme.scaleR       = diag(mme.scaleR/(mme.df.residual - 1))*(mme.df.residual-2)/mme.df.residual #diag(R_prior_mean)*(ν-2)/ν
-        if mme.M != 0
-            for Mi in mme.M
-                Mi.df        = Mi.df - mme.nModels
-                Mi.scale    = diag(Mi.scale/(Mi.df - 1))*(Mi.df-2)/Mi.df
-            end
-        end
+
+    # NNBayes mega trait: from multi-trait to multiple single-trait
+    if mme.MCMCinfo.mega_trait == true
+        printstyled(" - Bayesian Alphabet:                multiple independent single-trait Bayesian models are used to sample marker effect. \n",bold=false,color=:green)
+        printstyled(" - Multi-threaded parallelism:       $nThread threads are used to run single-trait models in parallel. \n",bold=false,color=:green)
+        nnbayes_mega_trait(mme)
+    elseif mme.nonlinear_function != false  #only print for NNBayes
+        printstyled(" - Bayesian Alphabet:                multi-trait Bayesian models are used to sample marker effect. \n",bold=false,color=:green)
+    end
+
+    # NNBayes: modify parameters for partial connected NN
+    if is_nnbayes_partial
+        nnbayes_partial_para_modify2(mme)
     end
     ############################################################################
     #Make incidence matrices and genotype covariates for training observations
@@ -328,18 +334,15 @@ function runMCMC(mme::MME,df;
     ############################################################################
     describe(mme)
     if mme.nModels ==1 && RRM != false
-        mme.output=MCMC_BayesianAlphabet_RRM(chain_length,mme,df,
+        mme.output=MCMC_BayesianAlphabet_RRM(mme,df,
                                   Φ                        = RRM,
                                   Pi                       = Pi,
                                   burnin                   = burnin,
-                                  methods                  = methods,
-                                  estimatePi               = estimatePi,
-                                  estimateScale            = estimateScale,
-                                  starting_value           = mme.MCMCinfo.starting_value,
+                                  starting_value           = mme.sol,
                                   outFreq                  = printout_frequency,
                                   output_samples_frequency = output_samples_frequency,
-                                  output_file              = output_samples_file,
-                                  update_priors_frequency  = update_priors_frequency)
+                                  nIter                    = chain_length,
+                                  output_folder            = mme.MCMCinfo.output_folder)
     else
         mme.output=MCMC_BayesianAlphabet(mme,df)
     end
@@ -362,7 +365,22 @@ function runMCMC(mme::MME,df;
     versioninfo()
     printstyled("\n\nThe analysis has finished. Results are saved in the returned ",bold=true)
     printstyled("variable and text files. MCMC samples are saved in text files.\n\n\n",bold=true)
+
+    # make MCMC samples for indirect marker effect
+    if causal_structure != false
+
+        #generate the MCMC sample file for indirect and direct effect file.
+        generate_indirect_marker_effect_sample(mme.lhsVec,output_folder,causal_structure,"structure_coefficient_MCMC_samples.txt")
+        generate_overall_marker_effect_sample(mme.lhsVec,output_folder,causal_structure)
+
+        # generate marker effet file for direct, indirect, and overall effect
+        generate_marker_effect(mme.lhsVec, output_folder,causal_structure,"direct")
+        generate_marker_effect(mme.lhsVec, output_folder,causal_structure,"indirect")
+        generate_marker_effect(mme.lhsVec, output_folder,causal_structure,"overall")
+    end
+
     return mme.output
+
 end
 ################################################################################
 # Print out Model or MCMC information
@@ -429,6 +447,7 @@ end
 * (internal function) Print out MCMC information.
 """
 function getMCMCinfo(mme)
+    is_nnbayes_partial       = mme.nonlinear_function != false && mme.is_fully_connected==false
     if mme.MCMCinfo == false
         printstyled("MCMC information is not available\n\n",bold=true)
         return
@@ -479,18 +498,26 @@ function getMCMCinfo(mme)
             @printf("%-30s %20s\n","Method",Mi.method)
             for Mi in mme.M
                 if Mi.genetic_variance != false
-                    if mme.nModels == 1
+                    if mme.nModels == 1 && mme.MCMCinfo.RRM == false || is_nnbayes_partial
                         @printf("%-30s %20.3f\n","genetic variances (genomic):",Mi.genetic_variance)
-                    else
+                    elseif mme.nModels==1 && mme.MCMCinfo.RRM != false
+                        @printf("%-30s\n","genetic variances (genomic):")
+                        Base.print_matrix(stdout,round.(Mi.genetic_variance,digits=3))
+                        println()
+                    elseif mme.nModels>1
                         @printf("%-30s\n","genetic variances (genomic):")
                         Base.print_matrix(stdout,round.(Mi.genetic_variance,digits=3))
                         println()
                     end
                 end
                 if !(Mi.method in ["GBLUP"])
-                    if mme.nModels == 1
+                    if mme.nModels == 1 && mme.MCMCinfo.RRM == false || is_nnbayes_partial
                         @printf("%-30s %20.3f\n","marker effect variances:",Mi.G)
-                    else
+                    elseif mme.nModels==1 && mme.MCMCinfo.RRM != false
+                        @printf("%-30s\n","marker effect variances:")
+                        Base.print_matrix(stdout,round.(Mi.G,digits=3))
+                        println()
+                    elseif mme.nModels>1
                         @printf("%-30s\n","marker effect variances:")
                         Base.print_matrix(stdout,round.(Mi.G,digits=3))
                         println()