Skip to content
GitLab
Commit eaac4f31 authored by mariefbourdon's avatar mariefbourdon
Browse files

220211 figures article

parent 186d1d08
Hide whitespace changes
Inline Side-by-side
.Rproj.user/9DAE6990/pcs/windowlayoutstate.pper
View file @ eaac4f31
{
"left" : {
"panelheight" : 632,
"splitterpos" : 164,
"panelheight" : 1402,
"splitterpos" : 363,
"topwindowstate" : "NORMAL",
"windowheight" : 670
"windowheight" : 1440
},
"right" : {
"panelheight" : 632,
"splitterpos" : 397,
"panelheight" : 1402,
"splitterpos" : 880,
"topwindowstate" : "NORMAL",
"windowheight" : 670
"windowheight" : 1440
}
}
\ No newline at end of file
.Rproj.user/9DAE6990/pcs/workbench-pane.pper
View file @ eaac4f31
{
"TabSet1" : 1,
"TabSet1" : 0,
"TabSet2" : 3,
"TabZoom" : {
}
......
article/.RData
View file @ eaac4f31
No preview for this file type
article/article_figures.Rmd
View file @ eaac4f31
......@@ -127,7 +127,7 @@ na_plot
Proportions of markers with more than 75% of missing genotypes:
```{r prop_missing}
tab2 %>% mutate(prop_NA=n_NA/176) %>% filter(prop_NA > 0.75) %>% summarise(p=n()/count(tab2)%>%pull())
tab2 %>% mutate(prop_NA=n_NA/176) %>% filter(prop_NA > 0.50) %>% summarise(p=n()/count(tab2)%>%pull())
```
......@@ -147,7 +147,15 @@ prop_plot <- tab2 %>% filter(n_NA<88) %>% filter(!chr %in% c("M","X","Y")) %>%
theme(aspect.ratio=0.8,
legend.position=c(0.8,0.8),
legend.title = element_blank()) +
theme(plot.title = element_text(hjust = 0.4,face="bold",size=14))
theme(plot.title = element_text(hjust = 0.4,face="bold",size=14)) +
geom_segment(x = 0.2, y = 0.82,
xend = 0.25, yend = 0.72,
color = "#1d91c0",
arrow = arrow(type="closed",length=unit(6,"points"))) +
geom_segment(x = 0.39, y = 0.62,
xend = 0.44, yend = 0.52,
color = "#1d91c0",
arrow = arrow(type="closed",length=unit(6,"points")))
prop_plot
```
### Grid
......@@ -165,7 +173,7 @@ grid <- plot_grid(na_plot,prop_plot,
grid
ggsave("figures/naprop.pdf",grid,width=5,height=8)
ggsave("figures/fig2.pdf",grid,width=5,height=8)
rm(na_plot,prop_plot)
```
......@@ -208,19 +216,16 @@ rm(allele,strains_allele)
```{r barplot}
none <- tab2 %>% nrow()
match <- tab2 %>% filter(exclude_match==0) %>% nrow()
poly <- tab2 %>% filter(exclude_match==0&exclude_poly==0) %>% nrow()
prop <- tab2 %>% filter(exclude_match==0&exclude_poly==0&exclude_prop==0) %>% nrow()
allele <- tab2 %>% filter(exclude_match==0&exclude_poly==0&exclude_prop==0&exclude_allele==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()
barplot_df <- tibble(
fct = c("none","match","poly","prop","allele"),
markers = c(none, match, poly, prop, allele)
)
functions_plot <- barplot_df %>% ggplot(aes(x=markers,y=fct)) +
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("allele", "prop", "poly","match","none")) +
scale_y_discrete(limits=c("prop","poly", "na", "allele","match","none")) +
theme(aspect.ratio=0.7) +
labs(title="Number of markers kept after each step",
x="Number of markers",
......@@ -229,8 +234,6 @@ functions_plot <- barplot_df %>% ggplot(aes(x=markers,y=fct)) +
theme(plot.title = element_text(hjust = 0.4,face="bold",size=14))
functions_plot
ggsave("figures/functions.pdf",functions_plot,width=5,height=3)
rm(none,allele,match,poly,prop,barplot_df)
```
......@@ -240,8 +243,8 @@ rm(none,allele,match,poly,prop,barplot_df)
```{r cross_before}
# filter at minima: remove non polymorphic and NA>0.5
tab_before <- mark_poly(stuart_tab)
tab_before <- mark_prop(tab_before,cross="F2",homo=0,hetero=0)
tab_before <- mark_na(stuart_tab)
tab_before <- mark_poly(tab_before)
# 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")
......@@ -263,7 +266,7 @@ load("files/cluster/newmap_before.rda")
# plot
plotMap(cross_before,newmap_before,shift=TRUE)
plotmap_before <- ~plotMap(cross_before,newmap_before,shift=TRUE,main="Before stuart")
plotmap_before <- ~plotMap(cross_before,newmap_before,shift=TRUE,main="Before stuart", ylab='')
```
### Before: plot genome scan
......@@ -295,7 +298,10 @@ pheno_before_plot <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05"
theme(legend.position = "none") +
ggtitle("")
pheno_before_plot
```
```{r before_scan}
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",
......@@ -305,8 +311,14 @@ pheno_before_zoom <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05"
theme(legend.position = "none") +
ggtitle("")
pheno_before_zoom
```
### Before figure
ggsave("figures/zoom_peak_13.pdf",pheno_before_zoom,width=3)
```{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)
```
......@@ -350,40 +362,39 @@ plotmap_after <- ~plotMap(cross_after,newmap_after,shift=TRUE,main="After stuart
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_after.csv")
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")
```
### After: plot estimated map 2
```{r after_map2}
# import cross
cross_after <- read.cross(format="csv",file="files/cluster2/cross_after.csv",
cross_after2 <- read.cross(format="csv",file="files/cluster2/cross_after2.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
# load rda with newmap
load("files/cluster2/newmap_after.rda")
load("files/cluster2/newmap_after2.rda")
# plot
plotMap(cross_after,newmap_after,shift=TRUE)
plotmap_after <- ~plotMap(cross_after,newmap_after,shift=TRUE,main="After stuart")
plotMap(cross_after2,newmap_after2,shift=TRUE)
plotmap_after <- ~plotMap(cross_after2,newmap_after2,shift=TRUE,main="After stuart")
```
### After: plot genome scan
```{r after_scan}
# load rda with perm
load("files/cluster2/after_1000p.rda")
# load("files/after_1000p.rda")
load("files/cluster2/after_1000p2.rda")
# calculate thresholds
threshold_after <- summary(after_1000p,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63
threshold_after <- summary(after_1000p2,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63
# scanone
cross_after <- calc.genoprob(cross_after, step=2.0, off.end=0.0,
cross_after <- calc.genoprob(cross_after2, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
pheno_after <- scanone(cross=cross_after, 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")
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")
summary(pheno_after)
# Plot
......@@ -397,13 +408,12 @@ pheno_after_plot <- qtl_plot(pheno_after,lod=data.frame(group = c("alpha=0.05",
ggtitle("")
pheno_after_plot
```
### Grid
### Grid after
```{r grid fig1, fig.height = 7, fig.width = 13, fig.align = "center"}
grid <- plot_grid(as_grob(plotmap_before),as_grob(plotmap_after),pheno_before_plot,pheno_after_plot,ncol=2,labels=c("A","B","C","D"),label_size=20)
grid
grid_after <- plot_grid(as_grob(plotmap_after),pheno_after_plot,ncol=1,labels=c("A","B"),label_size=20)
ggsave("figures/before_after.pdf",grid,width=11,height = 7)
ggsave("figures/fig3.pdf",grid_after,width=7,height = 10)
```
# Differences between real and reference genotypes for parental strains
......@@ -429,8 +439,7 @@ knitr::kable(table_dif)
# Pheno data format
```{r pheno}
format_pheno <- phenos[1:6,]
write_csv(format_pheno,file="sup/tableS1.csv")
print(xtable::xtable(format_pheno, type = "latex"), file = "tables/tab_alleles.tex",include.rownames=FALSE)
```
article/files/cluster/before_after.R
View file @ eaac4f31
##################################################################
########### before_after.R
########### R script for stuart article figures
##################################################################
## Goal and data
###########################
# The goal of this script is to calculate est.map and permutations of cross_before et cross_after
#
# It imports the csv files "cross_before.csv" and "cross_after.csv"
## Import packages
library(qtl)
##############################
## Before
##############################
## read cross
cross_before <- read.cross(format="csv",file="./cross_before.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
## calc geno prob
cross_before <- calc.genoprob(cross_before, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
# newmap_before <- est.map(cross=cross_before,error.prob=0.01)
# save(newmap_before,file="./newmap_before.rda")
## calculate 1000p
before_1000p <- 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", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE)
save(before_1000p,file="./before_1000p.rda")
##############################
## After
##############################
## read cross
cross_after <- read.cross(format="csv",file="./cross_after.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
## calc geno prob
cross_after <- calc.genoprob(cross_after, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
newmap_after <- est.map(cross=cross_after,error.prob=0.01)
save(newmap_after,file="./newmap_after.rda")
## calculate 1000p
after_1000p <- scanone(cross=cross_after,
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", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE)
##################################################################
########### before_after.R
########### R script for stuart article figures
##################################################################
## Goal and data
###########################
# The goal of this script is to calculate est.map and permutations of cross_before et cross_after
#
# It imports the csv files "cross_before.csv" and "cross_after.csv"
## Import packages
library(qtl)
##############################
## Before
##############################
## read cross
cross_before <- read.cross(format="csv",file="./cross_before.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
## calc geno prob
cross_before <- calc.genoprob(cross_before, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
newmap_before <- est.map(cross=cross_before,error.prob=0.01)
save(newmap_before,file="./newmap_before.rda")
## calculate 1000p
before_1000p <- 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", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE)
save(before_1000p,file="./before_1000p.rda")
##############################
## After
##############################
## read cross
cross_after <- read.cross(format="csv",file="./cross_after.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
## calc geno prob
cross_after <- calc.genoprob(cross_after, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
newmap_after <- est.map(cross=cross_after,error.prob=0.01)
save(newmap_after,file="./newmap_after.rda")
## calculate 1000p
after_1000p <- scanone(cross=cross_after,
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", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE)
save(after_1000p,file="./after_1000p.rda")
\ No newline at end of file
article/files/cluster/before_after.err
View file @ eaac4f31
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--Omitting 1120 male heterozygote genotypes on the X chromosome.
--Omitting 1198 male heterozygote genotypes on the X chromosome.
2: In fixXgeno.f2(cross, alleles) :
--There appear to be some individuals of each cross direction, but "pgm" is not provided.
Check the X chr genotype data and include a "pgm" column in the phenotype data.
......@@ -12,7 +12,7 @@ Warning message:
In matchchr(chr, names(x$geno)) : Chr Y not found
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--Omitting 236 male heterozygote genotypes on the X chromosome.
--Omitting 239 male heterozygote genotypes on the X chromosome.
2: In fixXgeno.f2(cross, alleles) :
--There appear to be some individuals of each cross direction, but "pgm" is not provided.
Check the X chr genotype data and include a "pgm" column in the phenotype data.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment