ngsJulia: population genetic analysis of next-generation DNA sequencing data with Julia language

A sound analysis of DNA sequencing data is important to extract meaningful information and infer quantities of interest. Sequencing and mapping errors coupled with low and variable coverage hamper the identification of genotypes and variants and the estimation of population genetic parameters. Methods and implementations to estimate population genetic parameters from sequencing data available nowadays either are suitable for the analysis of genomes from model organisms only, require moderate sequencing coverage, or are not easily adaptable to specific applications. To address these issues, we introduce ngsJulia, a collection of templates and functions in Julia language to process short-read sequencing data for population genetic analysis. We further describe two implementations, ngsPool and ngsPloidy, for the analysis of pooled sequencing data and polyploid genomes, respectively. Through simulations, we illustrate the performance of estimating various population genetic parameters using these implementations, using both established and novel statistical methods. These results inform on optimal experimental design and demonstrate the applicability of methods in ngsJulia to estimate parameters of interest even from low coverage sequencing data. ngsJulia provide users with a flexible and efficient framework for ad hoc analysis of sequencing data.ngsJulia is available from: https://github.com/mfumagalli/ngsJulia


Introduction
Population genetics, i.e., the study of genetic variation within and between groups, plays a central role in evolutionary inferences.The quantification of genetic diversity serves the basis for the inference of neutral 1 and adaptive 2 events that characterised the history of different populations.Additionally, the comparison of allele frequencies between groups (i.e.cases and controls) is an important aspect in biomedical and clinical sciences. 3 the last 20 years, next-generation sequencing (NGS) technologies allowed researchers to generate unprecedented amount of genomic data for a wide range of organisms. 4This revolution transformed population genetics (therefore also labelled as population genomics) to a data-driven discipline.Data produced by short-read sequencing machines (still the most accessible platform worldwide) consist of a collection of relatively short (approx.100 base pairs) fragments of DNA which are then mapped or de novo assembled to form a contiguous sequence. 4At each genomic position, all observed sequenced reads are used to infer the per-sample genotype (an operation called 'genotype calling') and the inter-samples variability, i.e., whether a particular site is polymorphic (an operation called 'single-nucleotide polymorphism (SNP) calling'). 5 this end, several software packages have been implemented to perform genotype and SNP calling from NGS data, the most popular ones being samtools/bcftools, 6 GATK, 7 and freeBayes. 8When, on average, few reads map at each genomic position (a scenario referred to as 'low-coverage' or 'low-depth'), genotypes and SNPs cannot be assigned with confidence due to the high data uncertainty. 9,10Under these circumstances, statistical methods that integrate data uncertainty into genotype likelihoods and propagate it to downstream analyses have been proposed. 5Software packages like ANGSD 11 and ngsTools, 12 among others reviewed by Lou et al., 13 implement a statistical framework to estimate population genetic metrics from low-coverage sequencing data.Similarly, an affordable generation of sequencing data from large sample sizes can be obtained via pooled sequencing experiments, where assignment of individual samples is typically not retained. 14Several new and popular software for the analysis of pooled sequencing have been proposed in recent years. 15,16spite these advances, most of these implementations are either tuned and suitable for model organisms only (e.g., with haploid or diploid genomes) or not easily adaptable to novel applications.Therefore, an accessible computational framework for building and testing ad hoc population genetic analyses from NGS data is in dire need.Among programming languages, Julia 17 has emerged as both a powerful and easy-to-use dynamically typed language that is widely used in many fields of data sciences, including genomics. 18While several Julia packages are currently available for both population genetic and bioinformatic analyses (e.g., BioJulia), to our knowledge, a suitable framework for custom population genetic analysis from NGS data is not available yet.
Here we present ngsJulia, a set of templates and functions in Julia language to process NGS data and create custom analyses in population genetics.To illustrate its applicability, we further introduce two implementations, ngsPool and ngsPloidy, for the analysis of pooled sequencing data and polyploid genomes, respectively.By extensive simulations, we show the performance of several methods implemented in these programs under various experimental conditions.We also introduce novel statistical methods to estimate population genetic parameters from NGS data and demonstrate their applicability and suggest optimal experimental design.We finally discuss further directions and purposes for ngsJulia and bioinformatics for NGS data analysis.

