1 Introduction

This tutorial demonstrates how to run Qploidy2’s copy number estimation using the nQuack EM-based BAF model selection route. It uses the same alfalfa mapping population dataset as the main Qploidy2 tutorial and assumes that you have already completed the Standardization step described there (standardize()).

Important: The nQuack route requires allele read counts (X, Y, R columns) and is therefore exclusive to sequencing-based data (WGS, targeted sequencing such as DArTag/GTseq/AgriSeq, etc.). It cannot be used with intensity-based genotyping arrays (Axiom, Illumina), because arrays do not produce read counts. For array data, use hmm_estimate_CN() (the grid-based BAF model selection route) instead.

1.1 hmm_estimate_CN_nQuack(): a self-contained nQuack route

The nQuack integration is a separate, self-contained function, hmm_estimate_CN_nQuack(). It implements the same HMM structure as hmm_estimate_CN but using nQuack EM models to estimate both the overall sample ploidy estimate and per-window BAF likelihoods/dosage calling. Because nQuack only models copy-number states 2 through 6, cn_grid must be a subset of 2:6; use hmm_estimate_CN() instead for CN=0/1/LOH detection or ploidies outside 2:6.

Its signature mirrors hmm_estimate_CN() for the shared windowing/HMM arguments (cn_grid, chr, segment_zscore, exp_ploidy, etc.), plus three arguments specific to nQuack, all forwarded to nQuack::bestquack():

Argument Default Description
nQuack_dists "normal" The single distribution family nQuack evaluates: one of "normal", "beta", or "beta-binomial". nQuack’s BIC/LL values are not comparable across distribution families, so only one family can be passed per call; providing more than one, or NULL, raises an error. To choose a family, follow nQuack’s guidelines and validate candidates on samples with known ploidy first.
nQuack_type "fixed" The nQuack EM type: one of "fixed" (recommended), "fixed_2", or "fixed_3". BIC/LL are not comparable across types, so only one can be evaluated per call.
nQuack_uniform 0 1 to include an additional uniform mixture component (<distribution>-uniform), or 0 to omit it.

hmm_estimate_CN_multi() auto-detects which function to run per sample: if any of nQuack_dists, nQuack_type, or nQuack_uniform is passed via ..., it dispatches to hmm_estimate_CN_nQuack(); otherwise it runs hmm_estimate_CN().

nQuack must be installed before using this route:

install.packages("nQuack")

For full details on the grid-based route, see the Qploidy2 alfalfa tutorial. For background on nQuack itself, see Gaynor et al. (2024) and the nQuack vignettes.

2 Data Preparation

Follow Steps 1–4 of the main tutorial to obtain the standardized object. The code below reproduces those steps for completeness.

library(Qploidy2)
packageVersion("Qploidy2")
## [1] '1.22.2'

2.1 Download the alfalfa dataset

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 = "alfalfa_F1_marker_panel_dataset_publicly_available.vcf.gz")
vcf_path_local <- "alfalfa_F1_marker_panel_dataset_publicly_available.vcf.gz"

2.2 Read input files

data     <- qploidy_read_vcf(vcf_path_local)
genos    <- qploidy_read_vcf(vcf_path_local, geno     = TRUE)
geno.pos <- qploidy_read_vcf(vcf_path_local, geno.pos = TRUE)

2.3 Standardize

The nQuack integration requires the columns X (reference allele count), Y (alternative allele count), and R (total depth) that are present in the $data slot of the standardized object. Make sure your standardization step preserves them (the default standardize() call does).

qploidy_standardization <- standardize(
  data                   = data,
  genos                  = genos,
  geno.pos               = geno.pos,
  ploidy.standardization = 4,
  threshold.n.clusters   = 4,
  threshold.missing.geno = 0.9,
  min.depth              = 30,
  max.depth              = 600,
  n.cores                = 2,
  out_filename           = "standardization.tsv.gz",
  verbose                = TRUE
)

