## ----message=FALSE, warning=FALSE, include=FALSE------------------------------
library(kableExtra)
library(dplyr)
library(nQuack)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# # Download all data files
# SetupBasicExample(overwrite = FALSE)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# library(nQuack)
# library(dplyr)
# library(kableExtra)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# # Set in and out paths of files
# inpath <- "../inst/extdata/01_raw/"
# outpath <- "../inst/extdata/02_prepared/"
# 
# # List files in the inpath and remove their ending
# filelist <- list.files(path = inpath, pattern = "*.bam" )
# filelist <- gsub(".bam", "", filelist)
# 
# for( i in 1:length(filelist)){
#   prepare_data(filelist[i], inpath, outpath)
# }
# 

## ----eval=FALSE, include=TRUE-------------------------------------------------
# inpathtext <- "../inst/extdata/02_prepared/"
# newfilelist <- list.files(path = inpathtext, pattern = "*.txt" )
# 
# for(i in 1:length(newfilelist)){
#     samp <- newfilelist[i]
#     temp <- process_data(paste0(inpathtext, samp),
#                          min.depth = 10,
#                          max.depth.quantile.prob = 1,
#                          error = 0.01,
#                          trunc = c(0.15,0.85))
# 
# 
#     write.csv(temp,
#               file = paste0("../inst/extdata/03_processed/", gsub(".txt", "", samp), ".csv"),
#               row.names = FALSE)
# }
# 

## ----eval=FALSE, include=TRUE-------------------------------------------------
# samples <- c("MLG013", "MLG014", "MLG015")
# 
# for(i in 1:length(samples)){
#   temp <- as.matrix(read.csv(paste0("../inst/extdata/03_processed/", samples[i], ".csv")))
#   out1 <- quackNormal(xm = temp,
#                       samplename = samples[i],
#                       cores = 10,
#                       parallel = TRUE)
#   out2 <- quackBeta(xm = temp,
#                     samplename = samples[i],
#                     cores = 10,
#                     parallel = TRUE)
#   out3 <- quackBetaBinom(xm = temp,
#                          samplename = samples[i],
#                          cores = 10,
#                          parallel = TRUE)
#   allout <- rbind(out1, out2, out3)
#   write.csv(allout,
#             file = paste0("../inst/extdata/04_output/", samples[i], ".csv"),
#             row.names = FALSE)
# }
# 

## ----echo=FALSE, message=FALSE, warning=FALSE, fig.align='center'-------------
recordtime <- read.csv(file = "../inst/extdata/04_output/modeltime.csv")
recordtime$normal <- round(abs(recordtime$normal), 2)
recordtime$beta <- round(abs(recordtime$beta), 2)
recordtime$beta.binomial <- round(abs(recordtime$beta.binomial), 2)



kableExtra::kbl(recordtime) %>%
  kable_paper("hover", full_width = F) %>%
   footnote("Time in seconds.")


## ----eval=FALSE, message=FALSE, warning=FALSE, include=TRUE-------------------
# inpathtext <- "../inst/extdata/04_output/"
# samples <- c("MLG013", "MLG014", "MLG015")
# 
# for(i in 1:length(samples)){
#   temp <- read.csv(paste0(inpathtext, samples[i], ".csv"))
#   summary <- quackit(model_out =  temp,
#                      summary_statistic = "BIC",
#                      mixtures = c("diploid", "triploid", "tetraploid"))
#   write.csv(summary,
#             file = paste0("../inst/extdata/05_interpret/", samples[i], ".csv"),
#             row.names = FALSE)
# }
# 

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.align='center'--------------
 # Create key
key <- data.frame(sample = c("MLG013", "MLG014", "MLG015"), 
           ploidal.level = c("diploid", "tetraploid", "triploid"))

# Read in quackit() output
dfs <- lapply(list.files("../inst/extdata/05_interpret/", full.names = TRUE  ), read.csv)
alloutput <- do.call(rbind, dfs)

# Combined
alloutputcombo <- dplyr::left_join(alloutput, key)

# Check the accuracy
alloutputcombo <- alloutputcombo %>%
                  dplyr::mutate(accuracy = ifelse(winnerBIC == ploidal.level, 1, 0))

## What distribution and model type should we use?
sumcheck <- alloutputcombo %>% 
            group_by(Distribution, Type) %>% 
            summarize(total = n(), correct = sum(accuracy))

kbl(sumcheck) %>%
  kable_paper("hover", full_width = F) 

## ----eval=FALSE, include=TRUE-------------------------------------------------
# samples <- c("MLG013", "MLG014", "MLG015")
# out <- c()
# 
# for(i in 1:length(samples)){
#   temp <- as.matrix(read.csv(paste0("../inst/extdata/03_processed/", samples[i], ".csv")))
#   out[[i]] <- bestquack(temp,
#                       distribution = "normal",
#                       type = "fixed",
#                       uniform = 1,
#                       mixtures = c("diploid", "triploid", "tetraploid"),
#                       samplename = samples[i])
# }
# 

## ----eval=FALSE, include=TRUE-------------------------------------------------
# samples <- c("MLG013", "MLG014", "MLG015")
# bout <- c()
# 
# for(i in 1:length(samples)){
#   temp <- as.matrix(read.csv(paste0("../inst/extdata/03_processed/", samples[i], ".csv")))
#   bout[[i]] <- quackNboots(temp,
#                         nboots = 100,
#                       distribution = "normal",
#                       type = "fixed",
#                       uniform = 1,
#                       mixtures = c("diploid", "triploid", "tetraploid"),
#                       samplename = samples[i])
# }
# write.csv(bout[[1]], file = "../inst/extdata/06_boots/MLG013-boots.csv", row.names = FALSE)
# write.csv(bout[[2]], file = "../inst/extdata/06_boots/MLG014-boots.csv", row.names = FALSE)
# write.csv(bout[[3]], file = "../inst/extdata/06_boots/MLG015-boots.csv", row.names = FALSE)

## -----------------------------------------------------------------------------
MLG013boot <- read.csv("../inst/extdata/06_boots/MLG013-boots.csv")
MLG013boot

## ----eval=FALSE, include=TRUE-------------------------------------------------
#   temp <- as.matrix(read.csv("../inst/extdata/06_boots/MLG129.csv"))
#   check  <- quackNboots(temp,
#                         nboots = 1000,
#                         distribution = "normal",
#                         type = "fixed",
#                         uniform = 1,
#                         mixtures = c("diploid", "triploid", "tetraploid"),
#                         samplename = "MLG129")
# 
#   write.csv(check, file = "../inst/extdata/06_boots/MLG129-boots.csv", row.names = FALSE)

## -----------------------------------------------------------------------------
MLG129boot <- read.csv("../inst/extdata/06_boots/MLG129-boots.csv")
MLG129boot