J Barros

D Appadoo

D Mcnaughton

E G Robertson

C Medcraft

R Plathe

P Roy

L Manceron
email: laurent.manceron@synchrotron-soleil.fr

The Rotational Spectrum of Propynal in the 250-700 GHz range using Coherent Synchrotron Radiation Fourier Transform Spectrometry

The Coherent Synchrotron radiation obtained at SOLEIL has allowed, for the first time, measurement of broad band absorption, Fourier Transform spectra in the sub-THz range at 0.001 cm -1 resolution, with a signal-to-noise ratio of 100 attained in a few hours. The technique has been applied to the observation of the pure rotational spectrum of propynal (HC≡CCHO) in the 250 -700 GHz range. The new data have allowed the determination of more precise values for rotational constants and for higher order rotational constants. The potential of the technique for rapid recording of rotational spectra in this energy range is outlined.

Introduction.

Propynal, also known as propiolic aldehyde, is a mildly reactive molecule, which has been detected in interstellar molecular clouds, using its rotational spectrum in the 15-40 GHz range [1][2][3]. The first microwave laboratory studies [4,5] allowed a first determination of rotational constants and some centrifugal distortion constants. Later, Winnewisser extended the analysis of the ground state (GS) to the millimeter wave region (95-200 GHz) and obtained a refined set of rotational and centrifugal constants [6]. More recently, McKellar and coworkers have recorded the IR spectrum in the 650 cm -1 region at 0.003 cm -1 resolution, allowing an analysis of the υ 7 , υ 8 and υ 11 bands. In the course of this latter investigation, the authors noted discrepancies between GS levels predicted using the reference 6 constants and combination differences from the observed infrared transitions, especially for higher K values. The data were thus refitted using Watson's a-reduced Hamiltonian to yield a slightly different set of rotational and quartic centrifugal distortion constants and improved values for the sextic constants including an H k term [7]. Even more recently, Jaman and coworkers obtained new millimeter wave data of propynal in a DC discharge containing propargyl alcohol as a precursor and measured new lines in the 36-180 GHz range. In their analysis, they refitted the rotational and centrifugal distortion constants, this time using the s-reduced Hamiltonian instead but not making use of the ground state combination differences from McKellar and coworkers. In view of these discrepancies, it seems desirable to extend the range of observation of the rotational spectrum at higher resolution than the IR work to help determine more precise higher order centrifugal distortion constants. Also, the direct observation of GS rotational lines should yield a higher accuracy than the use of combination difference of infrared lines measured at lower resolution. In view of the importance of the 100-900 GHz spectral range for new ground-based radioastronomy instruments [9], it is of importance to illustrate with a concrete example the capabilities of synchrotron radiation (SR)-based, broad band, interferometric techniques for recording survey absorption spectra in the sub 1 THz range.

The synchrotron provides a source of infrared radiation that combines the advantages of a high brilliance and a wide spectral coverage [10 and reference therein]. To enhance its intensity in the THz range, a number of synchrotron radiation facilities around the world have studied the production of Coherent Synchrotron Radiation (CSR) obtained by reducing the bunch length so that it becomes comparable to the emitted wavelength [11][START_REF] Müller | Proceedings of the 10th European Particle Accelerator Conference[END_REF][START_REF] Martin | [END_REF][14][15][16][17]. Although this operation mode allows a large flux enhancement in the THz range (•10 10 4 ), the interaction of the electrons with their own electromagnetic field causes a microstructure (microbunching) in the electron bunch, resulting in bursts of THz radiation and with accompanying source instabilities [14]. The reported applications of CSR mostly exploit the intensity increase of this source for low resolution studies, where the stability has less impact on the signal-to-noise (S/N) ratio. For high resolution measurements, the use of CSR can be greatly hampered by these micro-bunching instabilities [10]. Recent developments at SOLEIL, however, demonstrated the potential advantages of CSR (brightness, stability, small beam divergence) which improve the S/N allowing the observation of weak bands even at high resolution [16,17]. Even though the ultimate available resolution (20 -30 MHz) using this technique is far from the sub-MHz resolution of the millimeter wave monochromatic scanning set-up, such studies present complementary properties: the whole spectral range is measured simultaneously, which is useful for the study of reactive species, the line intensities are easier to extract and the large number of lines recorded is advantageous in obtaining accurate constants. Here, such a study is presented for which the CSR emission is ideally suited: the rotational structure of propynal in the sub-THz range.

Experimental.

The absorption measurements of the propynal molecule were performed on the AILES beamline with the CSR source. The flux of this emission between 250 and 750 GHz (about 0.4 to 1.2 mm wavelength) is advantageous for performing broad band interferometric spectroscopy measurements (here about nine hours). For high resolution measurements, (here •30 MHz resolution) the signal-to-noise has been optimized through minimizing the source instabilities [14,16,17]. These optimal conditions allowing for an intense emission over a wide frequency range are found for a momentum compaction factor: α = α 0 /25 (with α 0 the nominal compaction factor) , corresponding to a theoretical bunch length at zero current of around 3 ps (1 mm), and a current per bunch of 70µA. The Synchrotron was operating in the so-called hybrid mode with an isolated electron bunch (in one quarter of the orbit) and 312 equally small bunches. With these conditions, the CSR appears quasi-stable in the 7 to 21 cm -1 range, allowing for high resolution measurements.

The AILES beamline of SOLEIL has been designed to minimize known sources of instabilities, specifically in the far-infrared and THz ranges [18,19], these characteristics are briefly summarized here. Near the bending magnet output, the edge and constant field radiations of the bending magnet are extracted by a slotted mirror. The collected solid angle is 78 mrad horizontally and 20 mrad vertically, allowing 10 16 to 10 17 photons to be extracted per 0.1 % band pass of coherent radiation. The transfer optics placed in the ring vacuum chamber are fixed to the highly stable ground of the storage ring. The experimental beamline is under vacuum (pressure less than10 -6 mbar) and equipped with a Bruker IFS 125HR high-resolution interferometer (882 cm maximum optical path difference) coupled to two 4.2 K Infrared Laboratories composite bolometers with a specified theoretical noise equivalent power of about 10 -13 W/√Hz . The use of separate bolometers for sample and reference allows the non random noise generated in the CSR to be removed from the final spectrum. In the following, the synthesis of the propynal sample is summarized, the experimental measurement conditions are given and the results of this study presented.

To obtain propynal with high purity, we used a synthesis technique based on the retro-ene reaction of propargyl ether, described by McNab et al. [20] followed by a purification procedure. For this purpose 7 g of propargyl ether (98% from Sigma Aldrich) were placed in a flash pyrolysis chamber under vacuum allowing for thermal decomposition at 750°C into two products, the propynal, and allene (H 2 C=C=CH 2 ). After the pyrolysis, the reaction products were trapped at near liquid nitrogen temperature and gradually warmed under static vacuum. At around -34 ° C as the allene evaporates, the gas is pumped away and liquid propynal remains. This purification procedure was performed twice and the purity was checked using both mass spectrometry and infrared absorption, until only propynal was detected.