If you already have the standardized file from the main tutorial, reload it:

qploidy_standardization <- read_qploidy_standardization("standardization.tsv.gz")

3 Export nQuack input format from standardized BAF

hmm_estimate_CN_nQuack() builds the nQuack input matrix automatically (see the next section). If you instead want to inspect the data nQuack sees, or fit quackNormal(), quackBeta(), or quackBetaBinom() yourself outside of Qploidy2’s HMM, use convert2nQuack() to convert the standardized object into the list-of-matrices format these functions expect — one matrix per sample.

nquack_input <- convert2nQuack(
  qploidy_standardization_object = qploidy_standardization,
  min_depth              = 10,
  seed                   = 123
)

length(nquack_input)
## [1] 185
head(nquack_input[[1]])
##        R allele
## [1,] 111     31
## [2,] 100     47
## [3,] 182    136
## [4,] 320    236
## [5,] 102     55
## [6,] 583    292

By default (from_standardized_baf = TRUE), convert2nQuack() derives the allele counts fed to nQuack from the standardization-corrected baf column rather than the raw X/Y read counts. For each retained heterozygous site it inverts the standardized BAF back into reference/alternative counts:

  • alt_count <- round(baf * R)
  • ref_count <- R - alt_count

and then randomly keeps one of the two (reference or alternative) with equal probability — the same approach used internally by hmm_estimate_CN_nQuack(). This ensures nQuack sees allele frequencies that have already been corrected for marker/assay biases by standardize(), which matters most for data where raw X/Y ratios can be systematically skewed per marker.

If you would rather fit nQuack on the raw, uncorrected allele counts, set from_standardized_baf = FALSE:

nquack_input_raw <- convert2nQuack(
  qploidy_standardization_object = qploidy_standardization,
  min_depth              = 10,
  from_standardized_baf  = FALSE
)

Each element of the returned list is a two-column integer matrix (one row per qualifying site): column 1 is the total read depth R, column 2 is the randomly selected allele count. You can fit any of the nQuack models directly on this output, e.g.:

nQuack::quackBeta(xm = nquack_input[[1]], samplename = names(nquack_input)[1], cores = 1)

See ?convert2nQuack for the full filtering details (homozygous-site removal, minimum depth, genotype dosage filter).

4 How nQuack Input Is Built Inside Qploidy2

You do not need to prepare nQuack input yourself. hmm_estimate_CN_nQuack() builds the nQuack input matrix automatically from the $data slot of your standardized object.

For transparency, here is what happens internally for each sample:

Step 1 — Discard homozygous markers. Only markers where the sample is heterozygous are kept. A site is considered heterozygous when its allele ratio is neither 0 (all-reference) nor 1 (all-alternative), and its dosage call (geno) is neither 0 (homozygous reference) nor the maximum dosage for the dataset (homozygous alternative). In practice this means markers where the sample carries at least one copy of each allele. Homozygous markers are discarded because their allele frequency provides no information for distinguishing ploidy levels — only heterozygous sites generate the mixture of allele-frequency peaks that nQuack models.

Step 2 — Remove low-coverage sites. Sites with total read depth R < 10 are discarded. At very low depths, allele counts are too noisy to reliably estimate allele frequencies. The threshold of 10 reads is nQuack’s recommended minimum for meaningful allele-frequency estimation.

Step 3 — Randomly sample one allele per site. For each remaining heterozygous site, either the reference (X) or the alternative (Y) allele count is randomly drawn with equal probability (P = 0.5). This produces a matrix with two columns:

  • column 1: R, the total read depth at that site
  • column 2: sampled_count, the read count for the randomly chosen allele (either the reference count X or the alternative count Y)

The ratio sampled_count / R gives one allele frequency observation per site. Because the reference and alternative alleles are sampled randomly, the allele frequencies cluster symmetrically between 0 and 0.5 regardless of which allele is the reference — this is what nQuack’s mixture models expect.

