data 2-3

.Rproj.user/shared/notebooks/paths

-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/DESCRIPTION="5140B408"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/geno_strains.R="EBA82473"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/mark_allele.R="A0A182D3"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/mark_estmap.R="57E1825F"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/mark_match.R="FFA1BE52"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/mark_prop.R="30FA9E8C"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/tab_mark.R="F0B2417"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/R/write_rqtl.R="F07CE16C"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/README.Rmd="297B98D"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/man/mark_na.Rd="C8CA0D4F"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/man/mark_prop.Rd="893B273D"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/vignettes/stuaRt.Rmd="DA6206CB"
-/Users/mariebourdon/Documents/PhD/stuart_package/stuart/vignettes/stuart.Rmd="54D793B9"
+/home/marie/Documents/HB_CC011/F2/Geno/Maestro/F2_geno_estmap1.R="42B0AE06"
+/home/marie/Documents/stuart/stuart_package/stuart/R/mark_estmap.R="F759F000"
+/home/marie/Documents/stuart/stuart_package/stuart/vignettes/stuart.Rmd="01396A37"

article/article_figures.Rmd

@@ -207,7 +207,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", ylab='')
+plotmap_before <- ~plotMap(cross_before,newmap_before,shift=TRUE,main="Before stuart", ylab='^')
 ```
 
 ### Before: plot genome scan
@@ -235,7 +235,8 @@ pheno_before_plot <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05"
          ylab="LOD score",
          title="QTL mapping",
          x.label = element_blank(),
-         size=0.6) +
+         size=0.6,
+         strip.pos="bottom") +
     theme(legend.position = "none") +
     ggtitle("")
 pheno_before_plot

article/data2/cluster/before_after.R

+##################################################################
+########### before_after.R
+########### R script for stuart article
+##################################################################
+
+## 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/data2/cluster/before_after.err

+Warning messages:
+1: In fixXgeno.f2(cross, alleles) :
+   --Omitting 382 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!
+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.
+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().
+Warning message:
+In matchchr(chr, names(x$geno)) : Chr Y not found

article/data2/cluster/before_after.out

+ --Read the following data:
+	 94  individuals
+	 3408  markers
+	 4  phenotypes
+ --Cross type: f2 
+ --Read the following data:
+	 94  individuals
+	 1293  markers
+	 4  phenotypes
+ --Cross type: f2 

article/data2/cluster2/after.R

+##################################################################
+########### after.R
+########### R script for stuart article 
+##################################################################
+
+## Goal and data
+###########################
+
+# The goal of this script is to calculate est.map and permutations of cross_after
+# 
+# It imports the csv file "cross_after.csv"
+
+## Import packages
+library(qtl)
+
+##############################
+
+## read cross
+cross_after <- 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, 
+                              error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")
+
+## calculate est map
+newmap_after2 <- est.map(cross=cross_after,error.prob=0.01)
+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"), 
+                        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

article/data2/cluster2/after.sh

+#!/bin/sh
+
+#SBATCH -J "after"
+#SBATCH --qos=normal
+#SBATCH --mem=1000
+#SBATCH -o after.out -e after.err
+#SBATCH --mail-type=END  --mail-user=mbourdon@pasteur.fr
+
+Rscript after.R  || exit 1
+
+
+exit 0

article/data2/data2.Rmd

+---
+title: "data2.Rmd"
+author: "Marie Bourdon"
+date: "01/06/2022"
+output: html_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+library(stuart)
+library(magrittr)
+library(readr)
+library(stringr)
+library(qtl)
+```
+
+## Load
+```{r}
+genos <- read_csv("geno_data2.csv",show_col_types = FALSE) #genotypes of F2s
+annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) #annotation file for miniMUGA
+phenos <- read_csv("pheno_data2.csv",show_col_types = FALSE) #phenotypes of F2s
+parents <- read_csv("parents_data2.csv",show_col_types = FALSE) #genotypes of parental strains (genotyped with F2s)
+strns_ref <- read_csv("ref_geno_data2.csv",show_col_types = FALSE) #reference genotypes of parental strains
+```
+
+```{r}
+tab <- tab_mark(genos,annot_mini,"cM_cox")
+```
+
+## Before: creation of Rqtl csv file 
+
+```{r cross_before}
+# filter at minima: remove non polymorphic and NA>0.5 
+tab_before <- mark_na(tab)
+tab_before <- mark_poly(tab_before)
+
+#join with annotation file from miniMUGA
+strns_ref <- strns_ref %>% select(name,StrainA,StrainB) %>% right_join(annot_mini,.,by=c("marker"="name")) %>% select(c(marker,chr,cM_cox,StrainA,StrainB))
+
+# create rqtl csv file
+cross_before <- write_rqtl(geno=genos,pheno=phenos,tab=tab_before,ref=strns_ref,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster/cross_before.csv")
+
+# import cross
+cross_before <- read.cross(format="csv",file="cluster/cross_before.csv",
+                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
+
+
+load("cluster/newmap_before.rda")
+plotMap(cross_before,newmap_before,shift=TRUE)
+
+```
+
+## create file for parental strains genotyped
+
+```{r}
+#our genotypes
+
+#create tibble with individivual names
+parental_strains <- tibble::tibble(StrainA = c("StrainsA_1","StrainsA_2"),
+                                   StrainB = c("StrainsB_1","StrainsB_2"))
+
+
+#create data frame with geno_strains
+strains <- geno_strains(annot=annot_mini,geno=parents,
+                        strn=parental_strains,cols=c("chr","cM_cox"))
+rm(parental_strains)
+
+#summary
+summary(strains)
+```
+
+
+## Use of stuart's functions
+
+```{r}
+tab2 <- mark_match(stuart_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))
+tab2 <- mark_allele(tab2,ref=strains,cross="F2",par1="StrainA",par2="StrainB")
+```
+
+
+
+### estmap
+
+```{r}
+# create rqtl csv file
+write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster/cross_after.csv")
+
+# import cross
+cross_after <- read.cross(format="csv",file="cluster/cross_after.csv",
+                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)
+
+
+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")
+```
+
+```{r estmap}
+
+```
+
+
+## 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)
+```
\ No newline at end of file

article/data2/data2.Rproj

+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX