---
title: "Verify the pedigree entries"
author: "Marnin Wolfe"
date: "2021-May-12"
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, 
                      tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```

# Previous step

2.  [Get BLUPs combining all trial data](02-GetBLUPs.html): Combine data from all trait-trials to get BLUPs for downstream genomic prediction. Fit mixed-model to multi-trial dataset and extract BLUPs, de-regressed BLUPs and weights. Include two rounds of outlier removal.

# Read pedigree

Read the pedigree I downloaded for IITA

```{r}
library(tidyverse); library(magrittr)

ped<-read_delim(here::here("data/DatabaseDownload_2021May04","pedigree.txt"),
                delim = "\t")
ped_tms20<-readxl::read_xlsx(here::here("data/DatabaseDownload_2021Aug08","Pedigree.xlsx"))
ped %<>% 
  select(-Cross_Type) %>% 
  bind_rows(ped_tms20 %>% 
              select(Accession_name,Pedigree) %>% 
              rename(Accession=Accession_name) %>% 
              separate(Pedigree,c("Female_Parent","Male_Parent"),"/"))
#rm(ped_tms20)
```

Filter: Keep only complete pedigree records.

```{r}
ped %<>% 
  filter(!is.na(Female_Parent),
         !is.na(Male_Parent),
         Female_Parent!="?",
         Male_Parent!="?") %>% 
  distinct
```

Number of full-sib families?

```{r}
ped %>% distinct(Female_Parent,Male_Parent) %>% nrow()
```

Summarize distribution of full-sib family sizes

```{r}
ped %>%   
  count(Female_Parent,Male_Parent) %>% arrange(desc(n))  %>%  summary(.$n)
```

```{r}
ped %>% head
```

# Match SNP-to-pedigree names

Goal is to identify DNA samples names for listed accessions and parents in the pedigree. Important to choose same samples used in genomic predictions where possible, esp. phenotyped training clones.

However, there may be non-phenotyped clones that are genotyped genomic selection progeny, which I still want in my analysis.

First with the union of the parent and accession IDs in the pedigree.

```{r}
pednames<-union(ped$Accession,union(ped$Female_Parent,ped$Male_Parent)) %>% 
  tibble(germplasmName=.)
pednames$germplasmName %>% length # number of names in ped
```

Add a "Cohort" variable corresponding to the genetic groups or cycles in the germplasm.

```{r}
pednames %<>%
  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")))))))))
pednames %>% 
  count(Cohort)
```

The names in the pedigree downloaded from cassavabase should match the names in the **`germplasmName`** variable in the plot-basis pheno data, also downloaded from cassavabase. From that plot-basis data, make a data.frame of the unique **`germplasmName`-to-`FullSampleName`** matches, where **`FullSampleName`** is the column with names matching VCF files / DNA samples.

```{r}
dbdata<-readRDS(here::here("output","IITA_ExptDesignsDetected_2021Aug08.rds"))
phenos2genos<-dbdata %>% 
  distinct(GID,germplasmName,FullSampleName) %>% 
  filter(!is.na(FullSampleName))
phenos2genos %>% head
```
How many from each cohort in the pednames match a germplasmName with a FullSampleName in the plot-basis trial data?

```{r}
pednames %>%
  inner_join(phenos2genos) %>% 
  count(Cohort)
```

I want haplotypes of genotyped progeny in the pedigree *even* if they aren't phenotyped.

The `*.fam` file for the **RefPanelAndGSprogeny** VCF I will ultimately use to extract haplotypes and do predictions with is already in the `data/` directory because it's used in the [standard match-genos-to-phenos step](01-cleanTPdata.html#%5BUser_input%5D_Assign_genos_to_phenos).

```{r}
gids_in_fam<-read.table(here::here("data",
                                   "chr1_RefPanelAndGSprogeny_ReadyForGP_72719.fam"), 
                        stringsAsFactors = F, header = F)$V2
length(gids_in_fam)
```

Split the names based on a ":" for the GBS samples, and "\_A" seems to work (for IITA) to split the DArT sample ID from the germplasmName for DArTseqLD samples. The suffix / ID added to DNA sample names at DARt seems to vary depending on what gets submitted to them. GBS-era samples (almost) 100% had the ":" separator.

```{r}
gids_in_fam %<>%
  tibble(FullSampleName=.) %>%
  separate(FullSampleName,c("germplasmName","DNA_ID"),":|_A",remove = F) %>% 
  select(-DNA_ID)
gids_in_fam %>%
  filter(grepl("TMS18",FullSampleName)) %>% nrow()