Step 4 — Fit nQuack models. The two-column matrix is passed to quackNormal, quackBeta, and quackBetaBinom. Each function fits mixture models for all candidate ploidies (from cn_grid, restricted to the range 2–6 supported by nQuack). The model with the lowest BIC across all families and all ploidies is declared the winner, and its ploidy becomes exp_ploidy for the HMM.

The random allele sampling in Step 3 means results can vary slightly between runs. Set set.seed() before calling hmm_estimate_CN_nQuack() for reproducibility.

5 BAF Model Selection and HMM: the nQuack Route

5.1 Single sample

hmm_estimate_CN_nQuack() lets nQuack handle the entire BAF model selection step. It accepts exactly one (nQuack_dists, nQuack_type, nQuack_uniform) combination per call — it does not test multiple distributions, types, or uniform settings and pick the best by BIC, since nQuack’s BIC/LL values are not comparable across distribution families, EM types (fixed/fixed_2/fixed_3), or uniform/non-uniform mixtures (see Restricting distribution families with nQuack_dists below). If you are unsure which nQuack_dists/nQuack_type/nQuack_uniform to use, follow nQuack’s own guidance: evaluate candidate combinations on samples with known ploidy, as described in the nQuack tutorials, and then pass the validated combination here. Qploidy2 uses that distribution and the winning ploidy to build the HMM emission parameters.

set.seed(42)
res_nq <- hmm_estimate_CN_nQuack(
  qploidy_standarize_result = qploidy_standardization,
  sample_id                 = "85-126",
  cn_grid                   = c(2, 3, 4, 5, 6),
  nQuack_dists = "normal", 
  nQuack_uniform = 1,
  nQuack_type = "fixed",
  verbose                   = TRUE
)

res_nq
## hmm_CN result
##   Copy-number grid: 2, 3, 4, 5, 6 
##   Expected ploidy: 4 
##   Minimum SNPs per window: 31 
##   Initial state probabilities (pi0): 0.125, 0.125, 0.5, 0.125, 0.125 
##   Estimated z means per CN: -0.093, 0.299, 0.572, 1.347, 1.711 
##   Estimated z mean: 0.6914988 
##   Estimated z sigma: 0.215 
##   BAF Emission distribution: normal-uniform 
##   Final log-likelihood: 0 
##   Range of CN_call (window-level): 4 - 5

Inspect the karyotype:

plot_karyotype(res_nq, sample_name = "85-126", nrow = 2)

Color by HMM posterior confidence:

plot_karyotype(res_nq, sample_name = "85-126", nrow = 2, color_by = "post_max_CN")

Detailed diagnostic track:

p <- plot_cn_track(hmm_CN      = res_nq,
                   data        = data,
                   geno.pos    = geno.pos,
                   sample_id   = "85-126",
                   show_window_lines = FALSE,
                   z_by_mean   = TRUE)
p$arranged

5.2 Restricting distribution families with nQuack_dists

nQuack_dists accepts exactly one distribution family: "normal" (default), "beta", or "beta-binomial". nQuack’s own documentation states that BIC/LL values are not comparable across distribution families, so passing more than one family (or NULL to test all three) is rejected with an error. To choose a family, follow nQuack’s guidelines: validate candidate families on samples with known ploidy first, then pass that single family here. nQuack_type ("fixed", "fixed_2", or "fixed_3") and nQuack_uniform (0 or 1, to exclude/include a uniform mixture component) are similarly restricted to a single value per call:

# Normal family, fixed type, with a uniform mixture component
res_nq_normal <- hmm_estimate_CN_nQuack(
  qploidy_standarize_result = qploidy_standardization,
  sample_id                 = "85-126",
  cn_grid                   = c(2, 3, 4, 5, 6),
  nQuack_dists              = "normal",
  nQuack_type               = "fixed",
  nQuack_uniform            = 1
)

