Commit 1367ad7d authored by Marie Bourdon's avatar Marie Bourdon
Browse files

modify vignette

parent 03f045b2
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")) annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
data(genos) data(genos)
summary(genos) summary(genos)
...@@ -510,3 +345,168 @@ genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", ...@@ -510,3 +345,168 @@ genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2",
data(stuart_tab) data(stuart_tab)
summary(stuart_tab) summary(stuart_tab)
View(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")
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
/home/marie/Documents/stuart_package/stuart/R/mark_prop.R="19C6446D" /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/phenos-data.R="75BCF577"
/home/marie/Documents/stuart_package/stuart/R/stuart_cross-data.R="A0BD213E" /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/stuart_tab-data.R="9C18AF59"
/home/marie/Documents/stuart_package/stuart/R/tab_mark.R="60E85DC0" /home/marie/Documents/stuart_package/stuart/R/tab_mark.R="60E85DC0"
/home/marie/Documents/stuart_package/stuart/R/write_rqtl.R="D25FAC55" /home/marie/Documents/stuart_package/stuart/R/write_rqtl.R="D25FAC55"
......
#' 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"
% 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}
No preview for this file type
No preview for this file type
...@@ -174,30 +174,22 @@ library(qtl) ...@@ -174,30 +174,22 @@ library(qtl)
data(stuart_cross) data(stuart_cross)
summary(stuart_cross) summary(stuart_cross)
stuart_cross <- calc.genoprob(stuart_cross, step=2.0, off.end=0.0, data(stuart_newmap)
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)
``` ```
Then we can plot the newmap to compare it with the positions found in the annotation file. Then we can plot the newmap to compare it with the positions found in the annotation file.
```{r plot_map} ```{r plot_map}
# plotMap(stuart_cross,stuart_newmap,shift=TRUE) plotMap(stuart_cross,stuart_newmap,shift=TRUE)
``` ```
### mark_estmap ### mark_estmap
It seems to be still 2 problematic markers. We can identify them with the `mark_estmap()`function. It seems to be still 2 problematic markers. We can identify them with the `mark_estmap()`function.
```{r mark_estmap} ```{r mark_estmap}
# bad_markers <- mark_estmap(map=stuart_newmap,dist=100) bad_markers <- mark_estmap(map=stuart_newmap,dist=100)
# print(bad_markers) print(bad_markers)
``` ```
### Exclusion of last problematic 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_ ...@@ -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. Let's focus on these markers.
```{r bad_markers_tab} ```{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} ```{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. 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 ...@@ -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()`: 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} ```{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()`: 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} ```{r}
# tab2 <- tab2 %>% filter(!marker %in% 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")
```
good_rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,
ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
```
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