---
title: "analysis_inbredlines_and_matrilines"
author: "Patrice"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Preliminaries
```{r, results='hide', message=FALSE}
  # Package names 
  packages <- c("lme4", "dplyr", "ggplot2", "scales", "gridExtra", "ggpubr", "tidyverse", "patchwork")
  
  # Install packages not yet installed
  installed_packages <- packages %in% rownames(installed.packages())
  if (any(installed_packages == FALSE)) {
    install.packages(packages[!installed_packages])
  }
  
  # Packages loading
  invisible(lapply(packages, library, character.only = TRUE))
  
  options(max.print=1300)
```

``` {r, results='hide', message=FALSE}
# useful function
antilogit<-function(x){exp(x)/(1+exp(x))}
```

```{r, results='hide', message=FALSE}
# setting plot theme
theme_set(theme_bw())
theme_update(text = element_text(size = 30,"Roboto Condensed Light"),
             plot.title = element_text(size = 30, face = "bold", hjust = 0.5),
             legend.title = element_text(size=30),
             legend.text = element_text(size=30),
             axis.text.y = element_text(size=30),
             axis.text.x = element_text(size=30),
             axis.title.x = element_text(size=30),
             axis.title.y = element_text(size=30),
             panel.grid = element_blank()
             )
```

# Data
``` {r, results='hide', message=FALSE}
#dataset inbredlines
inbredlines<- read.table("inbredlines.csv", header=TRUE, sep=";", na.strings="NA", dec=".", strip.white=TRUE)
#we start by eliminating individuals that turned out not to be cross fertlized (conformity to expected genotype=no, we keep NA and yes)
inb<-subset(inbredlines, is.na(inbredlines$conformity_to_expected_genotype)|inbredlines$conformity_to_expected_genotype=="yes")
str(inb)
g0inb<- read.table("inbredlines_mothers.csv", header=TRUE, sep=";", na.strings="NA", dec=".", strip.white=TRUE)
str(g0inb)

#datasets g0 and g1 for HFR-LFR
g0<- read.table("G0HFRLFR.csv", header=TRUE, sep=";", na.strings="NA", dec=".", strip.white=TRUE)
g1<- read.table("G1HFRLFR.csv", header=TRUE, sep=";", na.strings="NA", dec=".", strip.white=TRUE)
str(g0)
```

## Compute semi-sterility and sterility status of mothers

```{r}

#we use a threshold of 10 outcrossed offspring to declare fertility

#HFR&LFR dataset
#the "g0" file contains all data of babies pigmented and albino produced by  "g0" mothers of the experimental individuals (G1)
#the mothers are not repeated 
#mothers have been crossed with fathers and the offspring laid by the FATHERS serve as indicators of male-sterility of the MOTHER !!!
# the fathers are all albino, mothers may be pigmented homozygous ("pig") or pigmented heterozygous ("pigh") 
# the sterility status of the mother is defined using a threshold of number of outcrossed babies produced by partners :  0 for sterility, between 1 and 9  for semisterility, at least 10 of or fertility. Beware that the corresponding numbers of pigmented babies (to which unhatched pigmented, if any, must be added) is 10 if the focal is pigmented homozygous ("pig") but 5 when the mother is heterozygous ("pigh").

#I first create a column for NHP (not hatched pigmented) with 0 instead of Na (the unhatched are not counted if there are already many hatched pigmented, ie largely enough to declare the individual male-fertile)
g0$nhpnum<-ifelse(is.na(g0$NHP),0,g0$NHP)
#then i compute the total number of pigmented, hatched or unhatched
g0$pigtot<-g0$nhpnum+g0$pignum
#then I compute a heterozygosity factor that multiplies nb pigmented by 2 to fiund outcrossed offspring if mother is heterozygous
g0$h<-ifelse(g0$morph_mother=="pigh",2,1)

#then I compute the mat_status10 variable
g0$mats10<-case_when(g0$pigtot*g0$h >= 10 ~ 'fertile'
          ,g0$pigtot*g0$h >0 ~ 'semisterile'
          ,g0$pigtot*g0$h ==0 ~ 'sterile'
          ,TRUE  ~ NA_character_)
g0clean<-subset(g0,is.na(g0$mats10)==F)


#inbredlines dataset

g0inb$nhpnum<-ifelse(is.na(g0inb$NHP_ftmother),0,g0inb$NHP_ftmother)
g0inb$pigtot<-g0inb$nhpnum+g0inb$num_pigmother

levels(as.factor(g0inb$mother_test__partner_no_or_new))
#partners can be albino (=alb-1) or albinew (=alb-2)
#compute the expected frequency of pigmented in the progeny
#the frequency is 100% if the partner is albinew (complementation), or if the mother is homozygous pigmented (pig)
#it is 50% if the partner is albino and the mother is heterozygous pigmented (pigh)

g0inb$exppigmum<-ifelse(g0inb$mother_test__partner_no_or_new=="new"|g0inb$morph_mother=="pig",1,
                      ifelse((g0inb$morph_mother=="pigh")&(g0inb$mother_test__partner_no_or_new=="no"),0.5,NA))
#g0inb$exppigmum

g0inb$mats10<-case_when(    g0inb$pigtot==0 ~ 'sterile'
          ,g0inb$pigtot < 10*g0inb$exppigmum  ~ 'semisterile'
          ,g0inb$pigtot >= 10*g0inb$exppigmum ~ 'fertile'
          ,TRUE~NA_character_ )
g0inbclean<-subset(g0inb,is.na(g0inb$mats10)==F)
#g0inbclean$mats10
```

