Commit 6a644002 authored by mariefbourdon's avatar mariefbourdon
Browse files

data2-3 ok

parent 158fb22c
No preview for this file type
......@@ -177,6 +177,7 @@ grid <- plot_grid(na_plot,prop_plot,
grid
ggsave("figures/fig2.pdf",grid,width=5,height=8)
ggsave("figures/fig2.png",grid,width=5,height=8)
rm(na_plot,prop_plot)
```
......@@ -266,6 +267,7 @@ pheno_before_zoom
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)
ggsave("figures/fig1.png",grid_before,width=7,height = 15)
```
......@@ -362,6 +364,7 @@ pheno_after_plot
grid_after <- plot_grid(as_grob(plotmap_after),pheno_after_plot,ncol=1,labels=c("A","B"),label_size=20)
ggsave("figures/fig3.pdf",grid_after,width=7,height = 10)
ggsave("figures/fig3.png",grid_after,width=7,height = 10)
```
# Differences between real and reference genotypes for parental strains
......
No preview for this file type
##################################################################
########### 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
Warning messages:
1: In fixXgeno.f2(cross, alleles) :
--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,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 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
3343 markers
4 phenotypes
--Cross type: f2
--Read the following data:
176 individuals
2241 markers
4 phenotypes
--Cross type: f2
#!/bin/sh
#SBATCH -J "before_after"
#SBATCH --qos=normal
#SBATCH --mem=1000
#SBATCH -o before_after.out -e before_after.err
#SBATCH --mail-type=END --mail-user=mbourdon@pasteur.fr
Rscript before_after.R || exit 1
exit 0
This diff is collapsed.
This diff is collapsed.
......@@ -130,39 +130,22 @@ load("cluster/newmap_after.rda")
plotMap(cross_after,newmap_after,shift=TRUE)
tab2 <- mark_estmap(tab=tab2,map=newmap_after,annot=annot_mini)
# create new rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster2/cross_after2.csv")
```
### 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")
load("cluster/after_1000p.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
threshold_after <- summary(after_1000p,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,
cross_after <- calc.genoprob(cross_after, 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")
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")
summary(pheno_after)
# Plot
......
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