Alignment-free methods for polyploid genomes: quick and reliable genetic distance estimation

Polyploid genomes pose several inherent challenges to population genetic analyses. While alignment-based methods are fundamentally limited in their applicability to polyploids, alignment-free methods bypass most of these limits. We investigated the use of Mash, a k-mer analysis tool that uses the MinHash method to reduce complexity in large genomic datasets, for basic population genetic analyses of polyploid sequences. We measured the degree to which Mash correctly estimated pairwise genetic distance in simulated diploid and polyploid short-read sequences with various levels of missing data. Mash-based estimates of genetic distance were comparable to alignment-based estimates, and were less impacted by missing data. We also used Mash to analyze publicly available short-read data for three polyploid and one diploid species, then compared Mash results to published results. For both simulated and real data, Mash accurately estimated pairwise genetic differences for polyploids as well as diploids as much as 476 times faster than alignment-based methods, though we found that Mash genetic distance estimates could be biased by per-sample read depth. Mash may be a particularly useful addition to the toolkit of polyploid geneticists for rapid confirmation of alignment-based results and for basic population genetics in reference-free systems with poor quality DNA.


10
Abstract 11 Polyploid genomes pose several inherent challenges to population genetic analyses. While alignment-12 based methods are fundamentally limited in their applicability to polyploids, alignment-free methods 13 bypass most of these limits. We investigated the use of Mash, a k-mer analysis tool that uses the MinHash 14 method to reduce complexity in large genomic datasets, for basic population genetic analyses of polyploid 15 sequences. We measured the degree to which Mash correctly estimated pairwise genetic distance in 16 simulated diploid and polyploid short-read sequences with various levels of missing data. Mash-based 17 estimates of genetic distance were comparable to alignment-based estimates, and were less impacted by 18 missing data. We also used Mash to analyze publicly available short-read data for three polyploid and one 19 diploid species, then compared Mash results to published results. For both simulated and real data, Mash 20 accurately estimated pairwise genetic differences for polyploids as well as diploids as much as 476 times 21 faster than alignment-based methods, though we found that Mash genetic distance estimates could be 22 biased by per-sample read depth. Mash may be a particularly useful addition to the toolkit of polyploid 23 geneticists for rapid confirmation of alignment-based results and for basic population genetics in 24 reference-free systems with poor quality DNA .  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43 alleles A and B could have a genotype ABBBBBBB, which should have a MAF of 0.125. The genotype 70 AABBBBBB should have a MAF of 0.25. If a sequenced locus' MAF were 0.18, for instance, it would be 71 challenging to determine which is the true genotype. This challenge is compounded by the fact that fewer 72 high-quality genomes have been assembled for polyploids, particularly for higher ploidy levels. where commonly used statistics may apply to polyploid datasets, there are practical concerns as well. 81 Most bioinformatic tools require biallelic SNPs, so polyploid data analysts must produce custom analysis 82 scripts that are challenging for reviewers to evaluate, or constrain analyses to the few tools that 83 accommodate polyploids. 84

85
Solutions are being continuously developed to overcome the challenges of polyploid genome analysis. 86 Alternatives to F-statistics have been developed that are less impacted by ploidy. One of the most useful 87 is ϱ (Ronfort et al. 1998), which is equivalent to the average relatedness of individuals within populations 88 versus the average relatedness. ϱ is comparable between ploidy levels (Meirmans et al. 2018), but has not 89 been applied widely in empirical studies. For certain applications, unsupervised machine learning 90 methods have shown great promise in identifying population genetic parameters using methods that are 91 applicable to polyploids (Flagel et al. 2019). These methods are agnostic to many of the assumptions of 92 other population genetic models, but the requirement of extensive training data mean they will not be 93 applicable for all systems. 94 Currently, there are few resources that are both robust to the large degree of variation inherent to 96 polyploid genomes, and repeatable enough to provide useful benchmarks between population genetic 97 analyses. Since many of the challenges with polyploids come from the alignment and SNP-calling steps, 98 an ideal solution would completely bypass these steps. Many alignment-free sequence analysis techniques 99 have been developed in recent years, mostly using "words" or k-mers, subsamples of sequences with 100 length k. Analyses that use k-mers add flexibility to sequence analysis and can better make use of and so are a promising tool for population genetics in polyploid systems. There are many ways to 104 compare k-mer decompositions of sequence data, but one of the most promising is the use of the 105 MinHash algorithm. MinHash is a type of locality-sensitive hashing, an algorithmic technique to collapse 106 similar input items into short strings, or hashes. MinHash then uses a Jaccard similarity coefficient to 107 quickly assess how similar sets of hashes are. For sequence data, raw sequence reads for a sample can be 108 decomposed into k-mer tables, which are then hashed and can be quickly compared between individuals. 109 In practice, this process is many times faster than sequence alignment -researchers were able to cluster all In this paper we evaluate the performance of Mash to estimate pairwise sample relationships within and 116 across populations and species in polyploids. First, we simulate haploid and polyploid DNA sequences 117 with known phylogenetic relationships, and show that Mash accurately recovers cladograms for these 118 samples. Next, we remove reads from these simulated datasets to show the impact of missing data on 119 Mash estimation relative to alignment-based analysts, and find that Mash is more robust to missing data 120 than alignment-based estimation at estimating pairwise relationships. Finally, in several case studies, we use Mash to estimate kinship in reduced-representation sequencing reads for polyploid genera Capsella 122 and Reynoutria, as well as whole-genome data in polyploid Panicum and diploid Oryza. We first evaluated the performance of Mash using simulated haploid and polyploid sequence data. We 129 simulated phylogenetic trees using toytree (Eaton 2020) in Python 3.7 for 50 individuals using a tree 130 height of 1e6, then simulated SNP loci following these trees by generating sequences in ipcoal 131 (McKenzie & Eaton 2020). For haploid data, we simulated 1000 loci with 100bp each for 50 individuals, 132 diverging with an effective population size of 1e5, and a recombination rate of 1e-9. For polyploid data, 133 we simulated an allopolyploid species with two subgenomes by running two simulations based on the 134 same tree and initial parameters, but different random seeds, then combining reads into a joint fasta file. 135