# Tables and graphs of maternal status
```{r}

#mothers in the inbred line dataset
#rename levels in conformity with the names used in the paper
levels(as.factor(g0inbclean$matriline_short))
g0inbclean$maternal_type<-case_when(g0inbclean$matriline_short == "KI" ~ 'KS'
          ,g0inbclean$matriline_short == "WI"~ 'DS'
          ,g0inbclean$matriline_short == "KF" ~ 'KF'
          ,TRUE  ~ 'NA')
#number of mothers per category, inbred lines dataset
table(g0inbclean$maternal_type, g0inbclean$mats10)
#proportions of mothers per category, inbred lines dataset
table(g0inbclean$maternal_type, g0inbclean$mats10)/as.vector(table(g0inbclean$maternal_type))

# Figure 2a
figure_2a <- ggplot(
  g0inbclean,aes(x = maternal_type,fill = mats10)) +
  geom_bar(position = "fill", color = "white", size = 1) +  # White borders
  scale_fill_manual(
    name = "Male-Fertility Status",  # Legend title
    values = c("#0073e6", "#c44601", "#f57600")  # Colorblind-friendly colors
  ) +  
  #guides(fill = guide_legend(title = NULL)) +
  labs(
    x = "Maternal Type", 
    y = "Proportions",
    title = "Inbred lines dataset"  # Adding the title
  )

figure_2a

# Figure 2b
#mothers in the HFRLFR dataset
levels(as.factor(g0clean$G0mother_short))
g0clean$maternal_type<-case_when(g0clean$G0mother_short == "KI" ~ 'KS'
          ,g0clean$G0mother_short == "DI"~ 'DS'
          ,g0clean$G0mother_short == "KF" ~ 'KF'
          ,TRUE  ~ 'NA')
#number of mothers per category, HFRLFR dataset
table(g0clean$maternal_type, g0clean$mats10)
#proportions
table(g0clean$maternal_type, g0clean$mats10)/as.vector(table(g0clean$maternal_type))

figure_2b <- ggplot(
  g0clean, aes(x = maternal_type, fill = mats10)) +
  geom_bar(position = "fill", color = "white", size = 1) +  # White borders
  scale_fill_manual(
    name = "Male-Fertility Status",  # Legend title
    values = c("#0073e6", "#c44601", "#f57600")  # Colorblind-friendly colors
  ) +  
  #guides(fill = guide_legend(title = NULL)) +
  labs(
    x = "Maternal Type", 
    y = "Proportions",
    title = "Selected lines dataset"  # Adding the title
  )

figure_2b

figure2_final <- figure_2a+figure_2b +
  plot_layout(guides = "collect")+
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 35),
        legend.position = 'bottom')

### to save it in an readable version
ggsave("figure2.png", figure2_final, width = 14, height = 10, dpi = 600)

######################

#table and proportions over the whole dataset
table(g0clean$maternal_type, g0clean$mats10)+table(g0inbclean$maternal_type, g0inbclean$mats10)

(table(g0clean$maternal_type, g0clean$mats10)+table(g0inbclean$maternal_type,g0inbclean$mats10))/(as.vector(table(g0clean$maternal_type)+table(g0inbclean$maternal_type)))

```

# Tests on frequency of sterility among maternal groups

```{r}

#inbred lines dataset
#create binomial variable male fertility in 1/0

g0inbclean$matsbinary<-cbind(ifelse(g0inbclean$mats10=="fertile",1,0),1-ifelse(g0inbclean$mats10=="fertile",1,0))
mod1inbmat<-glmer(matsbinary~1+maternal_type+(1|block),data=g0inbclean,family=binomial)
summary(mod1inbmat)
mod2inbmat<-glmer(matsbinary~1+(1|block),data=g0inbclean,family=binomial)
summary(mod2inbmat)
anova(mod1inbmat,mod2inbmat)

#HFRLFR lines dataset
#create binomial variable male fertility in 1/0

g0clean$matsbinary<-cbind(ifelse(g0clean$mats10=="fertile",1,0),1-ifelse(g0clean$mats10=="fertile",1,0))
mod1mat<-glmer(matsbinary~1+maternal_type+(1|block),data=g0clean,family=binomial)
summary(mod1mat)
mod2mat<-glmer(matsbinary~1+(1|block),data=g0clean,family=binomial)
summary(mod2mat)
anova(mod1mat,mod2mat)


```


# Compute semi-sterility an sterility status in the G1 dataset :LFR versus HFR

```{r}

#we use a threshold of 10 outcrossed offspring to declare male fertility

#first step is to compute the sterility status of mothers
#note that in the g1 dataset, data for each mother may be repeated as many times as she has offspring in the dataset

# first get a numerical version of the number of unhatched pigmented embryos in the progeny of the mother's partner with 0 instead of NA; note that the NAs represent counts that have not been done because they would not have changed the fertility status (enough had been counted in the live babies already);putting a zero instead has no consequence
g1$NHPmatnum<-ifelse(is.na(g1$NHP_ftmother),0,g1$NHP_ftmother)
#then sum all pigmented found in the mother's partner progeny
g1$pigmattot<-g1$num_pigmother+g1$NHPmatnum
#then compute the threshold for number of pigmented that applies to each particular mother; partner is always albino but mother may be homozygous pigmented ("pig", in which case the number of pigmented equals number of sired offspring) or heterozygous ("pigh", in which case the number of sired offspring is estimated by multiplying the pigmented by two)
g1$pigmaxmat<-ifelse(g1$morph_mother=="pig",10,ifelse(g1$morph_mother=="pigh",5,NA))
#now give to each mother its male sterility status
g1$mats10<-case_when(g1$pigmattot >= g1$pigmaxmat ~ 'fertile'
          ,g1$pigmattot >0 ~ 'semisterile'
          ,g1$pigmattot==0 ~ 'sterile'
          ,TRUE  ~ 'NA')

#g1$mats10
#note there is one missing maternal status because the mother's partner died before it could lay eggs..

#second step is to compute the sterility status of g1
# first get a numerical version of the number of unhatched pigmented embryos in the progeny of the g1's partner with 0 instead of NA; note that the NAs represent counts that have not been done because they would not have changed the fertility status (enough had been counted in the live babies already);putting a zero instead has no consequence
g1$NHPnum<-ifelse(is.na(g1$NHP_fertility_test),0,g1$NHP_fertility_test)
#then sum all pigmented found in the g1's partner progeny
g1$pigtot<-g1$numpig+g1$NHPnum
#there are a few NA's that we keep, these are cases were could not test that the partner was able to lay (eg died before etc)
#the threshold applicable to g1's depends on their own genotype which is usually "pigh" (pig heterozygous) but in a few occasions albino ("alb") and on the genotype of the partner that may be "no" (albino alb-1) or "new" (albinew alb-2). We used albino partners to test at the same time that the tested individuals were indeed heterozygous (as expected if the albino father HFR or LFR has indeed inseminated the K mother), or homozygous (in case the k mother self-fertilized and was itself homozygous pig, case of the fertile K from lyon). in all instances that could be tested, we never obtained ratios consistent with homozygous status (pig) of the mother (this would be 100% pigmented progeny with an albino partner) so we conclude that the crosses were ok.
#in total the threshold 10 should be divided by two when the mother is pigh and the partner alb, remain the same when the partner is albinew (or when the mother is pig homozygous but this is never the case in this dataset)
g1$pigmax<-ifelse(g1$expected_genotype_g1=="pigh"&g1$morph_partner=="no",5,ifelse(g1$morph_partner=="new",10,NA))

#now give to each g1 its male sterility status
g1$s10<-case_when(g1$pigtot >= g1$pigmax ~ 'fertile'
          ,g1$pigtot >0 ~ 'semisterile'
          ,g1$pigtot==0 ~ 'sterile'
          ,TRUE  ~ 'NA')


#third step we simplify down by collapsing the sterile and semisterile category
#there remains 1 for fertile 0 for sterile

g1$fertilemother<-ifelse(g1$mats10=="sterile"|g1$mats10=="semisterile",0,ifelse(g1$mats10=="fertile",1,NA))
g1$fertileoffspring<-ifelse(g1$s10=="sterile"|g1$s10=="semisterile",0,ifelse(g1$s10=="fertile",1,NA))

#finally we keep only offspring whose status is known
g1clean<-subset(g1,is.na(g1$fertileoffspring)==FALSE)


```
# Tables and graphs in the G1 dataset :LFR versus HFR

