dataesr/scientific-paragraphs-categorization huggingface dataset:
To test the performance of our method in such a scenario, we analyzed the moments when an animal slept when a CS+ was being presented in cued fear extinction sessions. Indeed, the unanticipated presentations of a strong sound can wake up the sleeping animal. In the example in Figure  3F ; Supplementary Video 5, the first couple of CS+ presentations induced freezing behavior, which was detected by our method because smoothed spindle power remained low during these immobility periods. After some minutes, the animal fell asleep, which was accompanied by increased smoothed spindle power. Each subsequent CS presentation briefly woke the animal up, resulting in sleep-wake-sleep transitions detected by our method. The smoothed spindle power can therefore serve as a precise tag for SWS. Finally, while ground truth data for the internal state of the animal is impossible to obtain, we compared the results of the algorithm to the scoring by an expert observer, blind to the algorithm results, who manually labeled 2 s video snippets (n = 600) classified as either freezing or sleep by the algorithm across five recording sessions from three different animals.
Data Acquisition and Processing An inertial measurement unit (IMU, non wireless version of the one described in  Pasquet et al. , 2016)  recorded 3D angular velocity and linear acceleration of the animals' heads (sampled at 300 Hz). Animal behavior was also recorded by lateral video cameras (50 Hz sampling rate) in contexts A and B (acA25000, Basler). Brain activity was recorded wideband at 20,000 Hz using a KJE-1001 data acquisition system (Amplipex, Szeged, Hungary). Neurophysiological and behavioral data were explored with NeuroScope  (Hazan et al. , 2006). LFPs were derived from wideband signals by downsampling all channels to 1,250 Hz. At the end of the experiments, recording sites were marked with small electrolytic lesions (∼20 µA for 20 s, one lesion per bundle). After a delay of at least 3 days to permit glial scarring, rats were deeply anesthetized with a lethal dose of pentobarbital and intracardially perfused with saline (0.
Slow Wave Sleep Detection Immobility was detected as above. LFP data was visualized using Neuroscope  (Hazan et al. , 2006 ) and a channel with distinct spindles during sleep (mPFC) was selected. The LFP signal was filtered in the spindle band (9-17 Hz), and the instantaneous amplitude was estimated using the Hilbert transform. This amplitude was then smoothed (Gaussian window of 14 s), and immobility periods where the smoothed amplitude exceeded a k-means identified threshold (Supplementary Figure  2 ) were classified as SWS. The minimum duration allowed for sleep bouts and the maximum duration of small movements are parameters that users are free to select. For the results presented here these parameters were set to 30 and 1 s, respectively.
Inter-species RNA-Seq datasets are increasingly common, and have the potential to answer new questions on gene expression patterns across the evolution. Single species differential expression analysis is a now 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 tools in various settings. Such benchmarks are still missing in the inter-species gene expression context. In this work, we take a first step in this direction by developing and implementing 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.
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 and 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; Gilad et al. , 2006; Khaitovich et al. , 2004; 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). As a compilation of a large number of datasets and atlases for gene expression of healthy wild-type individuals, the well maintained Bgee database  (Bastian et al.
, 2021)  is an important resource to ease the comparison of expression patterns across animal species. Since changes in expression may underlie complex phenotypes, across species gene expression datasets can be used to test a wide range of evolutionary scenarios  (Dunn et al. , 2013; Romero et al. , 2012). Tested hypotheses include for instance expression divergence  (Gu, 2004) ; strength of expression conservation  (Gu et al. , 2019) ; coevolution of gene expression  (Cope et al. , 2020) ; test of the orthologous conjecture  (Rogozin et al. , 2014; Dunn et al. , 2018) ; detection of the "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,b) ; or detection of duplication-specific effects in expression evolution  (Fukushima and Pollock, 2020).
The paper is structured as follows: first, we describe the main methods used to perform differential analysis or shift detection across multiple species. We then explain our simulation method to generate synthetic inter-species RNA-Seq data using a Poisson log-normal model. Our simulation tool is integrated in the Bioconductor package compcodeR. Finally, a targeted simulation study, that draws its parameters from a recent inter-species RNA-Seq study  (Stern and Crandall, 2018a) , allows us to compare the current statistical methods, and propose some recommendations. 2 Review of methods used to compare level of expression across species.
Normalization and Transformations As we will see below, instead of using a generalized linear model on raw count data, it is possible to use a simple linear model on normalized data. The normalization step is essential to transform count measurements into continuous values, and to unlock the use of linear models. The normalization should be designed to temper the sample and gene specific technical biases, as well as to render the data homoscedastic (i. e. with homogeneous variance across samples). Three main normalization scores are used in the literature. They all rely on the normalized library size M i for sample i, defined as: M i = g y gi m i , with m i the scaling normalization factor described above. The Count Per Million (CPM) score incorporates sample-specific normalization only: CPM gi = y gi M i /10 6.
Strategy 2: Linear Model on Normalized Data Assuming the data has been normalized and transformed properly, it can be modelled, for each gene g, using a simple linear regression: Y g = Xθ g + E g , ( 1 ) where Y g is the vector of the n normalized measurements for gene g, E g is a vector of Gaussian independent and identically distributed residuals, and, as previously, X is the design matrix and θ g the associated vector of coefficients. This model is implemented in the popular R package limma  (Smyth, 2004; Smyth et al. , 2005) , that uses an empirical Bayes moderated statistic to test whether the coefficient of θ g associated with the group segregation is significantly different from zero. This method is appropriate to analyze datasets with low sample size, but a large number of genes that are pooled in a hierarchical model to get a better estimation of the variance.
Probabilistic Models and Data Simulation Building on existing RNA-Seq methods  (Robles et al. , 2012; Soneson and Delorenzi, 2013; Soneson, 2014) , we developed a new inter-species 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.
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; Choi et al. , 2017) , which proves essential for inter-specific 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 Normalisation 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 (3) to: µ gi = λ gi gi p h=1 λ hi gi M i. ( 7 ) 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. Lengths Simulation. The lengths are simulated according to the pPLN model described above, with expectations and dispersions empirically estimated from the dataset at hand. (Stern et al. , 2017). The root was dated to 65 million years before the present  (Stern et al. , 2017) , but the tree was re-scaled to unit height for the analyses. The "sight" design (dark purple and light green squares) matches with the biological vision status of the species studied  (Stern and Crandall, 2018a).
Material and Methods Gene Expression Underlying Vision Loss in Cave Animals. We used our new simulation framework to generate realistic synthetic datasets, that were set to mimic the features of a recently published inter-species RNA-Seq dataset  (Stern and Crandall, 2018a) , while varying the level of evolutionary dependence. In this study  (Stern and Crandall, 2018a) , the authors analyzed the molecular mechanisms involved in vision loss in the North American family Cambaridae of crayfish species. They selected 8 blind and 6 sighted crayfish species, for which a time-calibrated maximum likelihood phylogeny is known  (Stern et al. , 2017). 3560 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.
2 ). Base Simulation Parameters. Following the methodology described in the previous section, we simulated a "base scenario" dataset using the estimated crayfish tree re-scaled to unit height (t = 1), with the observed vision status design ("sight" design, see Fig. 2 ), 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.
2 × σ 2 g , with σ 2 g = (σ 2 BM ) g t + s 2 g. We chose a base effect size of 3, with 150 differentially expressed genes out of the 3560 simulated ones. From this base scenario, we varied several parameters in order to study their impacts on the simulated data. Each scenario was replicated 50 times. Star Tree and NB Simulations. To check that our new pPLN framework produced datasets with properties similar to the well known NB framework, we replaced the crayfish tree with a star-tree, that mimics the NB situation where all species and replicates are independent. Tree Group Design. 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.
2 ), 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. Differential Analysis Phylogenetic Asymptotic Effective Sample Size. To quantify the intrinsic difficulty of a design compared to another, we propose a new "differential analysis phylogenetic asymptotic effective sample size" (dapaESS). Given a phylogenetic tree T , we first remove all replicates, so that there are no zero-length branches.
We studied the effect of these choices by testing combinations of other normalization methods (RPKM length normalization; or a simple CPM, i. e. no length normalization), and an other transformation function (square root, as advocated in  Musser and Wagner (2015) ). Scores Used. To assess the performance of the inference methods, based on the list of true (simulated) differentially expressed genes, we computed the number of True Positives (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN). We used the Matthews correlation coefficient  Chicco and Jurman (2020). 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).
Results PLN and NB Simulation Frameworks Produce Similar Datasets. When parametrized to produce the same moments, the pPLN framework on a star tree produces datasets that are similar in difficulty to the classical NB framework (Fig. 3 , first two columns). While limma controls the FDR to the nominal rate, DESeq2 fails to control the FDR in this case with a lot of variance (empirical dispersion range from 0. 1 to 5, see also Fig. 5 ). As it assumes a NB distribution of the counts, DESeq2 suffers from the deviation from this model, as opposed to limma, which performs equally well in both cases. As showed by the countsimQC analysis, the datasets simulated with the pPLN and the NB frameworks have similar features, and are comparable to the original empirical dataset (data not shown, comparison report   (Stern and Crandall, 2018a) , with an effect size of 3, a BM model of evolution with added intra-species variation accounting for 20% of the total variance, on the maximum likelihood tree, with the observed sight groups (see Fig.
2 ). It is compared to a pPLN model with the same parameters, but in a case where all samples are 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 are applied to each scenarios, and their FDR is compared. The black dashed line represents the nominal rate of 5% used to call positives. For limma, the counts are normalized using log 2 (TPM) values. Boxplots are on 50 replicates. available on the GitHub repository github. com/i2bc/InterspeciesDE). While preserving the univariate moments, the tree included in the framework (Fig. 3 , last column), introduces some phylogenetic correlations between the species, and leads to a spectacular loss of power of both methods, with a rate of false discoveries higher than three quarters.
Phylogenetic Data Requires Correlation Modeling. For data simulated according to the base scenario, methods that explicitly model sample correlations (limma cor and phylolm) perform best (Fig. 4 , dark purple line). limma cor exhibits the best behavior with the highest MCC, and a TPR reaching about 80%. Its FDR is still above the nominal rate (median around 10%). Tree group Design Matters. The alt designs produces datasets with the clearest signal, (Fig. 4 , light orange line). In this case, limma cor is able to correctly control for the FDR. Although phylolm methods have slightly higher FDR, they achieve a better TPR reaching almost 100%, leading to a better overall MCC score. At the opposite of the spectrum, the block designs produces datasets with a very weak signal, with differentially expressed genes counts M-A values strongly overlapping a very diffuse non-differentially expressed genes distribution (Fig.
002). This is consistent with the fact that the original dataset uses de novo assembled data, that naturally exclude any small counts, and hence the need for a mean-variance trend correction (see Discussion). (Stern and Crandall, 2018a) , with the observed sight groups, and added intra-species variation accounting for 20% of the total variance). The counts are length-normalized (x axis) using CPM (length not taken into account, none), RPKM or TPM, and transformed using the square root (light orange) or the log 2 (dark purple) functions. Boxplots are on 50 replicates.
Simulation Study Our targeted simulation study illustrates some of the specificities of inter-species 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. 4 ), that make the analysis unreliable and hard to exploit. Indeed, the limma method with added correlation seems to outperform other tools, including phylogenetic comparative methods, 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 analyse inter-species 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 intra-species 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, it produces datasets with comparable features (Fig. 2 ). 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.
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 this review, we focused our attention of the detection of shifts of expression between groups spanning across species. However, inter-species datasets are also used to address many other questions, such as equality of within-species variance, expression divergence, or detection of neutral versus 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.
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 lead to more difficult datasets, 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.
Gu et al. (2019)  recently proposed an empirical Bayes approach to deal with this parameter in an RNA-Seq setting. One possible direction could be to adapt this method to a differential analysis problem. More generally, our simulation studies illustrate the need for new statistical tools for inter-species differential analysis, that would combine the strengths of both the classical RNA-Seq literature, that can deal with the specificities of this noisy data, and the PCM literature, that takes into account the phylogeny, an information that can be crucial to correctly interpret inter-species data. 6 Key Points • Inter-species RNA-Seq datasets have a complex structure, and require a dedicated simulation tool that can generate count data with phylogeny induced correlations and that can take varying gene lengths into account.
• Differential analysis for inter-species RNA-Seq data requires a tool that can take at least within-species sample correlations into account and an adequate length normalisation procedure. • The experimental design of the group allocation on the phylogeny has a strong impact on the differential expression signal, and is well captured by the dapaESS score, that can be computed a priori before any statistical analysis or data collection. Marie-Laure Martin and Guillem Rigaill for useful discussions, and to Claire Ducos, Marie Michel and Sarah Jelassi for their work during their master internship. This work was partly funded by the I2BC and the MI CNRS through the MODELCOG (M. G. ) and X-TrEM projects (Sylvain Merlot). We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020.
Figure5: 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 3410 non-differentially expressed genes is shown as a tile plot, with deeper blues representing high probability values. The M-A values of the 150 differentially expressed genes are shown as red dots.