For the present measurement, 0.225 mbar of propynal was introduced into a White cell aligned to an 8.85 m optical path. The interferometer was equipped with a Si/Mylar composite beamsplitter and the moving mirror scanned at the highest possible velocity of 5.06 cm/s (80 kHz He-Ne laser fringe modulation), as the highest Fourier frequency in the wavenumber range of interest (about 125 Hz) falls well below the detector roll-off (full response here up to 200Hz). This enables averaging over a larger number of interferograms. The detector was fitted with a 200 cm -1 cold low pass filter and the spectrometer with a 12.5 mm entrance aperture, but the SR effective source diameter was measured at the interferometer entrance focal plane to be about 4.5mm FWHM with 100% of the incoming energy at 8mm diameter. This results in a source diameter fulfilling the resolution criterion. For the propynal spectrum and reference spectrum (empty cell), 186 scans at maximum resolution were recorded for about 9 hours. The calculation of the transmittance may ideally be obtained when measuring both empty and gas filled cell in strictly the same optical conditions, aiming at suppressing the fine structures discussed previously [14]. The S/N ratio at the peak of the bandpass (around 10 cm -1 ) can be estimated to about 57. It should be noted, however that, during the course of preliminary measurements, an increase in signal was observed, causing one detector to saturate. This gradual increase in coherent emission was obtained as the electron bunch length distribution narrowed as a result of the topping-up procedure and gave rise to a 20% general increase of the signal, causing the detectors to reach saturation point. In order to operate the detector in non-saturated conditions, we chose to slightly misalign the cell to record the gas filled cell spectra. The signal to noise ratio of the transmission, based on this ratio, was then about 45 at best. The spectrum calibration was performed using CH 3 Cl lines measured just prior to this study. Water lines were not exploited for this purpose as this molecule is present in both the spectrometer and the cell, to be avoided for a precise line determination [START_REF] Davies | Fourier Transform Spectrometry[END_REF]. Using CH 3 Cl data in the HITRAN [START_REF] Rothman | The HITRAN 2004 molecular spectroscopic database[END_REF] database, from the microwave measurements of Wlodarczak et al. [23], a recalibration shift of 3.5×10 -4 cm -1 (10.5 MHz) was determined. This adjusts for a small systematic error in the frequency axis, with random errors on line positions of the calibration gas of the order of 5×10 -5 cm -1 (1.5 MHz). This small systematic error is commonly due to slight non parallelism of the alignment laser and IR moving mirror axis and for a small part (estimated to about 2.2×10 -4 cm -1 here) to slight distortion of the line shape due to off-focus and off-axis radiation [START_REF] Davies | Fourier Transform Spectrometry[END_REF]24].

Results and Analysis.

The geometry of propynal makes it a near prolate asymmetric top with the A rotational constant much larger than the B and C constants, which are of the same order of magnitude (Figure 1) and a dipole moment with µ a and µ b components. As a result, the rotational spectrum recorded here in the 8 to 19 cm -1 range displays quite regularly spaced q R Ka stacks arising from the µ a selection rule, and also a few r Q Ka sub branches from the µ b selection rule, where the overlapping lines lead to strong absorption corresponding to the high intensity features, with the r Q 3 sub branch saturating in the spectrum. Propynal has several low lying fundamental vibrations (in particular ν 12 and ν 9 below 300 cm -1 ) with a substantial population at room temperature, responsible for the intricate pattern of weaker lines between the well-defined K stacks.

The analysis was carried out as follows: since very good starting parameters were available in the literature, we first carried out a simulation in the spectral region of interest, using the parameters from reference [7] and the same A-reduced asymmetric rotor Hamiltonian [27] and I r representation as in this earlier work. Simulations were performed using C.M. Western's PGopher program [START_REF] Pgopher | a Program for Simulating Rotational Structure[END_REF] and a first set of about 600 transitions could be assigned in a straightforward manner.

The RMS error produced when using the ground state parameters from literature to calculate assignable transitions in the 250-700 GHz region were relatively large (about 0.002 cm -1 for constants from ref [6], about 0.01 cm -1 for constants from ref [7] and about 0.0025 cm -1 for constants from ref [8]). As expected, examination of the residuals showed that although the results for the already measured microwave lines were, as expected, excellent, systematic deviations appeared in the higher quantum number regions. Thus, using our experimental data, a new fit of the molecular parameters was carried out using H.B. Pickett SPFIT and SPCAT program suites [START_REF] Pickett | [END_REF]. This produced a new set of parameters which allowed a better match of predicted line positions with higher quantum number lines. New sets of lines were progressively introduced and this procedure repeated until about 1200 transitions (941 a-type and 257 b-type lines) could be reproduced with an average RMS deviation of about 2×10 -4 cm -1 . A total of 140 pure rotational lines from refs [4][5][6]8] were also included in the combined fit. An example comparison of simulated and observed spectrum is given in Figure 3. The results are also compared to literature values using the residuals plots in figure 4. The new set of effective constants is presented in table 1 and compared with literature values. 3. In red, sub-200 GHz lines from references [6] and [8], in black, 250-700 GHz lines from this work.

On the whole, the precision of the rotational and centrifugal distortion constants is improved. Some quartic and all sextic constants now have improved values, well outside the defined error bars from earlier work and a H J sextic constant is added. Addition of other, higher order, constants did not significantly improve the fit and were not determined with enough certainty anyway, and these have not been retained. It should also be noted that the RMS error using these parameters is close to the precision of our frequency determination (about 10 -4 cm -1 ), and that the residuals are now evenly distributed around zero without systematic deviation.

In the simulation process, the relative intensities of the a-type and b-type transitions were manually adjusted to yield the best visual agreement. Although this procedure does not take into account the instrument function rigorously, it gives however, an acceptable approximation of the relative intensities for weak (unsaturated) lines. From the relative intensities, we could thus derive a µ a /µ b = 1.57± 0.1 dipole moment component ratio. This determination compares favourably with the µ a /µ b = 2.36 ±0.02/1.47±0.02 = 1.60 ± 0.04 ratio, determined previously with great precision from Stark effect measurements by Brown and Godfrey [25]. This also implies that, using with standard procedures with precise optical path, pressure and temperature measurements [30], line strength measurements could be carried out and from this, estimation of total dipole moments. Here, as an illustration, an approximate measurement of intensity over four isolated lines with moderate intensity (about 1.9×10 -22 cm 2 .molec. -1 cm -1 ) around 10 cm -1 leads to a dipole moment estimate of 2.7± 0.3 D (compared to 2.78D experimentally [25]). Clearly, measurements over several pressures are needed for an adequate measurement of the dipole moment or, in other cases, the vibrational partition function, but this illustrates this possibility.

As explained earlier, many lines pertaining to vibrationally excited states are observed in the spectrum, these could be extracted by subtracting out the ground state predicted lines [31], but in this paper no attempt was made to analyse these data.

Conclusions

This study shows that Coherent Synchrotron Radiation can be used to obtain, in a few hours, rotational spectra of a light, reactive molecule such as propynal in the sub-THz range using FT-based instrumentation. Although this technique has limited resolution capabilities compared to purely electronic or opto-electronic scanning devices, it is shown here that the large number of observable lines in this spectral range enables molecular constants to be retrieved with good precision, as a new, improved set of effective molecular constants could be derived. In a similar manner, improvement in the rotational constants was achieved using FIR data for ethylene oxide, a lighter molecule where the higher frequency lines are available using normal synchrotron radiation [32]. Also, the multiplex advantage of the technique may be of use when studying reactive or elusive molecules with lifetime of the order of an hour. Intensity measurements are also possible using the same methods as in the IR region, which can be an advantage for sub-mm wave/ IR cross calibration studies, as switching between these two spectral ranges can be easily done, using the very same instrumentation. a. Converted from the older set of centrifugal distortion constants. Here k= -2R5 see ref [33]. b. Not specified. From MW lines of ref [6] and combination differences in IR bands.

Highlights

-Improved molecular constants for Propynal -Use of Coherent Synchrotron radiation for broad band rotational spectroscopy -Fourier Transform spectra in the 250 -700 GHz range at 0.001 cm -1 resolution

Figure 1 .

 1 Figure 1. Molecular structure of propynal with the main axis and the dipole moment vector. The interatomic distances and angles are given in Table 1 (from references 5, 25 and 26).

Figure 2 .

 2 Figure 2. Survey of the rotational spectrum of propynal in the 8-19 cm -1 region, recorded at room temperature with 0.001 cm -1 resolution, using CSR at Beamline AILES (SOLEIL), 0.225 mbar sample pressure and 8.85 m path length.

Figure 3 .

 3 Figure 3. Detail of the rotational spectrum of propynal for the J = 40 ← 39 region. Only lines marked with an asterisk have been retained in the fit, most others being excited vibrational state transitions. Experimental conditions are the same as in Fig. 2. Gaussian line shapes with 4.5×10-4 cm -1 HWHM were assumed for the ground state line simulations.

