diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths
index a68c4f6117f52e4b4130c6af6e151b7733125f14..5bb19586df1d3deb5ca60eb400ec87138317759e 100644
--- a/.Rproj.user/shared/notebooks/paths
+++ b/.Rproj.user/shared/notebooks/paths
@@ -1,7 +1,11 @@
-/home/marie/Documents/HB_CC011/F2/Geno/Maestro/F2_geno_estmap1.R="42B0AE06"
/home/marie/Documents/stuart/stuart_package/stuart/DESCRIPTION="55556B0B"
+/home/marie/Documents/stuart/stuart_package/stuart/NAMESPACE="D52CF639"
/home/marie/Documents/stuart/stuart_package/stuart/R/geno_strains.R="68944065"
/home/marie/Documents/stuart/stuart_package/stuart/R/mark_estmap.R="F759F000"
+/home/marie/Documents/stuart/stuart_package/stuart/R/mark_poly.R="25EFE570"
+/home/marie/Documents/stuart/stuart_package/stuart/R/mark_prop.R="E4688428"
+/home/marie/Documents/stuart/stuart_package/stuart/R/tab_mark.R="37E604E4"
+/home/marie/Documents/stuart/stuart_package/stuart/R/write_rqtl.R="04457D88"
/home/marie/Documents/stuart/stuart_package/stuart/README.Rmd="4E06E7DC"
/home/marie/Documents/stuart/stuart_package/stuart/doc/stuart.R="A03ACF3D"
/home/marie/Documents/stuart/stuart_package/stuart/man/genos.Rd="75E7412F"
diff --git a/DESCRIPTION b/DESCRIPTION
index 88b43cabf627310d01654cd72213eee16cc610d8..f44252071ae2376c75d819588d19f3a1a1ac8b51 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: stuart
Title: Sort markers of lab strains genotyping results
-Version: 1.0.7
+Version: 1.1.0
Authors@R:
person(given = "Marie",
family = "Bourdon",
@@ -12,7 +12,7 @@ License: GPL-3
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.1
Depends:
R (>= 3.5.0),
dplyr,
diff --git a/R/mark_allele.R b/R/mark_allele.R
index 3abbcee969bff9ba6492cebf333569407f9b70f9..9a350010a4bf39324e3e920edeab67046fd4457b 100755
--- a/R/mark_allele.R
+++ b/R/mark_allele.R
@@ -6,6 +6,7 @@
#' @param cross cross type, either "F2" or "N2".
#' @param par1 first parental strain used in the cross, the name must be written as in the "ref" data frame. For a backcross the strain used to backcross F1 individuals must be indicated as par1.
#' @param par2 second parental strain used in the cross, the name must be written as in the "ref" data frame.
+#' @param parNH wether markers for which one of the parental strains has an undetermined or heterozygous genotype (default is TRUE)
#'
#' @import dplyr
#'
diff --git a/R/mark_prop.R b/R/mark_prop.R
index 0d756dfddafaae10d75b2cb2f11f383c04524015..77794e868db0bdb468495665d4ce08b4ec3b5535 100755
--- a/R/mark_prop.R
+++ b/R/mark_prop.R
@@ -16,7 +16,7 @@
#' @param homo1X a vector of two numbers. The lower and upper limits for the proportion of homozygous individuals for markers on X chromosome. This argument is for homozygous genotype with the highest expected proportion.
#' @param homo2X a vector of two numbers. The lower and upper limits for the proportion of homozygous individuals for markers on X chromosome. This argument is for homozygous genotype with the lowest expected proportion.
#' @param heteroX a vector of two numbers. The lower and upper limits for the proportion of heterozygous individuals for markers on X chromosome.
-#' @param na proportion of non-genotyped individuals above which the marker is excluded.
+#' @param pval chi2 p value threshold under which the marker is excluded. Must not be used if the homo and hetero arguments are used.
#'
#' @import dplyr
#' @import tidyselect
diff --git a/R/tab_mark.R b/R/tab_mark.R
index 3cc24c9c4c2011a9b6bd5029cb85fa10b0407560..b8a85dd60a0c31f32911f5c31b33b95d615055ed 100755
--- a/R/tab_mark.R
+++ b/R/tab_mark.R
@@ -23,137 +23,73 @@ tab_mark <- function(geno,annot,pos){
stop("No position data")
}
- #rename df columns
- geno <- geno %>% rename("marker"=1,
- "id"=2,
- "allele_1"=3,
- "allele_2"=4)
-
- #create geno column in geno df
- geno <- geno %>% unite(Geno,c("allele_1","allele_2"),sep="",remove=FALSE)
-
- #recode genotypes to have all heterozygous encoded the same way (ex: only "AT", no "TA")
- geno <- geno %>% mutate(Geno=recode(Geno,
- "TA" = "AT",
- "GA" = "AG",
- "CA" = "AC",
- "GT" = "TG",
- "CT" = "TC",
- "GC" = "CG"))
-
-
- #create df with counts for each genotype
- tab <- tibble(marker = as.character(unique(geno$marker)),
- allele_1 = NA,
- allele_2 = NA,
- n_HM1 = NA,
- n_HM2 = NA,
- n_HT = NA,
- n_NA = NA)
-
-
- ## loop to count genotype
- for(i in tab$marker){
- #extract alleles for each marker
- Alleles <- geno %>% filter(marker==i) %>%
- select(c(marker,id,Geno,allele_1,allele_2)) %>%
- pivot_longer(c(allele_1,allele_2),names_to="Allele_name",values_to="Allele") %>%
- distinct(Allele) %>% filter(Allele != "-")
- Alleles <- as.factor(paste(Alleles$Allele))
-
- #sort alleles
- Alleles <- factor(Alleles,levels=c("A","T","C","G"))
- Alleles <- sort(Alleles)
-
- #add alleles and counts, only for markers with alleles (not markers with no genotyped ind)
- if(all(rapportools::is.empty(Alleles))==FALSE){
-
- #add alleles to tab
- tab <- tab %>% mutate(allele_1 = ifelse(marker == i,
- paste(Alleles[1]), allele_1))
-
-
-
- #count for homozygous for allele 1
- n1 <- geno %>% filter(marker==i) %>%
- filter(Geno == paste(Alleles[1],Alleles[1],sep="")) %>%
- summarise(n=n())
-
-
- #add count for homozygous for allele 1 to tab
- tab <- tab %>% mutate(n_HM1 = ifelse(marker == i,
- n1$n, n_HM1))
-
-
- }
-
- #if marker not polymorphic
- if(is.na(Alleles[2])==TRUE){
- #NA as allele_2
- tab <- tab %>% mutate(allele_2 = ifelse(marker == i,
- NA, allele_2))
-
- #NA as n_HM2
- tab <- tab %>% mutate(n_HM2 = ifelse(marker == i,
- NA, n_HM2))
-
- #NA as n_HT
- tab <- tab %>% mutate(n_HT = ifelse(marker == i,
- NA, n_HT))
- } else {
- #add alleles to tab
- tab <- tab %>% mutate(allele_2 = ifelse(marker == i,
- paste(Alleles[2]), allele_2))
-
-
- #count for homozygous for allele 2
- n2 <- geno %>% filter(marker==i) %>%
- filter(Geno == paste(Alleles[2],Alleles[2],sep="")) %>%
- summarise(n=n())
-
- #add count for homozygous for allele 1 to tab
- tab <- tab %>% mutate(n_HM2 = ifelse(marker == i,
- n2$n, n_HM2))
-
-
- #count for heterozygous
- n3 <- geno %>% filter(marker==i) %>%
- filter(Geno == paste(Alleles[1],Alleles[2],sep="")) %>%
- summarise(n=n())
-
-
- #add count for homozygous for allele 1 to tab
- tab <- tab %>% mutate(n_HT = ifelse(marker == i,
- n3$n, n_HT))
-
-
- }
-
- #count for NA
- n4 <- geno %>% filter(marker==i) %>%
- filter(Geno == "--" |
- Geno == paste(Alleles[1],"-",sep="") | Geno == paste(Alleles[2],"-",sep="") |
- Geno == paste("-",Alleles[1],sep="") | Geno == paste("-",Alleles[2],sep="")) %>%
- summarise(n=n())
-
- #add count for NA to tab
- tab <- tab %>% mutate(n_NA = ifelse(marker == i,
- n4$n, n_NA))
+ geno <- geno %>% rename(marker = 1, id = 2, allele_1 = 3,
+ allele_2 = 4)
+
+ tab <- geno %>% group_by(marker) %>%
+ summarise(alleles=toString(unique(c(allele_1,allele_2))))
+
+ tab <- tab %>% separate(alleles,c("al1","al2","al3"),sep=", ",fill="right") %>%
+ mutate(allele_1=case_when(al1 == "A" | al2 == "A" | al3 == "A" ~ "A",
+ al1 == "T" | al2 == "T" | al3 == "T" ~ "T",
+ al1 == "C" | al2 == "C" | al3 == "C" ~ "C",
+ al1 == "G" | al2 == "G" | al3 == "G" ~ "G")) %>%
+ mutate(allele_2=case_when(#is.na(al2)==TRUE & is.na(al3)==TRUE ~ NA,
+ allele_1 != "T" & (al1=="T" | al2=="T" | al3=="T") ~ "T",
+ allele_1 != "C" & (al1=="C" | al2=="C" | al3=="C") ~ "C",
+ allele_1 != "G" & (al1=="G" | al2=="G" | al3=="G") ~ "G")) %>%
+ select(-c(al1,al2,al3))
+
+ # genotype per ind
+ geno <- geno %>% unite(Geno, c("allele_1", "allele_2"), sep = "",
+ remove = FALSE) %>% select(marker,id,Geno)
+ geno <- geno %>% mutate(Geno = recode(Geno, TA = "AT", GA = "AG",
+ CA = "AC", GT = "TG", CT = "TC", GC = "CG"))
+
+ #summarise
+ geno <- geno %>% group_by(Geno,marker) %>% summarise(n=n())
+ geno <- geno %>% pivot_wider(names_from=Geno,values_from=n)
+
+ # add columns if missing
+ genocols <- c(AA = NA_integer_, AT = NA_integer_, AG = NA_integer_, AC = NA_integer_,
+ TT = NA_integer_, TC = NA_integer_, TG = NA_integer_,
+ CC = NA_integer_, CG = NA_integer_,
+ GG = NA_integer_)
+
+ if(!rapportools::is.empty(setdiff(names(genocols), names(geno)))){
+ geno <- geno %>% add_column(!!!genocols[setdiff(names(genocols), names(geno))])
}
- #change class of counts as numeric :
- tab$n_HM1 <- tab$n_HM1 %>% as.numeric()
- tab$n_HM2 <- tab$n_HM2 %>% as.numeric()
- tab$n_HT <- tab$n_HT %>% as.numeric()
- tab$n_NA <- tab$n_NA %>% as.numeric()
- #add 0 for null counts
- tab <- tab %>% mutate_at(.vars=vars(n_HM1,n_HM2,n_HT,n_NA),~replace(., is.na(.), 0))
-
- #save useful columns in annot dataframe
+ #join tab and geno
+ tab <- full_join(tab,geno,by="marker")
+
+ #replace NA with 0
+ tab <- tab %>% mutate_at(c(4:14), ~replace(., is.na(.), 0))
+
+ #create n_HM1, n_HT, n_HM2, n_NA columns and suppress others
+ tab <- tab %>% mutate(n_HM1=case_when(allele_1=="A"~AA,
+ allele_1=="T"~TT,
+ allele_1=="C"~CC,
+ allele_1=="G"~GG,
+ T ~ 0)) %>%
+ mutate(n_HM2=case_when(allele_2=="T"~TT,
+ allele_2=="C"~CC,
+ allele_2=="G"~GG,
+ T ~ 0)) %>%
+ mutate(n_HT=case_when(allele_1=="A" & allele_2=="T" ~ AT,
+ allele_1=="A" & allele_2=="C" ~ AC,
+ allele_1=="A" & allele_2=="G" ~ AG,
+ allele_1=="T" & allele_2=="C" ~ TC,
+ allele_1=="T" & allele_2=="G" ~ TG,
+ allele_1=="C" & allele_2=="G" ~ CG,
+ T ~ 0)) %>%
+ rename("n_NA"="--") %>%
+ select(marker:allele_2,n_HM1:n_HT,n_NA)
+
+ #save useful columns in annot dataframe and merge
annot <- annot %>% select(marker,chr,!!sym(pos))
tab <- right_join(annot,tab,by="marker")
-
#return
return(tab)
}
diff --git a/R/write_rqtl.R b/R/write_rqtl.R
index 314e591bc2a1d8c9c041ddcddcbb847e7168c727..1c70ed21d4f1380ea449db7fe058ec3d4a5ec23b 100755
--- a/R/write_rqtl.R
+++ b/R/write_rqtl.R
@@ -20,7 +20,7 @@
#'
#### write_rqtl ####
## write data frame in rqtl format (csv), if path != NA writes the file in the path indicated
-write_rqtl <- function(geno,pheno,tab,ref,par1,par2,par_N=TRUE,prefix,pos,path=NA){
+write_rqtl <- function(geno,pheno,tab,ref,par1,par2,prefix,pos,path=NA){
#rename df columns
geno <- geno %>% rename("marker"=1,
"id"=2,
diff --git a/man/geno_strains.Rd b/man/geno_strains.Rd
index a8c4be973026aa65dbc9d3dd9ccb364b6d30a0d7..0faf914e6dd16afcb75d07fb6f70240144135304 100755
--- a/man/geno_strains.Rd
+++ b/man/geno_strains.Rd
@@ -2,7 +2,7 @@
% Please edit documentation in R/geno_strains.R
\name{geno_strains}
\alias{geno_strains}
-\title{Create haplotype for inbred strains into a dataframe}
+\title{Summarizes genotypes of multiple individuals}
\usage{
geno_strains(annot, geno, strn, cols)
}
@@ -16,7 +16,9 @@ geno_strains(annot, geno, strn, cols)
\item{cols}{name of the columns from the annot data frame to keep in the output of this function.}
}
\description{
-This functions adds columns for parental strains used in the cross in the annotation data frame, from the genotype data frame in which one or several animal of the parental strains were genotyped.
+This functions formats genotypes of genotyped individuals from a two letters encoding to a one letter encoding.
If several animals of one strain were genotyped, a consensus is created from these animals.
-The consensus is created as follow : if the indivuals carry the same allele, this allele is kept, otherwise, the allele is noted as "N". If individuals show residual heterozygosity, it is encoded as "H".
+The one letter encoding is made as follow: if the individual is homozygous, the letter of the allele (A, T, G or C) is kept. If the individual show residual heterozygosity, it is encoded as "H". If the genotype is missing, it is encoded as "N".
+The consensus is created as follow : if the indivuals carry the same allele, this allele is kept. If one or several individuals are heterozygous, the genotype is encoded as "H". If one individual has a missing genotypes but another was correctly genotyped, its genotype is kept.
+If the genotypes are completely discordant (i.e. if two individuals are both homozygous but for different alleles), the genotype is encoded as "D".
}
diff --git a/man/mark_allele.Rd b/man/mark_allele.Rd
index 8131b8618057be89f6d0e50c0f8b2dd184eebf45..8711aadcb9a89437ee9f393cc788041d08eb05a6 100755
--- a/man/mark_allele.Rd
+++ b/man/mark_allele.Rd
@@ -16,6 +16,8 @@ mark_allele(tab, ref, cross, par1, par2, parNH = TRUE)
\item{par1}{first parental strain used in the cross, the name must be written as in the "ref" data frame. For a backcross the strain used to backcross F1 individuals must be indicated as par1.}
\item{par2}{second parental strain used in the cross, the name must be written as in the "ref" data frame.}
+
+\item{parNH}{wether markers for which one of the parental strains has an undetermined or heterozygous genotype (default is TRUE)}
}
\description{
This functions uses the dataframe produced by the tab_mark function and fills the "exclude" column for all the markers which have alleles observed in the individuals of the cross that do not correspond to the alleles observed in the parental strains. For example, a marker which is not polymorphic between the two parental strains but which has two alleles in the cross individuals will be excluded.
diff --git a/man/mark_prop.Rd b/man/mark_prop.Rd
index b2ff00629645304859a9da155102f2e53f810b2a..beb797079e55f18ae118601c2a44fc30398c424f 100644
--- a/man/mark_prop.Rd
+++ b/man/mark_prop.Rd
@@ -12,8 +12,7 @@ mark_prop(
pval = NA,
homo1X = NULL,
homo2X = NULL,
- heteroX = NULL,
- na = 0.5
+ heteroX = NULL
)
}
\arguments{
@@ -25,13 +24,13 @@ mark_prop(
\item{hetero}{proportion of heterozygous individuals under which the marker is excluded.}
+\item{pval}{chi2 p value threshold under which the marker is excluded. Must not be used if the homo and hetero arguments are used.}
+
\item{homo1X}{a vector of two numbers. The lower and upper limits for the proportion of homozygous individuals for markers on X chromosome. This argument is for homozygous genotype with the highest expected proportion.}
\item{homo2X}{a vector of two numbers. The lower and upper limits for the proportion of homozygous individuals for markers on X chromosome. This argument is for homozygous genotype with the lowest expected proportion.}
\item{heteroX}{a vector of two numbers. The lower and upper limits for the proportion of heterozygous individuals for markers on X chromosome.}
-
-\item{na}{proportion of non-genotyped individuals above which the marker is excluded.}
}
\description{
This function uses the dataframe produced by the tab_mark function and fills the "exclude" column for all the markers that present too much missing genotypes or odd proportions of each genotype.
diff --git a/man/write_rqtl.Rd b/man/write_rqtl.Rd
index 0f1280d636698fcdec2c7c512a4fa921f764a668..0a52b05dbafb614aa89bd3209a8d3f86d255b50e 100755
--- a/man/write_rqtl.Rd
+++ b/man/write_rqtl.Rd
@@ -4,18 +4,7 @@
\alias{write_rqtl}
\title{Create data frame in Rqtl CSV format}
\usage{
-write_rqtl(
- geno,
- pheno,
- tab,
- ref,
- par1,
- par2,
- par_N = TRUE,
- prefix,
- pos,
- path = NA
-)
+write_rqtl(geno, pheno, tab, ref, par1, par2, prefix, pos, path = NA)
}
\arguments{
\item{geno}{data frame with the genotyping results for your cross}
diff --git a/stuart_1.0.7.pdf b/stuart_1.0.7.pdf
deleted file mode 100644
index 70522494431bc2abaf9cebe05fd9b604966ee7ef..0000000000000000000000000000000000000000
Binary files a/stuart_1.0.7.pdf and /dev/null differ
diff --git a/stuart_1.1.0.pdf b/stuart_1.1.0.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..1bf457d41368cb7fbce661c5b2d60c98fc887183
Binary files /dev/null and b/stuart_1.1.0.pdf differ