dataesr/scientific-paragraphs-categorization huggingface dataset:
Funding: This work has benefited from an " Investissement d'Avenir " grant (CEBA, reference:  ANR-10-LABX-0025 ) managed by  Agence Nationale de la Recherche. This project has received financial support from the  CNRS  through the  80|Prime program. The  UltrafleXtreme MALDI ToF/ToF  mass spectrometer used in this study was funded by a grant from the  Région Ile-de-France  and the  CNRS-ICSN (DIM Analytics Equipements Mi-Lourds 2012 ). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.
Microorganisms 2022 , 14 Figure 6. 2022146 Figure 6. Visualization of 230 lipid fingerprints of environmental fungi (230 isolates-BNSB's database) set by t-SNE. Each node represents the lipid fingerprint of an isolate and is colored according to the taxonomic order of the organism. Arrows indicate illogical attributions within clusters (annotated Q2, X, Y and Z).
Interspecies RNA-Seq datasets are increasingly common, and have the potential to answer new questions about the evolution of gene expression. Single-species differential expression analysis is now a well-studied problem that benefits from sound statistical methods. Extensive reviews on biological or synthetic datasets have provided the community with a clear picture on the relative performances of the available methods in various settings. However, synthetic dataset simulation tools are still missing in the interspecies gene expression context. In this work, we develop and implement a new simulation framework. This tool builds on both the RNA-Seq and the phylogenetic comparative methods literatures to generate realistic count datasets, while taking into account the phylogenetic relationships between the samples.
We illustrate the usefulness of this new framework through a targeted simulation study, that reproduces the features of a recently published dataset, containing gene expression data in adult eye tissue across blind and sighted freshwater crayfish species. Using our simulated datasets, we perform a fair comparison of several approaches used for differential expression analysis. This benchmark reveals some of the strengths and weaknesses of both the classical and phylogenetic approaches for interspecies differential expression analysis, and allows for a reanalysis of the crayfish dataset. The tool has been integrated in the R package compcodeR, freely available on Bioconductor.
Introduction The study and analysis of gene expression differences across species is a long-standing problem  (King and Wilson 1975). The development of microarray technologies led to the gathering of the first large scale across species gene expression datasets, that allowed for the formulation and study of various hypotheses regarding the link between gene expression and evolution  (Enard 2002; Khaitovich et al. 2004; Gilad et al. 2006; Whitehead and Crawford 2006). RNA-Sequencing technologies have changed the way to measure gene expression  (Wang et al. 2009) , making comparisons across several species easier, even for species with no reference genome available  (Perry et al. 2012; Romero et al. 2012). Interspecies gene expression data are increasingly common, with well-curated resources, such as the Bgee database  (Bastian et al.
2021) , that make it available to the community. Since changes in expression may underlie complex phenotypes, across species gene expression datasets can be used to test a wide range of evolutionary scenarios  (Romero et al. 2012; Dunn et al. 2013). Tested hypotheses include, for instance, expression divergence  (Gu 2004) ; the strength of expression conservation  (Gu et al. 2019) ; the coevolution of gene expression  (Cope et al. 2020) ; test of the orthology conjecture  (Rogozin et al. 2014; Dunn et al. 2018) ; the detection of "phylogenetic signal"  (Musser and Wagner 2015) ; equality of within-species variance  (Catalán et al. 2019) ; constant stabilizing selection, loss through drift, parallel, or divergent selection  (Stern and Crandall 2018a, 2018b) ; or the detection of duplication-specific effects in expression evolution  (Fukushima and Pollock 2020).
We consider the level of expression of genes as a quantitative trait evolving across several species, and we detect genes with a shift in the level of expression across species as performed in for example,  Brawand et al. (2011) ,  Perry et al. (2012) ,  Torres-Oliva et al. (2016) , and  Stern and Crandall (2018a). The specificities of interspecies RNA-Seq data are multifold. RNA-Seq data are counts, usually measured on a low number of samples. In addition, several technical biases affect the measured level of expression of a given gene in a given sample, either gene-specific (such as heterogeneity of gene length and GC content across genes and samples), or sample-specific (such as heterogeneity in library size across samples). Finally, since the level of expression of a gene is measured across several species, the phylogenetic relationships between species induce some correlations in the data.
Simulation Framework We designed and implemented a new simulation framework, that we used to benchmark and calibrate differential analysis methods on synthetic datasets that exhibit features specific to interspecies RNA-Seq dataset. Our approach is based on a phylogenetic Poisson log-normal (pPLN) model, that relies on two layers. First, latent variables are simulated using a Brownian motion (BM) or an Ornstein-Uhlenbeck (OU) stochastic process on the phylogenetic tree, using well-studied tools from the PCM literature. This layer can account for correlations induced by the phylogenetic relationships, including intraspecies-independent variations. Second, the log values of these latent variables are used to define the parameter of a Poisson distribution, from which the final simulated count values are drawn.
Built on top of tools tailored for RNA-Seq, the framework can simulate realistic count data that reproduce some of the important characteristics of a given empirical dataset, with genes FIG. 2. The base scenario [pPLN (real tree), right] had empirical moments drawn from  (Stern and Crandall 2018a) , with an effect size of 3, a BM model of evolution with added intraspecies variation accounting for 20% of the total variance, on the maximum likelihood tree, with the observed "sight" groups (see fig. 1 ). It is compared with a pPLN model with the same parameters, but in a case where all samples were independent [pPLN (star tree), middle], and to a NB model with the same moments and effect size (NB, left). The DESeq2 (black) and limma (light orange) inference methods were applied to each scenario.
Reanalysis of the Crayfish Dataset We applied our simulation framework with parameters empirically drawn from a recent study on the evolution of gene expression underlying vision loss in cave animals  (Stern and Crandall 2018a). The synthetic datasets had characteristics similar to the original data, while allowing us to control the actual simulation parameters, and hence assess the performance of various differential expression analysis tools. We tested some of the most popular methods, along with several normalization and transformation strategies, both issued from the RNA-Seq or the PCM literature. (See supplementary section A, Supplementary Material online for a survey of all the methods used. ) Enlightened by this focused benchmark, we proposed a reanalysis of the original dataset, resulting in a list of possibly differentially expressed genes that is much shorter than previously published one, and that we believe is more robust.
Results Simulation Studies  Stern and Crandall (2018a)  collected an interspecies RNA-Seq dataset to study the molecular mechanisms involved in vision loss in the North American family Cambaridae of crayfish species. We exploited our new simulation framework to generate realistic synthetic datasets following a base scenario, that was set to mimic the features of the crayfish dataset, using the estimated crayfish tree with the observed vision status design ("sight" design, see fig. 1 ), and matching the empirical counts and gene lengths moments. From this base scenario, we varied several parameters in order to study the impact of evolutionary dependence on the simulated data. In all simulated datasets, we could control exactly which genes were generated as differentially expressed, and which genes had a constant expression across clades (see section "A framework to simulate interspecies RNA-seq datasets" in the Materials and Methods).
This allowed us to compare the list of truly differentially expressed genes with the list of candidate genes found by the various statistical methods tested. We tested the performance of the following differential expression analysis methods: DESeq2  (Love et al. 2014) ; limma  (Ritchie et al. 2015) , or limma cor  (Smyth et al. 2005) ; and phylolm  (Ho and Ané 2014a)  with BM or OU processes. See Materials and Methods Section FIG. 4. M-A plots (log2 fold change as a function of the mean of normalized counts for each gene) of the datasets produced with base pPLN parameters (effect size of 3, BM model with added intraspecies variation accounting for 20% of the total variance), on the maximum likelihood tree  (Stern and Crandall 2018a) , with the "block" (left) and "alt" (right) designs. The M-A values distribution for the 3,410 nondifferentially expressed genes is shown as a tile plot, with deeper blues representing high probability values.
PLN and NB Simulation Frameworks Produce Similar Datasets To check that our new pPLN framework produced datasets with properties similar to the well known NB framework as implemented in compcodeR  (Soneson 2014) , we replaced the crayfish tree with a star-tree, that mimics the NB situation where all species and replicates are independent. When parametrized to produce the same moments, the pPLN framework on a star tree produced datasets that were similar in difficulty to the classical NB framework (fig. 2 , first two columns). While limma controlled the FDR to the nominal rate, DESeq2 failed to control the FDR. However, both methods had a better TPR under the pPLN model (see supplementary fig. S1 , Supplementary Material online), and DESeq2 had the best MCC score under the pPLN on a star tree.
Tree Group Design Matters The group design on the tree is known to strongly impact the properties of the data, in particular through its "phylogenetic effective sample size"  (Ané 2008; Bartoszek 2016 ). To study its effect in a gene expression context, we replaced the "sight" design with a "block" and "alt" design (see fig. 1 ), that were chosen to model two extreme situations. In the "block" design, all the species with a given group are nested within a single clade, so that the differential expression signal is redundant with the phylogenetic signal. At the other end of the spectrum, the "alt" design was chosen so that sister species are in different groups, in order to maximize the contrast between organisms that share a long common history. We expect the "alt" design to produce datasets with a stronger signal.
The "alt" designs produced datasets with the clearest signal (fig. 3 , light orange line). In this case, limma cor was able to correctly control for the FDR. Although phylolm methods had slightly higher FDR, they achieved a better TPR, reaching almost 100%, leading to a better overall MCC score. At the opposite of the spectrum, the "block" design produced datasets with a very weak signal, with differentially expressed genes counts M-A values strongly overlapping the nondifferentially expressed genes distribution, which was very diffuse (fig. 4 ). All methods applied to the "block" design had FDR higher or equal to about 50% (fig. 3 , black line). The BM phylolm tool had the least bad MCC score (about 0. 5), although with the worst TPR (around 50%). The relative difficulties of each design was correctly captured by the normalized dapaESS (dapaESSn = dapaESS/dapaESS(Ind), so that dapaESSn = 1 in the independent case).
Reanalysis of the Crayfish Dataset While  Stern and Crandall (2018a)  found a list of 93 differentially expressed genes (see supplementary table S2, Supplementary Material online in Stern and Crandall 2018a), the limma cor method found evidence for only 6, with only one gene that was not in the previous list. Among those, one gene was clearly associated with vision (coding for the Rhodopsin protein). The phylolm OU method found evidence for 17 differentially expressed genes, including the same 5 genes common to OUwie and limma cor. Raising the threshold from 0. 05 to 0. 1 A Phylogenetic Framework to Simulate Synthetic Interspecies RNA-Seq Data • https://doi. org/10. 1093/molbev/msac269MBE (respectively, 0. 2), the limma cor method selected another 4 (resp. , 25) genes, including 1 (resp. , 6), matching with the list from  Stern and Crandall (2018a).
Simulation Study Our targeted simulation study illustrates some of the specificities of interspecies RNA-Seq differential expression analysis. First, it is essential to take the correlation between replicates within a given species into account. Failure to do so leads to very high rates of false discoveries (fig. 3 ), that make the analysis unreliable and hard to exploit. Indeed, the limma method with added correlation seems to outperform other tools, including PCMs, in many settings. These results tend to indicate that, even if the full tree is not included in the analysis, incorporating these simple correlations between replicates might be sufficient to efficiently analyze interspecies datasets, at least for some simulation designs. The group design on the tree was indeed found to be extremely important (fig.
Simulation Design In this work, we proposed a method to simulate RNA-Seq gene expression across multiple species. Similar to intraspecies simulation tools  (Dillies et al. 2013; Soneson and Delorenzi 2013; Soneson 2014) , our simulation method can use empirical datasets to set the value of parameters such that the simulated datasets are as close as possible to the real ones, with matching empirical marginal expectation and variance. When applied to independent species, our pPLN model produces datasets with features comparable to the classical NB model (fig. 2  and supplementary section D, Supplementary Material online). In our specific simulation studies, we use the dataset from  Stern and Crandall (2018a). This dataset was obtained using de novo assembled data. In addition, we focused on genes with one-to-one orthologuous relationships across species.
As a consequence, this dataset had a low number of zeros and small counts, and a large variance across samples (empirical dispersion ranged from 0. 1 to 5, see fig. 4 ). The simulated datasets had similar characteristics, which could explain the low performance of DESeq2 in terms of FDR, even when the data were simulated without correlation (fig. 2 ), and the fact that the trend procedure did not add any power to the limma method. In addition, as DESeq2 explicitly assumes a NB distribution of the counts, it suffers from the deviation from this model, as opposed to limma, which controlled the FDR to the nominal rate in both the NB and pPLN models (fig. 2 ). Interspecies RNA-Seq gene expression datasets are very diverse, with specificities depending on the underlying biological question being studied.
Although it does not take into account the number of replicates or the specific evolution model, we showed that it could well represent the difficulty of an experimental design. The strength of this criterion is that it only depends on the timed species tree and the tips group allocation, and can be computed before any statistical inference or even data collection. It can hence be used as a practical guide on the expected power of the experimental design. In particular, if the normalized dapaESSn is lower than one, results from methods that do not take the phylogeny into account should be interpreted with particular caution. In this work, we focused our attention on the detection of shifts of expression between groups spanning across species. However, interspecies datasets are also used to address many other questions, such as equality of withinspecies variance, expression divergence, or detection of neutral vs.
directed evolution regimes. Several tools from the PCM literature have been used to this end, that rely on various models of trait evolution with appropriate parameter constraints. Since our simulation tool is modular, Bastide et al. • https://doi. org/10. 1093/molbev/msac269 MBE those various processes could be implemented, in order to produce realistic RNA-Seq datasets with the desired structure. More generally, our tool could be extended to take into account any correlation structure between samples, not necessarily deriving from a phylogenetic model. Such an extended framework could help researchers to test the statistical properties of these complex inference models.
Inference Tools In this study, we focused on a few inference tools, that come either from the RNA-Seq or the PCM literature, limiting ourselves to methods implemented in R and that can do differential analysis. Although a more comprehensive simulation design would be needed to draw stronger conclusions, our results show that simulations under the OU model could lead to more difficult datasets for some group designs, and that even methods that include the OU model in their framework fail to completely correct for this effect. This could be linked with the fact that the estimation of the selection strength in an OU model is a notoriously difficult question, especially on an ultrametric tree  (Ho and Ané 2014b; Cooper et al. 2016). Having to estimate this parameter for thousands of genes is bound to generate some instability, and to deteriorate the performance of those tools.
Reanalysis of the Crayfish Dataset In the setting that was most similar to the empirical Crayfish dataset, we found that the limma cor method worked best, with a similar TPR but smaller FDR compared with phylogenetic methods. When we applied this method to the Crayfish dataset, we found a list of only 6 differentially expressed genes between sighted vs. blind species across all clades. This list was robust to the choice of the detection threshold, as raising it from 0. 05 to 0. 1 only added 4 candidates. Allowing more false discoveries with a threshold of 0. 2 output a total of 31 genes including only approximately a third matching the original list of 93 genes found by  Stern and Crandall (2018a) , comforting our suspicion that the data only support a limited number of differentially expressed genes.
Other tools, designed to test for other patterns of evolution as mentioned above, might be able to detect other genes. Finally, when using RPKM instead of TPM normalization, we found a list of genes that only partially matched with the previous ones. Although TPM and RPKM seemed to performed similarly well on the simulated data (see fig. 7 ), this result shows that normalization can have a strong impact when analyzing biological datasets, and the robustness of these methods should be carefully assessed, when possible. Unfortunately, still little is known regarding the evolutionary genetics of vision loss in crayfishes, which makes the biological validation of the results difficult.
A Framework to Simulate Interspecies RNA-Seq Datasets Building on existing RNA-Seq methods  (Robles et al. 2012; Soneson and Delorenzi 2013; Soneson 2014) , we developed a new interspecies simulation framework that can generate realistic count datasets, and takes into account, first, the gene expression correlations induced by the phylogeny and, second, the different lengths a given gene can have in different species.
For each gene g, the dispersion parameter α g is the same for all samples, while the expression level λ gi can only take two values: λ gS 1 if i is in S 1 and λ gS 2 if i is in S 2. Given λ gS 1 , we take λ gS 2 as λ gS 2 = λ gS 1 if g is not differentially expressed; λ gS 1 × (e + X e g ) if g is up-regulated in S 2 ; λ gS 1 × (e + X e g ) -1 if g is down-regulated in S 2 ; ⎧ ⎪ ⎨ ⎪ ⎩ with e the minimal differential effect size, and X e g random variables independent identically distributed according to an exponential distribution with parameter 1. The values of the parameters are set to match the empirical counts expectation and dispersion of a real datasets.
Realistic Simulations using the Poisson Log-Normal Distribution The Poisson log-normal (PLN) distribution has been advocated as an alternative to the NB distribution for the analysis of RNA-Seq data. Being more flexible, it is particularly well suited in the presence of correlations  (Gallopin et al. 2013; Zhang et al. 2015) , which proves essential for interspecies datasets, as demonstrated in the next section. We show here how the parameters of a PLN model can be chosen to match first-and second-order moments of the NB model described above, making it possible to simulate realistic datasets under this more flexible framework. The PLN Distribution. Under the PLN model, for each gene g and sample i, we assume that the observed count random variable Y gi follows a Poisson distribution, with log parameter a Gaussian latent variable Z gi , such that Z gi ∼ N m gi , σ 2 g 􏼐 􏼑 Y gi | Z gi ∼ P exp (Z gi ) 􏼁.
Taking Differential Gene Lengths into Account Length Normalization of Counts. Let ℓ gi denote the length of the gene g for sample i. Following  Robinson and Oshlack (2010) , we take this length into account by changing equation (1) to μ gi = λ gi ℓ gi 􏽐 p h=1 λ hi ℓ hi M i. ( 5 ) Note that the same overall sequencing depth M i is attributed to each sample, but that, because of the weighted average, it is preferentially allocated to longer genes. A Phylogenetic Framework to Simulate Synthetic Interspecies RNA-Seq Data • https://doi. org/10. 1093/molbev/msac269MBE Lengths Simulation. The lengths are simulated according to the pPLN model described above, with expectations and dispersions empirically estimated from the dataset at hand.
Gene Expression Underlying Vision Loss in Cave Animals The real dataset used to set simulation parameters and for the real data analysis case study was extracted from  (Stern and Crandall 2018a). In this study, the authors selected eight blind and six sighted crayfish species, for which a time-calibrated maximum likelihood phylogeny is known  (Stern et al. 2017). 3,560 orthologous gene expressions were estimated using the method RNA-Seq by Expectation Maximization (RSEM)  (Li and Dewey 2011) , with one to three replicates per species (see fig. 1 ).
Base Simulation Parameters We used the real dataset to set the parameters of our simulations. We took the estimated crayfish tree rescaled to unit height (t = 1), with the observed vision status design ("sight" design, see fig. 1 ), and matching the empirical counts and gene lengths expectation and dispersion. The expression level λ gS 1 and the dispersion α g were estimated from the dataset for each gene g, while for each sample i the simulation sequencing depth M i was independently drawn from a uniform distribution with bounds M min and M max the observed empirical minimal and maximal values of the library size across all samples. We used a BM model of trait evolution, with an independent layer of individual variation s 2 g representing 20% of the total tip variance σ 2 g for each gene g: s 2 g = 0.
MBE We also computed the true positive rate (TPR = TP/ (TP + FN)) and the false discovery rate (FDR = FP/ (FP + TP)). In addition, we compared the features of the simulated datasets with the empirical one using the countsimQC R package  (Soneson and Robinson 2018). Reanalysis of the Crayfish Dataset  Stern and Crandall (2018a)  used the OUwie package  (Beaulieu et al. 2012)  on each gene to compare an OU model with a single optimal value to a model with two optimal values, one for the sighted species and one for the blind. This method is similar to the phylogenetic ANOVA, but differs on two important aspects. First, it takes as entry the within-species empirical means and variance instead of the individual values. Second, it uses a likelihood ratio test assuming a chi-square distribution with one degree of freedom, instead of the conditional t-test used in phylolm.