article_figures.Rmd 50.6 KB
Newer Older
Marie Bourdon's avatar
Marie Bourdon committed
1
2
3
---
title: "article_figures.Rmd"
author: "Marie Bourdon"
mariefbourdon's avatar
mariefbourdon committed
4
date: "June 2022"
Marie Bourdon's avatar
Marie Bourdon committed
5
6
7
output: html_document
---

Marie Bourdon's avatar
Marie Bourdon committed
8
# Goal and raw data
Marie Bourdon's avatar
Marie Bourdon committed
9
10
11

The goal of this script is to produce figure for the stuart package manuscript.

Marie Bourdon's avatar
Marie Bourdon committed
12
This scripts uses data from the package, and other files found in the "files" directory.
Marie Bourdon's avatar
Marie Bourdon committed
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

This script uses the qtl_plot() ggplot based function to plot QTL mapping results. This function is in the script "QTL_plot.R".

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(magrittr)
library(qtl)
library(cowplot)
library(grid)
library(gridExtra)
library(gridGraphics)

library(stuart)

source("files/QTL_plot.R")
```

Marie Bourdon's avatar
Marie Bourdon committed
31
# Data load and use of stuart functions
Marie Bourdon's avatar
Marie Bourdon committed
32

Marie Bourdon's avatar
Marie Bourdon committed
33
## genos dataframe
Marie Bourdon's avatar
Marie Bourdon committed
34
35
36
37
38
39

Data frame from stuart with genotypes of 176 F2 individuals and parental strains.
```{r genos}
data(genos)
summary(genos)
```
Marie Bourdon's avatar
Marie Bourdon committed
40
## phenos dataframe
Marie Bourdon's avatar
Marie Bourdon committed
41
42
43
44
45
46

Data frame from stuart with phenotypes of 176 F2 individuals for a quantitative trait.
```{r phenos}
data(phenos)
summary(phenos)
```
mariefbourdon's avatar
mariefbourdon committed
47

Marie Bourdon's avatar
Marie Bourdon committed
48
## annotation file
Marie Bourdon's avatar
Marie Bourdon committed
49

50
Annotation file from K.Broman: https://kbroman.org/MUGAarrays/mini_annotations.html
Marie Bourdon's avatar
Marie Bourdon committed
51
52
```{r annot}
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
Marie Bourdon's avatar
Marie Bourdon committed
53
summary(annot_mini)
Marie Bourdon's avatar
Marie Bourdon committed
54
```
Marie Bourdon's avatar
Marie Bourdon committed
55
## parental strains genotypes
Marie Bourdon's avatar
Marie Bourdon committed
56
57
58
59
60
61

2 genotypes tables will be used:

* One table with our genotypes of the strains used in the cross: "strains"

* One table with reference genotypes: "strns_ref"
Marie Bourdon's avatar
Marie Bourdon committed
62
63
64

```{r strains}
#our genotypes
Marie Bourdon's avatar
Marie Bourdon committed
65
66
67
68
69
70
71
72
73
74
75
76

#create tibble with individivual names
parental_strains <- tibble::tibble(parent1 = c("StrainsA_1","StrainsA_2"),
                                   parent2 = c("StrainsB_1","StrainsB_2"))


#create data frame with geno_strains
strains <- geno_strains(annot=annot_mini,geno=genos,
                        strn=parental_strains,cols=c("chr","cM_cox"))
rm(parental_strains)

#remove genotypes of parental strains from genos data frame
Marie Bourdon's avatar
Marie Bourdon committed
77
78
genos %<>% filter(!Sample.ID %in% c("StrainsA_1","StrainsA_2","StrainsB_1","StrainsB_2"))

Marie Bourdon's avatar
Marie Bourdon committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#summary
summary(strains)
```

```{r strns_ref}
#reference genotypes
#load parental strains genotype data from Neogen
strns_ref <- read.csv("files/ref_genotypes.csv")

#join with annotation file from miniMUGA
strns_ref <- strns_ref %>% select(name,parent1,parent2) %>% right_join(annot_mini,.,by=c("marker"="name")) %>% select(c(marker,chr,cM_cox,parent1,parent2))

#summary
summary(strns_ref)
```

Marie Bourdon's avatar
Marie Bourdon committed
95
## Use of stuart's functions
Marie Bourdon's avatar
Marie Bourdon committed
96
97
98
99
100
101
102
```{r}
data(stuart_tab)
summary(stuart_tab)

tab2 <- mark_match(stuart_tab,ref=strns_ref)
tab2 <- mark_poly(tab2)
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1))
Marie Bourdon's avatar
Marie Bourdon committed
103
tab2 <- mark_allele(tab2,ref=strns_ref,cross="F2",par1="parent1",par2="parent2")
Marie Bourdon's avatar
Marie Bourdon committed
104

Marie Bourdon's avatar
Marie Bourdon committed
105
data(stuart_newmap)
mariefbourdon's avatar
mariefbourdon committed
106
tab2 <- mark_estmap(tab=tab2,map=stuart_newmap,annot=annot_mini)
Marie Bourdon's avatar
Marie Bourdon committed
107
108
```

Marie Bourdon's avatar
Marie Bourdon committed
109
110
111
112
113
114
115
# Graphs

## Grid mark_prop

### Graph missing genotypes

```{r graph_NA}
mariefbourdon's avatar
mariefbourdon committed
116
na_plot <- tab2 %>% mutate(prop_NA=n_NA/176) %>% ggplot(aes(x=prop_NA)) +
Marie Bourdon's avatar
Marie Bourdon committed
117
118
  geom_histogram() +
  scale_y_log10() +
