dataesr/scientific-paragraphs-categorization huggingface dataset:
These are then counted at specific genomic locations to effectively serve as molecular counters. The genomic locations typically originate either from reference databases (for gene-centric analyses) or are identified de-novo through peak-calling algorithms (for epigenetic or regulatory non-coding events analyses). The measured sequencing signal's intensity can also be highly diverse, ranging from a few thousand UMIs per observation in single-cell analyses to hundreds of millions of reads in bulk sequencing. However, a common aspect across all these whole-genome sequencing protocols is that the sequencing signal per observation ultimately gets quantified through counts in a large number of features (e. g. , genes or genomic peaks). We demonstrate that a generic methodology can be employed to conduct between-observation comparisons in a wide variety of experimental designs, with low to high observation counts, and shallow to deep sequencing.
High resolution consensus peak identification in large datasets When working with sequencing data that target regulatory elements, the signal typically manifests in the form of peaks in sequencing protocols such as ChIP-seq, ATAC-seq, or CAGE. Various peak-callers have been developed to identify these signal regions on genomic tracks  [23, 29]. To integrate observation-level peaks into consensus peaks, which serve as key sampling points for sequencing reads and downstream analyses, we propose a simple tool. This tool, which we used in previous work, accepts per-experiment peak-calling results (from an external tool) as input and outputs consensus peaks based on the genome-wide peak density. While identifying consensus peaks based on simple overlaps is an acceptable approach when working with a small number of experiments, it becomes problematic with a larger number of experiments.
r a i,j = K 2/3 i,j × H(-α reg j × K i,j ) -µ 2/3 i,j × H(-α reg j × µ i,j ) (µ i,j + α reg j × µ 2 i,j ) 1/6 Where H (x) = H 2 F 1 (2/3, 1/3, 5/3, x), H 2 F 1 being the Gauss Hypergeometric function. We found Pearson residuals to be much more sensitive to outliers, and as they are still skewed and heavy-tailed due to the linear nature of the transformation, their values have to be clipped to prevent large values from driving a majority of the variance. The clipping value is a sensitive hyperparameter and has to be adjusted per dataset, with too large values being sensitive to outliers, and too small values removing biological signal. Additionally, SCTransform uses an arbitrary lower bound on the variance estimate to prevent weakly expressed genes to drive a majority of the variance. We found that Anscombe residuals do not require clipping (Figure  S2 ), give low weights to weakly expressed features (Figure  S4  A), as well as having better known theoretical statistical properties concerning normality.
Differential expression We implement DESeq2 with Log Fold Change (LFC) shrinkage  [30]  via an rpy2 wrapper to conduct differential expression between two conditions, as it uses a statistical model similar to ours, supporting a per feature, per observation normalisation factor. However, DESeq2 can be computationally intensive for large datasets. When working with more than 50 samples, which provides sufficient statistical power, we transition to a Welch's t-test performed on the NB residuals, which has been shown to work well for single-cell data  [2]. It should be noted, though, that the t-test offers less statistical power than DESeq2.
Optimal number of components for principal component analysis To automatically identify the optimal number of Principal Components, we utilize Horn's Parallel Permutation Analysis, which has proven to be one of the most effective methods to determine the number of components in factor analysis  [1]. This approach involves generating row-wise permutations for each feature, computing PCA on these permuted datasets, and selecting the number of components at the threshold where the eigenvalues from the randomized dataset exceed those from the actual dataset. By default, we use the residuals as input to the PCA, and due to the computational cost of this approach, we only perform three permutations by default. This is justifiable, as the randomized eigenvalues are very stable on large matrices (Figure  S4 D ).
Data visualisation For visualising the similarity among observations in datasets with a large number of features and observations, dimensionality reduction methods such as UMAP (Uniform Manifold Approximation and Projection) or t-SNE  [16]  (t-Distributed Stochastic Neighbor Embedding) have become standard, particularly for the analysis of single-cell data. Here, we provide a wrapper for UMAP, which is, by default, performed on the PCA representation. By default, Pearson correlation is used as the metric when there are more than 10 input dimensions; otherwise, the Euclidean distance is used.
Analysis of benchmark datasets We manually retrieved matching peak calling BED files, ChIP BAM files, and Input BAM files from ENCODE  [19]  (see Supplementary Tables). We selected only those H3K4Me3 ChIP-seq experiments that were mapped onto hg38, and removed samples with any audit errors. In cases where multiple replicates of the input file were present, we retained only the most deeply sequenced one and used it across all biological replicates. We used a coarse annotation (the same as in the ENCODE data portal) instead of the detailed cell types as there is a large number of precise cell types with no replicates. The list of genes associated with Gene Ontology terms was retrieved from the g:ProfileR website  [21]. MUFFIN covers multiple aspects of the analysis of several types of functional sequencing data, from the identification of consensus peaks, to differential expression and Gene Set Enrichment Analysis.
We obtained the filtered 10k PBMC scRNAseq count tables from 10X Genomics (https://www. 10xgenomics. com/resources/datasets/10khuman-pbmcs-3-ht-v3-1-chromium-x-3-1-high ). We also analyzed the dataset using Seurat with SCTransform V2 as a reference, with the default settings from its tutorial vignette. We compared the similarity between our method and SCTransform using the jaccard index of similarities between clusters, and using Adjusted Mutual Information (AMI) to compare similarities across all clusters. TCGA ATAC-seq  [3]  count tables were obtained from the NIH (https://gdc. cancer. gov/about-data/publications/ATACseq-AWG ). The list of candidate genes associated with specific cancer hallmarks was obtained from the CHG database  [28]. All analyses were conducted using the default settings, with the exception of the normalisation method: for the singlecell dataset, we employed scran normalisation; for the TCGA ATAC-seq dataset, we applied DESeq's median of ratios; and for the ChIP-seq dataset, we used our custom normalisation method.
Additionally, the suite offers high-resolution consensus peak identification for large datasets, integrating observation-level peaks into consensus peaks-a crucial feature when dealing with a larger number of experiments. Tools are also provided for linking genomic regions to genes and functional annotations, including utilities for renaming regions based on their nearest gene and performing Gene Set Enrichment Analysis. Furthermore, the suite encompasses a variety of tools for functional sequencing analyses, covering count modeling for between-observation comparisons, differential expression analysis, feature selection for downstream tasks, determination of the optimal number of components for principal component analysis, graph clustering, and data visualization  (1). A detailed description of the tools is available in the Methods.
To validate the suite's efficacy, it was tested in three largely different datasets : i) the ENCODE Epigenetic Atlas of Immune Cells using H3K4Me3 ChIP-seq, ii) 10,000 PBMCs sequenced with 10X scRNA-seq, and iii) 400 cancer samples sequenced via ATAC-seq originating from TCGA. These applications showcase the versatility and robustness of the MUFFIN suite across diverse biological and technical contexts.
ENCODE epigenetic atlas of Immune cells through H3K4Me3 ChIP-seq In order to evaluate the general applicability of our proposed methodology, we utilised the H3K4Me3 ChIP-seq dataset from the comprehensive ENCODE epigenetic atlas of immune cells. H3K4me3, an epigenetic modification known for marking active promoters, is a crucial mediator of transcriptional activity and plays a significant role in cellular differentiation and function. Here, multiple difficulties are linked to this dataset : there are no reference regions to sample sequencing signals from, the sample sizes are rather small, and the sequencing protocol uses a control track. Our first task was to identify which regions to consider for quantifying H3K4Me3 occupancy. For this purpose, we employed our density-based approach to identify consensus peaks from per-experiment peak calling results.
ENCODE cCREs ReMap ChIP-seq density identifies the consensus peaks at high resolution in dense regions, and correlates with the Transcription Factor binding density from the external ReMap database  [8]  (Figure  2A ). Subsequently, we counted reads in immuno-precipitated and input sample pairs at consensus peak locations. Ultimately, the counts are normalised using our centering and scaling approach (Methods), transformed using Anscombe residuals, and can be processed through a visualisation or dimensionality reduction routine. As shown in Figure  2  B, a heatmap view shows a clear distinction between immune cell types and reveals distinct H3K4Me3 occupancy across different cell types.
MUFFIN consensus peaks Next, we sought to identify regions differentially marked by H3K4Me3 between B-cells (n=5) and T-cells (n=26). We found that the regions most differentially bound in B-cells are indeed located near B-cell marker genes such as HLA-DOB, which is expressed in Antigen Presenting Cells, or PAX5, a transcription factor involved in B-Cell differentiation. Conversely, for T-cells, we also retrieved T-cell markers like CD247 (T-cell surface glycoprotein CD3 zeta chain) or CD3E (CD3 Epsilon Subunit Of T-Cell Receptor Complex). Using our tool of Genomic Regions Enrichment Analysis, we find that these differentially occupied regions are located near genes annotated with immune and lymphocyte-related Gene Ontology terms (Figure  2 C-D ). As a whole, this confirms that our methods are able to analyze ChIP-seq data and are likely to work well with other similar datasets that are using control sequencing.
TNFRSF4 ITGB1 GPR183 TRADD FHIT CD40LG ADTRP TSHZ2 SLC4A10 GZMK CXCR6 IL23R GZMH GNLY FGFBP2 TRGC2 GZMK CCL5 DUSP2 CD8B SMIM6 CLEC4C LRRC26 LILRA4 CD8B LINC02446 CD8A KLRK1 CLEC4E CFD CYP27A1 FTH1 MX1 ISG15 RSAD2 IFIT1 TREML1 CAVIN2 GP1BB SPARC CES1 PADI4 CYP1B1 KCTD12 CLEC10A FCER1A MRC1 PKIB APOBEC3A LIPA PID1 HLA-DPA1 CDKN1C FCGR3A HES4 MS4A7 KLRF1 SH2D1B SPON2 CLIC3 IGKC TCL1A IGHD AL139020. 1 IGLC3 IGLC2 TCL1A CD200 IGKC IGHA1 SSPN TNFRSF13B IGLC2 IGLC3 IGLC1 Cell-type clustering of 10,000 PBMCs sequenced by 10X scRNA-seq After investigating the differential occupancy of H3K4Me3 across various immune cell types, we sought to evaluate the robustness and versatility of our methodology using another widely adopted dataset. Hence, we turned our focus to a single-cell RNA-sequencing (scRNA-seq) dataset, which offers a distinct set of challenges due to the large number of observations, inherent sparsity, shallow sequencing depth and high dimensionality of the data.
The 10X Genomics scRNA-seq dataset of Peripheral Blood Mononuclear Cells (PBMCs), a mixed population of various immune cells, is often used as a benchmark case study. The PBMC dataset encompasses numerous immune cell types, allowing us to assess the efficacy of our approach in identifying and distinguishing between different cellular phenotypes. Here, count quantification per gene, per cell has already been performed. We used UMAP to visualise the transcriptomic similarities between individual cells and applied graph clustering to pinpoint cell types (Figure  3 A ). We clearly identify canonical cell types, such as Monocytes, B-Cells, T-Cells or platelets. Their markers, identified via the coefficients of a logistic regression, are also coherent (Figure  3   for platelets. Furthermore, we are able to identify cell types at a fine resolution, highlighting clusters of cells with similar transcriptomes, but with a few very specific markers, especially in the B-Cell clusters and the monocyte clusters.
The Cancer Genome Atlas open chromatin landscape of cancers via ATAC-seq Finally, we studied the open chromatin landscape in cancers using Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) in tissues. ATAC-seq is a technique used to assess the chromatin accessibility, thereby providing insights into the regulatory regions active in a particular cell type, tissue or condition. This sequencing protocol typically generates deeply sequenced datasets with an extremely large number of features. The read count table per ATAC peak, per sample was already provided, but a consensus peak identification as performed with the ChIP-seq dataset could have been performed. In Figure  4  A, using UMAP, we can clearly distinguish the different cancer types, and even distinguish between different breast cancer molecular subtypes.
Discussion In conclusion, our proposed framework, MUFFIN, represents a valuable and generic toolset for the processing and analysis of count data derived from a variety of high-throughput sequencing experiments. We have demonstrated its performance across diverse datasets, from epigenetic to transcriptomic, and from bulk to single-cell resolution. The underlying model is highly flexible and can adapt to a wide range of experimental designs, and relies on already well-tested, state-of-art statistical approaches. One of MUFFIN's key strengths lies in its ability to address multiple facets of count-based genomic assays in a generic way with minimal hand-tuning : normalization, count transformation, feature selection, dimensionality reduction, visualization, differential expression, and clustering.
Additionally, it offers specialized tools for analyzing data that do not rely on gene annotations. This includes a tool for identifying consensus peaks at a high resolution, and another one for performing functional enrichment of genes located near genomic regions of interest. Furthermore, MUFFIN seamlessly integrates in the existing Python Scanpy ecosystem and its diverse range of tools. The modular architecture of MUFFIN could also allow future packages to be built on top of it. However, it's worth noting a potential limitation of MUFFIN, which is its compatibility with extremely large, sparse datasets. As its count transformation does not support sparse data formats, this may restrict its use in situations where memory capacity is a significant concern, such as in the analysis of extremely large integrative single-cell datasets consisting of millions of cells.
Acknowledgements We would like to thank  Dr. Lionel Spinelli  for engaging in discussions that have contributed positively to our work. This work was supported with ; PhD Fellowship to P. D. L. from the  French Ministry of Higher Education and Research (MESR) ;  Institut National de la Santé et de la Recherche Médicale (INSERM) ; The  Core Cluster of the Institut Français de Bioinformatique (IFB)  ( ANR-11-INBS-0013 ) for granting access to its high performance computing resources. The results shown here are based upon data generated by the  TCGA Research Network , the  ENCODE Consortium  and the  ENCODE production laboratories.
Direct RNA sequencing data pre-processing Raw fast5 files were processed with the Master of Pores pipeline (version 1. 5,  (22) , https://github. com/biocorecrg/master_of_pores). Fast5 files were demultiplexed with DeePlexiCon using default parameters  (21)  and basecalled with Guppy basecaller v4. 0. (https://nanoporetech. com). Reads were mapped to D. melanogaster rRNA sequences, obtained from Ensembl: FBtr0346885 (28S), FBtr0346878 (18S), FBtr0346887 (5. 8S) and FBtr0086426 (5S). Mapping was performed using minimap2 v2. 17 (https://github. com/lh3/minimap2) with "-ax splice -k 14 -uf" options.
Analysis of RNA modifications in direct RNA sequencing data Basecalling features (base quality, mismatch frequency, insertion frequency and deletion frequency) of each 5-mer were extracted from 2 replicates of untreated and Fib KD S2R+ reads using Epinano (23), version 1. 2 (https://github. com/enovoa/EpiNano). Then, each position at the centre of the 5-mer was assigned a score that consisted of the difference between the scaled sum of the aforementioned frequencies in the untreated and in the Fib KD sample, as previously described  (24). Then, we calculated the median of scores for each rRNA transcript and each replicate. For de novo discovery of Nm sites, nucleotides with a score greater than 3x the median score of all positions in the same transcript in both replicates or 5x the median in one of the replicates were further kept.
Identification of Nm sites on Drosophila melanogaster rRNA Different rRNA annotations are available in Drosophila melanogaster including multiple ones from the Ensembl 103 assembly, one GenBank reference that was used for previous snoRNA predictions, and a PDB reference that was used for 3D representations. To facilitate the navigation between previously predicted and novel identified Nm positions (see below), we aligned annotated rRNA sequences with Sanger-sequenced data from both embryonic cell line S2R+ cells and wildtype Canton-S flies (Supplementary Table  1  and Supplementary Figure  1A ). The most similar sequences were FBtr0346885 (28S), FBtr0346878 (18S), FBtr0346887 (5. 8S) and FBtr0086426 (5S) from the Ensembl assembly. Consistent with previous observations  (33) , we confirmed that the 28S rRNA has a "hidden break" based on the lower read coverage that spans nucleotides 1814 to 1858, but also based on the presence on the electrogram of two RNA peaks just below 2 kb that reflects a mixture of 18S and two similarly sized fragments of 28S (Supplementary Figure  1B-C ).
Finally, as expected, the 3 Nm sites that were Fib-independent in the RMS dataset and missed by the de novo approach were also missed by the validation approach (Figure  2A ). Overall, RMS ΔMethScores and Epinano scores correlated significantly (Pearson r=0. 58, p<0. 001), although the amplitude of RMS ΔMethScores was larger than that of Epinano scores (Figure  2D  and 2E ). We further noticed that nanopore was more sensitive to differential methylation on Am sites in comparison to RMS, which had no visible bias towards any base (Figure  2E , Supplementary Figure  4C  and 4D ). In conclusion, while 77% of RMSidentified Nm sites were de novo discovered with nanopore sequencing, 95% of Fibdependent RMS-identified sites could be validated with nanopore, indicating that the nanopore approach is a bona fide validation tool of Nm modifications on rRNA.
Data availability Basecalled FAST5 of nanopore direct RNA sequencing runs have been deposited to ENA, under accession code PRJEB45722. Tracks for IGV visualization of the mapped reads (BAM files), predicted rRNA modifications in D. melanogaster (based on homology to known human rRNA modifications), EpiNano RNA modification scores, as well as the code to perform EpiNano score analysis per k-mer and de novo discovery of Nm sites, can be found in GitHub (https://github. com/novoalab/Nm_Nanopore_Drosophila).
Promoter methylation levels and RT-qPCR analysis Quantification of single cytosine percent methylation at specific CpG sites in five genes (cyp19a1a, foxl2a, foxl2b, dmrt1, pdgfrb) was performed using targeted bisulfite-pyrosequencing method, as described in  Pierron et al. (2019). Primers were designed from the Pyromark assay design software (Table  1 ). For cyp19a1a, foxl2a and dmrt1, we used the already published method  (Pierron et al. , 2021). For foxl2b and pdgfrb, the mean methylation percent of all CpG sites was used, corresponding to three or five sites respectively. The average methylation level of one CpA was measured to confirm great conversion after bisulfite treatment (Supplementary data 1). Gene transcription was measured by RT-qPCR using the ΔΔCt method  (Livak and Schmittgen, 2001).
Statistical analysis For all measured parameters data are presented as mean values ± SE, with significance taken as P < 0. 05. Irradiation effect on cortisol and gene transcription levels were assessed using Anova and with BoxCox transformation when normality and homoscedasticity of the error terms were not verified. When assumptions were not met after transformation, a non-parametric test was used (kruskal wallis). For parametric and nonparametric test, Tukey HSD and Dunnet test with Bonferroni correction were used respectively as post hoc test. Irradiation effect on F0 methylation and sex differentiation for F0 and F1 methylation were assessed using GLM with binomial distribution. Irradiation effect on F1 methylation was assessed using a GLMM (Generalized Linear Mixed Model, glmm function in R), with binomial distribution and with spawn as a random effect as explained by  Guirandy et al.
Introduction Most top-performing approaches for predictive uncertainty estimation with Deep Neural Networks (DNNs)  [39, 2, 45, 15]  are essentially based on ensembles, in particular Deep Ensembles (DE)  [39] , which have been shown to display many strengths: stability, mode diversity, good calibration, etc. [13]. In addition, through the Bayesian lens, ensembles enable a more straightforward separation and quantification of the sources and forms of uncertainty  [16, 39, 47] , which in turn allows a better communication of the decisions to humans  [5, 37]  or to connected modules in an autonomous system  [49]. This is crucial for real-world decision making systems. Although originally introduced as simple and scalable alternative to Bayesian Neural Networks In most BNNs all weights follow Gaussian distributions and are assumed to be mutually independent: each weight is factorized by a Gaussian distribution.
Preliminaries We consider a training dataset D = {(x i , y i )} n i=1 with n samples and labels, corresponding to two random variables X ∼ P X and Y ∼ P Y. Without loss of generality we represent x i ∈ R d as a vector, and y i as a scalar label. We process the input data x i with a neural network f Θ (•) with parameters Θ, that outputs a classification or regression prediction. We view the neural network as a probabilistic model with f Θ (x i ) = P (Y = y i | X = x i , Θ). In the following, when there are no ambiguities, we discard the random variable from notations. For classification, P (y i | x i , Θ) is a categorical distribution over the set of classes over the domain of Y , typically corresponding to the crossentropy loss function, while for regression P (y i | x i , Θ) is a Gaussian distribution of real values over the domain of Y when using the squared loss function.
For simplicity we unroll our reasoning for the classification task. In supervised learning, we leverage gradient descent for learning Θ that minimizes the cross-entropy loss, which is equivalent to finding the parameters that maximize the likelihood estimation (MLE) P (D | Θ) over the training set Θ MLE = arg max Θ (xi,yi)∈D log P (y i | x i , Θ), or equivalently minimize the following loss function: L MLE (Θ) = - (xi,yi)∈D log P (y i | x i , Θ). (1) The Bayesian approach enables adding prior information on the parameters Θ, by placing a prior distribution P(Θ) upon them. This prior represents some expert knowledge about the dataset and the model. Instead of maximizing the likelihood, we can now find the maximum a posteriori (MAP) weights for P(Θ | D) ∝ P(D | Θ)P(Θ) to compute Θ MAP = arg max Θ (xi,yi)∈D log P(y i | x i , Θ) + log P(Θ), i.