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()

rawdata %>% 
  count(TrialType,CompleteBlocks,IncompleteBlocks) %>% 
  spread(TrialType,n) %>% 
  rmarkdown::paged_table()

rawdata %>% 
  count(locationName,TrialType) %>% 
  spread(locationName,n) %>% 
  rmarkdown::paged_table()

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")

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()

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()

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")

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")

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))

mrks %>% count(Chr) %>% rmarkdown::paged_table()

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")))))

ped2check %>% 
  count(FemaleParent,MaleParent) %>% 
  mutate(Prop=round(n/sum(n),2))

ped2check %>% 
  count(Cohort,FemaleParent,MaleParent) %>% 
  spread(Cohort,n) %>% 
  mutate(across(is.numeric,~round(./sum(.),2))) %>% 
  rmarkdown::paged_table()

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")))))

ped %>% 
  count(SireID,DamID) %$% summary(n)

ped %>% nrow(.) # 4259 pedigree entries

ped %>% 
  count(Cohort,name = "Number of Verified Pedigree Entries")

ped %>% 
  distinct(Cohort,SireID,DamID) %>% 
  count(Cohort,name = "Number of Families per Cohort")

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()

summarized_parentfolds %>% summarize(across(is.numeric,median,.names = "median{.col}"))

c(logFYLD=20,
  HI=10,
  DM=15,
  MCMDS=-10,
  logRTNO=12,
  logDYLD=20,
  logTOPYLD=15,
  PLTHT=10)

# GBLUPs
### Two models AD and DirDom
#gpreds_ad<-readRDS(file = here::here("output","genomicPredictions_ModelAD.rds"))
gpreds_dirdom<-readRDS(file = here::here("output","genomicPredictions_ModelDirDom.rds"))
si_getgvs<-gpreds_dirdom$gblups[[1]] %>% 
#si_getgvs<-gpreds_ad$gblups[[1]] %>% 
  filter(predOf=="GETGV") %>% 
  select(GID,SELIND)

## IITA Germplasm Ages
ggcycletime<-readxl::read_xls(here::here("data","PedigreeGeneticGainCycleTime_aafolabi_01122020.xls")) %>% 
  mutate(Year_Accession=as.numeric(Year_Accession))

# Need germplasmName field from raw trial data to match GEBV and cycle time
rawdata<-readRDS(here::here("output","IITA_ExptDesignsDetected_2021May10.rds"))
si_getgvs %<>% 
  left_join(rawdata %>% 
              distinct(germplasmName,GID)) %>% 
  group_by(GID) %>% 
  slice(1) %>% 
  ungroup()
# table(ggcycletime$Accession %in% si_getgvs$germplasmName)
# FALSE  TRUE 
#   193   614 

si_getgvs %<>% 
  left_join(.,ggcycletime %>% 
              rename(germplasmName=Accession) %>% 
              mutate(Year_Accession=as.numeric(Year_Accession))) %>% 
  mutate(Year_Accession=case_when(grepl("2013_|TMS13",germplasmName)~2013,
                                  grepl("TMS14",germplasmName)~2014,
                                  grepl("TMS15",germplasmName)~2015,
                                  grepl("TMS18",germplasmName)~2018,
                                  !grepl("2013_|TMS13|TMS14|TMS15|TMS18",germplasmName)~Year_Accession))

# Declare the "eras" as PreGS\<2012 and GS\>=2013.
si_getgvs %<>% 
  filter(Year_Accession>2012 | Year_Accession<2005)
si_getgvs %<>% 
  mutate(GeneticGroup=ifelse(Year_Accession>=2013,"GS","PreGS"))

# , fig.height=10, fig.width=12
si_getgvs %>% 
  select(GeneticGroup,GID,Year_Accession,SELIND) %>% 
  ggplot(.,aes(x=Year_Accession,y=SELIND,color=GeneticGroup)) + 
  geom_point(size=1.25) + geom_smooth(method=lm, se=TRUE, size=1.5) + 
  # facet_wrap(~Trait,scales='free_y', ncol=2) + 
  theme_bw() +
  theme(axis.text = element_text(face = 'bold',angle = 0, size=14),
        axis.title = element_text(face = 'bold',size=16),
        strip.background.x = element_blank(),
        strip.text = element_text(face='bold',size=18)) + 
  scale_color_viridis_d() + 
  labs(title = "Selection Index GETGV vs. Accession Age by 'era' [GS vs. PreGS]",
       subtitle = "SI GETGV from modelType='DirDom'")

si_gebvs<-gpreds_dirdom$gblups[[1]] %>% 
  filter(predOf=="GEBV") %>% 
  select(GID,SELIND) %>% 
  mutate(GeneticGroup=case_when(grepl("2013_|TMS13",GID)~"C1",
                                  grepl("TMS14",GID)~"C2",
                                  grepl("TMS15",GID)~"C3",
                                  grepl("TMS18",GID)~"C4",
                                  !grepl("2013_|TMS13|TMS14|TMS15|TMS18",GID)~"PreGS"),
         GeneticGroup=factor(GeneticGroup,levels = c("PreGS","C1","C2","C3","C4")))
si_gebvs %>% 
  group_by(GeneticGroup) %>% 
  summarize(meanGEBV=mean(SELIND),
            stdErr=sd(SELIND)/sqrt(n()),
            upperSE=meanGEBV+stdErr,
            lowerSE=meanGEBV-stdErr) %>% 
  ggplot(.,aes(x=GeneticGroup,y=meanGEBV,fill=GeneticGroup)) + 
  geom_bar(stat = 'identity', color='gray60', size=1.25) + 
  geom_linerange(aes(ymax=upperSE,
                     ymin=lowerSE), color='gray60', size=1.25) + 
  #facet_wrap(~Trait,scales='free_y', ncol=1) + 
  theme_bw() +
  geom_hline(yintercept = 0, size=1.15, color='black') + 
  theme(axis.text.x = element_text(face = 'bold',angle = 0, size=12),
        axis.title.y = element_text(face = 'bold',size=14),
        legend.position = 'none',
        strip.background.x = element_blank(),
        strip.text = element_text(face='bold',size=14)) + 
  scale_fill_viridis_d() + 
  labs(x=NULL,y="Mean GEBVs",
       title="Mean +/- Std. Error Selection Index GEBV by Cycle")