mark_prop.R 3.12 KB
Newer Older
Marie  BOURDON's avatar
Marie BOURDON committed
1
2
#' @title Exclude markers depending on proportions of homo/hetorozygous
#'
Marie Bourdon's avatar
Marie Bourdon committed
3
#' @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.
Marie Bourdon's avatar
Marie Bourdon committed
4
5
6
7
8
#' @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.
Marie  BOURDON's avatar
Marie BOURDON committed
9
10
#'
#' @import dplyr
Marie Bourdon's avatar
Marie Bourdon committed
11
#' @import tidyselect
Marie  BOURDON's avatar
Marie BOURDON committed
12
13
14
15
16
17
#'
#' @export
#'

#### mark_prop ####
## excludes markers depending on proportions of homo/hetorozygous
mariefbourdon's avatar
mariefbourdon committed
18
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
Marie Bourdon's avatar
Marie Bourdon committed
19
20
  #calculate total number of individuals genotyped for each marker
  tab <- tab %>% mutate(n_geno = (n_HM1 + n_HM2 + n_HT))
mariefbourdon's avatar
mariefbourdon committed
21
22


Marie  BOURDON's avatar
Marie BOURDON committed
23
24
25
26
27
  #stop if cross != "F2" or "N2"
  if(!cross %in% c("F2","N2")){
    stop("Cross must be F2 or N2")
  }

mariefbourdon's avatar
mariefbourdon committed
28
29
30
  #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")
Marie  BOURDON's avatar
Marie BOURDON committed
31
32
  }

mariefbourdon's avatar
mariefbourdon committed
33
34
35


  #stop with prop of na
Marie Bourdon's avatar
Marie Bourdon committed
36
37
38
  #calculate proportion
  tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))

mariefbourdon's avatar
mariefbourdon committed
39
40
41
42
  tab <- tab %>%
    mutate(exclude_prop=case_when(p_NA > na ~ 1,
                                  T ~ 0))

Marie Bourdon's avatar
Marie Bourdon committed
43
44
45
  #stock colnames to join
  names <- colnames(tab)

Marie Bourdon's avatar
Marie Bourdon committed
46
  #exclude with prop of homo/hetero
mariefbourdon's avatar
mariefbourdon committed
47
  if(is.na(pval)==TRUE){
Marie Bourdon's avatar
Marie Bourdon committed
48
49
50
51
52
53
    #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)


mariefbourdon's avatar
mariefbourdon committed
54
55
56
57
58
59
    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
      ))
Marie  BOURDON's avatar
Marie BOURDON committed
60

Marie Bourdon's avatar
Marie Bourdon committed
61
  #exclude with pval chisq.test
Marie Bourdon's avatar
Marie Bourdon committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
  } 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))
    }
mariefbourdon's avatar
mariefbourdon committed
84
85
86


  }
Marie Bourdon's avatar
Marie Bourdon committed
87
88
  tab <- tab %>% select(all_of(names)) %>% select(-c(n_geno,p_NA))
  return(tab)
Marie  BOURDON's avatar
Marie BOURDON committed
89
}