mariefbourdon's avatar
mariefbourdon committed
119
  theme_classic() +
Marie Bourdon's avatar
Marie Bourdon committed
120
121
122
123
124
125
126
127
128
  labs(title="Proportion of missing genotyped",
       x="Proportion of NA",y="Number of markers") +
  theme(
    aspect.ratio=0.8,
    plot.title = element_text(hjust = 0.4,face="bold",size=14))

na_plot
```

mariefbourdon's avatar
mariefbourdon committed
129
130
Proportions of markers with more than 75% of missing genotypes:
```{r prop_missing}
mariefbourdon's avatar
mariefbourdon committed
131
tab2 %>% mutate(prop_NA=n_NA/176) %>% filter(prop_NA > 0.50) %>% summarise(p=n()/count(tab2)%>%pull())
mariefbourdon's avatar
mariefbourdon committed
132
133
134
```


Marie Bourdon's avatar
Marie Bourdon committed
135
136
137
### Graph proportion of genotypes

```{r graph_prop}
mariefbourdon's avatar
mariefbourdon committed
138
prop_plot <- tab2 %>% filter(n_NA<88) %>% filter(!chr %in% c("M","X","Y")) %>%
mariefbourdon's avatar
mariefbourdon committed
139
  ggplot(aes(x=n_HM1/(n_HM1+n_HM2+n_HT),y=n_HM2/(n_HM1+n_HM2+n_HT),color=as.factor(exclude_prop))) +
Marie Bourdon's avatar
Marie Bourdon committed
140
  geom_point() +
mariefbourdon's avatar
mariefbourdon committed
141
  scale_color_manual(values=c("#66bd63","#b2182b"),labels = c("Retained", "Excluded")) +
mariefbourdon's avatar
mariefbourdon committed
142
143
144
  geom_hline(yintercept = 0.1,linetype="dashed",size=.3) +
  geom_vline(xintercept = 0.1,linetype="dashed",size=.3) +
  geom_abline(intercept = 0.9, slope=-1,linetype="dashed",size=.3) +
Marie Bourdon's avatar
Marie Bourdon committed
145
146
147
148
149
150
  #geom_text(aes(label=ifelse(exclude_prop=="1",SNP.Name,'')),hjust=0, vjust=0,size=2) +
  labs(title="Exclusion of markers with mark_prop()",
       x="Proportion of homozygous individuals AA",
       y="Proportion of homozygous individuals BB",
       color="Exclusion") +
  theme_classic() +
mariefbourdon's avatar
mariefbourdon committed
151
152
153
  theme(aspect.ratio=0.8,
        legend.position=c(0.8,0.8),
        legend.title = element_blank()) +
mariefbourdon's avatar
mariefbourdon committed
154
  theme(plot.title = element_text(hjust = 0.4,face="bold",size=14)) +
mariefbourdon's avatar
mariefbourdon committed
155
  geom_segment(x = 0.2, y = 0.78,
mariefbourdon's avatar
mariefbourdon committed
156
157
               xend = 0.25, yend = 0.72,
               color = "#1d91c0",
mariefbourdon's avatar
mariefbourdon committed
158
159
               arrow = arrow(type="closed",length=unit(4,"points"))) +
  geom_segment(x = 0.39, y = 0.58,
mariefbourdon's avatar
mariefbourdon committed
160
161
               xend = 0.44, yend = 0.52,
               color = "#1d91c0",
mariefbourdon's avatar
mariefbourdon committed
162
               arrow = arrow(type="closed",length=unit(4,"points")))
Marie Bourdon's avatar
Marie Bourdon committed
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
prop_plot
```
### Grid

Figure with graphs proportion of NA and exclusion depending on genotype proportions


```{r grid_na_prop}
grid <- plot_grid(na_plot,prop_plot,
                  ncol=1,
                  labels = c('A', 'B'), 
                  label_size = 20,
                  byrow = TRUE
                  )

grid

mariefbourdon's avatar
mariefbourdon committed
180
ggsave("figures/fig2.pdf",grid,width=5,height=8)
mariefbourdon's avatar
mariefbourdon committed
181
ggsave("figures/fig2.png",grid,width=5,height=8)
Marie Bourdon's avatar
Marie Bourdon committed
182
183
184
185
186
187
188
189
190
191

rm(na_plot,prop_plot)
```

## Graph before after

### Before: creation of Rqtl csv file 

```{r cross_before}
# filter at minima: remove non polymorphic and NA>0.5 
mariefbourdon's avatar
mariefbourdon committed
192
193
tab_before <- mark_na(stuart_tab)
tab_before <- mark_poly(tab_before)
Marie Bourdon's avatar
Marie Bourdon committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214

# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab_before,ref=strns_ref,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster/cross_before.csv")
```

### Before: newmap and permutation

These objects were produced on our cluster with the following script: /files/cluster/before_after.R

### Before: plot estimated map

```{r before_map}
# import cross
cross_before <- read.cross(format="csv",file="files/cluster/cross_before.csv",
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
load("files/cluster/newmap_before.rda")

# plot
plotMap(cross_before,newmap_before,shift=TRUE)
mariefbourdon's avatar
mariefbourdon committed
215
plotmap_before <- ~plotMap(cross_before,newmap_before,shift=TRUE,main="Before stuart", ylab='')
Marie Bourdon's avatar
Marie Bourdon committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
```

### Before: plot genome scan

```{r before_scan}
# load rda with perm
load("files/cluster/before_1000p.rda")
# load("files/before_1000p.rda")


# calculate thresholds
threshold_before <- summary(before_1000p,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63

# scanone
cross_before <- calc.genoprob(cross_before, step=2.0, off.end=0.0, 
                                 error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")


pheno_before <- scanone(cross=cross_before, chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), pheno.col="Pheno", model="normal", method="em")
summary(pheno_before)

# Plot
pheno_before_plot <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
mariefbourdon's avatar
mariefbourdon committed
243
244
         size=0.6,
         strip.pos="bottom") +
Marie Bourdon's avatar
Marie Bourdon committed
245
246
247
    theme(legend.position = "none") +
    ggtitle("")
pheno_before_plot
mariefbourdon's avatar
mariefbourdon committed
248
249
```

mariefbourdon's avatar
mariefbourdon committed
250

mariefbourdon's avatar
mariefbourdon committed
251

mariefbourdon's avatar
data3-4    
mariefbourdon committed
252
```{r before_plot}
mariefbourdon's avatar
mariefbourdon committed
253
254
255
pheno_before_zoom <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
mariefbourdon's avatar
220615    
mariefbourdon committed
256
257
258
         chrs = "13",
         size=0.6,
         rug = TRUE) +
mariefbourdon's avatar
mariefbourdon committed
259
260
261
262
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
mariefbourdon's avatar
mariefbourdon committed
263
264
    ggtitle("")
pheno_before_zoom
mariefbourdon's avatar
mariefbourdon committed
265
266
267
```

### Before figure
mariefbourdon's avatar
mariefbourdon committed
268

mariefbourdon's avatar
mariefbourdon committed
269
270
271
272
```{r grid_before}
grid_before <- plot_grid(as_grob(plotmap_before),pheno_before_plot,pheno_before_zoom,ncol=1,labels=c("A","B","C"),label_size=20)

ggsave("figures/fig1.pdf",grid_before,width=7,height = 15)
mariefbourdon's avatar
mariefbourdon committed
273
ggsave("figures/fig1.png",grid_before,width=7,height = 15)
Marie Bourdon's avatar
Marie Bourdon committed
274
275
276
277
278
279
280
281
282
```


### After: creation of Rqtl csv file 

```{r cross_after}
# filter with stuart functions: use the good data for parental strains (strains df)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 <- mark_poly(tab2)
mariefbourdon's avatar
mariefbourdon committed
283
tab2 <- mark_na(tab2)
Marie Bourdon's avatar
Marie Bourdon committed
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1))
tab2 <- mark_allele(tab2,ref=strains,cross="F2",par1="parent1",par2="parent2")

# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster/cross_after.csv")
```

### After: newmap 

These objects were produced on our cluster with the following script: /files/cluster/before_after.R

### After: plot estimated map

```{r after_map}
# import cross
cross_after <- read.cross(format="csv",file="files/cluster/cross_after.csv",
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
load("files/cluster/newmap_after.rda")

# plot
plotMap(cross_after,newmap_after,shift=TRUE)
plotmap_after <- ~plotMap(cross_after,newmap_after,shift=TRUE,main="After stuart")
```
### Remove last problematic markers with mark_estmap

```{r after_estmap}
#filter with mark_estmap
313
tab2 <- mark_estmap(tab2,newmap_after,annot_mini)
Marie Bourdon's avatar
Marie Bourdon committed
314
315

# create rqtl csv file
mariefbourdon's avatar
mariefbourdon committed
316
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster2/cross_after2.csv")
Marie Bourdon's avatar
Marie Bourdon committed
317
318
319
320
321
322
```

### After: plot estimated map 2

```{r after_map2}
# import cross
mariefbourdon's avatar
mariefbourdon committed
323
cross_after2 <- read.cross(format="csv",file="files/cluster2/cross_after2.csv",
Marie Bourdon's avatar
Marie Bourdon committed
324
325
326
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
mariefbourdon's avatar
mariefbourdon committed
327
load("files/cluster2/newmap_after2.rda")
Marie Bourdon's avatar
Marie Bourdon committed
328
329

# plot
mariefbourdon's avatar
mariefbourdon committed
330
331
plotMap(cross_after2,newmap_after2,shift=TRUE)
plotmap_after <- ~plotMap(cross_after2,newmap_after2,shift=TRUE,main="After stuart")
Marie Bourdon's avatar
Marie Bourdon committed
332
333
334
335
336
337
```

### After: plot genome scan

```{r after_scan}
# load rda with perm
mariefbourdon's avatar
mariefbourdon committed
338
load("files/cluster2/after_1000p2.rda")
Marie Bourdon's avatar
Marie Bourdon committed
339
340

# calculate thresholds
mariefbourdon's avatar
mariefbourdon committed
341
threshold_after <- summary(after_1000p2,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63
Marie Bourdon's avatar
Marie Bourdon committed
342
343

# scanone
mariefbourdon's avatar
mariefbourdon committed
344
cross_after <- calc.genoprob(cross_after2, step=2.0, off.end=0.0, 
Marie Bourdon's avatar
Marie Bourdon committed
345
346
347
                                 error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")


mariefbourdon's avatar
mariefbourdon committed
348
pheno_after <- scanone(cross=cross_after2, chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), pheno.col="Pheno", model="normal", method="em")
Marie Bourdon's avatar
Marie Bourdon committed
349
350
351
352
353
354
355
356
summary(pheno_after)

# Plot
pheno_after_plot <- qtl_plot(pheno_after,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_after[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
mariefbourdon's avatar
mariefbourdon committed
357
358
         size=0.6,
         strip.pos="bottom") +
Marie Bourdon's avatar
Marie Bourdon committed
359
360
361
    theme(legend.position = "none") +
    ggtitle("")
pheno_after_plot
mariefbourdon's avatar
mariefbourdon committed
362
363
364
365
366
367
368
369
370
371
372

qtl_plot(pheno_after,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_after[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="10") +
    theme(legend.position = "none") +
    ggtitle("")
Marie Bourdon's avatar
Marie Bourdon committed
373
```
mariefbourdon's avatar
mariefbourdon committed
374

mariefbourdon's avatar
mariefbourdon committed
375
### Grid after
Marie Bourdon's avatar
Marie Bourdon committed
376
377

```{r grid fig1, fig.height = 7, fig.width = 13, fig.align = "center"}
mariefbourdon's avatar
mariefbourdon committed
378
grid_after <- plot_grid(as_grob(plotmap_after),pheno_after_plot,ncol=1,labels=c("A","B"),label_size=20)
Marie Bourdon's avatar
Marie Bourdon committed
379

mariefbourdon's avatar
mariefbourdon committed
380
ggsave("figures/fig3.pdf",grid_after,width=7,height = 10)
mariefbourdon's avatar
mariefbourdon committed
381
ggsave("figures/fig3.png",grid_after,width=7,height = 10)
Marie Bourdon's avatar
Marie Bourdon committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
```

# Differences between real and reference genotypes for parental strains

```{r dif_number}
dif <- full_join(strains,strns_ref,by=c("marker","chr","cM_cox")) %>% 
  mutate(dif=case_when((!parent1.x%in%c("N","H") & 
                          !parent1.y%in%c("N","H") & 
                          parent1.x!=parent1.y) ~ 1,
                       (!parent2.x%in%c("N","H") &
                          !parent2.y%in%c("N","H") & 
                          parent2.x!=parent2.y) ~ 1, 
                       T~0))

dif %>% filter(dif==1) %>% count()
```

mariefbourdon's avatar
mariefbourdon committed
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
## Table alleles different between parental strains and F2s

```{r allele}
#investigation of the role of mark_allele function
#prove that some marker with non corresponding alleles between parents and F2

#keep only markers that are exlcuded with mark_allele
allele <- tab2%>% 
  filter(exclude_allele==1&exclude_poly==0&exclude_prop==0&exclude_na==0&exclude_estmap==0)
strains_allele <- strains %>% filter(marker %in% allele$marker)

#number of markers not excluded by other functoins
tab2 %>% 
  filter(exclude_poly==0&exclude_prop==0&exclude_na==0&exclude_estmap==0) %>% nrow()

#join with strains genotypes to have parental strains
allele <- left_join(allele,strains_allele,by=c("marker"="marker"))


#most of markers excluded with mark_allele that were not excluded with other functions have N/H as genotype for parents
#keep only those with non missing/heterozygous genotypes
allele %<>% filter(parent1 != "N" & parent2 != "N")
allele %<>% select(marker,parent1,parent2,allele_1,allele_2)


print(allele)

print(xtable::xtable(allele, type = "latex"), file = "tables/tab_alleles.tex",include.rownames=FALSE)
rm(allele,strains_allele)
```

Marie Bourdon's avatar
Marie Bourdon committed
430
431
432
```{r dif_table}
table_dif <- dif %>% filter(dif==1) %>% select(marker,parent1_ref=parent1.y,parent1_geno=parent1.x,parent2_ref=parent2.y,parent2_geno=parent2.x) %>% head()
knitr::kable(table_dif)
mariefbourdon's avatar
mariefbourdon committed
433
434
```

mariefbourdon's avatar
mariefbourdon committed
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
## Number of markers kept after each function

```{r barplot}
none <- tab2 %>% nrow()
match <- tab2 %>% filter(exclude_match==0) %>% nrow()
allele <- tab2 %>% filter(exclude_match==0&exclude_allele==0) %>% nrow()
naf <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0) %>% nrow()
poly <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0) %>% nrow()
prop <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0&exclude_prop==0) %>% nrow()
estmap <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0&exclude_prop==0&exclude_estmap==0) %>% nrow()

functions_df <- tibble(fct=c("none","match","allele","na","poly","prop","estmap"),
                       markers=c(none,match,allele,naf,poly,prop,estmap))

functions_plot <- functions_df %>% ggplot(aes(x=markers,y=fct)) +
  geom_bar(stat="identity",width=0.6) +
  geom_text(aes(label=markers), hjust=1.3, color="white", size=3.5) +
  scale_y_discrete(limits=c("estmap","prop","poly", "na", "allele","match","none")) +
  theme(aspect.ratio=0.7) +
  labs(title="Number of markers kept after each step",
       x="Number of markers",
       y="Function used") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.4,face="bold",size=14))