# Beta-Binomial family, fixed type, with a uniform mixture component
res_nq_seq <- hmm_estimate_CN_nQuack(
  qploidy_standarize_result = qploidy_standardization,
  sample_id                 = "85-126",
  cn_grid                   = c(2, 3, 4, 5, 6),
  nQuack_dists              = "beta-binomial",
  nQuack_type               = "fixed",
  nQuack_uniform            = 1
)

5.3 Multiple samples

Run the nQuack route for all samples in parallel. Because nQuack fits a mixture model per sample, this takes longer than the grid route. hmm_estimate_CN_multi()’s n_cores argument controls how many samples are processed simultaneously.

hmm_estimate_CN_multi() detects the nQuack route automatically: passing any of nQuack_dists, nQuack_type, or nQuack_uniform (via ...) is enough to dispatch each sample to hmm_estimate_CN_nQuack() instead of hmm_estimate_CN().

set.seed(42)
multi_esti_nq <- hmm_estimate_CN_multi(
  qploidy_standarize_result = qploidy_standardization,
  sample_ids                = "all",
  n_cores                   = 3,
  cn_grid                   = c(2, 3, 4, 5, 6),
  nQuack_dists              = "normal",
  nQuack_type               = "fixed",
  nQuack_uniform            = 1
)
samples <- unique(qploidy_standardization$data$SampleName)
compare_cn_track(hmm_CN          = multi_esti_nq,
                 samples_to_plot = samples[1:30], 
                 gray_CN = 4,
                 facet_ncol      = 8)

6 Comparing Routes: nQuack vs. Grid

For a well-behaved tetraploid dataset the two routes typically agree. The table below summarises their trade-offs:

Feature Grid (hmm_estimate_CN()) nQuack (hmm_estimate_CN_nQuack())
BAF model selection Grid search over KDE bandwidth nQuack EM for ploidy + all BAF parameters
Speed Faster Slower
Mixture families tested Parametric templates Normal, Beta, Beta-Binomial
Platform compatibility Array (intensity-based BAF) and sequencing (read-count-derived ratio) Sequencing only (requires raw X/Y/R read counts)
CN range Any (including 0/1/LOH) 2–6 only

7 Ploidy type summary

types_df <- count_types(multi_esti_nq)
types_df
## 
## Sample classification summary
## ====================================================
##   2x  (n = 11)
##     2x - euploid            10
##     2x - segmental          1
##   3x  (n = 3)
##     3x - euploid            1
##     3x - aneuploid          1
##     3x - segmental          1
##   4x  (n = 168)
##     4x - euploid            118
##     4x - aneuploid          30
##     4x - segmental          20
##   5x  (n = 3)
##     5x - euploid            2
##     5x - aneuploid          1
## ----------------------------------------------------
##   Total euploid              131
##   Total aneuploid            32
##   Total segmental aneuploidy 22
##   Grand total                185
## ====================================================
compare_cn_track_summary(
  multi_esti_nq, 
  gray_CN = 4,
  samples_to_plot = c(types_df$aneuploid, types_df$segmental)
)

7.1 Per-sample, per-chromosome summaries

summ_chr    <- summarize_cn_mode(multi_esti_nq, level = "chromosome")
summ_sample <- summarize_cn_mode(multi_esti_nq, level = "sample")

head(summ_chr)
head(summ_sample)

7.2 Write results

write_hmm_CN(multi_esti_nq, prefix = "alfalfa_nQuack")

Reload in a later session:

multi_esti_nq <- read_hmm_CN(
  by_marker_file = "alfalfa_nQuack_by_marker.csv.gz",
  by_window_file = "alfalfa_nQuack_by_window.csv.gz",
  params_file    = "alfalfa_nQuack_params.rds"
)

7.3 MAPpoly export

for_mappoly <- export_mappoly(multi_esti_nq)
head(for_mappoly)

8 Troubleshooting nQuack Integration

8.1 nQuack is not installed

Error: nQuack is not installed. Install with:
  install.packages("nQuack")

Install nQuack and restart R before re-running.

