This vignette describes how to prepare the input files for use with
Qploidy2, an R package designed for ploidy, aneuploidy and
large-scale copy number variation estimation using genotyping platform
data.
Before describing the inputs, it is important to understand the
conditions under which the Qploidy2 methodology is
effective, as well as its limitations. This understanding will help you
determine whether Qploidy2 is the right tool for your
dataset and research objectives.
Qploidy2 methodology work?The Qploidy2 approach 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 (make sure you have sufficient coverage for accurate
copy number estimation and you filter your data appropriately). - All
DNA samples were prepared following the same library preparation
protocol.
- You know the ploidy of at least a subset of 60
samples or you know the most common ploidy in
the dataset.
- Your dataset includes heterozygous samples.
Qploidy2 methodology NOT work?The methodology will not be effective under the following
circumstances:
- Your marker data comes from RADseq or
GBS (Genotyping-by-Sequencing) platforms (sparse and
inconsistent marker coverage across samples makes the signal
standardization unreliable).
- You intend to combine datasets from different sequencing
batches.
- 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 Qploidy2. Combining
all 480 samples into a single analysis will lead to incorrect results.
See section “Check for batch effect” in the tutorial “Aneuploid
Detection in Alfalfa F1 Population” for more details. - You do
not have a subset of samples with known ploidy or lack
a predominant ploidy in your dataset.
- Qploidy2 may still work depending on the quality of the data. Give it
a try, but it is not guaranteed. - Your samples consist exclusively of
inbred lines (homozygous individuals), as dosage
calling relies on the segregation of allele intensities, which does not
occur in fully homozygous lines. - Your samples consist of
first-generation, colchicine-treated individuals in which genome
duplication was directly induced, as the complete duplication of the
genome doesn’t generate all possible dosages of a polyploid. Use
progenies of these individuals instead, as they will segregate for all
possible dosages.
To install the Qploidy2 package, you can use the
following commands:
Qploidy2 requires three main data.frames
that we call:
data): contains the allele intensities or read
counts for each marker and sample.genos): contains the genotype dosages for each
marker and sample of known ploidy or for all samples if overall ploidy
is being used.genos.pos): contains the genomic
positions for each marker.If you have standard format files, Qploidy2 provides
functions to read them and convert them into these three required
data.frames.
In this vignette, we demonstrate file format conversion using simulated data which mimics standard formats from: - VCF (Variant Call Format) file - Axiom array summary file - Illumina array summary file
Data
data.frame formatThe required FORMAT in a VCF for Qploidy2 is the AD (Allele Depth) field, which contains the read counts for the reference and alternative alleles. Let’s generate our simulated VCF file. This file contains 500 markers and 60 samples, with a mix of tetraploid, diploid, and triploid samples. The tetraploid samples are used as the reference for standardization because they are the most common ploidy in the dataset.
library(Qploidy2)
# Simulate a VCF file with 500 markers and 60 samples
simulate_vcf(
seed = 1234, file_path = "vcf_example_simulated.vcf",
n_tetraploid = 35, n_diploid = 5, n_triploid = 10,
n_markers = 500
)The resulting data data.frame contains the following
columns: MarkerName and SampleName
(identifiers), X and Y (allele signal
intensities for the reference and alternative alleles, respectively),
R (total signal intensity, i.e. X + Y), and
ratio (the B-allele frequency,
i.e. Y / R).
In case your samples were genotyped using the Axiom array platform,
you should have a summary file or .CEL files.
If you have .CEL files, you can convert them into a
summary file using the apt-probeset-genotype
plugin available in the Analysis Power Tools (APT) platform from Thermo
Fisher Scientific (Waltham, MA). A summary file is composed
of a header that can have many lines and a body that looks like
this:
In Qploidy2, use the read_axiom function to
read the summary file. It skips the header and reads the
body of the file. The function will return a data frame.
MarkerName SampleName X Y R ratio
1 AX-86752740 Lemon_Fiz 977.7548 793.4656 1771.220 0.4479768
2 AX-86752740 High_Voltage 1292.2205 582.7886 1875.009 0.3108191
3 AX-86752740 Lemon_Fiz.1 919.8055 870.9619 1790.767 0.4863624
4 AX-86752740 High_Voltage.1 1368.1266 645.9481 2014.075 0.3207170
5 AX-86752740 Brite_Eyes 238.5416 1242.7593 1481.301 0.8389648
6 AX-86752740 Morden_Blush.1 302.1575 1542.4037 1844.561 0.8361900
In case you have a Illumina array summary file, Qploidy2
expects the following format:
[Header]
GSGT Version 2.0.5
Processing Date 4/17/2023 9:31 AM
Content specieX.bpm
Num SNPs 26251
Total SNPs 26251
Num Samples 100
Total Samples 100
[Data]
SNP Name Sample ID GC Score Theta X Y X Raw Y Raw Log R Ratio
mk-1 SAMP-2 0.6814 0.247 0.023 0.009 504 100 0.5211668
mk-10 SAMP-2 0.9204 0.000 0.006 0.000 346 27 -1.021583
mk-100 SAMP-2 0.1695 0.367 0.030 0.020 576 158 1.177404
mk-101 SAMP-2 0.8655 0.222 0.026 0.010 538 102 0.3459635
mk-103 SAMP-2 0.8932 0.000 0.019 0.000 465 18 0.7578704
mk-104 SAMP-2 0.5146 0.175 0.041 0.012 675 114 1.016254
mk-105 SAMP-2 0.8060 0.450 0.018 0.015 458 132 0.4305636
mk-106 SAMP-2 0.0000 0.386 0.002 0.002 315 57 NaN
mk-107 SAMP-2 0.4804 0.495 0.026 0.025 534 189 1.0151
Verify if your file looks like this one. If confirmed, you can use
the read_illumina_array function to read the data. If you
include more than one file, the function will merge them into a single
data frame.
data <- read_illumina_array(
"data/illumina/Set1.txt", # This is a example code, data not available
"data/illumina/Set2.txt",
"data/illumina/Set3.txt"
)
head(data)MarkerName SampleName X Y R ratio
1 CT1 S1-2_1 0.023 0.009 0.032 0.2812500
2 CT10 S1-2_1 0.006 0.000 0.006 0.0000000
3 CT100 S1-2_1 0.030 0.020 0.050 0.4000000
4 CT101 S1-2_1 0.026 0.010 0.036 0.2777778
5 CT103 S1-2_1 0.019 0.000 0.019 0.0000000
6 CT104 S1-2_1 0.041 0.012 0.053 0.2264151
Genos
data.frame formatThe samples contained in the genos data frame will
depend on the chosen path below.
Qploidy2 applies a standardization method to allele
intensities or read counts to enable the comparison of copy numbers
across the genome of a sample. To achieve this, it requires a set of
reference samples. There are two ways to select these reference
samples:
genos data.frame should only contain the genotype
dosages for this subset of samples. You will use this subset as the
reference in the Qploidy2 standardization step.genos data.frame should contain the
genotype dosages for all samples. You can initially run
Qploidy2 using all samples as the reference. From this
first run, you can identify samples with the common ploidy, then select
this subset to use as the reference in a second run of
Qploidy2. Check section re-standardization in
the tutorial “Ploidy Estimated in a Mixed Population” for more
details.Below are the steps to proceed for each case, tagged as (1) and (2):
For this method, you will need to separate your subset of samples with known ploidy:
In our simulated example, the sample names tell us which ploidy they are. Let’s create a vector with the tetraploid samples for us to use later:
## [1] "Tetraploid1" "Tetraploid2" "Tetraploid3" "Tetraploid4" "Tetraploid5"
## [6] "Tetraploid6" "Tetraploid7" "Tetraploid8" "Tetraploid9" "Tetraploid10"
## [11] "Tetraploid11" "Tetraploid12" "Tetraploid13" "Tetraploid14" "Tetraploid15"
## [16] "Tetraploid16" "Tetraploid17" "Tetraploid18" "Tetraploid19" "Tetraploid20"
## [21] "Tetraploid21" "Tetraploid22" "Tetraploid23" "Tetraploid24" "Tetraploid25"
## [26] "Tetraploid26" "Tetraploid27" "Tetraploid28" "Tetraploid29" "Tetraploid30"
## [31] "Tetraploid31" "Tetraploid32" "Tetraploid33" "Tetraploid34" "Tetraploid35"
## [36] "Diploid1" "Diploid2" "Diploid3" "Diploid4" "Diploid5"
## [41] "Triploid1" "Triploid2" "Triploid3" "Triploid4" "Triploid5"
## [46] "Triploid6" "Triploid7" "Triploid8" "Triploid9" "Triploid10"
## [1] "Tetraploid1" "Tetraploid2" "Tetraploid3" "Tetraploid4" "Tetraploid5"
## [6] "Tetraploid6" "Tetraploid7" "Tetraploid8" "Tetraploid9" "Tetraploid10"
## [11] "Tetraploid11" "Tetraploid12" "Tetraploid13" "Tetraploid14" "Tetraploid15"
## [16] "Tetraploid16" "Tetraploid17" "Tetraploid18" "Tetraploid19" "Tetraploid20"
## [21] "Tetraploid21" "Tetraploid22" "Tetraploid23" "Tetraploid24" "Tetraploid25"
## [26] "Tetraploid26" "Tetraploid27" "Tetraploid28" "Tetraploid29" "Tetraploid30"
## [31] "Tetraploid31" "Tetraploid32" "Tetraploid33" "Tetraploid34" "Tetraploid35"
For this method, you will use the entire dataset as the reference
during the first round of the Qploidy2 run:
To proceed, you need to determine the dosage for all reference samples. If you haven’t done this yet, you can use any suitable dosage caller.
fitPoly
package.updog or polyRAD
packages.These packages determine dosages by analyzing the distribution of allele intensities or read counts across all samples for each marker, using an informed ploidy model. You will provide either the most common ploidy in your population or the known ploidy of your selected subset as the informed ploidy.
Warning:
Depending on the number of samples and markers, this step may take a
significant amount of time to complete. It is highly recommended to run
it on a high-performance computing system where you can utilize multiple
cores and avoid issues with excessive RAM usage.
## [1] 50000 4
If your VCF lacks the GT information for your reference samples,
perform dosage calling using the AD field as described in the next
section. After that, you can create the genos data frame
with the dosage information and the genos.pos data frame
with the genomic positions of the markers.
ID column of
the fixed portion, the MarkerName column will contain
NAs. To resolve this, you can use the following code to
concatenate the chromosome and position as marker names:library(tidyr)
# Prepare inputs for updog
## Approach (1)
# tobe_genotyped <- data[which(data$SampleName %in% tetraploid_samples),]
## Approach (2)
tobe_genotyped <- data
ref <- pivot_wider(tobe_genotyped[, 1:3], names_from = SampleName, values_from = X)
ref <- as.matrix(ref)
rownames_ref <- ref[, 1]
ref <- ref[, -1]
ref <- apply(ref, 2, as.numeric)
rownames(ref) <- rownames_ref
size <- pivot_wider(tobe_genotyped[, c(1, 2, 5)], names_from = SampleName, values_from = R)
size <- as.matrix(size)
rownames_size <- size[, 1]
size <- size[, -1]
size <- apply(size, 2, as.numeric)
rownames(size) <- rownames_size
size[1:10, 1:10]## Tetraploid1 Tetraploid2 Tetraploid3 Tetraploid4 Tetraploid5
## chr1_mk1 215 202 188 212 202
## chr1_mk2 184 194 208 205 213
## chr1_mk3 182 200 193 196 214
## chr1_mk4 189 202 202 203 203
## chr1_mk5 194 209 200 199 207
## chr1_mk6 186 206 178 195 205
## chr1_mk7 216 200 206 191 193
## chr1_mk8 196 215 201 209 199
## chr1_mk9 198 216 211 197 205
## chr1_mk10 204 195 205 187 205
## Tetraploid6 Tetraploid7 Tetraploid8 Tetraploid9 Tetraploid10
## chr1_mk1 195 195 205 210 196
## chr1_mk2 200 195 209 200 195
## chr1_mk3 205 199 205 194 220
## chr1_mk4 215 204 218 203 188
## chr1_mk5 213 208 191 191 208
## chr1_mk6 203 191 200 197 200
## chr1_mk7 214 198 204 203 198
## chr1_mk8 212 200 195 197 205
## chr1_mk9 204 198 200 189 208
## chr1_mk10 212 200 201 199 200
## Tetraploid1 Tetraploid2 Tetraploid3 Tetraploid4 Tetraploid5
## chr1_mk1 62 98 97 156 158
## chr1_mk2 93 92 56 105 158
## chr1_mk3 143 110 146 52 98
## chr1_mk4 43 102 149 158 95
## chr1_mk5 101 107 200 153 151
## chr1_mk6 148 102 40 49 56
## chr1_mk7 112 0 150 141 94
## chr1_mk8 50 56 54 155 101
## chr1_mk9 98 104 58 101 108
## chr1_mk10 145 148 151 145 106
## Tetraploid6 Tetraploid7 Tetraploid8 Tetraploid9 Tetraploid10
## chr1_mk1 147 149 154 99 145
## chr1_mk2 103 50 104 0 45
## chr1_mk3 150 53 100 105 109
## chr1_mk4 55 50 105 165 144
## chr1_mk5 56 154 154 98 101
## chr1_mk6 150 94 0 147 0
## chr1_mk7 70 99 50 98 154
## chr1_mk8 151 151 140 46 114
## chr1_mk9 54 44 101 86 156
## chr1_mk10 96 200 152 150 200
library(updog)
multidog_obj <- multidog(
refmat = ref,
sizemat = size,
model = "norm", # It depends on your population structure
ploidy = 4, # Most common ploidy in your population
nc = 2
) # Change the parameters accordingly!!## | *.#,%
## ||| *******/
## ||||||| (**..#**. */ **/
## ||||||||| */****************************/*%
## ||| &****..,*.************************/
## ||| (....,,,*,...****%********/(******
## ||| ,,****%////,,,,./.****/
## ||| /**// .*///....
## ||| .*/*/%# .,/ .,
## ||| , **/ #% .* ..
## ||| ,,,*
##
## Working on it...done!
genos <- data.frame(
MarkerName = multidog_obj$inddf$snp,
SampleName = multidog_obj$inddf$ind,
geno = multidog_obj$inddf$geno,
prob = multidog_obj$inddf$maxpostprob
)
head(genos)Notice that, if you are following approach (2), diploid and triploid sample genotypes will be incorrectly called as tetraploids. However, because they are not the most common samples in the dataset, they should not interfere significantly with the dosage calling of the tetraploids. You can go back to this step after a first round of ploidy evaluation and call the dosages of the tetraploids only.
Here we show an option of genotype calling for array data using fitpoly software:
library(fitPoly)
saveMarkerModels(
ploidy = 4, # Most common ploidy among Texas Garden Roses collection
data = data_reference, # output of the previous function
filePrefix = "fitpoly_output/texas_roses", # Define the path and prefix for the result files
ncores = 2
) # Change it accordingly to your system!!!!
library(vroom) # Package to speed up reading files
fitpoly_scores <- vroom("fitpoly_output/texas_roses_scores.dat")
genos <- data.frame(
MarkerName = fitpoly_scores$MarkerName,
SampleName = fitpoly_scores$SampleName,
geno = fitpoly_scores$maxgeno,
prob = fitpoly_scores$maxP
)
head(genos) MarkerName SampleName geno prob
1 AX-86752740 1000 Wishes 2 0.7831438
2 AX-86752740 10004-N008 2 0.9892338
3 AX-86752740 10037_N046 2 0.9844171
4 AX-86752740 10038-N001 2 0.9916813
5 AX-86752740 10043_N019 2 0.9871057
6 AX-86752740 10043_N049 2 0.9955792
Unlike a VCF file, a summary file does not contain the
genomic positions of the markers. You will need to provide this
information in a separate file, as shown in the following example.
# File containing markers genomic positions
genos.pos_file <- read.table("geno.pos_roses_texas.txt", header = T)
head(genos.pos_file) SNP probes chr pos
1 Affx-86843634 AX-86752747 6 255852
2 Affx-86842821 AX-86752763 6 62032267
3 Affx-86839613 AX-86752769 5 9169310
4 Affx-86838724 AX-86752790 1 44259488
5 Affx-86840823 AX-86752809 7 12991931
6 Affx-86842443 AX-86752817 2 53861719
# Edit for Qploidy2 input format
genos.pos <- data.frame(
MarkerName = genos.pos_file$probes,
Chromosome = genos.pos_file$chr,
Position = genos.pos_file$pos
)
head(genos.pos) MarkerName Chromosome Position
1 AX-86752747 6 255852
2 AX-86752763 6 62032267
3 AX-86752769 5 9169310
4 AX-86752790 1 44259488
5 AX-86752809 7 12991931
6 AX-86752817 2 53861719
By the end of this tutorial you should have the three required data
frames (data, genos, and
genos.pos) ready to be used as input for
Qploidy2. You are now ready to run the analysis. Choose the
tutorial that best fit your purpose to go next:
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
Manuscript in preparation. Please contact the author for more information.
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).
Supported by Breeding Insight.