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

## ----eval=FALSE, message=FALSE, warning=FALSE, include=TRUE-------------------
# # Load Package
# library(Qploidy2)
# 
# ## Download Data
# vcf_path_web <- "https://github.com/Breeding-Insight/BIGapp-PanelHub/raw/refs/heads/long_seq/alfalfa/GenoBrew_example/alfalfa_F1_marker_panel_dataset_publicly_available.vcf.gz"
# download.file(vcf_path_web, destfile = "inst/extdata/07_Qploidy2/alfalfa_F1_marker_panel_dataset_publicly_available.vcf.gz")
# vcf_path_local <- "inst/extdata/07_Qploidy2/alfalfa_F1_marker_panel_dataset_publicly_available.vcf.gz"
# 
# ## Read in data
# data <-  Qploidy2::qploidy_read_vcf(vcf_path_local)
# genos <-  Qploidy2::qploidy_read_vcf(vcf_path_local, geno = TRUE)
# geno.pos <-  Qploidy2::qploidy_read_vcf(vcf_path_local, geno.pos = TRUE)
# 
# ## Standardize
# qploidy_standardization <- Qploidy2::standardize(data = data,
#                                                  genos = genos,
#                                                  geno.pos = geno.pos,
#                                                  ploidy.standardization = 4,
#                                                  threshold.n.clusters = 4,
#                                                  n.cores = 2,
#                                                  out_filename = "../inst/extdata/07_Qploidy2/standardization.tsv.gz",
#                                                  verbose = TRUE)
# 

## ----eval=FALSE, message=FALSE, warning=FALSE, include=TRUE-------------------
# ## Read in the output
# qploidy_standardization <- Qploidy2::read_qploidy_standardization("../inst/extdata/07_Qploidy2/standardization.tsv.gz")
# 
# ## Subsample to just the data frame with coverage information
# qdata <- qploidy_standardization$data
# 
# ## Here we are interested in pulling out the information we need for nQuack - the total coverage and coverage for a randomly sampled allele for each sample.
# ### In this data frame:
# ### X = coverage of the reference
# ### Y = coverage of the alternative
# ### R = total coverage
# 
# ## Make a list of samples
# samples <- unique(qdata$SampleName)
# templist <- c()
# 
# for(i in 1:length(samples)){
#   temp <- qdata[which(qdata$SampleName == samples[i]), ]
#   ## Remove sites that are not biallelic for the individual
#   temp <- temp[which(temp$ratio != 1 & temp$ratio != 0),]
#   temp <- temp[which(temp$geno != 0 & temp$geno != 4), ]
#   ## R = total coverage, X = coverage of the reference, & Y = coverage of the alternative
#   xmdf <- data.frame(temp$R, temp$X, temp$Y)
#   xmr <- matrix(nrow = nrow(xmdf), ncol = 2)
#   ## Randomly select coverage of the reference or the alternative allele
#   for(y in 1:nrow(xmdf)){
#     xmr[y, 1] <- xmdf[y, 1]
#     xmr[y, 2] <- xmdf[y, sample(x = c(2,3), size = 1, prob = c(0.5, 0.5))]
#   }
#   ## Remove sites with coverage less than 10 and only keep biallelic sites
#   xmr <- xmr[which(xmr[,1] >= 10 & xmr[,2] != 0), ]
#   templist[[i]] <- as.matrix(xmr)
# }
# 

## ----eval=FALSE, message=FALSE, warning=FALSE, include=TRUE-------------------
# for(i in 1:length(samples)){
#   out1 <- quackNormal(xm = templist[[i]],
#                       samplename = samples[i],
#                       cores = 10,
#                       parallel = TRUE)
#   data.table::fwrite(out1, file = "../inst/extdata/07_Qploidy2/output.csv",
#                      append = TRUE)
#   out2 <- quackBeta(xm = templist[[i]],
#                     samplename = samples[i],
#                     cores = 10,
#                     parallel = TRUE)
#    data.table::fwrite(out2,
#                      file = "../inst/extdata/07_Qploidy2/output.csv",
#                      append = TRUE)
#   out3 <- quackBetaBinom(xm = templist[[i]],
#                          samplename = samples[i],
#                          cores = 10,
#                          parallel = TRUE)
#    data.table::fwrite(out3,
#                      file = "../inst/extdata/07_Qploidy2/output.csv",
#                      append = TRUE)
# }
# 

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.align='center'--------------
modoutput <-  read.csv("../inst/extdata/07_Qploidy2/output.csv")
summary <- c()
samples <- unique(modoutput$sample)
for(i in 1:5){
  temp <- modoutput[which(modoutput$sample == samples[i]), ]
  summary[[i]] <- nQuack::quackit(model_out =  temp, 
                     summary_statistic = "BIC", 
                     mixtures = c("diploid", "tetraploid"))
}
alloutputcombo <- do.call(rbind, summary)
alloutputcombo <- alloutputcombo %>%
                  dplyr::mutate(accuracy = ifelse(winnerBIC == "tetraploid", 1, 0))
sumcheck <- alloutputcombo %>% 
            group_by(Distribution, Type) %>% 
            summarize(total = n(), correct = sum(accuracy))

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