functions_plot
rm(none,allele,match,poly,prop,barplot_df)
```

mariefbourdon's avatar
mariefbourdon committed
464
465
## Pheno data format

mariefbourdon's avatar
mariefbourdon committed
466
467
```{r pheno}
format_pheno <- phenos[1:6,]
mariefbourdon's avatar
mariefbourdon committed
468
write_csv(format_pheno,file="sup/tableS1.csv")
mariefbourdon's avatar
mariefbourdon committed
469
470
471
print(xtable::xtable(format_pheno, type = "latex"), file = "tables/tab_alleles.tex",include.rownames=FALSE)
```

mariefbourdon's avatar
mariefbourdon committed
472
473
474
475
476
477
478
479
480
481
## est rf

```{r}
source("files/find_linked_markers.R")
estrf_matrix_after <- pull.rf(cross_after, what=c("lod"))

find_linked_markers(estrf_matrix_after,mark="S6J011498219",annot=annot_mini)
```


mariefbourdon's avatar
data3-4    
mariefbourdon committed
482
483
484
485
## Narrow peaks

Investigation of high lod score peaks 

mariefbourdon's avatar
mariefbourdon committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
```{r before_ann}
chrs <- c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X")
ann_dat_text<-data.frame(
    chr=factor(chrs,
               levels=chrs),
    lod=c(NA,7.5,NA,NA,15,NA,10.5,10.5,NA,5.5,NA,NA,18.3,NA,NA,NA,NA,NA,NA,NA),
    label=c(NA,"1",NA,NA,"2",NA,"3","4",NA,"5",NA,NA,"6",NA,NA,NA,NA,NA,NA,NA),
    x=c(NA,17,NA,NA,27,NA,70,10,NA,30,NA,NA,33,NA,NA,NA,NA,NA,NA,NA)
)

pheno_before_annot <- pheno_before_plot +  geom_text(
    # the new dataframe for annotating text
    data = ann_dat_text,
    mapping = aes(x = x,
                  y = lod,
                  label = label,
                  color="blue")
)
pheno_before_annot