Figure 4 .

 4 Figure 4. Plot of the fit residuals (observed minus calculated line positions) using (a) constants, reduction and representation as in reference[8], (b) constants from reference[7], (c) this work, constants from table3. In red, sub-200 GHz lines from references[6] and[8], in black, 250-700 GHz lines from this work.

Table 1 .

 1 Comparison of ground state parameters for Propynal using Watson's a-reduced Hamiltonian in the I r representation (1σ uncertainty in parenthesis).

	Constant	This work	Ref [6] a	Ref [7]
	(MHz)
	A	68035.308 (25)	68035.299 (43)	68035.263 (39)
	B	4826.28805 (97)	4826.3014 (73)	4826.293 (6)
	C	4499.52695 (95)	4499.5107 (69)	4499.51 (6)
	ΔJ	1.91792 (111) • 10 -3	1.917 (012) • 10 -3	1.913 (4) • 10 -3
	ΔJK	-0.1481185 (186)	-0.148102 (47)	-0.148044 (39)
	ΔK	9.0033 (57)	8.99124 (97)	8.9926 (87)
	δJ	3.3742 (15) • 10 -4	3.458 (14) • 10 -4	3.4564 (12) • 10 -4
	δK	0.03396 (26)	0.0390 (32)	0.0347 (28)
	HJ	7.549 (214) • 10 -9	--	--
	HJK	-0.8348 (082) • 10 -6	-7.41 (59) • 10 -7	-7.11 (48) • 10 -7
	HKJ	-8.245 (107) • 10 -6	-8.310 (130) • 10 -6	-8.264 (12) • 10 -6
	HK	3.190 (147) • 10 -3	--	2.36 (6) • 10 -3
	J max	77	23	<50-b
	Ka max	14	14	-b
	Kc max	75	21	-b
	Number of a-type Lines	941	91	-b
	Number of b-type Lines	257	44	-b

Acknowledgement.

The authors wish to thank the Australian synchrotron and the Soleil synchrotron for the collaboration and the Australian government for funding through the International Science Linkages (ISL), French-Australian science and technology (FAST) program.


Luca Ferretti
email: luca.ferretti@gmail.com

Sebástian E Ramos-Onsins

A generalized Watterson estimator for next-generation sequencing: from trios to autopolyploids

Keywords: Site frequency spectrum, Population genetics, Summary statistics, Maximum likelihood, Composite likelihood

Several variations of the Watterson estimator of variability for Next Generation Sequencing (NGS) data have been proposed in the literature. We present a unified framework for generalized Watterson estimators based on Maximum Composite Likelihood, which encompasses most of the existing estimators. We propose this class of unbiased estimators as generalized Watterson estimators for a large class of NGS data, including pools and trios. We also discuss the relation with the estimators that have been proposed in the literature and show that they admit two equivalent but seemingly different forms, deriving a set of combinatorial identities as a byproduct. Finally, we give a detailed treatment of Watterson estimators for single or multiple autopolyploid individuals.

Introduction

The rescaled mutation rate per base θ plays an important role in population genetics models. Its definition is θ = 2pN e µ where µ denotes the mutation rate per base, N e the effective population size and p the ploidy of the species. In the context of the Standard Neutral Model (SNM) at low mutation rate, commonly used estimator include the Watterson estimator [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], based on the number of segregating sites S in a sample from the population, and the pairwise nucleotide diversity Π [START_REF] Tajima | Evolutionary relationship of dna sequences in finite populations[END_REF], defined as the average number of differences per base between two individual sequences. Despite its simple interpretation and its popularity, Π is an inconsistent estimator of θ, while the Watterson estimator is a good estimator since it is unbiased and it corresponds to the Maximum Composite Likelihood estimator for θ. Furthermore, as shown in [START_REF] Roychoudhury | Sufficiency of the number of segregating sites in the limit under finite-sites mutation[END_REF], it is a sufficient statistics for θ for large sequences. For unequal mutation rates between different alleles, the Watterson estimator actually measures the net rate of mutation towards different alleles, rescaled by population size.

Today, variability analyses on a genome-wide scale can be easily done by Next Generation Sequencing (NGS) technologies. NGS technologies can sequence a single complete genome at relatively high redundancy, but the sequencing of many samples increases substantially the cost of the experiment. Several strategies have been used to obtain sequence data from many individuals, from restriction reduced libraries to pooled samples. These strategies may reduce substantially the coverage across the genome (i.e. how many regions are sequenced) and the sequence redundancy per nucleotide base (i.e. how many sequences cover each position). From this point of view, the study of variability on NGS data can be seen as a missing data problem, where information about several samples is missing at covered positions. The Watterson θ estimator has been generalized to missing data by Ferretti et al. [START_REF] Ferretti | Neutrality tests for sequences with missing data[END_REF], but it covers just a limited number of cases.

Estimators of variability based on NGS data should take into account the relatively high probability of sequencing errors in order not to overestimate the number of variants. Sequence redundancy at a given position is fundamental for detecting and removing errors but also for detecting correctly homozygote or heterozygote positions in diploid individuals. NGS estimators of variability were first proposed by Lynch [START_REF] Lynch | Estimation of nucleotide diversity, disequilibrium coefficients, and mutation rates from high-coverage genome-sequencing projects[END_REF] for a single diploid individual using the sequence error rate and the number of reads observed at each variant. Hellmann et al. [START_REF] Hellmann | Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals[END_REF] and Jiang et al. [START_REF] Jiang | Population genetic inference from resequencing data[END_REF] proposed Watterson estimators for multiple individuals taking into account the combinatorics of reads from different individuals and homologous chromosomes. Later, Futschik and Schlötterer [START_REF] Futschik | The next generation of molecular markers from massively parallel sequencing of pooled dna samples[END_REF] and Ferretti et al. [START_REF] Ferretti | Population genomics from pool sequencing[END_REF] developed variability estimators for pooled sample data using Method of Moments (MM) and Maximum Composite Likelihood (MCL) methods, respectively. We observe that for the analysis of a single diploid individual, they all reduce to the same estima-tor except for [START_REF] Futschik | The next generation of molecular markers from massively parallel sequencing of pooled dna samples[END_REF] (the differences between this estimator and the MCL have been detailed in [START_REF] Ferretti | Population genomics from pool sequencing[END_REF]). In this work, we study the relations between these estimators and present a common framework for estimators of variability. We develop new Watterson estimators for other kinds of NGS data, like trios or autopolyploids.

In our framework, we deal with all these cases as special instances of a general, unified approach. Data are represented by NGS or Sanger sequences coming from several units. Units can be haploid, diploid or polyploid individuals or pools, each one containing a different number of lineages. In section 2, we derive Maximum Composite Likelihood estimators for a generic class of data. These MCL estimators are not unbiased, but we show that they can be well approximated by unbiased estimators that share the same functional form as the Watterson estimator. In section 3, we propose these unbiased estimators as generalized Watterson estimators for a large class of NGS data, such as multiple haploid/diploid/polyploid individuals, pools, trios, inbred lines, and combinations of them. We present explicit formulae for these estimators in section 4. We discuss their relation with other estimators that have been proposed in the literature, showing in section 5 that many of our estimators admit two equivalent but seemingly different forms. As a byproduct, this equivalence implies a set of combinatorial identities. Finally, in section 6 we treat in details the case of single and multiple autopolyploid individuals and we provide the Watterson estimator for autopolyploids.

General Watterson estimators

Maximum Composite Likelihood estimators