136
To simulate missing data, we used custom scripts in R to randomly remove sequence data for 5-75% of 137 reads. We compared Mash and alignment-based analyses for each level of missing data. We first calculated the true distance between individuals as the proportion of nucleotide differences (π) in 142 the simulated sequence data. We then compared these true distances to distances inferred by alignment-143 based and alignment-free methods. 144

145
To provide the simplest comparison of methods, we ran a naive pipeline for both alignment-based and 146 Mash methods, using the default parameters throughout. We aligned reads using bwa and samtools, then For Mash, we first sketched the simulated reads using mash sketch, which uses minhashing to produce a 153 reduced-complexity representation of the read set. We then calculated Mash distances between 154 individuals using mash dist. We imported pairwise Mash distances into R and built trees again using 155 hclust. We calculated the true distance between individuals as the proportion of nucleotide differences 156 (π), and compared this to missing distance matrices using Mantel tests. We also compared the true tree to 157 each tree inferred from missing data using the comparePhylo function in ape (Paradis et al. 2004). we removed k-mers with counts under 2, used a k-mer length of 21, and increased the number of k-mers 166 per sketch to 1e6. We used the same analysis tools as for simulated data, and visualized sample 167 relationships using Principal Coordinates Analysis in R. We found that per-sample read count biased 168 Mash analyses, so we randomly removed reads from the raw sequence files to standardize inputs. 169

