From 7f522196e02480aaac8aa2cb0605ea59b69297d8 Mon Sep 17 00:00:00 2001 From: Marie Bourdon <mbourdon@pasteur.fr> Date: Wed, 23 Jun 2021 10:43:14 +0200 Subject: [PATCH] mark_prop chisq backcross --- R/mark_prop.R | 40 +++++++++++++++++++++++++++------------- vignettes/stuaRt.Rmd | 2 ++ 2 files changed, 29 insertions(+), 13 deletions(-) diff --git a/R/mark_prop.R b/R/mark_prop.R index 3b0ff5d..a7b2508 100755 --- a/R/mark_prop.R +++ b/R/mark_prop.R @@ -15,9 +15,6 @@ #### mark_prop #### ## excludes markers depending on proportions of homo/hetorozygous mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){ - #stock colnames to join - names <- colnames(tab) - #calculate total number of individuals genotyped for each marker tab <- tab %>% mutate(n_geno = (n_HM1 + n_HM2 + n_HT)) @@ -42,6 +39,9 @@ mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){ 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 @@ -56,19 +56,33 @@ mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){ cross=="N2" & (p_HM1 < homo | p_HT < hetero) ~ 1, T ~ exclude_prop )) - } #exclude with pval chisq.test - ## NEED TO ADD THIS FILTER IF CROSS = N2 - if(is.na(pval)==FALSE){ + } 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 %>% 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)) } - tab <- tab %>% select(all_of(names),exclude_prop) + tab <- tab %>% select(all_of(names)) %>% select(-c(n_geno,p_NA)) + return(tab) } diff --git a/vignettes/stuaRt.Rmd b/vignettes/stuaRt.Rmd index 48568f2..d156d33 100755 --- a/vignettes/stuaRt.Rmd +++ b/vignettes/stuaRt.Rmd @@ -15,6 +15,8 @@ knitr::opts_chunk$set( ) ``` +# stuart + Marie Bourdon June 2021 -- GitLab