Composite Likelihood Estimation of parameters has been extensively used in population and quantitative genetics for estimating the linkage disequilibrium among positions [START_REF] Devlin | Disequilibrium mapping: composite likelihood for pairwise disequilibrium[END_REF][START_REF] Weir | Inferences about linkage disequilibrium[END_REF] and for estimating evolutionary parameters as the level of variation [START_REF] Hellmann | Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals[END_REF], the recombination rate [START_REF] Hudson | Two-locus sampling distributions and their application[END_REF], [START_REF] Mcvean | A coalescent-based method for detecting and estimating recombination from gene sequences[END_REF], the strength of positive selection [START_REF] Kim | Detecting a local signature of genetic hitchhiking along a recombining chromosome[END_REF][START_REF] Zhu | A composite-likelihood approach for detecting directional selection from dna sequence data[END_REF] or demographic parameters [START_REF] Garrigan | Composite likelihood estimation of demographic parameters[END_REF]. Maximum Composite Likelihood is an appropriate method to estimate the nucleotide variability across large regions because it has minimum mean squared error for large recombining regions, since in this case the Composite Likelihood is a good approximation for the exact likelihood and therefore the estimator is approximately asymptotically efficient. Furthermore, it turns out to be based on the same statistics as the Watterson estimator (the total number of segregating sites S) and it actually reduces to the Watterson estimator if the data are represented by complete sequences.

In this paper, we consider allelic variants represented by segregating sites, or Single Nucleotide Polymorphisms (SNPs), but the methods can be applied to generic variants with low mutation rate.

For each site, there are features that depend on an eventual SNP (for example, allele frequencies) and features that do not depend on the allelic content (for example the read depth, i.e. the number of sequences that contain data for a given site). We summarize the SNP features with the index ξ ∈ Ξ and the features of each site that do not depend on the allelic content with the index ϕ ∈ Φ. Both indices indicate mutually exclusive, collectively exhaustive features.

We denote by p ϕ,ξ (θ) the probability that a site with features ϕ contains an observed SNP with features ξ for the sample studied. For small values of θ (that is, in the infinite site model), we can expand it in Taylor series and using the fact that there are no SNPs without mutations, i.e. p ϕ,ξ (0) = 0, we find that these probabilities are proportional to θ multiplied by a quantity that depends on the population model and the sequencing setup:

p ϕ,ξ (θ) θZ ϕ,ξ . (1)
We denote by S ϕ,ξ the number of segregating sites with features ξ in positions with features ϕ, and by L ϕ the number of sites with features ϕ. We also define the quantities S ϕ = ξ∈Ξ S ϕ,ξ and Z ϕ = ξ∈Ξ Z ϕ,ξ , the total number of segregating sites S = ϕ∈Φ S ϕ and total length L = ϕ∈Φ L ϕ . Under the composite approximation, all sites are independent. The Composite Likelihood is therefore the simple product of probabilities:

CL(θ) = ϕ∈Φ ξ∈Ξ p ϕ,ξ (θ) S ϕ,ξ •   ϕ∈Φ 1 - ξ∈Ξ p ϕ,ξ (θ) Lϕ-Sϕ   . (2)
Substituting eq. ( 1) for p ϕ,ξ (θ) and taking the log, we obtain

log(CL(θ)) = S log(θ) + ϕ∈Φ ξ∈Ξ S ϕ,ξ log(Z ϕ,ξ ) + ϕ∈Φ (L ϕ -S ϕ ) log(1 -θZ ϕ ) .
(3) For S > 0, the loglikelihood is always negative and tends to -∞ both for θ → 0 and θ → 1/ max ϕ∈Φ (Z ϕ ), so the maximum can be obtained from the zeros of the first derivative of the log(CL) in equation [START_REF] Roychoudhury | Sufficiency of the number of segregating sites in the limit under finite-sites mutation[END_REF]. After rearrangements, we obtain the equation that defines the Maximum Composite Likelihood Estimator (MCLE):

L = ϕ∈Φ L ϕ -S ϕ 1 -θMCLE Z ϕ , (4)
which is valid for all values of θ < 1/ max ϕ∈Φ (Z ϕ ) since the second derivative in

θ = θMCLE is negative if ϕ∈Φ (L ϕ -S ϕ )Z ϕ /(1 -θMCLE Z ϕ ) 2 > 0, which is always verified.
For the simplest case of the original Watterson estimator [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], all sites are equivalent and there is no site feature ϕ, so the above MCLE equation ( 4) can be easily rewritten as L = (L -S)/(1 -θMCLE Z), i.e. θMCLE = θW = S/(LZ). The SNP features Ξ correspond simply to the derived allele counts i = 1 . . . n -1 and the probability of a SNP of frequency i/n is related to the expected frequency spectrum ξ i -defined as the count of SNPs of frequency i/n in the sample -by p i (θ) = E(ξ i )/L. For the standard neutral model,

p i (θ) = θ/i, therefore Z = n-1 i=1 Z i = n-1 i=1 p i (θ)
/θ is given by the harmonic number a n = n-1 i=1 1/i. Then the MCLE in this case corresponds precisely to the unbiased Watterson estimator θW = S/(L n-1 i=1 1/i) = S/(La n ). The estimator (4) is defined implicitly, so it is not easy to use. An explicit, approximate MCL estimator can be derived in two equivalent ways: either (i) by expanding equation (4) at first order in the small parameters θ and S ϕ /L ϕ with constant ratio S ϕ /(θL ϕ ), or (ii) by taking the small θ, large L limit of the likelihood (2) with θL and L ϕ /L constant; in this limit, the S ϕ,ξ are Poisson distributed random variables with mean L ϕ θZ ϕ,ξ

CL(θ) ϕ∈Φ ξ∈Ξ (L ϕ θZ ϕ,ξ ) S ϕ,ξ S ϕ,ξ ! e -LϕθZ ϕ,ξ = θ S e -θ ϕ∈Φ LϕZϕ ϕ∈Φ ξ∈Ξ (L ϕ Z ϕ,ξ ) S ϕ,ξ S ϕ,ξ !
(5) and since the dependence on θ lies in the first term which is a function of the statistics S only, S is a sufficient statistics for θ by the Fisher-Neyman factorization theorem, as already observed in [START_REF] Roychoudhury | Sufficiency of the number of segregating sites in the limit under finite-sites mutation[END_REF].

Both ways lead to the same estimator. We define the resulting approximate MCL estimator as the generalized Watterson estimator:

θW = S ϕ∈Φ L ϕ Z ϕ . (6)
This estimator depends only on the total number of segregating sites, like the original Watterson estimator, since S is a sufficient statistics for small θ. Furthermore, it is unbiased. In fact, E(S) = ϕ∈Φ ξ∈Ξ E(S ϕ,ξ ) = ϕ∈Φ ξ∈Ξ L ϕ p ϕ,ξ (θ) = ϕ∈Φ θL ϕ Z ϕ . Both the implicit estimator in equation ( 4) and the formula [START_REF] Hellmann | Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals[END_REF], which is unbiased, can be used. The relative error of the Watterson estimator ( 6) with respect to the MCLE ( 4) is very small. By expanding eq. ( 4) at first order in θMCLE -θW and at first nonzero order in θW and S φ /(L φ Z φ ), we obtain an error estimate of order θ multiplied by a weighted covariance between Z and the relative fluctuations of θW for different ϕs: 2 . This error is usually negligible, since the r.h.s. is suppressed by a factor of θ 1; furthermore it has mean 0, since it correlates the fluctuations of the mutation process with the fluctuations of the sequencing process, but the two processes are independent.

θMCLE -θW θW -θW ϕ∈Φ L ϕ Z ϕ ϕ ∈Φ L ϕ Z ϕ Z ϕ - ϕ ∈Φ L ϕ Z ϕ L S ϕ /(L ϕ Z ϕ ) -θW θW (7) up to terms of order ( θW , S φ /(L φ Z φ ))
The only information needed to compute (4) or ( 6) are the factors Z ϕ = ξ∈Ξ p ϕ,ξ /θ, which depends both on the model and on the sequencing setup. In the rest of the paper, we will specialize these factors for several combinations of NGS data.

Application to NGS and sequence data

The data: sequences aligned to a reference genome

In this section we deal with combinations of complete sequences, genotypes and NGS data from different sources in a unified way. Our data are represented by reads, sequences or genotypes1 , aligned to a reference genome. Each read/sequence/genotype originates from a single unit: units can be individuals of different ploidy, or pools of individuals. Complete sequences are considered as sequences coming from an haploid unit, so the two complete sequences of the two homologous chromosomes from a diploid organism are equivalent to two different haploid units.

We denote by U the number of units. The features associated to the units are (i) the number of copies of homologous chromosomes present in each unit and (ii) the evolutionary relationships between the units. We denote the number of copies of homologous chromosomes in the ith unit by c i , i = 1 . . . U . Diploid individuals will have c i = 2, polyploid individuals will have c i equal to their ploidy and pools will have c i equal to the number of individuals in the pool multiplied by their ploidy. We denote the set of numbers {c i } i=1...U by {c}.

The evolutionary relationships (denoted here by the generic symbol χ) include all the available information relevant for the probabilities that some of the sequenced chromosome derived by the same lineage, either because they are actually from the same individual or because are identical by descent. For example, two different pools could contain two genetically identical individuals, or two individuals could be parent and offspring, or a single diploid individual could originate from a single inbred line with given inbreeding coefficient, and so on.

We assume that the number of NGS reads covering each position depends only on the sequencing process and not on the allelic composition of the sequence. In this case, we associate to each position x the read depth of the ith unit r i (x), i = 1 . . . U , i.e. the number of reads or sequences from the ith unit that cover position x. For Sanger sequences and genotyping, we define an "effective read depth' r i (x) = 0 for positions with missing data and r i (x) = +∞ otherwise. We denote the set of read depths {r i } i=1...U by {r}.

An example of the data is given in Table 1.

Estimators for Next Generation Sequencing

We consider NGS data like the ones described above. We can derive the general Watterson estimator for this case by using the definition (6) with ϕ = {r} = {r i } i=1...U , the set of read depths of the different units at a given site. We use the short form {r} = {r i } i=1...U and {c} = {c i } i=1...U for the information about the site features. Z can be computed by conditioning on the number of unrelated homologous chromosomes or lineages which actually contribute to the data, denoted by j, and then averaging over j:

ξ∈Ξ p {r},ξ = θZ {r} = ∞ j=2 P (SNP|χ, {c}, {r}, j) • P c (j|χ, {c}, {r}) (8)
where P (SNP| . . .) is the probability of observing a SNP among j independent lineages and P c (j| . . .) is the distribution of j, which depends also on {c}

Reference A C A C G T A A T C G C Sanger: (unit 1) A C A C G T T A A C G C (unit 2) A C A C C T -A T C G C Genotyping: (unit 3) A A C C - C C C G T T T T A A T T C C G G C G NGS reads: (unit 4) A C G T T A (unit 4) A C A C C (unit 4) C C T T A A (unit 5) T A A C G C Read depths: unit 1: c 1 = 1 r 1 = ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ unit 2: c 2 = 1 r 2 = ∞ ∞ ∞ ∞ ∞ ∞ 0 ∞ ∞ ∞ ∞ ∞ unit 3: c 3 = 2 r 3 = ∞ ∞ 0 ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ ∞ unit 4: c 4 = 2 r 4 = 1 1 2 3 3 2 2 2 1 0 0 0 unit 5: c 5 = 2 r 5 = 0 0 0 0 0 0 1 1 1 1 1 1
Table 1: Example (not real) of data from five sequencing units from different sources and technologies. The data come from four diploid individuals, of which two are parent and child. The data are aligned to the reference genome (the sequence at the top). There are two complete Sanger sequences (units 1 and 2) with a missing base in the 7th position, one sequence of genotypes (unit 3) with a missing genotype at the 3rd position, and four NGS read (three coming from unit 4 and one from unit 5). The allelic content of the positions with SNPs is shown in bold red. At the bottom, we report the "effective' read depths for all units. All units are diploid individuals (therefore c = 2). Units 1 and 2 are the two homologous sequences of the parent (hence c 1 = c 2 = 1) and unit 5 comes from the child (c 5 = 2), while units 3 and 4 come from unrelated individuals (c

3 = c 4 = 2).
and on their relationships χ. The probability of observing a SNP among j independent lineages depends just on j, not on the redundancy of each lineage in the sequences, and on the expected site frequency spectrum E(ξ k |j) of the model, being equal to P (SNP|j) = j-1 k=1 E(ξ k |j)/L. In the case of the standard neutral Wright-Fisher model we have E(ξ k |j) = θL/k and therefore

P (SNP|j) = θ j-1 k=1 1/k = θa j . The general Watterson estimator for the Standard Neutral Model (SNM) is then θW = S {r} L {r} ∞ j=2 P c (j|χ, {c}, {r}) • a j . ( 9
)
This form was found by [START_REF] Hellmann | Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals[END_REF] and [START_REF] Pérez-Enciso | Massive parallel sequencing in animal genetics: wherefroms and wheretos[END_REF][START_REF] Ferretti | Population genomics from pool sequencing[END_REF] in specific cases, but it holds for a very large class of data as shown. In the next section we will find the expression of P c (j| . . .) for the most common sequencing setups. Note that this estimator does not take into account sequencing errors. The treatment of sequencing error will be presented in the Discussion.

Complete sequences and genotyping data can be also analyzed by the above formulae, provided that an effective read depth r = +∞ is considered for positions with data, and r = 0 for positions with missing data. In fact, if there are missing data, the information is equivalent to the absence of NGS data (r = 0), while presence of data means that the genotype is known with certainty, as it would be with large read depth for NGS data (r = ∞).

Finally, the general Watterson estimator for an arbitrary model with spectrum E(ξ k |j) = θL ξk,j is

θW = S {r} L {r} ∞ j=2 P c (j|χ, {c}, {r}) j-1 i=1 ξi,j . (10)

A simple example: data from a single diploid individual

In this section we present the case of a single diploid individual [Lynch]. The single individual represents a single unit U = 1, and being diploid (assuming unrelated parents) we have two homologous sequences, therefore c 1 = 2. Assuming that each read is extracted at random from one of the two homologous sequences, the estimator is specified by the probabilities

P c (j = 1|c 1 = 2, r 1 ) = 2 -r 1 +1 (11) P c (j = 2|c 1 = 2, r 1 ) = 1 -2 -r 1 +1 (12)
Reference since the probability of extracting all r 1 reads from a given sequence is 2 -r 1 . Therefore the estimator from equation ( 9), taking into account the relevant harmonic factors a 1 = 0 and

A G A C C A T NGS reads: A C C T T G C C C A T A G C C Read depths: r 1 = 1 2 3 3 2 2 2
a 2 = 1, is θW = S ∞ r 1 =2 L r 1 (1 -2 -r 1 +1 ) . (13)
As an example, consider the sequences in Table 2. In the sequence of length L = 7 there are S = 2 segregating sites. There are bases with read depth 1, 2 and 3 and their numbers are L 1 = 1, L 2 = 4 and L 3 = 2 respectively. The value of the Watterson estimator for this sequence is therefore θW = 2/(1 • 0 + 4 • 1/2 + 2 • 3/4) = 4/7 0.57.

Distributions of the number of sequenced lineages

As discussed in the previous sections, the distribution of the number of sequenced lineages P c (j| . . .) is actually enough to define the Watterson estimator. Before deriving its expression for a number of cases, we introduce some notation.

We denote the Stirling numbers of second kind for j sets from r objects by S(r, j). We define the probability distribution

P * (j|c, r) = c! S(r, j) (c -j)! c r (14)
that corresponds to the probability of extracting exactly j different objects with r extractions (with repetitions) from a set of c objects [START_REF] Ferretti | Population genomics from pool sequencing[END_REF]. In fact, the number of possible extractions from a set of c objects is c r , while the number of extractions of precisely j objects is given by the product of the number of ordered choices of j objects out of c, that is c!/(cj)!, multiplied by the number of ways to distribute the j objects across r extractions S(r, j). Since all extractions are equiprobable, the ratio gives the the probability [START_REF] Kim | Detecting a local signature of genetic hitchhiking along a recombining chromosome[END_REF]. This equation will often appear in the formulae for the distribution of the number of lineages. We denote by I(x) the indicator function that takes the value 1 if x is true and 0 otherwise. We also denote by δ i,j the Kronecker delta, that is, the identity matrix δ i,j = I(i = j). Note that P * (j|c, 0) = δ j,0 .