```{r}

str(g1clean)
#rename maternal and paternal types in the HFRLFR dataset
levels(as.factor(g1clean$matriline_short))
g1clean$maternal_type<-case_when(g1clean$matriline_short == "KI" ~ 'KS'
          ,g1clean$matriline_short == "WI"~ 'DS'
          ,g1clean$matriline_short == "KF" ~ 'KF'
          ,TRUE  ~ 'NA')

table(g1clean$maternal_type, g1clean$mats10)

g1clean$paternal_type<-case_when(g1clean$paternal_line == "HFR_2" ~ 'HFR'
          ,g1clean$paternal_line == "LFR_2"~ 'LFR'
          ,TRUE  ~ 'NA')
#note that the number of fertile mothers in DS and of nonfertile mothers in KF are
#respectively 1 and 2 so there is no way we can test the effect of maternal phenotype
#within the DS and KF sub-datasets

#create a combined maternal and paternal type in the HFRLFR dataset
g1clean$matpat<-paste(g1clean$maternal_type,g1clean$paternal_type)
levels(as.factor(g1clean$matpat))

#tabulate fertility data for g1 offspring by maternal and paternal type
table(g1clean$matpat, g1clean$s10)

table(g1clean$matpat, g1clean$s10)/as.vector(table(g1clean$matpat))

#plot fertility data for g1 offspring by maternal and paternal type
library(tidyr)
# Split the 'matpat' variable into two columns (assuming space as a delimiter)


figure_3a<-ggplot(g1clean,aes(x = matpat,fill = s10)) +
  geom_bar(position = "fill", color = "white", size = 1) +  
  theme_update(text = element_text(size = 30,"Roboto Condensed Light"),
             plot.title = element_text(size = 30, face = "bold", hjust = 0.5),
             legend.title = element_text(size=30),
             legend.text = element_text(size=30),
             axis.text.y = element_text(size=30),
             axis.text.x = element_text(size=30),
             axis.title.x = element_text(size=30),
             axis.title.y = element_text(size=30),
             panel.grid = element_blank()
             )+
  scale_fill_manual(
    name = "Male-Fertility Status",  # Legend title
    values = c("#0073e6", "#c44601", "#f57600")  # Colorblind-friendly colors
  ) +  
  labs(
    x= "",
    y = "Proportions",
    title = "All data"  # Adding the title
  )

figure_3a


```
# GLMMs in the G1 dataset :LFR versus HFR

```{r}

#now let us study the difference between the different maternal/paternal origins, not taking maternal phenotype into consideration

g1clean$y<-cbind(g1clean$fertileoffspring,1-g1clean$fertileoffspring)
mod4<-glmer(y~1+maternal_type+paternal_type+maternal_type:paternal_type + (1|mother)+ (1|block),family="binomial",data=g1clean)
summary(mod4)

mod5<-glmer(y~maternal_type+paternal_type+ (1|mother)+ (1|block),family="binomial",data=g1clean)
summary(mod5)

anova(mod4,mod5)
#interaction not significant

mod6<-glmer(y~maternal_type+ (1|mother)+ (1|block),family="binomial",data=g1clean)
summary(mod6)

anova(mod5,mod6)
#the patriline effect is barely significant , this is due to the slight difference observed in the KI (more sterility in KI-HFR than in KI-LFR); in the other two categories (KF and W) there is no visible difference but the frequencies are so close to 0 or 100% that we cannot exclude it either

mod7<-glmer(y~paternal_type+ (1|mother)+ (1|block),family="binomial",control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)),data=g1clean)
summary(mod7)

anova(mod5,mod7)
#the matriline effect is indeed significant; the KFs are entirely restored and the W are never restored, which makes the model difficult to fit, as it mixes up very different categories... 

mod5_2<-glmer(y~maternal_type+paternal_type+ (1|block),family="binomial",data=g1clean)
summary(mod5_2)

anova(mod5,mod5_2)

mod5_3<-glmer(y~maternal_type+paternal_type+ (1|mother),family="binomial",data=g1clean)
summary(mod5_3)

anova(mod5,mod5_3)



#we  proceed by testing the possible difference among the two paternal lines within each maternal type separately
#in the KS

g1KS<-subset(g1clean,maternal_type=="KS")
mod7bis<-glmer(y~paternal_line+ (1|mother)+ (1|block),family="binomial",data=g1KS)
summary(mod7bis)
mod8bis<-glmer(y~1+ (1|mother)+ (1|block),family="binomial",data=g1KS)
summary(mod8bis)
anova(mod7bis,mod8bis)
mod9bis<-glmer(y~paternal_line+ (1|block),family="binomial",data=g1KS)
summary(mod9bis)
anova(mod7bis,mod9bis)

#in the KF

g1KF<-subset(g1clean,maternal_type=="KF")

mod7ter<-glmer(y~paternal_line+ (1|mother),family="binomial",data=g1KF)
summary(mod7ter)
mod8ter<-glmer(y~1+ (1|mother),family="binomial",data=g1KF)
summary(mod8ter)
anova(mod7ter,mod8ter)
mod9ter<-glm(y~paternal_line,family="binomial",data=g1KF)
summary(mod9ter)
anova(mod7ter,mod9ter)

#in the DS
g1D<-subset(g1clean,maternal_type=="DS")

mod7quat<-glmer(y~paternal_line+ (1|mother)+(1|block),family="binomial",data=g1D)
summary(mod7quat)
mod8quat<-glmer(y~1+ (1|mother)+(1|block),family="binomial",data=g1D)
summary(mod8quat)
anova(mod7quat,mod8quat)
mod9quat<-glmer(y~paternal_line+ (1|block),family="binomial",data=g1D)
summary(mod9quat)
anova(mod7quat,mod9quat)

#we now make a mode detailed model, only within the KS category, where we take into account potential effects of maternal phenotype

#we select only the data where the maternal phenotype is known
g1KSclean<-subset(g1KS,is.na(mats10)==F)
g1KSclean$matsbinary<-ifelse(g1KSclean$mats10=="fertile","mother\nmale-fertile","mother\nnot male-fertile")


g1KSclean$pat_matfert<-paste(g1KSclean$paternal_type,g1KSclean$matsbinary)

table(g1KSclean$pat_matfert,g1KSclean$s10)

table(g1KSclean$pat_matfert,g1KSclean$s10)/as.vector(table(g1KSclean$pat_matfert))

g1KSclean$pat_matfert <- gsub("HFR", "HFR\n", g1KSclean$pat_matfert)
g1KSclean$pat_matfert <- gsub("LFR", "LFR\n", g1KSclean$pat_matfert)

#plot fertility data within the KS category, as a function of paternaltype and maternal phenotype
figure_3b <- ggplot(g1KSclean,aes(x = pat_matfert,fill = s10)) +
  geom_bar(position = "fill", color = "white", size = 1) +  # White borders
  scale_fill_manual(
    name = "Male-Fertility Status",  # Legend title
    values = c("#0073e6", "#c44601", "#f57600")  # Colorblind-friendly colors
  ) +  
  labs(
    x= "",
    y = "Proportions",
    title = "KS Only"  # Adding the title
  )
figure_3b 

#### final figure 3
figure3_final <- figure_3a+figure_3b +
  plot_layout(guides = "collect")+
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 40),
        legend.position = 'bottom')


## to save it in an readable version
ggsave("figure3.png", figure3_final, width = 23, height = 12, dpi = 600)

##### 

#glmm on the KS dataset with maternal phenotype included in the predictors


g1KSclean$motherfertility<-as.factor(g1KSclean$matsbinary)

mod10<-glmer(y~motherfertility+paternal_type+motherfertility:paternal_type+ (1|mother)+(1|block),family="binomial",data=g1KSclean)
summary(mod10)

mod11<-glmer(y~motherfertility+paternal_type+ (1|mother)+(1|block),family="binomial",data=g1KSclean)
summary(mod11)
anova(mod10,mod11)


mod12<-glmer(y~motherfertility+(1|mother)+(1|block),family="binomial",data=g1KSclean)
summary(mod12)
anova(mod11,mod12)

mod13<-glmer(y~paternal_type+ (1|mother)+(1|block),family="binomial",data=g1KSclean)
summary(mod13)
anova(mod11,mod13)

mod10bis<-glmer(y~motherfertility+paternal_type+motherfertility:paternal_type+(1|block),family="binomial",data=g1KSclean)
summary(mod10bis)
anova(mod10,mod10bis)


```
# Compute semi-sterility an sterility status in the G1 dataset :inbred lines

