06-Results.Rmd

---
title: "Results"
site: workflowr::wflow_site
date: "2021-June-07"
output: 
  workflowr::wflow_html:
    code_folding: hide
editor_options:
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = F, 
                      tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```

# Raw data

Summary of the number of unique plots, locations, years, etc. in the cleaned plot-basis data. See [here](01-cleanTPdata.html) for details..

```{r}
library(tidyverse); library(magrittr); library(ragg)
rawdata<-readRDS(file=here::here("output","IITA_ExptDesignsDetected_2021May10.rds"))
rawdata %>% 
  summarise(Nplots=nrow(.),
            across(c(locationName,studyYear,studyName,TrialType,GID), ~length(unique(.)),.names = "N_{.col}")) %>% 
  rmarkdown::paged_table()
```
This _is not_ the same number of clones as are expected to be genotyped-and-phenotyped.

Break down the plots based on the trial design and TrialType (really a grouping of the population that is breeding program specific), captured by two logical variables, **CompleteBlocks** and **IncompleteBlocks**.  
```{r}
rawdata %>% 
  count(TrialType,CompleteBlocks,IncompleteBlocks) %>% 
  spread(TrialType,n) %>% 
  rmarkdown::paged_table()
```

Next, look at breakdown of plots by TrialType (rows) and locations (columns):
```{r, cols.print=16}
rawdata %>% 
  count(locationName,TrialType) %>% 
  spread(locationName,n) %>% 
  rmarkdown::paged_table()
```

```{r}
traits<-c("MCMDS","DM","PLTHT","BRNHT1","BRLVLS","HI",
          "logDYLD", # <-- logDYLD now included. 
          "logFYLD","logTOPYLD","logRTNO","TCHART","LCHROMO","ACHROMO","BCHROMO")
rawdata %>% 
  select(locationName,studyYear,studyName,TrialType,any_of(traits)) %>% 
  pivot_longer(cols = any_of(traits), values_to = "Value", names_to = "Trait") %>% 
  ggplot(.,aes(x=Value,fill=Trait)) + geom_histogram() + facet_wrap(~Trait, scales='free') + 
  theme_bw() + scale_fill_viridis_d() + 
  labs(title = "Distribution of Raw Phenotypic Values")
```
How many genotyped-and-phenotyped clones?

```{r}
rawdata %>% 
  select(locationName,studyYear,studyName,TrialType,germplasmName,FullSampleName,GID,any_of(traits)) %>% 
  pivot_longer(cols = any_of(traits), values_to = "Value", names_to = "Trait") %>%
  filter(!is.na(Value),!is.na(FullSampleName)) %>%
  distinct(germplasmName,FullSampleName,GID) %>% 
  rmarkdown::paged_table()
```
There are 8149 genotyped-and-phenotyped clones! 

# BLUPs

These are the BLUPs combining data for each clone across trials/locations _without_ genomic information, used as input for genomic prediction downstream.

```{r, message=F, rows.print=12}
blups<-readRDS(file=here::here("data","blups_forCrossVal.rds")) 
gidWithBLUPs<-blups %>% select(Trait,blups) %>% unnest(blups) %$% unique(GID)
rawdata %>% 
  select(observationUnitDbId,GID,any_of(blups$Trait)) %>% 
  pivot_longer(cols = any_of(blups$Trait), 
               names_to = "Trait", 
               values_to = "Value",values_drop_na = T) %>% 
  filter(GID %in% gidWithBLUPs) %>% 
  group_by(Trait) %>% 
  summarize(Nplots=n()) %>% 
  ungroup() %>% 
  left_join(blups %>% 
              mutate(Nclones=map_dbl(blups,~nrow(.)),
                     avgREL=map_dbl(blups,~mean(.$REL)),
                     Vg=map_dbl(varcomp,~.["GID!GID.var","component"]),
                     Ve=map_dbl(varcomp,~.["R!variance","component"]),
                     H2=Vg/(Vg+Ve)) %>% 
              select(-blups,-varcomp)) %>% 
  mutate(across(is.numeric,~round(.,3))) %>% arrange(desc(H2)) %>% 
  rmarkdown::paged_table()