General case: independent lineages

Assume that all lineages in these units are independent. This corresponds to sequencing many unrelated individuals in a population without inbreeding. If there are U units, c i is the number of lineages/homologous chromosomes in the ith unit, and r i is the number of reads/sequences coming from the ith unit, the probability P c (j|{c}, {r}) in the Watterson estimator is

P c (j|{c}, {r}) = c 1 i 1 =0 . . . c U i U =0 I j = U p=1 i p U q=1 P * (i q |c q , r q ) (15)
which is a product of probabilities for each unit of sequencing i 1 . . . i U chromosomes respectively, summed over all combinations resulting in j independent chromosomes.

In Section 5 we will present an alternative form for these Watterson estimators. In the rest of this section we specialize the expression [START_REF] Zhu | A composite-likelihood approach for detecting directional selection from dna sequence data[END_REF] to the most common scenarios.

Multiple haploid individuals

In this case all individuals have c i = 1, therefore P * (i q |c q = 1, r q ) = I(i q = I(r q > 0)) and the estimator reduces to the one proposed for missing data in [START_REF] Ferretti | Neutrality tests for sequences with missing data[END_REF], which is equivalent to

P c (j|c = n, {r}) = I j = n i=1 I(r i > 0) . ( 16
)
This choice was implicitly suggested also in [START_REF] Roychoudhury | Sufficiency of the number of segregating sites in the limit under finite-sites mutation[END_REF].

Multiple diploid individuals

In this case all individuals have c i = 2 and the MCLE was already derived by Hellmann et al. [START_REF] Hellmann | Population genetic analysis of shotgun assemblies of genomic sequences from multiple individuals[END_REF]:

P c (j|c = 2n, {r}) = 2 i 1 =0
. . .

2 in=0 I j = n p=1 i p n q=1 P * (i q |c q = 2, r q ) . (17)

Pools

In this case, there is a single unit of c chromosomes. The probability was derived by Ferretti et al. [START_REF] Ferretti | Population genomics from pool sequencing[END_REF] P c (j|c, r) = P * (j|c, r) [START_REF] Arnold | Extending coalescent theory to autotetraploids[END_REF] but see also Section 5 for a simpler formula. For multiple pools, the probability follows closely equation ( 15) where c 1 . . . c U are the numbers of (haploid) individuals inside each pool and r 1 . . . r U are the read depths of the pools at the position considered.

Related lineages

Sequencing unrelated individuals (either pooled together or sequenced separately) is the most common experimental setup for variability studies as described in the previous section, but it is not the only one. There are several cases where lineages in different units are related by identity (for example, the same individual sequenced both alone and in a pool with other individuals) or identity-by-descent, like for trios or inbred lines from a population. In this section we develop estimators for these cases.

Trios

A trio is a (diploid) family of mother, father and child that are sequenced separately. We assume that father and mother are two unrelated individuals from the same population. We restrict our analysis to autosomes, where the two alleles of the child are the copies of one paternal and one maternal allele. For complete sequences, the probability is just P (j) = δ j,4 since there are four independent lineages.

For NGS data, we denote by r M , r F and r C the read depths of mother, father and child respectively. We obtain P c (j|r M , r F , r C ) by conditioning on the number of lineages j sequenced in the parents. We can rewrite it in terms of the probability P c (j |c = 4, r M , r F ) of the parents alone (eq.( 17)):

P c (j|r M , r F , r S ) = 4 j =0 p t (j|j , r M , r F , r C )P c (j |c = 4, r M , r F ) ( 19
)
where p t (j|j , r M , r F , r C ) is the probability of sequencing j independent lineages in the trio given the number of independent lineages j sequenced from the parents. p t (j|j , r M , r F , r C ) is obtained case by case depending on the probability that the sequences of the child could contain new alleles with respect to the parental sequences and the probability to detect them. For example, consider the case j = 3. Then there is only a single allele in the parents that has not been sequenced. This allele is absent in the child with 50% probability (in this case j = 3 because no new alleles are sequenced in the child) or it could be present with 50% probability, but not sequenced (then j = 3 with probability 2 -rc ) or could be sequenced (then j = 4 with probability 1 -2 -rc ). The complete probability is

p t (j|j = 0, r M , r F , r C ) =P * (j|2, r S ) ( 20
)
p t (j|j = 1, r M , r F , r C ) = 1 2 P * (j -1|2, r S ) + 1 2 (δ j,1 2 -r S + δ j,2 (1 -2 -r S )) p t (j|j = 2, r M , r F , r C ) = 1 2 (1 + I(r M r F = 0))(δ j,2 2 -r S + δ j,3 (1 -2 -r S ))+ + 1 4 I(r M r F > 0) (δ j,2 + P * (j -2|2, r S )) p t (j|j = 3, r M , r F , r C ) = 1 2 δ j,3 + 1 2 (δ j,3 2 -r S + δ j,4 (1 -2 -r S )) p t (j|j = 4, r M , r F , r C ) =δ j,4 .
Multiple unrelated trios can be dealt with by replacing the probability P * (i q | . . .) in equation ( 15) with the probability [START_REF] Brenchley | Analysis of the bread wheat genome using whole-genome shotgun sequencing[END_REF] and replacing c q with 4.

Pooled trios

A pooled trio is a family of mother, father and child that are pooled together and sequenced. We consider n unrelated families, each family sequenced separately from the others, and denote by r i the total read depths.

P c (j|{r i } i=1...n ) is given by P c (j|{r i } i=1...n ) = 4 i 1 =0
. . .

4 in=0 I j = n p=1 i p n q=1 P pt (i q |r q ) ( 21
)
where P pt (i|r) is the probability of sequencing i homolog chromosomes for a single pooled trio. It can be derived case-by-case conditioning on the number i of sequenced chromosomes (identical or not) for a pool, obtaining

P pt (i|r) = i p pt (i|i )P * (i |c = 6, r) . ( 22
)
p pt (i|i ) can be found by conditioning on the number of non-inherited chromosomes sequenced, obtaining

p pt (i|i ) = δ i,i 4 2 i -2 + 4 2 i -1 + 2 i 6 i + δ i,i -1 4 2 i -3 + 2 2 i -2 6 i + δ i,i -2 2 i -4 6 i (23)
and finally

P pt (i|r) = 4 2 i -2 + 4 2 i -1 + 2 i P * (i|6, r) 6 i + (24) + 4 2 i -2 + 2 2 i -1 P * (i + 1|6, r) 6 i+1 + 2 i -2 P * (i + 2|6, r) 6 i+2
.

Pools and complete sequences with overlapping individuals

A potentially useful setup is the combination of complete sequences of few individuals and a pool of several individuals from the same population. In this situation, there could be individuals in the pool for which the complete sequence is also available.

Here we deal with the haploid case, but the results can be easily adapted to the diploid case by considering a diploid individual as a pair of haploids. Denote by m the number of individuals completely sequenced, by n the number of individuals pooled, and by o the overlapping between the two groups of individuals, i.e. the individuals in the pool that have also been sequenced separately. Denote by r the read depth of the pool. The distribution of j can be obtained by conditioning on the number l of pooled reads that come actually from the no individuals exclusive to the pool. The distribution of l is a binomial with probability o/n and r extractions, therefore

P c (j|r) = r l=0 P * (j -m|n -o, l) r l 1 - o n l o n r-l = (25) = (n -o)! n r (n -o -j + m)! r l=0 r l o r-l S(l, j -m) .
See also Section 5 for a simpler formula.

Inbred lines

