In this study we explored the mode of genetic inheritance of Varroa destructor, a parasitic mite of honeybees, so far known as a haplodiploid species with a arrhenotokous parthenogenesis reproductive mode. The mite shows contrasting phenomena of highly inbreeding life style (sib-mating), yet maintaining relatively high genetic variation.

The experimental setup consisted of a three-generational pedigree.
Sample size: 30 families, total of 223 individuals.
3 different colonies, from OIST apiary.
For details of mite collection and pedigree construction, please see the original manuscript.
All biosamples are available in Sequence Read Archive (SRA) under the accession PRJNA794941.
All data are available and reproducible from the GitHub page.

load libraries

library("tidyverse")
library("plyr")
library("dplyr")
library("ggplot2")
library("scales")
library("ggpubr")
library("gridExtra")
library("grid")
library("GGally")
library("vcfR") # for extracting genotype data from a vcf file
library("data.table")
library("stringr")
library("janitor")
library("gmodels")
library("rstatix")
library("freqtables")
library("broom")
library("cowplot") 

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE,
                     fig.width = 10,
                      fig.asp = 0.8,
                      out.width = "100%")
#fig.width = 6,fig.asp = 0.8,out.width = "100%"

the VCF was filtered using vcftools with the following filters:

vcftools --vcf snponly_freebayes.vcf --chr NW_019211454.1 --chr NW_019211455.1 --chr NW_019211456.1 --chr NW_019211457.1 --chr NW_019211458.1 --chr NW_019211459.1 --chr NW_019211460.1 --max-alleles 2 --minQ 15000 --minDP 16 --maxDP 40 --max-missing 0.5 --maf 0.2 --recode --recode-INFO-all --out Q15000BIALLDP16HDP40mis.5maf.2Chr7

Load Variant Call Format (VCF) file.

Extract genotypes for each site and individual. The metadata for all samples can be found in here.

vcf <- read.vcfR("/Users/nuriteliash/Documents/GitHub/varroa-linkage-map/data/vcf_filter/minQ_filter/Q15000BIALLDP16HDP40mis.5maf.2Chr7.recode.vcf", verbose = FALSE )
vcf
## ***** Object of Class vcfR *****
## 223 samples
## 7 CHROMs
## 33,925 variants
## Object size: 272.8 Mb
## 0 percent missing data
## *****        *****         *****
# extract the genotype for each site in each individual
gt <- extract.gt(vcf, element = "GT") 
gt <- as.data.frame(t(gt)) %>%
  rownames_to_column("ID")
#clean the ID names
gt$ID <- sub("_[^_]+$", "", gt$ID)

table <-  gt %>% 
  t() %>%
  as.data.frame() %>%
  row_to_names(row_number = 1) %>% 
  dplyr::select(contains(c("son", "dat", "fnd"))) # keep only adults of F0, F1 and F2 

# set the families (include only families with at least one adult F2)
family = grep("grndat|grnson",gt$ID, value=TRUE) %>%
  str_extract("[^_]+")  %>%
  unique()

# or, include all F2 samples, but indicate if they have an adult sister (may indicate if the F1 female was fertilized)
#family = grep("grn",gt$ID, value=TRUE) %>%
#  str_extract("[^_]+")  %>%
#  unique()

The individual ID name nomenclature is composed of 3 parts, separated by an underscore (_).
The first part is the family ID, the second is the unique name of the individual, and the the third part indicates its generation, sex, and its position in relation to the foundress mite (generation F0).

For example:
Individual ID 240_240a_son, belong to family 240,
has a unique name of 240a,
and is the son of the foundress mite of family 240, that is, its a male of F1 generation.

Individual ID 240_241b_grndat, belong to family 240,
has a unique name of 241b,
and is the granddaughter of the foundress mite of family 240, that is, its a female of F2 generation (and the daughter of 240_240a_son).

