Taehee Kang

Jia Zhang

Achintya Kundu

Klaus Reimann

Michael Woerner

Thomas Elsaesser
email: elsasser@mbi-berlin.de

Bernard Gil

Guillaume Cassabois

Christos Flytzanis

Giorgia Fugallo

Michele Lazzeri

Ryan Page

Debdeep Jena

Ultrafast nonlinear phonon response of few-layer hexagonal boron nitride

Nonequilibrium dynamics of TO phonons and low-frequency interlayer shear modes are mapped in femtosecond mid-infrared pump-probe experiments. Time-resolved changes of TO phonon absorption reveal a TO phonon lifetime of 1.2 ps, while the low-frequency excitations decay with a time constant of approximately 20 ps. The coupling of intralayer TO and interlayer motions manifests in a quasi-instantaneous redshift of the TO phonon resonance by some 10 cm -1 . Theoretical calculations account for the coupling scenario and underline the relevance of inter-phonon interactions for the nonlinear phonon response of few-layer materials.

The crystal structure of hexagonal boron nitride (hBN) consists of covalently bound 6-rings of boron and nitrogen atoms arranged in two-dimensional layers [1]. Adjacent layers are coupled via much weaker van der Waals interactions, making the material highly anisotropic. Mono-and few-layer hBN have received strong interest as a prototypical quasi-twodimensional quantum material [START_REF] Shi | Synthesis of few-layer hexagonal boron nitride thin film by chemical vapor deposition[END_REF], holding strong potential for synthesizing van der Waals heterostructures [START_REF] Geim | Van der Waals heterostructures[END_REF] and applications in optoelectronics [START_REF] Caldwell | Photonics with hexagonal boron nitride[END_REF], nanophotonics, and phonon polaritonics [START_REF] Giles | Ultralow-loss polaritons in isotopically pure boron nitride[END_REF]. The electronic structure [START_REF] Cassabois | Hexagonal boron nitride is an indirect bandgap semiconductor[END_REF][START_REF] Schué | Bright luminescence from indirect and strongly bound excitons in h-BN[END_REF], the linear and the nonlinear optical properties [START_REF] Li | Probing symmetry properties of few-layer MoS 2 and h-BN by optical second-harmonic generation[END_REF][START_REF] Kim | Stacking order dependent second harmonic generation and topological defects in h-BN bilayers[END_REF] of few-layer hBN have been studied in some detail, whereas insight in carrier and phonon dynamics has remained very limited [START_REF] Golla | Ultrafast relaxation of hot phonons in graphene-hBN heterostructures[END_REF].

The anisotropic structure of hBN results in distinctly different frequencies of phonons connected with intra-and interplane atomic motions. Phonon modes of bulk hBN have been studied by infrared and Raman spectroscopy [START_REF] Geick | Normal modes in hexagonal boron nitride[END_REF][START_REF] Reich | Resonant Raman scattering in cubic and hexagonal boron nitride[END_REF][START_REF] Serrano | Vibrational properties of hexagonal boron nitride: Inelastic x-ray scattering and ab-initio calculations[END_REF][START_REF] Cuscó | Temperature dependence of Raman-active phonons and anharmonic interactions in layered hexagonal BN[END_REF][START_REF] Cuscó | Isotope effects on phonon anharmonicity in layered van der Waals crystals: Isotopically pure hexagonal boron nitride[END_REF], inelastic x-ray scattering [START_REF] Serrano | Vibrational properties of hexagonal boron nitride: Inelastic x-ray scattering and ab-initio calculations[END_REF], and theoretical calculations [START_REF] Michel | Phonon dispersions and piezoelectricity in bulk and multilayers of hexagonal boron nitride[END_REF][START_REF] Michel | Theory of rigid-plane phonon modes in layered crystals[END_REF]. At the Γ point (q = 0), the infrared-active transverseoptical (TO) phonon of E 1u symmetry displays a frequency of 1367 cm -1 (E T O = 169.5 meV), nearly identical to a Ramanactive phonon of E 2g symmetry. Both modes are connected with different in-plane atomic motions, as shown for the TO phonon in Fig. 1(a) (inset). Second-order Raman spectra suggest a TO phonon anharmonicity on the order of 5 cm -1 [START_REF] Serrano | Vibrational properties of hexagonal boron nitride: Inelastic x-ray scattering and ab-initio calculations[END_REF].

Lifetimes of the TO phonon of bulk hBN have been estimated from the damping term in numerical line shape fits of the reststrahlen band, giving values from 0.8 to 1.3 ps for samples with the natural isotope variation of boron [START_REF] Giles | Ultralow-loss polaritons in isotopically pure boron nitride[END_REF][START_REF] Cuscó | Isotope effects on phonon anharmonicity in layered van der Waals crystals: Isotopically pure hexagonal boron nitride[END_REF][START_REF] Segura | Long lifetime of the E 1u in-plane infrared-active modes of h-BN[END_REF]. An analysis of the Raman line shapes of the E 2g mode at 1367 cm -1 suggests a similar lifetime of 0.7 to 1 ps. For both modes, similar decay channels mediated via three-and four-phonon anharmonic couplings have been proposed [START_REF] Segura | Long lifetime of the E 1u in-plane infrared-active modes of h-BN[END_REF].

Motions of the stacked BN planes relative to each other give rise to rigid-plane shear and compression modes at low energy. In bulk hBN, the Raman-active horizontally-polarized rigid-plane shear mode arises at 54 cm -1 (E SH = 6.74 meV). In hBN with a finite number N of BN layers, the discretization of the k-vector perpendicular to the layers results in an increasing number of optical phonon branches. Theory [START_REF] Michel | Phonon dispersions and piezoelectricity in bulk and multilayers of hexagonal boron nitride[END_REF][START_REF] Michel | Theory of rigid-plane phonon modes in layered crystals[END_REF] predicts a splitting of the bulk rigid-plane modes in N -1 compression modes and N -1 doubly degenerate shear modes, the latter at energies below E SH . So far, experimental results on such modes are very limited [START_REF] Stenger | Low frequency Raman spectroscopy of few-atomic-layer thick hBN crystals[END_REF] and couplings of intra-and interlayer modes unknown.

For a basic physical understanding of the low-frequency dielectric properties of hBN and applications in phonon polaritonics, insight in phonon-phonon interactions and the resulting nonequilibrium dynamics is mandatory. Moreover, coupled phonon excitations are relevant for energy relaxation and heat transport, including the interplay of in-plane and interplane motions. In this Letter, we report the first study of the nonlinear phonon response of few-layer hBN. Upon femtosecond mid-infrared excitation via the TO phonon absorption band, the resulting ultrafast dynamics exhibit both a pronounced coupling of inter-layer and TO phonon modes and population relaxation processes of the coupled modes on different time scales. The anharmonic phonon-phonon coupling is accounted for by in-depth theoretical calculations.

The few-layer hBN sample studied here was grown directly on a sapphire substrate and has an average thickness of 2.5 nm, corresponding to N = 8 to 9 hBN layers. For such N, the rigid-plane shear phonons are expected at frequencies below the bulk shear mode at 54 cm -1 [START_REF] Michel | Theory of rigid-plane phonon modes in layered crystals[END_REF][START_REF] Stenger | Low frequency Raman spectroscopy of few-atomic-layer thick hBN crystals[END_REF]. A temporally and spectrally resolved pump-probe approach with femtosecond mid-infrared pulses is employed to map nonequilibrium phonon dynamics. As depicted schematically in Fig. 1(b) (inset), the s-polarized (TE) pump and probe pulses were reflected from a gold-coated prism and hit the sample under slightly different angles of incidence. The probe beam reflected from the sample passed the second facet of the prism, was spectrally dispersed in a monochromator, and detected with a 64-element HgCdTe detector array (spectral resolution 2 cm -1 ). A reference pulse reflected from an unexcited part of the sample was detected by a second array. Details of sample preparation, interaction geometry and mid-infrared pulse generation are presented in the supplemental material (SM).