Consider a population from which n inbred lines are derived. We denote the initial heterozygosity by H and the final heterozygosity by H inbred . The degree of inbreeding is measured by the inbreeding coefficient F = (H -H inbred )/H, that is the relative decrease in heterozygosity H due to inbreeding, and is assumed to be known. F is also equal to the probability of identity by descent for the inbred line. For each line, a diploid individual is sequenced. Our aim is to estimate the heterozygosity of the initial population from the sequences of individuals from the inbred lines.

If complete sequences are available, since each individual has an independent probability F of being homozygote because of inbreeding, the distribution of the number of homozygotes is just a binomial. But the number of sequenced chromosomes is 2n minus the number of homozygotes, therefore

P c (j|c = 2n) = n 2n -j F 2n-j (1 -F ) j-n . (26)
If instead we have NGS reads, the distribution of sequenced chromosomes should account for the " effective homozygote probability" F + (1 -F )2 -rq+1 due to sampling:

P c (j|c = 2n, {r}) = 2 i 1 =0
. . .

2 in=0 I j = n p=1 i p • n q=1 I(r q = 0)δ iq,0 + (27) +I(r q i q > 0) F + (1 -F )2 -rq+1 + (i q -1) 2(1 -F )(1 -2 -rq+1 ) -1 .
Note that in this case, the formula ( 6) with ( 26), ( 27) can be inverted to give the expected variability E(S) for individuals from inbred lines with a given inbreeding coefficient F .

An equivalent form for Watterson estimators

Equivalence between the estimator of Jiang et al. and the Watterson estimator for pools

An unbiased estimator of θ based on S was proposed in [START_REF] Jiang | Population genetic inference from resequencing data[END_REF] for NGS data of multiple diploid individuals, even if this is not the most appropriate setup, as we will see immediately. The estimator is

θJ = S ∞ r=2 L r c-1 k=1 1 k 1 -k c r -1 -k c r ( 28
)
where c is twice the sample size (for diploids) and r is the total read depth. This estimator is unbiased since the mean of S is given by the probability θ/k of a SNP of frequency k in the sample multiplied by the probability of detecting it in a random extraction of r alleles, that is 1

-k c r -1 -k c r .
A first observation is that this estimator is not actually unbiased for reads coming from multiple individuals sequenced separately. In fact, it takes into account only the total number of reads, while an unbiased estimator would depend on how they are distributed among individuals. However, it is an unbiased estimator of θ for pooled sequences, since in that case information about the origin of the reads is lost.

Furthermore, there is only a single unbiased estimator proportional to S, since the proportionality constant is fixed by the bias of S. This means that the estimator θJ is actually the Watterson estimator θW for pools proposed in [START_REF] Ferretti | Population genomics from pool sequencing[END_REF]. The two different forms derive from different intermediate conditioning for Z: on the allele frequency k in the sample in the first case, on the number of lineages actually sequenced j in the second.

Note that in the light of this equivalence, the conclusions of [START_REF] Jiang | Population genetic inference from resequencing data[END_REF] about the differences between their estimator and Hellmann's one when applied to individual data are at least doubtful. They found both estimators to be biased and the variance of Hellmann's one to be significantly larger, but theory suggests that they are unbiased and the variance of Hellmann's one should be lower. In fact, numerical simulations performed in [START_REF] Ferretti | Population genomics from pool sequencing[END_REF] showed almost no bias, no sensible difference in variance and a very good correlation between them.

From the mathematical point of view, the equality between θJ and θW for pools and the related equalities that we will present in the next section depend on a family of combinatorial identities. The identity of the two estimators θJ and θW for pools implies identity of their denominators. The reasoning in this section is equivalent to a double counting proof of the combinatorial identity min(c,r)

j=1 c! (c -j)! S(r, j)a j = c-1 k=1 c r -k r -(c -k) r k ( 29
)
valid for all pairs of integers (c, r) such that c ≥ 1 and r ≥ 1. The identity involves Stirling numbers S(r, j) and harmonic numbers a j in a nontrivial way. Note that both sides of the identity are integers. This identity can also be proved directly [M.Mamino, persona communication]. This identity is a combination of a family of related identities for a general spectrum, presented in Appendix A.

General alternative form for the Watterson estimators

The above form of [START_REF] Jiang | Population genetic inference from resequencing data[END_REF] for the Watterson estimator for pools can be generalized to the whole family of estimators for units of independent lineages, described by equations ( 9) and [START_REF] Zhu | A composite-likelihood approach for detecting directional selection from dna sequence data[END_REF]. We follow the same notation as before, but we denote the total number of lineages by c = U i=1 c i . The general form for these estimators is

θW = S {r} L {r} c-1 k=1 1 k Π k ({c}, {r}) , (30)
Π k ({c}, {r}) = c 1 k 1 =0 . . . c U k U =0 I k = U i=1 k i U i=1 c i k i c k 1 - U i=1 k i c i r i - U i=1 1 - k i c i r i
where the multi-hypergeometric distribution U i=1 c i k i / c k describes how the alleles are assigned to the different units and the term

U i=1 k i c i r i + U i=1 1 -k i c i r i
is the probability of extracting just one of the two alleles. All the estimators of section 4.1 can be rewritten in this form. This form is often more convenient computationally than the combinatorics in equations ( 9), [START_REF] Zhu | A composite-likelihood approach for detecting directional selection from dna sequence data[END_REF].

We can also find an estimator similar to [START_REF] Futschik | On the inadmissibility of Watterson's estimator[END_REF] for a combination of a pool of n (haploid) individuals and m complete sequences, o of which are overlapping. In this case simple combinatorial reasoning on the probability of detecting a SNP of frequency k among the n + mo individuals (the SNP is detected unless all complete sequences share the same allele) leads to the unbiased estimator

θW = S ∞ r=2 L r n+m-o-1 k=1 1 - ( n-o k ) ( n+m-o k ) 1 -k n r - ( n-o k-m ) ( n+m-o k ) k-m+o n r 1 k
(32) that is equivalent to the case (25) of estimator (9).

Watterson estimators for autopolyploids

A particularly interesting and challenging set of data is represented by polyploid genomes. Species with ploidy greater than 2 are highly interesting from an evolutionary point of view, as well as economically in agrobiotech and breeding since it involves many commercial species of plants (e.g. potato, sugar cane) and fishes (e.g. Salmonidae).

Polyploid species are difficult both to sequence and to analyze, due to the complex homology/paralogy relation between the constituent genomes. However, some polyploids can be treated by the methods developed here. In particular, autopolyploids are polyploid organisms with different homologous chromosomes from the same species. Autotetraploid populations follow the standard coalescent as shown by [START_REF] Arnold | Extending coalescent theory to autotetraploids[END_REF], and this can be extended to autopolyploids that have similar transition probability matrices. Here we present estimators of variability for populations of autopolyploid species.

Polyploids can be considered as pools with number of lineages equal to their ploidy. Multiple polyploids can then be considered as combinations of pools, but they can be pooled themselves. We consider a species with ploidy p. The estimators for autopolyploids are given by θW

= S ∞ r=2 L r p-1 k=1 1 k 1 -k p r -1 -k p r (33)
for a single polyploid individual, where r is the read depth, and by the formula

θW = S {r} L {r} np-1 k=1 1 k Π k (n, p, {r}) , (34)
Π k (n, p, {r}) = p k 1 =0 . . . p kn=0 I k = n i=1 k i n i=1 p k i np k 1 - n i=1 k i p r i - n i=1 1 - k i p r i
for sequences from n polyploid individuals, where {r} = {r i } i=1...n are the read depths per individual. This is also equivalent to the formula (9) for θW with

P c (j|{r}) = p i 1 =0
. . .

p in=0 I j = n l=1 i l n q=1 P * (i q |p, r q ) . (35)
The estimator for a pool of polyploid individuals is the same as in the general case for pools with c = np, where n is the number of individuals in the pooled sample:

θW = S ∞ r=2 L r np-1 k=1 1 k 1 -k np r -1 -k np r . (36)

Discussion

In this paper we have presented a large family of generalized Watterson estimators that are suited for different types of NGS data, from haploids to polyploids, pools and trios, or a mix of NGS/Sanger data. These estimators are built on the Maximum Composite Likelihood approach; furthermore they are unbiased and depend linearly on S, which is a sufficient statistic for small θ. The general theory presented here includes all these estimators and many others. Existing estimators are assigned to the proper place in this unified framework.

We pay special attention to estimators for single and multiple autopolyploid individuals. Sequencing of these species has proved to be hard, but more and more projects will soon be devoted to some of the more interesting polyploid species from a commercial point of view, especially among domesticated plants [START_REF] Brenchley | Analysis of the bread wheat genome using whole-genome shotgun sequencing[END_REF][START_REF] Han | Genome-wide snp discovery in tetraploid alfalfa using 454 sequencing and high resolution melting analysis[END_REF][START_REF] Wang | The genome of the mesopolyploid crop species brassica rapa[END_REF]. Autopolyploids without a strong inbreeding follow the dynamics of the usual coalescent, so our theory is applicable to these species. On the other hand, allopolyploids (whose genome derives from different species) cannot be studied by the same technique since the differences between homologous chromosomes from different constituent species are much stronger and the divergence time between them is often of order of the divergence between species. Specific methods have to be developed for the analysis of variability in allopolyploids [START_REF] Page | Insights into the evolution of cotton diploids and polyploids from whole-genome re-sequencing[END_REF][START_REF] Brenchley | Analysis of the bread wheat genome using whole-genome shotgun sequencing[END_REF]. A simple approach could be the study of the variability of each constituent genome and, independently, the genetic differentiation between them.

We did not discuss an important issue with NGS data, that is, base errors. Sequencing errors and misalignments occur at an high rate in NGS data. The bases with lower quality can be removed from the reads or the sequences, however sequencing errors or similar effects can often generate false SNPs at low frequency and it could be difficult to distinguish them from true low frequency alleles. In this case, filtering or SNP calling is usually applied to the data, resulting in an unknown Z ϕ,ξ for these alleles. Denote by (Φ, Ξ) the set of features {ϕ, ξ} strongly affected by sequencing errors. In some cases it is possible to estimate Z ϕ,ξ or to correct S ϕ,ξ based on quality scores for called SNPs or known error rates. On the other hand, if it is not possible to estimate the contribution of the errors, a good practice is to discard the corresponding S ϕ,ξ , {ϕ, ξ} ∈ (Φ, Ξ) and to work with the approximate MCL estimator for this case, that is θW = {ϕ,ξ} / ∈(Φ,Ξ) S ϕ,ξ {ϕ,ξ} / ∈(Φ,Ξ) L ϕ Z ϕ,ξ (37) as proposed in [START_REF] Achaz | Testing for neutrality in samples with sequencing errors[END_REF] for sequence data and [START_REF] Futschik | The next generation of molecular markers from massively parallel sequencing of pooled dna samples[END_REF], [START_REF] Ferretti | Population genomics from pool sequencing[END_REF] for pooled reads. The only alternative is to estimate Z ϕ,ξ by heuristic methods.

Generalizing equation (37), it is also possible to extend the results of this paper to generic sums of the frequency spectrum, for example estimators of the form k i=1 ξ k / k i=1 1/k which consider only the lowest frequencies. The estimator proposed here assume the standard Wright-Fisher neutral model for the allele frequency spectrum. However, an arbitrary expected frequency spectrum E(ξ k |n) = θL ξk,n could be used in the place of the neutral spectrum θL/k. It is sufficient to replace a j by j-1 i=1 ξi,j in the denominator of equation [START_REF] Ferretti | Population genomics from pool sequencing[END_REF] or to replace 1/k by ξk,c in the denominator of equation [START_REF] Arratia | Logarithmic Combinatorial Structures: A Probabilistic Approach[END_REF]. This extends previous adaptations of the original Watterson estimator to null scenarios with demography or varying population size (e.g. [START_REF] Živković | Second-order moments of segregating sites under variable population size[END_REF], [START_REF] Ferretti | Optimal neutrality tests based on the frequency spectrum[END_REF]).

In this paper we used the composite likelihood approximation to derive the estimators of variability. However, the variance of these Watterson estimators depends on recombination. The usual formulae for ML work only for unlinked sites. In this case, in the limit θ → 0 and θL constant, S ϕ,ξ is Poisson distributed, i.e. Var(S ϕ,ξ ) = E(S ϕ,ξ ) and therefore Var( θW ) = θ/ ϕ∈Φ L ϕ Z ϕ . In the same limit, the variance for linked sites contain a term θ 2 L 2 coming from the covariances between sites [START_REF] Fu | Statistical properties of segregating sites[END_REF][START_REF] Ferretti | Neutrality tests for sequences with missing data[END_REF][START_REF] Ferretti | Population genomics from pool sequencing[END_REF]. An exact formula for this term of the variance is available only for a few cases: complete sequences, sequences with missing data [START_REF] Ferretti | Neutrality tests for sequences with missing data[END_REF] and pooled NGS reads [START_REF] Ferretti | Population genomics from pool sequencing[END_REF]. In the case of completely linked sites and known variance, these estimators could be improved [START_REF] Fu | Estimating effective population size or mutation rate using the frequencies of mutations of various classes in a sample of dna sequences[END_REF], also by shrinkage methods [START_REF] Futschik | On the inadmissibility of Watterson's estimator[END_REF]. Note that the variance of the MCLE could be estimated by the bootstrapping methods described in [START_REF] Roychoudhury | Composite likelihood-based inferences on genetic data from dependent loci[END_REF].

Finally, the theoretical framework developed in this paper allowed to obtain an interesting set of combinatorial identities. This is another example of the way research on theoretical population genetics is highly connected to some fields of mathematics, e.g. combinatorics [START_REF] Arratia | Logarithmic Combinatorial Structures: A Probabilistic Approach[END_REF] and could lead to further mathematical insights.

for integers (r, c, l) with r ≥ 1 and 1 ≤ l ≤ c. They involve Stirling numbers and multinomials. Any other identity in this family (including (29)) can be obtained as a linear combinations of these ones. Note that since the l.h.s. is 0 for l > r, these identities reduce to

  For a generic frequency spectrum E(ξ k |n) = θL ξk,n , equation (30) should be replaced by θW = S {r} L {r} c-1 k=1 ξk,n Π k ({c}, {r}).

  k, cl (c rk r -(ck) r ) = 0 (A.2)for r < l ≤ c.

Table 2 :

 2 Example (not real) of data from a single diploid individual. The reads are aligned to the reference genome (the sequence at the top). Positions with SNPs are shown in bold red.

Some genotyping methods (like DNA microarrays) preselect the possible SNPs or the positions containing a segregating site. These methods give biased estimates of variability and cannot be meaningfully combined with unbiased methods like Sanger or NGS.

Acknowledgments

We thank A. Fonseca Amaral, M. Pérez-Enciso, W. Burgos, B. Nevado and G. Achaz for useful discussions, M. Mamino for providing an explicit proof of the main combinatorial identity and two anonymous referees for their constructive comments. LF acknowledges support from ANR-12-JSV7-0007 (ANR, France). The project was funded by Grants CGL2009-09346 (MICINN, Spain) and AGL2013-41834-R (MEC, Spain) to SERO and by a Consolider Grant from Spanish Ministry of Research, CSD2007-00036 "Centre for Research in Agrigenomics".

Appendix A. Combinatorial identities

We can extend the previous identity [START_REF] Roychoudhury | Composite likelihood-based inferences on genetic data from dependent loci[END_REF] to a set of identities derived from the same equivalence of estimators but for an arbitrary frequency spectrum. The fundamental identities are obtained by double counting technique.

We considering a frequency spectrum concentrated around a single frequency f in the population (i.e. ξ(f

for the sample spectrum). By double counting, the two Watterson estimators of the form ( 10) and (31) should be equal, and therefore we can equal their denominators. By computing the Taylor expansion in the variable f of both sides and equating the coefficients of the lth power, we obtain: