dataesr/scientific-paragraphs-categorization huggingface dataset:
Introduction Massively parallel reporter assays (MPRAs) have transformed our ability to directly investigate the effects of genetic perturbations on cellular phenotypes such as gene expression, growth, and chemical resistance. They allow researchers to test tens of thousands of sequences at once, which provides vast amounts of data to study complex phenomena while minimizing cost and the impact of batch effects. MPRAs are now a fundamental tool and have been used to identify enhancer activity  [1] , determine the influence of chromatin on cis-regulatory sequences  [2] , identify sources of variation in red blood cell traits  [3] , and identify functional annotations in 3 untranslated regions  [4]. Related assays such as CRISPR screens have identified loss-of-function mutations affecting tumor growth and metastasis in mice  [5]  and genes that regulate T cell activation in the genome  [6, 7].
1 We collectively term these approaches as bottom-up: the analysis objective is to identify which individual perturbations have a causal impact on the phenotype of interest. In this work we develop a top-down MPRA analysis framework: the analysis objective is to determine the contribution of all sequence variation to experimental outcomes across the entire experiment. MPRAudit takes as input annotated subsets of the MPRA experiment and outputs estimates of the biological and technical variance across these subsets. These parameters have a range of utilities across MRPA use cases, and can be applied even when there is not sufficient power to identify individual causal perturbations with bottom-up approaches. For studies that consider perturbations to individual elements, MPRAudit provides an estimate of the quality of the dataset and sets an upper limit to the performance of prediction assays.
For studies that consider pairs of elements with and without a mutation, MPRAudit can additionally determine the contribution of the flanking sequence of a mutation to the effect size of that mutation. For studies that consider groups of elements, MPRAudit tells researchers how much novel biology is left once known biological classes are accounted for; and what subsets contain unknown biological sources of variation waiting to be discovered. We evaluate MPRAudit through application to simulated and real data sets. In realistic simulations, we show that MPRAudit gives accurate estimates of the variance explained by technical factors, and that MPRAudit's uncertainty estimates accurately match the empirical variation across simulations. We then apply MPRAudit to a simulated MPRA experiment, where we demonstrate the use of MPRAudit on simulated data of groups of gene expression measurements, paired effect size measurements, or groups of measurements.
Finally, we apply MPRAudit to experimental data, including our recent MPRA on the 3 UTR in human cell lines, where we demonstrate the use of MPRAudit on categorical data and show that the cis-regulatory effects of the 3 UTR sequences on gene expression remain largely unexplained by available models. MPRAudit is implemented as a free and open source software package. MPRAudit can be used in any barcoded molecular assay, as it is based on holding out or resampling over barcoded clones within a dataset and it does not require parametric model assumptions.
If A i varies significantly with sequence (i. e. i), then the sequence-based variation will be nonzero and the total variation will differ from the technical variation. In real experiments there might be additional forms of sequence-based variation from the PCR process (due to GC content, for instance), which we do not include in this simulation, but we do account for in the real analysis. Given these numbers of simulated RNA and DNA molecules, simulated counts of RNA and DNA for each clone and sequence are drawn from a negative binomial distribution to simulate overdispersion. Further details are given in the methods section. The heavy tails, truncated distributions, and complex relationships between clonal counts and aggregated outcomes in these types of datasets would pose problems for parametric methods, which is why we favor the flexibility and accuracy of a resampling-based approach.
Simulation Results for Single Perturbations One of the primary use cases that we envision for MPRAudit is to quantify the total fraction of sequence variation across the entire experiment. The explainability b 2 sets an upper limit on the accuracy of prediction methods, and provides a measurement of the reproducibility of the dataset and experiment as a whole. In addition, when a dataset is divided into groups (for instance by the presence or absence of a binding site motif), a comparison of b 2 across groups determines the remaining biological signal left in each group, which could be informative for followup investigations. We first simulate the null by setting A i = 2 to be constant for all sequences, guaranteeing that b 2 = 0 and that sequence will have no effect. The results shown in Figure  2  demonstrate that the bias is uniformly small, and rapidly decreases as the number of clones per sequence increases.
For the dataset simulated here, we expect the b 2 within the 5,000 "TATACAG" sequences to be lower than the b 2 within the 5,000 completely random sequences, since the TATACAG region has no variation. This is what we observe in Figure  3. In Figure  3e , the data is simulated under the null and there is no b 2 in either group. In Figure  3f , the data is simulated under the alternate with A i = 2 -0. 5f 1i + 0. 15f 2i , and b 2 is shown to be higher for the group with completely random nucleotides. In Figure  3g  the bias within each subgroup is shown to be similarly small.
Overview To demonstrate the utility and performance of MPRAudit in realistic settings, we looked at 12 real MPRA data sets from two different MPRA technologies. In each case we require at least 5 clones or barcodes per sequence. We also require a minimum number of counts per clone to filter out sequences with low statistical power and we follow the analysis guidelines of the published work when given. The first several datasets are from our fast-UTR study of 3 UTR transcripts in human cell lines. We analyze data from two cell lines (Jurkat and Beas2B)  [4]. In the design of this study there are 41,288 distinct 3 UTR sequences (including common variants and deliberate mutations) from 4,653 genes, which are each 160 nt long. Each sequence has twenty clones on average, which differ only by an 8 nt barcode.
Not all sequences were observed in all sequencing runs. For this experiment we require each measured sequence to have at least five clones with more than five counts of DNA, and for at least one clone of each measured sequence to have at least one count of RNA. Any sequence that does not satisfy these conditions is excluded from the study. We then analyze sequences according to the statistic RNA/(RNA+DNA): Ôi = ni c=1 Rc ni c=1 Rc + ni c=1 Dc. ( 8 ) The study is designed to investigate the effects of several classes of active motifs and introduces mutations to determine the impact of those motifs  [10, 17, 18, 35, 24]  This study considers single sequences, pairs of sequences, and groups of sequences in the course of its investigation. The second dataset is from a publicly available MPRA  [23]  where the authors compare the functional activities of 2,236 candidate liver enhancers in episomal (EPI) and chromosomally (CHR) integrated contexts.
Three replicates of each experiment were made available to the public online, so we have compared the datasets individually and to each other according to the statistic used in their study: Ôi = log 2     ni c=1 Rc ni c=1 Dc    . (9) For this dataset we did not place a filter on the minimum number of counts per barcode, as the smallest number is already seven. We combine the replicates at the clonal level by first normalizing each replicate by dividing by the total number of counts of RNA or DNA, then combining replicates by concatenating the normalized data files. Thus, the combined dataset effectively behaves as if it has a greater number of measured clones.
Single Perturbations: Contribution of individual sequences to variation in outcomes. Experiment-wide measures of noise. For each data set we estiamte the total b 2 and provide all results in Table  1. For the fast-UTR data, the results show that a large fraction of the variance in these datasets (≈50%) is explained by sequence variation, which means accurate predictors of gene expression should be possible to generate for these datasets. However the noise might be too high to establish accurate estimates of activity at the single-sequence level, and deeper sequencing could improve the signal. For the Inoue data sets, the b 2 is approximately 90%, implying a very low contribution from technical noise and only a modest benefit could be obtained form deeper sequencing. GC-content biases aggregated outcomes.
We also estimated the fraction of variation caused by GC content alone, by residualizing at the sequence level and determining the change in variance. To do this we fit a 5th-order polynomial to the RNA/(RNA+DNA) ratio as a function of GCcontent and subtract it from the gene expression (see supplement for further details). In general, we find that more GC leads to a decrease in gene expression, which runs counter to the trend of the AU-rich elements that this dataset is enriched for, indicating that GC-content is a separate effect from AU-mediated decay. Performance of computational prediction methods. Several methods have been proposed to predict the outcome of MPRA experiments directly from sequence  [25, 26, 27]. The performance (squared correlation between predicted and measured values) is bounded above by b 2.
We trained a k-mer lasso model  [26]  in the fast-UTR data set, which takes sequence as its input, and outputs a prediction of the RNA/DNA ratio (see Methods). The squared-correlation between the out of sample prediction and measured outcome was at most 0. 26. Comparing these percentages to Table  1 , we expect that this approach can be substantially improved with better methods and/or less noisy data. Similarly, in the Inoue data sets, we find that the b 2 of the datasets are much higher than the squared correlation of their prediction methods (Pearson R 2 of 0. 34 on sequence data and 0. 36 when combining chromatin annotations with sequence information), which indicates that most of the sources of sequence-related variation remain unexplained by the current set of predictors.
Fraction of Expression Variance Explained Dataset b 2 Fraction Explained by GC fast-UTR Jurkat @ T 0 0. 56250 ± 0. 00003 0. 165 fast-UTR Jurkat @ T 4 /(T 4 + T 0 ) 0. 5677 ± 0. 0001 0. 057 fast-UTR Beas2B @ T 0 0. 41592 ± 0. 00007 0. 095 fast-UTR Beas2B @ T 4 /(T 4 + T 0 ) 0. 5641 ± 0. 0002 0. 119 with stem length produces a slope that differs from zero with p = 7. 5 × 10 -7 , r 2 = 0. 10, and slope = -0. 063. Using MPRAudit, we can examine the full set of sequences with CDEs and determine that b 2 = 0. 8598 ± 0. 0002, so 86% of the total variance is determined by differences in sequence. Accounting for categories of CDE stem lengths, MPRAudit determines that the within-group explainability is b 2 within = 0. 8433 ± 0. 0009. The fraction of within-group variation explained by sequence variation has barely changed as a result of grouping by CDE stem length, because the number of sequences affected by these groupings are small in number.
This suggests that more novel biology is likely to be uncovered within the former group of sequences than the latter, even though the latter have the more active elements. The fact that b 2 happens to be much higher within these particular groups than the dataset average given in Table  1  suggests that these are particularly active regions and that there are further rules for CDE activity that are waiting to be uncovered. In the fast-UTR data we define the effect size of a mutation to be the signed difference in gene expression between a sequence and its mutant. To estimate the variance contributed from the mutations we run MPRAudit over the signed differences of paired sequences. If the mutations have no effect then the estimate of b 2 will be 0. Of the 41,288 distinctly designed sequences in the dataset, 17,066 sequences are deliberately designed mutations of a known motif, which are present in 9,754 controls (many of the controls have multiple mutations).
Discussion In this work we demonstrate a method to determine the fraction of variance in a molecular assay that is due to sequence variation as opposed to technical artifacts or statistical noise. The purpose of this method is to help researchers who use MPRAs to know where to target their analysis efforts in different types of datasets as well as characterize the results of their existing analyses. Using realistic simulations of two types of MPRA data, we showed that MPRAudit is unbiased and calibrated provided a modest minimum clonal and sequencing depth is achieved. Application to real data demonstrates how MPRAudit can determine additional biological and experimental information about the sources of variation within a dataset. In concurrent work  [37] , MPRAudit helped us to discover novel rules for sequence motifs.
[38]  Applications may be found in relation to pathway analysis or genes or annotations of interest. MPRAudit does not rely on parametric model assumptions and is broadly flexible with respect to the quantity being measured. We rely on jackknife resampling rather than some sort of weighted linear model because the paired counts of RNA and DNA make it difficult to calculate the proper weights for each clone, whereas this information is automatically incorporated in the jackknife. MPRAudit relies on only a single replicate experiment, which mitigates confounding batch effects, saves time and money, and is better at handling sparse datasets than measuring correlations between replicates. In the results section we note that there is a relationship between the correlation between replicates and the explainability: the squared correlation between replicates should approximately equal the product of the explainabilities; a short proof is given in the supplement.
If the correlation between two replicates is low, the correlation between them is not capable of determining whether one replicate is less reproducible than the other. It should be kept in mind that the correlation between signal and noise may be nonzero when the counts of RNA or DNA are very low, leading to nonsensical explainability estimates and the breakdown of our additive model in equation 2. It should also be noted that since the technical and total variances are estimated separately, when the explainability is truly zero the estimate of the explainability intuitively should be negative for roughly half of all measurements. Resampling is most effective when the number of samples is large. This can be a problem for both the number of clones per sequence for V ar( ), and the number of sequences in the dataset for V ar( Ô).
To address the former issue, we suggest placing a filter on the minimum number of clones. In our analysis we required at least 5 clones per sequence, but a higher cutoff would give greater accuracy. With regards to the number of sequences in the dataset, we acknowledge that MPRAudit performs best on a massive dataset such as MPRAs. Sparse datasets can also cause problems with subsampling or resampling, since each of the jackknifed subsamples must be well-defined. If O i is given by log( RN A DN A ), each jackknifed subsample must have a nonzero number of RNA counts and DNA counts. For sparse datasets a statistic such as log( RNA+1 DNA+1 ) might be worth considering. This will not be an issue for RN A RN A+DN A , which is well-defined as long as the total counts are nonzero, but it might become a problem in a ratio of ratios such as RN A1/(RN A1+DN A1) RN A0/(RN A0+DN A0)+RN A1/(RN A1+DN A1).
Missing data is also an important quality control issue. While having very few counts per clone is not in theory a limitation of MPRAudit as long as the number of clones is high, we find that requiring at least one count of RNA and DNA per sequence gives us far greater repeatability across replicate datasets in practice due to dropout  [39] ; the set of missing sequences has a low correlation between replicates. These cases produce erroneously low predictions of technical variance, since resampling will repeatedly produce the same RN A RN A+DN A ratio, predicting zero technical variance while creating the largest possible amount of variance in practice. The convergence and bias of various resampling techniques will differ for different methods. For instance, in our analysis of simulated data, using log( RN A DN A led to greater bias than RN A RN A+DN A.
These considerations should be weighed by the user and investigated for their particular use case. Since MPRAudit estimates the variance of the noise of each measurement, we anticipate that these techniques will be directly applicable to many related methods that can be applied to MPRAs. For example, sampling-based variances can be used to provide the weights to regressions, ANOVAs, Bayesian analyses  [14] , and other methods that accept precision weights in heteroskedastic models, improving the power of bottom-up MPRA analyses. The code that we have made freely and easily accessible gives some implementations of MPRAudit, and can be easily extended by altering just a few lines. Users can make extensions of MPRAudit to consider arbitrary choices of statistics, such as complex ratios of ratios or regressions of time series data.
Acknowledgments DAS, AB and NZ were funded by  NIH  grants  K25HL121295 ,  U01HG009080 ,  R01HG006399 ,  R01CA227237 ,  R03DE025665 ,  R01CA227466 ,  R01ES029929 ,  R01GM110251 ,  R01HL124285 , and  DoD  grant  W81XWH-16-2-0018  DJE and OLT were funded by  NIH  grants  R35 HL145235 ,  R01 GM110251 , and  R01 HL124285. DAS would like to thank  Andrew Dahl ,  Joel Mefford , and  Nadav Ahituv  for useful discussions.
Methods Modeling the fragmentation phase. Numerical simulations of asteroid disruptions, including both the fragmentation phase during which the asteroid is broken up into small pieces and the gravitational phase during which fragments may reaccumulate due to their mutual attractions and form rubble piles, were first conducted in the early 2000 s and successfully reproduced the size distributions of asteroid families  18. The simulations, which used a smooth particle hydrodynamics (SPH) hydrocode including a model of brittle failure  43  to compute the fragmentation phase and the parallel N-body code pkdgrav  44  to compute the gravitational phase, showed that all fragments larger than typically 200 m are rubble piles formed by reaccumulation of smaller pieces. Other simulations using a porous material model in the SPH simulations to simulate the disruption of microporous asteroids  45  showed the same role of the reaccumulation phase in systematically forming rubble piles for bodies larger than 200 m  46.
Another improvement in the modeling of the reaccumulation came from the implementation of the SSDEM in pkdgrav  21. It was first used to simulate successfully the formation of the bilobate shape of the comet 67P by disruption and reaccumulation  57. SSDEM allows the explicit modeling of the contact forces between particles, parameterized by various friction parameters adjusted to represent realistic material types. SSDEM is thus the most realistic method to compute the reaccumulation, but it comes with a higher cost in computational time. A technique was recently developed and demonstrated to execute the handoff of the output of the fragmentation phase computed by SPH simulations to the input of the reaccumulation phase computed by pkdgrav  20. With this technique, which is the one used in this paper, it is possible to reconfigure the output of SPH simulations in order to use the more computationally expensive SSDEM code to model more accurately particle collisions and compute the reaccumulation process more efficiently.
Acknowledgements This material is based upon work supported by  NASA  under Contract  NNM10AA11C  issued through the  New Frontiers Program  and is also supported by  JSPS International Planetary Network. P. M. acknowledges funding support from the  French space agency CNES , from  Academies of Excellence: Complex systems and Space , environment, risk, and resilience, part of the  IDEX JEDI of the UniversitéCote d'Azur  and, with M. J. , from the  European Union 's  Horizon 2020 research and innovation program  under grant agreement No. 870377  (project  NEO-MAPP ). R. L. B. acknowledges  NASA  grant  NNH12ZDA001N , which supported development of some of the computational Methods used in this research. M. H. acknowledges  NASA  grant  NNH17ZDA001N. We are grateful to the entire OSIRIS-REx and Hayabusa2 teams for making the encounters with Bennu and Ryugu possible.
Code availability The simulations of the fragmentation phase are performed with a custom SPH code for which the details can be found in refs. 19, 45, 47 (and references therein). The simulations of the reaccumulation phase are performed by the code pkdgrav in its custom version that includes the Soft Sphere Discrete Element Method, for which the details can be found in refs. 20, 56 (and references therein). We have provided a compiled version of pkdgrav as well as the necessary files required to reproduce the outcome of Case 1, as an example. The input data for other cases presented in this paper are available on reasonable request.
The modulation of gene expression is a complex process occurring at various levels from which the transcriptional regulation is the core control code  (Petrillo et al. , 2014). The transcription machinery is regulated by an interplay between DNA-binding proteins called transcription factors (TFs) and cis-regulatory elements (CREs). TFs bind short sequences known as TF binding sites (TFBS) or motifs located at CREs (e. g. , promoters, enhancers, silencers). TFs may act as either activators or repressors of gene expression, leading to dynamic changes of the cellular pathways. For peach, annotation of TFs is available in the plant transcription factor database (plantTFDB)  (Tian et al. , 2019). As of February 2020, plantTFDB v5. 0 stores 2780 peach TFs classified into 58 families (http://planttfdb.
The proposed approach relies on three fundaments, i) an accurate definition of coexpressed gene modules, ii) an assessment of the effect of upstream region length regarding the effectiveness of motif discovery and, finally iii) disclosing the effect of splitting the analysis around the TSS site in discovering potential cis-elements. All together, we demonstrate the utility of our strategy in analyzing genome wide data to provide insights on gene regulation dynamics across tissues and specific conditions. To the best of our knowledge, no work has been reported on cis-elements present in P. persica on genome wide level, hence the originality of our survey. Additionally, the predicted motifs from this study can be browsed at (https://eead-csiccompbio. github. io/coexpression_motif_discovery/peach/), where we provide readers with direct links to the results, source code and a Docker container to reproduce the analysis on any other plant species.
Detailed scanning results can be accessed at https://eead-csiccompbio. github. io/coexpression_motif_discovery/peach. On the other side, (AT) repetitive elements were also scanned to check their relevance, e. g. , whether they correspond to the TATA box. The underlying data included in Figure S 2 , showed that TFBSs of these motifs were remarkably distant to the TSS and were distributed across the whole proximal region.