Methods
Implementation ngsJulia was built in Julia language (Julia Programming Language, RRID:SCR_021666) and requires the packages 'GZip', 'Combinatorics', and 'ArgParse'.Auxiliary scripts to process output files were built in R (R Project for Statistical Computing, RRID:SCR_001905) version 3.6.3and require the package 'getopt'.ngsJulia receives REVISED Amendments from Version 2 In the revised version, we provide users with a comprehensive online documentation of APIs and CLIs at https://ngsjulia.readthedocs.io.Additionally, all functions have docstring documentation.The methodology behind some of the implementations has been clarified and references on the applicability of said methods have been added.We also rephrased few sentences for clarity as suggested by the reviewer and fixed few typos.We believe that these changes significantly enhance the user's experience and simplify the interpretation of results obtained using ngsJulia.
gzipped mpileup input files which can be generated using samtools. 19Output files are in text file format and can be easily parsed for producing summary plots and for further analyses.Scripts for some downstream analyses are provided in ngsJulia.
Operation ngsJulia is compatible with all major operating systems and is maintained at https://github.com/mfumagalli/ngsJulia.Documentation and tutorials are available via this GitHub repository and archived at Zenodo 20 at the time of writing.All analyses in the manuscript are performed in Julia language and R. ngsJulia implements functions to read and parse gzipped mpileup files and to output gzipped text files on various calculations (e.g., genotype and allele frequency likelihoods) and estimations (e.g., allele frequencies), as requested by the user.ngsJulia also allows for several data filtering options, including on global and per-sample depth, proportion or count of minor allele, and base quality.Finally, several options for SNP and biallelic and triallelic polymorphisms calling are available.
Nucleotide, genotype and allele frequency likelihoods ngsJulia provides utilities to calculate nucleotide and genotype likelihoods, i.e., the probability of observed sequencing data given a specific nucleotide or genotype, 21 for an arbitrary ploidy level, as in Soraggi et al. 22 We now describe how such quantities are calculated in ngsJulia.
Following the notation in Soraggi et al., 22 for one sample and one site, we let O be the observed NGS data, Y the ploidy, and G the genotype.Therefore, G has values in 0, 1,…, Y f g , i.e., the number of derived (or alternate) alleles.
In the simplest form, genotype likelihoods can be calculated by considering individual base qualities as probabilities of observing an incorrect nucleotide. 21We adopt the calculation of genotype likelihoods for an arbitrary ploidy level P OjG, Y ð Þas proposed in Soraggi et al. 22 From the genotype likelihoods with P OjG,Y ¼ 1 ð Þ , the two most likely alleles are identified by sorting P OjG, Y ¼ 1 ð Þvalues after pooling all sequencing reads together across all samples.This operation will restrict the range of possible genotypes to biallelic variation only.Note that this calculation is still valid for monomorphic sites, although the actual assignment of the minor allele is meaningless.
We now describe how to estimate F a,n , the frequency of allele a∈ A, C,G,T f gat site n.Similarly to Kim et al., 23 the loglikelihood function for F a,n is given by: where c i,n is the i À th read at site n, and C n is the total depth across all samples.Function 1 is maximised to obtain a maximum likelihood estimate (MLE) of the sample allele frequency, b F n , either with a grid-or golden-section-search algorithm.

SNP calling
To perform SNP calling, we implement a likelihood-ratio test (LRT) with one degree of freedom with null hypothesis H 0 : F n , as described by Kim et al. 23 Additionally, we develop a test for a site being biallelic or triallelic.The former can be interpreted as a further evidence of polymorphism, while the latter as a condition not to be met for the site being included in further estimations, as our models assume at most two alleles.The log-likelihood of site n being biallelic is equal to Þ , with i, j, and z being the most, second most, and third most likely allele with Y ¼ 1 (i.e.haploid genotype G), respectively.An LRT with one degree of freedom can be conducted to assess whether Allele frequency, site frequency spectrum and association test from pooled sequencing data Several estimators of population parameters from pooled sequencing data are implemented in ngsPool, a separate program which uses functions in ngsJulia.We now describe the statistical framework for the analysis of pooled sequencing data.
In case of data with unknown sample size, the MLE of the population allele frequency F n is calculated as in Equation 1with F∈ 0, 1 ð Þ.With known sample size, the same equation is used to calculate sample allele frequency likelihoods 24 for instance with f ∈ 0, 1,…, M Â Y f g for M samples of equal ploidy Y. From these likelihoods, we can calculate both the MLE and the expected value with uniform prior probability as estimators of F n .
A simple estimator of the site frequency spectrum (SFS) from pooled sequencing data is obtained by counting point-estimates of F n across all sites.We propose a novel estimator of the SFS implemented in ngsPool.Under the standard coalescent model with infinite sites mutations, we let the probability of derived allele frequency F in a sample of 25 The parameter K determines whether the population is deviating from a model of constant effective population size.For instance, K ¼ 1 is equal to the expected distribution of P F ð Þ under constant population size, while K > 1 models a population shrinking and K < 1 population growth.
We optimise the value of K to minimise the Kullback-Liebrel divergence between the expected distribution of P FjK ð Þ and the observed SFS.The latter can be obtained by either counting b F n across all sites or by integrating over the sample allele frequency probabilities with a uniform prior distribution).A threshold can be set to ignore allele frequencies with low probability to improve computing efficiency and reduce noise.Within this framework, folding spectra can be generated in case of unknown allelic polarisation.
Finally, we introduce a strategy to perform association tests from pooled sequencing data.Similarly to Kim et al., 23 we propose an LRT with one degree of freedom for null hypothesis H 0 : f cases ¼ f controls and alternate hypothesis H 1 : f cases 6 ¼ f controls .The likelihood of each hypothesis is calculated from P O n jF n ¼ f ð Þand, therefore, this strategy avoids the assignment of counts or per-site allele frequencies.A statistically significant LRT with one degree of freedom suggests a difference in allele frequencies between cases and controls, and possible association between the tested phenotype and alleles.

Ploidy levels and test for multiploidy
We now describe the statistical framework implemented in the program ngsPloidy to estimate ploidy levels and test for multiploidy.When multiple samples are available, two scenarios can be envisaged: (i) all samples have the same ploidy, (ii) each sample can have a different ploidy (multiploidy, i.e. as in tumor genomes). 26e log-likelihood function for a vector of ploidy levels f g for M samples and N sites is defined as: F n is the genotype probability given the ploidy and allele frequency at site n.
Once b F n is estimated and the inbreeding is known (or under the assumption of Hardy-Weinberg Equilibrium (HWE)), then genotype probabilities are fully defined.Equation 2 is optimised by maximising the marginal likelihood of each sample separately, assuming independence among samples and sites.

