Commit ae7929d3 authored by mariefbourdon's avatar mariefbourdon
Browse files

220627 figures 3

parent 3c831c51
No preview for this file type
......@@ -283,25 +283,21 @@ colnames(peak1_tab) <- c("marker","chr","pos","allele\n1","allele\n2",
"n\nHM1","n\nHM2","n\nHT","n\nNA")
colnames(peak3_tab) <- c("marker","chr","pos","allele\n1","allele\n2",
"n\nHM1","n\nHM2","n\nHT","n\nNA")
colnames(peak5_tab) <- c("marker","chr","pos","allele\n1","allele\n2",
"n\nHM1","n\nHM2","n\nHT","n\nNA")
narrow_grid <- ggdraw() +
draw_plot(pheno_before_annot_data2,0,.75,.5,.25) +
draw_plot(pheno_before_annot_data3,.5,.75,.5,.25) +
draw_plot(peak1,0,.5,.3,.25) +
draw_plot(tableGrob(peak1_tab,rows=NULL),.3,.5,.7,.25) +
draw_plot(peak3,0,.25,.3,.25) +
draw_plot(tableGrob(peak3_tab,rows=NULL),.3,.25,.7,.25) +
draw_plot(peak5,0,0,.3,.25) +
draw_plot(tableGrob(peak5_tab,rows=NULL),.3,0,.7,.25) +
draw_plot_label(c("A","B","C","D","E","F","G","H"),
c(0,.5,0,.3,0,.3,.5,.7),c(1,1,.75,.75,.5,.5,.25,.25))
draw_plot(pheno_before_annot_data2,0,.7,.5,.3) +
draw_plot(pheno_before_annot_data3,.5,.7,.5,.3) +
draw_plot(peak1,0,.35,.3,.3) +
draw_plot(tableGrob(peak1_tab,rows=NULL),.3,.35,.7,.3) +
draw_plot(peak3,0,0,.3,.3) +
draw_plot(tableGrob(peak3_tab,rows=NULL),.3,0,.7,.3) +
draw_plot_label(c("A","B","C","D","E","F"),
c(0,.5,0,.3,0,.3),c(1,1,.65,.65,.3,.3))
narrow_grid
ggsave(narrow_grid,file="figures/figS3prev.pdf",width=10,height=17)
ggsave(narrow_grid,file="figures/figS3prev.png",width=10,height=17)
ggsave(narrow_grid,file="figures/figS3prev.pdf",width=10,height=14)
ggsave(narrow_grid,file="figures/figS3prev.png",width=10,height=14)
rm(peak1,peak1_tab,peak2,peak2_tab,peak3,peak3_tab,peak5,peak5_tab)
```
......@@ -336,22 +332,18 @@ ggsave("figures/figS5.pdf",grid_scans,width=14,height = 5)
## Comparision of non mendelian/wrong alleles markers between the 4 datasets
```{r}
load("data2/tab2_data2.rda")
load("data3/tab2_data3.rda")
load("data4/tab2_data4.rda")
A REFAIRE!!!!
```{r}
colnames(tab2_data2) <- c(colnames(tab2_data2)[1:3],paste0(colnames(tab2_data2)[4:ncol(tab2_data2)],".2"))
colnames(tab2_data3) <- c(colnames(tab2_data3)[1:3],paste0(colnames(tab2_data3)[4:ncol(tab2_data3)],".3"))
colnames(tab2_data4) <- c(colnames(tab2_data4)[1:3],paste0(colnames(tab2_data4)[4:ncol(tab2_data4)],".4"))
tab_alldata <- full_join(tab2,tab2_data2) %>% full_join(.,tab2_data3) %>% full_join(.,tab2_data4)
tab_alldata <- full_join(tab2,tab2_data2) %>% full_join(.,tab2_data3)
#keep all markers with non mendelian prop in at least 1 dataset
prop_rem <- tab_alldata %>% filter((exclude_poly==0&exclude_prop==1) |
(exclude_poly.2==0&exclude_prop.2==1) |
(exclude_poly.3==0&exclude_prop.3==1) |
(exclude_poly.4==0&exclude_prop.4==1))
prop_rem <- tab_alldata %>% filter((exclude_poly==0&exclude_prop==1) | (exclude_poly==0&exclude_allele==1) |
(exclude_poly.2==0&exclude_prop.2==1) | (exclude_poly.2==0&exclude_allele.2==1) |
(exclude_poly.3==0&exclude_prop.3==1) | (exclude_poly.3==0&exclude_allele.3==1))
nrow(prop_rem)
......@@ -359,25 +351,24 @@ nrow(prop_rem)
prop_rem %<>% mutate(exclpr_data1=case_when(exclude_poly==0&exclude_prop==1 ~ 1, T ~ 0)) %>%
mutate(exclpr_data2=case_when(exclude_poly.2==0&exclude_prop.2==1 ~ 1, T ~ 0)) %>%
mutate(exclpr_data3=case_when(exclude_poly.3==0&exclude_prop.3==1 ~ 1, T ~ 0)) %>%
mutate(exclpr_data4=case_when(exclude_poly.4==0&exclude_prop.4==1 ~ 1, T ~ 0)) %>%
mutate(sompr=exclpr_data1+exclpr_data2+exclpr_data3+exclpr_data4)
mutate(sompr=exclpr_data1+exclpr_data2+exclpr_data3)
#note if marker retained in a data set
prop_rem %<>% mutate(excl_data1=case_when((exclude_match+exclude_poly+exclude_na+exclude_prop+exclude_allele+exclude_estmap)==0 ~ 0, T ~ 1)) %>%
mutate(excl_data2=case_when((exclude_match.2+exclude_poly.2+exclude_na.2+exclude_prop.2+exclude_allele.2+exclude_estmap.2)==0 ~ 0, T ~ 1)) %>%
mutate(excl_data3=case_when((exclude_match.3+exclude_poly.3+exclude_na.3+exclude_prop.3+exclude_allele.3+exclude_estmap.3)==0 ~ 0, T ~ 1)) %>%
mutate(excl_data4=case_when((exclude_match.4+exclude_poly.4+exclude_na.4+exclude_prop.4+exclude_allele.4+exclude_estmap.4)==0 ~ 0, T ~ 1)) %>%
mutate(som=excl_data1+excl_data2+excl_data3+excl_data4)
mutate(som=excl_data1+excl_data2+excl_data3)
#number of markers ok in a least 1 dataset
prop_rem %>% filter(som<4) %>% nrow()
prop_rem %>% filter(som<3) %>% nrow()
#number of markers always removed
prop_rem %<>% mutate(notpoly=case_when(som==4 & (exclude_poly+exclude_poly.2+exclude_poly.3+exclude_poly.4==3) ~ 0,
som==4 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3+exclude_poly.4==2) ~ 0,
som==4 & (exclude_poly+exclude_poly.2+exclude_poly.3+exclude_poly.4==0) ~ 2,
som==4 & (exclude_poly+exclude_poly.2+exclude_poly.3+exclude_poly.4<3) & (exclude_poly+exclude_poly.2+exclude_poly.3+exclude_poly.4>0) ~ 1,
som==4 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3+exclude_poly.4<2) & (exclude_poly+exclude_poly.3+exclude_poly.4>0) ~ 1))
prop_rem %<>% mutate(notpoly=case_when(som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3==2) ~ 0,
som==3 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3==1) ~ 0,
som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3==0) ~ 2,
som==3 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3==0) ~ 2,
som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3<2) & (exclude_poly+exclude_poly.2+exclude_poly.3>0) ~ 1,
som==3 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3<1) & (exclude_poly+exclude_poly.3>0) ~ 1))
prop_rem %>% group_by(notpoly) %>% summarise(n=n())
```
......@@ -388,12 +379,14 @@ prop_rem %>% group_by(notpoly) %>% summarise(n=n())
#keep all markers with non mendelian prop in at least 1 dataset
map_rem <- tab_alldata %>% filter(exclude_estmap==1 |
exclude_estmap.2==1 |
exclude_estmap.3==1 |
exclude_estmap.4==1)
exclude_estmap.3==1)
map_rem %<>% mutate(som=case_when(exclude_match+exclude_poly+exclude_na+exclude_prop+exclude_allele+exclude_estmap==0 ~ 0,
exclude_match.2+exclude_poly.2+exclude_na.2+exclude_prop.2+exclude_allele.2+exclude_estmap.2==0 ~ 0,
exclude_match.3+exclude_poly.3+exclude_na.3+exclude_prop.3+exclude_allele.3+exclude_estmap.3==0 ~ 0,))
```
11/15 are okay in another cross.
5/9 are okay in another cross.
No preview for this file type
Warning messages:
1: In getsex(cross) :
3 individuals with missing sex; assuming they're female
2: In fixXgeno.f2(cross, alleles) :
1: In fixXgeno.f2(cross, alleles) :
--Omitting 332 male heterozygote genotypes on the X chromosome.
3: In fixXgeno.f2(cross, alleles) :
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.
"pgm" was inferred (probably poorly).
4: In getsex(object) :
3 individuals with missing sex; assuming they're female
5: In getsex(object) :
3 individuals with missing sex; assuming they're female
6: In summary.cross(cross) :
3: In summary.cross(cross) :
Some markers at the same position on chr 1,2,3,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,X; use jittermap().
Warning message:
In est.map(cross = cross_before, error.prob = 0.01) : Didn't converge!
Warning messages:
1: In matchchr(chr, names(x$geno)) : Chr Y not found
2: In getsex(cross) :
3 individuals with missing sex; assuming they're female
3: In checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, :
Dropping 3 individuals with missing phenotypes.
In matchchr(chr, names(x$geno)) : Chr Y not found
Warning messages:
1: In getsex(cross) :
3 individuals with missing sex; assuming they're female
2: In fixXgeno.f2(cross, alleles) :
1: In fixXgeno.f2(cross, alleles) :
--Omitting 27 male heterozygote genotypes on the X chromosome.
3: 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.
"pgm" was inferred (probably poorly).
4: In getsex(object) :
3 individuals with missing sex; assuming they're female
5: In getsex(object) :
3 individuals with missing sex; assuming they're female
6: In summary.cross(cross) :
2: In summary.cross(cross) :
Some markers at the same position on chr 1,2,3,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,X; use jittermap().
Warning messages:
1: In matchchr(chr, names(x$geno)) : Chr Y not found
2: In getsex(cross) :
3 individuals with missing sex; assuming they're female
3: In checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, :
Dropping 3 individuals with missing phenotypes.
Warning message:
In est.map(cross = cross_after, error.prob = 0.01) : Didn't converge!
Warning message:
In matchchr(chr, names(x$geno)) : Chr Y not found
--Read the following data:
92 individuals
2774 markers
89 individuals
2762 markers
5 phenotypes
--Cross type: f2
--Read the following data:
92 individuals
2630 markers
89 individuals
2627 markers
5 phenotypes
--Cross type: f2
......@@ -215,8 +215,10 @@ write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="StrainA",par2="Str
cross_after2 <- read.cross(format="csv",file="cluster2/cross_after2.csv",
genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
# load rda with newmap
load("cluster2/newmap_after2.rda")
cross_after2 <- calc.genoprob(cross_after2, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
newmap_after2 <- est.map(cross=cross_after2,error.prob=0.01)
# plot
plotMap(cross_after2,newmap_after2,shift=TRUE)
......@@ -229,7 +231,9 @@ plotmap_after_data2 <- ~plotMap(cross_after_data3,newmap_after_data3,shift=TRUE,
```{r after_scan}
# load rda with perm
load("cluster2/after_1000p2.rda")
after_1000p2 <- 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"),
pheno.col="Pheno", model="normal", method="em", n.perm=1000, perm.Xsp=FALSE, verbose=FALSE)
# calculate thresholds
threshold_after <- summary(after_1000p2,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63
......@@ -361,9 +365,9 @@ chrs <- c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "
ann_dat_text<-data.frame(
chr=factor(chrs,
levels=chrs),
lod=c(NA,NA,NA,NA,11,NA,NA,NA,NA,13,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
label=c(NA,NA,NA,NA,"p4",NA,NA,NA,NA,"p5",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
x=c(NA,NA,NA,NA,20,NA,NA,NA,NA,27,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
lod=c(NA,NA,NA,NA,11,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
label=c(NA,NA,NA,NA,"p4",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
x=c(NA,NA,NA,NA,20,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
)
pheno_before_annot_data3 <- pheno_before_plot + geom_text(
......@@ -438,67 +442,6 @@ geno_plot4 <- pgmap %>% filter(pos > 1 & pos < 30) %>%
geno_plot4
```
### Peak 5
```{r peak5_zoom}
peak5 <- 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="10",
rug=TRUE) +
theme(legend.position = "none",
strip.background = element_blank(),
strip.text.x = element_blank()) +
xlab("Position (cM)") +
xlim(c(22,42)) +
ggtitle("Peak 5: 1 marker and 1 pseudomarker")
peak5
```
1 peak on 1 marker and one pseudomarker : S2C102505843 and c10.loc30 postionned between SAH102097335 and S2C102505843.
Here are the infos on genotype counts for these markers:
```{r summary_geno_peak5}
peak5_tab <- tab_before %>% filter(marker %in% c("SAH102097335","S2C102505843")) %>% select(marker:n_NA)
peak5_tab
```
For S2C102505843, all individuals except 1 are homozygous so this marker should be removed. The proportions for SAH102097335 seem correct.
Graph:
```{r geno_plot_peak5}
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:S6J105101847,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")
geno_plot5 <- 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,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 = 50)))
geno_plot5
```
## Phenotype distributions
```{r pheno_distrib}
......@@ -551,7 +494,7 @@ prop_plot
```{r}
tab2_data3 <- tab2
save(pheno_before_annot_data3,pheno_data3,peak5,peak5_tab,tab2_data3,
save(pheno_before_annot_data3,pheno_data3,tab2_data3,
rec_ratios_before_data3,rec_ratios_after_data3,
plotmap_before_data3,plotmap_after_data3,
cross_before_data3,newmap_before_data3,cross_after_data3,newmap_after_data3,
......
No preview for this file type
No preview for this file type
No preview for this file type
This diff is collapsed.
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