The linear reflectivity spectrum [Fig. 1(a)] of the hBN sample exhibits the TO phonon resonance at 1366 cm -1 , in agreement with an attenuated total reflection spectrum of the sample (Fig. S2 of the SM). The width of the resonance has a value ∆ν T O ≈ 15 cm -1 , which is about twice as large as values derived from the reststrahlen band of high-quality bulk hBN. We assign this discrepancy to structural imperfections of the few-layer sample. A more detailed analysis of the reflectivity spectrum is presented in the SM.

In Experiments with pump pulses of a reduced bandwidth ∆ν p give additional insight in the nonlinear phonon response. The pump energy E p was adjusted to implement a similar excitation level for ∆ν p = 20, 32, 54, and 200 cm -1 (FWHM). Pump-probe spectra for a delay time t D = 3 ps after the fast signal decay [Fig. 3(a)] display identical envelopes for the different ∆ν p . The time evolution of reflectivity changes at a fixed probe frequency (not shown) is independent from ∆ν p as well. In Fig. 3(b), the frequency positions of the positive and negative peaks in the spectra for t D = 3 ps are plotted vs ∆ν p . Figure 3(c) shows the frequency position of the positive peak (solid symbols) and the first moment M = ( νg(ν)dν)/ g(ν)dν of the redshifted positive component g(ν) of the pump-probe spectra as a function of E p (open symbols). Both quantities exhibit a redshift increasing roughly linearly with E p .

We now discuss the mechanisms behind the ultrafast nonlinear phonon response with the support of first principles evaluations [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF] of phonon lifetimes, anharmonic frequency shifts and inter-phonon anharmonic couplings [26, 32, 33]. The theoretical methods are accounted for in the SM. We recall that the TO phonon potential energy surface shows a small anharmonicity along the q T O coordinate of 5 cm -1 [START_REF] Serrano | Vibrational properties of hexagonal boron nitride: Inelastic x-ray scattering and ab-initio calculations[END_REF].

As schematically shown in Fig. 3(d), this anharmonicity leads to a redshift of the v = 1 to 2 compared to the fundamental v = 0 to 1 transition. Resonant excitation by the pump pulses [blue arrow in Fig. 3(d)] results in a population decrease in the v = 0 and an increase in the v = 1 TO phonon state. Concomitantly, the absorption/reflectivity on the v = 0 to 1 transition decreases, due to depopulation of the v = 0 state and stimulated emission from the v = 1 state. In parallel, a redshifted absorption/reflectivity increase arises on the v = 1 to 2 transition. The probe pulse [red arrows in Fig. 3(d)] maps such changes.

The dispersive shape of the transient spectra in Figs. 1(b,c) is in line with this picture. For a pump energy E p = 70 nJ, all reflectivity changes occur within the steady-state TO phonon peak [thick solid line in Fig. 1(b)], reflecting the small diagonal anharmonicity. The time resolved transients in Fig. 2(a) display a first decay time of 1.2 ps, which we assign to the population decay of the v = 1 state of the TO phonon. The 1.2-ps decay measured here directly supports indirect lifetime estimates from vibrational linewidths [START_REF] Giles | Ultralow-loss polaritons in isotopically pure boron nitride[END_REF][START_REF] Segura | Long lifetime of the E 1u in-plane infrared-active modes of h-BN[END_REF]. For E p = 1 µJ, one observes a similar 1.2-ps decay at ν pr = 1342 and 1366 cm -1 .

Relaxation of the v = 1 TO phonon population is connected with a population transfer to anharmonically coupled phonons at lower frequency [START_REF] Cuscó | Temperature dependence of Raman-active phonons and anharmonic interactions in layered hexagonal BN[END_REF][START_REF] Cuscó | Isotope effects on phonon anharmonicity in layered van der Waals crystals: Isotopically pure hexagonal boron nitride[END_REF][START_REF] Segura | Long lifetime of the E 1u in-plane infrared-active modes of h-BN[END_REF]. The theoretical calculations of Ref. 18 suggest a predominant 4-phonon decay channel, invoking the zone-center A 2u phonon at 820 cm -1 and ZO/ZA phonons at the M, K points of the Brillouin zone with frequencies around 300 cm -1 . A relaxation scheme and more detailed discussion based on our calculations is presented in the SM. The calculated TO-phonon lifetimes of 1 ps [START_REF] Segura | Long lifetime of the E 1u in-plane infrared-active modes of h-BN[END_REF] and 1.78 ps (SM) are close to the experimental value of 1.2 ps.

The pump-probe spectra (Fig. 1) display a redshift of the fundamental (v = 0 to 1) TO phonon absorption even after the decay of the TO v = 1 population, i.e., with the TO phonons in their v = 0 state. The non-exponential time-resolved transients measured with E p = 1 µJ at ν pr = 1349 and 1357 cm -1 [Fig. 2(b)] demonstrate the onset of this redshift at early delays during the decay of the v = 1 state, and, together with all other transients in Fig. 2, their decay with a time constant of 22 ps. The redshift of the TO phonon transition is due to the intermode anharmonic coupling to a subset of low-frequency phonons excited by the pump pulse. As sketched in Fig. 3(e), a Raman process which is resonantly enhanced by the TO transition dipole, excites low-frequency Raman active phonons q SH,i . This excitation mechanism is nonselective and generates excess populations of all Raman-active modes within the pump bandwidth. The nonequilibrium populations of lowfrequency modes redshift the TO phonon transition via the anharmonic coupling V anh . The redshift at a delay t D = 3 ps after the TO phonon decay scales with the pump energy E p [Fig. 3(c)], i.e., the generated excess population, and decays with a time constant of 22 ps by population relaxation of the low-frequency phonon states.

The experiments with reduced pump bandwidth ∆ν p (Fig. 3) show that the redshift at t D = 3 ps scales with E p but not with ∆ν p . The data for ∆ν p = 20 cm -1 (FWHM), however, set an upper limit of 30-40 cm -1 for the frequency of the phonons inducing the redshift. While optically accessible phonons below 50 cm -1 are absent in bulk hBN, a subset of Raman active interlayer shear modes i and 1-2 compression modes of few-layer hBN shows sub-40 cm -1 frequencies [cf. SM, [START_REF] Michel | Theory of rigid-plane phonon modes in layered crystals[END_REF]]. The observed TO redshift is assigned to a coupling with such low-frequency phonons in which nonequilibrium populations are generated by the Raman pump mechanism. Independent evidence for the exis-tence of low-frequency shear and compression modes comes from the Raman spectra of Ref. [START_REF] Stenger | Low frequency Raman spectroscopy of few-atomic-layer thick hBN crystals[END_REF], displaying a broad lowfrequency phonon signal superimposed by a sharp line due to the strongest Raman line at the highest shear-mode frequency [START_REF] Liang | Interlayer bond polarizability model for stacking-dependent lowfrequency Raman scattering in layered materials[END_REF].