With limited sample size, b
F n is not a good estimator of the population allele frequency and therefore genotype probabilities may not be well defined.In the simplest scenario, genotype probabilities can be set as uniformly distributed, with all genotypes being equally probable.However, the assignment of alleles into ancestral (e.g., wild-type) and derived (e.g., mutant) states is particularly useful to inform on genotype probabilities.Recalling Equation 2, we can substitute Note that E FjK ½ does not depend on the sequencing data for each specific site n.
Note that, in practice, Equation 2 is a composite likelihood function, as samples and sites are not independent observations due to shared population history and linkage disequilibrium, respectively.A solution to circumvent this issue is to perform a bootstrapping procedure, by sampling with replacement segments of the chromosome and estimate ploidy for each bootstrapped chromosome.The distribution of inferred ploidy levels from bootstrapped chromosomes provides a quantitative measurement of confidence in determining the chromosomal ploidy.It is not possible to calculate the likelihood of ploidy equal to one after SNP calling.Nevertheless, the identification of haploid genomes from sequencing data is typically trivial.

Unknown or uncertain ancestral allelic state
So far, we assumed to know which allele can be assigned to an ancestral state, and which one to a derived state.However, in some cases, such assignment is either not possible or associated with a certain level of uncertainty due to, for instance, ancestral polymorphisms or outgroup sequence genome from a closely related species not being available.Under these circumstances, we extend our formulation by adding a parameter underlying the probability that the assigned ancestral state is incorrectly identified.
Let us define R as the ancestral state and a as any possible allele in A,C, G,T f g .In practice, a can take only two possible values as we select only the two most likely alleles.We label this set of the two most common alleles as A and we assume that the true ancestral state is included in such set.The log-likelihood function of ploidy for a single sample m under unknown ancestral state is: Where P R ¼ a ð Þ indicates the probability that allele a is the ancestral state, and it is invariant across sites.If P R ¼ a ð Þ¼0:5, then the equation refers to the scenario of folded allele frequencies, where each allele is equally probable to be the ancestral state.
Finally, note that in a sufficiently large sample size, the major allele is more probable to represent the ancestral state.This probability depends on the shape of the site frequency spectrum, and it is equal to the cumulative distribution of P FjK ð Þ evaluated at F ¼ N=2.We can extend Equation 3 to reflect this parameter uncertainty with a being the major allele in P R ¼ a ð Þ.

Test for multiploidy
We introduce a novel test for multiploidy.If all samples have the same ploidy y∈Y, then Y i ¼ Y j ¼ y is true for all i, j ð Þ∈1,2, …,M.We propose an LRT for multiploidy with null hypothesis H 0 : supY

Data simulation
To benchmark the performance of methods implemented in ngsJulia, we simulated NGS data following a strategy previously proposed by Fumagalli et al. 27 available as a stand-alone R script.Briefly, individual genotypes are drawn according to probabilities depending on input parameters.The number of mapped reads at each position is modelled with a Poisson distribution and sequenced bases are sampled with replacement with a probability given by the quality score.As an illustration, the following code For the analysis of pooled sequencing data, we simulated 100,000 independent sites at sample sizes 20, 50, and 100 from a diploid population with a constant effective population size of 10,000.We imposed the average per-sample sequencing depth to be 0:5, 1, 2, or 5 with an average base quality of 20 in Phred score.To assess the performance of ngsPool, we calculated bias and root mean squared error (RMSE) between the estimated value of the true value, either from the sample or the whole population.While both metrics measure the distance with the true value, the bias retains the direction of the error (i.e.over-or under-estimation).To quantify the accuracy of SNP calling, we calculated F1 scores (the harmonic mean of precision and recall rates).
To simulate data for association test, we assumed an equal number of cases and controls (equal to 150) and 200 SNPs, either causal or non causal.For non causal sites, cases and controls have the same population allele frequency of 0:10.For causal sites, cases and controls have a population allele frequency of 0:09 and 0:04, respectively.These conditions are derived assuming a high risk allele frequency of 0:1, prevalence of 0:2, genotypic relative risk for the heterozygote of 2, genotypic relative risk for the homozygous state of 4. The sample size simulated guarantees at least 80% power with a false positive rate of 0:10.
To illustrate the usage of ngsPloidy, we simulated NGS data of genomes with different ploidy levels (one haploid, eight triploids, one tetraploid) at 1000 sites.NGS data was simulated assuming an average depth of 10 at the haploid level.
Code and simulated data sets analysed are available in ngsJulia GitHub repository.

Empirical data
To showcase the applicability of ngsJulia to real data, we deployed to two distinct empirical data sests.To test ngsPool, we used mapped reads of Arabidopsis lyrata produced in a previous study. 28Specifically, We used data belonging to the population B where 25 individuals of A. lyrata have been sequenced by both pool-sequencing and genotyping-by-sequencing.Since GBS harbours data on a fraction of the genome and Pool-seq is whole genome, We restricted the analysis to the base pairs (bp) sequenced by both techniques, for a total of 15202 bp.The mpileup files have been generated with samtools 6 (v1.14) and SNP calling in GBS data has been performed with Varscan 29 (v2.4.2), resulting in 333 SNPs.
To test ngsPloidy, we reanalysed whole-genome sequencing data of the human fungal pathogen Candida auris from a previous study. 30We processed data for 21 samples, each one consisting of 17 contigs.As the null hypotheses is that these samples are haploid, we tested ploidy levels from 1 to 4. We filtered out sites with a minimum depth lower than 5 and a proportion of minor alleles lower than 0.15.The number of sites analysed per contig ranged from 555 to 754.

Results
ngsJulia implements data structures and functions for an easy calculation of nucleotide and genotype likelihoods (of arbitrary ploidy) which serve the basis of genotype and SNP calling and for the estimation of allele frequencies and other summary statistics.It is particularly suitable for low-coverage sequencing data and for cases when there is high data uncertainty.To demonstrate the use of ngsJulia, we provide two custom applications from its templates and functions.

