dataesr/scientific-paragraphs-categorization huggingface dataset:
Transcriptomics Signatures We found that 315, 2,783, and 3,246 transcripts were differentially expressed at D013h, D3, and D6 respectively. The DEG at time D013h, D3, D6 were subsequently mapped onto IPA (Ingenuity Pathway Analysis), leading to the identification of 205, 1,496, and 2,101 transcripts, respectively, that were successfully mapped in the Ingenuity Knowledge Base, with approximately 600 transcripts shared between the 3 time points (Supplementary Text, Supplementary Figure  5 ). There were 13, 109, and 207 pathways significantly modulated (P , 0. 05) at D013h, D3, and D6, respectively. A heat map of the comparison analysis for the top 40 canonical pathways is presented in Figure  2. Canonical pathways are ranked from the highest total score (the negative log of the P-value derived from the Fisher's exact test) across the set of observations.
ACKNOWLEDGMENTS The authors wish to acknowledge  Asma a Fritah-Lafont  for the administrative management of the project all along its maturation and during the execution of the main tasks of the study. Her feedback on the manuscript also allowed improving its understanding. They also wish to acknowledge the technical assistance of  P. -Y. Durand  (sample preparation during the clinical study). They also acknowledge  Laurent Beloeil  for its feedback after reading the manuscript and  Jos ephine Abi-Ghanem  for its contribution to the development of the BioTracs-Mimosa software. All data and metadata will be provided upon request. The BioTracs-Mimosa software and code source are available on GitHub (https://github. com/bioaster/biotracs-m-mimosa).
Modular response analysis (MRA) is a widely used modeling technique to uncover coupling strengths in molecular networks under a steady-state condition by means of perturbation experiments. We propose an extension of this methodology to search genomic data for new associations with a network modeled by MRA and to improve the predictive accuracy of MRA models. These extensions are illustrated by exploring the cross talk between estrogen and retinoic acid receptors, two nuclear receptors implicated in several hormone-driven cancers such as breast. We also present a novel, rigorous and elegant mathematical derivation of MRA equations, which is the foundation of this work and of an R package that is freely available at https://github. com/bioinfo-ircm/aiMeRA/. This mathematical analysis should facilitate MRA understanding by newcomers.
Introduction Modular response analysis (MRA) was introduced to infer the coupling between components of a biological system in a steady-state  [1]. It can be applied to components at different levels of details, e. g. , individual genes or subsystems such as pathways or processes. It relies on the perturbation of individual components, the so-called modules. Various developments of MRA and related methods were recently reviewed  [2]  but, despite its success, MRA mathematical derivation was not provided in a systematic and rigorous manner. We thus reasoned that such an analysis was needed and it would facilitate the understanding of the methodology for newcomers. It is presented as a result and is the basis of the development of an open source R library (aiMeRA) that should make MRA accessible to biology labs without mathematical expertise.
We illustrate the use of the aiMeRA package by investigating the crosstalk between nuclear receptors (NRs) in a breast cancer (BC) cell line. A new extension of the method is also introduced to perform inferences at the genomic scale. The Blüthgen Lab recently released another R package to perform MRA computations  [3] , although with a specific focus on their particular edge-pruning and associated maximum likelihood extension of MRA  [4]  that is not our interest in this study. Estrogen receptors (ERs) belong to the NR superfamily, which act as transcription factors activated upon ligand binding. The two isoforms of ERs (ERα and ERβ) are involved in the control of cell proliferation and exhibit essential functions in tissue development and homeostasis, in particular in organs related to reproduction  [5].
Transcriptional data ERβ and RARβ expression could not be quantified in MELN cells. We hence learned networks involving an ERα module, which transcriptional activity was reported by the ERE/luciferase construct. That is, luciferase mRNA abundance measured ERα activity. ERα mRNA abundance would combine ligand-bound and free amounts of the receptor, but only the ligand-bound ones matter in the model. We did not try to distinguish between RARα and RARγ. We estimated their combined transcriptional activity by the mRNA abundance of the HOXA5 gene and the corresponding MRA module was named RARs. NRIP1 and LCoR activity were determined by their gene mRNA abundance. Since MELN cells are BC cells, we considered the E2-, RA-, or E2+RA-stimulated conditions as basal. That is, perturbations at ERα and RARs were negative (switch to ethanol).
The performance of this new type of MRA model (udMRA) was assessed by its ability to predict module  + 1 activity under the dual siNRIP1/siLCoR condition, as we did above for the other models. To avoid trivially successful predictions on genes that would not vary, we limited the benchmark to the 884 genes above that were also significantly regulated upon siNRIP1 or siLCoR in the E2-stimulated condition. That left us with 60 genes. In Fig. 4E , we report the relative errors observed applying udMRA and comparing with naïve predictions. udMRA yielded significantly better estimates of the added module activity. One could wonder whether perturbation magnitudes during double siRNA interference on the same biological system remain identical. That is, whether filling the vector  in Eq. (  11 ) with 1's at the perturbed module indices (what we did so far) is the best option.
aiMeRA usage The R package was designed to be generally applicable; it relies on the formulae presented here. It is able to work with any quantitative input, including biological and technical replicates. We included functionality to facilitate the definition of MRA model topologies (Fig. 5A ). Model construction only involves the execution of a few generic R functions and network plots can be generated within R directly (Fig. 5B ). It is also possible to export such graphs in the graphML format for loading into Cytoscape  [21]. More details are provided in the package documentation. aiMeRA is available from GitHub; submission to Bioconductor is pending.
Discussion MRA modeling  [1, 2]  is a widely used technique to learn the coupling between the modules of a biological system from perturbation data sets (Fig. 1 ). We introduced a new mathematical derivation of the method and implemented a generic R package called aiMeRA (Fig. 5 ). Application thereof was illustrated on an unpublished data set combining specific qPCR and broad RNA-seq data to explore crosstalk between ER and RAR, two important NRs involved in several tumors such as BC. Data analysis reconstituted some known interactions (Fig. 2 ) and supported a novel hypothesis that reciprocal negative coupling could be mediated by shared corepressors (Fig. 3 ). We showed that MRA transcriptional models trained from RNA-seq data are close to those trained from qPCR (Fig. 4A-C ). Which lead us to introduce an innovative application of MRA to probe genes genomewide, searching for replacement reporters of module activity.
That is, new genes that could be functionally related to an MRA module. In particular, we found that the progesterone receptor gene (PGR) reported on ligand-bound ERα transcriptional activity accurately. That observation indicates the potential value of this new use of MRA models since PGR is a widely used reporter of estrogen activity in BC in the clinic. We additionally investigated the possibility to build hybrid MRA models (udMRA & udMRA. ab) including unidirectional coupling to add modules that were not perturbed in the training data set (Fig. 4D ). In the context of the application presented in this report, i. e. , the transcriptional activity of a system of transcription factors, we showed that the ERα-RARs-NRIP1-LCoR udMRA and udMRA. ab models could outperform naïve predictors (Fig.
4E ). Other biological systems might be amenable to such modified models in the absence of strong back-coupling. The aiMeRA package methods support both RNA-seq data and udMRA models. human genome (hg38) with TopHat v2. 10 [23]  and read counts extracted with HTSeq-Count  [24]. The read count matrix was normalized with edgeR  [25]  TMM algorithm. Data are available from GEO under GSE143956.
Acknowledgments We thank  Simon Cabello-Aguilar  and  Meriem Mekedem  for useful discussions during the development of the project. JC was supported by the  Agence Nationale de la Recherche  (grant  ANR-16-CE17-0002-02 ), the Groupement des entreprises françaises dans la lutte contre le cancer (GEFLUC), and the  Fondation ARC  ( PJA 20151203332 ).
The signals were digitized at 24 kHz and band pass filtered between 0. 5 and 10 kHz using a Tucker-Davis Technologies RX5 Pentusa. Signals were thresholded at 3× standard deviation to isolate spikes, and spikes were sorted after a principal component extraction by KlustaKwik  (Harris et al. , 2000)  and custom-written Matlab (Mathworks) scripts.
Visual Stimuli Stimuli were back projected by a gamma-corrected Plus U2-X1130 85 Hz DLP projector onto a screen (Macada Innovision) placed 18 cm in front of the animal. Full screen size was 60 × 42 cm. Stimuli were produced by scripts using Psychophysics Toolbox 3  (Kleiner et al. , 2007)  running on Matlab. We first mapped the receptive fields of the units at the recording sites by presenting a 5 min movie (5 frames per second) of small white squares (approximately 5 degrees wide) in random positions on black background (ratio of white to black area: 1:30)  (Ahmadlou et al. , 2018). These receptive field positions were used to ascertain that we were recording in binocular V1. The next visual stimuli were full-screen, sine-wave, drifting gratings of 0. 05 cycles per degree. Drift frequencies were  0.
Data Analysis Analysis was done using Matlab scripts. For all stimuli, we computed the evoked visual responses, averaged over the duration of the stimulus, minus the spontaneous rate. The spontaneous rate was defined as the mean rate in the last 0. 5 s before stimulus onset. We averaged the response for each combination of contrast and temporal frequency over all drift directions. Only units were included that had a minimum response (i. e. , the evoked responses minus the spontaneous rate) of 1 spikes per second for at least one combination of temporal frequency and contrast. The response dependence on the temporal frequency was fitted with a difference of Gaussians (d. o. G. ), i. e. , R(f ) = R e exp(-1/2 f 2 /w e 2 ) -R i exp(-1/2 f 2 /w i 2 ), where f is the temporal frequency, R e and w e are the gain and width of the positive Gaussian and R i and w i of the wider negative Gaussian.
The fits were made by minimizing the summed squared error of the fit to the mean responses for all temporal frequencies, using the Matlab fminsearch implementation of the Nelder-Mead simplex algorithm. The fit was rejected if the optimal fit was found for R e < 10 -4 or R i < 10 -4 or w e > 10 × 16 Hz. The d. o. G. fit was used to calculate the optimal temporal frequency and the low and high half-maximum temporal frequencies at which the responses were half of the interpolated maximum response. The response dependence on the Michelson contrast of the stimulus was fitted with a Naka-Rushton function, i. e. , R(c) = R m c n /(σ n + c n ), where c is the contrast and R m , σ, and n are fitting parameters  (Albrecht and Hamilton, 1982). The fits were made by minimizing the summed squared error of the fit to all responses for all contrasts plus very small contributions of σ 2 and (n -2) 2 to reduce the degeneracy in fitting sometimes nearly linear data, using Matlab fminsearch.
From the fit, the C50 value was interpolated as the contrast at which the response would be half of the response at 100% contrast. The explained variance per unit for the temporal frequency and contrast response fits were calculated as 1 -{ i (F(i) -R(i)) 2 }/{ i (R(i) -R m ) 2 } , where F(i) are the fitted values for each frequency or contrast i, and R m is the mean of all R(i)'s. The population models in Section "Divisive Normalization" were fit by minimizing the norm of the difference between all the measured values and the fit values over all parameters, using Matlab fminsearch. For the normalization model, the optimal parameters were σ = 6. 9 and n = 0. 87. For the shunting-extended model, the optimal parameters were σ = 0. 5, τ = 0. 11/Hz, and n = 0. 99. An approximation for the LGN population tuning in the mouse was made by taking the values reported in the literature for the optimal temporal frequency, and low and high halfmaximum responses for the LGN population  (Tang et al.
Genotypic Categories Assignation and Reconstruction of Families We used the NewHybrids 1. 1 software  [47]  to reconstruct the structure of gene flow in our transplantation experiment for each population. In essence, NewHybrids attempts to assign individuals to a mixture of genotypic categories representing the possible structure of the gene flow. Because the C1 samples could be only either from pure origin (local or foreign) or F1 hybrids, we ran a first analysis using only C0 genotypes and C1 genotypes. C0 genotypes contributed to improve allelic frequencies estimation but were not used to estimate the π mixture, which was only assessed using C1 offspring genotypes (see  [47] ). Using these first assignments for the C1 offspring, we then ran a second analysis, integrating this time the C2 offspring samples, so as to determine their own specific π mixture, and to assign them to all possible genetic classes (pure local, pure foreign, F1, F2, backcross with local, backcross with foreign).
0. 6. 5 software  [48]. In particular, COLONY uses multi-locus genotypes to infer sibship among samples. Potential parental genotypes were included in the analyses (84 females and 116 males for Val Travers, 51 males and 67 females for Clarée). We performed long runs, using weak priors and full likelihood method, assuming polygamy for males and females, with inbreeding, for diploid dioecious species. These tests were repeated three times to validate results manually. We then tagged families as either local, foreign, F1, F2 or backcrosses by simply matching the individual assignations obtained from NewHybrids with the family structure obtained from COLONY. In some cases, some families could not safely be assigned because they were composed of more than one type of offspring (for instance, both local individual and hybrid individual were found in the same family, for a same run of the analysis, or between runs).
, 2020)  4 , on the (prior) assumption that they will self-isolate, when they realise they could be contagious. The priors for the non-susceptible and non-contagious proportion of the population are based upon clinical and serological studies reported over the past few weeks; e. g. ,  (Ing et al. , 2020; Stringhini et al. , 2020). Please see the code base for a detailed explanation of the role of these parameters in transition probabilities among states. Although the (scale) parameters are implemented as probabilities or rates, they are estimated as log parameters, denoted by ln  =. 4 Shedding of COVID-19 viral RNA from sputum can outlast the end of symptoms. Seroconversion occurs after 6-12 days but is not necessarily followed by a rapid decline of viral load. However, RNA shedding usually lasts longer than the shedding of infectious virus: many viruses produced are defective in some way, but still present RNA.
Of this effective population, a certain proportion are susceptible to infection: Heterogeneity in susceptibility: This was modelled in terms of a portion of the effective population that are not susceptible to infection. For example, they may have pre-existing immunity, e. g. , via cross-reactivity 5 Available from https://github. com/CSSEGISandData/COVID-19. These timeseries were smoothed with a Gaussian kernel to suppress spurious fluctuations at the weekends. 6 And of those 57%, some will be more exposed than others, conferring a further degree of heterogeneity, e. g. , people working in care homes and hospitals, whose staff show a high seroprevalence: Houlihan, C. , Vora, N. , Byrne, T. , Lewer, D. , Heaney, J. , Moore, D. A. , Matthews, R. , Adam, S. , Enfield, L. , Severn, A. , McBride, A.
Software note The figures in this report can be reproduced using annotated (MATLAB/Octave) code available as part of the free and open source academic software SPM (https://www. fil. ion. ucl. ac. uk/spm/), released under the terms of the GNU General Public License version 2 or later. The routines are called by a demonstration script that can be invoked by typing >> DEM_COVID_I at the MATLAB prompt. At the time of writing, these routines are undergoing software validation in our internal source version control system-that will be released in the next public release of SPM (and via GitHub at https://github. com/spm/). In the interim, please see https://www. fil. ion. ucl. ac. uk/spm/covid-19/. The data used in this technical report are available for academic research purposes from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University, hosted on GitHub at https://github.
Acknowledgements The  Wellcome Centre for Human Neuroimaging  is supported by core funding from  Wellcome  [ 203147/Z/16/Z ]. A. R. is funded by the  Australian Research Council  (Refs:  DE170100128  and  DP200100757 ). A. J. B. is supported by a  Wellcome Trust  grant  WT091681MA. CL is supported by an  MRC   Clinician Scientist award  ( MR/R006504/1 ). The authors declare no conflicts of interest.
In addition, newly available data on mask efficacy enabled the exploration of a third, 'universal mask-use' scenario to investigate the potential population-level benefits of increased mask use in addition to the same threshold-driven reinstatement of SDMs. In this best-case scenario model, 'universal' was defined as 95% of people wearing masks in public, based on the highest observed coverage of mask use globally (in Singapore) during the COVID-19 pandemic to date (Supplementary Information). Two derivative scenarios were also included to assist understanding, nuance and policy resolution around the three boundary scenarios. The first scenario, termed 'plausible reference + 85% mask use' , modeled less than universal mask use in public (85%) in the presence of reinstatement of SDMs. The second was a scenario of universal mask use (95%) in the absence of any NPIs (termed 'mandate easing + universal mask use').
We have therefore tried to benchmark our model against alternative models of the COVID-19 pandemic and fully document our predictive performance with a range of measures  21. In addition, we have provided all the data and model code to enable full reproducibility and increased transparency, provided sensitivity analyses to some of our core assumptions; and presented a range of likely futures  36  in the form of mandate-easing, plausible reference and universal mask-use scenarios (as well as two derivative scenarios thereof) for decision-makers to review. In addition, triangulation of other outputs of the SEIR model, such as the proportion of the population that are affected, are also provided and tested against independent data, in this case seroprevalence surveys (Extended Data Fig. 9 ). Finally, because uncertainty compounds with increased distance into the future predicted, the data, model and its assumptions will be iteratively updated as the pandemic continues to unfold (https://www.
Data identification and processing. IHME forecasts include data from local and national governments, hospital networks and associations, the World Health Organization, third-party aggregators and a range of other sources. Data sources and corrections are described in detail in the Supplementary Information and in the data availability statement. Briefly, daily confirmed case and death numbers due to COVID-19 are collated from the Johns Hopkins University data repository; we supplement and correct this dataset as needed to improve the accuracy of our projections and adjust for reporting-day biases (Supplementary Information). Testing data are obtained from Our World in Data (https://ourworldindata. org/), The COVID Tracking Project (https://covidtracking. com/) and supplemented with data from additional government websites (Supplementary Information).
Social distancing data are obtained from a number of different official and open sources, which vary by state (Supplementary Information). Mobility data are obtained from Facebook Data for Good (https://dataforgood. fb. com/docs/covid19/), Google (https://www. google. com/covid19/mobility/), SafeGraph (https://www. safegraph. com/dashboard/covid19-shelter-in-place/) and Descartes Labs (https://www. descarteslabs. com/mobility/; Supplementary Information). Mask-use data are obtained from the Facebook Global Symptom Survey (in collaboration with the University of Maryland Social Data Science Center), the Kaiser Family Foundation, YouGov COVID-19 Behavioural Tracker survey (https://today. yougov. com/ covid-19/) and PREMISE (https://www. premise. com/covid-19/; Supplementary Information). Specific sources for data on licensed bed and ICU capacity and average annual utilization in the United States are detailed in the Supplementary Information.
Before modeling, observed cumulative deaths are smoothed using a spline-based smoothing algorithm with randomly placed knots 37. Uncertainty is derived from bootstrapping and resampling of the observed deaths. The time series of case data is used as a leading indicator of death based on an infection fatality ratio (IFR) and a lag from infection to death. These smoothed estimates of observed deaths by location are then used to create estimated infections based on an age distribution of infections and on age-specific IFRs. The age-specific infections were collapsed into total infections by day and state and used as data inputs in the SEIR model. Detailed descriptions of data smoothing and transformation steps are provided in the Supplementary Information.
These covariates are representative at a population level and are time invariant. 9] [40]  ). Time-varying covariates include pneumonia excess mortality seasonality, diagnostic tests administered per capita, population-level mobility and personal mask use. These are described below. Pneumonia seasonality. We used weekly pneumonia mortality data from the National Center for Health Statistics Mortality Surveillance System (https://gis. cdc. gov/grasp/fluview/mortality. html) from 2013 to 2019 by US state. Pneumonia deaths included all deaths classified by the full range of the International Classification of Disease codes in J12-J18. 9. We pooled data over available years for each state and found the weekly deviation from the annual, state-specific mean mortality due to pneumonia. We then fit a seasonal pattern using a Bayesian meta-regression model with a flexible spline and assumed annual periodicity (Supplementary Information).