This picture relies on the concept that due to optical selection rule the pump pulse excites only phonons with zero wavevector, i.e., in an infinitesimally small volume in the Brillouin zone (BZ). If we indicate the phonon population as n(q), where q is the wavevector, immediately after the pump excitation we have n(q) = n 0 (q) + Aδ (q) where n 0 (Q) is the thermal (T = 300 K) population and δ is the Dirac distribution [START_REF]A is basically a phonon population resulting from a Bose-Einstein[END_REF]. Such a distribution does not alter the phonon lifetimes, which are determined by the population distribution over all the BZ (and not by what happens only at q = 0). However, it can induce a frequency shift of the TO mode through anharmonic coupling with the low-energy shear modes at the BZ center. Interestingly, in bulk hBN the TO and the lowenergy shear modes are weakly coupled, with a predominant contribution of 4-phonon processes and a coupling term of the order of -0.02 cm -1 . With this weak coupling in order to experimentally observe a TO shift ∆ν T O = 10 cm -1 , we should assume A ≈ ∆ν T O /(0.02 cm -1 ) ∼ 500. For the bulk E SH = 6.7 meV, this value of A corresponds to an effective energy per unit cell of ∼ 1 eV to be provided by the pump pulse, a value prohibitively high under the present experimental conditions.

On the contrary, the generation of a TO shift ∆ν T O = 10 cm -1 becomes possible when passing from bulk to few-layer hBN. Our calculations show that the coupling term displays a strong dispersion along the z-axis and increases up to two orders of magnitude. For a maximum shear-mode frequency of 30 cm -1 (cf. Fig. 3), such strong coupling shrinks the excitation energy required per unit cell to roughly 20 meV, well in the range of the experiment. Our picture is further supported by the calculated lifetime of the low-energy shear modes of τ SH = 25.81 ps, close to the experimental decay time of 22 ps. The large difference between τ SH and the 1.2-ps TO phonon lifetime assures, as experimentally observed, that the TO nonlinear signal decreases while the redshift remains the same during the TO decay time.

In conclusion, the ultrafast nonlinear phonon response of few-layer hexagonal boron nitride reveals a TO phonon lifetime of 1.2 ps and a 22-ps decay time of low-frequency phonons. Nonequilibrium populations of rigid-plane shear modes below 30 cm -1 induce a transient redshift of the TO phonon resonance, a direct manifestation of the anharmonic coupling between the TO-intralayer and interlayer degrees of freedom. This coupling is highly relevant for future extensions of phonon polaritonics into the nonlinear regime.

FIG. 2 .

 2 FIG. 1. (a)Reflecticity spectrum of the few-layer hBN sample in the range of the TO phonon resonance at 1366 cm -1 (symbols). Dash-dotted line: experimental reflectivity spectrum of the sapphire substrate. Inset: Atomic elongations connected with the TO phonon of hBN. (b) Transient reflectivity spectra (symbols) measured with a pump energy E p = 70 nJ. The change of reflectivity ∆R/R 0 = (R -R 0 )/R 0 is plotted as a function of probe frequency for different delay times (R, R 0 : reflectivity of the sample with and without excitation). Solid lines are to guide the eye. The thick solid line gives the inverted and scaled reflectivity spectrum from panel (a) after subtraction of the substrate reflectivity. Inset: Pump-probe geometry applied in the experiments. (c) Same as (b) for a pump energy E p = 1.0 µJ. The arrows mark the probe frequencies at which the transients in Fig.2were recorded.

FIG. 3 .

 3 FIG.3. Pump-probe spectra at a delay t D = 3 ps recorded with pump bandwidths ∆ν p . The spectra are normalized to the maximum tive ∆R/R 0 . The pump energy was E p = 400 nJ for ∆ν p = 200 cm -1 . For a reduced ∆ν p , E p was adjusted to make the amount of absorbed pump energy identical to the data for E p = 400 nJ. (b) Frequency positions of the maximum negative (1) and positive (2) reflectivity changes as a function of ∆ν p . (c) Peak position and first moment of the positive reflectivity change as a function of E p for excitation with ∆ν p = 200 cm -1 (d) Schematic of the anharmonic TO phonon potential along the coordinate q T O with optical transitions relevant for excitation (blue arrow) and probing (red arrows). (e) Anharmonic TO phonon potential along q T O and vibrational levels of low-frequency shear modes i with a Raman excitation process (blue arrows). The TO phonon and shear modes are anharmonically coupled via V anh .


Nishant Thakur

Nathalie Pujol

Jacques Van Helden

Robert H Waterston

Ladeana W Hillier

Laurent Tichit

Jonathan J Ewbank
email: ewbank@ciml.univ-mrs.fr

An integrated analysis tool reveals intrinsic biases in gene set enrichment

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In recent years, there has been a boom in genomic, transcriptomic and epigenomic studies, largely fuelled by advances in sequencing technologies and the attendant reduction in costs. They often result in the production of long lists of candidate genes.

Various publicly-available resources classify genes on the basis of their structure, interactions or function. A common way to interpret a new gene list is to exploit this available knowledge, looking for enrichment of genes in defined classes. Despite inherent limitations and the tendency for results to be moulded to fit "a sensible biological narrative" [START_REF] Pavlidis | A critical assessment of storytelling: gene ontology categories and the importance of validating genomic scans[END_REF], particularly in hypothesis-driven research [START_REF] Yanai | A hypothesis is a liability[END_REF], this technique of gene functional enrichment has been used for more than a decade [START_REF] Huang Da | Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists[END_REF] and there are currently dozens of available tool, many listed at http://omictools.com/.

There are, however, concerns about these tools, even for the most highly cited ones. A major problem with many of them is that they are not regularly updated. Among the 21 most popular tools in 2016, 12 had not been updated since 2011; 84% of citations to gene enrichment tools refer to those that are outdated (Wadi, et al., 2016). For example, the data in DAVID was not updated for 7 years, between version 6.7, made public in 2009, and version 6.8, released in 2016, and according to the DAVID website has not been updated since. Despite this, the tool is still being cited several thousand times each year. This is a problem for 2 reasons. Firstly, gene predictions change. For the model organism Caenorhabditis elegans, between the WS220 release of the database Wormbase in [START_REF] Sleumer | Caenorhabditis elegans cisRED: a catalogue of conserved genomic elements[END_REF]WS255 (released in 2016), for instance, more than 1700 new genes were defined, and the predicted structure of hundreds of existing genes modified.

If tools do not update their source data, an increasing number of genes may be either incorrectly associated with an annotation, or simply absent. Secondly, these tools lack up-to-date information about gene function. Recently, gene ontology (GO) annotations have increased on average by 12.5% every year (Wadi, et al., 2016). There is a similar trend for other annotations, like those from Reactome [START_REF] Fabregat | The Reactome pathway Knowledgebase[END_REF] and KEGG [START_REF] Kanehisa | KEGG for representation and analysis of molecular networks involving diseases and drugs[END_REF]). Thus, for example, fully 80% of functional annotations available in mid-2016 were absent from DAVID 6.7. This is clearly a severe limitation, and can substantially bias analyses of gene lists (Wadi, et al., 2016). Some tools that perform gene enrichment analyses are regularly updated. One example is g:Profiler [START_REF] Reimand | Profiler-a web server for functional interpretation of gene lists (2016 update)[END_REF]. Despite its power, in common with similar tools, g:Profiler captures annotations for multiple species from centralized databases, rather than leveraging the annotations available in species-specific databases. The GO data for C. elegans, for example, is derived from GO Consortium's compendium, so that it almost always lags behind Wormbase releases. Further, Wormbase contains more-orless detailed descriptions of the phenotypes associated with mutation or RNAiknockdown for ~40% of genes. Although much of this information in time reaches generic databases, there is always a delay, so at a given moment, for some annotations, Wormbase is the sole source. Wormbase also harbours more than 1600 transcriptome expression datasets that are manually curated to a standard machinereadable format, allowing automated retrieval and data integration. This represents a powerful resource for functional enrichment analysis that is more readily exploited than the raw data available through the main transcriptome databases such as the Gene Expression Omnibus (GEO).