ngsPool: analysis of pooled sequencing data
We used ngsJulia to implement a separate program, called ngsPool, to perform population genetic analysis from pooled sequencing data.Specifically, ngsPool implements established and novel statistical methods to estimate allele frequencies and site frequency spectra (SFS) and perform association tests from pooled sequencing data.
ngsPool uses functions in ngsJulia to parse mpileup files as input.As an illustration, the following code julia ngsPool.jl--fin test.mpileup.gz--fout test.out.gz--lrtSnp 6.64 will parse test.mpileup.gzfile and write estimates of allele frequencies in test.out.gzfile from unknwon sample size after performing SNP calling with an LRT value of 6.64 (--lrtSnp 6.64, equivalent to a p-value of 0:01).
Depending on the options selected by the user, output files contain, various results are printed on the screen, including • inferred major allele, • inferred minor allele, • LRT statistic for SNP calling, • LRT for bi-and tri-allelic calling, • three estimators of the minor allele frequency at each site.

Estimation of allele frequencies
To illustrate the usage of ngsPool, we estimated allele frequencies based on simulated data at different experimental conditions.We sought to compare the performance among different estimators implemented in the program.If the sample size is not provided, ngsPool provides a simple MLE of the population allele frequency assuming haploidy.Alternatively, if the sample size is provided by the user, ngsPool calculates sample allele frequency likelihoods and returns both the MLE and the expected value of the allele frequency using a uniform prior probability.
Results show that the error of estimating allele frequencies decreases with increasing depth and sample size, as expected (Figure 1).Likewise, the error for estimating the population allele frequency is more pronounced for lower sample sizes.MLE values tend to be less biased than expected values and sample estimates appear to be unbiased even at depth 1 for moderate sample size (Figure 1).
Furthermore, we simulated NGS data at fixed population allele frequency and compared the distribution of true sample allele frequencies and estimated values using a MLE approach from known sample size.Figure 2 shows that most of the deviation from the true population allele frequencies occurs at intermediate frequencies (F equal to 0.5).This effect is more evident for low depth and low sample size (Figure 2).
We then assessed the effect of low-frequency variants on SNP calling.Specifically, we calculated F1 scores for SNP calling when the population allele frequency is 0 (not a SNP) or greater than 0 (a SNP). Figure 3 shows how the prediction performance increases with the population allele frequency (F), sample size, and depth (D).For instance, an F1 score greater than 0.75 is obtained with 20 samples only for F ¼ 0:05 and D > 0:5.On the other hand, the same F1 score is achieved with 50 samples even with F ¼ 0:025 and also at F ¼ 0:02 but only if D > 1.
SNP calling under-performs when there is no variation in the population (F ¼ 0), with sequencing errors and sampling statistical uncertainty generating estimates of F greater than 0.  ngsPool implements several methods to estimate the SFS from sample allele frequencies.As described in the methods, a simpler estimator is based on assigning the most likely sample allele frequency at each site (labelled count).ngsPool implements novel estimators of SFS from pooled sequencing data as described in the method section.A script implements an algorithm to fit the theoretical SFS to the observed SFS.The latter can be calculated either by assigning per-site MLE of allele frequencies (labelled Fit_count) or by integrating the uncertainties across all sample allele frequency likelihoods (labelled Fit_afl).
Figure 4 shows the error in estimating either the population or sample SFS at various settings with different methods.The error decreases with increasing depth and sample size.Estimating SFS by fitting the theoretical SFS without assignment of allele frequencies generally outperforms other tested strategies (Figure 4).
Notably, the novel approach implemented in ngsPool to estimate the parameter K of the SFS distribution (see Methods) allow us to directly quantify the error in inferring demographic events.In fact, all simulations assumed constant population size, equivalent to K ¼ 1. Figure 5 shows the estimated values of K by fitting the SFS either by assigning (counting) allele frequencies (Fit_count) or by using allele frequency likelihoods (Fit_afl).For low-to-moderate depth and sample size, estimates of K tend to suggest population expansion (K < 1), possibly due to an over-estimation of the abundance of low-frequency alleles.However, the error is reduced when integrating the data uncertainty with sample allele frequency likelihoods, as estimates of K values tend to be closer to the true simulated value of 1 (Figure 5).
ngsPool implements a script to perform association tests from pooled sequencing data.Specifically, the script calculates an LRT statistic, with null hypothesis being that allele frequencies of cases and controls (or any two groups) are the same, as used by Kim et al. 23 It uses sample allele frequency likelihoods and, therefore, it maintains data uncertainty and avoids the assignment of counts or per-site allele frequencies.An LRT statistics significantly greater than 0 indicates a difference in allele frequencies between cases and controls.
Figure 6 compares the distribution of LRT statistics between causal and non causal sites at different experimental scenarios.The distribution of LRT at causal SNPs is skewed towards higher values for increasing depth, indicating more support to find phenptype-SNP association.Nevertheless, a clear separation between the distributions of causal and non causal SNPs is observed at low depth (Figure 6).To illustrate the application of ngsPool to real data, we reanalysed genomic data from Arabidopsis lyrata. 28In this study, authors generated data both by genotype-by-sequencing and by pooled-sequencing, with the former providing groundtruth values for genotype and allele calls.Estimates of minor allele frequencies using ngsPool yield a lower RMSE (0.09372) than estimates from Varscan (0.09783) across all SNPs analysed.

