dataesr/scientific-paragraphs-categorization huggingface dataset:
Since we assume that the cells are independently distributed, we can average the probabilities for the hypotheses across a large set of randomly sampled cell pairs, one from each subpopulation. The Bayes factor from the averaged probability will provide an estimate of whether cells from one subpopulation tend to express g at a higher frequency. Notably, we draw here data from the prior distribution and not the posterior for given cells. As a consequence, these Bayes factors can be approximated in a unbiased fashion using a naive Monte Carlo estimator. We noticed in the case of the real dataset that the aggregate posterior on u might not perfectly match the prior for rare cell types. Consequently, we replaced the prior by the aggregate posterior for all the analysis in this manuscript. Scalability From a scalability perspective, our methods scales to a million cells in less than an our as reported in our original manuscript  [29].
Remarkably, our inference procedure is scalable with respect to the number of batches (with number of floating point operations linear per dataset per iterations) while even the most efficiently implemented method such as Scanorama needs to compare each dataset pairs (quadratic). Our runtime on the large dataset from Scanorama is 42 min for 250 iterations over 100K cells from 26 datasets, still 30 min longer than their algorithm. We expect that further work on implementation details such as GPU usage efficiency, memory loading, data format and iteration monitoring will help close the computational gap between these algorithms.
Datasets We report an extensive list of datasets at Supplementary Table  1. For all UMI based datasets we took the raw counts without any normalization as input to scVI. Gene Selection A common practice in data harmonization is to perform gene selection prior to harmonization. This assumption is critical when the number of genes that can be taken into account by the algorithm is small and potentially biological signal could be lost. scVI is however designed for large datasets which do not fall into the high-dimensional statistics data regime  [29]. Remarkably, there is no need for crude gene filtering as part of our pipeline and we adopt it as part of this publication only for concerns of fairness in benchmarking. For real datasets, we calculated the dispersion (variance to mean ratio) for all genes using Seurat in each dataset and selected g =1, 000 genes with the highest dispersion from each.
The performance of scVI is not as affected by gene set and we use the same gene selection scheme as in  [13]  to ensure fairness in our comparison. We then took the union of these gene list as input to Seurat Alignment, MNN and scANVI. One exception is the differential expression study for which we kept the gene set (g =3, 346) to have it match the bulk reference as in  [29]. Cell type labeling for the Tabula Muris Dataset For the Tabula Muris dataset, cell types are defined by first reducing the dimensions of the data by principal component analysis and then performing nearest-neighbor-graph-based clustering. The labels for Smart-Seq2 and 10x data are derived independently. All cells in both dataset are labeled, but there is also a possibility that they are mislabelled since the labels are computationally derived.
Since cells used in Smart-Seq2 are first FACS sorted into each plate, some cell types might have been lost during the sorting process, resulting in incomplete overlap in cell types between the two datasets. Hierarchical cell type labeling for the mouse nervous system dataset The multilevel labels are generated through an iterative process that is described in detail in the original publication  [60]. The clustering was performed with strict quality filters, takes into account anatomical information and were validated at different levels using existing scRNAseq dataset, osmFISH, RNAscope and others. The cell types taxonomy is derived differently for each level and the details can be found in the original publication. Cell type clusters were obtained by Louvain clustering on a multiscale kNN graph and DBSCAN.
The first level separates neurons and non-neuronal cells. The second level separates peripheral neuronal system from central neuronal system. The third layer separate anterior posterior domain, and the fourth layer is split by excitatory versus inhibitory neurotransmitter. At this level, all cells are divided into 39 subsets, each corresponding to a coarser cell type definition. Then, within each subset the authors defined N=28 enriched genes and used linkage (correlation distance and Ward method) to construct the dendrogram. Normalization of SmartSeq2 data For the MarrowMT-ss2 dataset, we normalized the read counts per gene by relative transcript length (average transcript lengths of a gene divided by average gene length over all genes), and subsequently took the integer part of the normalized count.
Maximum Posterior Probability We evaluate the performance of the scANVI classifier at transferring labels from an annotated dataset to an unannotated dataset by looking at the maximum posterior probability for the observed classes. By default scANVI classifier sets the number of classes to the same number of cell types in the merged dataset. In the case of N observed labels from the annotated dataset and one unannotated dataset (thus the cell type label is "Unlabeled") scANVI assumes N + 1 classes. For each cell, scANVI assigns a posterior probability for each of the N + 1 classes. The maximum posterior probability for the observed classes is the highest probability of a cell being assigned to one of the N observed classes.
Signature calculus To compute the signature of a cell, we followed the normalization procedure from  [24]  which consists in dividing by total numbers of UMIs, applying a entry-wise transformation x 7 ! log(1 + 10 4 x) and z-score normalization for each gene. Then, we aggregated over the genes of interest for each cell by applying the sign from the gene-set and averaging. that there is a common signal exactly reproducible between the two datasets and that CCA capture most of the biological variation, which is not obvious true considering limited suitability of linear Gaussian model for scRNA-seq. Seurat anchors  [20]  relies on CCA and suffer from the same problems as MNNs with its specific normalization scheme. Finally, the recent LIGER  [18]  method at first sight seems like a non-probabilistic of scVI since it also learns a degenerate conditional distribution via its integrative non-negative matrix factorization  [18]  (NMF is a noisy-less version of a Gamma-Poisson generative model).
Figure8: Supplementary study of harmonizing datasets with different cellular composition. We show here the case where each of the two datasets has a unique cell types and share all the others. For each box plot, we report over all the possible combinations of left-out cell types. (a) Batch entropy mixing on the unique populations (lower is better). (b) Conservation of the original structure across all populations (unique and non-unique; higher is better). (c) Batch Entropy mixing on the non-unique populations (higher is better). The boxplots are standard Tukey boxplots where the box is delineated by the first and third quartile and the whisker lines are the first and third quartile plus minus 1. 5 times the box height. The dots are outliers that fall above or below the whisker lines.
Additional metric for retainment of structure via k-means. For scANVI we perform semi-supervision using the cell type label (not k-means cluster labels) from only one of the two datasets. Thus we train two separate models SCANVI1 and SCANVI2. To obtain a measure of clustering conservation, we first run k-means clustering in the latent space of dataset 1, then in the harmonized latent space using only cells from dataset 1. We compute the adjusted Rand Index of the two clustering results. We then do the same for dataset 2 and the final score is the average for both datasets.
This may complicate analysis leading to spurious results. Some methods seek to correct the noise, sometimes also leveraging combinations between sets  (20). In particular, matrix factorization methods are quite effective although costful on such noisy data  (21). However, using itemset mining to find combinations of interest based on a criterion and assessing their enrichment are two different approaches, which are worthwhile to now be combined. This paper proposes a method named OLOGRAM-MODL to leverage both, by calculating the significance of mined combinations of overlaps of interest. OLOGRAM is a statistical framework that evaluates the statistical enrichment of a given combination of sets overlapping, while MODL is an itemset mining algorithm used to select the most representative combinations.
However, querying overlaps has a complexity of O( n i ) where n i is the number of regions in the i th set. This is in contrast with interval trees (used in BEDTOOLS fisher) whose complexity is O(log n) per queried region, simplifying to O(n 2 log n 1 ) for only two sets. Computing the combination enrichment. For any combination ␥ , we fit a Negative Binomial model on each measure N(␥ ) and S(␥ ) (Definition 4), and give the P-value for the actual observed occurrences to happen by chance. By default, the approach computes the enrichment of all combinations observed in the real data. To perform a selection among those instead, the MODL algorithm (see MODL itemset mining algorithm section) is proposed. The user can also provide a custom selection, as the intersection matrix itself can be accessed in our Python API to help make that selection.
In OLOGRAM-MODL, it can be used to pre-select combinations of sets which are of interest, hence reducing afterward the total number of enrichment computations and making the results easier to interpret. Algorithm 1 Multiple Overlap Dictionary Learning (MODL) algorithm for combination mining Extraction of relevant sub-complexes with dictionary learning. In the OLOGRAM-MODL approach, the input matrix of MODL is the matrix of overlap flags provided by the algorithm in Figure  2. It has one row per overlap and one column per set in the real, non-shuffled data. However, any matrix matching this format can be used. Principle of itemset mining via dictionary learning. Dictionary learning is a special case of matrix factorization, where the goal is to learn the sparse code U and the dictionary V from the data X under certain constraints, minimizing the reconstruction error.
Implementation and availability of the method As an update of OLOGRAM, the code is written in Python 3, with some performance-critical tasks in Cython and C++. To preserve RAM when working with large files, the shuffles are divided into batches. However, OLOGRAM must still remember all computed intersections within a run; hence, we also permit merging different runs as superbatches if necessary, although fitting is not assessed then. The examples take a timescale of minutes to run, up to hours for very large cases (see Running times section) on a 4-core 3. 5 Ghz processor. For the MODL subroutines of dictionary learning and sparse coding, the Scikit-Learn implementation is used  (27). OLOGRAM-MODL is accessible through the command line interface of pygtftk  (28)  which is available on Bioconda.
Data and results We use two main types of data to produce the results of this study. Two different types are used: either a set of genomic intervals as BED files, or a representation of the overlaps between those sets of intervals as a matrix X where each row is a combination. The details about the creation and the detailed sources of this data is presented in Supplementary Note S-IV. Artificial datasets are generated with a known ground truth, so we can use them to demonstrate and quantify the behavior of our algorithms. Those are used in the results presented in section Artificial data and comparison to existing approaches section. For the artificial matrices, they have a large number of rows to perform a large scale experiment and ensure the results are not due to random chance, and the itemsets selected represent biologically interesting cases, with a diversity of complexes and two overlapping ones.
Conversely, we also use three biological datasets to produce the results of Biological results section: human and murine transcription factor binding sites ChIP-Seq data, as well as human sc-ATAC-seq data. The two first datasets serve to demonstrate OLOGRAM's utility in discovering complexes of regulator, and the third demonstrates how it can scale to larger problems. When selecting the transcrip- tional regulators to study, we tried to avoid cherry picking by using a wide variety of TRs, and including many commonly studied ones. We present the full list in Supplementary Data. The complete workflow used to produce the Results of this study, along with the full data used, is available as a Snakefile (Snakemake file, (  29 )) at https://github. com/qferre/ologram-modl supp mat. Users can re-use and adapt the commands presented there for their own workflow.
RESULTS In summary, OLOGRAM-MODL takes an input a set of genomic regions (usually in BED format) and will determine for each encountered combination of genomic regions whether it is enriched, meaning whether it is encountered more than expected by chance. The MODL algorithm can be used to restrict the combinations for which this is calculated to combinations of interest (custom selections are also possible). The first goal of the experiments presented in this section is to validate our statistical model and our itemset mining algorithm. We do so using artificial data for which the ground truth is known, and the results can be compared to. We also compare to existing approaches. In the second part, we show that OLOGRAM-MODL is not only able to deal with true biological data but also a fair way to discover relevant complex of genomics regions and get novel insights.
Artificial data and comparison to existing approaches First, we want to establish that our statistical model is valid in the context of combination enrichment. We also want to show that the MODL algorithm can find relevant itemsets in situations representative of real use cases. To do so, we benchmark OLOGRAM-MODL against data with a known ground truth to which the results can be compared. The datasets used will refer to the data defined in Supplementary Note S-IV, see this note for more details. S1  and in the rest of this section. We show that OLOGRAM-MODL is easy-to-use and brings novel insights compared to existing approaches.
Running times. A full run of 1000 shuffles in the large 'coarse' artificial dataset took 40 s on a laptop (i7-7820HQ processor, single thread), against 10 s for MULTOVL. Typical runtimes are in the order of a few minutes, such as for the example in Regulatory complexes in murine promoters section. In extreme cases, such as the very large data use case of Overlap between single-cell ATAC-seq in PBMC data section with 75+ large files, the runtime was 1. 5 h per 4 shuffles, mostly due to the analysis in steps (ii) and (iii) (as defined in Time scaling section) and not the shuffles themselves. In general our running times are of similar orders of magnitude compared to existing tools. For all these cases, detailed runtimes, data characteristics and computer specifications are presented in the legends of the Supplementary Figures.
Other algorithms give them a lower ranking, because if the rule ABC is true the rules AB and ABC are equally true. Itemset miners in general are vulnerable to noise. Approximate itemset miners (such MODL purports to be) are proposed as a solution  (25). Indeed, when using a clusteringbased approximate itemset miner (CL-MAX, (  30 )), results that are close to MODL's. However, CL-MAX is more vulnerable to a poor choice of hyperparameters. We also compare OLOGRAM with GINOM on their demonstration data. There is a link between the fold changes in OLOGRAM and in GINOM, as shown in Supplementary Table  S4. The combinations selected by the two approaches differ somewhat since MODL is not based on enrichment, but MODL's selection is still relevant, and there are commonalities between the two. MODL behaves as would be expected of a frequent itemset mining algorithm, in that it seeks to rebuild the original data matrix and does not select sets alone but with their co-localized sets (i.
Biological results Having demonstrated that OLOGRAM-MODL's behavior matches our expectations, we now apply the tool to a variety of biological problems to demonstrate that it can present relevant insights in different use cases, while remaining easy to use. Nevertheless, the rigorous validation of our approach rests on the results obtained with artificial data, where the ground truth is known. This section on biological data should rather be seen as illustrative of the possibilities of the approach when applied to biological problems. Regulatory complexes in murine promoters. We begin with a fundamental use case for OLOGRAM: the analysis of regulatory complexes. We perform a search in murine promoters, by looking for the transcriptional regulators (TRs) among our selected TRs that are associated with three different query TRs: CTCF, NANOG and IRF1.
DISCUSSION The OLOGRAM-MODL approach consists of two steps. We combine an optional itemset mining algorithm to find interesting combinations, with a statistical model to determine the enrichment of the relevant combinations, asserting whether this combination occurs in the real data across more base pairs that would be expected by chance. The integration with the pygtftk toolkit facilitates complex queries in bioinformatic worfklows. As a command line tool whose dependencies can be handled by conda, it is convenient to run on clusters with reproducible workflow managers such as Snakemake. Thanks to the use of a Monte Carlo approach and to C++ optimization, OLOGRAM-MODL can even be run on a laptop.
Perspectives OLOGRAM-MODL is applicable to any problem that can be reduced to quantifying the significance of overlaps between n sets of position intervals. Besides epigenetic marks binding sites, associations between sets of regions such as 'promoters of housekeeping genes' or 'Binding sites for the regulator X in the experimental condition Y' can also be integrated. This is especially important in our 'big data' and multi-omics era, as combinations can also be combinations of datasets. The core idea of MODL, selecting itemsets according to how well they rebuild the ensemble of all itemsets, is precedented  (42). The MODL algorithm can be applied to any submodular problem, as the API supports custom error functions. For example, since variable selection in Naive Bayes classifiers is indeed submodular  (31) , MODL could select combinations that help such a classifier predict the query, as we demonstrate in Results and comparison section.
An individual's brown chroma score was obtained as the average of these four points. This method correlates well with observational assessments using colour chips (r = -0. 78, p <. 0001, Dreiss & Roulin, 2010)  and has high repeatability (97. 6%; Antoniazza et al. , 2010). Brown chroma data from  Burri et al. (2016)  were used to complete the data set, using the same individuals as for the MC1R analysis (Appendix S1). Given the marked non-normality of the data, a nonparametric Kruskal-Wallis test was performed to detect differences in coloration between the six populations. Further, a pairwise Wilcoxon rank sum test was used to identify significant differences between pairs of populations using a Bonferroni correction.
| Population structure and genetic diversity To investigate population structure among our samples, sNMF v. 1. 2 (Frichot et al. , 2014)  was run for K 2 to 6 in 25 replicates to infer individual clustering and admixture proportions. For this analysis, singletons were excluded and the remaining SNPs were pruned for linkage disequilibrium (LD) with PLINK v1. 946  (Purcell et al. , 2007;  parameters -indep-pairwise 50 10 0. 1) as recommended by the authors, yielding 319,801 SNP. The same data set was used to perform a principal component analysis (PCA) with the R package SNPRelate  (Zheng et al. , 2012). Treemix  (Pickrell & Pritchard, 2012 ) was used to infer population splits in our data, using the LD-pruned data set further filtered to include no missing data (180,764 SNP). To detect meaningful admixture between populations, 10 replicates were run for 0 to 8 migration events, with the tree rooted on the PT population, representative of the glacial refugium.
An extra run without migration events was conducted with a North American owl as an outgroup in the data set to verify that the root did not affect the topology of the tree. To estimate population statistics, individuals found to be mis-assigned to their given population based on genetic structure analyses (PCA and sNMF) were removed so as not to bias allelic frequencies (n = 3 individuals from Ireland). Individual expected and observed heterozygosity and population-specific private alleles were estimated from the full data set (6,721,999 SNP) using custom R scripts for each genetic lineage identified by sNMF with K = 4 (Figure  1c ). To account for differences in sample sizes, private alleles were calculated by randomly sampling 9 individuals from the larger populations (GB and western Europe) 10 times in a bootstrap-fashion and estimating the mean.
| Migration surface estimate The estimated effective migration surface (EEMS) v. 0. 0. 9 software  (Petkova et al. , 2016)  was used to visualize geographic regions with higher or lower than average levels of gene flow within our data set. The provided tool bed2diff was used to compute the matrix of genetic dissimilarities, from the data set pruned for LD produced above. The free google Maps api v. 3 tool (http://www. birdt heme. org/usefu l/v3tool. html) was used to draw the polygon outlining the study area in western Europe (Figure  2a ). EEMS was run with 750 demes in three independent chains of 5 million MCMC iterations with a 1 million iterations burnin. Results were checked for MCMC chain convergence visually and through the linear relation between the observed and fitted values for within-and between-demes estimates using the accompanying R package ReeMsplots v.
| Treatment and analyses of capturerecapture data In addition to genomic data, recapture data of ringed barn owls across Europe were obtained from the EURING database  (obtained in March 2020; du Feu et al. , 2016; Speek et al. , 2001). Specifically, we estimated the frequency of crossings over open water (flights between two land masses over a water body) between GB and central and western Europe, as well as between GB and Ireland. To do so, we kept records of birds that had been recaptured at least once after ringing (n = 94,797 recaptures, n = 80,083 individuals, from 1910 to 2019) and filtered the accuracy of the "time of capture" parameter to a period of within 6 weeks of the reported date to exclude potentially unreliable data points. We extracted the number of birds ringed and recaptured in GB and Ireland, as well as in the countries that produced or received migrant birds from these islands and western Europe (Belgium, Denmark, France, Spain, Germany, Switzerland and The Netherlands).
| Post-glacial species distribution To support the demographic scenarios tested in the following section, we modelled the past spatial distribution of barn owls in western Europe, in order to identify the regions of high habitat suitability at the Last Glacial Maximum (LGM, 20,000 years BP). A complete description of the models can be found in Methods S1. Briefly, using maximum entropy modelling (MaxEnt), a presence-only based modelling tool, we built species distribution models (SDM) for the Western Palearctic (Figure  S1 ) based on climatic variables extracted from the WorldClim database (Figure  S2 ;  Hijmans et al. , 2005)  at 5 arc min resolution. The best combination of feature and regularization multiplier based on the corrected AIC (as recommended by  Warren & Seifert, 2011)  was achieved with a quadratic model with 1 as regularization multiplier (Table  S6 ).