Commit 2054450f authored by mariefbourdon's avatar mariefbourdon
Browse files

modif mark_prop chi

parent 716a1393
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
}
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
#recode genotypes depending on parents' genotypes
geno <- geno %>% mutate(Geno = case_when(
#if one allele not genotyped:
Allele1...Forward=="N" | Allele2...Forward=="N" ~ "NA",
#if both alleles genotyped
##homozygous 0
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par1) ~ "0",
##homozygous 2
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par2) ~ "2",
##heterozygous
Allele1...Forward!=Allele2...Forward ~ "1",
#if parental strains are N/H
##homozygous for parent that is N/H
###homozygous 0
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par1)%in%c("H","N") ~ "0",
###homozygous 2
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par2)%in%c("H","N") ~ "2",
)
)
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
#recode genotypes depending on parents' genotypes
geno <- geno %>% mutate(Geno = case_when(
#if one allele not genotyped:
Allele1...Forward=="N" | Allele2...Forward=="N" ~ "NA",
#if both alleles genotyped
##homozygous 0
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par1) ~ "0",
##homozygous 2
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par2) ~ "2",
##heterozygous
Allele1...Forward!=Allele2...Forward ~ "1",
#if parental strains are N/H
##homozygous for parent that is N/H
###homozygous 0
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par1)%in%c("H","N") ~ "0",
###homozygous 2
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par2)%in%c("H","N") ~ "2",
)
)
#keep positions of markers
markers <- geno %>% select(SNP.Name,chr,!!sym(pos)) %>% distinct()
markers <- markers %>% arrange(chr,bp_mm10)
#keep only interesting columns in geno file
geno <- geno %>% arrange(chr,bp_mm10)
geno <- geno %>% select(SNP.Name,Sample.ID,Geno)
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
#recode genotypes depending on parents' genotypes
geno <- geno %>% mutate(Geno = case_when(
#if one allele not genotyped:
Allele1...Forward=="N" | Allele2...Forward=="N" ~ "NA",
#if both alleles genotyped
##homozygous 0
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par1) ~ "0",
##homozygous 2
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par2) ~ "2",
##heterozygous
Allele1...Forward!=Allele2...Forward ~ "1",
#if parental strains are N/H
##homozygous for parent that is N/H
###homozygous 0
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par1)%in%c("H","N") ~ "0",
###homozygous 2
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par2)%in%c("H","N") ~ "2",
)
)
#keep positions of markers
markers <- geno %>% select(SNP.Name,chr,!!sym(pos)) %>% distinct()
markers <- markers %>% arrange(chr,bp_mm10)
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
#recode genotypes depending on parents' genotypes
geno <- geno %>% mutate(Geno = case_when(
#if one allele not genotyped:
Allele1...Forward=="N" | Allele2...Forward=="N" ~ "NA",
#if both alleles genotyped
##homozygous 0
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par1) ~ "0",
##homozygous 2
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par2) ~ "2",
##heterozygous
Allele1...Forward!=Allele2...Forward ~ "1",
#if parental strains are N/H
##homozygous for parent that is N/H
###homozygous 0
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par1)%in%c("H","N") ~ "0",
###homozygous 2
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par2)%in%c("H","N") ~ "2",
)
)
#keep positions of markers
markers <- geno %>% select(SNP.Name,chr,!!sym(pos)) %>% distinct()
markers <- markers %>% arrange(chr,bp_mm10)
print(markers)
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
#keep positions of markers
markers <- geno %>% select(SNP.Name,chr,bp_mm10,!!sym(pos)) %>% distinct()
write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#extract snps non excluded
if("exclude_match" %in% colnames(tab)){
tab <- tab %>% filter(exclude_match==0)
}
if("exclude_poly" %in% colnames(tab)){
tab <- tab %>% filter(exclude_poly==0)
}
if("exclude_prop" %in% colnames(tab)){
tab <- tab %>% filter(exclude_prop==0)
}
if("exclude_allele" %in% colnames(tab)){
tab <- tab %>% filter(exclude_allele==0)
}
#filter genotypes for non excluded markers in geno file
geno <- geno %>% select(c(SNP.Name,Sample.ID,Allele1...Forward,Allele2...Forward)) %>% filter(SNP.Name %in% tab$SNP.Name)
#recode parents' names to match column names nomenclature
par1 <- make.names(par1)
par2 <- make.names(par2)
#keep parental lines genotypes
colnames(ref) <- make.names(colnames(ref))
ref <- ref %>% select(marker,chr,bp_mm10,!!sym(pos),!!sym(par1),!!sym(par2))
#merge genotypes with parents
geno <- left_join(geno,ref,by=c("SNP.Name"="marker"))
#recode "-" in "N" in geno file
geno <- geno %>% mutate(Allele1...Forward = recode(Allele1...Forward,
"-" = "N"))
geno <- geno %>% mutate(Allele2...Forward = recode(Allele2...Forward,
"-" = "N"))
#recode geno in factors with same levels
geno <- geno %>% mutate(Allele1...Forward = factor(Allele1...Forward,levels=c("A","C","G","H","N","T")))
geno <- geno %>% mutate(Allele2...Forward = factor(Allele2...Forward,levels=c("A","C","G","H","N","T")))
#recode genotypes depending on parents' genotypes
geno <- geno %>% mutate(Geno = case_when(
#if one allele not genotyped:
Allele1...Forward=="N" | Allele2...Forward=="N" ~ "NA",
#if both alleles genotyped
##homozygous 0
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par1) ~ "0",
##homozygous 2
Allele1...Forward==Allele2...Forward & Allele1...Forward==!!sym(par2) ~ "2",
##heterozygous
Allele1...Forward!=Allele2...Forward ~ "1",
#if parental strains are N/H
##homozygous for parent that is N/H
###homozygous 0
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par1)%in%c("H","N") ~ "0",
###homozygous 2
(!!sym(par1)%in%c("H","N") | !!sym(par2)%in%c("H","N")) &
!!sym(par2)%in%c("H","N") ~ "2",
)
)
#keep positions of markers
markers <- geno %>% select(SNP.Name,chr,bp_mm10,!!sym(pos)) %>% distinct()
print(markers)
markers <- markers %>% arrange(chr,bp_mm10)
}
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
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)
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", "StrainsB_1","StrainsB_2"))
data(stuart_tab)
summary(stuart_tab)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 %>% filter(exclude_match==1)
tab2 <- mark_poly(tab2)
head(tab2)
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1)
head(tab2)
tab2 <- mark_allele(tab=tab2,ref=strains,par1="parent1",par2="parent2")
tab2 %>% arrange(desc(exclude_allele)) %>% head()
strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354","gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% select(marker,parent1,parent2)
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
devtools::build(path="/Users/mariebourdon/stuart_package/stuart",vignettes = FALSE)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
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)
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", "StrainsB_1","StrainsB_2"))
data(stuart_tab)
summary(stuart_tab)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 %>% filter(exclude_match==1)
tab2 <- mark_poly(tab2)
head(tab2)
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1)
head(tab2)
tab2 <- mark_allele(tab=tab2,ref=strains,par1="parent1",par2="parent2")
tab2 %>% arrange(desc(exclude_allele)) %>% head()
strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354","gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% select(marker,parent1,parent2)
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
library(dplyr)
library(stuart)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(dplyr)
library(stuart)
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)
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", "StrainsB_1","StrainsB_2"))
data(stuart_tab)
summary(stuart_tab)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 %>% filter(exclude_match==1)
tab2 <- mark_poly(tab2)
head(tab2)
tab <- tab %>% mutate(n_geno = tab %>% select(n_HM1,n_HM2,n_HT) %>% rowSums(na.rm=TRUE))
#calculate proportion of each genotype
tab <- tab %>% mutate(p_HM1 = n_HM1/n_geno)
tab <- tab %>% mutate(p_HM2 = n_HM2/n_geno)
tab <- tab %>% mutate(p_HT = n_HT/n_geno)
tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))
#stop if cross != "F2" or "N2"
if(!cross %in% c("F2","N2")){
stop("Cross must be F2 or N2")
}
#stop of homo&hetero or pval not specified
if((is.na(homo)==TRUE & is.na(hetero)==TRUE) | is.na(pval)==TRUE){
stop("Arguments homo and hetero or argument pval must be specified")
}
# #stop if no value for "homo" for F2 cross
# if(cross=="F2" & (is.na(homo)==TRUE | is.na(hetero)==TRUE)){
# stop("Arguments homo and hetero must be specified for F2 crosses")
# }
#
# #stop if no value for "homo" and "hetero" for N2 cross
# if(cross=="N2" & (is.na(homo)==TRUE | is.na(hetero)==TRUE)){
# stop("Arguments homo and hetero must be specified for N2 crosses")
# }
#exclude markers according to proportion of na
tab <- tab %>% mutate(exclude_prop=case_when(p_NA > na ~ 1, #exclude markers according to proportion of na
cross=="F2" & (p_HM1 < homo | p_HM2 < homo | p_HT < hetero) ~ 1, #exclude markers according to proportion of homozygous individuals for F2 cross
cross=="N2" & (p_HM1 < homo | p_HT < hetero) ~ 1, #exclude markers according to proportion of homozygous and heterozygous individuals for N2 cross
T ~ 0))
tab <- tab %>% select(-c(p_HM1,p_HM2,p_HT,p_NA,n_geno))
return(tab)
}
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1)
head(tab2)
tab2 <- mark_allele(tab=tab2,ref=strains,par1="parent1",par2="parent2")
tab2 %>% arrange(desc(exclude_allele)) %>% head()
strains %>% filter(marker %in% c("gJAX00038569","gJAX00425031","gUNC12245354","gUNC15530876","gUNC21555204","gUNC21596600")) %>% arrange(marker) %>% select(marker,parent1,parent2)
rqtl_file <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")
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)
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", "StrainsB_1","StrainsB_2"))
data(stuart_tab)
summary(stuart_tab)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 %>% filter(exclude_match==1)
tab2 <- mark_poly(tab2)
head(tab2)
mark_prop <- function(tab,cross,homo=NA,hetero=NA,pval=NA,na=0.5){
#calculate total number of individuals genotyped for each marker
tab <- tab %>% mutate(n_geno = tab %>% select(n_HM1,n_HM2,n_HT) %>% rowSums(na.rm=TRUE))
#calculate proportion of each genotype
tab <- tab %>% mutate(p_HM1 = n_HM1/n_geno)
tab <- tab %>% mutate(p_HM2 = n_HM2/n_geno)
tab <- tab %>% mutate(p_HT = n_HT/n_geno)
tab <- tab %>% mutate(p_NA = n_NA/(n_geno+n_NA))
#stop if cross != "F2" or "N2"
if(!cross %in% c("F2","N2")){
stop("Cross must be F2 or N2")
}
#stop of homo&hetero or pval not specified
if((is.na(homo)==TRUE | is.na(hetero)==TRUE) & is.na(pval)==TRUE){
stop("Arguments homo and hetero or argument pval must be specified")
}
# #stop if no value for "homo" for F2 cross
# if(cross=="F2" & (is.na(homo)==TRUE | is.na(hetero)==TRUE)){
# stop("Arguments homo and hetero must be specified for F2 crosses")
# }
#
# #stop if no value for "homo" and "hetero" for N2 cross
# if(cross=="N2" & (is.na(homo)==TRUE | is.na(hetero)==TRUE)){
# stop("Arguments homo and hetero must be specified for N2 crosses")
# }
#exclude markers according to proportion of na
tab <- tab %>% mutate(exclude_prop=case_when(p_NA > na ~ 1, #exclude markers according to proportion of na
cross=="F2" & (p_HM1 < homo | p_HM2 < homo | p_HT < hetero) ~ 1, #exclude markers according to proportion of homozygous individuals for F2 cross
cross=="N2" & (p_HM1 < homo | p_HT < hetero) ~ 1, #exclude markers according to proportion of homozygous and heterozygous individuals for N2 cross
T ~ 0))
tab <- tab %>% select(-c(p_HM1,p_HM2,p_HT,p_NA,n_geno))
return(tab)