```

```{r}
blups %>% 
  select(Trait,blups) %>% 
  unnest(blups) %>% 
  ggplot(.,aes(x=drgBLUP,fill=Trait)) + geom_histogram() + facet_wrap(~Trait, scales='free') + 
  theme_bw() + scale_fill_viridis_d() + theme(legend.position = 'none') + 
  labs(title = "Distribution of de-regressed BLUP Values")
```

```{r}
blups %>% 
  select(Trait,blups) %>% 
  unnest(blups) %>% 
  ggplot(.,aes(x=Trait,y=REL,fill=Trait)) + geom_boxplot(notch=T) + #facet_wrap(~Trait, scales='free') + 
  theme_bw() + scale_fill_viridis_d() + 
  theme(axis.text.x = element_text(angle=90),
        legend.position = 'none') + 
  labs(title = "Distribution of BLUP Reliabilities")
```
# Marker density and distribution
```{r}
library(tidyverse); library(magrittr);
snps<-readRDS(file=here::here("data","dosages_IITA_filtered_2021May13.rds"))
mrks<-colnames(snps) %>% 
  tibble(SNP_ID=.) %>% 
  separate(SNP_ID,c("Chr","Pos","Allele"),"_") %>% 
  mutate(Chr=as.integer(gsub("S","",Chr)),
         Pos=as.numeric(Pos))
mrks %>% 
  ggplot(.,aes(x=Pos,fill=as.character(Chr))) + geom_histogram() + 
  facet_wrap(~Chr,scales = 'free') + theme_bw() + 
  scale_fill_viridis_d() + theme(legend.position = 'none',
                                 axis.text.x = element_text(angle=90))
```
```{r, rows.print=18}
mrks %>% count(Chr) %>% rmarkdown::paged_table()
```

# Pedigree

Summarize the pedigree and the [verification results described here](03-validatePedigree.html).

```{r}
library(tidyverse); library(magrittr);
ped2check_genome<-readRDS(file=here::here("output","ped2check_genome.rds"))
ped2check_genome %<>% 
  select(IID1,IID2,Z0,Z1,Z2,PI_HAT)
ped2check<-read.table(file=here::here("output","ped2genos.txt"),
                      header = F, stringsAsFactors = F) %>% 
  rename(FullSampleName=V1,DamID=V2,SireID=V3)
ped2check %<>% 
  select(FullSampleName,DamID,SireID) %>% 
  inner_join(ped2check_genome %>% 
               rename(FullSampleName=IID1,DamID=IID2) %>% 
               bind_rows(ped2check_genome %>% 
                           rename(FullSampleName=IID2,DamID=IID1))) %>% 
  distinct %>% 
  mutate(FemaleParent=case_when(Z0<0.32 & Z1>0.67~"Confirm",
                                       SireID==DamID & PI_HAT>0.6 & Z0<0.3 & Z2>0.32~"Confirm",
                                       TRUE~"Reject")) %>% 
  select(-Z0,-Z1,-Z2,-PI_HAT) %>% 
  inner_join(ped2check_genome %>% 
               rename(FullSampleName=IID1,SireID=IID2) %>% 
               bind_rows(ped2check_genome %>% 
                           rename(FullSampleName=IID2,SireID=IID1))) %>% 
  distinct %>% 
  mutate(MaleParent=case_when(Z0<0.32 & Z1>0.67~"Confirm",
                                       SireID==DamID & PI_HAT>0.6 & Z0<0.3 & Z2>0.32~"Confirm",
                                       TRUE~"Reject")) %>% 
  select(-Z0,-Z1,-Z2,-PI_HAT)
rm(ped2check_genome)

ped2check %<>% 
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS18",FullSampleName,ignore.case = T),"TMS18",
                       ifelse(grepl("TMS15",FullSampleName,ignore.case = T),"TMS15",
                              ifelse(grepl("TMS14",FullSampleName,ignore.case = T),"TMS14",
                                     ifelse(grepl("TMS13|2013_",FullSampleName,ignore.case = T),"TMS13","GGetc")))))
```

Proportion of accessions with male, female or both parents in pedigree confirm-vs-rejected?
```{r}
ped2check %>% 
  count(FemaleParent,MaleParent) %>% 
  mutate(Prop=round(n/sum(n),2))
```

Proportion of accessions _within_ each Cohort with pedigree records confirmed-vs-rejected?
```{r}
ped2check %>% 
  count(Cohort,FemaleParent,MaleParent) %>% 
  spread(Cohort,n) %>% 
  mutate(across(is.numeric,~round(./sum(.),2))) %>% 
  rmarkdown::paged_table()
```

Use only fully-confirmed families / trios. 
Remove any without both parents confirmed.


```{r}
ped<-read.table(here::here("output","verified_ped.txt"),
                header = T, stringsAsFactors = F) %>% 
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS18",FullSampleName,ignore.case = T),"TMS18",
                       ifelse(grepl("TMS15",FullSampleName,ignore.case = T),"TMS15",
                              ifelse(grepl("TMS14",FullSampleName,ignore.case = T),"TMS14",
                                     ifelse(grepl("TMS13|2013_",FullSampleName,ignore.case = T),"TMS13","GGetc")))))
```
Summary of family sizes
```{r}
ped %>% 
  count(SireID,DamID) %$% summary(n)
```
```{r}
ped %>% nrow(.) # 4259 pedigree entries
ped %>% 
  count(Cohort,name = "Number of Verified Pedigree Entries")
```
```{r}
ped %>% 
  distinct(Cohort,SireID,DamID) %>% 
  count(Cohort,name = "Number of Families per Cohort")
```
730 families. Mean size 5.85, range 1-77.

# Parent-wise Cross-validation

```{r}
parentfolds<-readRDS(file=here::here("output","parentfolds.rds"))
summarized_parentfolds<-parentfolds %>% 
  mutate(Ntestparents=map_dbl(testparents,length),
         Ntrainset=map_dbl(trainset,length),
         Ntestset=map_dbl(testset,length),
         NcrossesToPredict=map_dbl(CrossesToPredict,nrow)) %>% 
  select(Repeat,Fold,starts_with("N"))
summarized_parentfolds %>% 
  rmarkdown::paged_table()
```
```{r}
summarized_parentfolds %>% summarize(across(is.numeric,median,.names = "median{.col}"))
```

Selection index weights
```{r}
c(logFYLD=20,
  HI=10,
  DM=15,
  MCMDS=-10,
  logRTNO=12,
  logDYLD=20,
  logTOPYLD=15,
  PLTHT=10)
```


## Prediction accuracy


3. [Check prediction accuracy](03-CrossValidation.html): Evaluate prediction accuracy with cross-validation.

### Selection Index Accuracy

```{r}
library(ggdist)
cvmeans<-readRDS(here::here("output","cvMeanPredAccuracyAD.rds"))
cvvars<-readRDS(here::here("output","cvVarPredAccuracyAD.rds"))
acc<-cvmeans %>% 
  filter(Trait=="SELIND") %>% 
  mutate(VarComp=gsub("Mean","",predOf),
         predOf="Mean") %>% 
  bind_rows(cvvars %>% 
              filter(Trait1=="SELIND") %>% 
              rename(Trait=Trait1) %>% 
              select(-Trait2) %>% 
              mutate(VarComp=gsub("Var","",predOf),
         predOf="Var")) 

colors<-viridis::viridis(4)[1:2]

acc %>% 
  mutate(predOf=factor(predOf,levels=c("Mean","Var")),
         VarComp=factor(VarComp,levels=c("BV","TGV"))) %>% 
  ggplot(.,aes(x=VarComp,y=AccuracyEst,fill=VarComp)) +
  ggdist::stat_halfeye(adjust=0.5,.width = 0,fill='gray',width=0.75) +  
  geom_boxplot(width=0.12,notch = TRUE) +
  ggdist::stat_dots(side = "left",justification = 1.1,
                    binwidth = 0.03,dotsize=0.6) +
  theme_bw() + 
  scale_fill_manual(values = colors) + 
  geom_hline(yintercept = 0, color='black', size=0.9) +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title = element_text(face='bold',color = 'black'),
        strip.text.x = element_text(face='bold',color='black',size=14),
        axis.text.y = element_text(face = 'bold',color='black'),
        legend.text = element_text(face='bold'),
        legend.position = 'bottom') +
  facet_grid(.~predOf) + 
  labs(title="Selection Index Prediction Accuracy") + 
  coord_cartesian(xlim = c(1.2, NA))
