Merge pull request #57 from GregFa/main

Update example and corrected subset functions

README.md

@@ -98,20 +98,27 @@ Load bxd data using the function `get_geneticstudydata()`:
 data = get_geneticstudydata(file);
 ```
 
+```julia
+# The current version of `BigRiverQTL` does not have imputation functions.
+# Remove the  missing data
+data = get_data_completecases(data);
+```
 
 ```julia
 # Data types
+# gmap contains 
 # makers info 
 gInfo = data.gmap;
-# pehnotype info 
+
+# phenotype info 
 pInfo = data.phenocov;
 # phenotype values 
 pheno = data.pheno.val;
 
-# We can get the genotype matrix using the following command:
-geno = reduce(hcat, data.geno.val);
-# For computing reasons, we need to convert the geno matrix in Float64
-geno_processed = convert(Array{Float64}, geno);
+# We can get the genotype matrix using the following command.
+# For computing reasons, we need to convert the geno matrix in Float64.
+# One way to do it is to multiply by 1.0
+geno = reduce(hcat, data.geno.val).*1.0;
 ```
 
 #### Preprocessing
@@ -123,9 +130,9 @@ geno_processed = convert(Array{Float64}, geno);
 #################
 traitID = 1112;
 pheno_y = pheno[:, traitID];
-pheno_y2=ones(length(pheno_y));
-idx_nothing = findall(x->x!=nothing,pheno_y)
-pheno_y2[idx_nothing]=pheno_y[idx_nothing];
+pheno_y2 = ones(length(pheno_y));
+idx_not_missing = findall(!ismissing, pheno_y)
+pheno_y2[idx_not_missing] = pheno_y[idx_not_missing];
 ```
 
 #### Kinship
@@ -135,7 +142,7 @@ pheno_y2[idx_nothing]=pheno_y[idx_nothing];
 ###########
 # Kinship #
 ###########
-kinship = kinship_gs(geno_processed,.99);
+kinship = kinship_gs(geno,.99);
 ```
 
 #### Scan
@@ -148,7 +155,7 @@ kinship = kinship_gs(geno_processed,.99);
 
 single_results_perms = scan(
 	pheno_y2,
-	geno_processed,
+	geno,
 	kinship;
 	permutation_test = true,
 	nperms = 1000,
@@ -169,10 +176,10 @@ single_results_perms = scan(
 plot_QTL(single_results_perms, gInfo, mbColname = "Pos")
 
 ```
-![image](images/QTL_example.png)
+![image](images/QTL_thrs_example.png)
 
 ```julia
 # Manhattan plots
 plot_manhattan(single_results_perms, gInfo, mbColname = "Pos")
 ```
-![image](images/manhattan_example.png)
+![image](images/manhattan_thrs_example.png)

example/example_qtl.ipynb

@@ -107,6 +107,18 @@
     "data = get_geneticstudydata(file);"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "168ea215-ff1d-41a7-92fe-6f75b74f56dd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The current version of `BigRiverQTL` does not have imputation functions.\n",
+    "# Remove the  missing data\n",
+    "data = get_data_completecases(data);"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -118,19 +130,16 @@
     "# gmap contains \n",
     "# makers info \n",
     "gInfo = data.gmap;\n",
-    "# pehnotype info \n",
+    "\n",
+    "# phenotype info \n",
     "pInfo = data.phenocov;\n",
     "# phenotype values \n",
     "pheno = data.pheno.val;\n",
     "\n",
-    "\n",
-    "# We can get the genotype matrix using the following command:\n",
-    "geno = reduce(hcat, data.geno.val);\n",
-    "# For computing reasons, we need to convert the geno matrix in Float64\n",
-    "geno_processed = geno .- 1.0 |>\n",
-    "    x -> replace(x, missing => 0.5) |>\n",
-    "\n",
-    "# geno_processed = convert(Matrix{Float64}, geno_processed);\n"
+    "# We can get the genotype matrix using the following command.\n",
+    "# For computing reasons, we need to convert the geno matrix in Float64.\n",
+    "# One way to do it is to multiply by 1.0\n",
+    "geno = reduce(hcat, data.geno.val).*1.0;"
    ]
   },
   {
@@ -153,9 +162,9 @@
     "#################\n",
     "traitID = 1112;\n",
     "pheno_y = pheno[:, traitID];\n",
-    "pheno_y2=ones(length(pheno_y));\n",
-    "idx_nothing = findall(x->x!=nothing,pheno_y)\n",
-    "pheno_y2[idx_nothing]=pheno_y[idx_nothing];"
+    "pheno_y2 = ones(length(pheno_y));\n",
+    "idx_not_missing = findall(!ismissing, pheno_y)\n",
+    "pheno_y2[idx_not_missing] = pheno_y[idx_not_missing];"
    ]
   },
   {
@@ -176,7 +185,7 @@
     "###########\n",
     "# Kinship #\n",
     "###########\n",
-    "kinship = kinship_gs(geno_processed,.99);"
+    "kinship = kinship_gs(geno,.99);"
    ]
   },
   {
@@ -200,7 +209,7 @@
     "\n",
     "single_results_perms = scan(\n",
     "\tpheno_y2,\n",
-    "\tgeno_processed,\n",
+    "\tgeno,\n",
     "\tkinship;\n",
     "\tpermutation_test = true,\n",
     "\tnperms = 1000,\n",
@@ -235,7 +244,8 @@
     "#########\n",
     "\n",
     "# QTL plots\n",
-    "plot_QTL(single_results_perms, gInfo, mbColname = \"Pos\")\n"
+    "plot_QTL(single_results_perms, gInfo, mbColname = \"Pos\")\n",
+    "savefig(joinpath(@__DIR__, \"..\", \"images\", \"QTL_thrs_example.png\"))"
    ]
   },
   {
@@ -246,7 +256,8 @@
    "outputs": [],
    "source": [
     "# Manhattan plots\n",
-    "plot_manhattan(single_results_perms, gInfo, mbColname = \"Pos\")"
+    "plot_manhattan(single_results_perms, gInfo, mbColname = \"Pos\")\n",
+    "savefig(joinpath(@__DIR__, \"..\", \"images\", \"manhattan_thrs_example.png\"))"
    ]
   }
  ],

example/example_qtl.jl

@@ -48,25 +48,25 @@ file = joinpath(data_dir, "bxd.json");
 data = get_geneticstudydata(file);
 
 
+# The current version of `BigRiverQTL` does not have imputation functions.
+# Remove the  missing data
+data = get_data_completecases(data);
+
 # +
 # Data types
 # gmap contains 
 # makers info 
 gInfo = data.gmap;
-# pehnotype info 
+
+# phenotype info 
 pInfo = data.phenocov;
 # phenotype values 
 pheno = data.pheno.val;
 
-
-# We can get the genotype matrix using the following command:
-geno = reduce(hcat, data.geno.val);
-# For computing reasons, we need to convert the geno matrix in Float64
-geno_processed = geno .- 1.0 |>
-    x -> replace(x, missing => 0.5) |>
-
-# geno_processed = convert(Matrix{Float64}, geno_processed);
-
+# We can get the genotype matrix using the following command.
+# For computing reasons, we need to convert the geno matrix in Float64.
+# One way to do it is to multiply by 1.0
+geno = reduce(hcat, data.geno.val).*1.0;
 # -
 
 # #### Preprocessing
@@ -76,16 +76,16 @@ geno_processed = geno .- 1.0 |>
 #################
 traitID = 1112;
 pheno_y = pheno[:, traitID];
-pheno_y2=ones(length(pheno_y));
-idx_nothing = findall(x->x!=nothing,pheno_y)
-pheno_y2[idx_nothing]=pheno_y[idx_nothing];
+pheno_y2 = ones(length(pheno_y));
+idx_not_missing = findall(!ismissing, pheno_y)
+pheno_y2[idx_not_missing] = pheno_y[idx_not_missing];
 
 # #### Kinship
 
 ###########
 # Kinship #
 ###########
-kinship = kinship_gs(geno_processed,.99);
+kinship = kinship_gs(geno,.99);
 
 # #### Scan
 
@@ -96,7 +96,7 @@ kinship = kinship_gs(geno_processed,.99);
 
 single_results_perms = scan(
 	pheno_y2,
-	geno_processed,
+	geno,
 	kinship;
 	permutation_test = true,
 	nperms = 1000,
@@ -114,8 +114,9 @@ single_results_perms = scan(
 
 # QTL plots
 plot_QTL(single_results_perms, gInfo, mbColname = "Pos")
-
+savefig(joinpath(@__DIR__, "..", "images", "QTL_thrs_example.png"))
 # -
 
 # Manhattan plots
 plot_manhattan(single_results_perms, gInfo, mbColname = "Pos")
+savefig(joinpath(@__DIR__, "..", "images", "manhattan_thrs_example.png"))

src/BigRiverQTL.jl

@@ -89,7 +89,9 @@ module BigRiverQTL
     #########
     # Utils #
     #########
-    include("utils/wrangling_utils.jl")
-    export get_geno_completecases, summary_missing, select_sample, select_marker
+    include("utils/subset_utils.jl")
+    export  select_sample, select_marker
 
+    include("utils/missing_utils.jl")
+    export get_geno_completecases, get_data_completecases, summary_missing
 end