```{r}
#let us choose a threshold of 10 sired offspring, including unhatched ones, to declare fertility

#sterility test for mothers

#we first compute the number of pigmented offspring laid by the mother's partner, including unhatched ones if there are some

inb$pigmatemum<-inb$num_pigmother+ifelse(is.na(inb$NHP_ftmother_num),0,inb$NHP_ftmother_num)
levels(as.factor(inb$morph_mother))
#mothers can be alb,pig or pigh
levels(as.factor(inb$mother_test__partner_no_or_new))
#partners can be albino (=alb-1) or albinew (=alb-2)
#compute the expected frequency of pigmented in the outcrossed progeny
#the frequency is 100% if the partner is albinew (complementation), or if the mother is homozygous pigmented (pig)
#it is 50% if the partner is albino and the mother is heterozygous pigmented (pigh)

inb$exppigmum<-ifelse(inb$mother_test__partner_no_or_new=="new"|inb$morph_mother=="pig",1,
                      ifelse((inb$morph_mother=="pigh")&(inb$mother_test__partner_no_or_new=="no"),0.5,NA))
inb$mats10<-case_when(    inb$pigmatemum==0 ~ 'sterile'
          ,inb$pigmatemum < 10*inb$exppigmum  ~ 'semisterile'
          ,inb$pigmatemum >= 10*inb$exppigmum ~ 'fertile'
          ,TRUE~NA_character_ )


#sterility test for G1 focal offspring

#we first compute the number of pigmented offspring laid by the G1's partner, including unhatched ones if there are some

inb$pigmateg1<- inb$P_babies_num +ifelse(is.na(inb$NHP_num),0,inb$NHP_num)
levels(as.factor(inb$G1_morph))
#G1s can be alb,pig or pigh
levels(as.factor(inb$morph_partner_albino_or_albinew))
#partners can be albino (=alb-1) or albinew (=alb-2)
#compute the expected frequency of pigmented in the outcrossed progeny
#the frequency is 100% if the partner is albinew (complementation), or if the mother is homozygous pigmented (pig)
#it is 50% if the partner is albino and the mother is heterozygous pigmented (pigh)

inb$exppigg1<-ifelse(inb$morph_partner_albino_or_albinew=="new"|inb$G1_morph=="pig",1,
                      ifelse((inb$G1_morph=="pigh")&(inb$morph_partner_albino_or_albinew=="no"),0.5,NA))
inb$s10<-case_when(    inb$pigmateg1==0 ~ 'sterile'
          ,inb$pigmateg1 < 10*inb$exppigg1  ~ 'semisterile'
          ,inb$pigmateg1 >= 10*inb$exppigg1 ~ 'fertile'
          ,TRUE~NA_character_ )


#prepare the binomial variable to be analysed
inb$y<-cbind(ifelse(inb$s10=="fertile",1,0),1-ifelse(inb$g1sterility=="fertile",1,0))

#clean dataset without NA in the offspring fertility status
inbclean<-subset(inb, is.na(inb$s10)==F)


```
# Tables and graphs in the G1 dataset : inbred lines

