{ "cells": [ { "cell_type": "code", "execution_count": 66, "id": "7fadb1ac-3e38-4019-ac26-4af39ec227b9", "metadata": {}, "outputs": [], "source": [ "# Install Julia from https://julialang.org/\n", "\n", "# Install JWAS from https://reworkhow.github.io/JWAS.jl/latest/\n" ] }, { "cell_type": "code", "execution_count": 67, "id": "137ee049-269a-4184-9dcd-0046ac50d71e", "metadata": {}, "outputs": [], "source": [ "# Single-trait analysis - https://github.com/reworkhow/JWAS.jl/wiki/Single-Trait-Analysis" ] }, { "cell_type": "code", "execution_count": 68, "id": "cb8647c3-fe3b-4687-9cb9-4f851cf5b21e", "metadata": {}, "outputs": [], "source": [ "# Step 1: Load packages" ] }, { "cell_type": "code", "execution_count": 69, "id": "2eb0830b-ee75-4912-b275-54b0090fc9b1", "metadata": {}, "outputs": [], "source": [ "using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets" ] }, { "cell_type": "code", "execution_count": 70, "id": "5bb92302-5e91-4bf9-b145-b098ac95a55d", "metadata": {}, "outputs": [], "source": [ "# Step 2: Read data " ] }, { "cell_type": "code", "execution_count": 71, "id": "7a13b0a9-0f7f-4439-b589-9708ac9a4189", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"/Users/gota/.julia/packages/JWAS/6mjQ8/src/4.Datasets/src/../data/genotypes.csv\"" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phenofile = dataset(\"phenotypes.csv\")\n", "pedfile = dataset(\"pedigree.csv\")\n", "genofile = dataset(\"genotypes.csv\")" ] }, { "cell_type": "code", "execution_count": 72, "id": "8bcfa745-898c-4d70-880d-975776ad1999", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32mThe delimiter in pedigree.csv is ','.\u001b[39m\n", "Pedigree information:\n", "#individuals: 100\n", "#sires: 39\n", "#dams: 38\n", "#founders: 20\n", "\u001b[32mThe delimiterd in genotypes.csv is ','. \u001b[39m\u001b[32mThe header (marker IDs) is provided in genotypes.csv.\u001b[39m\n", "Genotype informatin:\n", "#markers: 1000; #individuals: 60\n" ] } ], "source": [ "phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=[\"NA\"])\n", "pedigree = get_pedigree(pedfile,separator=\",\",header=true);\n", "genotypes = get_genotypes(genofile,separator=',',method=\"BayesC\");" ] }, { "cell_type": "code", "execution_count": 73, "id": "094c60f8-c92e-4cc3-a9a5-46bab492aadf", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

5 rows × 11 columns (omitted printing of 2 columns)

IDy1y2y3x1x2x3dambv1
StringFloat64Float64Float64Float64StringStringString?Float64
1a1-0.579682-0.05820891.04560.77g1mmissing-1.07018
2a2-2.02246-2.3843-1.57196-1.02g2mmissing-1.03328
3a3-1.48071-0.421258-0.2722450.52g1fmissing-1.2916
4a4-3.03031-2.0634-0.604634-1.05g4mmissing-1.39308
5a52.188192.224251.723062.06g3mmissing0.267546
" ], "text/latex": [ "\\begin{tabular}{r|cccccccccc}\n", "\t& ID & y1 & y2 & y3 & x1 & x2 & x3 & dam & bv1 & \\\\\n", "\t\\hline\n", "\t& String & Float64 & Float64 & Float64 & Float64 & String & String & String? & Float64 & \\\\\n", "\t\\hline\n", "\t1 & a1 & -0.579682 & -0.0582089 & 1.0456 & 0.77 & g1 & m & \\emph{missing} & -1.07018 & $\\dots$ \\\\\n", "\t2 & a2 & -2.02246 & -2.3843 & -1.57196 & -1.02 & g2 & m & \\emph{missing} & -1.03328 & $\\dots$ \\\\\n", "\t3 & a3 & -1.48071 & -0.421258 & -0.272245 & 0.52 & g1 & f & \\emph{missing} & -1.2916 & $\\dots$ \\\\\n", "\t4 & a4 & -3.03031 & -2.0634 & -0.604634 & -1.05 & g4 & m & \\emph{missing} & -1.39308 & $\\dots$ \\\\\n", "\t5 & a5 & 2.18819 & 2.22425 & 1.72306 & 2.06 & g3 & m & \\emph{missing} & 0.267546 & $\\dots$ \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m5×11 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m ID \u001b[0m\u001b[1m y1 \u001b[0m\u001b[1m y2 \u001b[0m\u001b[1m y3 \u001b[0m\u001b[1m x1 \u001b[0m\u001b[1m x2 \u001b[0m\u001b[1m x3 \u001b[0m\u001b[1m dam \u001b[0m ⋯\n", "\u001b[1m \u001b[0m│\u001b[90m String \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m String \u001b[0m\u001b[90m String \u001b[0m\u001b[90m Stri\u001b[0m ⋯\n", "─────┼──────────────────────────────────────────────────────────────────────────\n", " 1 │ a1 -0.579682 -0.0582089 1.0456 0.77 g1 m \u001b[90m miss\u001b[0m ⋯\n", " 2 │ a2 -2.02246 -2.3843 -1.57196 -1.02 g2 m \u001b[90m miss\u001b[0m\n", " 3 │ a3 -1.48071 -0.421258 -0.272245 0.52 g1 f \u001b[90m miss\u001b[0m\n", " 4 │ a4 -3.03031 -2.0634 -0.604634 -1.05 g4 m \u001b[90m miss\u001b[0m\n", " 5 │ a5 2.18819 2.22425 1.72306 2.06 g3 m \u001b[90m miss\u001b[0m ⋯\n", "\u001b[36m 4 columns omitted\u001b[0m" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "first(phenotypes,5)" ] }, { "cell_type": "code", "execution_count": 74, "id": "c40dc157-731c-4570-80f1-755016a2ccb4", "metadata": {}, "outputs": [], "source": [ "# Step 3: Build Model Equations" ] }, { "cell_type": "code", "execution_count": 75, "id": "3316aeda-e615-4f39-b1c8-f4380b4545e8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes\"" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_equation =\"y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes\"" ] }, { "cell_type": "code", "execution_count": 76, "id": "59f07cee-a52a-47da-8e23-f515ec879541", "metadata": {}, "outputs": [], "source": [ "model = build_model(model_equation);" ] }, { "cell_type": "code", "execution_count": 77, "id": "1abd66a8-a528-423f-b587-5045bc8b77e6", "metadata": {}, "outputs": [], "source": [ "# Step 4: Set Factors or Covariates" ] }, { "cell_type": "code", "execution_count": 78, "id": "cfb60ded-35c9-4ea2-a952-f953d8e6a264", "metadata": {}, "outputs": [], "source": [ "set_covariate(model,\"x1\");" ] }, { "cell_type": "code", "execution_count": 79, "id": "8481056c-dcca-42eb-829a-555ea2239464", "metadata": {}, "outputs": [], "source": [ "# Step 5: Set Random or Fixed Effects" ] }, { "cell_type": "code", "execution_count": 80, "id": "918afc48-b356-48c6-aca9-84299a8cd4e9", "metadata": {}, "outputs": [], "source": [ "set_random(model,\"x2\");\n", "# set_random(model,\"ID dam\",pedigree);" ] }, { "cell_type": "code", "execution_count": 81, "id": "0d099ef4-e8b7-447d-8dca-f77ce1518536", "metadata": {}, "outputs": [], "source": [ "# # Step 6: Run Analysis" ] }, { "cell_type": "code", "execution_count": 82, "id": "b2f84a94-ff0f-4abb-800e-2bbba5c8e90a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32mThe folder results is created to save results.\u001b[39m\n", "\u001b[32mChecking genotypes...\u001b[39m\n", "\u001b[32mChecking phenotypes...\u001b[39m\n", "\u001b[32mIndividual IDs (strings) are provided in the first column of the phenotypic data.\u001b[39m\n", "\u001b[31mIn this complete genomic data (non-single-step) analyis, 40 phenotyped individuals are not genotyped. These are removed from the analysis.\u001b[39m\n", "\u001b[32mPredicted values for individuals of interest will be obtained as the summation of Any[] (Note that genomic data is always included for now).\u001b[39m\u001b[32mPhenotypes for 60 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.\u001b[39m\n", "\u001b[32mPrior information for genomic variance is not provided and is generated from the data.\u001b[39m\n", "\u001b[32mPrior information for residual variance is not provided and is generated from the data.\u001b[39m\n", "\u001b[32mPrior information for random effect variance is not provided and is generated from the data.\u001b[39m\n", "\n", "The prior for marker effects variance is calculated from the genetic variance and π.\n", "The mean of the prior for the marker effects variance is: 0.003622\n", "\n", "\n", "\n", "\u001b[0m\u001b[1mA Linear Mixed Model was build using model equations:\u001b[22m\n", "\n", "y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes\n", "\n", "\u001b[0m\u001b[1mModel Information:\u001b[22m\n", "\n", "Term C/F F/R nLevels\n", "intercept factor fixed 1\n", "x1 covariate fixed 1\n", "x2 factor random 5\n", "x2*x3 interaction fixed 10\n", "ID factor fixed 60\n", "dam factor fixed 29\n", "\n", "\u001b[0m\u001b[1mMCMC Information:\u001b[22m\n", "\n", "chain_length 100\n", "burnin 0\n", "starting_value true\n", "printout_frequency 101\n", "output_samples_frequency 1\n", "constraint false\n", "missing_phenotypes true\n", "update_priors_frequency 0\n", "seed false\n", "\n", "\u001b[0m\u001b[1mHyper-parameters Information:\u001b[22m\n", "\n", "random effect variances (y1:x2): [0.623]\n", "residual variances: 0.623\n", "\n", "\u001b[0m\u001b[1mGenomic Information:\u001b[22m\n", "\n", "complete genomic data (i.e., non-single-step analysis)\n", "\n", "Genomic Category genotypes\n", "Method BayesC\n", "genetic variances (genomic): 1.247\n", "marker effect variances: 0.004\n", "π 0.0\n", "estimatePi true\n", "estimateScale false\n", "\n", "\u001b[0m\u001b[1mDegree of freedom for hyper-parameters:\u001b[22m\n", "\n", "residual variances: 4.000\n", "random effect variances: 5.000\n", "marker effect variances: 4.000\n", "\n", "\n", "\n", "\u001b[32mThe file results/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_marker_effects_genotypes_y1.txt is created to save MCMC samples for marker_effects_genotypes_y1.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_marker_effects_variances_genotypes.txt is created to save MCMC samples for marker_effects_variances_genotypes.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_pi_genotypes.txt is created to save MCMC samples for pi_genotypes.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_y1.x2_variances.txt is created to save MCMC samples for y1:x2_variances.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_EBV_y1.txt is created to save MCMC samples for EBV_y1.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_genetic_variance.txt is created to save MCMC samples for genetic_variance.\u001b[39m\n", "\u001b[32mThe file results/MCMC_samples_heritability.txt is created to save MCMC samples for heritability.\u001b[39m\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mrunning MCMC ...100%|███████████████████████████████████| Time: 0:00:00\u001b[39m\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "\u001b[0m\u001b[1mThe version of Julia and Platform in use:\u001b[22m\n", "\n", "Julia Version 1.6.3\n", "Commit ae8452a9e0 (2021-09-23 17:34 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin19.5.0)\n", " CPU: Intel(R) Core(TM) i5-4278U CPU @ 2.60GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-11.0.1 (ORCJIT, haswell)\n", "\n", "\n", "\u001b[0m\u001b[1mThe analysis has finished. Results are saved in the returned \u001b[22m\u001b[0m\u001b[1mvariable and text files. MCMC samples are saved in text files.\u001b[22m\n", "\n", "\n" ] } ], "source": [ "out=runMCMC(model,phenotypes);" ] }, { "cell_type": "code", "execution_count": 83, "id": "ea0e0ff6-1209-4460-96bb-4868b8dff131", "metadata": {}, "outputs": [], "source": [ "# Check the results" ] }, { "cell_type": "code", "execution_count": 84, "id": "9447e66f-3994-4481-a1a8-33664dc28b23", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Dict{Any, Any} with 7 entries:\n", " \"EBV_y1\" => \u001b[1m60×3 DataFrame\u001b[0m…\n", " \"heritability\" => \u001b[1m1×3 DataFrame\u001b[0m…\n", " \"location parameters\" => \u001b[1m106×5 DataFrame\u001b[0m…\n", " \"residual variance\" => \u001b[1m1×3 DataFrame\u001b[0m…\n", " \"pi_genotypes\" => \u001b[1m1×3 DataFrame\u001b[0m…\n", " \"genetic_variance\" => \u001b[1m1×3 DataFrame\u001b[0m…\n", " \"marker effects genotypes\" => \u001b[1m1000×5 DataFrame\u001b[0m…" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out" ] }, { "cell_type": "code", "execution_count": 85, "id": "9fb96034-5c28-4441-811a-40169a22cada", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "KeySet for a Dict{Any, Any} with 7 entries. Keys:\n", " \"EBV_y1\"\n", " \"heritability\"\n", " \"location parameters\"\n", " \"residual variance\"\n", " \"pi_genotypes\"\n", " \"genetic_variance\"\n", " \"marker effects genotypes\"" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "keys(out)" ] }, { "cell_type": "code", "execution_count": 86, "id": "7e92dc90-a4cd-4638-a215-b1fbfa846d0d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

1 rows × 3 columns

CovarianceEstimateSD
AnyAnyAny
1y10.648410.141073
" ], "text/latex": [ "\\begin{tabular}{r|ccc}\n", "\t& Covariance & Estimate & SD\\\\\n", "\t\\hline\n", "\t& Any & Any & Any\\\\\n", "\t\\hline\n", "\t1 & y1 & 0.64841 & 0.141073 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m1×3 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m Covariance \u001b[0m\u001b[1m Estimate \u001b[0m\u001b[1m SD \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\n", "─────┼────────────────────────────────\n", " 1 │ y1 0.64841 0.141073" ] }, "execution_count": 86, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out[\"heritability\"]" ] }, { "cell_type": "code", "execution_count": 87, "id": "4411950c-8e69-41dd-a7ce-7e11f57f5df7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

1 rows × 3 columns

CovarianceEstimateSD
AnyAnyAny
1y1_y10.4149340.283966
" ], "text/latex": [ "\\begin{tabular}{r|ccc}\n", "\t& Covariance & Estimate & SD\\\\\n", "\t\\hline\n", "\t& Any & Any & Any\\\\\n", "\t\\hline\n", "\t1 & y1\\_y1 & 0.414934 & 0.283966 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m1×3 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m Covariance \u001b[0m\u001b[1m Estimate \u001b[0m\u001b[1m SD \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\n", "─────┼────────────────────────────────\n", " 1 │ y1_y1 0.414934 0.283966" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out[\"residual variance\"]" ] }, { "cell_type": "code", "execution_count": 88, "id": "fcbb88fd-817b-44d6-9b1a-ac0b81d9e43d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

1,000 rows × 5 columns

TraitMarker_IDEstimateSDModel_Frequency
AnyAnyAnyAnyAny
1y1m1-0.004063940.03648310.8
2y1m20.006340280.04698960.84
3y1m3-0.008253310.04364170.81
4y1m4-0.005313350.0426850.87
5y1m50.003256270.04590290.83
6y1m6-0.004930910.04112990.83
7y1m70.00824310.04353720.81
8y1m8-0.001290830.04568030.79
9y1m9-0.01038260.04091480.87
10y1m100.004375310.0439690.83
11y1m110.004331190.04424280.78
12y1m12-0.0009376530.04020710.8
13y1m13-0.005450290.04431240.84
14y1m140.003281230.04591520.9
15y1m15-0.006966940.04567810.85
16y1m165.13691e-50.04196890.83
17y1m170.00306850.04371310.83
18y1m180.005252230.04518650.83
19y1m19-0.008288770.04519340.79
20y1m200.00784350.04680860.85
21y1m210.001744290.04473790.87
22y1m22-0.007310890.04192740.81
23y1m23-0.014030.04232160.91
24y1m240.003640110.04904320.88
25y1m250.005946730.0543510.8
26y1m260.003592380.04269030.79
27y1m27-0.001752650.0426880.82
28y1m28-0.001658770.04196330.77
29y1m29-0.004368270.04280240.8
30y1m300.005067340.03983620.76
⋮⋮⋮⋮⋮⋮
" ], "text/latex": [ "\\begin{tabular}{r|ccccc}\n", "\t& Trait & Marker\\_ID & Estimate & SD & Model\\_Frequency\\\\\n", "\t\\hline\n", "\t& Any & Any & Any & Any & Any\\\\\n", "\t\\hline\n", "\t1 & y1 & m1 & -0.00406394 & 0.0364831 & 0.8 \\\\\n", "\t2 & y1 & m2 & 0.00634028 & 0.0469896 & 0.84 \\\\\n", "\t3 & y1 & m3 & -0.00825331 & 0.0436417 & 0.81 \\\\\n", "\t4 & y1 & m4 & -0.00531335 & 0.042685 & 0.87 \\\\\n", "\t5 & y1 & m5 & 0.00325627 & 0.0459029 & 0.83 \\\\\n", "\t6 & y1 & m6 & -0.00493091 & 0.0411299 & 0.83 \\\\\n", "\t7 & y1 & m7 & 0.0082431 & 0.0435372 & 0.81 \\\\\n", "\t8 & y1 & m8 & -0.00129083 & 0.0456803 & 0.79 \\\\\n", "\t9 & y1 & m9 & -0.0103826 & 0.0409148 & 0.87 \\\\\n", "\t10 & y1 & m10 & 0.00437531 & 0.043969 & 0.83 \\\\\n", "\t11 & y1 & m11 & 0.00433119 & 0.0442428 & 0.78 \\\\\n", "\t12 & y1 & m12 & -0.000937653 & 0.0402071 & 0.8 \\\\\n", "\t13 & y1 & m13 & -0.00545029 & 0.0443124 & 0.84 \\\\\n", "\t14 & y1 & m14 & 0.00328123 & 0.0459152 & 0.9 \\\\\n", "\t15 & y1 & m15 & -0.00696694 & 0.0456781 & 0.85 \\\\\n", "\t16 & y1 & m16 & 5.13691e-5 & 0.0419689 & 0.83 \\\\\n", "\t17 & y1 & m17 & 0.0030685 & 0.0437131 & 0.83 \\\\\n", "\t18 & y1 & m18 & 0.00525223 & 0.0451865 & 0.83 \\\\\n", "\t19 & y1 & m19 & -0.00828877 & 0.0451934 & 0.79 \\\\\n", "\t20 & y1 & m20 & 0.0078435 & 0.0468086 & 0.85 \\\\\n", "\t21 & y1 & m21 & 0.00174429 & 0.0447379 & 0.87 \\\\\n", "\t22 & y1 & m22 & -0.00731089 & 0.0419274 & 0.81 \\\\\n", "\t23 & y1 & m23 & -0.01403 & 0.0423216 & 0.91 \\\\\n", "\t24 & y1 & m24 & 0.00364011 & 0.0490432 & 0.88 \\\\\n", "\t25 & y1 & m25 & 0.00594673 & 0.054351 & 0.8 \\\\\n", "\t26 & y1 & m26 & 0.00359238 & 0.0426903 & 0.79 \\\\\n", "\t27 & y1 & m27 & -0.00175265 & 0.042688 & 0.82 \\\\\n", "\t28 & y1 & m28 & -0.00165877 & 0.0419633 & 0.77 \\\\\n", "\t29 & y1 & m29 & -0.00436827 & 0.0428024 & 0.8 \\\\\n", "\t30 & y1 & m30 & 0.00506734 & 0.0398362 & 0.76 \\\\\n", "\t$\\dots$ & $\\dots$ & $\\dots$ & $\\dots$ & $\\dots$ & $\\dots$ \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m1000×5 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m Trait \u001b[0m\u001b[1m Marker_ID \u001b[0m\u001b[1m Estimate \u001b[0m\u001b[1m SD \u001b[0m\u001b[1m Model_Frequency \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\u001b[90m Any \u001b[0m\n", "──────┼────────────────────────────────────────────────────────────\n", " 1 │ y1 m1 -0.00406394 0.0364831 0.8\n", " 2 │ y1 m2 0.00634028 0.0469896 0.84\n", " 3 │ y1 m3 -0.00825331 0.0436417 0.81\n", " 4 │ y1 m4 -0.00531335 0.042685 0.87\n", " 5 │ y1 m5 0.00325627 0.0459029 0.83\n", " 6 │ y1 m6 -0.00493091 0.0411299 0.83\n", " 7 │ y1 m7 0.0082431 0.0435372 0.81\n", " 8 │ y1 m8 -0.00129083 0.0456803 0.79\n", " 9 │ y1 m9 -0.0103826 0.0409148 0.87\n", " 10 │ y1 m10 0.00437531 0.043969 0.83\n", " 11 │ y1 m11 0.00433119 0.0442428 0.78\n", " ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮\n", " 991 │ y1 m991 -0.00874085 0.0377388 0.83\n", " 992 │ y1 m992 -0.0110955 0.0391157 0.78\n", " 993 │ y1 m993 -0.00174481 0.0371269 0.78\n", " 994 │ y1 m994 -0.00463736 0.0412571 0.79\n", " 995 │ y1 m995 -0.00120571 0.0414865 0.81\n", " 996 │ y1 m996 0.00929234 0.0472444 0.8\n", " 997 │ y1 m997 -0.00156771 0.0442376 0.89\n", " 998 │ y1 m998 -0.000405303 0.0472067 0.88\n", " 999 │ y1 m999 -0.00260627 0.0475014 0.82\n", " 1000 │ y1 m1000 0.00648095 0.0479945 0.85\n", "\u001b[36m 979 rows omitted\u001b[0m" ] }, "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ "out[\"marker effects genotypes\"]" ] }, { "cell_type": "code", "execution_count": 89, "id": "e86fe35d-40ea-4745-ac1d-6804d932bbff", "metadata": {}, "outputs": [], "source": [ "# Bayesian marker effect SEM-GWAS - https://github.com/reworkhow/JWAS.jl/wiki/Integrating-Phenotypic-Causal-Networks-in-GWAS" ] }, { "cell_type": "code", "execution_count": 90, "id": "d91db5e0-e345-4ebf-9dc2-2ccb77d019d5", "metadata": {}, "outputs": [], "source": [ "using JWAS,DataFrames,CSV,Statistics,JWAS.Datasets" ] }, { "cell_type": "code", "execution_count": 91, "id": "b8617842-b1d5-4844-ae60-9534b98a7129", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32mThe delimiter in pedigree.csv is ','.\u001b[39m\n", "Pedigree information:\n", "#individuals: 100\n", "#sires: 39\n", "#dams: 38\n", "#founders: 20\n", "\u001b[32mThe delimiterd in genotypes.csv is ','. \u001b[39m\u001b[32mThe header (marker IDs) is provided in genotypes.csv.\u001b[39m\n", "Genotype informatin:\n", "#markers: 1000; #individuals: 60\n" ] }, { "data": { "text/html": [ "

5 rows × 11 columns (omitted printing of 2 columns)

IDy1y2y3x1x2x3dambv1
StringFloat64Float64Float64Float64StringStringString?Float64
1a1-0.579682-0.05820891.04560.77g1mmissing-1.07018
2a2-2.02246-2.3843-1.57196-1.02g2mmissing-1.03328
3a3-1.48071-0.421258-0.2722450.52g1fmissing-1.2916
4a4-3.03031-2.0634-0.604634-1.05g4mmissing-1.39308
5a52.188192.224251.723062.06g3mmissing0.267546
" ], "text/latex": [ "\\begin{tabular}{r|cccccccccc}\n", "\t& ID & y1 & y2 & y3 & x1 & x2 & x3 & dam & bv1 & \\\\\n", "\t\\hline\n", "\t& String & Float64 & Float64 & Float64 & Float64 & String & String & String? & Float64 & \\\\\n", "\t\\hline\n", "\t1 & a1 & -0.579682 & -0.0582089 & 1.0456 & 0.77 & g1 & m & \\emph{missing} & -1.07018 & $\\dots$ \\\\\n", "\t2 & a2 & -2.02246 & -2.3843 & -1.57196 & -1.02 & g2 & m & \\emph{missing} & -1.03328 & $\\dots$ \\\\\n", "\t3 & a3 & -1.48071 & -0.421258 & -0.272245 & 0.52 & g1 & f & \\emph{missing} & -1.2916 & $\\dots$ \\\\\n", "\t4 & a4 & -3.03031 & -2.0634 & -0.604634 & -1.05 & g4 & m & \\emph{missing} & -1.39308 & $\\dots$ \\\\\n", "\t5 & a5 & 2.18819 & 2.22425 & 1.72306 & 2.06 & g3 & m & \\emph{missing} & 0.267546 & $\\dots$ \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m5×11 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m ID \u001b[0m\u001b[1m y1 \u001b[0m\u001b[1m y2 \u001b[0m\u001b[1m y3 \u001b[0m\u001b[1m x1 \u001b[0m\u001b[1m x2 \u001b[0m\u001b[1m x3 \u001b[0m\u001b[1m dam \u001b[0m ⋯\n", "\u001b[1m \u001b[0m│\u001b[90m String \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m String \u001b[0m\u001b[90m String \u001b[0m\u001b[90m Stri\u001b[0m ⋯\n", "─────┼──────────────────────────────────────────────────────────────────────────\n", " 1 │ a1 -0.579682 -0.0582089 1.0456 0.77 g1 m \u001b[90m miss\u001b[0m ⋯\n", " 2 │ a2 -2.02246 -2.3843 -1.57196 -1.02 g2 m \u001b[90m miss\u001b[0m\n", " 3 │ a3 -1.48071 -0.421258 -0.272245 0.52 g1 f \u001b[90m miss\u001b[0m\n", " 4 │ a4 -3.03031 -2.0634 -0.604634 -1.05 g4 m \u001b[90m miss\u001b[0m\n", " 5 │ a5 2.18819 2.22425 1.72306 2.06 g3 m \u001b[90m miss\u001b[0m ⋯\n", "\u001b[36m 4 columns omitted\u001b[0m" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Step 2: Read data \n", "phenofile = dataset(\"phenotypes.csv\")\n", "pedfile = dataset(\"pedigree.csv\")\n", "genofile = dataset(\"genotypes.csv\")\n", "phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=[\"NA\"])\n", "pedigree = get_pedigree(pedfile,separator=\",\",header=true);\n", "genotypes = get_genotypes(genofile,separator=',',method=\"BayesC\");\n", "first(phenotypes,5)" ] }, { "cell_type": "code", "execution_count": 92, "id": "06e5ab4c-e474-466d-9822-d3068e031504", "metadata": {}, "outputs": [], "source": [ "# Step 3: Build Model Equations\n", "\n", "model_equation =\"y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes\n", " y2 = intercept + x1 + x2 + ID + genotypes\n", " y3 = intercept + x1 + ID + genotypes\";\n", "model = build_model(model_equation);" ] }, { "cell_type": "code", "execution_count": 93, "id": "ae37f7e7-982b-4e28-b0d9-6fb381c53aca", "metadata": {}, "outputs": [], "source": [ "# Step 4: Set Factors or Covariates\n", "set_covariate(model,\"x1\");" ] }, { "cell_type": "code", "execution_count": 94, "id": "14a06068-ce2c-46a0-bbe3-3e2217a02713", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[32mx2 is not found in model equation 3.\u001b[39m\n", "\u001b[32mdam is not found in model equation 2.\u001b[39m\n", "\u001b[32mdam is not found in model equation 3.\u001b[39m\n" ] } ], "source": [ "# Step 5: Set Random or Fixed Effects\n", "set_random(model,\"x2\");\n", "set_random(model,\"ID dam\",pedigree);" ] }, { "cell_type": "code", "execution_count": 95, "id": "1fc0b077-5814-4897-bcf4-b9c43d734e7a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[31mThe folder results already exists.\u001b[39m\n", "\u001b[32mThe folder results1 is created to save results.\u001b[39m\n", "\u001b[32mChecking pedigree...\u001b[39m\n", "\u001b[32mChecking genotypes...\u001b[39m\n", "\u001b[32mChecking phenotypes...\u001b[39m\n", "\u001b[32mIndividual IDs (strings) are provided in the first column of the phenotypic data.\u001b[39m\n", "\u001b[31mIn this complete genomic data (non-single-step) analyis, 40 phenotyped individuals are not genotyped. These are removed from the analysis.\u001b[39m\n", "\u001b[32mPredicted values for individuals of interest will be obtained as the summation of Any[\"y1:ID\", \"y2:ID\", \"y3:ID\", \"y1:dam\"] (Note that genomic data is always included for now).\u001b[39m\u001b[32mPhenotypes for 60 observations are used in the analysis.These individual IDs are saved in the file IDs_for_individuals_with_phenotypes.txt.\u001b[39m\n", "\u001b[32mPrior information for genomic variance is not provided and is generated from the data.\u001b[39m\n", "\u001b[32mPrior information for residual variance is not provided and is generated from the data.\u001b[39m\n", "\u001b[32mPrior information for random effect variance is not provided and is generated from the data.\u001b[39m\n", "\u001b[32mPrior information for random effect variance is not provided and is generated from the data.\u001b[39m\n", "\n", "\u001b[0mPi (Π) is not provided.\n", "\u001b[0mPi (Π) is generated assuming all markers have effects on all traits.\n", "\n", "The prior for marker effects covariance matrix is calculated from genetic covariance matrix and Π.\n", "The mean of the prior for the marker effects covariance matrix is:\n", " 0.001811 0.0 0.0\n", " 0.0 0.001849 0.0\n", " 0.0 0.0 0.001429\n", "\n", "\n", "\n", "\u001b[0m\u001b[1mA Linear Mixed Model was build using model equations:\u001b[22m\n", "\n", "y1 = intercept + x1 + x2 + x2*x3 + ID + dam + genotypes\n", "y2 = intercept + x1 + x2 + ID + genotypes\n", "y3 = intercept + x1 + ID + genotypes\n", "\n", "\u001b[0m\u001b[1mModel Information:\u001b[22m\n", "\n", "Term C/F F/R nLevels\n", "intercept factor fixed 1\n", "x1 covariate fixed 1\n", "x2 factor random 5\n", "x2*x3 interaction fixed 10\n", "ID factor random 100\n", "dam factor random 100\n", "\n", "\u001b[0m\u001b[1mMCMC Information:\u001b[22m\n", "\n", "chain_length 100\n", "burnin 0\n", "starting_value true\n", "printout_frequency 101\n", "output_samples_frequency 1\n", "constraint true\n", "missing_phenotypes false\n", "update_priors_frequency 0\n", "seed false\n", "\n", "\u001b[0m\u001b[1mHyper-parameters Information:\u001b[22m\n", "\n", "random effect variances (y1:x2,y2:x2):\n", " 0.623 0.0\n", " 0.0 0.637\n", "random effect variances (y1:ID,y2:ID,y3:ID,y1:dam):\n", " 0.623 0.0 0.0 0.0\n", " 0.0 0.637 0.0 0.0\n", " 0.0 0.0 0.492 0.0\n", " 0.0 0.0 0.0 0.623\n", "genetic variances (polygenic):\n", " 0.623 0.0 0.0 0.0\n", " 0.0 0.637 0.0 0.0\n", " 0.0 0.0 0.492 0.0\n", " 0.0 0.0 0.0 0.623\n", "residual variances: \n", " 0.623 0.0 0.0\n", " 0.0 0.637 0.0\n", " 0.0 0.0 0.492\n", "\n", "\u001b[0m\u001b[1mGenomic Information:\u001b[22m\n", "\n", "complete genomic data (i.e., non-single-step analysis)\n", "\n", "Genomic Category genotypes\n", "Method BayesC\n", "genetic variances (genomic): \n", " 0.623 0.0 0.0\n", " 0.0 0.637 0.0\n", " 0.0 0.0 0.492\n", "marker effect variances: \n", " 0.002 0.0 0.0\n", " 0.0 0.002 0.0\n", " 0.0 0.0 0.001\n", "\n", "Π: (Y(yes):included; N(no):excluded)\n", "\n", "[\"y1\", \"y2\", \"y3\"] probability\n", "[\"Y\", \"Y\", \"Y\"] 1.0\n", "[\"N\", \"N\", \"N\"] 0.0\n", "[\"N\", \"N\", \"Y\"] 0.0\n", "[\"Y\", \"Y\", \"N\"] 0.0\n", "[\"N\", \"Y\", \"Y\"] 0.0\n", "[\"Y\", \"N\", \"N\"] 0.0\n", "[\"Y\", \"N\", \"Y\"] 0.0\n", "[\"N\", \"Y\", \"N\"] 0.0\n", "\n", "estimatePi true\n", "estimateScale false\n", "\n", "\u001b[0m\u001b[1mDegree of freedom for hyper-parameters:\u001b[22m\n", "\n", "residual variances: 7.000\n", "random effect variances: 6.000\n", "polygenic effect variances: 8.000\n", "marker effect variances: 4.000\n", "\n", "\n", "\n", "\u001b[32mThe file results1/MCMC_samples_residual_variance.txt is created to save MCMC samples for residual_variance.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_polygenic_effects_variance.txt is created to save MCMC samples for polygenic_effects_variance.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_marker_effects_genotypes_y1.txt is created to save MCMC samples for marker_effects_genotypes_y1.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_marker_effects_genotypes_y2.txt is created to save MCMC samples for marker_effects_genotypes_y2.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_marker_effects_genotypes_y3.txt is created to save MCMC samples for marker_effects_genotypes_y3.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_marker_effects_variances_genotypes.txt is created to save MCMC samples for marker_effects_variances_genotypes.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_pi_genotypes.txt is created to save MCMC samples for pi_genotypes.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_y1.x2_y2.x2_variances.txt is created to save MCMC samples for y1:x2_y2:x2_variances.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_y1.ID_y2.ID_y3.ID_y1.dam_variances.txt is created to save MCMC samples for y1:ID_y2:ID_y3:ID_y1:dam_variances.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_EBV_y1.txt is created to save MCMC samples for EBV_y1.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_EBV_y2.txt is created to save MCMC samples for EBV_y2.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_EBV_y3.txt is created to save MCMC samples for EBV_y3.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_genetic_variance.txt is created to save MCMC samples for genetic_variance.\u001b[39m\n", "\u001b[32mThe file results1/MCMC_samples_heritability.txt is created to save MCMC samples for heritability.\u001b[39m\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mrunning MCMC ...100%|███████████████████████████████████| Time: 0:00:02\u001b[39m\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "\u001b[0m\u001b[1mThe version of Julia and Platform in use:\u001b[22m\n", "\n", "Julia Version 1.6.3\n", "Commit ae8452a9e0 (2021-09-23 17:34 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin19.5.0)\n", " CPU: Intel(R) Core(TM) i5-4278U CPU @ 2.60GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-11.0.1 (ORCJIT, haswell)\n", "\n", "\n", "\u001b[0m\u001b[1mThe analysis has finished. Results are saved in the returned \u001b[22m\u001b[0m\u001b[1mvariable and text files. MCMC samples are saved in text files.\u001b[22m\n", "\n", "\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n", "Compute the model frequency for each marker (the probability the marker is included in the model).\n" ] } ], "source": [ "# Step 6: Run Analysis\n", "# 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] (column index affects row index, and a lower triangular matrix is required) for\n", "# trait 1 -> trait 2 and trait 1 -> trait 3, phenotypic causal networks will be incorporated using structure equation models.\n", "my_structure = [0.0 0.0 0.0\n", " 1.0 0.0 0.0\n", " 1.0 0.0 0.0]\n", "out=runMCMC(model,phenotypes,causal_structure=my_structure);" ] }, { "cell_type": "code", "execution_count": 96, "id": "6de3f9d7-c551-49f7-8a4b-7bb246c20cf2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compute the model frequency for each marker (the probability the marker is included in the model).\n" ] }, { "data": { "text/html": [ "

1,000 rows × 2 columns

marker_IDmodelfrequency
Abstrac…Float64
1m10.91
2m20.94
3m30.88
4m40.96
5m50.94
6m61.0
7m70.98
8m80.96
9m90.8
10m100.96
11m110.95
12m120.99
13m130.97
14m140.89
15m150.96
16m160.81
17m170.99
18m180.96
19m190.96
20m200.98
21m211.0
22m220.96
23m231.0
24m240.93
25m250.96
26m260.94
27m270.85
28m280.96
29m290.98
30m300.96
⋮⋮⋮
" ], "text/latex": [ "\\begin{tabular}{r|cc}\n", "\t& marker\\_ID & modelfrequency\\\\\n", "\t\\hline\n", "\t& Abstrac… & Float64\\\\\n", "\t\\hline\n", "\t1 & m1 & 0.91 \\\\\n", "\t2 & m2 & 0.94 \\\\\n", "\t3 & m3 & 0.88 \\\\\n", "\t4 & m4 & 0.96 \\\\\n", "\t5 & m5 & 0.94 \\\\\n", "\t6 & m6 & 1.0 \\\\\n", "\t7 & m7 & 0.98 \\\\\n", "\t8 & m8 & 0.96 \\\\\n", "\t9 & m9 & 0.8 \\\\\n", "\t10 & m10 & 0.96 \\\\\n", "\t11 & m11 & 0.95 \\\\\n", "\t12 & m12 & 0.99 \\\\\n", "\t13 & m13 & 0.97 \\\\\n", "\t14 & m14 & 0.89 \\\\\n", "\t15 & m15 & 0.96 \\\\\n", "\t16 & m16 & 0.81 \\\\\n", "\t17 & m17 & 0.99 \\\\\n", "\t18 & m18 & 0.96 \\\\\n", "\t19 & m19 & 0.96 \\\\\n", "\t20 & m20 & 0.98 \\\\\n", "\t21 & m21 & 1.0 \\\\\n", "\t22 & m22 & 0.96 \\\\\n", "\t23 & m23 & 1.0 \\\\\n", "\t24 & m24 & 0.93 \\\\\n", "\t25 & m25 & 0.96 \\\\\n", "\t26 & m26 & 0.94 \\\\\n", "\t27 & m27 & 0.85 \\\\\n", "\t28 & m28 & 0.96 \\\\\n", "\t29 & m29 & 0.98 \\\\\n", "\t30 & m30 & 0.96 \\\\\n", "\t$\\dots$ & $\\dots$ & $\\dots$ \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m1000×2 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m marker_ID \u001b[0m\u001b[1m modelfrequency \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Abstract… \u001b[0m\u001b[90m Float64 \u001b[0m\n", "──────┼───────────────────────────\n", " 1 │ m1 0.91\n", " 2 │ m2 0.94\n", " 3 │ m3 0.88\n", " 4 │ m4 0.96\n", " 5 │ m5 0.94\n", " 6 │ m6 1.0\n", " 7 │ m7 0.98\n", " 8 │ m8 0.96\n", " 9 │ m9 0.8\n", " 10 │ m10 0.96\n", " 11 │ m11 0.95\n", " ⋮ │ ⋮ ⋮\n", " 991 │ m991 1.0\n", " 992 │ m992 0.97\n", " 993 │ m993 0.89\n", " 994 │ m994 0.87\n", " 995 │ m995 0.95\n", " 996 │ m996 0.93\n", " 997 │ m997 0.98\n", " 998 │ m998 0.97\n", " 999 │ m999 0.99\n", " 1000 │ m1000 0.97\n", "\u001b[36m 979 rows omitted\u001b[0m" ] }, "execution_count": 96, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# GWAS of Direct Marker Effects on Trait y2\n", "\n", "#Compute the model frequency for each marker (the probability the marker is included in the model).\n", "marker_effects_file=\"results1/MCMC_samples_marker_effects_genotypes_y2.txt\"\n", "out = GWAS(marker_effects_file,header=true)\n", "\n" ] }, { "cell_type": "code", "execution_count": 97, "id": "1ed78d85-cba8-43b1-84cb-bf57da27bf66", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compute the model frequency for each marker (the probability the marker is included in the model).\n" ] }, { "data": { "text/html": [ "

1,000 rows × 2 columns

marker_IDmodelfrequency
Abstrac…Float64
1m10.92
2 m20.92
3 m30.95
4 m41.0
5 m50.91
6 m60.96
7 m70.99
8 m80.98
9 m90.9
10 m100.99
11 m110.94
12 m120.98
13 m130.95
14 m140.91
15 m150.94
16 m160.88
17 m170.99
18 m180.99
19 m190.96
20 m200.99
21 m210.95
22 m220.96
23 m231.0
24 m240.92
25 m250.94
26 m260.96
27 m270.93
28 m280.95
29 m290.98
30 m300.94
⋮⋮⋮
" ], "text/latex": [ "\\begin{tabular}{r|cc}\n", "\t& marker\\_ID & modelfrequency\\\\\n", "\t\\hline\n", "\t& Abstrac… & Float64\\\\\n", "\t\\hline\n", "\t1 & m1 & 0.92 \\\\\n", "\t2 & m2 & 0.92 \\\\\n", "\t3 & m3 & 0.95 \\\\\n", "\t4 & m4 & 1.0 \\\\\n", "\t5 & m5 & 0.91 \\\\\n", "\t6 & m6 & 0.96 \\\\\n", "\t7 & m7 & 0.99 \\\\\n", "\t8 & m8 & 0.98 \\\\\n", "\t9 & m9 & 0.9 \\\\\n", "\t10 & m10 & 0.99 \\\\\n", "\t11 & m11 & 0.94 \\\\\n", "\t12 & m12 & 0.98 \\\\\n", "\t13 & m13 & 0.95 \\\\\n", "\t14 & m14 & 0.91 \\\\\n", "\t15 & m15 & 0.94 \\\\\n", "\t16 & m16 & 0.88 \\\\\n", "\t17 & m17 & 0.99 \\\\\n", "\t18 & m18 & 0.99 \\\\\n", "\t19 & m19 & 0.96 \\\\\n", "\t20 & m20 & 0.99 \\\\\n", "\t21 & m21 & 0.95 \\\\\n", "\t22 & m22 & 0.96 \\\\\n", "\t23 & m23 & 1.0 \\\\\n", "\t24 & m24 & 0.92 \\\\\n", "\t25 & m25 & 0.94 \\\\\n", "\t26 & m26 & 0.96 \\\\\n", "\t27 & m27 & 0.93 \\\\\n", "\t28 & m28 & 0.95 \\\\\n", "\t29 & m29 & 0.98 \\\\\n", "\t30 & m30 & 0.94 \\\\\n", "\t$\\dots$ & $\\dots$ & $\\dots$ \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m1000×2 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m marker_ID \u001b[0m\u001b[1m modelfrequency \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Abstract… \u001b[0m\u001b[90m Float64 \u001b[0m\n", "──────┼───────────────────────────\n", " 1 │ m1 0.92\n", " 2 │ m2 0.92\n", " 3 │ m3 0.95\n", " 4 │ m4 1.0\n", " 5 │ m5 0.91\n", " 6 │ m6 0.96\n", " 7 │ m7 0.99\n", " 8 │ m8 0.98\n", " 9 │ m9 0.9\n", " 10 │ m10 0.99\n", " 11 │ m11 0.94\n", " ⋮ │ ⋮ ⋮\n", " 991 │ m991 0.98\n", " 992 │ m992 0.99\n", " 993 │ m993 0.91\n", " 994 │ m994 0.91\n", " 995 │ m995 0.97\n", " 996 │ m996 0.92\n", " 997 │ m997 0.98\n", " 998 │ m998 0.97\n", " 999 │ m999 0.99\n", " 1000 │ m1000 0.98\n", "\u001b[36m 979 rows omitted\u001b[0m" ] }, "execution_count": 97, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# GWAS of Indirect Marker Effects on Trait y2\n", "\n", "#Compute the model frequency for each marker (the probability the marker is included in the model).\n", "marker_effects_file=\"results1/MCMC_samples_indirect_marker_effects_genotypes_y2.txt\"\n", "out=GWAS(marker_effects_file,header=true)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1d82b75c-0dc7-484d-8905-10fae9527464", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.3", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.3" } }, "nbformat": 4, "nbformat_minor": 5 }