Last updated: 2021-08-19
Checks: 7 0
Knit directory: IITA_2021GS/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20210504) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 182b83c. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: code/.DS_Store
Ignored: data/.DS_Store
Untracked files:
Untracked: data/DatabaseDownload_2021Aug08/
Untracked: data/DatabaseDownload_2021May04/
Untracked: data/GBSdataMasterList_31818.csv
Untracked: data/IITA_GBStoPhenoMaster_33018.csv
Untracked: data/NRCRI_GBStoPhenoMaster_40318.csv
Untracked: data/PedigreeGeneticGainCycleTime_aafolabi_01122020.xls
Untracked: data/Report-DCas21-6038/
Untracked: data/blups_forGP.rds
Untracked: data/chr1_RefPanelAndGSprogeny_ReadyForGP_72719.fam
Untracked: data/dosages_IITA_2021Aug09.rds
Untracked: data/haps_IITA_2021Aug09.rds
Untracked: data/recombFreqMat_1minus2c_2021Aug02.qs
Untracked: output/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/04-PreprocessDataFiles.Rmd) and HTML (docs/04-PreprocessDataFiles.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 1c03315 | wolfemd | 2021-08-11 | Build site. |
| Rmd | e4df79f | wolfemd | 2021-08-11 | Completed IITA_2021GS pipeline including imputation and genomic prediction. Last bit of cross-validation and cross-prediction finishes in 24 hrs. |
| html | 934141c | wolfemd | 2021-07-14 | Build site. |
| Rmd | 04a2ca8 | wolfemd | 2021-07-05 | Add genotypic dominance matrix. |
| html | e66bdad | wolfemd | 2021-06-10 | Build site. |
| Rmd | a8452ba | wolfemd | 2021-06-10 | Initial build of the entire page upon completion of all |
Impute DCas21-6038: with West Africa reference panel merged with additional GS progeny (IITA TMS18)
Validate the pedigree obtained from cassavabase: Before setting up a cross-validation scheme for predictions that depend on a correct pedigree, add a basic verification step to the pipeline. Not trying to fill unknown relationships or otherwise correct the pedigree. Assess evidence that relationship is correct, remove if incorrect.
Extract haps from VCF with bcftools
library(tidyverse); library(magrittr)
pathIn<-"~/IITA_2021GS/output/"
pathOut<-pathIn
vcfName<-"AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08"
system(paste0("bcftools convert --hapsample ",
pathOut,vcfName," ",
pathIn,vcfName,".vcf.gz "))Read haps to R
library(data.table)
haps<-fread(paste0(pathIn,vcfName,".hap.gz"),
stringsAsFactors = F,header = F) %>%
as.data.frame
sampleids<-fread(paste0(pathIn,vcfName,".sample"),
stringsAsFactors = F,header = F,skip = 2) %>%
as.data.frameExtract needed GIDs from BLUPs and pedigree: Subset to: (1) genotyped-plus-phenotyped and/or (2) in verified pedigree.
blups<-readRDS(file=here::here("output",
"IITA_blupsForModelTraining_twostage_asreml_2021Aug09.rds"))
blups %>%
select(Trait,blups) %>%
unnest(blups) %>%
distinct(GID) %$% GID -> gidWithBLUPs
genotypedWithBLUPs<-gidWithBLUPs[gidWithBLUPs %in% sampleids$V1]
length(genotypedWithBLUPs)
# [1] 8669
ped<-read.table(here::here("output","verified_ped.txt"),
header = T, stringsAsFactors = F)
pednames<-union(ped$FullSampleName,
union(ped$SireID,ped$DamID))
length(pednames)
# [1] 4615
samples2keep<-union(genotypedWithBLUPs,pednames)
length(samples2keep)
# [1] 8857
# write a sample list to disk for downstream purposes
# format suitable for subsetting with --keep in plink
write.table(tibble(FID=0,IID=samples2keep),
file=here::here("output","samples2keep_IITA_2021Aug09.txt"),
row.names = F, col.names = F, quote = F)Add sample ID’s
hapids<-sampleids %>%
select(V1,V2) %>%
mutate(SampleIndex=1:nrow(.)) %>%
rename(HapA=V1,HapB=V2) %>%
pivot_longer(cols=c(HapA,HapB),
names_to = "Haplo",values_to = "SampleID") %>%
mutate(HapID=paste0(SampleID,"_",Haplo)) %>%
arrange(SampleIndex)
colnames(haps)<-c("Chr","HAP_ID","Pos","REF","ALT",hapids$HapID)Subset haps
hapids2keep<-hapids %>% filter(SampleID %in% samples2keep)
dim(haps) # [1] 61239 46669
haps<-haps[,c("Chr","HAP_ID","Pos","REF","ALT",hapids2keep$HapID)]
dim(haps) # [1] 61239 17719Format, transpose, convert to matrix and save!
haps %<>%
mutate(HAP_ID=gsub(":","_",HAP_ID)) %>%
column_to_rownames(var = "HAP_ID") %>%
select(-Chr,-Pos,-REF,-ALT)
haps %<>% t(.) %>% as.matrix(.)
saveRDS(haps,file=here::here("data","haps_IITA_2021Aug09.rds"))To ensure consistency in allele counting, create dosage from haps manually.
dosages<-haps %>%
as.data.frame(.) %>%
rownames_to_column(var = "GID") %>%
separate(GID,c("SampleID","Haplo"),"_Hap",remove = T) %>%
select(-Haplo) %>%
group_by(SampleID) %>%
summarise(across(everything(),~sum(.))) %>%
ungroup() %>%
column_to_rownames(var = "SampleID") %>%
as.matrix
dim(dosages)
# [1] 8857 61239
saveRDS(dosages,file=here::here("data","dosages_IITA_2021Aug09.rds"))Apply two variant filters: According to the results of preliminary experiments I will work with two marker sets downstream, both filtered to MAF>=1%.
plink --indep-pairwise 50 25 0.98 (considered the “full_set”)plink --indep-pairwise 1000 'kb' 50 0.6 (considered the “reduced_set”)cd ~/IITA_2021GS/
export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;
# Lightly pruned "full_set"
plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
--keep output/samples2keep_IITA_2021Aug09.txt \
--maf 0.01 \
--indep-pairwise 50 25 0.98 \
--out output/samples2keep_IITA_MAFpt01_prune50_25_pt98;
# 60896 variants and 8857 people pass filters and QC.
# Note: No phenotypes present.
# Pruned 4504 variants from chromosome 1, leaving 3064.
# Pruned 1552 variants from chromosome 2, leaving 2033.
# Pruned 1578 variants from chromosome 3, leaving 2049.
# Pruned 1931 variants from chromosome 4, leaving 1438.
# Pruned 1805 variants from chromosome 5, leaving 1865.
# Pruned 1600 variants from chromosome 6, leaving 1739.
# Pruned 725 variants from chromosome 7, leaving 970.
# Pruned 1440 variants from chromosome 8, leaving 1696.
# Pruned 1599 variants from chromosome 9, leaving 1621.
# Pruned 1082 variants from chromosome 10, leaving 1521.
# Pruned 1261 variants from chromosome 11, leaving 1622.
# Pruned 1253 variants from chromosome 12, leaving 1477.
# Pruned 1257 variants from chromosome 13, leaving 1353.
# Pruned 2456 variants from chromosome 14, leaving 2716.
# Pruned 1733 variants from chromosome 15, leaving 1772.
# Pruned 1390 variants from chromosome 16, leaving 1329.
# Pruned 1108 variants from chromosome 17, leaving 1485.
# Pruned 1306 variants from chromosome 18, leaving 1566.
# Pruning complete. 29580 of 60896 variants removed.
# Stronger "reduced_set"
plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
--keep output/samples2keep_IITA_2021Aug09.txt \
--maf 0.01 \
--indep-pairwise 1000 'kb' 50 0.6 \
--out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt6;
# 60896 variants and 8857 people pass filters and QC.
# Note: No phenotypes present.
# Pruned 7035 variants from chromosome 1, leaving 533.
# Pruned 2949 variants from chromosome 2, leaving 636.
# Pruned 3033 variants from chromosome 3, leaving 594.
# Pruned 2931 variants from chromosome 4, leaving 438.
# Pruned 3172 variants from chromosome 5, leaving 498.
# Pruned 2832 variants from chromosome 6, leaving 507.
# Pruned 1336 variants from chromosome 7, leaving 359.
# Pruned 2635 variants from chromosome 8, leaving 501.
# Pruned 2842 variants from chromosome 9, leaving 378.
# Pruned 2148 variants from chromosome 10, leaving 455.
# Pruned 2455 variants from chromosome 11, leaving 428.
# Pruned 2292 variants from chromosome 12, leaving 438.
# Pruned 2210 variants from chromosome 13, leaving 400.
# Pruned 4583 variants from chromosome 14, leaving 589.
# Pruned 3137 variants from chromosome 15, leaving 368.
# Pruned 2266 variants from chromosome 16, leaving 453.
# Pruned 2251 variants from chromosome 17, leaving 342.
# Pruned 2442 variants from chromosome 18, leaving 430.
# Pruning complete. 52549 of 60896 variants removed.After initial test results, the “reduced_set” was too harsh. Try a “medium_set”:
cd ~/IITA_2021GS/
export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;
plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
--keep output/samples2keep_IITA_2021Aug09.txt \
--maf 0.01 \
--indep-pairwise 1000 'kb' 50 0.8 \
--out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt8;
# Before main variant filters, 8857 founders and 0 nonfounders present.
# Calculating allele frequencies... done.
# 343 variants removed due to minor allele threshold(s)
# (--maf/--max-maf/--mac/--max-mac).
# 60896 variants and 8857 people pass filters and QC.
# Note: No phenotypes present.
# Pruned 6645 variants from chromosome 1, leaving 923.
# Pruned 2568 variants from chromosome 2, leaving 1017.
# Pruned 2722 variants from chromosome 3, leaving 905.
# Pruned 2751 variants from chromosome 4, leaving 618.
# Pruned 2862 variants from chromosome 5, leaving 808.
# Pruned 2515 variants from chromosome 6, leaving 824.
# Pruned 1168 variants from chromosome 7, leaving 527.
# Pruned 2336 variants from chromosome 8, leaving 800.
# Pruned 2618 variants from chromosome 9, leaving 602.
# Pruned 1862 variants from chromosome 10, leaving 741.
# Pruned 2156 variants from chromosome 11, leaving 727.
# Pruned 2056 variants from chromosome 12, leaving 674.
# Pruned 2010 variants from chromosome 13, leaving 600.
# Pruned 4123 variants from chromosome 14, leaving 1049.
# Pruned 2856 variants from chromosome 15, leaving 649.
# Pruned 2088 variants from chromosome 16, leaving 631.
# Pruned 2029 variants from chromosome 17, leaving 564.
# Pruned 2201 variants from chromosome 18, leaving 671.
# Pruning complete. 47566 of 60896 variants removed.
# Marker lists written to
# output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt8.prune.in and
# output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt8.prune.outUsed plink to output a list of pruned SNPs.
# 1) start a screen shell
screen; # or screen -r if re-attaching...
# 2) start the singularity Linux shell inside that
#singularity shell /workdir/$USER/rocker.sif;
singularity shell ~/rocker2.sif;
# Project directory, so R will use as working dir.
cd /home/mw489/IITA_2021GS/;
# 3) Start R
Rlibrary(tidyverse); library(magrittr); library(genomicMateSelectR)
RhpcBLASctl::blas_set_num_threads(56)
dosages<-readRDS(file=here::here("data","dosages_IITA_2021Aug09.rds"))
full_set<-read.table(here::here("output",
"samples2keep_IITA_MAFpt01_prune50_25_pt98.prune.in"),
header = F, stringsAsFactors = F)
medium_set<-read.table(here::here("output",
"samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt8.prune.in"),
header = F, stringsAsFactors = F)
reduced_set<-read.table(here::here("output",
"samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt6.prune.in"),
header = F, stringsAsFactors = F)
full_set<-tibble(FULL_SNP_ID=colnames(dosages)) %>%
separate(FULL_SNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>%
mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>%
filter(SNP_ID %in% full_set$V1)
medium_set<-tibble(FULL_SNP_ID=colnames(dosages)) %>%
separate(FULL_SNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>%
mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>%
filter(SNP_ID %in% medium_set$V1)
reduced_set<-tibble(FULL_SNP_ID=colnames(dosages)) %>%
separate(FULL_SNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>%
mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>%
filter(SNP_ID %in% reduced_set$V1)
# > dim(reduced_set)
# [1] 8347 6
# > dim(medium_set)
# [1] 13330 6
# > dim(full_set)
# [1] 31316 6
snpsets<-full_set %>%
mutate(Set="full_set") %>%
bind_rows(medium_set %>%
mutate(Set="medium_set")) %>%
bind_rows(reduced_set %>%
mutate(Set="reduced_set")) %>%
nest(snps2keep=-Set)
saveRDS(snpsets,file = here::here("data","snpsets.rds"))# Kinships from full snp set
A<-kinship(dosages[,full_set$FULL_SNP_ID],type="add")
D<-kinship(dosages[,full_set$FULL_SNP_ID],type="domGenotypic")
saveRDS(A,file=here::here("output","kinship_A_IITA_2021Aug09.rds"))
saveRDS(D,file=here::here("output","kinship_Dgeno_IITA_2021Aug09.rds"))
# kinships from reduced snp set
A<-kinship(dosages[,medium_set$FULL_SNP_ID],type="add")
D<-kinship(dosages[,medium_set$FULL_SNP_ID],type="domGenotypic")
saveRDS(A,file=here::here("output","kinship_A_MediumSNPset_IITA_2021Aug09.rds"))
saveRDS(D,file=here::here("output","kinship_Dgeno_MediumSNPset_IITA_2021Aug09.rds"))
# kinships from reduced snp set
A<-kinship(dosages[,reduced_set$FULL_SNP_ID],type="add")
D<-kinship(dosages[,reduced_set$FULL_SNP_ID],type="domGenotypic")
saveRDS(A,file=here::here("output","kinship_A_ReducedSNPset_IITA_2021Aug09.rds"))
saveRDS(D,file=here::here("output","kinship_Dgeno_ReducedSNPset_IITA_2021Aug09.rds"))The genetic map and recombination frequency matrix including all the SNPs that will be in any downstream analysis, were created previously. See here.
Download recombFreqMat from Cassavabase FTP server: Click Here* to download recombFreqMat_1minus2c_2021Aug02.qs.
Copy to project directory.
cp ~/implementGMSinCassava/data/recombFreqMat_1minus2c_2021Aug02.qs ~/IITA_2021GS/data/ped<-read.table(here::here("output","verified_ped.txt"),
header = T, stringsAsFactors = F)Select traits and data to be analyzed.
# This list from Dec. 2020 GeneticGain rate estimation
# These were what Ismail/IITA/BMGF wanted to see
# Will cross-validate these traits
traits<-c("logDYLD","logFYLD","logRTNO","logTOPYLD","MCMDS","DM","BCHROMO",
"PLTHT","BRLVLS","BRNHT1","HI")
# only 8 of these on the SELIND
# Full trait list = 14:
## traits<-c("MCMDS","DM","PLTHT","BRNHT1","BRLVLS","HI",
## "logDYLD", # <-- logDYLD now included.
## "logFYLD","logTOPYLD","logRTNO","TCHART","LCHROMO","ACHROMO","BCHROMO")library(tidyverse); library(magrittr);
blups<-readRDS(file=here::here("output",
"IITA_blupsForModelTraining_twostage_asreml_2021Aug09.rds"))
dosages<-readRDS(file=here::here("data","dosages_IITA_2021Aug09.rds"))
blups %>%
select(Trait,blups) %>%
unnest(blups) %>%
distinct(GID) %$% GID -> gidWithBLUPs
# [1] 77854
length(gidWithBLUPs)
genotypedWithBLUPs<-gidWithBLUPs[gidWithBLUPs %in% rownames(dosages)]
length(genotypedWithBLUPs)
# [1] 8669
blups %<>%
filter(Trait %in% traits) %>%
select(Trait,blups,varcomp) %>%
mutate(blups=map(blups,~filter(.,GID %in% genotypedWithBLUPs)))
saveRDS(blups,file=here::here("data","blups_forGP.rds"))# 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)
SIwts logFYLD HI DM MCMDS logRTNO logDYLD logTOPYLD PLTHT
20 10 15 -10 12 20 15 10 Soon, will get a precise list from IITA team. For now: went to cassavabase wizard, made a list of all IITA accessions in field trials dated 2020 and 2021 at 4 locations (IBA, IKN, MOK, UBJ), as these can likely source for crossing nurseries. Cross reference that list with my genos-2-phenos matches to get a list of clones I can make genomic predictions for.
library(tidyverse); library(magrittr);
dbdata<-readRDS(here::here("output","IITA_ExptDesignsDetected_2021Aug08.rds"))
accessions_possibly_infield<-read.table(here::here("data","Accessions_IITA_LikelyInField_IbaIknMokUbj_2020to2021.txt"),
header = F, stringsAsFactors = F) %>%
as_tibble %>%
rename(germplasmName=V1) %>%
inner_join(dbdata %>%
distinct(germplasmName,GID,FullSampleName)) %>%
filter(!is.na(FullSampleName)) %>%
mutate(Cohort=NA,
Cohort=ifelse(grepl("TMS20",germplasmName,ignore.case = T),"TMS20",
ifelse(grepl("TMS19",germplasmName,ignore.case = T),"TMS19",
ifelse(grepl("TMS18",germplasmName,ignore.case = T),"TMS18",
ifelse(grepl("TMS17",germplasmName,ignore.case = T),"TMS17",
ifelse(grepl("TMS16",germplasmName,ignore.case = T),"TMS16",
ifelse(grepl("TMS15",germplasmName,ignore.case = T),"TMS15",
ifelse(grepl("TMS14",germplasmName,ignore.case = T),"TMS14",
ifelse(grepl("TMS13|2013_",germplasmName,ignore.case = T),"TMS13","GGetc")))))))))
saveRDS(accessions_possibly_infield,
here::here("data","accessions_possibly_infield_2021Aug10.rds"))
# accessions_possibly_infield %>%
# count(Cohort)
# Cohort n
# GGetc 1107
# TMS13 85
# TMS14 74
# TMS15 19
# TMS16 1
# TMS17 38
# TMS18 671
# TMS19 186
# TMS20 674 # 1) start a screen shell
screen; # or screen -r if re-attaching...
# 2) Pull a singularity impage from rocker/
singularity pull ~/rocker2.sif docker://rocker/tidyverse:latest;
# 3) start the singularity Linux shell inside that
singularity shell ~/rocker2.sif;
# Project directory, so R will use as working dir.
cd /home/mw489/IITA_2021GS/;
# 3) Start R
R# Install genomicMateSelectR to user-accessible libPath
libPath<-"/home/mw489/R/x86_64-pc-linux-gnu-library/4.1"
withr::with_libpaths(new=libPath, devtools::install_github("wolfemd/genomicMateSelectR", ref = 'master'))Parent-wise and standard cross-validation: estimate selection index (and component trait) prediction accuracies using the direction-dominance (DirDom) model.
Genomic predictions: First, predict of individual GEBV/GETGV for all selection candidates using all available data and return marker effects for use downstream. Next, Select a top set of candidate parents, for whom we would like to predict cross performances. Finally, predict all pairwise crosses of candidate parents and evaluate them for genomic mate selection. Select the top crosses and plant a crossing nursery with the parents indicated.
sessionInfo()R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] workflowr_1.6.2
loaded via a namespace (and not attached):
[1] Rcpp_1.0.7 whisker_0.4 knitr_1.33 magrittr_2.0.1
[5] R6_2.5.0 rlang_0.4.11 fansi_0.5.0 stringr_1.4.0
[9] tools_4.1.0 xfun_0.25 utf8_1.2.2 git2r_0.28.0
[13] jquerylib_0.1.4 htmltools_0.5.1.1 ellipsis_0.3.2 rprojroot_2.0.2
[17] yaml_2.2.1 digest_0.6.27 tibble_3.1.3 lifecycle_1.0.0
[21] crayon_1.4.1 later_1.2.0 sass_0.4.0 vctrs_0.3.8
[25] promises_1.2.0.1 fs_1.5.0 glue_1.4.2 evaluate_0.14
[29] rmarkdown_2.10 stringi_1.7.3 bslib_0.2.5.1 compiler_4.1.0
[33] pillar_1.6.2 jsonlite_1.7.2 httpuv_1.6.1 pkgconfig_2.0.3