Commit d95c71ca authored by mariefbourdon's avatar mariefbourdon
Browse files

data 2-3-4

parent 5c5d1e1b
No preview for this file type
No preview for this file type
##################################################################
########### before_after.R
########### R script for stuart article
########### R script for stuart article figures
##################################################################
## Goal and data
......@@ -23,7 +23,7 @@ cross_before <- read.cross(format="csv",file="./cross_before.csv",
## 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")
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
newmap_before <- est.map(cross=cross_before,error.prob=0.01)
......@@ -32,8 +32,8 @@ 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)
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")
##############################
......@@ -46,7 +46,7 @@ cross_after <- read.cross(format="csv",file="./cross_after.csv",
## 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")
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
newmap_after <- est.map(cross=cross_after,error.prob=0.01)
......@@ -55,6 +55,6 @@ 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)
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
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--Omitting 382 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.
"pgm" was inferred (probably poorly).
3: In summary.cross(cross) :
Some markers at the same position on chr 1,2,3,5,6,7,8,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!
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 matchchr(chr, names(x$geno)) : Chr Y not found
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--Omitting 57 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.
"pgm" was inferred (probably poorly).
3: In summary.cross(cross) :
Some markers at the same position on chr 5,7,8,10,11,12,13,14,15,17,19,X; use jittermap().
Some markers at the same position on chr 1,2,5,6,7,8,9,10,11,12,13,14,15,17,18,19,X; use jittermap().
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:
94 individuals
3408 markers
176 individuals
3343 markers
4 phenotypes
--Cross type: f2
--Read the following data:
94 individuals
1293 markers
176 individuals
2241 markers
4 phenotypes
--Cross type: f2
This diff is collapsed.
##################################################################
########### after.R
########### R script for stuart article
########### R script for stuart article figures
##################################################################
## Goal and data
......@@ -16,11 +16,11 @@ library(qtl)
##############################
## read cross
cross_after <- read.cross(format="csv",file="./cross_after2.csv",
cross_after2 <- read.cross(format="csv",file="./cross_after2.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,
cross_after2 <- calc.genoprob(cross_after2, step=2.0, off.end=0.0,
error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
## calculate est map
......@@ -29,7 +29,7 @@ save(newmap_after2,file="./newmap_after2.rda")
## calculate 1000p
after_1000p2 <- 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"),
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)
save(after_1000p2,file="./after_1000p2.rda")
\ No newline at end of file
save(after_1000p2,file="./after_1000p2.rda")
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--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.
"pgm" was inferred (probably poorly).
3: In summary.cross(cross) :
Some markers at the same position on chr 1,2,5,6,7,8,9,10,11,12,13,14,15,17,18,19,X; use jittermap().
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:
176 individuals
2238 markers
4 phenotypes
--Cross type: f2
This diff is collapsed.
......@@ -12,6 +12,8 @@ library(magrittr)
library(readr)
library(stringr)
library(qtl)
library(ggplot2)
source("../files/QTL_plot.R")
```
## Load
......@@ -47,7 +49,38 @@ cross_before <- read.cross(format="csv",file="cluster/cross_before.csv",
load("cluster/newmap_before.rda")
plotMap(cross_before,newmap_before,shift=TRUE)
```
### Before: plot genome scan
```{r before_scan}
# load rda with perm
load("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(),
size=0.6,
strip.pos="bottom") +
theme(legend.position = "none") +
ggtitle("")
pheno_before_plot
```
## create file for parental strains genotyped
......@@ -73,7 +106,7 @@ summary(strains)
## Use of stuart's functions
```{r}
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 <- mark_match(tab,ref=strains)
tab2 <- mark_poly(tab2)
tab2 <- mark_na(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))
......@@ -102,37 +135,74 @@ tab2 <- mark_estmap(tab=tab2,map=newmap_after,annot=annot_mini)
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster2/cross_after2.csv")
```
```{r estmap}
### After: plot estimated map 2
```{r after_map2}
# import cross
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")
# plot
plotMap(cross_after2,newmap_after2,shift=TRUE)
plotmap_after <- ~plotMap(cross_after2,newmap_after2,shift=TRUE,main="After stuart")
```
```{r after_scan}
# load rda with perm
load("cluster2/after_1000p2.rda")
# 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
# scanone
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_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
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(),
size=0.6) +
theme(legend.position = "none") +
ggtitle("")
pheno_after_plot
```
## Number of markers kept after each function
```{r barplot}
## TO DO WHEN ESTMAP OK
# 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)
# TO DO WHEN ESTMAP OK
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,functions_df)
```
\ No newline at end of file
No preview for this file type
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