```

There are also TMS20 in the latest DArT report and the DB sourced ped contains matches to the DARt sample names.
```{r}
gids_in_fam<-gids_in_fam %>% 
  bind_rows(ped_tms20 %>% 
              select(Accession_name,Lab_ID) %>% 
              rename(germplasmName=Accession_name,
                     FullSampleName=Lab_ID) %>% 
              mutate(FullSampleName=gsub("\\.","_",FullSampleName)))
gids_in_fam %>% head
```

```{r}
gids_in_fam %>%
  filter(grepl("TMS16|TMS17|TMS18|TMS19|TMS20",FullSampleName)) %>% nrow()
```

There are 3713 "TMS16" though "TMS20" clones genotyped.

```{r}
pednames2genos<-pednames %>% 
  inner_join(phenos2genos) %>% 
  bind_rows(pednames %>%
              anti_join(phenos2genos) %>% 
              inner_join(gids_in_fam))
pednames2genos %>% count(Cohort)  
```
Are there germplasmName in the `pednames2genos` match table with multiple DNA samples?
```{r}
pednames2genos %>% count(germplasmName) %>% arrange(desc(n))
```
Of course there are. Will need to pick.

Which pednames have BLUPs?
```{r}
blups<-readRDS(file=here::here("output","IITA_blupsForModelTraining_twostage_asreml_2021Aug09.rds"))
blups %>% 
  select(Trait,blups) %>% 
  unnest(blups) %>% 
  distinct(GID) %$% GID -> gidWithBLUPs
pednames %>%
  inner_join(phenos2genos) %>% 
  filter(FullSampleName %in% gidWithBLUPs) %>% 
  count(Cohort)
```


```{r}
# if there are any blups for a germplasmName
# keep only the FullSampleName/GID associated
# else keep all
pednames2genos %<>% 
  mutate(HasBLUPs=ifelse(GID %in% gidWithBLUPs,T,F)) %>% 
  nest(DNAsamples=-c(germplasmName,Cohort)) %>% 
  mutate(AnyBLUPs=map_lgl(DNAsamples,~any(.$HasBLUPs)),
         DNAsamples=ifelse(AnyBLUPs==T,
                           map(DNAsamples,~filter(.,HasBLUPs==TRUE)),
                           DNAsamples)) %>% 
  select(-AnyBLUPs) %>% 
  unnest(DNAsamples)
# Among all remaining
# Select only one GID to use for each germplasmName
pednames2genos %<>% 
  group_by(germplasmName) %>% 
  slice(1) %>% 
  ungroup()
pednames2genos %>% 
  count(Cohort)
```
Now make a pedigree with both Accession and parent names matching the genos (FullSampleName) rather than phenos (germplasmName).

```{r}
ped2genos<-ped %>% 
  rename(germplasmName=Accession) %>% 
  inner_join(pednames2genos %>% 
               select(-GID,-HasBLUPs)) %>% 
  left_join(pednames2genos %>% 
              select(-GID,-Cohort,-HasBLUPs) %>% 
              rename(Female_Parent=germplasmName,
                     DamID=FullSampleName)) %>% 
  left_join(pednames2genos %>% 
              select(-GID,-Cohort,-HasBLUPs) %>% 
              rename(Male_Parent=germplasmName,
                     SireID=FullSampleName)) 
ped2genos %<>% 
  filter(!is.na(FullSampleName),
         !is.na(DamID),
         !is.na(SireID))
ped2genos %>% distinct# %>% filter(grepl("TMS20",germplasmName))
```

In the end, considering only pedigree entries where the entire trio (offspring + both parents) are genotyped, the pedigree has 6257 entries to check.

```{r}
ped2genos %>% count(Cohort,DamID,SireID) %>% 
  ggplot(.,aes(x=Cohort,y=n,fill=Cohort)) + 
  geom_boxplot(notch = T) + theme_bw() +
  ggtitle("Distribution of family sizes (genotyped only)")
```

```{r}
ped2genos %>% 
  count(Cohort,DamID,SireID) %$% summary(n)
```

Number of families with at least 10 genotyped members, by cohort:

```{r}
ped2genos %>% 
  count(Cohort,DamID,SireID) %>% 
  filter(n>=10) %>% 
  count(Cohort)
```

```{r}
ped2genos %>% 
  select(FullSampleName,DamID,SireID) %>% 
  write.table(.,file=here::here("output","ped2genos.txt"),row.names=F, col.names=F, quote=F)
```





# Run plink IBD (--genome)

**Alternative to the below:** Could compute everything manually based on mendelian rules. Kinship coefficients directly from the relationship matrix used for prediction would also be useful, for example, the estimated inbreeding coefficient of an individual is 1/2 the relationship of its parents.

PLINK1.9 pipeline to use:

1.  Subset whole-pop. plink file (`AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08`) to only lines in the pedigree.
2.  LD-prune `--indep-pairwise 100 25 0.25` stringent, but somewhat arbitrary
3.  Compute IBD-relationships `--genome`
4.  Parent-offspring relationships determination (see below)

**Determine parent-offspring relationship status based on `plink` IBD:**

-   should have a kinship $\hat{\pi} \approx 0.5$.

-   Three standard IBD probabilities are defined for each pair; the probability of sharing zero (Z0), one (Z1) or two (Z2) alleles at a randomly chosen locus IBD.

-   The expectation for siblings in terms of these probabilities is Z0=0.25, Z1=0.5 and Z2=0.25.

-   The expectation for parent-offspring pairs is Z0=0, Z1=1 and Z2=0.

-   Based on work I did in *2016* (never published), declare a parent-offspring pair where: Z0\<0.313 and Z1\>0.668.

```{r, eval=F}
ped2check<-read.table(file=here::here("output","ped2genos.txt"),
                      header = F, stringsAsFactors = F)

pednames<-union(ped2check$V1,union(ped2check$V2,ped2check$V3)) %>% 
  tibble(FID=0,IID=.)
write.table(pednames,file=here::here("output","pednames2keep.txt"), 
            row.names = F, col.names = F, quote = F)
```

Checked plink's order-of-operations and combining --keep and --indep-pairwise 
in the same filter call should result in the correct ordering: 
first subset samples, then LD prune.

```{c, eval=F}
cd ~/IITA_2021GS/
export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;
plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
  --keep output/pednames2keep.txt \
  --indep-pairwise 100 25 0.25 \
  --genome \
  --out output/pednames_Prune100_25_pt25;
```

Creates a 2GB `*.genome`, >6000 samples samples worth of pairwise relationships. 

```{bash, eval=F}
#cd /home/jj332_cas/marnin/implementGMSinCassava/
#export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;
#plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
#  --indep-pairwise 100 25 0.25 --out output/Prune100_25_pt25;
#plink --bfile output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08 \
#  --extract output/Prune100_25_pt25.prune.in --genome \
#  --out output/AllChrom_RefPanelAndGSprogeny_ReadyForGP_2021Aug08_Prune100_25_pt25

# That wastefully creates a >40GB `*.genome` file with all pairwise relationships.

# Brute force solution is to read that, grab the needed relationships, and delete it....
```

```{r, eval=F}
library(tidyverse); library(magrittr); library(data.table)

genome<-fread(here::here("output/", 
                         "pednames_Prune100_25_pt25.genome"),
              stringsAsFactors = F,header = T) %>% 
  as_tibble

ped2check<-read.table(file=here::here("output","ped2genos.txt"),
                      header = F, stringsAsFactors = F)

head(genome)
# A tibble: 6 × 14
#    FID1 IID1      FID2 IID2    RT       EZ    Z0    Z1     Z2 PI_HAT   PHE   DST
#   <int> <chr>    <int> <chr>   <chr> <int> <dbl> <dbl>  <dbl>  <dbl> <int> <dbl>
# 1     0 TMS15F1…     0 TMS15F… OT        0 1     0     0      0         -1 0.720
# 2     0 TMS15F1…     0 TMS15F… OT        0 0.854 0.103 0.0434 0.0948    -1 0.763
# 3     0 TMS15F1…     0 TMS15F… OT        0 0.868 0.132 0      0.0661    -1 0.718
# 4     0 TMS15F1…     0 TMS15F… OT        0 1     0     0      0         -1 0.726
# 5     0 TMS15F1…     0 TMS15F… OT        0 0.407 0.593 0      0.296     -1 0.775
# 6     0 TMS15F1…     0 TMS15F… OT        0 1     0     0      0         -1 0.715
dim(genome)
# [1] 20470401       14
```

```{r, eval=F}
ped2check_genome<-genome %>% 
  semi_join(ped2check %>% rename(IID1=V1,IID2=V2)) %>% 
  bind_rows(genome %>% semi_join(ped2check %>% rename(IID1=V2,IID2=V1))) %>% 
  bind_rows(genome %>% semi_join(ped2check %>% rename(IID1=V1,IID2=V3))) %>% 
  bind_rows(genome %>% semi_join(ped2check %>% rename(IID1=V3,IID2=V1)))
