mark_prop.R

#' @title Exclude markers depending on proportions of homo/hetorozygous
#'
#' @description This function uses the dataframe produced by the tab_mark function and fills the "exclude" column for all the markers that present odd proportions of each genotype. You can define these proportions thanks to the arguments of the function.
#' @param tab data frame obtained with tab_mark function.
#' @param cross F2 or N2.
#' @param homo proportion of homozygous individuals under which the marker is excluded. Will apply on both homozygous genotypes for a F2, but only on one for N2.
#' @param hetero proportion of heterozygous individuals under which the marker is excluded.
#' @param na proportion of non-genotyped individuals above which the marker is excluded.
#'
#' @import dplyr
#' @import tidyselect
#'
#' @export
#'

#### mark_prop ####
## excludes markers depending on proportions of homo/hetorozygous
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
  #calculate total number of individuals genotyped for each marker
  tab <- tab %>% mutate(n_geno = (n_HM1 + n_HM2 + n_HT))


  #stop if cross != "F2" or "N2"
  if(!cross %in% c("F2","N2")){
    stop("Cross must be F2 or N2")
  }

  #stop of homo&hetero or pval not specified
  if((is.na(homo)==TRUE | is.na(hetero)==TRUE) & is.na(pval)==TRUE){
    stop("Arguments homo and hetero or argument pval must be specified")
  }



  #stop with prop of na
  #calculate proportion
  tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))

  tab <- tab %>%
    mutate(exclude_prop=case_when(p_NA > na ~ 1,
                                  T ~ 0))

  #stock colnames to join
  names <- colnames(tab)

  #exclude with prop of homo/hetero
  if(is.na(pval)==TRUE){
    #calculate proportion of each genotype
    tab <- tab %>% mutate(p_HM1 = n_HM1/n_geno)
    tab <- tab %>% mutate(p_HM2 = n_HM2/n_geno)
    tab <- tab %>% mutate(p_HT = n_HT/n_geno)


    tab <- tab %>%
      mutate(exclude_prop=case_when(p_NA > na ~ 1,
                                    cross=="F2" & (p_HM1 < homo | p_HM2 < homo | p_HT < hetero) ~ 1,
                                    cross=="N2" & (p_HM1 < homo | p_HT < hetero) ~ 1,
                                    T ~ exclude_prop
      ))

  #exclude with pval chisq.test
  } else if(is.na(pval)==FALSE){

    #if cross F2
    if(cross=="F2"){
      tab <- tab %>% filter(p_NA != 1) %>% rowwise() %>%
        mutate(.,chi_pval = tibble(n_HM1,n_HM2,n_HT) %>%
                 chisq.test(p=c(0.25,0.25,0.5)) %>% .$p.value) %>%
        full_join(.,tab,by=all_of(names))

      tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
                                                   T ~ exclude_prop))

      #if cross N2
    } else if(cross=="N2"){
      tab <- tab %>% filter(p_NA != 1) %>% rowwise() %>%
        mutate(.,chi_pval = tibble(n_HM1,n_HT) %>%
                 chisq.test(p=c(0.5,0.5)) %>% .$p.value) %>%
        full_join(.,tab,by=all_of(names))

      tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
                                                   T ~ exclude_prop))
    }


  }
  tab <- tab %>% select(all_of(names)) %>% select(-c(n_geno,p_NA))
  return(tab)
}