dataesr/scientific-paragraphs-categorization huggingface dataset:
4,5) and explored theoretically  6  , the forecasting of gene expression dynamics over time in relation to environmental conditions remains an open challenge  7. Although previous work  8, 9  have approached the problem in marine systems with relatively stable environmental conditions, such as those associated with marine oxygen minimum zones, using meta-omic (DNA, RNA, protein) and environmental parameter information to model biogeochemical cycles, a generalized framework for time-variable integration of meta-omic datasets into models of community ecology remains to be established. The surface community of biological wastewater treatment plants (BWWTPs) represents an excellent candidate to become a model system to establish such a modelling framework for the following three reasons  10. Firstly, BWWTPs share the challenges linked to most environments, as it is an open system with a constant influx of new populations  11  and exchange of matter and energy with the environment (that is, access to open air and sun irradiation).
However, these challenges can be mitigated by keeping operational parameters (for example, Article https://doi. org/10. 1038/s41559-023-02241-3 being all-purpose powerful (neural network). However, ARIMA assumes that the parameters behind the process are constant, while Prophet can model a time-dependent evolution of the ARIMA parameters. On the other hand, ARIMA is the only model (in the current R package fable  29  implementation) where it is possible to obtain information about the contribution of the individual variables to the forecasting, making it suitable as an explicative model. This allowed us to correctly predict the gene abundance and expression of the populations in the community.
Characterization of the microbial community From the experimental period between 21 March 2011 and 3 May 2012 (ref. 27), we previously obtained and analysed 51 weekly samples, to which we added 21 samples collected in the month of June during the years 2012-2016. The 72 samples were submitted to the same meta-omic analyses-MG, metatranscriptomics (MT) and metaproteomics (MP)-and processed individually to obtain 72 metagenomic assemblies, collections of metagenome-assembled genomes (MAGs), plasmids, viruses, unbinned prokaryotic chromosomal contigs and the corresponding gene expression at the transcriptional and proteomic levels. The combined datasets of the previous time series alongside the new samples were analysed together with updated bioinformatic workflows to allow a coherent comparison between samples along the time series while addressing the batch effect arising from combining the two sets (see Methods and Extended Data Fig.
2a ). To form coherent sets spanning the whole time series, we individually clustered the bins (prokaryotic and eukaryotic) and the contigs (viral, plasmid and unbinned) according to their sequence (see Methods), which led to a total of 144 representative MAGs (rMAGs) and 1,681,736 representative contigs (rContigs), yielding 4,711,952 open reading frames (ORFs) (Supplementary Tables  1  and 2 ). A KEGG Orthology group (KO term) was assigned to 55% of the total retrieved ORFs, while taxonomic affiliations were assigned to 38. 5%. The number of ORF copies as well as their detected gene expression and protein abundances were determined over the extended dataset (see Methods). We found on average 2. 2 × 10 6 ± 4. 8 × 10 5 (s. d. ) ORFs, 9. 1 × 10 5 ± 1. 7 × 10 5 transcripts and 2. 4 × 10 5 ± 2.
5 × 10 4 protein groups per sample. However, most of the genes were not found to be expressed over the entire dataset or were only detected in a few samples. This suggests that an important fraction of the gene pool in the LAOs is not specifically required for community function or their expression levels are below the detection limit, hinting that their cumulative functional effort may be compartmentalized. This finding supports previous results  30  showing how a large portion of the community is redundant and only few functions are keystone. Read recruitment (on the ORF level) per sample was on average 59 ± 9% for the MG and 82 ± 3% for the MT, and peptide matching was 27 ± 4%. The recruitment improved from the previous work on the same datasets, which reported that 26 ± 3% and 27 ± 3% of the MG and MT reads mapped against the MAGs, respectively  27.
The seasonal version of the forecasting method, autoregressive integrated moving average (ARIMA), computes cyclical (seasonality), autoregressive (temporal self-dependence), differencing (difference between consecutive timepoints) and moving-average (averaging of consecutive timepoints) components of a time series  25. It thereby offers a very flexible framework for time-series modelling  25. Here we combine SVD and several time-series algorithms into a generalizable framework for modelling the temporal dynamics of multilayered meta-omics data (Extended Data Fig. 1 ). We demonstrate the power of this framework through analysis of integrated meta-omic and environmental parameter datasets from a microbial community enriched in lipid-accumulating organisms (LAOs) on the surface of the anaerobic tank of the BWWTP in Schifflange (Luxembourg).
The sample set comprises 51 time-resolved samples collected between March 2011 and May 2012 for training, with 21 additional samples collected between 2012 and 2016 for testing and validation. For both sets, the biomolecules were co-extracted  26  and the data for the training set were presented in a previous study  27. We reconstructed the metagenomic (MG) structure of the community, alongside its taxonomy, genetic potential, transcript and protein levels. We employed SVD to extract relevant temporal patterns, which were then clustered into 17 fundamental signals. These were integrated with collected environmental parameters to build an ARIMA model, augmented with seasonal components that could explain the observed signals. Multiple models (ARIMA, Prophet  28  and NNETAR neural network model  29  ) were trained to forecast the signals for the subsequent 5 yr.
Abundance Fig. 1  | Diversity and quality of the rMAGs and their representativeness in the meta-omic dataset. a, Phylogenetic tree of the rMAGs in the LAO community (generated using GTDB-Tk 79 ) contains the 126 bacterial rMAGs in the system (the 18 archaeal MAGs were not included). The heat map ring contains the CheckM quality measures per rMAG (completeness, contamination and MAG-originally strain heterogeneity), which were filtered to be at least 75% complete and at a maximum 25% contaminated (median: 2%). The violin plots contain the timeaveraged (train time series) depth profiles over the contigs forming the rMAG. The two sections of the tree noted as * and ** highlight the strains of M. parvicella and Moraxella sp. , respectively. b. The cumulative length of the contigs (longer than 1,000 nt; see Methods) for the 25 most abundant phyla displayed for the rMAGs and unbinned contigs.
Article https://doi. org/10. 1038/s41559-023-02241-3 Hence, we computed one matrix per omic layer for the six formed taxonomic descriptors (phylum, class, order, family, genus and species) and two functional ones (KO terms and pathways) (Extended Data Fig. 2b ). The resulting 27 matrices (3 original and 24 summary) were used to compute the system's eigengenes (EGs)  22. In previous work, the first EG in a time series was shown to represent 'steady state' gene expression, encapsulating the largest explained variance (EV). Therefore, the first EG (average EV 50 ± 22% in all the datasets) was removed. We screened the subsequent EGs for time dependency (see Methods), selecting a set of 210 EGs, and assessed how much of the data variation they explained (Extended Data Fig. 2c ). To reduce potential redundancy associated with the time-resolved EGs identified across multiple data layers and bring together the same temporal behaviours, we clustered the set of 210 EGs into 17 representative EGs (see Methods).
Forecasting of future timepoints From the analysis of the signals identified in the training datasets, it is already possible to identify five signal groups: (1) alternative basal states, for example, two alternative stable states of abundance/expression (S5, S14); (2) perturbation, that is, standing wave with varying amplitude and frequency (S4, S8, S15); (3) cyclical, that is, standing wave with constant amplitude and frequency (S6, S11, S12); (4) 'crashes', that is, quick shifts in the state and reversion to basal states (S1, S10, S16); and (5) mixed, that is, the other factors (Extended Data Fig. 4c-f ). Alternative stable states, perturbations and crashes (groups 1, 2 and 4) are hard to model without observing multiple times the shift and the perturbation events, respectively. In addition, these scenarios may include permanent shifts into a new community equilibrium or transitory signals in the community that will be eventually resolved (for example, a viral infection).
To forecast such events, experimental information (such as one derived from co-culturing) on microbial interactions would be required, which is beyond the scope of this study. The 17 signals were used to train three models (with various parameters) from the package fable  29  , and the best-performing model on the training set was selected for each of them (see Methods and Extended Data Fig. 9 ). In detail, ARIMA, Prophet and neural network models (with up to four Fourier terms for ARIMA and Prophet) were trained for S1-17 using the environmental variables as external regressors. The 51 weeks spanning the 2011-2012 data were used as a training set as well as to select the three best-scoring models to build a combined Article https://doi. org/10. 1038/s41559-023-02241-3 one (see Methods). In the end, the model with the smallest root mean square error (RMSE) was selected for forecasting.
A total of 21 new samples were collected in the month of June of the subsequent 5 yr to validate the model for the MG and MT data. The month of June was chosen because it is far from disruptive events (such as rain, snow and very cold temperatures) that occur in autumn and winter. To predict the behaviours of the community in these cases, we would have needed a longer training set spanning multiple yearly cycles. To assess the accuracy of the forecasting, we computed the residues of the model and checked whether they were consistent with a white noise distribution. Therefore, we showed in 16 out of 17 cases that the modelling was sufficient to reproduce the training data (Fig. 3 ). There were six cases in which the modelling was fully successful: S1, S2, S4, S5, S10 and S16. The six correctly forecast signals account for 34.
Collection of environmental variables The environmental variables were collected on site by the researcher(s) while they were performing the sampling. These include dry matter, https://doi. org/10. 1038/s41559-023-02241-3 phosphate, nitrate, ammonium, oxygen, conductivity, pH, temperature and oxygen (Supplementary Table  5 ), following previously established protocol  27. The other variables were retrieved from the automated data collection routine of the Schifflange BWWTP, which measures these values online and aggregates them as 2 h averages starting at 1:00. These recordings include the same variables for different parts of the plant (inflow, both vats, outflow) with the addition of other measurements such as the in/outflow volume. For simplicity, we used exclusively the variable pertaining to the inflow, both vats and outflow in this study (Supplementary Table  6 ).
Co-assembly of metagenomics and metatranscriptomics reads All the samples from the training and the test datasets followed the same bioinformatic pipeline. Sample-wise preprocessing of the MG and MT data was performed using IMP (v. 3. 0)  34  (https://git-r3lab. uni. lu/IMP/imp3) with custom parameters, that is, (1) Illumina Truseq2 adapters were trimmed, and (2) the step involving the filtering of reads of human origin step was omitted for the preprocessing. The reads were corrected using BayesHammer  35  per sample, per omic. The resulting MG and MT reads were assembled with metaSPAdes (v. 3. 13. 1)  36  and rnaSPAdes (v. 3. 13. 1)  37  , respectively. The MG and MT reads of each sample were re-assembled together using the contigs and 'highly filtered' transcripts from the first assemblies as trusted contigs.
MG and MT quantification and filtering The filtered MG and MT reads were aligned to the ORF reference set using bwa  58  and sorted using samtools (v. 1. 11)  59. The resulting sorted bam files were processed using bam2hits (v. 1. 0. 9)  60  and the output split with a maximum number of 100,000 ORFs per subset while respecting the bam2hits read groups. Each subset was quantified with mmseq (v. 1. 0. 9)  60  and mmcollapse  61  , then the quantifications per sample were the normalized form of fragments per kilobase million, merged and re-normalized to fragments per kilobase million. Values of gene abundance and expression inferior to 10 -7 were considered equal to 0, and ORFs and transcripts that were not present in at least 20% of the training set were discarded from further analysis.
MP quantification and filtering Raw MP data were retrieved from the PRIDE repository with accession number PXD013655 (ref. 27); samples were processed as previously described  32  and re-analysed. The complete set of predicted ORFs was subsetted to obtain smaller sample-specific databases. The MG alignment files generated in the previous step were processed with featurecounts  62  , and all the ORFs with a count greater than for the given sample were included in the appropriate sample. Each sample-specific database was concatenated with a cRAP database of contaminants (https://thegpm. org/cRAP; downloaded in July 2019) and the human UniProtKB Reference Proteome (UniProt Consortium, 2021), and decoys were generated by adding the reversed sequences of all protein entries to the databases for the estimation of false discovery rates.
The search was performed using SearchGUI (v. 3. 3. 20) with X!Tandem 64 , MS-GF+  65  and Comet 66 as search engines and the following parameters: trypsin was used as the digestion enzyme, and a maximum of two missed cleavages was allowed. The tolerance levels for matching to the database were 10 ppm for MS1 and 15 ppm for MS2. Carbamidomethylation of cysteine residues and oxidation of methionines were set as fixed and variable modifications, respectively. Peptides with length between 7 and 60 amino acids and with a charge state between +2 and +4 were considered for identification. The results from SearchGUI were merged using PeptideShaker-1. 16. 45  (ref. 67), and all identifications were filtered to achieve a protein false discovery rate of 1%. The sample-specific peptide-spectrum matches obtained for each analysis were then used to calculate dataset-wide protein groups using the Occam subgroup method from the Pout2Prot algorithm  68.
Eigengenes clustering and Granger causality network Considering that we required a clustering approach that is independent of scale, we computed the Pearson correlations between pairs of EGs, the output was made absolute and the Minkowski distance was computed. The clusters were retrieved using the 'cutreeDynamic' function (deepSplit=0, pamRespectsDendro=FALSE, minCluster-Size=3) from the dynamicTreeCut package  74  (because it can accommodate a complex structuring of the data), resulting in 17 groups (Extended Data Fig. 2d ). From each of the 17 groups, a representative EG was selected according to the following criteria: (1) MG or MT (because MP data do not exist beyond the training set) and (2) smoothest profile (minimal median of the absolute de-trended time series). The resulting EGs are S1-17 in Fig.
2a. The signals were tested two at a time with the Granger causality test (grangertest) from the lmtest package (v. 0. 9-38)  75  , and if P < 0. 05, the two signals were considered connected. The signals were screened for nonlinearity via empirical dynamic modelling as implemented in the R package rEDM (v. 1. 14. 0)  76. We first identified the best number of lags (embedding value) to analyse the signals using the simplex function and default parameters. The signals were screened with the S-map method  77  , and only three signals appeared to be putatively nonlinear: S7, S8 and S17. All the causal links identified with the Granger causality test were also tested with the convergent cross-mapping method (Extended Data Fig. 5 ) using the function ccm with library size=c(20,50,1) and default parameters.
Modelling the signals and model selection For each signal, we trained multiple models using three techniques (ARIMA, Prophet and neural network using the functions ARIMA, Prophet and NNETAR, respectively, all implemented in the R package fable  29  ), alongside a range of values for the parameters accounting for seasonal components. Each signal was modelled as a separate process whereby the signal itself was the target of the model and the environmental parameters were the only exogenous variables. Therefore, we did not use any information transfer among the signals in the modelling. We fitted ARIMA with up to four Fourier components (see 'Eigengenes and their analysis' section), while the parameters P, d and q were automatically optimized by the function, leading to five ARIMA models (one for each increment of Fourier transform terms, starting with 0).
For the Prophet modelling, we specified seasonality (period = 'year', type = 'additive' and order = from 0 to 4, analogously to the Fourier transform terms of the ARIMA) and growth (type = 'logistic'), resulting in five Prophet models. The neural network function was used whereby the number of nodes in the hidden layer was set to 10, 20 and 30. The 13 models for each signal were scored according to their RMSE, and the 3 models with the lowest RMSE were combined (weighted by 1 - RMSE), as a 14th ensemble model. The RMSE was calculated for the 14 models as well (Extended Data Fig. 9 ). For each signal, the model with the lowest RMSE was selected for the putative generative process and used to forecast the test set with the function 'forecast' of the fable package and supplying environmental parameter readings.
Acknowledgements We thank the  Luxembourg National Research Fund (FNR)  for supporting this work through various funding instruments. Specifically, a  PRIDE doctoral training  unit grant ( PRIDE/15/10907093 ),  CORE grants  ( CORE/17/SM/11689322 ), a  European Union   ERASysAPP grant  ( INTER/SYSAPP/14/05 ) and an  ATTRACT grant  ( A09/03 ) all awarded to P. W. , as well as a  CORE Junior  ( C15/SR/10404839 ) grant to E. E. L. M. The project received financial support from the  Integrated Biobank of Luxembourg  with funds from the  Luxembourg Ministry of Higher Education and Research. This work was also supported by the  European Research Council (ERC)  under the  European Union 's  Horizon 2020 research and innovation programme  (Grant Agreement No. 863664 ). The work of P. M. was funded by the ' Plan Technologies de la Santé du Gouvernement du Grand-Duché de Luxembourg ' through the  Luxembourg Centre for Systems Biomedicine (LCSB) ,  University of Luxembourg.
Code availability The meta-omics pipeline IMP (v. 3. 0)  34  is maintained and developed at the GitLab page: https://git-r3lab. uni. lu/IMP/imp3. The code used in the analysis is available at https://git-r3lab. uni. lu/ESB/lao/lao_ts and https://github. com/fdelogu/microforecast, while the data required to start the analysis are available on Zenodo at https://doi. org/10. 5281/ zenodo. 7225349. The full list of software and R package versions are listed in the Git pages.
Accurate phylogenies are fundamental to our understanding of the pattern and process of evolution. Yet, phylogenies at deep evolutionary timescales, with correspondingly long branches, have been fraught with controversy resulting from conflicting estimates from models with varying complexity and goodness of fit. Analyses of historical as well as current empirical datasets, such as alignments including Microsporidia, Nematoda or Platyhelminthes, have demonstrated that inadequate modeling of across-site compositional heterogeneity, which is the result of biochemical constraints that lead to varying patterns of accepted amino acid along sequences, can lead to erroneous topologies that are strongly supported. Unfortunately, models that adequately account for across-site compositional heterogeneity remain computationally challenging or intractable for an increasing fraction of contemporary datasets.
, 1994; Bergsten, 2005; Brinkmann et al. , 2005; Philippe et al. , 2011b). LBA is a bias in the inferred topology, as a result of which long branches are more likely to branch together than with short ones, independent of the true evolutionary history. Substitution models  (Jukes and Cantor, 1969)  account for the possibility of multiple substitutions per site, and reduce LBA  (Felsenstein, 1973 ; but see  Farris, 1999)  compared to models based on maximum parsimony. Probabilistic substitution models describe the evolution of a site as a series of transitions between states. Every possible transition occurs at a specific 2 rate which is the expected number of transitions from one state to another per unit time. Standard substitution models perform well on most datasets  (Ripplinger and Sullivan, 2010).
Models ignoring across-site compositional heterogeneity are prone to LBA because they underestimate the probability of convergent substitutions at compositionally constrained sites. In particular, the probability of independent substitutions to the same state depends on the number of acceptable amino-acids, which differs across sites. Models ignoring across-site compositional heterogeneity pool all sites, and ignore the variation of the evolutionary process across sites. Indeed, analyses of a series of datasets exhibiting previously contentious evolutionary relationships provide evidence that ignoring across-site compositional heterogeneity can lead to LBA  (Phillips et al. , 2004; Brinkmann et al. , 2005; Philippe et al. , 2005b,a; Delsuc et al. , 2006; Lartillot et al. , 2007; Philippe et al.
, 2008; Schrempf et al. , 2020)  We distinguish between general distribution mixture models estimated from curated training databases, and distribution mixture models directly estimated from the datasets at hand. For example,  Wang et al. (2008)  directly estimate mixture model components using principal component analysis. Susko et al. (2018)  use a composite likelihood approach and additional methods such as taxon weighing. In contrast, the rationale behind providing and using general mixture models is the assumption that the underlying evolutionary processes share universal features. Quang et al. (2008)  use the expectation maximization algorithm to infer general mixture models consisting of 10, 20,. , 60 components (C10, C20,. , C60, collectively CXX models). Schrempf et al. (2020)  used a clustering approach together with different compositional transformations to provide a set of general mixture models, termed universal distribution mixtures (UDM), with the number of components ranging from four up to several thousand.
They also provide the clustering method EDCluster to infer dataset specific distribution mixture models. Finite mixtures can be used in the context of maximum likelihood (ML) phylogenetic inference. However, statistical analyses of model fit and investigation of known cases of LBA indicate that a large number of components are necessary for robustness against LBA  (Schrempf et al. , 2020) , which is computationally intensive and memory intense. Bayesian approaches can more easily accommodate richer mixtures. In partic-4 ular, nonparametric Bayesian methods do not require explicit specification of the number of mixture components nor their stationary distributions. In particular, the CAT model  (Lartillot and Philippe, 2004 ) uses a Dirichlet process prior to approximate an arbitrary mixture of stationary distributions across sites.