```{r}

str(inbclean)
#rename maternal and paternal types in the HFRLFR dataset
levels(as.factor(inbclean$matriline_short))
inbclean$maternal_type<-case_when(inbclean$matriline_short == "KI" ~ 'KS'
          ,inbclean$matriline_short == "WI"~ 'DS'
          ,inbclean$matriline_short == "KF" ~ 'KF'
          ,TRUE  ~ 'NA')

table(inbclean$maternal_type, inbclean$mats10)

levels(as.factor(inbclean$paternal_origin))
inbclean$paternal_pop<-case_when(inbclean$paternal_origin == "Lyon" ~ 'Lyon'
          ,inbclean$paternal_origin == "montpellier"~ 'Montpellier'
          ,TRUE  ~ 'NA')
levels(as.factor(inbclean$paternal_pop))
#note that the number of fertile mothers in DS and of nonfertile mothers in KF are
#respectively 1 and 2 so there is no way we can test the effect of maternal phenotype
#within the DS and KF sub-datasets

#create a combined maternal and paternal type in the inbred lines dataset
inbclean$matpat<-paste(inbclean$maternal_type,inbclean$paternal_pop)
levels(as.factor(inbclean$matpat))

#tabulate fertility data for g1 offspring by maternal and paternal type
table(inbclean$matpat, inbclean$s10)

table(inbclean$matpat, inbclean$s10)/as.vector(table(inbclean$matpat))

#plot fertility data for g1 offspring by maternal and paternal type

inbclean$matpat <- gsub("KS Montpellier", "KS\nMontpellier", inbclean$matpat)
inbclean$matpat <- gsub("KF Montpellier", "KF\nMontpellier", inbclean$matpat)
inbclean$matpat <- gsub("KF Lyon", "KF\nLyon", inbclean$matpat)
inbclean$matpat <- gsub("DS Montpellier", "DS\nMontpellier", inbclean$matpat)
inbclean$matpat <- gsub("DS Lyon", "DS\nLyon", inbclean$matpat)
inbclean$matpat <- gsub("KS Lyon", "KS\nLyon", inbclean$matpat)

figure_4 <- ggplot(inbclean,aes(x = matpat,fill = s10)) +
  geom_bar(position = "fill", color = "white", size = 1) +  # White borders
  scale_fill_manual(
    name = "Male-Fertility Status",  # Legend title
    values = c("#0073e6", "#c44601", "#f57600")  # Colorblind-friendly colors
  ) +  
  labs(
    x= "",
    y = "Proportions") +
  theme(legend.position = "bottom")
figure_4

## to save it in an readable version
ggsave("figure4.png", figure_4, width = 10, height = 8, dpi = 600)

```
# GLMMs in the G1 dataset :inbred lines