8.2 Array data — nQuack route not applicable

If your data comes from a hybridization-based genotyping array (Axiom, Illumina), the X, Y, and R columns will contain probe intensities rather than read counts, and the nQuack model assumptions will not hold. Use hmm_estimate_CN() (the grid-based BAF model selection route) for array data instead.

8.3 Too few heterozygous sites

nQuack needs a minimum number of heterozygous sites (at least 10 after filtering) to fit mixture models. If a sample has too few informative sites, hmm_estimate_CN_nQuack() stops with an error instead of estimating ploidy — there is no template fallback on this route. Either supply exp_ploidy explicitly, or use hmm_estimate_CN() (which does support a grid-based template fallback) for such samples. Low heterozygosity samples (< 5 %) are typically consistent with haploid (CN = 1) or highly inbred lines.

8.4 nQuack calls an unexpected ploidy (high-variance samples)

On samples with very wide BAF distributions (high allele-frequency variance, often due to low sequencing depth or complex rearrangements), nQuack’s fixed-mean normal mixture model may favour a higher ploidy by BIC simply because more components absorb variance better. There is currently no reliable automated safeguard against this: a one-component “diploid sigma” check was previously used to flag and force such calls to diploid, but that check could not distinguish genuine noise from real separation between dosage peaks (a noiseless tetraploid sample, for example, can produce a sigma well above the override’s threshold purely from its three distinct peaks at 0.25/0.5/0.75) and has been removed.

If you suspect an inflated ploidy call, inspect the sample’s BAF distribution directly (e.g. with plot_baf_hist()) and consider improving sequencing depth or reviewing the standardization quality for that sample.

8.5 Reproducibility

nQuack’s input preparation randomly samples one allele per heterozygous site (Step 3 in the input-building pipeline above). Set a seed before calling hmm_estimate_CN_nQuack() or hmm_estimate_CN_multi() for reproducible results:

set.seed(42)
res <- hmm_estimate_CN_nQuack(...)

8.6 nQuack supports copy numbers 2–6 only

hmm_estimate_CN_nQuack() only models CN states 2 through 6 (nQuack’s internal emstep* functions have no representation for CN=1/LOH or CN > 6), and it has no template fallback for out-of-range states: cn_grid (excluding 0, which is also rejected) must be a subset of 2:6, or the function stops with an informative error — including the default cn_grid, which includes CN=1 and must be restricted explicitly.

If you need CN=1/LOH detection or ploidies outside 2:6, use hmm_estimate_CN() (the grid/template route) instead — it has no such restriction.

9 How to Cite

9.0.1 Qploidy — Standardization

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. https://doi.org/10.1002/tpg2.70044

9.0.2 Qploidy2 — HMM and BAF model selection

Manuscript in preparation. Please contact the author for more information.

9.0.3 nQuack — EM-based BAF model selection

Gaynor, M., Landis, J., O’Connor, T., Laport, R., Doyle, J., Soltis, D., Ponciano, J., & Soltis, P. (2024). nQuack: An R package for predicting ploidal level from sequence data using site-based heterozygosity. Applications in Plant Sciences, 12(4), e11606. https://doi.org/10.1002/aps3.11606

10 Acknowledgments

10.0.1 Qploidy version 1.0.0

Funded in part by the Robert E. Basye Endowment in Rose Genetics, Dept. of Horticultural Sciences, Texas A&M University, and USDA’s National Institute of Food and Agriculture (NIFA), Specialty Crop Research Initiative (SCRI) projects: “Tools for Genomics-Assisted Breeding in Polyploids: Development of a Community Resource” (Award No. 2020-51181-32156); and “Developing Sustainable Rose Landscapes via Rose Rosette Disease Education, Socioeconomic Assessments, and Breeding RRD-Resistant Roses with Stable Black Spot Resistance” (Award No. 2022-51181-38330).

10.0.2 Qploidy versions > 1.0.0 and Qploidy2

Supported by Breeding Insight.