Thus while generic databases tools have a broad appeal, there is a clear demand for model organism databases (MODs) [START_REF] Oliver | Model organism databases: essential resources that need the support of both funders and users[END_REF] and tools that are speciesspecific. WormExp is one of many C. elegans-specific bioinformatic tools. It overcomes some of the limitations of the more generic ones, but is deliberately limited in its scope of data, principally from transcriptome studies [START_REF] Yang | WormExp: a web-based application for a Caenorhabditis elegans-specific gene expression enrichment analysis[END_REF]. Others include a database of time-resolved expression data [START_REF] Grun | Conservation of mRNA and protein expression during development of C. elegans[END_REF][START_REF] Stoeckius | Global characterization of the oocyte-to-embryo transition in Caenorhabditis elegans uncovers a novel mRNA clearance mechanism[END_REF], catalogues of potential transcription factor binding sites, established on the basis of the DNA conservation (cisRED [START_REF] Sleumer | Caenorhabditis elegans cisRED: a catalogue of conserved genomic elements[END_REF] or ChiPseq experiments (motif-disc [START_REF] Araya | Regulatory analysis of the C. elegans genome with spatiotemporal resolution[END_REF]), tools such as Worm-Cat, GExplore and WormMine for large-scale data mining related to gene or protein function [START_REF] Ding | Stress-responsive and metabolic gene regulation are altered in low Sadenosylmethionine[END_REF][START_REF] Harris | new views of curated biology[END_REF][START_REF] Hutter | GExplore: a web server for integrated queries of protein domains, gene expression and mutant phenotypes[END_REF], Tissue Enrichment Analysis among gene sets [START_REF]Tissue enrichment analysis for C. elegans genomics[END_REF], WormNet that can generate new members for a pathway, infer functions from network neighbours, or predict "hub" genes, on the basis of annotations and expression data [START_REF] Cho | WormNet v3: a network-assisted hypothesis-generating server for Caenorhabditis elegans[END_REF] and GeneModules to identify active "modules" in gene expression data [START_REF] Cary | Application of Transcriptional Gene Modules to Analysis of Caenorhabditis elegans[END_REF]. On the other hand, to the best of our knowledge, there is no single tool for enrichment analysis that exploits the diverse sets of data available for C. elegans. As reported here, we created a tool YAAT (for "yet another analysis tool") that performs enrichment analysis using a very extensive dataset (> 10000 classes) gathered from multiple sources. To increase its utility, we supplemented it with several complementary analytical methods, enriched class clustering, evaluation of conservation, phylogenetic and co-expression profiling, as well as tools for visual exploration of links between classes.

In the course of benchmarking to explore the robustness of test results, we revealed biases that have the potential to skew the results of broad-category enrichment analyses. Our results show how caution is needed when interpreting any result from gene enrichment analyses in C. elegans. These concerns are likely to be true across species.

SYSTEM AND METHODS

Data collection

To constitute the database for functional enrichment analysis, we expanded our previous collection of expression and phenotypic data from diverse resources [START_REF] Zugasti | A quantitative genome-wide RNAi screen in C. elegans for antifungal innate immunity genes[END_REF]. Firstly, for expression data we extracted ~1600 Serial Pattern of Expression Levels Locator (SPELL; [START_REF] Hibbs | Exploring the functional landscape of gene expression: directed search of large microarray compendia[END_REF]) expression clusters, curated from several hundred articles, from the Wormbase FTP site ftp://caltech.wormbase.org/pub/wormbase/spell_download/tables/. We supplemented this data with non-redundant datasets from WormExp [START_REF] Yang | WormExp: a web-based application for a Caenorhabditis elegans-specific gene expression enrichment analysis[END_REF]. As neither Wormbase expression clusters nor WormExp are comprehensive, we manually curated a further 188 datasets from 71 articles, building on our previously described in-house set derived from 69 articles [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF]. Most of the phenotypic data was extracted from Wormbase. Thus 1907 phenotypic classes, associated with a total of 8077 genes, were downloaded from Wormbase release WS255

(ftp://ftp.wormbase.org/pub/wormbase/releases/WS255/ONTOLOGY/phenotype_association.W S255.wb). Further, 101 phenotypic classes were extracted from DRSC FlyRNAi database (http://www.flyrnai.org/RNAi_all_hits.txt) and the correspondence between fly genes and their nematode homologues established as previously described [START_REF] Zugasti | A quantitative genome-wide RNAi screen in C. elegans for antifungal innate immunity genes[END_REF]. We also automatically extracted the list of genes associated with each of the 132 C. elegans pathways present in KEGG. Finally, we added regulatory information into YAAT by including 202 datasets for the putative gene targets for 94 transcription factors from the last modENCODE release [START_REF] Araya | Regulatory analysis of the C. elegans genome with spatiotemporal resolution[END_REF]. In addition, we pulled into YAAT the gene associations for 5668 GO terms.

Maintaining data up to date

The database for functional enrichment analysis includes categories of "static" data from sources that are not updated, such as from the DRSC FlyRNAi database, modENCODE and any transcriptome data referenced to gene not a specific sequence, and "dynamic" data, for example RNAi phenotypes from Wormbase that are referenced to specific genomic coordinates. The static data was kept up to date by tracking changes in gene structure prediction for each set, using an R script based on the algorithms underlying Wormbase Converter [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF]. The different sources of dynamic data have their particular release schedules. Our aim was to match YAAT updates to the Wormbase release cycle, with an interval of 5 releases (http://www.wormbase.org/about/release_schedule). We were successful in completing one update cycle and therefore data were generated for WS255 and WS260. The latter is a very extensive resource of 10395 datasets.

Enrichment analysis

For functional enrichment statistics, we used the hypergeometric distribution using the R "stats" package. The p-value is calculated as

phyper(q = x-1, m = m, n = n, k = k, lower.tail = FALSE)
Where k is the number of genes in the user's query.

N is total number of genes annotated in the catalogue of reference classes (e.g. 20129

genes for the analyses of data from WS255 reported here).

m is the number of genes annotated in the functional class of interest.

x is number of genes in the user list that are annotated in the functional class of interest.

n= N -m

For multiple comparisons, p-values are corrected with the p.adjust function in the R "stats" package, with various types of correctional statistics: "holm" [START_REF] Holm | A Simple Sequentially Rejective Multiple Test Procedure[END_REF], "hochberg" [START_REF] Hochberg | A sharper Bonferroni procedure for multiple tests of significance[END_REF], "hommel" [START_REF] Hommel | A stagewise rejective multiple test procedure based on a modified Bonferroni test[END_REF], "BH" or its alias "fdr" and "BY" [START_REF] Benjamini | The Control of the False Discovery Rate in Multiple Testing under Dependency[END_REF].

Clustering of enriched classes.

Having defined the enriched classes for the user's list of genes, a binary matrix of enriched classes and the overlapping genes in each class is generated, where rows represent genes found in at least one enriched class and columns represents enriched classes. Genes common to the user list and the enriched class are set to be "1", otherwise genes are assigned the value "0". The distance between all the genes and all the enriched classes are calculated and based on the distances, hierarchical clustering is performed using "hclust" function in R.

Phylogenetic profiles

For phylogenetic profile analysis, using WS255 we made a reference set of 101 Genome release 31). The complete predicted protein sets for a further 57 metazoan species, chosen for their phylogenetic diversity, but including all 5 available Caenorhabditis species, were downloaded from Ensembl (http://metazoa.ensembl.org; release 32). Each C. elegans protein was compared to the complete set of proteins in each of the helminths or diverse proteomes (including C. elegans) using BLASTP.

When multiple isoforms existed for a single protein, the longest one was selected. For each C. elegans query protein, only the best hit in each species was considered. Based on this analysis, for each of the 2 datasets, we generated a "BinMat" matrix (n proteins × m species), where for each entry C ij , if there was a BLASTP hit in species "j" for C. elegans protein "i" with a P-value less than 10 -5 , then C ij was set to "1", and otherwise to "0". Users provide lists of proteins or genes. In the latter case, for each gene, the corresponding protein (the longest isoform, if applicable) is identified. Starting with this list of "N" genes, the tool will extract a "QuBinMat" matrix (N proteins × m species) from "BinMat". Next, the distances between the "N" proteins within the QuBinMat matrix are calculated and based on these distances, clustering is performed with the "hclust"

clustering algorithm in R.

To calculate an average conservation score across the N species for a given list of n genes, first, for each gene, the number of species predicted to contain an ortholog, as defined above, was determined (o 1 , o 2 ,..o n ). Next, the sum of the number of orthologues was divided by the maximum possible score (∑ 1-n o /Nn).

Co-expression dataset

A coherent set of the results from 240 independent RNAseq experiments performed by the modENCODE consortium was assembled. In addition to transcriptome profiling of worms at different developmental stages, the set included data from different mutant strains and worms exposed to a variety of pathogens. The raw sequencing reads for each experiment were treated in an identical fashion, with expression level reported as average depth of coverage per base per million reads (dcpm). The data, at a gene level in WS220, was converted to WS255 using WormBase Converter. Pearson correlations are calculated for the genes of interest within the dataset, the calculated correlations clustered as above and plotted using the function heatmap.2 within the R package ggplot.

Fungal infection dataset

To define a test set of genes strongly induced by D. coniospora ("Dc_Up" genes), using published data (Table S7 in [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF])), we selected those with a >2.5fold increase in infected versus control worms and with a minimum expression level of 0.1 dcpm in infected worms. The list of 145 genes was converted from WS170 to WS260 using WormBase Converter [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF], giving a list of 146 genes, including 2 that are non-coding.

IMPLEMENTATION Constitution of reference gene sets for YAAT

In order to analyse a large dataset of genes derived from a genome-wide RNAi screen, we developed a functional class enrichment and clustering tool [START_REF] Zugasti | A quantitative genome-wide RNAi screen in C. elegans for antifungal innate immunity genes[END_REF].

Given the interest generated by this tool, we decided to improve it and adapt it with the aim of providing a web-based resource, YAAT. The first improvement was to expand the range of underlying data. Our primary data source was Wormbase. From there, we extracted, for example, the genes associated with each of ca. 2000 phenotypic classes.

A single gene can give rise to multiple transcripts, and hence to different protein isoforms. Some functional annotations are associated with a specific transcript, but we did not attempt to capture this level of detail in YAAT. Thus the reference object within YAAT is a gene (with unique Wormbase gene identifier; WBGene). Indeed, for annotations associated with proteins (e.g. for enzyme activities or derived from proteomic studies), the corresponding gene name was used. Thus for the sake of simplicity, here we use the term "gene" in an indiscriminate way when describing classes that were defined on the gene, transcript or protein level.

Wormbase has a regular release schedule, with updates available in principal every 2 months. Some data sources for YAAT are static. For example, the association of transcription factors and putative target genes established by the modENCODE consortium was generated using WS220 [START_REF] Araya | Regulatory analysis of the C. elegans genome with spatiotemporal resolution[END_REF] and although this largescale study has been extended in the context of modERN [START_REF] Kudron | The ModERN Resource: Genome-Wide Binding Profiles for Hundreds of Drosophila and Caenorhabditis elegans Transcription Factors[END_REF], the modENCODE mapping of transcription factor targets has not been updated, despite subsequent changes in gene predictions. Other data sources, including KEGG [START_REF] Tanabe | Using the KEGG database resource[END_REF], are updated regularly. We previously developed a tool, Wormbase Converter, to deal specifically with this issue [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF]. We used a script based on the Wormbase Converter algorithm to ensure that we had a homogeneous set of data, using a single reference release of Wormbase genes, irrespective of the data's source.

Although we successfully performed one semi-automatic update of the YAAT database (from WS255 to WS260), with the passage of time, this became unfeasible for a number of reasons. One barrier to updating YAAT was changes in the location and availability of third party data. Common problems included out-dated URLs (e.g. the broken link to "a tab-delimited file of all public 'hits' (positive results) from DRSC screens" from https://fgr.hms.harvard.edu/utility-tools) and re-organisation of ftp data repositories. The use of database-specific gene identifiers (e.g. in KEGG) also presented a challenge as the inevitable occasional parsing errors required manual correction. Further, the release cycles of the different databases are not synchronised with Wormbase updates. For example, the Ensembl database is a variable number of releases behind Wormbase.

This was particularly problematic for the establishment of phylogenetic profiles (see below) and could only be resolved by fastidious and very time-consuming crossvalidation and error checking. Thus although YAAT was designed at the outset to be readily updated, we failed in that aim and have not yet been able to implement a robust pipeline for semi-automatic updating. Nevertheless, compared to many other tools,

YAAT is based on relatively recent data and provides functions that they cannot.

Two major sources of transcriptome data for YAAT were the pre-computed expression clusters from the C. elegans implementation of Serial Pattern of Expression Levels

Locator (SPELL [START_REF] Hibbs | Exploring the functional landscape of gene expression: directed search of large microarray compendia[END_REF]) and WormExp [START_REF] Yang | WormExp: a web-based application for a Caenorhabditis elegans-specific gene expression enrichment analysis[END_REF]. SPELL includes only published datasets, while the underlying data in WormExp comes from the many public transcriptome databases and includes unpublished data. These sources thus provided a partially redundant coverage of transcriptome data. We also expanded our pre-existing in-house functional class database [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF] by manually extracting data from transcriptome studies. Since this dataset only references gene identities, not expression values, it includes some studies for which the expression values were never made publicly available and that are therefore absent from SPELL and WormExp (Figure 1A). This transcriptome information was combined with functional and structural categories from different sources to constitute the YAAT database (Figure 1B).

This set covers essentially all (>99%) protein-coding genes in C. elegans. The average class size was 272, with a median of 28 and a maximum of 10,705 genes (Figure 1C).

Classes from different sources had markedly different size distributions. Thus the 1907 classes derived from Wormbase phenotypes were biased towards small classes (median 6, average 34), while those derived from SPELL (1662 classes) were biased towards large classes (median 104, average 484; Figure 1C).

Less than 130 genes (0.6%) were associated with a single category, in most cases (67/127 genes), the class "C. elegans-specific genes" [START_REF] Zhou | Genome-wide identification of lineage-specific genes within Caenorhabditis elegans[END_REF]. On the other hand, 95% of genes were associated with at least 10 classes, most frequently 21, with the median number of classes for each gene being 58 (Figure 1D). The complete dataset therefore has good gene coverage and the potential to provide insight into the functional relationships between the members of most gene lists.

Gene sets analysis using YAAT

In its most basic implementation, a user inputs a list of genes of interest and sets analysis parameters. Users can choose between different methods to determine the statistical significance of each gene set enrichment; the default is by false discovery rate (FDR). Thresholds for FDR or P-values and for the maximum reference class size can be defined. In the context of enrichment analysis, including gene classes that contain as many as 50% of all protein-coding genes makes relatively little sense; we generally, and arbitrarily set this parameter to 1500 genes (i.e. ca. 7.5% of all protein-coding genes).

Enriched functional classes are returned in a hyperlinked table that gives access to the underlying data. The results are also available as a downloadable text file and, as detailed further below, in different graphical formats.

Analysis of the composition of the YAAT dataset reveals intrinsic biases

To characterise the dataset and the tool, we took the list of genes for each and every gene class and input this in turn into YAAT, to investigate the relationship between class sizes and number of enriched gene classes. We observed a positive correlatory trend between the number of genes in a class and the number of enriched functional classes for the entire set, and also for subsets corresponding to the Wormbase phenotype, WormExp or SPELL classes (Figure 2). There were, however, numerous classes that did not follow this trend. There were, for example, 4 large SPELL classes (of >5000 genes) associated with an unexpectedly low number of enriched categories (≤120, compared to a median of 416). It is not clear whether this reflects a technical artefact, but they come from a single study that transcriptionally profiled different cell types [START_REF] Spencer | A spatial and temporal map of C. elegans gene expression[END_REF]. At the other end of the scale, the class "Down <-1.5 by tunicamycin in N2", with 126 genes returned more than 400 enriched classes (Table S1).

During the analysis, we noticed an apparent bias in the categories of classes returned for different gene lists. When the gene list used as a query was derived from a SPELL class, overall there appeared to be an over-representation of classes comprised of transcriptionally regulated genes (i.e. SPELL or WormExp). This was not seen with gene lists derived from WBPhenotype classes (Table S1). To probe this further, and to avoid the confounding effect of class size, we took the results for the 22 WBPhenotype classes constituted of 200-300 genes, and selected the 25 SPELL classes (only one per published study, chosen at random) of the same size that gave at least 300 enriched classes when run through YAAT (Figure 3A). Within these groups of gene lists, there was a striking difference in distribution of the enriched gene classes. While groups of genes defined by transcriptome studies (i.e. SPELL classes) were found to be enriched almost exclusively for other SPELL or WormExp classes, the groups defined on the basis of a phenotype were enriched for a far more diverse range of categories (Figure 3B), in proportions that approached their overall distribution (Figure 1B). Part of the reason for this result, as discussed further below, is that the respective gene "universes", i.e. the non-redundant list of genes covered by all the classes of the various categories, are very different. The universe of SPELL or WormExp classes includes essentially all protein coding genes, while the other categories of classes only provide a partial coverage of the entire gene set. For example, most of the genes (57 ± 8 %, n = 13) in the SPELL-derived lists of 200-300 genes that returned the fewest (<5)

WBPhenotype classes (Figure 3B), are not contained within the WBPhenotype universe of 8077 genes. The consistent bias that was observed is a caveat for the interpretation of this type of global enrichment analysis that draws its results from multiple types of gene annotations. In order to circumvent this potential issue, YAAT provides users with the possibility of restricting the enrichment analysis to defined categories of data.

Enriched class clustering to investigate functional connections

When a list of genes is entered, YAAT performs a standard enrichment analysis. To provide a graphical representation of the relationship between the members of the enriched classes, YAAT returns the results as a hierarchical cluster plot. Users can choose between 7 standard clustering methods, with Euclidean-distance based clustering being the default. The source data used for YAAT contains many related data sets. As a single example, from SPELL, there are 4 gene sets, "B.thuringiensis_0.5mix_upregulated_12h", "B.thuringiensis_0.1mix_upregulated_12h", "B.thuringiensis_0.5mix_upregulated_6h" and "B.thuringiensis_0.1mix_upregulated_6h" [START_REF] Yang | Overlapping and unique signatures in the proteomic and transcriptomic responses of the nematode Caenorhabditis elegans toward pathogenic Bacillus thuringiensis[END_REF] for which >70% of all the constituent genes are in more than one class, with 24% common to all 4. Such groups of related genes can strongly influence the clustering analysis and mask relationships that would otherwise be apparent between less-related classes. To counter this, YAAT offers the user the possibility to exclude classes that have too great an overlap from the analysis. Additionally, YAAT offers users the possibility to exclude any number of defined gene lists, which is useful, for example, when querying YAAT with part, or all, of a gene list that is in the database.

Although hierarchical clustering is very powerful, and the utility of these analyses has been demonstrated previously [START_REF] Kim | Coordinated inhibition of C/EBP by Tribbles in multiple tissues is essential for Caenorhabditis elegans development[END_REF][START_REF] Zugasti | A quantitative genome-wide RNAi screen in C. elegans for antifungal innate immunity genes[END_REF], it can be biased because of the hierarchical nature of tree building. This can be avoided by calculating nto-n distances between all classes and displaying the result as a network, where nodes represent enriched classes and edges represent the similarity between the different classes. YAAT therefore includes such a complementary network representation of these results, illustrating the degree of overlap of genes between different significantly enriched classes.

As an example, a set of 144 C. elegans protein-coding genes that are strongly upregulated upon infection with the fungal pathogen D. coniospora (Dc_Up genes; see Methods), returns 39 highly enriched (p<10 -15 ) gene classes. As expected from the global analysis of the dataset, there was a very strong bias of enriched classes; all were derived from transcriptome analyses (Table S2). In part this reflects the lack of representation of these genes in the other categories; only 34/144 were present in the WBPhenotype universe. Relaxing the threshold (P-value < 0.01) revealed the common functional annotation, "pathogen susceptibility increased" for 5 of these genes (fmo-2;

thn-1; thn-2; ilys-3; Y65B4BR.1), since their individual or compound knockdown was associated with loss of resistance to one or more pathogens [START_REF] Evans | Pseudomonas aeruginosa suppresses host immunity by activating the DAF-2 insulin-like signaling pathway in Caenorhabditis elegans[END_REF][START_REF] Luhachack | EGL-9 controls C. elegans host defense specificity through prolyl hydroxylation-dependent and -independent HIF-1 pathways[END_REF][START_REF] O'rourke | Genomic clusters, putative pathogen recognition molecules, and antimicrobial genes are induced by infection of C. elegans with M. nematophilum[END_REF][START_REF] Sahu | Genomic analysis of immune response against Vibrio cholerae hemolysin in Caenorhabditis elegans[END_REF][START_REF] Shapira | A conserved role for a GATA transcription factor in regulating epithelial innate immune responses[END_REF]. They are likely to form part of a general but relatively unexplored pathogen defence strategy that is the subject of current investigation [START_REF] Dasgupta | NHR-49 Transcription Factor Regulates Immunometabolic Response and Survival of Caenorhabditis elegans during Enterococcus faecalis Infection[END_REF][START_REF] Naim | Cell nonautonomous roles of NHR-49 in promoting longevity and innate immunity[END_REF]. This illustrates one way in which YAAT can be used to generate new knowledge at the gene level.

YAAT provides a link to the sets of genes shared between the query list and each enriched class, allowing iterative rounds of analysis, to explore in more depth the relationships between classes. Since the results are available as a downloadable text file, users can represent the results graphically, such as with volcano plots. To provide a graphical representation of the relationship between the members of the enriched classes, YAAT returns the results as a hierarchical cluster plot. With the Dc_Up set of query genes, there were several distinct clusters. Within them, a number of classes clustered as expected. For example, the genes regulated upon D. coniospora infection and in the osmotic stress-related mutants osm-7 and osm-11 overlap significantly [START_REF] Rohlfing | Genetic and physiological activation of osmosensitive gene expression mimics transcriptional signatures of pathogen infection in C. elegans[END_REF] and so co-clustered (Figure 4A; Table S2).

When the similarity between the enriched classes is calculated and then represented in the form of a class linkage network, the connections between a more select group of classes becomes evident (Figure 4B). In this case, the analysis highlights 3 known connections with the anti-fungal innate immunity, again, the overlap with the response to osmotic stress [START_REF] Dodd | A Damage Sensor Associated with the Cuticle Coordinates Three Core Environmental Stress Responses in Caenorhabditis elegans[END_REF][START_REF] Pujol | Anti-fungal innate immunity in C. elegans is enhanced by evolutionary diversification of antimicrobial peptides[END_REF][START_REF] Rohlfing | Genetic and physiological activation of osmosensitive gene expression mimics transcriptional signatures of pathogen infection in C. elegans[END_REF][START_REF] Zugasti | A quantitative genome-wide RNAi screen in C. elegans for antifungal innate immunity genes[END_REF]; the reciprocal relationship with the response to bacterial pathogens that colonize the intestine [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF]; and the link to acs-3/nhr-25 signalling [START_REF] Ward | Defects in the C. elegans acyl-CoA Synthase, acs-3, and Nuclear Hormone Receptor, nhr-25, Cause Sensitivity to Distinct, but Overlapping Stresses[END_REF]. These last 2 categories were not in the same sub-cluster in the dendrogram (Figure 4A, Table S2), illustrating the utility of this complementary analysis.

There was also an overlap with lists corresponding to genes induced upon exposure to 2 concentrations of the bioactive polyphenol quercetin acid [START_REF] Pietsch | Meta-Analysis of Global Transcriptomics Suggests that Conserved Genetic Pathways are Responsible for Quercetin and Tannic Acid Mediated Longevity in C. elegans[END_REF], in part explicable by the >90% overlap between these 2 latter classes. Perhaps more interestingly, this representation revealed an overlap between the Drechmeria-induced genes and genes up-regulated by Burkholderia pseudomallei infection [START_REF] Lee | Burkholderia pseudomallei suppresses Caenorhabditis elegans immunity by specific degradation of a GATA transcription factor[END_REF] but repressed by allantoin [START_REF] Calvert | A network pharmacology approach reveals new candidate caloric restriction mimetics in C. elegans[END_REF]. The commonly regulated genes include a hallmark of epidermal antifungal immunity, 8 antimicrobial peptide genes (Figure 4B). B. pseudomallei degrades ELT-2, an intestinal GATA transcription factor, and so switches off defence gene expression in the gut [START_REF] Lee | Burkholderia pseudomallei suppresses Caenorhabditis elegans immunity by specific degradation of a GATA transcription factor[END_REF]. How this would switch on immune gene expression in the epidermis is not known. Allantoin extends C. elegans lifespan in a DAF-16-independent manner through a mechanism that is not fully understood [START_REF] Calvert | A network pharmacology approach reveals new candidate caloric restriction mimetics in C. elegans[END_REF]. The potential for allantoin to modulate innate immunity merits further investigation. While the results with this test dataset illustrate the confounding effect of including highly-related gene lists in the database, more importantly they suggest that YAAT has the potential to identify functional connections linking different regulatory mechanisms.

Establishing functional connections using co-expression analysis

Genes that function together often are expressed together. Therefore, as well as calculating and returning measures of gene class enrichment, YAAT generates, from a set of ModENCODE RNAseq data [START_REF] Agarwal | Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays[END_REF], information about gene coexpression. Taking the same set of 144 Dc_Up genes, a number of separate clusters can be resolved, including one that includes several antimicrobial peptides genes (Figure 4C). In addition to the annotated cnc, fip and nlp genes [START_REF] Pujol | Anti-fungal innate immunity in C. elegans is enhanced by evolutionary diversification of antimicrobial peptides[END_REF], this clustering draws attention not only to F48C1.9, recently demonstrated to be upregulated specifically in the epidermis upon D. coniospora infection [START_REF] Omi | Inducible expression of F48C1.9 encoding a nematode specific secreted peptide in the adult epidermis upon Drechmeria fungal infection[END_REF] but also to B0563.9 predicted to encode a secreted peptide of around 50 amino acids structural unrelated to previously characterised antimicrobial peptides, and F07C3.9, predicted to encode a 53 amino acid peptide, without a signal peptide, with homologues only in some Caenorhabditis species. Both of these merit further study.

Another clearly discernable cluster contained 5 ttr (TransThyretin-Related family domain) genes (Figure 4C). Several members of this large family of 59 genes have been implicated in C. elegans host defence (e.g. [START_REF] Simonsen | Quantitative proteomics identifies ferritin in the innate immune response of C. elegans[END_REF][START_REF] Treitz | Quantitative proteome analysis of Caenorhabditis elegans upon exposure to nematicidal Bacillus thuringiensis[END_REF]), but their role in the interaction with D. coniospora has not yet been investigated.

These results clearly illustrate how YAAT can help prioritise genes for in-depth study.

Establishing functional connections using phylogenetic analysis

Genes that function together are also often observed to co-evolve [START_REF] Pellegrini | Assigning protein functions by comparative genome analysis: protein phylogenetic profiles[END_REF]. To provide a comparative measure of gene conservation, YAAT displays the average level of conservation across a broad range of species for the entire query list. It also includes in the table of results a comparison between the average level of conservation for the entire list of genes in an enriched class and for the genes that are shared between that list and the genes input as a query. The figures for overall conservation for the members of the different classes varied widely. There was a strong bias towards conservation among the "phenotype" classes but not among those derived from the SPELL expression studies (Table S3). Not surprisingly, among the larger classes (>200 genes), the lowest score was associated with the class of lineage-specific genes [START_REF] Zhou | Genome-wide identification of lineage-specific genes within Caenorhabditis elegans[END_REF]. It was notable that 6 of 30 classes with the least conserved genes were lists of infection-regulated genes, presumably reflecting the lineage-specific adaptation to selection pressure from natural pathogens. At the other end of the scale, with high scores, apart from the classes transposed from Drosophila, which were by their very nature conserved, those related to fundamental cellular processes, including KEGG gene classes, were prominent. The observed skew in the conservation of genes found in the different types of functional classes could lead to another type of bias when analysing datasets.

In addition to giving a global measure of conservation, to provide insight into the pattern of conservation of the members of a gene list, YAAT performs phylogenetic profiling. By default, the results of a broad analysis against 113 diverse species are presented.

Given the increasing interest in translating research from C. elegans to parasitic worms (e.g. [START_REF] Keiser | Is Caenorhabditis elegans the Magic Bullet for Anthelminthic Drug Discovery?[END_REF]), users can opt for a more targeted one, using data from 127 nematode and helminth species. With the list of 144 Dc_Up, most genes were either found only in nematode species or just in C. elegans (Figure 4D). The 5 ttr genes mentioned above were found in the same sub-branch of the phylogenetic tree, increasing the likelihood that they function together in an unexplored aspect of innate defence.

Simultaneous clustering of class enrichment, co-expression and phylogenetic data

Genes that have the same pattern of class membership, co-expression and conservation may very well function together. To facilitate the identification of such groups of genes, YAAT uses a powerful new form of representation. For each set of class enrichment, co-expression and phylogenetic data, the order of genes in the clustering is enforced for the results of the other two sets. To demonstrate the utility of the approach, we queried YAAT with a combined list of genes corresponding to 2 large protein complexes with no common members (Table S4). Clustering by class enrichment resulted in a very good separation of the 2 lists. When this clustering was enforced on the co-expression data, clear groups were discernable (Figure 5A).

Similarly, clustering on the basis of co-expression separated the genes from the 2 complexes in a very satisfactory manner, while projection of this clustering on the class enrichment data also resulted in clear groups of genes (Figure 5B). With this data set, clustering on the basis of evolutionary conservation separated the 2 gene sets relatively poorly, presumably as a consequence of the fact that almost all the genes, regardless of the set, exhibit the same pattern of conservation. There was one marked exception, a group of genes specific to only one complex with a unique pattern of conservation (Figure 5C) that could represent the ancestral core complex.

As another example, we took a published list of genes determined to be induced by the immunostimulatory xenobiotic RPW-24 in an nhr-86-dependent manner [START_REF] Peterson | The nuclear hormone receptor NHR-86 controls anti-pathogen responses in C. elegans[END_REF]. In this case, there was a broader phylogenetic distribution. Looking at the clustering based on evolutionary conservation, one can see a coherent group of genes that were specific to the 5 Caenorhabditis species that also gave a clear cluster of 26 genes across the other 2 graphs (Figure 5D). Despite representing less than 1/5 of the input list, these genes were the majority of those found in 2 classes, "N2_UV_upregulated" and "PA14_vs_OP50_upregulated_8hr" (Table S5). As could be expected, at the gene level, 8 had annotations linked to the response to a range of bacterial pathogens, and 4 to abiotic stress. Interestingly, several correspond to proteins of the same family, in close proximity to each other in the genome (Table S5).

These are all indications for recently evolved genes potentially in cellular defence, and again illustrate how YAAT can help identify groups of genes for in-depth functional study. successive Wormbase releases using an established methodology [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF] proved unworkable after just one update cycle. We are well aware of the material and theoretical difficulties in maintaining any biological database, but have no concrete solutions to the multiple problems we encountered trying to keep YAAT up to date.

DISCUSSION

Regarding the analyses themselves, rigorous statistical enrichment analysis requires knowledge of the gene "universe" sampled when a class was defined. It matters whether a class was defined after sampling all genes or only a subset. For example, in our previous transcriptome studies to determine the response of C. elegans to infection, we have used cDNA [START_REF] Mallo | Inducible antibacterial defense system in C. elegans[END_REF], long oligonucleotide [START_REF] Wong | Genome-wide investigation reveals pathogen-specific and shared signatures in the response of Caenorhabditis elegans to infection[END_REF] and tiling microarrays [START_REF] Engelmann | A comprehensive analysis of gene expression changes provoked by bacterial and fungal infection in C. elegans[END_REF]. These theoretically cover between 35%

and close to 100% of the predicted protein coding genes. In reality, because of the sensitivity of the techniques, these figures are upper limits. Similarly, for analyses by RNA-seq, in principal all genes are assayed, but in fact the number is determined by the depth of sequencing [START_REF] Tarazona | Differential expression in RNA-seq: a matter of depth[END_REF], and in some studies may be restricted to a particular category. For example, a recent study looked specifically for non-coding transcripts regulated by CEP-1 [START_REF] Xu | Analysis of the p53/CEP-1 regulated non-coding transcriptome in C. elegans by an NSR-seq strategy[END_REF]. All of the genes identified in this study that are associated with a Wormbase SPELL class included in YAAT have both proteincoding and non-coding transcripts. Currently <550 such genes are predicted, and this would be the relevant universe for these classes, not the >25,000 non-protein-coding genes. For other types of class, determining the universe is problematic, be it because of the high rate of false negatives in RNAi screens, or lack of relevant information (i.e. degree of saturation) for genetic screens. Indeed, it is very rare for the sampled universe to be reported, and even when it is, this information is generally not available in a machine-readable standardised format. In common with most current functional enrichment analysis tools, in YAAT we assumed that the universe for each set is the entire complement of protein coding genes. This will introduce inevitable inaccuracy in the calculation of the statistical significance of any given enrichment, especially if the universe of a class of interest is small.

YAAT goes beyond the functionality of standard enrichment analysis since it provides class clustering, allowing the inference of functional links between groups of genes.

Further, we have included the option to perform co-expression and phylogenetic profiling, complementary methods for establishing groups of functionally-related genes.

YAAT was originally provided as a web-based service on an academic web server.

Following alteration to the host site's security protocols, external access was lost.

Despite our best efforts over many months, we have not succeeded in reinstalling YAAT elsewhere, in part because of the complexity of the underlying programming environment with the attendant multiplicity of dependencies. Unfortunately, YAAT was conceived before the widespread implementation of containers like Singularity that are specifically designed to circumvent such constraints. We will continue to try to restore and online version of YAAT and implement methods to allow it to be updated. Even if we do not succeed, we hope that our work with YAAT will serve as a model for the type of integrated analysis tool that could be useful for any large-scale exploration of biological function as well as illustrating the various pitfalls that can plague meaningful interpretation of gene lists. One gene subcluster is highlighted, with the positions of the indicated genes marked by the vertical line on the right. The query parameters and full results for all panels are in Table S2. The query parameters, the full species colour code and full results for all panels are in Tables S4 andS5.

  Johnson2 similarity scores (sim = (a/(a+b))+(a/(a+c)); a = number of shared genes, b = number of genes only found in one class, c = number of genes only found in the other class) were calculated between each pair of enriched classes. Nodes represent the enriched functional classes and edges represent the similarity between different classes. The user sets a threshold for the maximum distance between nodes, thereby defining the number of edges. Nodes representing functional classes found in the same hierarchical cluster are given the same colour. This representation complements the heatmap display that suffers from the usual drawback of occasional incorrect clustering seen with agglomerative methods.

  Various stand-alone or web-based functional enrichment analysis tools have been developed over the last 20 years. The results of our tests with YAAT indicate that enrichment analyses should be based on comparisons with the broadest possible collection of functional classes to avoid the biases inherent to each type of gene set. Indeed, by collecting transcriptomic, functional and phenotypic data from the literature, Wormbase, KEGG and other resources, we amassed more than 10,000 gene classes, covering essentially all protein coding genes in Wormbase. Since data evolves, from the outset, we included a pipeline to update data from the respective resources automatically. For the static data such as modENCODE assigned transcription factor targets, the most accurate updating would involve reassigning genes to ChIP-seq peaks on the basis of revised gene structure predictions. Even for RNA seq data, each revision of gene predictions should be accompanied by a re-calling of alignments and read counts. Such practices are today unfeasible. Indeed, even our intermediate strategy, tracking changes in gene identity (and associated annotations) over

Figure 1 .

 1 Figure 1. Sources and types of data in YAAT (WS255). A. Venn diagram showing

Figure 2 .

 2 Figure 2. Distribution of the number of classes returned as a function of the size of the

Figure 3 .

 3 Figure 3. Sets of genes defined by transcriptome, but not phenotype, are preferentially

Figure 4 .

 4 Figure 4. A. Clustering of enriched classes for 144 D. coniospora up-regulated genes.

Figure 5 .

 5 Figure 5. A. Clustering of enriched classes for genes from 2 independent protein

Acknowledgements

Supported by institutional grants from the Institut national de la santé et de la recherche médicale, Centre national de la recherche scientifique and Aix-Marseille University to the CIML, and the Agence Nationale de la Recherche program grants (ANR-16-CE15-0001-01, ANR-11-LABX-0054 (Labex INFORM) and ANR-11-IDEX-0001-02 (A*MIDEX)). We thank Damien Courtine and Lionel Spinelli for their input.

Enriched gene classes

Genes in class