@@ -18,6 +18,7 @@ library("sommer")
18
18
library("rrBLUP")
19
19
library("ggplot2")
20
20
library("tidyverse")
21
+ library("factoextra")
21
22
library("data.table")
22
23
23
24
knitr::opts_chunk$set(echo = FALSE)
@@ -93,7 +94,7 @@ print(g)
93
94
### ANOVA table
94
95
95
96
``` {r echo=TRUE}
96
- kable( anova(fit) )
97
+ anova(fit)
97
98
```
98
99
99
100
### Add population structure
@@ -118,11 +119,11 @@ Xscaled = scale(X, center = TRUE, scale = TRUE)
118
119
119
120
## uncomment if you want to recalculate PCs from SNP data
120
121
## (this may take a while)
121
- # res <- prcomp(X, rank. = 3)
122
+ res <- prcomp(X, rank. = 3)
122
123
123
124
## we can load precomputed PCs from the RData object in /data
124
- load("../data/pca.RData")
125
- names(res)
125
+ # load("../data/pca.RData")
126
+ # names(res)
126
127
127
128
npcs = length(res$sdev)
128
129
print(paste("QUESTION: The total number of PCs is", npcs, ": can you guess why?"))
@@ -144,14 +145,15 @@ sum(variance[c(1:3)])
144
145
```
145
146
146
147
``` {r echo=TRUE}
147
- library("factoextra")
148
+ # library("factoextra")
148
149
fviz_eig(res, ncp = 25)
149
150
```
150
151
151
152
#### Use principal components to model population structure
152
153
153
154
``` {r}
154
- head(res$x[,c(1,2,3,123,124,125)])
155
+ head(res$x[,c(1,2,3)])
156
+
155
157
## we use the first 3 PCs
156
158
df <- cbind.data.frame(df,res$x[,c(1:3)])
157
159
```
@@ -255,7 +257,6 @@ model1_x <- rrBLUP::GWAS(
255
257
)
256
258
```
257
259
258
-
259
260
``` {r}
260
261
# model1_x <- sommer::GWAS(
261
262
# fixed = phenotype~1,
0 commit comments