170
While these examples are useful for ensuring that Mash can at least match the performance of alignment-171 based methods, it is not clear in these reference-free studies how well either method estimates true 172 population relationships. We downloaded whole genome-sequencing data for 16 samples in Panicum distance between individuals using 11 million SNPs. We then simulated reduced-representation sequence 175 data by trimming whole-genome sequences to reads that contain a typical RAD (Restriction-Associated 176 Digest) restriction enzyme cut site (EcoRI), then dropped low-quality reads to standardize read count 177 across samples. Using the whole-genome SNP data as "true" population relationships, we then evaluated 178 how well Mash and alignment-based methods approximate these relationships. For sequence alignments, 179 we used bwa (Li 2013) to align to the switchgrass v5 reference genome, then called SNPs using the 180 Stacks ref_map.pl RAD-processing pipeline (Catchen et al. 2013). We called 783 SNPs from this 181 pipeline, and estimated the genetic distance between aligned samples again using the dist.genpop function 182 in adegenet. For Mash analyses, we calculated genetic distances from whole-genome sequences as well as 183 simulated RAD sequences. We used the same parameters as for the case studies: removing k-mers with 184 counts under 2, a k-mer length of 21, and increasing the number of k-mers per sketch to 1e6. Mash analyses and alignment-based analyses differed greatly when we introduced missing data (Figures 1  214   & 2). Alignment-based analyses were robust to lower levels of missing data, at 5% missing data, there 215 was a 97.2% concordance between missing and true distances. However, alignment-based distances 216 deviated more greatly from true distances at higher levels ( Figure 1). Indeed, Mash-based distance 217 estimates never showed less than 50% discordance from true distances, whereas alignment-based analyses 218 performed more poorly at levels of missing data >40%, showing a marked decrease (Figure 1. Mash analyses in samples with greatly divergent read counts, but we were able to adjust by randomly 246 trimming reads from raw sequence data. For Capsella, Mash showed five population groups, rather than 247 the three found by ADMIXTURE, but showed high similarity to reconstructed population distributions 248 despite using different methods ( Figure 3A). We compared compute time for Reynoutria between Mash 249 and alignment-based methods. To generate distance matrices from raw sequence reads for 25Gb of data 250 took 0.233 compute hours in Mash and 108 compute hours for alignment. We could not benchmark 251 Capsella alignment-based analyses, but Mash analyzed 44.3Gb of data in 0.410 compute hours. Since 252 sketching is easily multithreaded, most computers can do this in a fraction of the time. 253

254
To evaluate the performance of Mash against "known" population relationships in real sequence data, we 255 compared the performance of Mash distances against alignment-based methods for simulated RAD data. 256 We found that Mash effectively reproduced population relationships in Panicum virgatum for both whole-257 genome and simulated RAD data. Comparatively, the alignment-based method performed poorly, 258 showing much less difference between known subpopulations than Mash of RAD sequences ( Figure 3C).   This study showed that Mash performs well at assessing genetic relationships between polyploid 310 individuals in a fast, naive, and repeatable manner for both simulated and real data. For simulated data, 311 we showed that Mash outperforms alignment-based methods in distance estimation when missing data is 312 40% or greater (a reasonable amount of missing data, even in whole genome sequencing projects) and is 313 orders of magnitude faster. For real data, we have shown that Mash can reproduce population structure 314 estimates in both polyploids and diploids, though per-sample read count should be standardized before 315 analysis. Mash should therefore be considered a useful tool for researchers investigating population 316 relationships in their samples, particularly when missing data is high, or when rapid population 317 relationship assessment is necessary. Mash for the simulated data in this study, but neither sourmash nor kmer-db showed advantage over 341 Mash. 342 343 II. Potential applications for k-mer hashing based methods 344

345
As the quantity of available sequence data continues to increase exponentially, new tools that can rapidly 346 process these data are needed (Barone et al. 2017). Locality-sensitive hashing has been used for data-347 intensive jobs like clustering metagenomes (Rasheed et al. 2013) and genome assembly (Berlin et al. 348 2015). The use of k-mer methods in molecular ecology has been comparatively limited, perhaps since 349 historically this field has used less data-intensive approaches. However, now that genomics are more 350 feasible for many molecular ecology researchers, k-mer hashing may solve some of the scaling issues that 351 arise. In particular, three main directions will be facilitated by Mash and similar k-mer hashing methods: 352 population genetics for polyploids, reference-free genome-wide association studies (GWAS), and 353 genomics of ancient and otherwise degraded samples. Although we do not describe their use here, the accuracy of k-mer based methods suggests that they may 377 also be useful for generating Φst and other genome-wide estimates between populations, particularly in 378 polyploid systems. We suggest that future efforts be dedicated to assessing the accuracy of k-mers for this 379 task. Traditionally, population genetics has used a small sample of loci from a genome to infer the history of 391 the organism. A major issue with this method has always been that genetic changes that accompany 392 demographic changes are missed when using only relatively few neutral loci, particularly when these variants, k-mers will also show this variation even in short-read data through presence and absence 397 (Voichek & Weigel 2020). As molecular ecology studies increasingly use whole-genome coverage to 398 assess population-level relationships, k-mer hashing can allow for vastly greater compute times and 399 improved detection of both small and structural variation. So far, this method has mostly been applied to 400 population genetics of microbial taxa (Ondov et al. 2016;Kachroo et al. 2019), but has the potential to be 401 useful as more population geneticists adopt whole-genomic methods. In this study, we showed that Mash 402 can accurately reconstruct population relationships from whole-genome data (Figure 3c), though we did 403 not specifically investigate the importance of structural variants. 404 405

Reference-free GWAS 406
A breakthrough study recently showed that k-mer hashes can be used in the place of SNPs to perform 407 reference-free GWA (Voichek & Wiegel 2020). The authors showed that k-mer based GWAS 408 outperformed SNP-based methods, particularly in confidence of detection of causal variants and in the 409 number of phenotypes caused by structural variants (Voichek & Wiegel 2020). Reference-free GWA is 410 not completely new, but early attempts used computationally inefficient methods that are not scalable to 411 large, complex genomes (Lees et al. 2016). As k-mer counts accurately represent complex polyploid genotypes in an efficient way, they pave the way for large-scale GWA analyses in polyploid systems. For 413 example, just as k-mer counts replace genotypes in the linear mixed models commonly used for GWAS, 414 k-mer-derived distance matrices can replace the relatedness matrix that is used to control for population 415 structure in the same analyses. Critically, neither of these methods rely on an assumption of diploidy. 416 417

Genomics of ancient and degraded samples 418
A major finding in this study is that Mash greatly outperforms alignment-based methods when samples 419 have missing data. This leads to similar issues as those posed in polyploid genomes. Since most 420 phylogenetic analyses use Bayesian methods, not distance-based hierarchical clustering, this may not be 421 an issue in that field. Phylogenetic analyses have indeed been found to be largely robust to missing data 422 systems with a high level of missing data, large genomic arrangements, or data-heavy whole-genome 436 comparisons, the simplicity and speed of these methods may offer at minimum a confirmation of 437 alignment-based results. Other methods will need to be tested against Mash to determine which applies for particular population genetic studies. Similarly, combining machine learning with k-mer based 439 analyses may hold great promise for fast and reliable population genetic estimation.