diff --git a/.Rhistory b/.Rhistory index 647405bfe5b9c74f5ed909725945acdcdb2bd40f..338cca9ae0df8999c2e5d3491375f5a2f3275512 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,168 +1,3 @@ -bad_markers <- tab_map %>% filter(exclude == 1) %>% pull(marker) -} -mark_estmap <- function(map,dist=100){ -#initialize variables -marker <- c() -chr <- c() -pos <- c() -place <- c() -previous <- c() -follow <- c() -#get information in newmap -for(i in names(map)){ -marker <- c(marker,names(map[[i]])) -chr <- c(chr,rep(i,times=length(map[[i]]))) -pos <- c(pos,unname(map[[i]])) -place <- c(place,"first",rep("middle",times=(length(map[[i]])-2)),"last") -prev <- c(NA,unname(map[[i]])[1:length(map[[i]])-1]) -previous <- c(previous,prev) -fol <- c(unname(map[[i]])[2:length(map[[i]])],NA) -follow <- c(follow,fol) -} -tab_map <- tibble(marker = marker, -chr = chr, -place = place, -pos = pos, -previous = pos-previous, -follow = follow-pos) -tab_map <- tab_map %>% mutate(exclude=case_when(place == "first" & follow > dist ~ 1, -place == "middle" & previous > dist & follow > dist ~ 1, -place == "last" & previous > dist ~ 1, -T ~ 0)) -bad_markers <- tab_map %>% filter(exclude == 1) %>% pull(marker) -} -mark_estmap() -mark_estmap(map=newmap_after) -mark_estmap <- function(map,dist=100){ -#initialize variables -marker <- c() -chr <- c() -pos <- c() -place <- c() -previous <- c() -follow <- c() -#get information in newmap -for(i in names(map)){ -marker <- c(marker,names(map[[i]])) -chr <- c(chr,rep(i,times=length(map[[i]]))) -pos <- c(pos,unname(map[[i]])) -place <- c(place,"first",rep("middle",times=(length(map[[i]])-2)),"last") -prev <- c(NA,unname(map[[i]])[1:length(map[[i]])-1]) -previous <- c(previous,prev) -fol <- c(unname(map[[i]])[2:length(map[[i]])],NA) -follow <- c(follow,fol) -} -tab_map <- tibble(marker = marker, -chr = chr, -place = place, -pos = pos, -previous = pos-previous, -follow = follow-pos) -tab_map <- tab_map %>% mutate(exclude=case_when(place == "first" & follow > dist ~ 1, -place == "middle" & previous > dist & follow > dist ~ 1, -place == "last" & previous > dist ~ 1, -T ~ 0)) -bad_markers <- tab_map %>% filter(exclude == 1) %>% pull(marker) -return(bad_markers) -} -mark_estmap(map=newmap_after) -mark_estmap <- function(map,dist=100){ -#initialize variables -marker <- c() -chr <- c() -pos <- c() -place <- c() -previous <- c() -follow <- c() -#get information in newmap -for(i in names(map)){ -marker <- c(marker,names(map[[i]])) -chr <- c(chr,rep(i,times=length(map[[i]]))) -pos <- c(pos,unname(map[[i]])) -place <- c(place,"first",rep("middle",times=(length(map[[i]])-2)),"last") -prev <- c(NA,unname(map[[i]])[1:length(map[[i]])-1]) -previous <- c(previous,prev) -fol <- c(unname(map[[i]])[2:length(map[[i]])],NA) -follow <- c(follow,fol) -} -tab_map <- tibble(marker = marker, -chr = chr, -place = place, -pos = pos, -previous = pos-previous, -follow = follow-pos) -tab_map <- tab_map %>% mutate(exclude=case_when(place == "first" & follow > dist ~ 1, -place == "middle" & previous > dist & follow > dist ~ 1, -place == "last" & previous > dist ~ 1, -T ~ 0)) -bad_markers <- tab_map %>% filter(exclude == 1) %>% pull(marker) -return(bad_markers) -} -mark_estmap(map=newmap_after,dist=1) -mark_estmap <- function(map,dist=100){ -#initialize variables -marker <- c() -chr <- c() -pos <- c() -place <- c() -previous <- c() -follow <- c() -#get information in newmap -for(i in names(map)){ -marker <- c(marker,names(map[[i]])) -chr <- c(chr,rep(i,times=length(map[[i]]))) -pos <- c(pos,unname(map[[i]])) -place <- c(place,"first",rep("middle",times=(length(map[[i]])-2)),"last") -prev <- c(NA,unname(map[[i]])[1:length(map[[i]])-1]) -previous <- c(previous,prev) -fol <- c(unname(map[[i]])[2:length(map[[i]])],NA) -follow <- c(follow,fol) -} -tab_map <- tibble(marker = marker, -chr = chr, -place = place, -pos = pos, -previous = pos-previous, -follow = follow-pos) -tab_map <- tab_map %>% mutate(exclude=case_when(place == "first" & follow > dist ~ 1, -place == "middle" & previous > dist & follow > dist ~ 1, -place == "last" & previous > dist ~ 1, -T ~ 0)) -bad_markers <- tab_map %>% filter(exclude == 1) %>% pull(marker) -return(bad_markers) -} -mark_estmap(map=newmap_after,dist=100) -calc.errorlod(crosss) -calc.errorlod(cross) -dplyr::tibble -dplyr::select() -roxygen2::roxygenise() -rm(list = c("mark_estmap")) -roxygen2::roxygenise() -stuart_cross <- read.cross(format="csv",file="rqtl_file.csv", -genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE) -stuart_cross <- cross -View(stuart_cross) -View(cross) -View(cross) -View(cross) -usethis::use_data(stuart_cross) -data("stuart_cross") -View(stuart_cross) -roxygen2::roxygenise() -library(dplyr) -library(stuart) -data(stuart_cross) -data(stuart_cross) -summary(stuart_cross) -str(stuart_cross) -summary(stuart_cross) -knitr::opts_chunk$set( -collapse = TRUE, -comment = "#>" -) -library(dplyr) -library(stuart) annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) data(genos) summary(genos) @@ -510,3 +345,168 @@ genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", data(stuart_tab) summary(stuart_tab) View(stuart_tab) +knitr::opts_chunk$set( +collapse = TRUE, +comment = "#>" +) +library(dplyr) +library(stuart) +annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) +data(genos) +summary(genos) +data(phenos) +summary(phenos) +strains <- geno_strains(ref=annot_mini,geno=genos,par1=c("StrainsA_1","StrainsA_2"), +par2=c("StrainsB_1","StrainsB_2"),name1="parent1",name2="parent2") +head(strains) %>% print.data.frame() +genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", +"StrainsB_1","StrainsB_2")) +# how to use the function: +# stuart_tab <- tab_mark(genos) +data(stuart_tab) +summary(stuart_tab) +tab2 <- mark_match(stuart_tab,ref=strains) +tab2 %>% filter(exclude_match==1) %>% print.data.frame() +tab2 <- mark_poly(tab2) +head(tab2) %>% print.data.frame() +tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1) +head(tab2) %>% print.data.frame() +mark_prop(tab2,cross="F2",pval=0.05) %>% head() %>% print.data.frame() +tab2 <- mark_allele(tab=tab2,ref=strains,par1="parent1",par2="parent2") +tab2 %>% arrange(desc(exclude_allele)) %>% head() %>% print.data.frame() +strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354", +"gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% +select(marker,parent1,parent2) %>% print.data.frame() +stuart_cross <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, +ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox") +stuart_cross[1:10,1:7] %>% print.data.frame() +library(qtl) +data(stuart_cross) +summary(stuart_cross) +library(qtl) +data(stuart_cross) +summary(stuart_cross) +save(stuart_cross,file="stuart_cross.rda") +library(qtl) +data(stuart_cross) +summary(stuart_cross) +stuart_cross <- calc.genoprob(stuart_cross, step=2.0, off.end=0.0, +error.prob=1.0E-4, map.function="haldane", stepwidth="fixed") +save(stuart_cross,file="stuart_cross.rda") +devtools::build(path=".",vignettes=FALSE) +roxygen2::roxygenise() +devtools::build_manual(path=".") +devtools::build_vignettes() +load("/mnt/zeus/zeus_mouselab/marie/stuart_estmap/newmap.rda") +usethis::use_data(newmap,stuart_newmap) +stuart_newmap <- newmap +usethis::use_data(stuart_newmap,stuart_newmap) +roxygen2::roxygenise() +devools::build(path=".",vignettes=FALSE) +devtools::build(path=".",vignettes=FALSE) +library(stuart) +devtools::build_manual() +devtools::build_manual(path=".") +devtools::build_vignettes() +library(dplyr) +library(stuart) +data(stuart_estmap) +data(stuart_newmap) +knitr::opts_chunk$set( +collapse = TRUE, +comment = "#>" +) +library(dplyr) +library(stuart) +annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) +data(genos) +summary(genos) +data(phenos) +summary(phenos) +strains <- geno_strains(ref=annot_mini,geno=genos,par1=c("StrainsA_1","StrainsA_2"), +par2=c("StrainsB_1","StrainsB_2"),name1="parent1",name2="parent2") +head(strains) %>% print.data.frame() +genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", +"StrainsB_1","StrainsB_2")) +# how to use the function: +# stuart_tab <- tab_mark(genos) +data(stuart_tab) +summary(stuart_tab) +tab2 <- mark_match(stuart_tab,ref=strains) +tab2 %>% filter(exclude_match==1) %>% print.data.frame() +tab2 <- mark_poly(tab2) +head(tab2) %>% print.data.frame() +tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1) +head(tab2) %>% print.data.frame() +mark_prop(tab2,cross="F2",pval=0.05) %>% head() %>% print.data.frame() +tab2 <- mark_allele(tab=tab2,ref=strains,cross="F2",par1="parent1",par2="parent2") +tab2 %>% arrange(desc(exclude_allele)) %>% head() %>% print.data.frame() +strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354", +"gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% +select(marker,parent1,parent2) %>% print.data.frame() +stuart_cross <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, +ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox") +stuart_cross[1:10,1:7] %>% print.data.frame() +library(qtl) +data(stuart_cross) +summary(stuart_cross) +data(stuart_newmap) +plotMap(stuart_cross,stuart_newmap,shift=TRUE) +bad_markers <- mark_estmap(map=stuart_newmap,dist=100) +print(bad_markers) +tab2 %>% filter(marker %in% bad_markers) +strains %>% filter(marker %in% bad_markers) +newcross <- drop.markers(cross=stuart_cross,markers=bad_markers) +tab2 <- tab2 %>% filter(!marker %in% bad_markers) +good_rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, +ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path = "rqtl_file.csv") +#goodnewmap <- est.map(newcross) +save(good_rqtl_file,file="good_rqtl_file.rda") +load("/mnt/zeus/zeus_mouselab/marie/stuart_goodestmap/good_rqtl_file.rda") +write.table(good_rqtl_file,file="goodqtlfile.csv") +knitr::opts_chunk$set( +collapse = TRUE, +comment = "#>" +) +library(dplyr) +library(stuart) +annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) +data(genos) +summary(genos) +data(phenos) +summary(phenos) +strains <- geno_strains(ref=annot_mini,geno=genos,par1=c("StrainsA_1","StrainsA_2"), +par2=c("StrainsB_1","StrainsB_2"),name1="parent1",name2="parent2") +head(strains) %>% print.data.frame() +genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", +"StrainsB_1","StrainsB_2")) +# how to use the function: +# stuart_tab <- tab_mark(genos) +data(stuart_tab) +summary(stuart_tab) +tab2 <- mark_match(stuart_tab,ref=strains) +tab2 %>% filter(exclude_match==1) %>% print.data.frame() +tab2 <- mark_poly(tab2) +head(tab2) %>% print.data.frame() +tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1) +head(tab2) %>% print.data.frame() +mark_prop(tab2,cross="F2",pval=0.05) %>% head() %>% print.data.frame() +tab2 <- mark_allele(tab=tab2,ref=strains,cross="F2",par1="parent1",par2="parent2") +tab2 %>% arrange(desc(exclude_allele)) %>% head() %>% print.data.frame() +strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354", +"gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% +select(marker,parent1,parent2) %>% print.data.frame() +stuart_cross <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, +ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox") +stuart_cross[1:10,1:7] %>% print.data.frame() +library(qtl) +data(stuart_cross) +summary(stuart_cross) +data(stuart_newmap) +plotMap(stuart_cross,stuart_newmap,shift=TRUE) +bad_markers <- mark_estmap(map=stuart_newmap,dist=100) +print(bad_markers) +tab2 %>% filter(marker %in% bad_markers) +strains %>% filter(marker %in% bad_markers) +newcross <- drop.markers(cross=stuart_cross,markers=bad_markers) +save(newcross,file="newcross.rda") diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths index 7e87f342bc406ce58dc45ffce61dc48bb30c04eb..2532f969ec5b78787b271919d9f85c4e82d00a58 100644 --- a/.Rproj.user/shared/notebooks/paths +++ b/.Rproj.user/shared/notebooks/paths @@ -12,6 +12,7 @@ /home/marie/Documents/stuart_package/stuart/R/mark_prop.R="19C6446D" /home/marie/Documents/stuart_package/stuart/R/phenos-data.R="75BCF577" /home/marie/Documents/stuart_package/stuart/R/stuart_cross-data.R="A0BD213E" +/home/marie/Documents/stuart_package/stuart/R/stuart_newmap-data.R="35360211" /home/marie/Documents/stuart_package/stuart/R/stuart_tab-data.R="9C18AF59" /home/marie/Documents/stuart_package/stuart/R/tab_mark.R="60E85DC0" /home/marie/Documents/stuart_package/stuart/R/write_rqtl.R="D25FAC55" diff --git a/R/stuart_newmap-data.R b/R/stuart_newmap-data.R new file mode 100755 index 0000000000000000000000000000000000000000..a8b0039e38fd54ead69f2da71c3341ec4981eebf --- /dev/null +++ b/R/stuart_newmap-data.R @@ -0,0 +1,8 @@ +#' Newmap for stuart data +#' +#' Result of qtl::est.map() function with stuart data after exclusion of problematic markers +#' +#' @format A list of 20 elements + + +"stuart_newmap" diff --git a/data/stuart_newmap.rda b/data/stuart_newmap.rda new file mode 100644 index 0000000000000000000000000000000000000000..a94010c626c64ecbba11c564881edab5b3dc8ce1 Binary files /dev/null and b/data/stuart_newmap.rda differ diff --git a/man/stuart_newmap.Rd b/man/stuart_newmap.Rd new file mode 100644 index 0000000000000000000000000000000000000000..4fe7820ef888a853dd24431541ff23a81747014f --- /dev/null +++ b/man/stuart_newmap.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stuart_newmap-data.R +\docType{data} +\name{stuart_newmap} +\alias{stuart_newmap} +\title{Newmap for stuart data} +\format{ +A list of 20 elements +} +\usage{ +stuart_newmap +} +\description{ +Result of qtl::est.map() function with stuart data after exclusion of problematic markers +} +\keyword{datasets} diff --git a/stuart_1.0.3.9000.pdf b/stuart_1.0.3.9000.pdf index 5959a609b5698a37725340268b59e96b37c87d35..671f27e63e1fd385a42f93255b2f2664a852b1a3 100644 Binary files a/stuart_1.0.3.9000.pdf and b/stuart_1.0.3.9000.pdf differ diff --git a/stuart_1.0.3.9000.tar.gz b/stuart_1.0.3.9000.tar.gz index 56ccea5f23abd2d56a2d0b388d52050d090b4c1b..bbcc31e8824882d227ad9f7c0ba262b10a9d84cc 100644 Binary files a/stuart_1.0.3.9000.tar.gz and b/stuart_1.0.3.9000.tar.gz differ diff --git a/vignettes/stuart.Rmd b/vignettes/stuart.Rmd index 127a4ea3d8707df69921e7fdca3a72092ca07fd4..463314b3f0e440c94d351fc9f55267eaad432408 100755 --- a/vignettes/stuart.Rmd +++ b/vignettes/stuart.Rmd @@ -174,30 +174,22 @@ library(qtl) data(stuart_cross) summary(stuart_cross) -stuart_cross <- calc.genoprob(stuart_cross, step=2.0, off.end=0.0, - error.prob=1.0E-4, map.function="haldane", stepwidth="fixed") - -##### TO RUN ##### -# stuart_newmap <- est.map(cross=stuart_cross,error.prob=0.01) -# save the map as data in the package - -# data(stuart_newmap) -# summary(stuart_newmap) +data(stuart_newmap) ``` Then we can plot the newmap to compare it with the positions found in the annotation file. ```{r plot_map} -# plotMap(stuart_cross,stuart_newmap,shift=TRUE) +plotMap(stuart_cross,stuart_newmap,shift=TRUE) ``` ### mark_estmap It seems to be still 2 problematic markers. We can identify them with the `mark_estmap()`function. ```{r mark_estmap} -# bad_markers <- mark_estmap(map=stuart_newmap,dist=100) -# print(bad_markers) +bad_markers <- mark_estmap(map=stuart_newmap,dist=100) +print(bad_markers) ``` ### Exclusion of last problematic markers @@ -205,15 +197,13 @@ It seems to be still 2 problematic markers. We can identify them with the `mark_ Let's focus on these markers. ```{r bad_markers_tab} -tab2 %>% filter(marker %in% c("S6J010381992","S6J205609960")) -# tab2 %>% filter(marker %in% bad_markers) +tab2 %>% filter(marker %in% bad_markers) ``` -We can see that these two markers have correct proportions of each genotypes so there not excluded. Let's see the alleles of the parental strains for these markers. +We can see that these two markers have correct proportions of each genotypes so they were not excluded. Let's see the alleles of the parental strains for these markers. ```{r bad_markers_annot} -strains %>% filter(marker %in% c("S6J010381992","S6J205609960")) -# strains %>% filter(marker %in% bad_markers) +strains %>% filter(marker %in% bad_markers) ``` We can see that for these markers, one of the two parents has a missing genotype. In the cases where one parent has a missing of a heterozygous genotype, the `mark_allele()` function does not by default exclude the marker if the observed allele of the genotyped parent correspond to one of the alleles observed in the crossed individuals. @@ -223,16 +213,15 @@ To avoid the issues with these situations, you can use the `parNH=FALSE` argumen We rather advise you to keep these markers and identify possible problematic markers with `mark_estmap()` and exlude them. You have two possibilities to do so. First, you can drop these markers in the cross object with `qtl::drop.markers()`: ```{r drop_markers} -# newcross <- drop.markers(cross=stuart_cross,markers=bad_markers) +newcross <- drop.markers(cross=stuart_cross,markers=bad_markers) ``` Or if you prefere to have a correct CSV object that can always be used for QTL analysis, you can delete these markers in the marker tab and create a new cross file with `write_rqtl()`: ```{r} -# tab2 <- tab2 %>% filter(!marker %in% bad_markers) -# -# good_rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, -# ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path = "rqtl_file.csv") -``` +tab2 <- tab2 %>% filter(!marker %in% bad_markers) +good_rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2, + ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox") +```