mark_prop.R 3.76 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
    tab <- tab %>%
Marie Bourdon's avatar
Marie Bourdon committed
55
56
57
      mutate(exclude_prop=case_when(!chr %in% c("X","Y","M") & p_NA > na ~ 1,
                                    !chr %in% c("X","Y","M") & cross=="F2" & (p_HM1 < homo | p_HM2 < homo | p_HT < hetero) ~ 1,
                                    !chr %in% c("X","Y","M") & cross=="N2" & ((p_HM2 == 0 & p_HM1 < homo) |
Marie Bourdon's avatar
Marie Bourdon committed
58
59
60
                                                     (p_HM1 == 0 & p_HM2 < homo) |
                                                     (p_HT < hetero) |
                                                     (p_HM1 !=0 & p_HM2 != 0)) ~ 1,
mariefbourdon's avatar
mariefbourdon committed
61
62
                                    T ~ exclude_prop
      ))
Marie  BOURDON's avatar
Marie BOURDON committed
63

Marie Bourdon's avatar
Marie Bourdon committed
64
    #exclude with pval chisq.test
Marie Bourdon's avatar
Marie Bourdon committed
65
66
67
68
  } else if(is.na(pval)==FALSE){

    #if cross F2
    if(cross=="F2"){
Marie Bourdon's avatar
Marie Bourdon committed
69
      #chisq test
Marie Bourdon's avatar
Marie Bourdon committed
70
71
72
73
74
      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))

Marie Bourdon's avatar
Marie Bourdon committed
75
      #convert in score
Marie Bourdon's avatar
Marie Bourdon committed
76
77
78
79
80
      tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
                                                   T ~ exclude_prop))

      #if cross N2
    } else if(cross=="N2"){
Marie Bourdon's avatar
Marie Bourdon committed
81
82
83
84
85
86

      #stock homozygous data in a new column (either n_HM1 or n-<hm2)
      tab <- tab %>% mutate(n_hm = case_when(n_HM1>n_HM2 ~ n_HM1,
                                             T ~ n_HM2))

      #chisq test
Marie Bourdon's avatar
Marie Bourdon committed
87
      tab <- tab %>% filter(p_NA != 1) %>% rowwise() %>%
Marie Bourdon's avatar
Marie Bourdon committed
88
89
        mutate(.,chi_pval =chisq.test(x=c(n_hm,n_HT),p=c(0.5,0.5)) %>% .$p.value) %>%
        full_join(.,tab,by=c(all_of(names),"n_hm"))
Marie Bourdon's avatar
Marie Bourdon committed
90

Marie Bourdon's avatar
Marie Bourdon committed
91
      #convert in score
Marie Bourdon's avatar
Marie Bourdon committed
92
      tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
Marie Bourdon's avatar
Marie Bourdon committed
93
                                                   (n_HM1 !=0 & n_HM2 != 0) ~ 1,
Marie Bourdon's avatar
Marie Bourdon committed
94
95
                                                   T ~ exclude_prop))
    }
mariefbourdon's avatar
mariefbourdon committed
96
97
98


  }
Marie Bourdon's avatar
Marie Bourdon committed
99
100
  tab <- tab %>% select(all_of(names)) %>% select(-c(n_geno,p_NA))
  return(tab)
Marie  BOURDON's avatar
Marie BOURDON committed
101
}