Qploidy2 Tutorial
Large-scale copy number estimation in alfalfa mapping population using DArTag genotyping data
Breeding Insight
2026-07-28
1 Introduction
This tutorial demonstrates a typical workflow for ploidy and aneuploidy estimation using the Qploidy2 R package. We will cover data import, standardization, model selection, Hidden Markov Model (HMM)-based copy number estimation, visualization, and VCF export.
The dataset used here for our demonstration is an alfalfa mapping population genotyped using DArTag technology. Refer to this publication for more information about the dataset.
To start, download the example dataset to your local directory:
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 Remarks
The Qploidy standardization methodology is effective under the following conditions:
- Your marker data originates from Axiom or Illumina genotyping arrays
- Your marker data is derived from targeted sequencing platforms (e.g., DArTag, GTseq, AgriSeq)
- Your marker data is derived from whole-genome sequencing platforms
- All DNA samples were prepared following the same library preparation protocol (check this section for more details) For example: If you extracted DNA and sequenced two plates (192 samples) as one batch, and later sequenced an additional three plates (288 samples) as a second batch, you would need to analyze the two batches separately in Qploidy. Combining all 480 samples into a single analysis will lead to incorrect results due to batch effects.
- You know the ploidy of at least a subset (rough minimum of 60 samples) or you know the most common ploidy in the dataset. The higher the number of samples with known or common ploidy, the better the results
- Your dataset includes heterozygous samples
The Qploidy2 B-allele frequency (BAF) model selection and HMM methodologies are effective under the following conditions:
- The ratio or BAF distributions of your samples present clear peaks corresponding to different copy numbers
- The total read depth or z-scores variations within a sample follows a normal distribution, with the mean shifting according to the copy number
- The sample includes heterozygous loci, which are essential for BAF-based ploidy estimation. Samples with heterozygosity lower than 5% will be considered haploid (CN = 1) and they cannot be distinguished from inbred lines with higher CN
- If allele counts are based on a sequencing platform, the sequencing should have enough depth to capture dosages of the highest possible copy number in the dataset
- The number of markers is sufficient to provide a good resolution for the BAF distributions
- The larger the amount of markers, the smaller the copy number variations that can be detected
3 Install
- It requires R version superior to 3.6.0
Install Qploidy2 (if not already installed):
4 Input data
Qploidy2 receives three main objects named: -
marker-level data (data): should contain sample names,
marker ID, allele counts/intensities ratio and total for all samples -
marker positions (geno.pos): should contain marker IDs, and
respective genomic chromosome and positions - genotype dosages
(genos) - only required for standardization: contain
samples names, marker ID, and dosage/genotype considering expected
ploidy. The samples in this file should be the ones with a known ploidy
or the full dataset with dosage called considering the predominant
ploidy.
Qploidy2 has functions to facilitate the conversion of
standard file formats to these three objects. In this tutorial we will
show the qploidy_read_vcf function for VCF conversion. The
VCF must contain the AD field for allele counts. Check function
?read_axiom and ?read_illumina_array if your
data comes from genotyping arrays.
library(Qploidy2)
# Read marker-level data from VCF
data <- qploidy_read_vcf(vcf_path_local)
head(data)# Read genotype dosages (optionally, restrict to samples with known ploidy)
genos <- qploidy_read_vcf(vcf_path_local, geno = TRUE)
head(genos)# Read marker positions
geno.pos <- qploidy_read_vcf(vcf_path_local, geno.pos = TRUE)
head(geno.pos)5 Check data initial status
The plot_raw function provides diagnostic plots for
metrics that will be considered by Qploidy2 copy number
estimation. Here is a brief description of each from top to bottom:
- Heterozygosity: y-axis = proportion of heterozygous loci in a
genomic window with size defined by argument
window_size(default = 2Mb); x-axis = genomic positions - Allele counts/intensities ratio (B/(A+B)): y-axis = ratio value; x-axis = genomic positions
- Ratio histogram split by chromosome: y-axis = observation counts; x-axis = ratio value
- Ratio histogram considering all markers of the samples: y-axis = observation counts; x-axis = ratio value
- Read depth (sequencing platforms) or total intensities (array) in y-axis and genomic position in x-axis
Check the plots for a few samples. If the histograms do not show
clear peaks, consider applying Qploidy standardization. If
you already have clear peaks and do not want to try to improve them,
skip to item # BAF Model Selection.
5.1 Check for batch effect before standardization
Qploidy2 standardization assumes that each marker behaves consistently across samples, while allowing markers to differ from one another. Depending on the genotyping platform, however, this assumption may not hold — particularly if DNA extraction protocols were modified or samples were processed in different plates or batches.
For example, you may have extracted DNA at one time point and sequenced a plate with 96 samples in a given year, then used the stored DNA two years later to sequence an additional 96 samples. Changes in DNA quality over time can affect probe performance, leading to systematic differences between batches. In such cases, standardization should be performed separately for each batch.
You can assess whether batch effects are present in your dataset using the following function:
# provide a passport file (containing year of DNA extraction/plate number/ or any possible factor for a batch effect)
# This example dataset doesn't have a passport available, we will simulate one just to exemplify the function usage:
set.seed(123)
samples <- unique(data$SampleName)
passport_file <- data.frame(SampleName = samples, year_seq = sample(rep(c(2020, 2022, 2025), length(samples))))
head(passport_file)p <- pca_plot(input = data,
geno.pos = geno.pos,
passport_file = passport_file,
group_column = "year_seq",
sampleID_column = "SampleName",
col2use = "R")
p$PC1xPC2
The function draws a PCA using the total depth/intensities values of each marker x sample to identify batch effect between samples.
In this example, we did not detect evidence of batch effects, as samples did not cluster according to genotyping year. Therefore, we can proceed with standardization using the combined dataset.
If your dataset formed clear groups like the figure below referring to another dataset not covered in this tutorial:
In this case, perform the standardization separately for each and
merge the result using merge_qploidy_datas function.
You can also verify the effect of batches considering a single marker and all samples:
Again, for this dataset, no difference is observed but here are some markers examples from another dataset that has batch effect:
Dosage call methods such as fitPoly, updog,
and polyRAD use these type of ratio distributions (one
marker and all samples) to estimate each dosage cluster. If you have
batch effect between your runs (as the figures above), they will also
affect the quality of dosage call using these software. Therefore, you
should conduct dosage/genotype call individually for each batch.
5.2 Initial screening for chromosome abnormalities (optional)
For some genotyping platforms, total read depth is consistent across samples. When this assumption holds for your dataset, deviations in a sample’s total depth relative to the population average can provide an initial indication of copy number variation.
You can rapidly screen for candidate samples using the function below:
chromosome_level_test_plot_qploidy(input = data,
geno.pos = geno.pos,
selected_samples = samples[1:12], # plot the first 12
col2test = "R")
Before proceeding to the standardization step, consider removing abnormal chromosomes from your reference dataset (genos).
5.3 Overview of samples heterozygosity
# interactive graphic if plotly package is installed
p <- plot_heterozygosity(data, col2use = "ratio")
pThis is an interactive graphic; hover your mouse over the cells to get information on sample name, heterozygosity percentage, and the number of markers used for this estimation. Heterozygosity is an important metric for Qploidy2, as it uses the allele ratio heterozygous peaks as part of the input data for estimating ploidy. Samples with heterozygosity lower than 5% will be considered haploid (CN = 1) and cannot be distinguished from inbred lines with higher CN.
6 Standardization
Standardization of markers raw read counts or intensity data is essential for reducing technical noise and improving the accuracy of ploidy estimation. Markers from diverse platforms have different characteristics, and standardization helps to make them comparable.
The function has several parameters that can be adjusted based on the
characteristics of your dataset. You can access a description of all
parameters in ?standardize. Here we will focus only a few
of them:
Parameters: - genos: This data.frame
defines which samples are your reference for the standardization. If you
have samples with known ploidy, filter this object to keep only them.
You can use the command:
# Demonstration code, not applied to current dataset
# known_tetraploids <- paste0("Sample",1:100)
# str(known_tetraploids)
# genos <- genos[which(genos$SampleName %in% known_tetraploids),]If you are not sure about ploidy of at least 60 samples, just keep
all them and define the most common ploidy as
ploidy.standardization.
ploidy.standardization: Sets the predominant ploidy for the dataset defined in the objectgenos(e.g., 2 for diploid, 4 for tetraploid). This guides clustering and normalization, and should match the known or most frequent ploidy in your samples for best results.threshold.n.clusters: Minimum number of genotype clusters required for a marker to be retained. The maximum is defined as the possible number of the dosages according to the ploidy defined onploidy.standardization(ploidy + 1). Lower values requires imputation of dosage clusters for missing genotypes. The imputation is performed based on distance to other clusters, so it is recommended to have at least 3 clusters for a good imputation. If too many markers are filtered, relax this threshold (e.g., to 4 for tetraploids) to retain more markers with imputed clusters.
Other the the minimum number of clusters, the standardization step also apply some other filters:
threshold.missing.geno: Maximum fraction of missing datapoints allowed per marker. Markers with a higher fraction will be removed.min.depth: Maximum fraction of missing datapoints allowed per marker. Markers with a higher fraction will be removed.max.depth: Maximum read depth (R) threshold. Datapoints with R above this value will have their allele ratio set to missing and be excluded from BAF and Z-score computations.
If your data comes from sequencing (and not hybridization), filtering by depth is important for quality control. To decide on the minimum and maximum of depth, you may want to take a quick look on the total depth distribution:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 60.0 139.0 184.6 253.0 2635.0The threshold should be around the mean but giving enough space for the variations:
qploidy_standardization5 <- standardize(
data = data,
genos = genos,
geno.pos = geno.pos,
ploidy.standardization = 4, # Set to known or most frequent ploidy
threshold.n.clusters = 5, # Minimum clusters for marker retention
threshold.missing.geno = 0.9,
min.depth = 30,
max.depth = 600,
n.cores = 2,
out_filename = "standardization.tsv.gz",
verbose = TRUE
)
# Check filtering summary
qploidy_standardization5## This is an object of class 'qploidy_standardization'
## --------------------------------------------------------------------
## Parameters
##
## 1 Standardization type: intensities
## 2 Ploidy: 4
## 3 Min # of heterozygous classes (clusters) present: 5
## 4 Max proportion of missing genotype by marker: 0.9
## 5 Min genotype probability: 0.8
## 6 Min read depth (min.depth): 30
## 7 Max read depth (max.depth): 600
## --------------------------------------------------------------------
## Filters
##
## 1 # markers at raw data: 2766 (100%)
## 2 % datapoints filtered by low probability: - (0 %)
## 3 # datapoints filtered by depth (min/max): 84309 -
## 4 # markers filtered by missing data: 184 (6.65 %)
## 5 # markers filtered by min number of clusters: 2323 (83.98 %)
## 6 # markers filtered by lack of genomic information: 0 (0 %)
## 7 # markers remaining with estimated BAF: 259 (9.36 %)
## 8 # markers remaining with estimated z-score: 2664 (96.31 %)The graphics below provides a visual explanation of how the standardization worked for 6 different markers:
Observe three different situations:
- The second, fourth, fifth and sixth markers: these markers had at
least two samples included in the
genosdata.frame above (reference samples) for each dosage expected considering the ploidy defined atploidy.standardization = 4. Therefore, standardization could be performed for all the dosages. It re-position the centroids of ratios dosage clusters (red dots) to the expected position (blue dots). After standardization every markers have dosage clusters centroids in same positions, making possible to combine them for sample, chromosome or chromosome region level BAF distribution evaluation. - The first marker don’t have samples representing three of the dosages (2,3,4) and therefore the marker was filtered out and will not contain a BAF for downstream analysis (but kept for z-score).
- The third marker in the figure don’t have samples with dosage 4,
missing a single cluster but also discarded given that the defined
argument
threshold.n.clusters = 5requires at least two samples for all five clusters.
PS: We use B-allele frequency BAF terminology to refer to standardized ratios.
If too many markers are filtered due to number of clusters, you relax the cluster threshold:
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
)
qploidy_standardization## This is an object of class 'qploidy_standardization'
## --------------------------------------------------------------------
## Parameters
##
## 1 Standardization type: intensities
## 2 Ploidy: 4
## 3 Min # of heterozygous classes (clusters) present: 4
## 4 Max proportion of missing genotype by marker: 0.9
## 5 Min genotype probability: 0.8
## 6 Min read depth (min.depth): 30
## 7 Max read depth (max.depth): 600
## --------------------------------------------------------------------
## Filters
##
## 1 # markers at raw data: 2766 (100%)
## 2 % datapoints filtered by low probability: - (0 %)
## 3 # datapoints filtered by depth (min/max): 84309 -
## 4 # markers filtered by missing data: 184 (6.65 %)
## 5 # markers filtered by min number of clusters: 1801 (65.11 %)
## 6 # markers filtered by lack of genomic information: 0 (0 %)
## 7 # markers remaining with estimated BAF: 781 (28.24 %)
## 8 # markers remaining with estimated z-score: 2664 (96.31 %)In this case, Qploidy inputs the position of the missing cluster based on the mean distance between the existing ones making possible to perform the standardization. In our visual example, the third marker is kept after the change:
The function also generates the Qploidy file
standardization.tsv.gz that you can reload in another
session to recover the qploidy_standardization object:
This file can also be uploaded to GenoBrew app for
results curation. See section
Interactive curation of results.
7 Quality Assessment and Visualization
Visualize the effect of standardization for each sample and explore BAF and Z-score distributions:
plot_qploidy_standardization(x = qploidy_standardization,
sample = "J432",
type = c("Ratio_hist_overall", "BAF_hist_overall"))
See ?plot_qploidy_standardization for more plot
types.
7.1 Chromosome-level screening for abnormalities based on Z-score (optional)
You can use the same quick check we used before on the raw total depth (R) with the standardized z-scores. This is useful to identify samples with potential aneuploidies or large CNVs before running the HMM, and also to select samples with known ploidy to guide the BAF model selection.
chromosome_level_test_plot_qploidy(input = qploidy_standardization,
selected_samples = samples[1:12],
col2test = "z")
8 BAF Model Selection
Select the best BAF model for a sample. The approach implemented here
for this step is based on nQuack
methods. It tests different statistical distribution mixture models
fitting for each tested ploidy on the BAF distributions and selects the
best based on BIC. There are two ways of performing it, one “quick and
dirt” based on comparing parameters in a grid and another using
nQuack EM functions:
- EM approach BAF model selection: has a robust EM algorithm able to return tailored parameters for the selected distribution to users
- Grid approach BAF model selection: does not use an EM algorithm but a grid of parameters to be tested, which is less time- and computationally demanding but may not provide the best parameters for this step
This current tutorial apply the “grid approach”. Check tutorial Qploidy2 and nQuack integration for running the EM approach using nQuack and a comparison between the methods results.
This step helps determine the most likely ploidy and mixture model
parameters for BAF. The selected approach and model will be used in the
Qploidy2 HMM for copy number estimation.
You can select a sample with known ploidy (e.g., from the previous step) to guide the model selection, or test multiple samples and compare results.
8.1 Applying to raw data
sample_ratios <- data$ratio[data$SampleName == "J432"]
selected_model_raw <- select_best_baf_model(baf_vec = sample_ratios,
cn_grid= c(2,3,4,5,6),
bw = c(0.02, 0.04,0.1),
plot = TRUE)
selected_model_raw$plot
See that despite selecting as the best a model that inferred the ploidy as 4, the second place inferred as 3 and the likelihood between models don’t differ much. For this dataset, standardization is required for better results. This step and other downstream steps only make this move evident.
8.2 Applying to standardized data
selected_model <- select_best_baf_model(qploidy_standardization = qploidy_standardization,
sample = "J432",
cn_grid= c(2,3,4,5,6),
bw = c(0.02, 0.04,0.1),
plot = TRUE)
selected_model$plot
Using standardized model, the four best models returned ploidy 4 proving to be more stable than without standardization.
9 HMM Copy Number Estimation (Single Sample)
Estimate copy number for a sample using a Hidden Markov Model. This
step provides a robust approach that integrates raw or standardized
ratios (BAF) and total depth (Z-score) information across genomic
regions to infer CNV segments. Find all parameters description in
?hmm_estimate_CN.
9.1 Applying to raw data
res_raw <- hmm_estimate_CN(data = data,
geno.pos = geno.pos,
selected_model = selected_model_raw,
use_values = c("ratio","R"),
sample_id = "J432",
cn_grid = c(2,3,4,5,6))
# View results
res_raw## hmm_CN result
## Copy-number grid: 2, 4, 5, 6
## Expected ploidy: 4
## Minimum SNPs per window: 27
## Initial state probabilities (pi0): 0, 1, 0, 0
## Estimated z means per CN: 130.635, 200.22, 236.76, 357.835
## Estimated z mean: 200.222
## Estimated z sigma: 48.638
## BAF Emission distribution: negative_binomial
## Final log-likelihood: 0
## Range of CN_call (window-level): 4 - 4Not that the HMM multipoint approach consider the marker position order based on the reference genome used. Therefore, changing the reference genome may change the results.
Warning: Avoid to use unlikely copy numbers for your specie to decrease the number of tested scenarios.
Qploidy2 has different ways to visualize the estimated
karyotype:
You can add color by the CN call final likelihood:
The color show how confident the model was on calling the CN for each of the markers based on the combination of R/z-score and ratio/BAF.
You can also add color by the BAF weight:
This metric shows if the CN for the window was inferred based on the BAF/ratio given the presence of heterozygous genotypes within it, or if it was just based on the total depth/z-score given the lack of heterozygosity. The most confident estimations are the ones that are based on both BAF/ratio and R/z-score - they have BAF weight of 0.5 (and 0.5 for z-score weight).
The function plot_cn_track provides a detailed
diagnostic of the process:
p <- plot_cn_track(hmm_CN = res_raw, data = data,
geno.pos = geno.pos,
sample_id = "J432",
show_window_lines = FALSE,
z_by_mean = TRUE)
# All panels together
p$arranged
Despite not having the same peak resolution as the standardized data, the HMM was robust enough to retun correct result (CN 4) for all chromosomes. The standardization will be more critical for more complex scenarios, such as messy ratios, and/or presence of loss or gain of chromosomes. We have some example further in this tutorial.
9.2 Applying to standardized data
res <- hmm_estimate_CN(qploidy_standarize_result = qploidy_standardization,
sample_id ="J432",
cn_grid = c(2,3,4,5,6))
# View results
res## hmm_CN result
## Copy-number grid: 2, 3, 4, 6
## Expected ploidy: 4
## Minimum SNPs per window: 24
## Initial state probabilities (pi0): 0, 0, 1, 0
## Estimated z means per CN: 0.206, 0.228, 0.254, 0.297
## Estimated z mean: 0.2543716
## Estimated z sigma: 0.16
## BAF Emission distribution: gaussian
## Final log-likelihood: 0
## Range of CN_call (window-level): 4 - 4Visualize with the different plots:
p <- plot_cn_track(hmm_CN = res,
qploidy_standarize_result = qploidy_standardization,
sample_id = "J432",
show_window_lines = FALSE,
z_by_mean = TRUE)
# All panels together
p$arranged
See the results in a data.frame format:
## $cn_grid
## [1] 2 3 4 6
##
## $distribution
## [1] "gaussian"
##
## $mu
## 2 3 4 6
## 0.2055931 0.2281309 0.2543715 0.2965014
##
## $sigma
## [1] 0.1596447
##
## $A
## [,1] [,2] [,3] [,4]
## [1,] 3.587895e-03 4.083188e-04 0.9591808 3.682300e-02
## [2,] 1.291517e-06 1.763386e-01 0.8142277 9.432460e-03
## [3,] 1.745587e-09 1.517963e-07 0.9999971 2.751808e-06
## [4,] 7.563872e-07 9.434018e-05 0.9438000 5.610489e-02
##
## $pi0
## [1] 2.553166e-11 8.902838e-10 9.999999e-01 1.228514e-0710 HMM Copy Number Estimation (Multiple Samples)
Run HMM for all samples in parallel. The model will be selected
internally using select_best_baf_model; you just need to
define the grids (bw, add_uniform_grid, and
uniform_weight_grid) according to what you observed when
running it on a single sample.
Note: Running HMM for many samples can be
computationally intensive. Adjust n_cores based on your
resources and dataset size.
10.1 Applying to raw data
multi_esti_raw <- hmm_estimate_CN_multi(data = data,
geno.pos = geno.pos,
use_values = c("ratio","R"),
sample_ids = "all",
n_cores = 2,
cn_grid = c(2,3,4,5,6),
bw = c(0.02, 0.04,0.1))Visualize summary of CNV for multiple samples:
Again, not the best results for this dataset. Standardization is required.
10.2 Applying to standardized data
multi_esti <- hmm_estimate_CN_multi(qploidy_standarize_result = qploidy_standardization,
sample_ids = "all",
n_cores = 2,
cn_grid = c(2,3,4,5,6),
bw = c(0.02, 0.04,0.1))Visualize summary of CNV for multiple samples:
Not that sample 85-103 pointed in section Chromosome-level screening based on z-score
having a higher depth compared to the rest of the samples in chromosome
6 suggesting a possible aneuploidy now is presented as four copies. The
combination of mixed models and hidden markov model eliminated false
variations of the expected copy number. Applying the approach to the
standardized data also eliminated several false CNV identified when
running using the raw data. See comparison below between raw and
standardized data for sample 85-113.
# Plots for a specific sample within the multi-sample object
plot_karyotype(multi_esti_raw, sample_name = "85-113",
nrow = 2, color_by = "w_baf")
p_stand <- plot_cn_track(hmm_CN = multi_esti_raw,
qploidy_standarize_result = qploidy_standardization,
sample_id = "85-113",
show_window_lines = TRUE)
p_stand$arranged
p_stand <- plot_cn_track(hmm_CN = multi_esti,
data = data,
sample_id = "85-113",
show_window_lines = TRUE)
p_stand$arranged
On the other hand, sample 85-126 highlight a false
negative from raw data analysis:
# For a specific sample
plot_karyotype(multi_esti_raw, sample_name = "85-126", nrow = 2, color_by = "w_baf")
p_stand <- plot_cn_track(hmm_CN = multi_esti_raw,
qploidy_standarize_result = qploidy_standardization,
sample_id = "85-126",
show_window_lines = TRUE)
p_stand$arranged
p_stand <- plot_cn_track(hmm_CN = multi_esti,
data = data,
sample_id = "85-126",
show_window_lines = TRUE)
p_stand$arranged
11 Re-standardization (Optional)
If you want to improve BAF resolution (e.g., after identifying a subset of samples with known ploidy via HMM), re-standardize:
re_qploidy_standardization <- re_standardize(
data = data,
geno.pos = geno.pos,
hmm_CN_multi = multi_esti,
use_estimated_dosages = TRUE,
ploidy.standardization = 4,
threshold.n.clusters = 4,
threshold.missing.geno = 0.9,
min.depth = 30,
max.depth = 600,
n.cores = 2,
out_filename = "re_standardization.tsv.gz",
verbose = TRUE
)
re_qploidy_standardization## This is an object of class 'qploidy_standardization'
## --------------------------------------------------------------------
## Parameters
##
## 1 Standardization type: intensities
## 2 Ploidy: 4
## 3 Min # of heterozygous classes (clusters) present: 4
## 4 Max proportion of missing genotype by marker: 0.9
## 5 Min genotype probability: 0.5
## 6 Min read depth (min.depth): 30
## 7 Max read depth (max.depth): 600
## --------------------------------------------------------------------
## Filters
##
## 1 # markers at raw data: 2766 (100%)
## 2 % datapoints filtered by low probability: - (0.2 %)
## 3 # datapoints filtered by depth (min/max): 84309 -
## 4 # markers filtered by missing data: 1883 (68.08 %)
## 5 # markers filtered by min number of clusters: 2 (0.07 %)
## 6 # markers filtered by lack of genomic information: 0 (0 %)
## 7 # markers remaining with estimated BAF: 779 (28.16 %)
## 8 # markers remaining with estimated z-score: 2664 (96.31 %)If argument use_estimated_dosages is set to TRUE,
Qploidy2 estimates the dosage based on its own algorithm. If FALSE, it
uses the dosages contained in the initial genos object.
# Compare before and after
plot_qploidy_standardization(x = qploidy_standardization,
sample = "85-295",
type = c("Ratio_hist_overall", "BAF_hist_overall"))
plot_qploidy_standardization(x = re_qploidy_standardization,
sample = "85-295",
type = c("Ratio_hist_overall", "BAF_hist_overall"))
12 Re-run the HMM with the new standardization
re_multi_esti <- hmm_estimate_CN_multi(
qploidy_standarize_result = re_qploidy_standardization,
sample_ids = "all",
n_cores = 2,
cn_grid = c(2,3,4,5,6),
bw = c(0.02, 0.04,0.1))
# CNV profiles before
compare_cn_track(hmm_CN = multi_esti,
samples_to_plot = samples[1:50],
facet_ncol = 8)
# CNV profiles after
compare_cn_track(hmm_CN = re_multi_esti,
samples_to_plot = samples[1:50],
facet_ncol = 8)
Among the first 50 samples, sample 85-117 differ between
previous standardization and re-standardization results.
# Same plots can be generated
plot_karyotype(re_multi_esti, sample_name = "85-2", nrow = 2, color_by = "post_max_CN")
In this case, there are a small fragment with lower copy number (CN
3) on part of chromosome 7. And a higher CN in chromosome 8. If you
prefer to remove the dashed line showing the missing parts. Just set the
add_ghost argument to FALSE:
# Same plots can be generated
plot_karyotype(re_multi_esti, sample_name = "85-2", nrow = 2, color_by = "post_max_CN", add_ghost = FALSE)
p_re <- plot_cn_track(hmm_CN = re_multi_esti,
qploidy_standarize_result = re_qploidy_standardization,
sample_id = "85-2",
show_window_lines = TRUE)
plot_karyotype(multi_esti, sample_name = "85-2", nrow = 2, color_by = "post_max_CN")
p <- plot_cn_track(hmm_CN = multi_esti,
qploidy_standarize_result = qploidy_standardization,
sample_id = "85-2",
show_window_lines = TRUE)
p$arranged #before
For this dataset, there are a few changes for some samples like the one showed above. The predominant yellow color mean that BAF distributions likelihoods were not included in the HMM due to high dispersion. If this instability still leaves doubts about the CNV. You can further explore tweaking the other parameters, for example, increasing the number of markers in the windows to better fit the mixed models on the BAF distributions:
res_tweak <- hmm_estimate_CN(qploidy_standarize_result = re_qploidy_standardization,
sample_id = "85-2",
min_snps_per_window = 60,
cn_grid = c(2,3,4,5,6))
plot_karyotype(multi_esti, sample_name = "85-2", nrow = 2, color_by = "post_max_CN")
p <- plot_cn_track(hmm_CN = res_tweak,
qploidy_standarize_result = re_qploidy_standardization,
sample_id = "85-2",
show_window_lines = TRUE)
p$arranged
Explore more parameters consulting ?hmm_estimate_CN.
13 Refine results for a single sample
Once you get to a conclusion with results better explain the patterns, you can update your multi-sample HMM object by using the functions:
If you decide that the results for this sample is inconclusive, you
can set the rm_sample parameter to TRUE to
remove it from your final results.
14 Plotting all samples in batches (Optional)
Generate CNV plots for all samples in batches:
n_plots_by_page <- c(seq(1, length(samples), by = 40), length(samples))
for(i in 1:(length(n_plots_by_page) -1)){
comp_p <- compare_cn_track(hmm_CN = final_res,
samples_to_plot = samples[n_plots_by_page[[i]]:n_plots_by_page[i+1]],
facet_ncol = 12)
ggsave(comp_p, filename = paste0("Samples_", n_plots_by_page[i], "_", n_plots_by_page[i +1], ".png"), height = 10, width = 11)
}15 Interactive curation of results
WARNING: Qploidy2 HMM results accuracy depends on
many factors that sometimes the call likelihoods are not capable to
detect. We strongly suggest to review the diagnostic plots for
non-expected results. The GenoBrew tool described below is
a great resource to make the graphics evaluation and results curation
process quick.
15.1 Export HMM results
The resulted objects from function update_hmm_CN_multi,
hmm_estimate_CN, and hmm_estimate_CN_multi are
all classified as
## [1] "hmm_CN"## [1] "hmm_CN"## [1] "hmm_CN"## [1] "hmm_CN"Therefore, all of them can be exported by the function:
You will see three generated files:
- “[prefix]_by_marker.csv.gz”
- “[prefix]_by_window.csv.gz”
- “[prefix]_params.rds”
You can read them back to R using the read_hmm_CN
function:
15.2 Deploy
GenoBrew
The resulted files from write_hmm_CN can be uploaded on
GenoBrew Shiny interface. The tab CNV profiles
allow you to interactively inspect the BAF, z-score and CN calls plots,
re-run the HMM adjusting parameters for specific samples and update the
results.
Visit GenoBrew GitHub repository and
tutorial
for more information.
16 Call dosages
Warning: Beta version - Dosage call requires more testing and review.
Qploidy2 uses the BAF distributions to call dosages. See an example:
# res was generated by hmm_estimate_CN above, it contains only one sample
sample_name <- "85-17"
baf <- final_res$by_marker$baf[final_res$by_marker$SampleName == sample_name]
selected_model <- select_best_baf_model(baf_vec = baf,
sample = sample_name,
cn_grid= c(2,3,4),
bw = c(0.02, 0.04,0.1),
plot = TRUE)
dosages_one_sample <- call_BAF_dosages(baf, selected_model = selected_model, plot=TRUE)
head(dosages_one_sample$data)
Call dosages for all samples considering the copy number estimated:
Once you run the dosage call, you can color the karyotype by the dosage or dosage likelihoods:
17 Export results
17.1 Summarized by sample, chromossome, window or marker
# By chromosome
summ_hmm_chr <- summarize_cn_mode(final_res, level = "chromosome")
head(summ_hmm_chr)You can use any of your favorite R table export functions to generate a file from those data.frames:
17.2 karyotype notation
Function karyotype_notation summarises per-marker
copy-number calls into a compact karyotype table: - Ploidy: reports
overall ploidy - ChrID: chromosome ID - ChrCN: reports most common CN
for the chromosome - Segment: reports position of segment with CN
different than estimated overall ploidy (ISCN-ispired notation:
g.start_end×CN). If not CNV segments were estimated, column reports
NA
17.3 VCF format
You can read the VCF with vcfR and extract fields:
## Scanning file to determine attributes.
## File attributes:
## meta lines: 12
## header_line: 13
## variant count: 2664
## column count: 194
## Meta line 12 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 2664
## Character matrix gt cols: 194
## skip: 0
## nrows: 2664
## row_num: 0
## Processed variant 1000Processed variant 2000Processed variant: 2664
## All variants processed## FORMAT
## [1,] "GT:CN:AD:BAF:Z:PMC:PMD:PLC:PLD:WBF"
## [2,] "GT:CN:AD:BAF:Z:PMC:PMD:PLC:PLD:WBF"
## [3,] "GT:CN:AD:BAF:Z:PMC:PMD:PLC:PLD:WBF"
## [4,] "GT:CN:AD:BAF:Z:PMC:PMD:PLC:PLD:WBF"
## [5,] "GT:CN:AD:BAF:Z:PMC:PMD:PLC:PLD:WBF"
## 85-1
## [1,] "./.:4:1,339:.: 1.38:1.000:.:.,110.0,0.0,87.6,79.0:.,.,.,.,.,.,.:0.5000"
## [2,] "./.:4:128,0:.:-0.31:1.000:.:.,110.0,0.0,87.6,79.0:.,.,.,.,.,.,.:0.5000"
## [3,] "0/1/1/1:4:24,90:0.78: 1.51:1.000:1.000:.,110.0,0.0,87.6,79.0:814.5,373.4,101.9,0.0,67.8,.,.:0.5000"
## [4,] "./.:4:158,96:.: 1.26:1.000:.:.,110.0,0.0,87.6,79.0:.,.,.,.,.,.,.:0.5000"
## [5,] "./.:4:36,1:.:-0.96:1.000:.:.,110.0,0.0,87.6,79.0:.,.,.,.,.,.,.:0.5000"
## 85-100
## [1,] "./.:4:0,239:.:-0.49:1.000:.:.,.,-0.0,461.3,.:.,.,.,.,.,.,.:0.4170"
## [2,] "./.:4:104,0:.:-0.73:1.000:.:.,.,-0.0,461.3,.:.,.,.,.,.,.,.:0.4170"
## [3,] "0/0/1/1:4:28,37:0.55:-0.83:1.000:1.000:.,.,-0.0,461.3,.:403.1,116.7,0.0,52.9,275.5,.,.:0.4170"
## [4,] "./.:4:211,34:.: 1.08:1.000:.:.,.,-0.0,461.3,.:.,.,.,.,.,.,.:0.4170"
## [5,] "./.:4:47,0:.:-0.53:1.000:.:.,.,-0.0,461.3,.:.,.,.,.,.,.,.:0.4170"
## 85-102
## [1,] "./.:4:55,204:.:-0.12:1.000:.:186.8,104.1,0.0,.,58.0:.,.,.,.,.,.,.:0.5000"
## [2,] "./.:4:104,0:.:-0.73:1.000:.:186.8,104.1,0.0,.,58.0:.,.,.,.,.,.,.:0.5000"
## [3,] "0/0/1/1:4:44,49:0.50: 0.51:1.000:1.000:186.8,104.1,0.0,.,58.0:343.8,87.1,0.0,82.6,334.8,.,.:0.5000"
## [4,] "./.:4:228,0:.: 0.74:1.000:.:186.8,104.1,0.0,.,58.0:.,.,.,.,.,.,.:0.5000"
## [5,] "./.:4:33,1:.:-1.09:1.000:.:186.8,104.1,0.0,.,58.0:.,.,.,.,.,.,.:0.5000"
## 85-103
## [1,] "./.:4:53,211:.:-0.03:1.000:.:138.1,111.7,0.0,.,41.4:.,.,.,.,.,.,.:0.3495"
## [2,] "./:.:.,.:.:.:.:.:.:."
## [3,] "0/0/1/1:4:33,56:0.61: 0.31:1.000:0.925:138.1,111.7,0.0,.,41.4:487.4,159.1,0.3,11.3,191.8,.,.:0.3495"
## [4,] "./.:4:132,48:.:-0.22:1.000:.:138.1,111.7,0.0,.,41.4:.,.,.,.,.,.,.:0.3495"
## [5,] "./.:4:33,0:.:-1.14:1.000:.:138.1,111.7,0.0,.,41.4:.,.,.,.,.,.,.:0.3495"## 85-1 85-100 85-102 85-103 85-104
## chr1.1_000194324 "4" "4" "4" "4" "4"
## chr1.1_000532584 "4" "4" "4" NA "4"
## chr1.1_000735393 "4" "4" "4" "4" "4"
## chr1.1_000837330 "4" "4" "4" "4" "4"
## chr1.1_000915014 "4" "4" "4" "4" "4"The VCF format dosages are adapted to the CN estimated by Qploidy2. For example, if the CN is 3, the dosages can be 0,1,2,3 and are represented as 0/0/0, 0/0/1, 0/1/1, 1/1/1. If the CN is 4, the dosages can be 0,1,2,3,4 and are represented as 0/0/0/0, 0/0/0/1, 0/0/1/1, 0/1/1/1, 1/1/1/1.
## 85-157 85-158 85-168 85-169 85-17 85-171
## chr7.1_078445464 "4" "4" "4" "4" "4" "4"
## chr7.1_078682539 NA "4" "4" "4" "4" "4"
## chr7.1_079844136 "4" "4" "4" "4" "4" "4"
## chr7.1_080085927 NA "4" "4" "4" "4" "4"
## chr7.1_080411685 "4" "4" "4" "4" "4" "4"## 85-157 85-158 85-168 85-169 85-17 85-171
## chr7.1_078445464 NA NA NA NA NA NA
## chr7.1_078682539 "./" "0/0/0/1" "0/0/0/1" "0/0/1/1" "0/0/0/1" "0/1/1/1"
## chr7.1_079844136 NA NA NA NA NA NA
## chr7.1_080085927 "./" NA NA NA NA NA
## chr7.1_080411685 "0/0/0/1" "0/0/1/1" "0/0/1/1" "0/0/1/1" "0/0/0/1" "0/1/1/1"The VCF contains the likelihood of each dosage call in the
PMD field. The PMD is the probability of the most likely
dosage call.
## 85-1 85-100 85-102 85-103 85-104
## chr1.1_000194324 NA NA NA NA NA
## chr1.1_000532584 NA NA NA NA NA
## chr1.1_000735393 "1.000" "1.000" "1.000" "0.925" "1.000"
## chr1.1_000837330 NA NA NA NA NA
## chr1.1_000915014 NA NA NA NA NA## 85-157 85-158 85-168 85-169 85-17 85-171
## chr7.1_078445464 NA NA NA NA NA NA
## chr7.1_078682539 NA "1.000" "1.000" "1.000" "0.880" "1.000"
## chr7.1_079844136 NA NA NA NA NA NA
## chr7.1_080085927 NA NA NA NA NA NA
## chr7.1_080411685 "0.914" "1.000" "1.000" "1.000" "1.000" "1.000"18 How to cite
18.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
18.0.2 Qploidy2 - HMM and grid approach BAF model selection
Manuscript in preparation. Please contact the author for more information.
18.0.3 nQuack - EM approach 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. doi:10.1002/aps3.11606
19 Acknowledgments
19.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).
19.0.2 Qploidy versions > 1.0.0 and Qploidy2
Supported by Breeding Insight.