---
title: "Qploidy2 to nQuack"
output: rmarkdown::html_vignette
description: >
  How to use a standardized VCF from Qploidy2 in nQuack.
vignette: >
  %\VignetteIndexEntry{Qploidy2 to nQuack}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
resource_files:
  - ../inst/extdata/07_Qploidy2/output.csv
---

```{r message=FALSE, warning=FALSE, include=FALSE}
library(kableExtra)
library(dplyr)
library(nQuack)
library(data.table)
```

## Introduction to Qploidy2

Standardization and data cleaning are very important for determining ploidal level from sequence data. In 2025, Taniguti et al.  published their new tool that helps with this - [Qploidy](https://github.com/Cristianetaniguti/Qploidy/). Here we provide a tutorial for using data standardized with the extension of this approach, [Qploidy2](https://github.com/Breeding-Insight/Qploidy2/). 

Here, I followed their [tutorial on Alfalfa](https://breeding-insight.github.io/Qploidy2/Qploidy_alfalfa_tutorial.html#1_Introduction) and am using the output of `Qploidy2::standardize()`in nQuack. 

Find out more about this cool tool in their publication:

Taniguti, C. H., Lau, J., Hochhaus, T., Arias, D. C. L., Hokanson, S. C., Zlesak, D. C., Byrne, D. H., Klein, P. E., & Riera-Lizarazu, O. (2025). Exploring chromosomal variations in garden roses: Insights from high-density SNP array data and a new tool, Qploidy. The Plant Genome, e70044. \doi{10.1002/tpg2.70044}.




## Using data from Qploidy2 with nQuack 

## Predicting individual's ploidal level

Here, I followed the [tutorial on Alfalfa](https://breeding-insight.github.io/Qploidy2/Qploidy_alfalfa_tutorial.html#1_Introduction) to generate a file using `Qploidy2::standardize()`. To find out more about these steps, please see Qploidy2's [tutorial on Alfalfa](https://breeding-insight.github.io/Qploidy2/Qploidy_alfalfa_tutorial.html#1_Introduction).

```{r 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)

```

### Step 1: Modifying input from Qploidy2 to nQuack

Below, I subsample a data frame from an object of the class 'qploidy_standardization', generated above with `Qploidy2::standardize()`. I then split this data frame by SampleName and subset the information needed to use nQuack - the total coverage and coverage for a randomly sampled allele at every site which is biallelic for that individual. 

```{r 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)
}

```

### Step 2: Model inference 

Here we are following the [Basic Example](https://mlgaynor.com/nQuack/articles/BasicExample.html) and inferring ploidal level for the complete sample. If you are interested in a sliding window approach - we suggest identifying the most accurate model for your sample and then applying only this model with `bestquack()` on each sample, but [subsampled into windows](https://mlgaynor.com/nQuack/articles/FAQ.html#should-i-subsample-my-data). The sliding  to identify if there are regional differences in ploidal level.

Note, I wrote the output as I looped through the samples by using `data.table::fwrite()`.


```{r 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)
}

```

#### Identify the most accurate model

Using our function `quackit()`, you can easily interpret model output. Here we are selecting models based on the BIC score and only considering diploid and tetraploid mixtures. I only ran 5 samples for this example and assumed all samples were tetraploid - we sadly cannot identify the most accurate approach (distribution and type) here since many have 100% accuracy. 

```{r 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) 
```


