.Rhistory 17.4 KB
Newer Older
Marie Bourdon's avatar
Marie Bourdon committed
1
2
3
4
5
6
7
8
9
10
devtools::build_manual(path=".")
remove.packages("stuart")
devtools::install_gitlab(repo="mouselab/stuart",host="gitlab.pasteur.fr")
library(stuart)
View(annot_mini)
annot_mini %>% filter(is.na(chr)==FALSE) %>% filter(is.na(cM_cox)==TRUE)
devtools::build(path=".",vignettes=FALSE)
library(stuart)
write_rqtl
chisq.test(c(20,23),p=c(0.5,0.5))
mariefbourdon's avatar
mariefbourdon committed
11
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
12
13
#stock colnames to join
names <- colnames(tab)
Marie Bourdon's avatar
Marie Bourdon committed
14
15
#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
16
17
18
19
20
21
22
23
24
#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
25
26
#calculate proportion
tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))
mariefbourdon's avatar
mariefbourdon committed
27
28
29
tab <- tab %>%
mutate(exclude_prop=case_when(p_NA > na ~ 1,
T ~ 0))
30
#exclude with prop of homo/hetero
mariefbourdon's avatar
mariefbourdon committed
31
if(is.na(pval)==TRUE){
32
33
34
35
#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
36
37
38
39
40
41
42
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
))
}
43
44
#exclude with pval chisq.test
## NEED TO ADD THIS FILTER IF CROSS = N2
mariefbourdon's avatar
mariefbourdon committed
45
if(is.na(pval)==FALSE){
Marie Bourdon's avatar
Marie Bourdon committed
46
47
#if cross F2
if(cross=="F2"){
mariefbourdon's avatar
mariefbourdon committed
48
49
50
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) %>%
51
full_join(.,tab,by=all_of(names))
mariefbourdon's avatar
mariefbourdon committed
52
53
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
Marie Bourdon's avatar
Marie Bourdon committed
54
55
56
57
58
59
60
61
} else if(cross=="N2"){
ab <- 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
62
}
63
tab <- tab %>% select(all_of(names),exclude_prop)
Marie Bourdon's avatar
Marie Bourdon committed
64
return(tab)
65
}
Marie Bourdon's avatar
Marie Bourdon committed
66
67
68
69
70
71
72
73
test <- tibble(marker=c("mark1","mark2","mark3","mark4"),)
library(tidyverse)
test <- tibble(marker=c("mark1","mark2","mark3","mark4"),allele_1=c(A,G,T,C),allele_2=c(T,T,A,G),n_HM1=c(48,45,1,26),n_HM2=c(0,0,0,0),n_HT=c(52,55,49,74))
test <- tibble(marker=c("mark1","mark2","mark3","mark4"),allele_1=c("A","G","T","C"),allele_2=c("T","T","A","G"),n_HM1=c(48,45,1,26),n_HM2=c(0,0,0,0),n_HT=c(52,55,49,74))
View(test)
mark_prop(test,cross="N2",pval=0.05)
test <- tibble(marker=c("mark1","mark2","mark3","mark4"),allele_1=c("A","G","T","C"),allele_2=c("T","T","A","G"),n_HM1=c(48,45,1,26),n_HM2=c(0,0,0,0),n_HT=c(52,55,49,74),n_NA=c(0,0,0,0))
mark_prop(test,cross="N2",pval=0.05)
mariefbourdon's avatar
mariefbourdon committed
74
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
75
76
#stock colnames to join
names <- colnames(tab)
mariefbourdon's avatar
mariefbourdon committed
77
#calculate total number of individuals genotyped for each marker
mariefbourdon's avatar
mariefbourdon committed
78
tab <- tab %>% mutate(n_geno = (n_HM1 + n_HM2 + n_HT))
mariefbourdon's avatar
mariefbourdon committed
79
80
81
82
83
84
85
86
87
#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
88
89
#calculate proportion
tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))
mariefbourdon's avatar
mariefbourdon committed
90
91
92
tab <- tab %>%
mutate(exclude_prop=case_when(p_NA > na ~ 1,
T ~ 0))
93
#exclude with prop of homo/hetero
mariefbourdon's avatar
mariefbourdon committed
94
if(is.na(pval)==TRUE){
95
96
97
98
#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
99
100
101
102
103
104
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#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),exclude_prop)
return(tab)
}
mark_prop(test,cross="N2",pval=0.05)
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))
#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")
mariefbourdon's avatar
mariefbourdon committed
141
}
Marie Bourdon's avatar
Marie Bourdon committed
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#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))
print(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
))
161
#exclude with pval chisq.test
Marie Bourdon's avatar
Marie Bourdon committed
162
163
164
} else if(is.na(pval)==FALSE){
#if cross F2
if(cross=="F2"){
mariefbourdon's avatar
mariefbourdon committed
165
166
167
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) %>%
168
full_join(.,tab,by=all_of(names))
mariefbourdon's avatar
mariefbourdon committed
169
170
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
Marie Bourdon's avatar
Marie Bourdon committed
171
172
173
174
175
176
177
178
179
#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
180
}
Marie Bourdon's avatar
Marie Bourdon committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
tab <- tab %>% select(all_of(names),exclude_prop)
return(tab)
}
mark_prop(test,cross="N2",pval=0.05)
mark_prop(test,cross="N2",homo=0.1,hetero=0.1)
mark_prop(test,cross="F2",homo=0.1,hetero=0.1)
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))
#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))
206
print(tab)
Marie Bourdon's avatar
Marie Bourdon committed
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#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))
print(tab)
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
}
}
240
tab <- tab %>% select(all_of(names),exclude_prop)
Marie Bourdon's avatar
Marie Bourdon committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
return(tab)
}
mark_prop(test,cross="F2",pval=0.05)
mark_prop(test,cross="N2",pval=0.05)
View(test)
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))
#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))
265
print(tab)
Marie Bourdon's avatar
Marie Bourdon committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
#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
))
print(names)
#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))
print(tab)
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
298
}
Marie Bourdon's avatar
Marie Bourdon committed
299
300
301
302
303
}
tab <- tab %>% select(all_of(names),exclude_prop)
return(tab)
}
mark_prop(test,cross="N2",pval=0.05)
mariefbourdon's avatar
mariefbourdon committed
304
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
305
306
307
#stock colnames to join
names <- colnames(tab)
print(names)
mariefbourdon's avatar
mariefbourdon committed
308
#calculate total number of individuals genotyped for each marker
mariefbourdon's avatar
mariefbourdon committed
309
tab <- tab %>% mutate(n_geno = (n_HM1 + n_HM2 + n_HT))
mariefbourdon's avatar
mariefbourdon committed
310
311
312
313
314
315
316
317
318
#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
319
320
#calculate proportion
tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))
mariefbourdon's avatar
mariefbourdon committed
321
322
323
tab <- tab %>%
mutate(exclude_prop=case_when(p_NA > na ~ 1,
T ~ 0))
Marie Bourdon's avatar
Marie Bourdon committed
324
print(tab)
325
#exclude with prop of homo/hetero
mariefbourdon's avatar
mariefbourdon committed
326
if(is.na(pval)==TRUE){
327
328
329
330
#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
331
332
333
334
335
336
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
))
337
#exclude with pval chisq.test
Marie Bourdon's avatar
Marie Bourdon committed
338
339
340
} else if(is.na(pval)==FALSE){
#if cross F2
if(cross=="F2"){
mariefbourdon's avatar
mariefbourdon committed
341
342
343
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) %>%
344
full_join(.,tab,by=all_of(names))
mariefbourdon's avatar
mariefbourdon committed
345
346
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
Marie Bourdon's avatar
Marie Bourdon committed
347
348
349
350
351
352
#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))
353
print(tab)
Marie Bourdon's avatar
Marie Bourdon committed
354
355
356
357
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
}
}
358
tab <- tab %>% select(all_of(names),exclude_prop)
Marie Bourdon's avatar
Marie Bourdon committed
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
return(tab)
}
mark_prop(test,cross="N2",pval=0.05)
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))
print(tab)
#stock colnames to join
names <- colnames(tab)
print(names)
#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))
411
print(tab)
Marie Bourdon's avatar
Marie Bourdon committed
412
413
tab <- tab %>% mutate(exclude_prop=case_when(chi_pval < pval ~ 1,
T ~ exclude_prop))
mariefbourdon's avatar
mariefbourdon committed
414
}
Marie Bourdon's avatar
Marie Bourdon committed
415
}
Marie Bourdon's avatar
Marie Bourdon committed
416
417
tab <- tab %>% select(all_of(names),exclude_prop)
return(tab)
Marie Bourdon's avatar
Marie Bourdon committed
418
}
Marie Bourdon's avatar
Marie Bourdon committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
mark_prop(test,cross="N2",pval=0.05)
mark_prop(test,cross="N2",pval=0.05)
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)
}
mark_prop(test,cross="N2",pval=0.05)
test <- tibble(marker=c("mark1","mark2","mark3","mark4"),allele_1=c("A","G","T","C"),allele_2=c("T","T","A","G"),n_HM1=c(46,43,1,26),n_HM2=c(0,0,0,0),n_HT=c(50,55,48,74),n_NA=c(4,2,1,0))
mark_prop(test,cross="N2",pval=0.05)
Marie Bourdon's avatar
Marie Bourdon committed
478
479
480
481
482
483
484
485
486
487
488
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(dplyr)
library(stuart)
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
data(genos)
summary(genos)
data(phenos)
summary(phenos)
Marie Bourdon's avatar
Marie Bourdon committed
489
490
491
492
493
strains <- geno_strains(ref=annot_mini,geno=genos,par1=c("StrainsA_1","StrainsA_2"),
par2=c("StrainsB_1","StrainsB_2"),name1="parent1",name2="parent2")
head(strains) %>% print.data.frame()
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2",
"StrainsB_1","StrainsB_2"))
Marie Bourdon's avatar
Marie Bourdon committed
494
495
496
data(stuart_tab)
summary(stuart_tab)
tab2 <- mark_match(stuart_tab,ref=strains)
Marie Bourdon's avatar
Marie Bourdon committed
497
tab2 %>% filter(exclude_match==1) %>% print.data.frame()
Marie Bourdon's avatar
Marie Bourdon committed
498
tab2 <- mark_poly(tab2)
Marie Bourdon's avatar
Marie Bourdon committed
499
head(tab2) %>% print.data.frame()
Marie Bourdon's avatar
Marie Bourdon committed
500
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1)
Marie Bourdon's avatar
Marie Bourdon committed
501
head(tab2) %>% print.data.frame()
Marie Bourdon's avatar
Marie Bourdon committed
502
503
tab2 <- mark_allele(tab=tab2,ref=strains,par1="parent1",par2="parent2")
tab2 %>% arrange(desc(exclude_allele)) %>% head() %>% print.data.frame()
Marie Bourdon's avatar
Marie Bourdon committed
504
505
506
507
508
509
strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354",
"gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>%
select(marker,parent1,parent2) %>% print.data.frame()
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,
ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
rqtl_file[1:10,1:7] %>% print.data.frame()
Marie Bourdon's avatar
Marie Bourdon committed
510
511
512
devtools::build(path=".",vignettes=FALSE)
devtools::build_vignettes()
devtools::build_manual(path=".")