```{r}

#now let us study the difference between the different maternal /paternal origins, not taking maternal phenotype into consideration

model1<-glmer(
  inbclean$y~1+paternal_pop+maternal_type+paternal_pop:maternal_type+(1+maternal_type|paternal_line)+(1|mother)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model1)

model2<-glmer(
  inbclean$y~1+paternal_pop+maternal_type+(1+maternal_type|paternal_line)+(1|mother)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model2)

anova(model1,model2)

model3<-glmer(
  inbclean$y~1+maternal_type+(1+maternal_type|paternal_line)+(1|mother)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model3)
anova(model2,model3)

model4inb<-glmer(
  inbclean$y~1+paternal_pop+(1+maternal_type|paternal_line)+(1|mother)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model4inb)
anova(model2,model4inb)

model1bis<-glmer(
  inbclean$y~1+paternal_pop+maternal_type+paternal_pop:maternal_type+(1+maternal_type|paternal_line)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model1bis)
anova(model1,model1bis)

model1ter<-glmer(
  inbclean$y~1+paternal_pop+maternal_type+paternal_pop:maternal_type+(1|mother)+(1|paternal_line)+(1|block),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model1ter)
anova(model1,model1ter)

model1quat<-glmer(
  inbclean$y~1+paternal_pop+maternal_type+paternal_pop:maternal_type+(1+maternal_type|paternal_line)+(1|mother),  family="binomial", data=inbclean,
  control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model1quat)
anova(model1,model1quat)
```
# Models within the DS type, inbred lines
```{r}

#separate models in different maternal contexts : DS type
#what is happening in the DS maternal type


inbcleanW<-subset(inbclean,inbclean$maternal_type=="DS")


model16<-glmer(y~1+paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block)
              ,family="binomial", data=inbcleanW
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model16)

model17<-glmer(y~1+
              (1|paternal_line)+(1|mother)+(1|block)
              ,family="binomial", data=inbcleanW
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model17)

anova(model16,model17)

model16bis<-glmer(y~1+paternal_pop
              +(1|paternal_line)+(1|block)
              ,family="binomial", data=inbcleanW
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model16bis)
anova(model16,model16bis)

model16ter<-glmer(y~1+paternal_pop
              +(1|mother)+(1|block)
              ,family="binomial", data=inbcleanW
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model16ter)
anova(model16,model16ter)


#############################################################
##are the differences among patrilines within the montpellier or the Lyon set of lines or both?
#let us first evaluate in MOntpellier
inb2WM<-subset(inbcleanW,inbcleanW$paternal_pop=="Montpellier")
model11M<-glmer(y~1+(1|paternal_line)
              +(1|mother)+(1|block),data=inb2WM
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model11M)
model12M<-glmer(y~1+
              (1|mother)+(1|block),data=inb2WM
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model12M)
anova(model11M,model12M)


#are there differences within the Lyon set of lines?
inb2WL<-subset(inbcleanW,inbcleanW$paternal_pop=="Lyon")
model13L<-glmer(y~1+(1|paternal_line)
              +(1|mother)+(1|block),data=inb2WL
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model13L)
model14L<-glmer(y~1+
              (1|mother)+(1|block),data=inb2WL
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model14L)
anova(model13L,model14L)

##look for patrilines with above-average male-fertility potential in W context
##below a script to plot estimates of paternal line means with SD and 95 CI
#we work with model 17
aa<-ranef(model17) #extract random effects
bb<-as.data.frame(aa)#convert them to data frame
print(bb)#print them
bbpat<-subset(bb,bb$grpvar=="paternal_line") #keep only paternal estimates

summary(model17)
#untransform estimates to put them into the natural scale
#we untransform also estimate+-SD and estimate+-1.95SD to plot sd and confidence intervals
#we use the estimate of the grand intercept -1.8434 from model17
bbpat$y<-antilogit(-1.8434+bbpat$condval)
bbpat$yplussd<-antilogit(-1.8434+bbpat$condval+bbpat$condsd)
bbpat$yminussd<-antilogit(-1.8434+bbpat$condval-bbpat$condsd)
bbpat$yplus95<-antilogit(-1.8434+bbpat$condval+1.95*bbpat$condsd)
bbpat$yminus95<-antilogit(-1.8434+bbpat$condval-1.95*bbpat$condsd)
bbpat$paternal_pop<-ifelse(substr(bbpat$grp, 1,1)=="N","Lyon","Montpellier")

linevaluesDS<-bbpat

#a plot of paternal line means +/-1SD, with a horizontal bar for the grand mean
figure_5a <- ggplot(bbpat, aes(x=grp, y=y, color=paternal_pop, shape=paternal_pop)) +
  guides(color=guide_legend(title=NULL), shape=guide_legend(title=NULL)) +
  ylim(0, 1) + 
  geom_pointrange(aes(ymin=yminussd, ymax=yplussd), size = 3) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  geom_hline(yintercept=antilogit(-1.8434)) +
  scale_color_manual(values = c( "#0072B2", "#E69F00")) + # Colorblind-friendly colors
  scale_shape_manual(values = c(16, 15))+  # 16 = Circle, 15 = Square
  labs(
    x= "Paternal Inbred Lines",
    y = "Estimates Male Fertility")

figure_5a

#a plot of paternal line means with 95CI, with a horizontal bar for the grand mean
figure_5a_95 <- ggplot(bbpat, aes(x=grp, y=y, color=paternal_pop, shape=paternal_pop)) +
  guides(color=guide_legend(title=NULL), shape=guide_legend(title=NULL)) +
  ylim(0, 1) + 
  geom_pointrange(aes(ymin=yminus95, ymax=yplus95), size = 1,show.legend = FALSE)  +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=15)) +
  geom_hline(yintercept=antilogit(-1.8434), size = 1, show.legend = FALSE) +
  scale_color_manual(values = c( "#0072B2", "#E69F00")) + # Colorblind-friendly colors
  scale_shape_manual(values = c(16, 15))+  # 16 = Circle, 15 = Square
  labs(
    x= "Paternal inbred lines",
    y = "Estimated male-fertility",
    title = "DS Mothers")

figure_5a_95


#in the plots paternal lines are ordered by effect
#we see that some lines seem to have higher frequencies of fertility than others
#these are lines A24, N34, T4 (N80, A3, A14, N55 to a lesser extent)
#we must keep these lines in mind and investigate them further
#howewer none of them is compatible with a dominant restorer because the frequency of fertility is always <50%

#test of maternal phenotype in the WS dataset

#how many of the individuals have a male-fertile mother
inbcleanW$mumst_binary<-as.factor(ifelse(inbcleanW$mats10=="fertile","fertile","not fertile"))
inb2WnoNA<-subset(inbcleanW,is.na(inbcleanW$mumst_binary)==FALSE)
#inb2WnoNA$mumst_binary
str(inb2WnoNA)
table(inb2WnoNA$mumst_binary,inb2WnoNA$s10)

table(inb2WnoNA$mumst_binary,inb2WnoNA$s10)/as.vector(table(inb2WnoNA$mumst_binary))

#is there a difference in the probability to be fertile depending on whether the mother is sterile versus fertile
model91<-glmer(inb2WnoNA$y~1+mumst_binary+paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block),
              data=inb2WnoNA
              ,family="binomial"
              ,control= glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model91)
model92<-glmer(inb2WnoNA$y~1+paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block),
              data=inb2WnoNA
              ,family="binomial"
              ,control= glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model92)
anova(model91,model92)
#the effect of maternal sterility is barely significant, there might be a few fertile mothers that have a more fertile background or mitotype; 


#plot fertility data within the DS category, as a function of paternal pop and maternal phenotype
inb2WnoNA$pat_matfert<-paste(inb2WnoNA$paternal_pop,inb2WnoNA$mumst_binary)
table(inb2WnoNA$pat_matfert,inb2WnoNA$s10)
#numbers are too low to represent montpellier and Lyon separately
figureS4A_a <- ggplot(inb2WnoNA,aes(x = mumst_binary,fill = s10)) + 
  geom_bar(position = "fill")+ 
  scale_fill_manual(values=c("#0073e6", "#c44601", "#f57600")) +
    labs(title = "Maternal type = DS",
          x = NULL, y = NULL) +
  guides(fill = guide_legend(title = NULL)) +
  theme(axis.text.x = element_text(size = 30))
figureS4A_a
```
# Models within the KS type, inbred lines
```{r}

#separate models in different maternal contexts : DS type
#what is happening in the DS maternal type


inbcleanKS<-subset(inbclean,inbclean$maternal_type=="KS")


model26<-glmer(y~1+paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block)
              ,family="binomial", data=inbcleanKS
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model26)

model27<-glmer(y~1+
              (1|paternal_line)+(1|mother)+(1|block)
              ,family="binomial", data=inbcleanKS
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model27)

anova(model26,model27)

model26bis<-glmer(y~1+paternal_pop
              +(1|paternal_line)+(1|block)
              ,family="binomial", data=inbcleanKS
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model26bis)
anova(model26,model26bis)

model26ter<-glmer(y~1+paternal_pop
              +(1|mother)+(1|block)
              ,family="binomial", data=inbcleanKS
              ,control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))

summary(model26ter)
anova(model26,model26ter)


#############################################################
##are the differences among patrilines within the montpellier or the Lyon set of lines or both?
#let us first evaluate in MOntpellier
inb2KM<-subset(inbcleanKS,inbcleanKS$paternal_pop=="Montpellier")
model21M<-glmer(y~1+(1|paternal_line)
              +(1|mother)+(1|block),data=inb2KM
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model21M)
model22M<-glmer(y~1+
              (1|mother)+(1|block),data=inb2KM
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model22M)
anova(model21M,model22M)


#are there differences within the Lyon set of lines?
inb2KL<-subset(inbcleanKS,inbcleanKS$paternal_pop=="Lyon")
model23L<-glmer(y~1+(1|paternal_line)
              +(1|mother)+(1|block),data=inb2KL
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model23L)
model24L<-glmer(y~1+
              (1|mother)+(1|block),data=inb2KL
              ,family="binomial"
              , control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model24L)
anova(model23L,model24L)

##look for patrilines with above-average male-fertility potential in KS context
##below a script to plot estimates of paternal line means with SD and 95 CI
#we work with models 39 and 41 (for Lyon and montpellier respectively)
#better to use the two separate models because we don't known whether the paternal variance ought to be the same in the two places
#start with lyon
summary(model23L)
aa<-ranef(model23L) #extract random effects
bb<-as.data.frame(aa)#convert them to data frame
bbpat<-subset(bb,bb$grpvar=="paternal_line") #keep only paternal estimates

#untransform estimates to put them into the natural scale
#we untransform also estimate+-SD and estimate+-1.95SD to plot sd and confidence intervals
#we use the estimate of the grand intercept 0.4905 from model39
bbpat$y<-antilogit(0.4905+bbpat$condval)
bbpat$yplussd<-antilogit(0.4905+bbpat$condval+bbpat$condsd)
bbpat$yminussd<-antilogit(0.4905+bbpat$condval-bbpat$condsd)
bbpat$yplus95<-antilogit(0.4905+bbpat$condval+1.95*bbpat$condsd)
bbpat$yminus95<-antilogit(0.4905+bbpat$condval-1.95*bbpat$condsd)


#a plot of paternal line means +/-1SD, with a horizontal bar for the grand mean
ggplot(bbpat, aes(x=grp, y=y)) + 
geom_pointrange(aes(ymin=yminussd, ymax=yplussd))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ geom_hline(yintercept=antilogit(0.4905))

#a plot of paternal line means with 95CI, with a horizontal bar for the grand mean
ggplot(bbpat, aes(x=grp, y=y)) + 
geom_pointrange(aes(ymin=yminus95, ymax=yplus95))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ geom_hline(yintercept=antilogit(0.4905))


#same with montpellier
summary(model21M)
aa2<-ranef(model21M) #extract random effects
bb2<-as.data.frame(aa2)#convert them to data frame
bbpat2<-subset(bb2,bb2$grpvar=="paternal_line") #keep only paternal estimates

#untransform estimates to put them into the natural scale
#we untransform also estimate+-SD and estimate+-1.95SD to plot sd and confidence intervals
#we use the estimate of the grand intercept -1.7188 from model39
bbpat2$y<-antilogit(-1.7188+bbpat2$condval)
bbpat2$yplussd<-antilogit(-1.7188+bbpat2$condval+bbpat2$condsd)
bbpat2$yminussd<-antilogit(-1.7188+bbpat2$condval-bbpat2$condsd)
bbpat2$yplus95<-antilogit(-1.7188+bbpat2$condval+1.95*bbpat2$condsd)
bbpat2$yminus95<-antilogit(-1.7188+bbpat2$condval-1.95*bbpat2$condsd)

#a plot of paternal line means +/-1SD, with a horizontal bar for the grand mean

ggplot(bbpat2, aes(x=grp, y=y)) + 
geom_pointrange(aes(ymin=yminussd, ymax=yplussd))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ geom_hline(yintercept=antilogit(-1.7188))

#a plot of paternal line means with 95CI, with a horizontal bar for the grand mean
ggplot(bbpat2, aes(x=grp, y=y)) + 
geom_pointrange(aes(ymin=yminus95, ymax=yplus95))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ geom_hline(yintercept=antilogit(-1.7188))

#plot both locations in the same graph

bbpat$paternal_pop<-"Lyon"
bbpat2$paternal_pop<-"Montpellier"

bbpat3<-rbind(bbpat2,bbpat)
str(bbpat3)
linevaluesKS<-bbpat3

figure_5b<-ggplot(bbpat3, aes(x=grp, y=y,color=paternal_pop, shape=paternal_pop))+
  guides(color=guide_legend(title=NULL), shape=guide_legend(title=NULL)) +
  ylim(0, 1)+ 
  geom_pointrange(aes(ymin=yminus95, ymax=yplus95), size = 1, show.legend = FALSE)+ 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=12))+ 
  geom_hline(aes(yintercept=antilogit(-1.7188), color = 'Montpellier'), size = 1, show.legend = FALSE) + 
  geom_hline(aes(yintercept=antilogit(0.4905), color = 'Lyon'), size = 1, show.legend = FALSE) +
  scale_color_manual(values = c( "#0072B2", "#E69F00")) +
  scale_shape_manual(values = c(16, 15))+
  labs(x= "Paternal inbred lines",
    y = "Estimated male-fertility",
    title = "KS Mothers")

figure_5b

#from the plots it seems clear that most of the lines from montpellier have no restoration potential; ie rate of fertility is around 20-30% max
#however a few of them have around 50% fertility (A22, A27, M6, a11, A18)
#In Lyon three of the lines behave as most lines from montpellier (fertility around 30% max, N80, N68, N17)
#Five of the lines restore fertility in more than half of the KI individuals (N34, N51, N67,N55, N2)
#one is "indecise" (N1)
#the results make it clear that neither maintainer and restorer genes have a 100% penetrant phenotype (ie either 100% sterile or 100% fertile)
#although some lines have close to 0 or 100% fertility
#beware however that it is an inherent property of the model to assume a gaussian distribution on the logit scale, therefore it would give an apparent continuous distribution even if the true distribution would be discontinuous

#let us draw the raw data to make some sort of a comparison
#I reorder the names of the lines to match the same order as in the plot of BLUPs

as.vector(bbpat3$grp[order(bbpat3$y)])

ordlines<- c("T3" , "T5" , "M8" , "A10" ,"A19", "A21" , "A8" , "A16","A23", "A24", "A12" ,"A20","A2" , "A7" ,"T1" , "A14", "A4",  "A13", "A9",  "A29", "A28", "A1",  "A15", "M7",  "A25", "A5",  "A17", "T7", "T10" , "M2" , "T2",  "T4",  "M1",  "T8",  "A3",  "A22" ,"A27", "M6",  "A11","A18","N80", "N68", "N17", "N1"  ,"N2" , "N55", "N67" ,"N51", "N34")
inbcleanKS$patriline<-factor(inbcleanKS$paternal_line,levels=ordlines)
inbcleanKS$g1st_binary<-ifelse(inbcleanKS$s10=="fertile","fertile","not fertile")
inbcleanKS$stbin<-factor(inbcleanKS$g1st_binary,levels=c('not fertile','fertile'))
                         
                         
                         
levels(inbcleanKS$patriline)<-ordlines

#then I plot the proportion of fertile individuals for each patriline
q<-ggplot(bbpat3, aes(x=grp, y=y,color=paternal_pop))+
theme_classic()+guides(color=guide_legend(title=NULL))+theme(plot.title = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=14),
        legend.text=element_text(size=14))+
  ylim(0, 1)+ 
geom_pointrange(aes(ymin=yminus95, ymax=yplus95))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=9))+ geom_hline(    aes(yintercept=antilogit(-1.7188),color = 'Montpellier' )  )+ geom_hline(    aes(yintercept=antilogit(0.4905),color = 'Lyon' )  )

figureS4B_b <-ggplot(inbcleanKS, aes(x=patriline, fill = stbin))+
  ylim(0, 1)+ 
  geom_bar(position = "fill")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=18))+
  scale_fill_manual(values=c("#c44601","#0073e6"))+
  labs(y = "Observed proportion of male-fertile individuals",
       x = NULL) +
  guides(fill = guide_legend(title = NULL))
  
figureS4B_b 


#relationship between line values in KS context and in DS context
linevaluesDS$grp
linevaluesKS$grp
levels(inbcleanKS$patriline)
levels(as.factor(inbcleanW$paternal_line))
listpatr<-intersect(levels(inbcleanKS$patriline),levels(as.factor(inbcleanW$paternal_line)))
jointable<-matrix(NA,nrow=length(listpatr),ncol=3)

for (i in 1:length(listpatr)) {
    jointable[i,1]<-listpatr[i]
    jointable[i,2]<-as.numeric(linevaluesDS$y[linevaluesDS$grp==listpatr[i]])
    jointable[i,3]<-linevaluesKS$y[linevaluesKS$grp==listpatr[i]]
}

latable<-as.data.frame(jointable)
names(latable)=c("paternal_line","value_in_DS","value_in_KS")
latable$value_in_DS<-as.numeric(latable$value_in_DS)
latable$value_in_KS<-as.numeric(latable$value_in_KS)
latable$paternal_pop<-ifelse(substr(latable$paternal_line, 1,1)=="N","Lyon","Montpellier")
str(latable)

#plot of predicted values 
figure5_c <-ggplot(latable, aes(x=value_in_KS,y=value_in_DS, color=paternal_pop, shape=paternal_pop))+
  ylim(0, 0.5)+
  xlim(0, 1)+ 
  theme_update(axis.title.x = element_text(size=30)) +
  geom_point(size=4)+ 
  scale_color_manual(values = c( "#0072B2", "#E69F00")) +
  scale_shape_manual(values = c(16, 15))+
  labs(
    x= "Fertility in KS Context",
    y = "Fertility in DS Context") +
  guides(
    color = guide_legend(title = NULL),
    shape = guide_legend(title = NULL)
  )

figure5_c

cor.test(latable$value_in_KS,latable$value_in_DS,method = 'pearson')
cor.test(latable$value_in_KS,latable$value_in_DS,method = 'spearman')
cor.test(latable$value_in_KS,latable$value_in_DS,method = 'kendall')


#test of maternal phenotype in the KS dataset

#how many of the individuals have a male-fertile mother
inbcleanKS$mumst_binary<-as.factor(ifelse(inbcleanKS$mats10=="fertile","fertile","not fertile"))
inb2KnoNA<-subset(inbcleanKS,is.na(inbcleanKS$mumst_binary)==FALSE)

table(inb2KnoNA$mumst_binary,inb2KnoNA$s10)

table(inb2KnoNA$mumst_binary,inb2KnoNA$s10)/as.vector(table(inb2KnoNA$mumst_binary))

#is there a difference in the probability to be fertile depending on whether the mother is sterile versus fertile
model91k<-glmer(y~1+ mumst_binary + paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block),
              data=inb2KnoNA
              ,family="binomial"
              ,control= glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model91k)
model92k<-glmer(y~1+paternal_pop
              +(1|paternal_line)+(1|mother)+(1|block),
              data=inb2KnoNA
              ,family="binomial"
              ,control= glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
summary(model92k)
anova(model91k,model92k)
#the effect of maternal sterility is not significant 


#plot fertility data within the KS category, as a function of paternal pop and maternal phenotype
inb2KnoNA$pat_matfert<-paste(inb2KnoNA$paternal_pop,inb2KnoNA$mumst_binary)
table(inb2KnoNA$pat_matfert,inb2KnoNA$s10)

inb2KnoNA$pat_matfert <- gsub("Montpellier fertile", "Montpellier\nFertile", inb2KnoNA$pat_matfert)
inb2KnoNA$pat_matfert <- gsub("Montpellier not fertile", "Montpellier\nnot fertile", inb2KnoNA$pat_matfert)
inb2KnoNA$pat_matfert <- gsub("Lyon not fertile", "Lyon\nnot fertile", inb2KnoNA$pat_matfert)
inb2KnoNA$pat_matfert <- gsub("Lyon fertile", "Lyon\nfertile", inb2KnoNA$pat_matfert)

figureS4A_b <- ggplot(inb2KnoNA,aes(x = pat_matfert,fill = s10)) +
    geom_bar(position = "fill") + 
  scale_fill_manual(values=c("#0073e6", "#c44601", "#f57600")) +
  guides(fill = guide_legend(title = NULL)) +
  labs(x = NULL,
       y = NULL,
       title = "Maternale type = KS") +
  theme(axis.text.x = element_text(size = 25))
figureS4A_b

#### for figure 5 

figure5_final<- figure_5a_95+ figure_5b + figure5_c +
  plot_layout(guides = "collect")+
  plot_annotation(tag_levels = 'A') &
  theme(plot.tag = element_text(size = 45),
        legend.position = 'bottom') 

figure5_final

ggsave("figure5.png", figure5_final, width = 25, height = 10)

#### figure S4

figureS4A_final<- figureS4A_a + figureS4A_b +
  plot_layout(guides = "collect") &
  theme(legend.position = 'bottom') 

figureS4A_final
ggsave("figureS4A_final.png", figureS4A_final, width = 16, height = 10, dpi = 600)


figure_5bbis<-ggplot(bbpat3, aes(x=grp, y=y,color=paternal_pop, shape=paternal_pop))+
  guides(color=guide_legend(title=NULL), shape=guide_legend(title=NULL)) +
  ylim(0, 1)+ 
  geom_pointrange(aes(ymin=yminus95, ymax=yplus95), size = 1)+ 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 15))+ 
  geom_hline(aes(yintercept=antilogit(-1.7188), color = 'Montpellier'), size = 1, show.legend = FALSE) + 
  geom_hline(aes(yintercept=antilogit(0.4905), color = 'Lyon'), size = 1, show.legend = FALSE) +
  scale_color_manual(values = c( "#0072B2", "#E69F00")) +
  scale_shape_manual(values = c(16, 15))+
  labs(
    x= NULL,
    y = "Estimated proportion of\nmale-fertile individuals",
    title = "KS Mothers")

figure_5bbis

figureS4B_final<- figure_5bbis / figureS4B_b +
  plot_layout(guides = "collect") &
  theme(legend.position = 'bottom') 

figureS4B_final
ggsave("figureS4B_final.png", figureS4B_final, width = 12, height = 16, dpi = 600)

##### 

```