rm(ann_dat_text)
rm(chrs)
```


mariefbourdon's avatar
mariefbourdon committed
511
### Peak 1
mariefbourdon's avatar
data3-4    
mariefbourdon committed
512

mariefbourdon's avatar
mariefbourdon committed
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
```{r peak1_zoom}
peak1 <- qtl_plot(pheno_before,
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="2",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("Peak 1")
peak1
```


1 peak on chromosome 2 1 pseudomarker : c2.loc10, postionned between gUNC2731905 and gUNCHS004244.
mariefbourdon's avatar
data3-4    
mariefbourdon committed
532
533
534
535

Here are the infos on genotype counts for these markers:

```{r summary_geno_chr2}
mariefbourdon's avatar
mariefbourdon committed
536
537
peak1_tab <- tab_before %>% filter(marker %in% c("gUNC2731905","gUNCHS004244")) %>% select(marker:n_NA)
peak1_tab
mariefbourdon's avatar
data3-4    
mariefbourdon committed
538
539
540
541
542
543
544
```

For gUNC2731905, all individuals except 1 are homozygous so this marker should be removed. The proportions for gUNCHS004244 seem correct.


Graph:

mariefbourdon's avatar
mariefbourdon committed
545
```{r geno_plot_peak1}
mariefbourdon's avatar
data3-4    
mariefbourdon committed
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["2"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["2"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(SNJ020344053:mUNC4608754,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

geno_plot2 <- pgmap %>% filter(pos > 10 & pos < 20) %>%
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1,size=3,
           x = map %>% filter(pos > 10 & pos < 20) %>% pull(pos),
           label = map %>% filter(pos > 10 & pos < 20) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
geno_plot2
mariefbourdon's avatar
mariefbourdon committed
569
570

rm(pgmap,phenotypes,map,genotypes,phenogeno)
mariefbourdon's avatar
data3-4    
mariefbourdon committed
571
572
```

mariefbourdon's avatar
mariefbourdon committed
573
This region contains many markers with non Mendelian proportions (many have an excess of homozygous). This leads to the creation of a pseudomarker between  gUNC2731905 and gUNCHS004244 with incorrect genotypic probabilities which leads to a spurious peak.
mariefbourdon's avatar
data3-4    
mariefbourdon committed
574
575
576


```{r}
mariefbourdon's avatar
mariefbourdon committed
577
newmap_before[["2"]][8:20]
mariefbourdon's avatar
data3-4    
mariefbourdon committed
578
579
```

mariefbourdon's avatar
mariefbourdon committed
580
581
582
Recombination fractions between adjacent markers in this regions are indeed high.

### Peak 2
mariefbourdon's avatar
data3-4    
mariefbourdon committed
583

mariefbourdon's avatar
mariefbourdon committed
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
```{r peak7_zoom}
peak2 <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="5",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("")
peak2
```

1 peak on chromosome 2 on 1 marker : mUNC050096588
mariefbourdon's avatar
data3-4    
mariefbourdon committed
603
604
605
606

Here are the infos on genotype counts for these markers:

```{r summary_geno_chr5}
mariefbourdon's avatar
mariefbourdon committed
607
608
peak2_tab <- tab_before %>% filter(marker %in% c("mUNC050096588")) %>% select(marker:n_NA)
peak2_tab
mariefbourdon's avatar
data3-4    
mariefbourdon committed
609
610
```

mariefbourdon's avatar
mariefbourdon committed
611
All individuals heterozygous so this marker should be removed.
mariefbourdon's avatar
data3-4    
mariefbourdon committed
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644

```{r}
newmap_before[["5"]][18:23]
```

This leads to enormous calculated distance with adjacent markers: 1186.0871-184.5811 = 2187.5930-1186.0871=1001.506 cM with the previous and the following marker.

Graph:

```{r geno_plot_chr5}
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["5"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["5"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(gUNCHS013469:SAC056009450,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

test_plot <- pgmap %>% filter(pos > 20 & pos < 30) %>% 
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1,size=3,
           x = map %>% filter(pos > 20 & pos < 30) %>% pull(pos),
           label = map %>% filter(pos > 20 & pos < 30) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
test_plot 
mariefbourdon's avatar
mariefbourdon committed
645
rm(pgmap,phenotypes,map,genotypes,phenogeno)
mariefbourdon's avatar
data3-4    
mariefbourdon committed
646
647
```

mariefbourdon's avatar
mariefbourdon committed
648
### Peak 3
mariefbourdon's avatar
data3-4    
mariefbourdon committed
649

mariefbourdon's avatar
mariefbourdon committed
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
```{r peak3_zoom}
peak3 <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="7",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("")
peak3
```

2 peaks on chromosome 7, one on 1 pseudomarker : c7.loc74, between UNC13823755 and S3J075374098; and one on 1 marker and 1 pseudomarker : c7.loc82 and gUNC13998623, c7.loc82 being between gUNC13979374 and gUNC13998623.
mariefbourdon's avatar
data3-4    
mariefbourdon committed
669
670
671

Here are the infos on genotype counts for these markers:

mariefbourdon's avatar
mariefbourdon committed
672
673
674
```{r summary_geno_peak3}
peak3_tab <- tab_before %>% filter(marker %in% c("UNC13823755", "S3J075374098","gUNC13979374","gUNC13998623")) %>% select(marker:n_NA)
peak3_tab 
mariefbourdon's avatar
data3-4    
mariefbourdon committed
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
```

There are no homozygous individuals for one allele for S3J075374098, gUNC13979374 and gUNC13998623 so these markers should be removed. Proportions seem correct for UNC13823755.

```{r}
newmap_before[["7"]][206:212]
```

Indeed, there is not an increased in the calculated distance with the markers before UNC13823755, but there is 29066.20-28864.88=201.32 cM between UNC13823755 and S3J075374098 and 29116.71-29076.90=39.81 cM between gUNC13979374 and the previous marker.

Graph:

```{r geno_plot_chr7}
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["7"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["7"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(SBJ070191318:gUNCHS021959,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

test_plot <- pgmap %>% filter(pos > 65 & pos < 90) %>% 
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1.5,size=3,
           x = map %>% filter(pos > 65 & pos < 90) %>% pull(pos),
           label = map %>% filter(pos > 65 & pos < 90) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
test_plot 
```

mariefbourdon's avatar
mariefbourdon committed
713
### Peak 4
mariefbourdon's avatar
data3-4    
mariefbourdon committed
714

mariefbourdon's avatar
mariefbourdon committed
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
```{r peak4_zoom}
peak4 <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="8",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("")
peak4
```

1 peak on chromosome 8 one on 1 marker : mbackupJAX00158395
mariefbourdon's avatar
data3-4    
mariefbourdon committed
734
735
736

Here are the infos on genotype counts for this marker:

mariefbourdon's avatar
mariefbourdon committed
737
738
739
```{r summary_geno_peak4}
peak4_tab <- tab_before %>% filter(marker %in% c("mbackupJAX00158395")) %>% select(marker:n_NA)
peak4_tab 
mariefbourdon's avatar
data3-4    
mariefbourdon committed
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
```

All individuals are homozygous except one so this marker should be removed.

```{r}
newmap_before[["8"]][1:5]
```

Indeed there are 189.3831-2.1670=187.2161 cM between mbackupJAX00158395 and the following marker.

Graph:

```{r geno_plot_chr8}
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["8"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["8"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(mbackupJAX00158395:S6J085066393,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

test_plot <- pgmap %>% filter(pos > 0 & pos < 10) %>% 
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1.5,size=3,
           x = map %>% filter(pos > 0 & pos < 10) %>% pull(pos),
           label = map %>% filter(pos > 0 & pos < 10) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
test_plot 
```


mariefbourdon's avatar
mariefbourdon committed
779
### Peak 5
mariefbourdon's avatar
data3-4    
mariefbourdon committed
780

mariefbourdon's avatar
mariefbourdon committed
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
```{r peak5_zoom}
peak5 <- qtl_plot(pheno_before,
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="10",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("Peak 5")
peak5
```

1 peak on chromosome 10 on 1 marker : S6J102311553
mariefbourdon's avatar
data3-4    
mariefbourdon committed
799
800
801

Here are the infos on genotype counts for this marker:

mariefbourdon's avatar
mariefbourdon committed
802
```{r summary_geno_peak5}
mariefbourdon's avatar
data3-4    
mariefbourdon committed
803
804
805
806
807
tab_before %>% filter(marker %in% c("S6J102311553")) %>% select(marker:n_NA)
```

The genotypic proportions for this marker seem correct.

mariefbourdon's avatar
mariefbourdon committed
808
```{r parents_geno_peak5}
mariefbourdon's avatar
data3-4    
mariefbourdon committed
809
810
811
812
813
814
815
816
strns_ref %>% filter(marker %in% c("S6J102311553"))
strains %>% filter(marker %in% c("S6J102311553"))
```

The alleles of parental strains seemed correct too and are the same in the reference alleles data and in our genotypes.

Graph:

mariefbourdon's avatar
mariefbourdon committed
817
```{r geno_plot_peak5}
mariefbourdon's avatar
data3-4    
mariefbourdon committed
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["10"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["10"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(gICR258:gUNC18984159,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

test_plot <- pgmap %>% filter(pos > 22 & pos < 36.5) %>% 
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1.5,size=3,
           x = map %>% filter(pos > 22 & pos < 36.5) %>% pull(pos),
           label = map %>% filter(pos > 22 & pos < 36.5) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
test_plot 
```

S6J102311553 seem to be surrounded by markers with an increased proportion of homozygous individuals.

mariefbourdon's avatar
mariefbourdon committed
845
846
847
848
849
```{r summary_geno_peak5}
peak5_tab <-tab_before %>% filter(marker %in% c("SSR102275149","gUNC18022023","S6J102311553","gUNC18048671","B10100061261")) %>% select(marker:n_NA)
peak5_tab
```

mariefbourdon's avatar
data3-4    
mariefbourdon committed
850
851
852
853
854
855
856

```{r}
newmap_before[["10"]][80:110]
```

Indeed, there is not an increased distance between S6J102311553 and the adjacent markers; but S6J102311553 is in the middle of a region with enormous calculated with its adjacent markers, between S2T101968115 (showing 19084.57-18083.06=1001.51 cM with the previous marker) and S2T102642258 (showing 20373.93-19372.42=1001.51 cM with the following marker)

mariefbourdon's avatar
mariefbourdon committed
857
### Peak 6
mariefbourdon's avatar
data3-4    
mariefbourdon committed
858

mariefbourdon's avatar
mariefbourdon committed
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
```{r peak6_zoom}
peak6 <- qtl_plot(pheno_before,
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom",
         chr="13",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("Peak 6")
peak6
```

1 peak on chromosome 13 on 1 marker : SAC132487883
mariefbourdon's avatar
data3-4    
mariefbourdon committed
877
878
879
880

Here are the infos on genotype counts for this marker:

```{r summary_geno_chr13}
mariefbourdon's avatar
mariefbourdon committed
881
882
peak6_tab <- tab_before %>% filter(marker %in% c("SAC132487883")) %>% select(marker:n_NA)
peak6_tab
mariefbourdon's avatar
data3-4    
mariefbourdon committed
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
```

All individuals are heterozygous at this loci so this marker should be removed.

```{r}
newmap_before[["13"]][48:50]
```

This leads to enormous calculated distance with adjacent markers: 9521.842-8520.336 = 8520.336-7518.830 = 1001.506 cM with the previous and the following marker.


Graph:

```{r geno_plot_chr13}
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["13"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["13"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(gUNCHS034900:S3H134792711,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

test_plot <- pgmap %>% filter(pos > 25 & pos < 35) %>% 
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1.7,size=3,
           x = map %>% filter(pos > 25 & pos < 35) %>% pull(pos),
           label = map %>% filter(pos > 25 & pos < 35) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 55)))
test_plot 
```

## Fold change distance between adjacent markers

```{r fold_change}
#before
names_mark <- c(names(newmap_before[["1"]]),names(newmap_before[["2"]]),names(newmap_before[["3"]]),names(newmap_before[["4"]]),
  names(newmap_before[["5"]]),names(newmap_before[["6"]]),names(newmap_before[["7"]]),names(newmap_before[["8"]]),
  names(newmap_before[["9"]]),names(newmap_before[["10"]]),names(newmap_before[["11"]]),names(newmap_before[["12"]]),
  names(newmap_before[["13"]]),names(newmap_before[["14"]]),names(newmap_before[["15"]]),names(newmap_before[["16"]]),
  names(newmap_before[["17"]]),names(newmap_before[["18"]]),names(newmap_before[["19"]]),names(newmap_before[["X"]]))
pos_mark <- c(newmap_before[["1"]],newmap_before[["2"]],newmap_before[["3"]],newmap_before[["4"]],
  newmap_before[["5"]],newmap_before[["6"]],newmap_before[["7"]],newmap_before[["8"]],
  newmap_before[["9"]],newmap_before[["10"]],newmap_before[["11"]],newmap_before[["12"]],
  newmap_before[["13"]],newmap_before[["14"]],newmap_before[["15"]],newmap_before[["16"]],
  newmap_before[["17"]],newmap_before[["18"]],newmap_before[["19"]],newmap_before[["X"]])
tibble_newmap_before <- tibble(marker=names_mark,
                               cM_calc=pos_mark)

compar_pos_before <- full_join(tibble_newmap_before,annot_mini) %>% select(marker,chr,cM_calc,cM_cox)
know <- compar_pos_before$cM_cox
calc <- compar_pos_before$cM_calc
compar_pos_before <- tibble(marker=compar_pos_before$marker,
       chr=compar_pos_before$chr,
       cM_cox=compar_pos_before$cM_cox,
       cox_prev=c(NA,compar_pos_before$cM_cox[1:11124]),
       cox_fol=c(compar_pos_before$cM_cox[2:11125],NA),
       cM_calc=compar_pos_before$cM_calc,
       calc_prev=c(NA,compar_pos_before$cM_calc[1:11124]),
       calc_fol=c(compar_pos_before$cM_calc[2:11125],NA)) %>%
  mutate(dif_prev=calc_prev/cox_prev,
         dif_fol=calc_fol/cox_fol)

#after
names_mark <- c(names(newmap_after2[["1"]]),names(newmap_after2[["2"]]),names(newmap_after2[["3"]]),names(newmap_after2[["4"]]),
  names(newmap_after2[["5"]]),names(newmap_after2[["6"]]),names(newmap_after2[["7"]]),names(newmap_after2[["8"]]),
  names(newmap_after2[["9"]]),names(newmap_after2[["10"]]),names(newmap_after2[["11"]]),names(newmap_after2[["12"]]),
  names(newmap_after2[["13"]]),names(newmap_after2[["14"]]),names(newmap_after2[["15"]]),names(newmap_after2[["16"]]),
  names(newmap_after2[["17"]]),names(newmap_after2[["18"]]),names(newmap_after2[["19"]]),names(newmap_after2[["X"]]))
pos_mark <- c(newmap_after2[["1"]],newmap_after2[["2"]],newmap_after2[["3"]],newmap_after2[["4"]],
  newmap_after2[["5"]],newmap_after2[["6"]],newmap_after2[["7"]],newmap_after2[["8"]],
  newmap_after2[["9"]],newmap_after2[["10"]],newmap_after2[["11"]],newmap_after2[["12"]],
  newmap_after2[["13"]],newmap_after2[["14"]],newmap_after2[["15"]],newmap_after2[["16"]],
  newmap_after2[["17"]],newmap_after2[["18"]],newmap_after2[["19"]],newmap_after2[["X"]])
tibble_newmap_after <- tibble(marker=names_mark,
                               cM_calc=pos_mark)

compar_pos_after <- full_join(tibble_newmap_after,annot_mini) %>% select(marker,chr,cM_calc,cM_cox)
know <- compar_pos_after$cM_cox
calc <- compar_pos_after$cM_calc
compar_pos_after <- tibble(marker=compar_pos_after$marker,
       chr=compar_pos_after$chr,
       cM_cox=compar_pos_after$cM_cox,
       cox_prev=c(NA,compar_pos_after$cM_cox[1:11124]),
       cox_fol=c(compar_pos_after$cM_cox[2:11125],NA),
       cM_calc=compar_pos_after$cM_calc,
       calc_prev=c(NA,compar_pos_after$cM_calc[1:11124]),
       calc_fol=c(compar_pos_after$cM_calc[2:11125],NA)) %>%
  mutate(dif_prev=calc_prev/cox_prev,
         dif_fol=calc_fol/cox_fol)

mean(compar_pos_before$dif_prev,na.rm=TRUE)
mariefbourdon's avatar
mariefbourdon committed
982
sd(compar_pos_before$dif_prev,na.rm=TRUE)
mariefbourdon's avatar
data3-4    
mariefbourdon committed
983
mean(compar_pos_after$dif_prev,na.rm=TRUE)
mariefbourdon's avatar
mariefbourdon committed
984
sd(compar_pos_after$dif_prev,na.rm=TRUE)
mariefbourdon's avatar
data3-4    
mariefbourdon committed
985
986
```

mariefbourdon's avatar
mariefbourdon committed
987
988
989
990
991
992
993
994
995
996
997
998
999
1000


```{r}
# #pgm: non
# cross_test <- cross_before
# cross_test[["pheno"]][["pgm"]] <- cross_after2[["pheno"]][["pgm"]]
# identical(cross_test[["pheno"]][["pgm"]],cross_after2[["pheno"]][["pgm"]])
# 
# 
# pheno_test <- scanone(cross=cross_test, chr=c("5"), pheno.col="Pheno", model="normal", method="em")
# pheno_before5 <- scanone(cross=cross_before, chr=c("5"), pheno.col="Pheno", model="normal", method="em")
# 
# 
# identical(pheno_test,pheno_before5)
For faster browsing, not all history is shown. View entire blame