saveRDS(ped2check_genome,file=here::here("output","ped2check_genome.rds"))
```

```{bash, eval=F}
cd ~/IITA_2021GS/output/;
rm pednames_Prune100_25_pt25.genome
```

# Verify and Correct Pedigree

```{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(ConfirmFemaleParent=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(ConfirmMaleParent=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)
```


```{r}
ped2check %>% 
  count(ConfirmFemaleParent,ConfirmMaleParent) %>% mutate(Prop=round(n/sum(n),2))
```
```{r}
ped2check %>% 
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS20",FullSampleName,ignore.case = T),"TMS20",
                       ifelse(grepl("TMS19",FullSampleName,ignore.case = T),"TMS19",
                              ifelse(grepl("TMS18",FullSampleName,ignore.case = T),"TMS18",
                                     ifelse(grepl("TMS17",FullSampleName,ignore.case = T),"TMS17",
                                            ifelse(grepl("TMS16",FullSampleName,ignore.case = T),"TMS16",
                                                   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"))))))))) %>% 
  filter(ConfirmFemaleParent=="Confirm",
         ConfirmMaleParent=="Confirm") %>% 
  count(Cohort,name = "BothParentsConfirmed")
```
I'm only interested in families / trios that are confirmed. Remove any without both parents confirmed.
```{r}
correctedped<-ped2check %>% 
  filter(ConfirmFemaleParent=="Confirm",
         ConfirmMaleParent=="Confirm") %>% 
  select(-contains("Confirm"))
```
```{r}
correctedped %>% 
  count(SireID,DamID) %>% arrange(desc(n))
```

```{r}
correctedped %>% 
  count(SireID,DamID) %$% summary(n)
```

```{r}
correctedped %>% 
  count(SireID,DamID) %>% 
  filter(n>=5) %>% arrange(desc(n)) %$% union(SireID,DamID) %>% 
  tibble(ParentsMoreThanFiveProg=.) %>%  
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS20",ParentsMoreThanFiveProg,ignore.case = T),"TMS20",
                       ifelse(grepl("TMS19",ParentsMoreThanFiveProg,ignore.case = T),"TMS19",
                              ifelse(grepl("TMS18",ParentsMoreThanFiveProg,ignore.case = T),"TMS18",
                                     ifelse(grepl("TMS17",ParentsMoreThanFiveProg,ignore.case = T),"TMS17",
                                            ifelse(grepl("TMS16",ParentsMoreThanFiveProg,ignore.case = T),"TMS16",
                                                   ifelse(grepl("TMS15",ParentsMoreThanFiveProg,ignore.case = T),"TMS15",
                                                          ifelse(grepl("TMS14",ParentsMoreThanFiveProg,ignore.case = T),"TMS14",
                                                                 ifelse(grepl("TMS13|2013_",ParentsMoreThanFiveProg,
                                                                              ignore.case = T),"TMS13","GGetc"))))))))) %>%
  count(Cohort, name="ParentsWithLeast5offspring")
```

```{r}
correctedped %>% 
  count(SireID,DamID) %>% 
  arrange(desc(n)) %$% union(SireID,DamID) %>% 
  tibble(Parents=.) %>%  
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS20",Parents,ignore.case = T),"TMS20",
                       ifelse(grepl("TMS19",Parents,ignore.case = T),"TMS19",
                              ifelse(grepl("TMS18",Parents,ignore.case = T),"TMS18",
                                     ifelse(grepl("TMS17",Parents,ignore.case = T),"TMS17",
                                            ifelse(grepl("TMS16",Parents,ignore.case = T),"TMS16",
                                                   ifelse(grepl("TMS15",Parents,ignore.case = T),"TMS15",
                                                          ifelse(grepl("TMS14",Parents,ignore.case = T),"TMS14",
                                                                 ifelse(grepl("TMS13|2013_",Parents,
                                                                              ignore.case = T),"TMS13","GGetc"))))))))) %>%
  count(Cohort, name="Parents")
```
```{r}
correctedped %>% 
  mutate(Cohort=NA,
         Cohort=ifelse(grepl("TMS20",FullSampleName,ignore.case = T),"TMS20",
                       ifelse(grepl("TMS19",FullSampleName,ignore.case = T),"TMS19",
                              ifelse(grepl("TMS18",FullSampleName,ignore.case = T),"TMS18",
                                     ifelse(grepl("TMS17",FullSampleName,ignore.case = T),"TMS17",
                                            ifelse(grepl("TMS16",FullSampleName,ignore.case = T),"TMS16",
                                                   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"))))))))) %>% 
  count(Cohort, name="BothParentsConfirmed")
```
### --> Write verified pedigree
```{r}
correctedped %>% 
  write.table(.,here::here("output","verified_ped.txt"),
              row.names = F, col.names = T, quote = F)
```

# Next Step

4.  [Preprocess data files](04-PreprocessDataFiles.html): Prepare haplotype and dosage matrices, GRMs, pedigree and BLUPs, genetic map *and* recombination frequency matrix, for use in predictions.
