diff --git a/article/article_figures.Rmd b/article/article_figures.Rmd
index 047b893c1371be614febd7b6bf64d7d6c7f873ff..98c9b3e9325b098cce621b0e5c99e2a08ca5774a 100644
--- a/article/article_figures.Rmd
+++ b/article/article_figures.Rmd
@@ -22,10 +22,12 @@ library(cowplot)
library(grid)
library(gridExtra)
library(gridGraphics)
+library(scales)
library(stuart)
source("files/QTL_plot.R")
+source("files/find_linked_markers.R")
```
# Data load and use of stuart functions
@@ -306,6 +308,7 @@ load("files/cluster/newmap_after.rda")
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}
@@ -314,6 +317,17 @@ tab2 <- mark_estmap(tab2,newmap_after,annot_mini)
# 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/cluster2/cross_after2.csv")
+
+
+
+# are these markers misplaced ?
+tab2 %>% filter(marker %in% c("S6J010381992","SS6071602326","S6J205609960"))
+
+estrf_matrix_after <- pull.rf(cross_after, what=c("lod"))
+
+find_linked_markers(estrf_matrix_after,mark="S6J010381992",annot=annot_mini)
+find_linked_markers(estrf_matrix_after,mark="SS6071602326",annot=annot_mini)
+find_linked_markers(estrf_matrix_after,mark="S6J205609960",annot=annot_mini)
```
### After: plot estimated map 2
@@ -472,10 +486,9 @@ print(xtable::xtable(format_pheno, type = "latex"), file = "tables/tab_alleles.t
## est rf
```{r}
-source("files/find_linked_markers.R")
-estrf_matrix_after <- pull.rf(cross_after, what=c("lod"))
+estrf_matrix_after2 <- pull.rf(cross_after2, what=c("lod"))
-find_linked_markers(estrf_matrix_after,mark="S6J011498219",annot=annot_mini)
+find_linked_markers(estrf_matrix_after2,mark="S6J011498219",annot=annot_mini)
```
@@ -919,72 +932,6 @@ test_plot <- pgmap %>% filter(pos > 25 & pos < 35) %>%
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)
-sd(compar_pos_before$dif_prev,na.rm=TRUE)
-mean(compar_pos_after$dif_prev,na.rm=TRUE)
-sd(compar_pos_after$dif_prev,na.rm=TRUE)
-```
-
-
```{r}
# #pgm: non
@@ -1275,6 +1222,7 @@ write_csv(allele_rec,"sup/tableS2.csv")
rm(map_1,map_2,map_3,map_4,pos_1,pos_2,pos_3,pos_4,df)
```
+## Grid with all 4 data sets
```{r}
load("data2/data2_peaks.rda")
@@ -1311,10 +1259,10 @@ narrow_grid
ggsave(narrow_grid,file="sup/figureS1.pdf",width=10,height=17)
# Figure S1 : Analysis of the F2 cross data illustrating the identification of spurious narrow non-significant peaks in QTL mappig in 4 data sets.
-# A: output of the scanone function of rqtl in a F2 between CC001 and Ifnar KO C57BL/6J showing the identification of 6 narrow non-significant peaks.
-# B: output of the scanone function of rqtl in a F2 between Ifnar KO C57BL/6J and Ifnar KO 129S2/SvPas showing the identification of 1 narrow non-significant peak.
-# C: output of the scanone function of rqtl in a (CC001xCC071)xCC071 backcross showing no narrow non-significant peak.
-# D: output of the scanone function of rqtl in a F2 between Ifnar KO C57BL/6N and CC021 showing the identification of 5 narrow non-significant peaks.
+# A: output of the scanone function of rqtl in a F2 between CC001 and Ifnar KO C57BL/6J showing the identification of 6 narrow non-significant peaks.
+# B: output of the scanone function of rqtl in a F2 between Ifnar KO C57BL/6J and Ifnar KO 129S2/SvPas showing the identification of 1 narrow non-significant peak.
+# C: output of the scanone function of rqtl in a (CC001xCC071)xCC071 backcross showing no narrow non-significant peak.
+# D: output of the scanone function of rqtl in a F2 between Ifnar KO C57BL/6N and CC021 showing the identification of 5 narrow non-significant peaks.
# E: Zoom on Peak 1. Peak 1 is located on a pseudomarker next to a marker with non mendelian proportions (gUNC2731905). F: alleles and genotype proportions for the adjacent markers of Peak 1. Peaks 3, 7, 8 and 9 are also due to a pseudomarkers located next to a markers with non mendelian proportions
# G: Zoom on Peak 6. Peak 6 is located on a marker with non mendelian proportions (SAC132487883). H: alleles and genotype proportions for SAC132487883. Peaks 2, 4 and 12 are also due to a markers with non mendelian proportions.
# I: Zoom on Peak 11. Peak 11 is located on one marker with non mendelian proportions (SAC132487883) and one adjacent pseudomarker. J: alleles and genotype proportions for SAC132487883 and the other adjacent marker to the pseudomarker on the peak. Peak 10 is also due do one marker with non mendelian proportions and one adjacent pseudomarker.
@@ -1347,16 +1295,35 @@ for(i in names(newmap_before)){
follow <- c(follow,fol)
}
+annot <- annot_mini %>% filter(marker %in% mark)
+kn_pos <- annot$cM_cox
+kn_prev <- c(NA, annot[1:(nrow(annot) - 1), "cM_cox"])
+kn_previous <- c(kn_previous, kn_prev)
+kn_fol <- c(annot[2:nrow(annot), "cM_cox"], NA)
+kn_follow <- c(kn_follow, kn_fol)
#create tab with positions
tab_map1 <- tibble(marker = mark,
chr = chr,
place = place,
pos = pos,
- previous = pos-previous,
- follow = follow-pos)
+ previous = previous,
+ prev_dif = pos-previous,
+ follow = follow,
+ fol_dif = follow-pos,
+ kn_pos = kn_pos,
+ kn_previous = kn_previous,
+ kn_prev_dif = kn_pos - kn_previous,
+ kn_follow = kn_follow,
+ kn_fol_dif = kn_follow - kn_pos)
+
+
+tab_map1 <- tab_map1 %>% mutate(kn_prev_df = case_when(is.na(previous) == TRUE ~ NA_real_, T ~ kn_previous))
+tab_map1 <- tab_map1 %>% mutate(kn_fol_dof = case_when(is.na(follow) == TRUE ~ NA_real_, T ~ kn_follow))
+
+tab_map1 <- tab_map1 %>% mutate(rat_prev = prev_dif/kn_prev_dif)
-tab_map1 %>% ggplot(aes(x=follow)) +
+tab_map1 %>% filter(is.na(rat_prev)==FALSE & rat_prev != Inf) %>% ggplot(aes(x=rat_prev)) +
geom_density() +
scale_x_log10()
@@ -1382,25 +1349,71 @@ for(i in names(newmap_after2)){
follow <- c(follow,fol)
}
+annot <- annot_mini %>% filter(marker %in% mark)
+kn_pos <- annot$cM_cox
+kn_prev <- c(NA, annot[1:(nrow(annot) - 1), "cM_cox"])
+kn_previous <- c(kn_previous, kn_prev)
+kn_fol <- c(annot[2:nrow(annot), "cM_cox"], NA)
+kn_follow <- c(kn_follow, kn_fol)
+
#create tab with positions
tab_map2 <- tibble(marker = mark,
chr = chr,
place = place,
pos = pos,
- previous = pos-previous,
- follow = follow-pos)
+ previous = previous,
+ prev_dif = pos-previous,
+ follow = follow,
+ fol_dif = follow-pos,
+ kn_pos = kn_pos,
+ kn_previous = kn_previous,
+ kn_prev_dif = kn_pos - kn_previous,
+ kn_follow = kn_follow,
+ kn_fol_dif = kn_follow - kn_pos)
+
+tab_map2 <- tab_map2 %>% mutate(kn_prev_df = case_when(is.na(previous) == TRUE ~ NA_real_, T ~ kn_previous))
+tab_map2 <- tab_map2 %>% mutate(kn_fol_dof = case_when(is.na(follow) == TRUE ~ NA_real_, T ~ kn_follow))
+
+tab_map2 <- tab_map2 %>% mutate(rat_prev = prev_dif/kn_prev_dif)
+
+
+# save in new df
+
+rec_ratios <- bind_rows(tab_map1 %>% select(marker,chr,kn_pos,kn_previous,kn_prev_dif,pos,previous,prev_dif,rat_prev) %>% mutate(state="before"),
+ tab_map2 %>% select(marker,chr,kn_pos,kn_previous,kn_prev_dif,pos,previous,prev_dif,rat_prev) %>% mutate(state="after"))
+
+fancy_scientific <- function(l) {
+ # turn in to character string in scientific notation
+ l <- format(l, scientific = TRUE)
+ # quote the part before the exponent to keep all the digits
+ l <- gsub("^(.*)e", "'\\1'e", l)
+ # turn the 'e+' into plotmath format
+ l <- gsub("e", "%*%10^", l)
+ # return this as an expression
+ parse(text=l)
+}
-ggplot() +
- geom_density(data = tab_map1,aes(x=follow),color="red") +
- geom_density(data = tab_map2,aes(x=follow),color="blue") +
- scale_x_log10()
+rec_ratios %>% ggplot(aes(x=rat_prev,color=state)) +
+ geom_density() +
+ scale_color_manual(values=c("blue","red"),labels=c("After filtering","Before filtering")) +
+ scale_x_log10(labels=fancy_scientific) +
+ labs(x="Ratio between the calculated and the known distance with the previous marker",
+ y="Density",
+ color="",
+ title="Fold change distance between adjacent markers") +
+ ggpubr::theme_classic2()
ggsave("distrib_freq_rec1.png")
-ggplot() +
- geom_density(data = tab_map1,aes(x=follow),color="red") +
- geom_density(data = tab_map2,aes(x=follow),color="blue") +
- scale_x_log10(limits=c(1e-03,1e+05))
+rec_ratios %>% ggplot(aes(x=rat_prev,color=state)) +
+ geom_density() +
+ scale_color_manual(values=c("blue","red"),labels=c("After filtering","Before filtering")) +
+ scale_x_log10(labels=fancy_scientific,limits=c(1e-03,1e+05)) +
+ labs(x="Ratio between the calculated and the known distance with the previous marker",
+ y="Density",
+ color="",
+ title="Fold change distance between adjacent markers") +
+ ggpubr::theme_classic2()
ggsave("distrib_freq_rec2.png")
```
diff --git a/article/distrib_freq_rec1.png b/article/distrib_freq_rec1.png
index a867c67600edd03866a373c18ec91b329dd2962d..473752e4ef34fde05ba9082865f7dd238ddba1be 100644
Binary files a/article/distrib_freq_rec1.png and b/article/distrib_freq_rec1.png differ
diff --git a/article/distrib_freq_rec2.png b/article/distrib_freq_rec2.png
index 2341e56126ce4873c68268d326b3c3e048b2b6d9..bea97aaa82e6ef5df65888457fda5215258841af 100644
Binary files a/article/distrib_freq_rec2.png and b/article/distrib_freq_rec2.png differ