```


### Accuracy predicting cross-means
```{r,fig.width = 7.5, fig.height=3.5}
cvmeans<-readRDS(here::here("output","cvMeanPredAccuracyAD.rds"))
#  bind_rows(readRDS(here::here("output","cvMeanPredAccuracyA.rds")))

cvmeans %<>% 
  mutate(Trait=factor(Trait,levels=c("SELIND",blups$Trait)),
         predOf=factor(predOf,levels=c("MeanBV","MeanTGV")))

cvmeans %>% 
  ggplot(.,aes(x=predOf,y=AccuracyEst,fill=predOf,color=predOf)) + 
  geom_boxplot(width=0.4,notch = TRUE, color='gray40') +
  ggdist::stat_dots(side = "left", justification = 1.3,
                    binwidth = 0.03,dotsize=0.5,layout="swarm") +
  scale_fill_manual(values = viridis::viridis(4)[1:2]) + 
  scale_color_manual(values = viridis::viridis(4)[1:2]) + 
  geom_hline(yintercept = 0, color='black', size=0.8) +
  facet_grid(.~Trait) + 
  labs(title="Accuracy predicting cross-means") + 
  theme_bw() + 
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_text(face='bold'),
        strip.text.x = element_text(face='bold'),
        legend.position='bottom',
        legend.text = element_text(face='bold'),
        panel.spacing = unit(0.05, "lines"))
```

### Accuracy predicting cross-variance params
```{r,fig.width = 7.5, fig.height=3.5}
cvvars<-readRDS(here::here("output","cvVarPredAccuracyAD.rds"))
# cvvars<-readRDS(here::here("output","cvVarPredAccuracyA.rds")) %>% 
#   bind_rows(readRDS(here::here("output","cvVarPredAccuracyAD.rds")))
cvvars %<>% 
  mutate(Trait1=factor(Trait1,levels=c("SELIND",blups$Trait)),
         Trait2=factor(Trait2,levels=c("SELIND",blups$Trait)),
         predOf=factor(predOf,levels=c("VarBV","VarTGV")))

cvvars %>% 
  filter(Trait1==Trait2) %>% 
  ggplot(.,aes(x=predOf,y=AccuracyEst,fill=predOf, color=predOf)) + 
  geom_boxplot(width=0.4,notch = TRUE,color='gray40') +
  ggdist::stat_dots(side = "left", justification = 1.3,layout='swarm',
                    binwidth = 0.03,dotsize=0.4) +
  scale_fill_manual(values = viridis::viridis(4)[1:2]) + 
  scale_color_manual(values = viridis::viridis(4)[1:2]) + 
  geom_hline(yintercept = 0, color='black', size=0.8) +
  facet_grid(.~Trait1) + 
  labs(title="Accuracy predicting cross-variances") + 
  theme_bw() + 
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_text(face='bold'),
        strip.text.x = element_text(face='bold'),
        legend.position='bottom',
        legend.text = element_text(face='bold'),
        panel.spacing = unit(0.05, "lines"))
```


```{r, fig.width = 7.5, fig.height = 7.5}
cvvars %>% 
  filter(Trait1!=Trait2) %>% 
  ggplot(.,aes(x=predOf,y=AccuracyEst,fill=predOf)) + 
  geom_boxplot(notch = TRUE) +
  scale_fill_manual(values = viridis::viridis(4)[1:2]) + 
  scale_color_manual(values = viridis::viridis(4)[1:2]) + 
  geom_hline(yintercept = 0, color='gray40', size=0.6) +
  facet_grid(Trait1~Trait2) + 
  labs(title="Accuracy predicting cross-covariances") + 
  theme_bw() + 
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_text(face='bold'),
        strip.text.x = element_text(size=8),
        strip.text.y = element_text(size=8,angle = 0),
        legend.position='bottom',
        legend.text = element_text(face='bold'),
        panel.spacing = unit(0.05, "lines"))

```


# Genomic Mate Selection

Placeholder for summaries of the full-model predictions of cross means and variances to-be-used for selection.

# Genetic Gain Estimates