From the 223 individuals of 30 families, we excluded non adults mites (for which sex could not be determined for sure).
Eventually, we kept 144 individuals (adult males and females in F0, F1 and F2 generations) of 26 families, and all 33,925 biallelic sites.

For example, these are the first 6 sites (13565-25513) of chromosome (NW_019211454.1), in family number 240.

table %>% select(contains(c("240_"))) %>% head()
##                      240_240a_son 240_241c_grnson 240_241b_grndat
## NW_019211454.1_18037          0/0            <NA>            <NA>
## NW_019211454.1_20998         <NA>            <NA>            <NA>
## NW_019211454.1_24935         <NA>             0/0            <NA>
## NW_019211454.1_25470         <NA>            <NA>            <NA>
## NW_019211454.1_25513         <NA>            <NA>            <NA>
## NW_019211454.1_35412         <NA>            <NA>            <NA>
##                      240_241a_grndat 240_241_dat 240_240_fnd
## NW_019211454.1_18037             0/0         0/0         0/0
## NW_019211454.1_20998             0/0         0/0         0/0
## NW_019211454.1_24935             0/0         0/0         0/0
## NW_019211454.1_25470             0/0         0/0         0/0
## NW_019211454.1_25513             0/0        <NA>         0/0
## NW_019211454.1_35412            <NA>         0/0         0/0

The family members include:

  • F0 generation: foundress female mite (240_240_fnd).
  • F1 generation: adult son (240_240a_son) and adult daughter (240_241_dat) of the foundress F0 mite. Because in varroa mite reproduction is via sib-mating, these brother and sister will also mate, to produce the F2 generation.
  • F2 generation: adult grandson (240_241c_grnson) and two adult granddaughters (240_241a_grndat and 240_241b_grndat), of the foundress F0 mite.

Each individual, is genotyped for each site with one of the three genotypes:

  • Homozygote for the reference allele = 0/0
  • Heterozygote = 0/1
  • Homozygote for the alternate allele = 1/1
  • “NA” = site genotype not determined

For more information about the mapping, variant calling and variant filtration, please see the Snakemake pipeline in here.

In the following code we viewed the F2 generation genotype of all possible nine crosses between F1 male and females.

  • Control crosses of homozygotic sites: The first four crosses were aimed mainly to detect false sites:
    1. F1 male (0/0) x female (0/0)
    2. F1 male (1/1) x female (1/1)
    3. F1 male (0/0) x female (1/1)
    4. F1 male (1/1) x female (0/0)

Then we crossed heterozygotic sites, to explore the mode of genetic inheritance in varroa mite:

  • Homozygotic male, with heterozygotic female:
    1. F1 male (0/0) x female (0/1)
    2. F1 male (1/1) x female (0/1)
  • Heterozygotic male, with homozygotic female:
    1. F1 male (0/1) x female (0/0)
    2. F1 male (0/1) x female (1/1)
  • Heterozygotic male, with heterozygotic female:
    1. F1 male (0/1) x female (0/1)

In addition to the fixed F1 genotypes, we also “fixed” the F0 female genotype to match that of her son (F1 male), so we can phase the alleles of the F1 generations.
“Phasing alleles” is the process of determine the parental origin of each allele.
for heterozygotic genotype phasing is critical, as it allows the tracking of the allele from one generation tot he next.

Abnormal and Normal males

While most F2 males exhibit a parthenogenetic mode of inheritance, a few males carry also paternal genotypes.
Plotting the distribution of F2 males (proportion of maternal sites), we set the threshold for “normal” males at 0.75. That is in “normal” males, at least 75% of the genotypes are maternal. The rest of the males, are considered “abnormal”, and were not included in the final analysis.

The “abnormal” males are: 412_413a_grnson, 400_401a_grnson, 458_459a_grnson, 240_241c_grnson, 426_427b_grnson and 300_301a_grnson

# define the group of abnormal males
abnorm_males = tibble(sample = c("240_241c_grnson", "400_401a_grnson","412_413a_grnson", "426_427b_grnson", "458_459a_grnson", "46_47d_grnson"), sex = "male", normality = "abnormal")

‘Sanity check’: Crossing homozygotic sites (crosses 1-4)

(1) F1 male (0/0) x female (0/0)

we expect all F2 offspring to be homozygotic (0/0) like their parents (F1)

obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/0")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/0")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/0")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) %>%
    left_join(abnorm_males, by = "sample") %>% 
    dplyr::select(-sex.y) %>% 
    replace(is.na(.), "normal") %>%
    dplyr::rename(sex = sex.x) %>%
    unite("type", sex,normality, remove = FALSE)
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/0") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
   mutate(count_sample = n_distinct(sample)) %>%
   select(c(sex,sites,count_sample)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 3
## # Groups:   sex [2]
##   sex    sites count_sample
##   <chr>  <dbl>        <int>
## 1 male   35620           26
## 2 female 35003           25
p_00_00 =  samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
          axis.title.x = element_blank()) +
      ylab("F2 genotype proportion") +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 cross: male (0/0) x female (0/0)") +
  scale_x_discrete(breaks=c("female","male"),
        labels=c(paste(filter(pooled_obs_count, sex == 'female')$count_sample, " females,\n",
filter(pooled_obs_count, sex == 'female')$sites, "sites"), 
paste(filter(pooled_obs_count, sex == 'male')$count_sample, " males,\n",
filter(pooled_obs_count, sex == 'male')$sites, "sites")))

(2) F1 male (1/1) x female (1/1)

we expect all F2 offspring to be homozygotic (1/1) like their parents (F1)

# define a list to put all the data frames in
# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "1/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "1/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "1/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="1/1") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male   24520
## 2 female 23344
p_11_11 =  samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (1/1) x F1 female (1/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

(3) F1 male (0/0) x female (1/1)

We expect zero sites for this cross, because the F1 are siblings.
Indeed, there are only 23 sites in the F2 pooled females, and 14 for the F2 pooled males:

# define a list to put all the data frames in
# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/0")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/0")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "1/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/0") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 0 × 2
## # Groups:   sex [0]
## # … with 2 variables: sex <chr>, sites <dbl>
p_00_11 =  samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (0/0) x F1 female (1/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

(4) F1 male (1/1) x female (0/0)

We expect zero sites for this cross, since there is no paternal inheritance to the males.

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "1/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "1/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/0")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="1/1") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 0 × 2
## # Groups:   sex [0]
## # … with 2 variables: sex <chr>, sites <dbl>
p_11_00 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (1/1) x F1 female (0/0)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

‘Informative crosses’: at least one of the parents is heterozygous (crosses 5-9)

Crossing homozygotic male, with heterozygotic female

(5) F1 male (0/0) x female (0/1)

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/0")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/0")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/0") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male    5691
## 2 female  4942
p_00_01 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (0/0) x F1 female (0/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

(6) F1 male (1/1) x female (0/1)

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "1/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "1/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/1") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male    6101
## 2 female  5389
p_11_01 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (1/1) x F1 female (0/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

Crossing heterozygotic male, with homozygotic female

The former crosses of heterozygotic females (5 and 6) show that F2 males can be heterozygotic and carry two alleles, like their mother.
But are these sites real? and if they are, can they be transmitted to their daughters?

(7) F1 male (0/1) x female (0/0)

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/0")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/0") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male    1978
## 2 female  2465
p_01_00 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (0/1) x F1 female (0/0)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

(8) F1 male (0/1) x female (1/1)

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "1/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/1") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male    1573
## 2 female  2391
p_01_11 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (0/1) x F1 female (1/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

Crossing heterozygotic male and female

(9) F1 male (0/1) x female (0/1)

# define a list to put all the data frames in
obs <- list()

for (fam in family) {
 
obs[[fam]] <- table %>%
  dplyr::select(starts_with(fam)) %>%
    dplyr::filter_at(vars(matches("_fnd")), all_vars(. == "0/1")) %>% # force F0 female to be homo, like her son
  dplyr::filter_at(vars(matches("_son")), all_vars(. == "0/1")) %>%
  dplyr::filter_at(vars(matches("_dat")), all_vars(. == "0/1")) %>% 
  dplyr::select(contains("grn")) %>%
  tidyr::pivot_longer(everything())  %>% 
  dplyr::rename(sample = name, gt = value) %>%
  dplyr::count(sample, gt, .drop = FALSE) %>% 
  dplyr::filter(gt %in% c("0/0", "1/1", "0/1")) %>%
  mutate(n = as.numeric(n)) %>%
  group_by(sample) %>%
  mutate(total = as.numeric(sum(n))) %>%
  mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female"))
}

# bind all families together, to a final data frame containing all observed counts
observed <- do.call("rbind", obs) %>% mutate(sample = as.character(sample))
samples <- data.frame(sample = rep(unique(observed$sample), each =3), gt = rep(c("0/0", "0/1", "1/1")))

samples_obs <- left_join(samples, observed, by=c("sample","gt")) %>% mutate(sex = case_when(
    grepl("son", sample) ~ "male",
    grepl("dat", sample) ~ "female")) %>%
    group_by(sample) %>%
    replace(is.na(.), 0) %>%
    mutate(total = as.numeric(sum(n))) %>%
    dplyr::mutate(prop = n/total) 
  
pooled_obs_count =  samples_obs %>% 
    dplyr::filter(total>=10) %>%
    filter(gt =="0/1") %>%
    group_by(sex) %>% 
    mutate(sites = sum(total)) %>% 
    select(c(sex,sites)) %>% 
    unique()

pooled_obs_count
## # A tibble: 2 × 2
## # Groups:   sex [2]
##   sex    sites
##   <chr>  <dbl>
## 1 male    3657
## 2 female  4618
p_01_01 = samples_obs %>% dplyr::filter(total>=10) %>%
    ggplot(aes(fill=gt, y=prop, x=sex)) + 
    geom_bar(position="fill", stat="identity", ) +
    labs(fill = "Genotype") + 
    theme_classic() +
    theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
     scale_fill_manual(values=c("#ffbf00", "#66b032","#1982c4"))+
     theme(legend.position = "none")+
  ggtitle("F1 male (0/1) x F1 female (0/1)",
          subtitle = paste("F2 female = ", filter(pooled_obs_count, sex == 'female')$sites, ", F2 male =", filter(pooled_obs_count, sex == 'male')$sites))

Plot sanity check crosses

legend <- get_legend(p_00_00)   # get the legend of the first one plot

# here the plots in a grid
prow <- plot_grid( p_00_00 + theme(legend.position="none"),
           # here you add the percentage
           p_11_11 + theme(legend.position="none")+ scale_y_continuous(),
           p_11_00 + theme(legend.position="none")+ scale_y_continuous(),
          p_00_11 + theme(legend.position="none")+ scale_y_continuous(),
         align = 'v',
           #labels = c("A", "B"),
           hjust = -1,
           nrow = 2)

# here you add the legend
plot_grid( prow, legend, rel_widths = c(1, .2))

Plot informative crosses

legend <- get_legend(p_01_01)   # get the legend of the first one plot

# here the plots in a grid
prow <- plot_grid( p_00_01 + theme(legend.position="none"),
           # here you add the percentage
           p_11_01 + theme(legend.position="none")+ scale_y_continuous(),
           p_01_00 + theme(legend.position="none")+ scale_y_continuous(),
          p_01_11 + theme(legend.position="none")+ scale_y_continuous(),
         align = 'v',
           #labels = c("A", "B"),
           hjust = -1,
           nrow = 2)

# here you add the legend
plot_grid( prow, legend, rel_widths = c(1, .2))