ngsPloidy: analysis of sequencing data from polyploid genomes
We further utilised ngsJulia to implement an additional program, ngsPloidy, for the estimation of ploidy from unknown genotypes.The method implemented is similar to the one proposed by Soraggi et al. 22 with some notable differences on the calculation of genotype probabilities (see Methods   Following Equation 2, the genotype probabilities for each tested ploidy are pre-calculated using a script provided in ngsPloidy.This script takes as input the value of parameter K (the shape of the expected SFS), the effective population size, and the probability that the major allele is the ancestral allele.The latter can be either set by the user (e.g., a value of 0.5 would be equivalent to unknown polarisation, as in a folded SFS) or be calculated from the expected population SFS itself.If genotype probabilities are not set, a uniform distribution is assigned.Further options allow for estimation only on called SNPs and/or genotypes.
ngsPloidy uses functions in ngsJulia to parse mpileup files as input.At the end of the computation, various results are printed on the screen, including • the number of analysed sites that passed filtering for each sample, • a matrix of ploidy log-likelihoods for each sample, • the log-likelihood and MLE of the ploidy vector (i.e., the individually estimated ploidy for each sample), • LRT scores for the test of multiploidy against all tested ploidy levels.
Additionally, if requested by the user, ngsPloidy can generate output files with several statistics for each site (e.g., estimate allele frequency), and all per-site genotype likelihoods for each sample and the tested ploidy.
To illustrate the usage of ngsPloidy, we deployed it to simulated data of an multiploid sample consisting of one diploid, eight triploid and one tetraploid genomes.We compared the performance of ploidy and multiploidy inference among different choices of genotype probabilities.The latter were derived either from the expected folded population site frequency, from the estimated ancestral population allele frequency, or from the calculated sample allele frequency.
In all tested cases, we inferred the correct vector of marginal ploidy levels.We therefore assessed the confidence in such inference by calculating the LRT statistics of ploidy and multiploidy inferences.Both were calculated by comparing the most against the second most likely vector of ploidy levels or the most likely vector of equal ploidy levels, respectively.
Results show that using the per-site estimate sampled allele frequency yields higher LRT statistics (and therefore confidence) than using expected population allele frequencies (Table 1).We reiterate that for all three cases we correctly identified the patterns of multiploidy.However, we should caution that with lower sample sizes we do not expect inferences using estimated sample allele frequencies to perform well.
To illustrate the applicability of ngsPloidy to real empirical data, we inferred ploidy levels on 21 isolates from an outbreak of Candida auris. 30In all contigs analysed, haploid was inferred as the most likely ploidy and multiploidy was rejected in all samples.These results are consistent with findings from the original study, although the identification of polyploidy in other isolates may be associated with multidrug-resistant phenotypes. 31

Discussion
Analyses presented here provide further support for the use of genotype and allele frequency likelihoods in the analysis of NGS data. 5Notably, we demonstrated how probabilistic estimates of population genetic parameters can be obtained in case of pooled sequencing data and short-read data from polyploid genomes.Additionally, we motivated the inference of SFS from allele frequency likelihoods as a direct way to infer demography from raw sequencing data.
ngsJulia offers new possibilities of software prototyping for custom analyses of NGS data for population genetic applications.Furthermore, it facilitates experimental design as it provides a platform to benchmark the efficacy of population genetic analysis from competing sequencing experiments.Finally, ngsJulia has accessible documentation and tutorials to inform users on the theory underpinning the implemented methods.We envisage that further improvements in ngsJulia will include the expansion of suitable input formats and data file types, and the compatibility with additional NGS data type, including from long-read sequencing experiments. 32

Conclusions
In this study, we introduce ngsJulia, a series of templates and functions in Julia language to analyse NGS data for population genetic purposes.We present two implementations for the analysis of pooled sequencing data and polyploid genomes, with the inclusion of novel methods.ngsJulia is a suitable framework for prototyping new software and for custom population genetic analyses from NGS data.
frequency with and without this term.Besides, how would someone with a heterogeneous population set this value?
It would be helpful to illustrate the actual distribution of reads across the genome.It is hard to imagine that with a uniform coverage of 0.5, which is not uncommon these days to receive from a WGS run, a good recall rate of the allele frequency can be achieved.It would, therefore, be helpful to understand the assumed input data better and provide the lambda value of the Poisson model for read simulations.
There is no information provided about the recall rate of simulated SNPs.How many SNPs were detected on the different sequencing depth levels, and how many were falsely discovered?This is a major problem in the other named variant callers at very low sequencing depth levels, as discussed in this manuscript, and it would be of high value to compare the accuracy and sensitivity of the new ngsPool approach and the "old" once, like bcftools and GATK.
Real Arabidopsis data comparing ngsPool to VarScan.As previously mentioned, why not compare also to GATK, FreeBayes, and bcftools?VarScan to my knowledge is not one of the top tools in terms of variant detection.
To give this manuscript more relevance and make the decision to use it or not more accessible for the reader, I would encourage the authors to provide the information if this new package performs more precisely than common tools, if not, or under which circumstances it does.If the accuracy is not improved, what are the other benefits?Time-saving, user-friendly output.1: the order in the graphics is not the same as the one in the description.That can be changed easily and helps the reader.

Is the rationale for developing the new software tool clearly explained? Yes
Is the description of the software tool technically sound?Partly

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Partly
Are the conclusions about the tool and its performance adequately supported by the findings presented in the article?Yes I believe, it is potentially interesting to approach some bioinformatics analyses (such as population genomics) with Julia.However, I am still missing a clear motivation for re-writing e.g.genotype likelihoods and SFS computation in Julia.There are many popgen tools that are fast and accurate enough, e.g.ANGSD, so there should be a clear advantage of ngsJulia in some way compared to other tools.If the motivation is the outstanding speed or accuracy of ngsJulia, it should be clearly demonstrated.So my main concern is the lack of comparison of ngsJulia with other popgen tools.
Next, I would not advertise ngsJulia as a "population genetic analysis of next-generation sequencing data", as from the text I got an impression that this was a quite pool-seq-specific tool.If ngsJulia is more general than a pool-seq data tool, it should be demonstrated and compared with other tools.Otherwise, I would suggest to change the title of the manuscript and make it less general and more specific.For example, accounting for polyploidy is a strength of ngsJulia, perhaps this should be emphasized more in the manuscript.
I found it interesting and novel to determine the optimal K in the 1/f^K approximation of SFS profile with Kullback-Leibler optimization.However, as it seems to me from the Figure 5, this implementation in ngsJulia only works for D ~ 5 and large sample sizes.Perhaps other tools will not do a great job either, but as I mentioned above, I am missing the comparison with other tools.
Minor comment: Page 7 I would encourage the authors to make the command line for ngsPool.jlmore self-explanatory.More specifically, the "-lrtSnp 6.64" flag is hard to grasp.Perhaps re-writing it in a way so that it becomes clear that 6.64 corresponds to p-value = 0.01, and replacing the flag "-lrtSnp 6.64" with "--pvalue 0.01" would be possible?

Stéphane De Mita
Institut national de la recherche agronomique (INRAE), Paris, France The article "ngsJulia: population genetic analysis of next-generation DNA sequencing data with Julia language" by Alex Mas-Sandoval and colleagues provides a library for the Julia addressing management of NGS data (provided in the PileUp format), SNP and genotype calling allowing for arbritary and dataset-varying ploidy, and maximum-likelihood parameter estimation.
I've really discovered the Julia language while reviewing this article and I agree that it makes sense to use this language for developing tools for the analysis of NGS data.Providing a toolkit to perfom SNP and genotype calling under a probabilistic models makes much sense and can be useful for the community.The ngsPool application provided can be of widespread use while ngsPloidy can be a very precious tool for specific systems which lack proper software right now.
Here is a list of random typos or detail-level suggestions: "polidy" in the "Revision" box.
○ Methods: you should add Combinatorics as a requirement.
○ I am not sure to understand the header "Operation".

○
One or two sentences, presenting biological systems where the study of multiploidy (or in general MLE of ploidy) is relevant, would help.
○ I am not sure giving an example of a R command line for simulating data is a good choice.For a Julia package I'd like to see some Julia code example, although it is not absolutely necessary.
○ Discussion: rewrite: "is highly applicable in educational contexts" ○ I do not get this sentence at all: "[ngsJulia] allows for efficient testing of experimental designs and, therefore, would be beneficial for initial planning of any sequencing experiments." My big issue with the paper, which is confirmed by a first look at the package, is that the content of the ngsJulia library isn't described, even succinctly.There should be an API and the article should point to it.
In my view there is still need for more documentation.Tutorials are great and (speaking by experience) they are often essential, but at some point you need for some more rigourous reference documentation (API for the library, CLI for the applications).The Application Programmer Interface should be documented so that the reading planning to write a script using the library can know what functions are available, what input they expect, and if possible what errors can be expected.This remark also applies to the Command Line Interface of the two applications ngsPloidy and ngsPool.The help page produces an option summary, which is fine, but the meaning of the options isn't documented.
Below I report the results of my testing of the three components of the software.
*** ngsPloidy *** The tutorial is a good resource to start, but some points need clarification, and a more rigourous reference manual of the command-line interface would help a lot.
I am not sure that it is a good idea to start by defining environment variables, which make the commands much less friendly.At least, do mention that you are defining environment variables so that unaware users can know where to look at if they need more information.Anyway the julia command should be available is it is installed in the environment which should be considered to be the default.By the way, "Julia language" is a poor choice of words.
I am sorry, but I didn't get at all what "genotype probabilities" where and whether and why would I need to to create them.I have no idea what writePars is doing and -p and -s arguments aren't clearly documented.
There is a "variabile" (instead of "variable") in the in simulMpileup.R manual.The simulation script needs its own documentation.
Personally, I would think that either bam files containing aligned reads, or VCF, would be a better choice than PileUp, because users are more likely to already have those.
Suggestions: you could print progress information while running, and displays error messages in stdout instead of stderr (especially useful since you advice users to redirect results which are printed in stdout.
I tested ngsPloidy using data from my own read simulator.I thought it was a good idea to bring data from an independent source, and also I am more comfortable using it.It generated data from diploid and tetraploid organisms, using a coalescent model to draw allelic frequencies.
Then I ran this type of commands: $ julia ../ngsJulia/ngsPloidy/ngsPloidy.jl --fin diplo.plp.gz--nSamples 20 > diplo.outRunning times were (tens of minutes up to an hour for rather small simulated datasets).I suspect that I should have dropped non-varying sites from the PileUp (since I didn't incorporate errors, 90% of my sites are completely fixed).
When I simulated only diploid or tetraploid individuals, they are completely inferred (although in some parameter settings I had tetraploid individuals estimated as a mixture of tetra-and pentaploid individuals).When I merge diploid and tetraploid individuals in a single PileUp, the software tended to estimate a mixture of diploid and triploid individuals.So I suspect that I am doing something wrong.Most likely I miss an important step of the procedure and the documentation might be to blame.
I got an exception when the pileup file didn't contain reference bases.I understand that this should be an error but it could be displayed in a more accessible manner.
Some test results are reported as "Inf" or "-Inf" which is also not necessarily clear.
I was not totally clear when --lrt* arguments should be used.
I got a problem to understand the difference between maf, saf_MLE and saf_E?It is also a problem with the article, where the logic is mentioned by not really explained.
I ran a simulation with a 20,000 bp genome, 20 diploid individuals, 10% of sites with diallelic polymorphism, 10,000 reads of 50 bp (so an average depth of 25X

Lindsay V. Clark
Roy J. Carver Biotechnology Center, University of Illinois at Urbana-Champaign, Urbana, IL, USA The manuscript by Mas-Sandoval et al. describes ngsJulia, a collection of functions and scripts in Julia and R for estimating genotype likelihoods from short read DNA sequencing data, as well as estimating allele frequencies in pooled samples, and estimating ploidy where it is unknown.It seems that the goal is also to enable biologists proficient in Julia to develop their own custom genotype calling and population genetic analysis scripts starting with the provided functions for estimating genotype likelihoods.As it is, the manuscript and software suffer from two major flaws: (1) the software and its outputs are not well documented for users or developers, and (2) testing is only performed on simulated data, not empirical data.
The software documentation is not very approachable as it currently is.End-to-end tutorials with a real dataset (or something closely resembling a real dataset) are needed for all applications.How are the output files formatted and how might the user further process those files to address biological questions?For example, it seems that one could use ngsJulia to generate a matrix of genotype calls across samples and loci, starting from a set of mpileup files, but from the documentation I have no idea how to accomplish this task.
There are eleven different functions for estimating genotype log likelihoods, depending on ploidy, with four different functions for haploids.This strikes me as poor software design.At minimum the eight calcGenoLogLikeX_MajorMinor functions could be collapsed into one, with ploidy as an added argument to the function.Having a single function would make it much easier to update the code and ensure that there aren't any bugs that prevent it from working consistently across ploidy levels.The hard-coding of tuples for every possible genotype, and the repetitive if-else statements, could be replaced with a more programmatic approach.The programmatic approach with ploidy as a function argument would also enable the software to be used for organisms with ploidy greater than octoploid.

Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Partly
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Bioinformatics, polyploidy, population genetics I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 22 Nov 2022

Matteo Fumagalli
The manuscript by Mas-Sandoval et al. describes ngsJulia, a collection of functions and scripts in Julia and R for estimating genotype likelihoods from short read DNA sequencing data, as well as estimating allele frequencies in pooled samples, and estimating ploidy where it is unknown.It seems that the goal is also to enable biologists proficient in Julia to develop their own custom genotype calling and population genetic analysis scripts starting with the provided functions for estimating genotype likelihoods.As it is, the manuscript and software suffer from two major flaws: (1) the software and its outputs are not well documented for users or developers, and (2) testing is only performed on simulated data, not empirical data.
Thank you for reviewing our manuscript.We addressed all concerns by improving the applicability and documentation of our software and by applying it to empirical data.
The software documentation is not very approachable as it currently is.End-to-end tutorials with a real dataset (or something closely resembling a real dataset) are needed for all applications.How are the output files formatted and how might the user further process those files to address biological questions?For example, it seems that one could use ngsJulia to generate a matrix of genotype calls across samples and loci, starting from a set of mpileup files, but from the documentation I have no idea how to accomplish this task.
We improved the documentation as requested in the specific comments (see below).
Examples from simulated mpileup files under various conditions are provided and now clarified.In our response, we clarified that ngsPloidy does not provide a matrix of genotype calls but rather estimates ploidy levels.
There are eleven different functions for estimating genotype log likelihoods, depending on ploidy, with four different functions for haploids.This strikes me as poor software design.At minimum the eight calcGenoLogLikeX_MajorMinor functions could be collapsed into one, with ploidy as an added argument to the function.Having a single function would make it much easier to update the code and ensure that there We replaced all occurrences of "aneuploidy" in the text and documentation with "multiploidy".
How should the user estimate the shape of the site frequency spectrum and the effective population size, which are required inputs for ngsPloidy?
The need to specify the shape of the site frequency spectrum and the effective population size is only required when the sample size is too low to have a meaningful estimate of the population allele frequency.Whilst these parameters are not known a prior for nonmodel species, in practice we found that, as long as values are within a reasonable range and the coverage is sufficient, there is not a significant bias associated with inaccurate settings of these parameters.In fact, they simply determine the prior probability of sampling allele frequencies and, for most species and populations, this distribution tends to have more mass at lower frequency values.We now indicate the contents of all columns of the simulated file.The -callGeno flag does not output called genotypes but rather estimate ploidy by first assigning (calling) genotypes, and thus disabling the integrating over unknown genotypes.This is now clarified in the documentation.

I tried "Case
Similarly to ngsPloidy, ngsPool is only tested on simulated data.
We now also present an application of ngsPool to empirical data.We reanalysed genomic data from Arabidopsis lyrata (Fracassetti et al. 2015).In this study, authors generated data both by genotype-by-sequencing and by pooled-sequencing, with the former providing ground-truth values for genotype and allele calls.We found that estimates of minor allele frequencies using ngsPool yield a lower RMSE (root mean squared error) than estimates from Varscan across all SNPs analysed.We added a paragraph both in the methods and results sections to describe these new analyses.

Minor comments:
There seems to be a missing right bracket on the left side of equation 2.

Fixed.
The Figure 2 caption needs to more explicitly define the abbreviations MLE and SAF.
We added the acronyms SAF and MLE in the caption.
All of the code to reproduce the analysis in the paper is provided, although there is no documentation explaining the code or indicating which scripts go with which figures.
In https://github.com/mfumagalli/ngsJulia/tree/master/paperwe provide all code to replicate the analyses in the paper.As now indicated, ploidy/do.shwould replicate the simulations and estimation of ploidy variation.We now clarify that in the pool and pool/plots folder each folder/script corresponds to the specific figure in the paper.
Competing Interests: No competing interests were disclosed.
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias • You can publish traditional articles, null/negative results, case reports, data notes and more • The peer review process is transparent and collaborative • Your article is indexed in PubMed after passing peer review • Dedicated customer support at every stage • pre-submission enquiries, contact research@f1000.com large value of LRT is suggestive of multiploidy.Statistical significance can be assessed with the LRT and M À 1 ð Þdegrees of freedom.

Figure 1 .
Figure1.Estimation of minor allele frequencies from pooled sequencing data.Bias and root mean square error (RMSE) against the reference true value from either the population or sample are provided for each value of depth D (the average number of sequenced reads per base pair) and sample size (on column panels) tested.Three estimators of allele frequencies are considered: a population maximum likelihood estimate (MLE) from unknown sample size and the expected value and MLE from known sample size.

Figure 2 .
Figure 2. Distribution of estimated minor allele frequencies from pooled sequencing data.True sample allele frequencies (SAF) and maximum likelihood estimates (MLE) are shown at different fixed population allele frequencies F and sample sizes (on rows, 20 and 50) and depths (the average number of sequenced reads per base pair, on columns, 0.5, 1, 2, and 5).

Figure 3 .
Figure3.Performance of single nucleotide polymorphism (SNP) calling from pooled sequencing data.F1 scores (harmonic means of precision and recall rates) for predicting either a SNP (true population allele frequency F greater than 0) or not (F ¼ 0) are reported are various values of F, sample sizes (on columns, 20 and 50) and depths D (the average number of sequenced reads per base pair).

Figure 4 .
Figure 4. Estimation of site frequency spectrum (SFS) from pooled sequencing data.Estimates based on counting allele frequencies (Count) or fitting from counted allele frequencies (Fit Count) or from allele frequency likelihoods (Fit afl) are calculated.Root mean square error (RMSE) values between true (either population or sample, on rows) and estimated SFS are reported at various depths D and sample sizes (on columns).

Figure 5 .
Figure 5. Estimation of parameter K (used to determine if the population is deviating from constant effective population size) of the site frequency spectrum at different depths D and sample sizes (on rows) from pooled sequencing data.Estimates based on fitting from counted allele frequencies (Fit count) and from allele frequency likelihoods (Fit afl) are reported.Note that the true simulated value of K is 1.

Figure 6 .
Figure 6.Performance of association test from pooled sequencing data.The distribution of likelihood-ratio test (LRT) statistics for casual and non causal single-nucleotide polymorphism (SNP) is reported at different depths (the average number of sequenced reads per base pair, on columns).

2 Reviewer
the rationale for developing the new software tool clearly explained?Partly Is the description of the software tool technically sound?Yes Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?Yes Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool?Yes Are the conclusions about the tool and its performance adequately supported by the findings presented in the article?Partly Competing Interests: No competing interests were disclosed.Reviewer Expertise: bioinformatics, genomics, medical genetics, ancient DNA I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.Reviewer Report 14 August 2023 https://doi.org/10.5256/f1000research.150492.r187509© 2023 De Mita S. This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.Stéphane De Mita Institut national de la recherche agronomique (INRAE), Paris, France After considering the updated version and using the new documentation, I managed to run some tests successfully.However, I couldn't understand the return value of calcGenoLike without extrapolating from what it done in the triploid case in the tutorial.Is the rationale for developing the new software tool clearly explained?Partly Is the description of the software tool technically sound?Partly Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?Partly Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool?Partly Are the conclusions about the tool and its performance adequately supported by the findings presented in the article?Partly Competing Interests: No competing interests were disclosed.Reviewer Expertise: Population genomics, phytopathology, programming (Python, C, C++).I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.Version Report 21 February 2023 https://doi.org/10.5256/f1000research.141471.r160519© 2023 De Mita S. This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
A" in the ngsPloidy tutorial.It took me a few minutes to figure out the (undocumented) contents of test.A.txt.I assume that columns 6-25 are the true genotypes, column 5 is the allele frequency, and column 26 is the sample allele frequency.The tutorial otherwise worked.However, I didn't get anything useful when I added the --callGeno flag.
Additionally, ngsPool can output a file with per-site sample allele frequency likelihoods.The following code julia ngsPool.jl--fintest.mpileup.gz--fouttest.out.gz\\--nChroms 20 --fsaf test.saf.gzwillproduce estimates of allele frequencies from the known sample size (specified by --nChroms 20) and allele frequency likelihoods in test.saf.gzfile.These files can then be used to estimate the SFS and perform an association test using two scripts provided in ngsPool.For instance, the code Rscript poolSFS.R test.saf.gz> sfs.txt will estimate the SFS, while the code Rscript poolAssoc.R test.cases.saf.gztest.controls.saf.gz> assoc.txtwill perform an association test assuming two sets of allele frequency likelihood files, one from cases (test.cases. saf.gz) and one from controls (test.controls.saf.gz).

Table 1 .
Confidence values to the correctly assigned ploidy vector and test for multiploidy.Likelihood-ratio test (LRT) statistics using different methods of incorporating allele frequencies are reported.

Is the rationale for developing the new software tool clearly explained? Yes Is the description of the software tool technically sound? Yes Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others? Partly Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Partly Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes
).It took ~ 4:30 minutes and the results are fitting the simulation parameters.All sites called as SNPs are correct and only occasionnaly a few number (~1/1000) of "real" SNPs aren't called, and whose are always singletons.There is also a good agreement between the 'real' and estimated maf.An API reference was badly needed.I tried to start writing a Julia script analysing one of the datasets used for previous testing and didn't get very far.My lack of knowledge of Julia didn't help, for sure.The last paragraph of the readme "In general, ngsJulia provides templates..." should come first.Is it normal that we need to import Combinatorics even if we are not directly using it (i.e.I reckon that ngsJulia should import it itself).Population genomics, phytopathology, programming (Python, C, C++).
The error message is: ERROR: LoadError: BoundsError: attempt to access SubString{String} at index[2]I thought it was because I was trying to pass a site (represented by a Reads object, I assume), containing a single read base.But in the end it is not this (sites with more bases are not accepted either).Obviously I am making a very basic error but the lack of documentation doesn't help me.Competing Interests: No competing interests were disclosed.Reviewer Expertise:

have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above. Version 1
© 2022 Clark L. This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.