/
05-CrossValidation.Rmd
372 lines (322 loc) · 17 KB
/
05-CrossValidation.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
---
title: "Parent-wise cross-validation to check the accuracy of predicting cross (co)-variances"
author: "Marnin Wolfe"
date: "2021-May-14"
output:
workflowr::wflow_html:
toc: true
editor_options:
chunk_output_type: inline
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = F,
eval = FALSE, # <- NOTE THAT EVAL SET TO FALSE!
tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```
# Previous step
4. [Preprocess data files](04-PreprocessDataFiles.html): Prepare haplotype and dosage matrices, pedigree and BLUPs, genetic map *and* recombination frequency matrix, for use in predictions.
# Automating cross-validation
In the manuscript, the cross-validation is documented many pages and scripts, [documented here](https://wolfemd.github.io/PredictOutbredCrossVar/).
For ongoing GS, I have a function `runCrossVal()` that manages all inputs and outputs easy to work with pre-computed accuracies.
Goal here is to make a function: `runParentWiseCrossVal()`, or at least make progress towards developing one.
*However*, for computational reasons, I imagine it might still be best to separate the task into a few functions.
My goal is to simplify and integrate into the pipeline used for NextGen Cassava. In the paper, used multi-trait Bayesian ridge-regression (MtBRR) to obtain marker effects, and also stored posterior matrices on disk to later compute posterior mean variances. This was computationally expensive and different from my standard univariate REML approach. I think MtBRR and PMV are probably the least biased way to go... but...
For the sake of testing a simple integration into the in-use pipeline, I want to try univariate REML to get the marker effects, which I'll subsequently use for the cross-validation.
Revised the functions in **`package:predCrossVar`** to increase the computational efficiency. Not yet included into the actual R package but instead sourced from `code/predCrossVar.R`. Additional speed increases were achieved by extra testing to optimize balance of `OMP_NUM_THREADS` setting (multi-core BLAS) and parallel processing of the crosses-being-predicted. Improvements will benefit users predicting with REML / Bayesian-VPM, but probably worse for Bayesian-PMV.
# Set-up server computing env.
Use a a singularity image from the rocker project, as recommended by Qi Sun to get an OpenBLAS-linked R environment that packages can easily be installed on.
**This first chunk is one-time only and doesn't take long. Saves a 650Mb `*.sif` file to server's /workdir/**
```{bash, eval=F}
# copy the project data
cd /home/jj332_cas/marnin/;
cp -R implementGMSinCassava /home/$USER;
# the project directory can be in my networked folder for 2 reasons:
# 1) singularity will automatically recognize and be able to access it
# 2) My analyses not read/write intensive; don't break server rules/etiquette
# set up a working directory on the remote machine
mkdir /workdir/$USER
cd /workdir/$USER/;
# pull a singularity image and save in the file rocker.sif
# next time you use the rocker.sif file to start the container
singularity pull rocker.sif docker://rocker/tidyverse:latest;
```
For analysis, operate each R session within a singularity Linux shell within a screen shell.
```{bash, eval=F}
# 1) start a screen shell
screen;
# 2) start the singularity Linux shell inside that
singularity shell /workdir/$USER/rocker.sif;
# Project directory, so R will use as working dir.
cd /home/mw489/implementGMSinCassava/;
# 3) Start R
R
```
# Parent-wise cross-validation
Fully-tested `runParentWiseCrossVal()` and component functions are in the `code/parentWiseCrossVal.R` script.
Below, source it and use it for a full cross-validation run.
```{r install packages if needed}
# install.packages(c("RhpcBLASctl","here","rsample","sommer","psych","future.callr","furrr","lme4"))
# install.packages('future.callr')
```
```{r primary inputs}
require(tidyverse); require(magrittr);
# 5 threads per Rsession for matrix math (openblas)
RhpcBLASctl::blas_set_num_threads(5)
# SOURCE CORE FUNCTIONS
source(here::here("code","parentWiseCrossVal.R"))
source(here::here("code","predCrossVar.R"))
# PEDIGREE
ped<-read.table(here::here("output","verified_ped.txt"),
header = T, stringsAsFactors = F) %>%
rename(GID=FullSampleName,
damID=DamID,
sireID=SireID) %>%
dplyr::select(GID,sireID,damID)
# Keep only families with _at least_ 2 offspring
ped %<>%
semi_join(ped %>% count(sireID,damID) %>% filter(n>1) %>% ungroup())
# BLUPs
blups<-readRDS(file=here::here("data","blups_forCrossVal.rds")) %>%
dplyr::select(-varcomp)
# GENOMIC RELATIONSHIP MATRICES (GRMS)
grms<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
D=readRDS(file=here::here("output",
"kinship_domGenotypic_IITA_2021July5.rds")))
## using A+domGenotypic (instead of domClassic used previously)
## will achieve appropriate dom effects for predicting family mean TGV
## but resulting add effects WILL NOT represent allele sub. effects and thus
## predictions won't equal GEBV, allele sub. effects will be post-computed
## as alpha = a + d(q-p)
# DOSAGE MATRIX
dosages<-readRDS(file=here::here("data",
"dosages_IITA_filtered_2021May13.rds"))
# RECOMBINATION FREQUENCY MATRIX
recombFreqMat<-readRDS(file=here::here("data",
"recombFreqMat_1minus2c_2021May13.rds"))
# HAPLOTYPE MATRIX
## keep only haplos for parents-in-the-pedigree
## those which will be used in prediction, saves memory
haploMat<-readRDS(file=here::here("data","haps_IITA_filtered_2021May13.rds"))
parents<-union(ped$sireID,ped$damID)
parenthaps<-sort(c(paste0(parents,"_HapA"),
paste0(parents,"_HapB")))
haploMat<-haploMat[parenthaps,colnames(recombFreqMat)]
# SELECTION INDEX WEIGHTS
## from IYR+IK
## note that not ALL predicted traits are on index
SIwts<-c(logFYLD=20,
HI=10,
DM=15,
MCMDS=-10,
logRTNO=12,
logDYLD=20,
logTOPYLD=15,
PLTHT=10)
```
## model=DirDom
Server 1: modelType="DirDom"
cbsulm17 - 112 cores, 512 GB RAM
```{r model DirDom - 5 full reps - parent-wise CV}
cvDirDom_5rep5fold<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=84,
modelType="DirDom",
ncores=20,nBLASthreads=5,
outName="output/cvDirDom_5rep5fold",
ped=ped,
blups=blups,
dosages=dosages,
haploMat=haploMat,
grms=grms,
recombFreqMat = recombFreqMat,
selInd = TRUE, SIwts = SIwts)
saveRDS(cvDirDom_5rep5fold,here::here("output","cvDirDom_5rep5fold_predAccuracy.rds"))
# [1] "Marker-effects Computed. Took 2.3851 hrs"
# [1] "Predicting cross variances and covariances"
# Joining, by = c("Repeat", "Fold")
# [1] "Done predicting fam vars. Took 59.08 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 18.63 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 64.82 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 20.41 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 46.42 mins for 156 crosses"
# [1] "Done predicting fam vars. Took 14.94 mins for 156 crosses"
# [1] "Done predicting fam vars. Took 63.45 mins for 210 crosses"
# [1] "Done predicting fam vars. Took 19.8 mins for 210 crosses"
# [1] "Done predicting fam vars. Took 50.62 mins for 171 crosses"
# [1] "Done predicting fam vars. Took 16.26 mins for 171 crosses"
# [1] "Done predicting fam vars. Took 49.87 mins for 163 crosses"
# [1] "Done predicting fam vars. Took 16.2 mins for 163 crosses"
# [1] "Done predicting fam vars. Took 73.37 mins for 253 crosses"
# [1] "Done predicting fam vars. Took 23.59 mins for 253 crosses"
# [1] "Done predicting fam vars. Took 56.32 mins for 190 crosses"
# [1] "Done predicting fam vars. Took 18.44 mins for 190 crosses"
# [1] "Done predicting fam vars. Took 47.33 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 15.79 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 59.18 mins for 189 crosses"
# [1] "Done predicting fam vars. Took 18.67 mins for 189 crosses"
# [1] "Done predicting fam vars. Took 64.72 mins for 205 crosses"
# [1] "Done predicting fam vars. Took 21.17 mins for 205 crosses"
# [1] "Done predicting fam vars. Took 63.97 mins for 213 crosses"
# [1] "Done predicting fam vars. Took 20.04 mins for 213 crosses"
# [1] "Done predicting fam vars. Took 53.03 mins for 180 crosses"
# [1] "Done predicting fam vars. Took 17.28 mins for 180 crosses"
# [1] "Done predicting fam vars. Took 58.67 mins for 199 crosses"
# [1] "Done predicting fam vars. Took 19.03 mins for 199 crosses"
# ....
# estimate 20 more hours, complete on July 12 very early AM?
# [1] "Accuracies predicted. Took 34.37369 hrs total.Goodbye!"
# Warning message:
# In for (ii in 1L:length(res)) { : closing unused connection 3 (localhost)
# > saveRDS(cvDirDom_5rep5fold,here::here("output","cvDirDom_5rep5fold_predAccuracy.rds"))
```
## model=AD
Server 2: modelType="AD"
cbsulm29 - 104 cores, 512 GB RAM
```{r model AD - 5 full reps - parent-wise CV}
grmsAD<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
D=readRDS(file=here::here("output",
"kinship_D_IITA_2021May13.rds")))
rm(grms)
cvAD_5rep5fold<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=84,
modelType="AD",
ncores=20,
outName="output/cvAD_5rep5fold",
ped=ped,
blups=blups,
dosages=dosages,
haploMat=haploMat,
grms=grmsAD,
recombFreqMat = recombFreqMat,
selInd = TRUE, SIwts = SIwts)
saveRDS(cvAD_5rep5fold,here::here("output","cvAD_5rep5fold_predAccuracy.rds"))
# [1] "Marker-effects Computed. Took 1.81086 hrs"
# [1] "Done predicting fam vars. Took 43.11 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 47.04 mins for 216 crosses"
# .....
# [1] "Accuracies predicted. Took 19.68694 hrs total.\n Goodbye!"
# [1] "Accuracies predicted. Took 19.73242 hrs total.Goodbye!"
# > saveRDS(cvAD_5rep5fold,here::here("output","cvAD_5rep5fold_predAccuracy.rds"))
```
# [TO DO] Standard clone-wise cross-validation
## Add DirDom to runCrossVal func
```{r runCrossVal}
#' @param byGroup logical, if TRUE, assumes a column named "Group" is present which unique classifies each GID into some genetic grouping.
#' @param modelType string, A, AD or ADE representing model with Additive-only, Add. plus Dominance, and Add. plus Dom. plus. AxD Epistasis (AD), respectively.
#' @param grms list of GRMs where each element is named either A, D, or, AD. Matrices supplied must match required by A, AD and ADE models. For ADE grms=list(A=A,D=D,AD=AD)...
#' @param augmentTP option to supply an additional set of training data, which will be added to each training model but never included in the test set.
#' @param TrainTestData data.frame with de-regressed BLUPs, BLUPs and weights (WT) for training and test. If byGroup==TRUE, a column with Group as the header uniquely classifying GIDs into genetic groups, is expected.
runCrossVal<-function(TrainTestData,modelType,grms,nrepeats,nfolds,ncores=1,
byGroup=FALSE,augmentTP=NULL,gid="GID",...){
require(sommer); require(rsample)
# Set-up replicated cross-validation folds
# splitting by clone (if clone in training dataset, it can't be in testing)
if(byGroup){
cvsamples<-tibble(GroupName=unique(TrainTestData$Group))
} else { cvsamples<-tibble(GroupName="None") }
cvsamples<-cvsamples %>%
mutate(Splits=map(GroupName,function(GroupName){
if(GroupName!="None"){
thisgroup<-TrainTestData %>%
filter(Group==GroupName) } else { thisgroup<-TrainTestData }
out<-tibble(repeats=1:nrepeats,
splits=rerun(nrepeats,group_vfold_cv(thisgroup, group = gid, v = nfolds))) %>%
unnest(splits)
return(out)
})) %>%
unnest(Splits)
## Internal function
## fits prediction model and calcs. accuracy for each train-test split
fitModel<-possibly(function(splits,modelType,augmentTP,TrainTestData,GroupName,grms){
starttime<-proc.time()[3]
# Set-up training set
trainingdata<-training(splits)
## Make sure, if there is an augmentTP, no GIDs in test-sets
if(!is.null(augmentTP)){
## remove any test-set members from augment TP before adding to training data
training_augment<-augmentTP %>% filter(!(!!sym(gid) %in% testing(splits)[[gid]]))
trainingdata<-bind_rows(trainingdata,training_augment) }
if(GroupName!="None"){ trainingdata<-bind_rows(trainingdata,
TrainTestData %>%
filter(Group!=GroupName,
!(!!sym(gid) %in% testing(splits)[[gid]]))) }
# Subset kinship matrices
traintestgids<-union(trainingdata[[gid]],testing(splits)[[gid]])
A1<-grms[["A"]][traintestgids,traintestgids]
trainingdata[[paste0(gid,"a")]]<-factor(trainingdata[[gid]],levels=rownames(A1))
if(modelType %in% c("AD","ADE")){
D1<-grms[["D"]][traintestgids,traintestgids]
trainingdata[[paste0(gid,"d")]]<-factor(trainingdata[[gid]],levels=rownames(D1))
if(modelType=="ADE"){
#AA1<-grms[["AA"]][traintestgids,traintestgids]
AD1<-grms[["AD"]][traintestgids,traintestgids]
diag(AD1)<-diag(AD1)+1e-06
#DD1<-grms[["DD"]][traintestgids,traintestgids]
#trainingdata[[paste0(gid,"aa")]]<-factor(trainingdata[[gid]],levels=rownames(AA1))
trainingdata[[paste0(gid,"ad")]]<-factor(trainingdata[[gid]],levels=rownames(AD1))
#trainingdata[[paste0(gid,"dd")]]<-factor(trainingdata[[gid]],levels=rownames(DD1))
}
}
# Set-up random model statements
randFormula<-paste0("~vs(",gid,"a,Gu=A1)")
if(modelType %in% c("AD","ADE")){
randFormula<-paste0(randFormula,"+vs(",gid,"d,Gu=D1)")
if(modelType=="ADE"){
randFormula<-paste0(randFormula,"+vs(",gid,"ad,Gu=AD1)")
#"+vs(",gid,"aa,Gu=AA1)",
#"+vs(",gid,"ad,Gu=AD1)")
#"+vs(",gid,"dd,Gu=DD1)")
}
}
# Fit genomic prediction model
fit <- mmer(fixed = drgBLUP ~1,
random = as.formula(randFormula),
weights = WT,
data=trainingdata)
# Gather the BLUPs
gblups<-tibble(GID=as.character(names(fit$U[[paste0("u:",gid,"a")]]$drgBLUP)),
GEBV=as.numeric(fit$U[[paste0("u:",gid,"a")]]$drgBLUP))
if(modelType %in% c("AD","ADE")){
gblups %<>% mutate(GEDD=as.numeric(fit$U[[paste0("u:",gid,"d")]]$drgBLUP))
if(modelType=="ADE"){
gblups %<>% mutate(#GEEDaa=as.numeric(fit$U[[paste0("u:",gid,"aa")]]$drgBLUP),
GEEDad=as.numeric(fit$U[[paste0("u:",gid,"ad")]]$drgBLUP))
#GEEDdd=as.numeric(fit$U[[paste0("u:",gid,"dd")]]$drgBLUP))
}
}
# Calc GETGVs
## Note that for modelType=="A", GEBV==GETGV
gblups %<>%
mutate(GETGV=rowSums(.[,grepl("GE",colnames(.))]))
# Test set validation data
validationData<-TrainTestData %>%
dplyr::select(gid,BLUP) %>%
filter(GID %in% testing(splits)[[gid]])
# Measure accuracy in test set
## cor(GEBV,BLUP)
## cor(GETGV,BLUP)
accuracy<-gblups %>%
mutate(GETGV=rowSums(.[,grepl("GE",colnames(.))])) %>%
filter(GID %in% testing(splits)[[gid]]) %>%
left_join(validationData) %>%
summarize(accGEBV=cor(GEBV,BLUP, use = 'complete.obs'),
accGETGV=cor(GETGV,BLUP, use = 'complete.obs'))
computeTime<-proc.time()[3]-starttime
accuracy %<>% mutate(computeTime=computeTime)
return(accuracy)
},otherwise = NA)
## Run models across all train-test splits
## Parallelize
require(furrr); plan(multicore, workers = ncores)
options(future.globals.maxSize=+Inf); options(future.rng.onMisuse="ignore")
cvsamples<-cvsamples %>%
mutate(accuracy=future_map2(splits,GroupName,
~fitModel(splits=.x,GroupName=.y,
modelType=modelType,augmentTP=NULL,TrainTestData=TrainTestData,grms=grms),
.progress = FALSE)) %>%
unnest(accuracy)
return(cvsamples)
}
```
# Next step / Results
6. [Genomic predictions](06-GenomicPredictions.html):
- A. Standard genomic prediction of individual GEBV and GETGV for all selection candidates using all available data.
- B. Predict cross means and variances for genomic mate selection
See [Results](07-Results.html): Home for plots and summary tables.