dataesr/scientific-paragraphs-categorization huggingface dataset:
We consider PAM clustering at two resolution levels: First, we search for two clusters and check how well these correspond to the known cell types (i. e. , glial and neuronal); next, we search for 12 clusters and check how well these correspond to the known regions within each cell type. We evaluate the clusterings using the adjusted Rand index (ARI)  36, 37  , comparing the derived partitions with the ground truth. Interestingly, all methods typically cluster the majority or all samples correctly according to cell type, except for cqn and no normalization (Supplementary Figure  12 ). However, when checking how well the di↵erent brain regions within each cell type are recovered by clustering the normalized datasets into 12 clusters (Figure  4b ), for each selected number of PCs, GC-FQ, smooth GC-FQ, and FQ-FQ perform best, while cqn and no normalization perform worst.
The good performance of (smooth) GC-FQ and FQ-FQ was already noticeable in the PCA plots in Supplementary Figure  11 , since these are the only methods where the STR and MDT regions are clearly separated from other brain regions for the neuronal cells in 2 dimensions. Asides from clustering, researchers often focus on discovering peaks that are di↵erentially accessible between biological groups. Here, we use edgeR (or DESeq2 for DESeq2 normalization) to fit a negative binomial generalized linear model for each peak in each normalized dataset. For each peak, we test for di↵erences in average accessibility between neuronal and non-neuronal cells, across all brain regions. Mean-di↵erence plots (MD-plots)  38  show a prominent GC-content bias in the fold-changes for all GC-unaware normalization methods (Figure  4c  and Supplementary Figure  13 ).
Discussion The evaluations in this manuscript highlight the importance of accounting for GC-content e↵ects in ATAC-seq datasets. Because of the sample-specificity of GC-content e↵ects, failing to adjust for GC-content using an appropriate normalization method can bias downstream analyses such as clustering and di↵erential accessibility analysis. We have proposed GC-aware normalization procedures and benchmarked these against state-of-the-art procedures using eight public ATAC-seq datasets. While GC-aware procedures perform better than GC-unaware procedures, none uniformly outperforms all others, although smooth GC-FQ generally performs well on average. The choice of an appropriate normalization procedure is dataset-specific and exploratory data analysis is essential to guide this choice. Similar GC-content e↵ects have also been noted in DNA-seq  39  , RNA-seq  11, 21  , and ChIP-seq 9 , amongst others.
For ChIP-seq datasets, Teng and Irizarry 9 recently developed a negative binomial mixture model to correct for GC-content e↵ects in both background and binding signal regions at the peak-calling stage, by accounting for GCcontent in the abundance estimation for a particular genomic region. While their method has been evaluated using ChIP-seq data, it may also be useful for ATAC-seq data. However, none of the publicly-available and processed datasets we have looked into account for GC-content e↵ects during the peak-calling stage. While in this manuscript we have focused on correcting GC-content bias at the level of called peaks, other approaches are possible. For example, Benjamini and Speed 39 argue that it is the GC-content of the full DNA fragment (vs. only the sequenced read) that most influences read counts.
A comparison of the peak-level normalization approaches discussed here with fragment-level approaches would be an interesting avenue for further research on how to best correct for GC-content e↵ects. Our work has focused on normalization of bulk ATAC-seq datasets. While FQ-based normalization procedures were found to perform favorably in this setting, it remains to be seen whether they perform equally well on single-cell ATAC-seq (scATAC-seq) datasets. The sparsity associated with scATAC-seq data suggests that their application could be limited and alternative normalization procedures may be needed.
Datasets Philip et al. 17  study CD8 T-cell dysfunction in acutely infected and chronic tumoral tissue, over several time points. We only focus on the mouse samples and consider time and treatment as the biological variables of interest. We did not find metadata on quality control (QC) or batch variables, so we do not use any in the scone evaluation. The count matrix corresponds to 75, 689 peaks for 41 samples and was downloaded from the Gene Expression Omnibus (GEO) with accession number GSE89308. Bryois et al. 28  study the adult human prefrontal cortex brain. We remove samples that are not schizophrenic or control samples, leaving a total of 272 samples, consisting of 135 individuals with schizophrenia and 137 controls. We consider the disease status as the biological variable of interest.
In the evaluation, we use the 32 QC variables that were available in the metadata, along with the top 10 principal components derived from the patients' genotypes. The sequencing index in the metadata is used as batch variable. The count matrix corresponds to 118, 152 peaks and was obtained through personal communication with the authors. It was not relevant to correct for the width of the peaks using cqn in this dataset, since all peaks have a length of 301bp. de la Torre-Ubieta et al. 14  study human cortical neurogenesis in the germinal zone and cortical plate of the developing cerebral cortex. Samples were derived from three individual donors and each donor was handled and processed separately, so we treat each donor as a batch and the brain region as the biological variable of interest.
The count matrix corresponds of 62, 005 peaks across 19 samples and was downloaded from the GEO with accession number GSE95023. Note that the replication in this dataset is technical, i. e. , consists of samples from the same human donor. Calderon et al. 23  study a repertoire of 32 immune cell types under resting and activated conditions in humans. The metadata include three QC variables (number of peaks called, number of sequenced reads, and transcription start site enrichment for each sample), which we use in the scone evaluation. Most donors are processed and sequenced separately, and therefore each donor represents a di↵erent batch. However, some samples underwent a second round of sequencing and the accessibility counts from these two sequencing rounds were summed, so we also treat these additionally sequenced samples as a separate batch.
The biological variables of interest are cell type and treatment. The count matrix corresponds to 829, 942 peaks across 175 samples and was downloaded from the GEO with accession number GSE118189. This dataset was filtered to retain peaks with at least 2 counts per million in at least 10 samples, reducing the dataset to 203, 448 peaks. Murphy et al. 31  study photoreceptors and bipolar cells in the mouse retina. The dataset combines two experiments and we define each experiment as a batch. We consider the cell type as the biological variable of interest. The count matrix corresponds to 110, 715 peaks across 12 samples and was downloaded from the GEO with accession number GSE131625. It was not relevant to correct for the width of the peaks using cqn in this dataset, since all peaks have a length of 201bp.
Rizzardi et al. 16  study neuronal and non-neuronal cell populations in the prefrontal cortex and nucleus accumbens in humans. We consider the combination of brain region and cell type as the biological variable of interest. We define a batch variable as the combination of the donor, flow cytometry run, and sequencing date variables. The count matrix corresponds to 961, 916 peaks across 22 samples and was downloaded from the GEO with accession number GSE96614. Brain Open Chromatin Atlas (case study). Fullard et al. 29  developed a human brain atlas of neuronal and non-neuronal cells across 14 distinct brain regions from 5 human donors. We define a batch variable as the flow cytometry date. Note that while the sequencing date is nested within the flow cytometry date, there are other variables in the metadata that might also be considered to define batches, e.
g. , PCR date. A total of 49 variables corresponding to potential technical e↵ects were included as QC measures. The biological variable of interest is defined as the combination of cell type and brain region. The count matrix corresponds to 300, 444 peaks across 115 samples and was downloaded from the Brain Open Chromatin Atlas (BOCA) website, at https://bendlj01. u. hpc. mssm. edu/multireg/. Mouse Tissue Atlas (case study). Liu et al. 30  created an ATAC-seq atlas of mouse tissues, spanning a total of 20 tissues for both male and female mice. We use the "Slide lane of sequencer" variable recorded in the metadata as batch variable. Four variables (mitochondrial reads, usable reads, transcription start site (TSS) enrichment, and number of reproducible peaks, as identified using the irreproducible discovery rate (IDR)  40  ) are used as QC measures.
In practice, we implement full-quantile normalization using the following procedure 1. given a data matrix Y J⇥n for J peaks (rows) and n samples (columns), 2. sort each column to get Y S , 3. replace all elements of each row by the median (or average) for that row, 4. obtain the normalized counts Ỹ by re-arranging (i. e. , unsorting) each column. Smooth-quantile (SQ) normalization (qsmooth). Full-quantile normalization assumes that the read count distribution has a similar shape for each sample and that the observed variability in global distributional properties corresponds to technical e↵ects. However, this may not always be the case in practice. To tackle this, Hicks et al. 42  developed smooth-quantile normalization, a variant of full-quantile normalization that is able to deal with datasets where there are large global di↵erences between biological conditions of interest.
It provides a balance between (a) full-quantile normalization between samples of each condition separately, and (b) full-quantile normalization on the full dataset. This balance is struck by calculating data-driven weights for each quantile, that specify which of the two normalization options is more appropriate. The weights are estimated in a smooth way across the quantiles, by contrasting the within-condition with the between-condition variability for each quantile. If the within-condition variability is significantly smaller than the between-condition variability, then the weights will favor normalization for each condition separately. Within-and-between-sample full-quantile (FQ-FQ) normalization. The FQ-FQ method, implemented in the EDASeq package 10 , accounts for GC-content e↵ects by performing two rounds of full-quantile normalization.
First, the features of each sample are grouped into (by default, 10) GC-content bins and full-quantile normalization is performed across bins within each sample (referred to as 'within-lane normalization'). Next, the data are normalized using full-quantile normalization across all samples. Conditional-quantile normalization (cqn). The cqn method  43  starts by assuming a Poisson model for the accessibility counts Y ji. Median regression is used to model, for each sample, the log-transformed accessibility count as a smooth function of GC-content as well as peak width, focusing on peaks with high average count (above 50 by default). Note that for the datasets from Bryois et al. 28  , Murphy et al. 31  , and Rizzardi et al. 16  all peaks have the same width and hence there is no peak width normalization.
Normalization performance: scone benchmark We use the Bioconductor R package scone  33  to implement and evaluate di↵erent normalization procedures. The first step in the scone workflow is to normalize the data using all normalization methods of interest. The normalized data are then used to calculate a range of evaluation measures (see Supplementary Methods for a description and evaluation of each of these). Since some measures tend to be biased towards particular normalization methods, we rely on a subset, selected based on our evaluation as described in Supplementary Results. As part of the evaluation, the log-count principal components of the normalized data are correlated to factors of unwanted variation as well as quality control variables (if available). The factors of unwanted variation are inferred using RUVSeq 45 , based on negative control features which we here set to be the peaks within known housekeeping genes.
4. 4. 2 Di↵erential accessibility analysis performance in a synthetic null and signal scenario Our approach to evaluate the impact of normalization on DA analysis is two-fold: First, we perform synthetic null comparisons for each real dataset; second, we generate synthetic signal datasets by simulating DA peaks from each real dataset (see Methods, Datasets for a descriptions of each dataset). Synthetic null scenario. In the null scenario, for each dataset, we create a two-group mock variable so that we expect no systematic di↵erences between the groups. Specifically, for each dataset, we perform stratified random sampling, where the samples for each biological condition are randomly split into two approximately equally-sized groups (e. g. , for a biological condition with 4 samples, there are 2 samples per group, and for a biological condition with 3 samples, one group comprises 1 sample and the other 2 samples).
For each dataset, following the assignment of samples to groups, we evaluate the performance of normalization methods based on a di↵erential accessibility analysis using this mock variable. Sythetic signal scenario. Additionally, we also evaluate the performance of normalization and DA analysis methods on synthetic signal datasets created from each of the real datasest and where 10% of all peaks are DA. Each synthetic dataset comprises 12 samples (based on 6 randomly selected samples from each group in the mock variable). For each selected sample i, we calculate its accessibility fraction for each peak j as F ji = Y ji P J j=1 Y ji. (5) A random subset comprising 10% of all peaks is simulated to be di↵erentially accessible, with equal up-/downregulation between groups, via independent binary random variables S j , equal to either 1 or 1, each with 1/2 probability.
The S j 's define the group g 2 {1, 2} of samples for which the accessibility fractions will be altered as follows ji = ( F ji exp I(Sj = 1)Zj , if g(i) = 1 F ji exp I(Sj =1)Zj , if g(i) = 2 , (6) where g(i) 2 {1, 2} denotes the group to which sample i belongs and Z j are independent Gaussian random variables with mean 0. 8 and standard deviation 0. 1. That is, the log-fold-change in accessibility between groups 2 and 1 is S j Z j. The choice of the mean and standard deviation for the log-fold-changes corresponds to fold-changes being on average 2. 25, with a minimum of ⇠ 1. 5 and a max of ⇠ 3. 5, a reasonable scenario. Sequencing depths N ⇤ i are simulated from a uniform distribution N ⇤ i ⇠ U (N min , N max ), where N min = min i P j Y ji and N max = max i P j Y ji denote, respectively, the minimum and maximum library sizes across all samples in the dataset.
Given ji and N ⇤ i , accessibility counts Y ⇤ ji are then simulated using a Multinomial distribution Y ⇤ ji | ji , N ⇤ i ⇠ Mult( ji , N ⇤ i ). (7) For each combination of sample size and DA signal strength, we evaluate 14 simulated datasets (a greater number of simulations was not feasible due to local memory limitations). Di↵erential accessibility analysis is performed using edgeR for all normalization methods, except for DESeq2 normalization where we rely on the native DESeq2 pipeline. For RUVg normalization, we use peaks overlapping with housekeeping genes as negative control features, as in the scone evaluation, excluding features that were simulated to contain signal between the groups. For each simulated dataset, we calculate the true positive rate (TPR) and false discovery rate (FDR), defined as T P R = T P T P + F N (8) F DR = F P T P + F P , where F N, F P , and T P denote, respectively, the numbers of false negatives, false positives, and true positives.
Mouse Tissue Atlas The raw count matrix from Liu et al. 30  is obtained as described in the Datasets section and is normalized using each of the twelve evaluated normalization procedures, with no peaks filtered out. Euclidean distances between samples are calculated on the log-normalized counts, adding an o↵set of 1 to avoid taking the log of zero. Hierarchical clustering trees are derived using complete linkage. Di↵erential accessibility analysis is performed using edgeR for all normalization methods, except for DESeq2 normalization where we rely on the native DESeq2 pipeline. Normalized counts are used directly as input for FQ, FQ-FQ, and (smooth) GC-FQ normalization, while normalization o↵sets are used for TMM and cqn. For RUVSeq normalization, we incorporate 4 inferred factors of unwanted variation in the design matrix.
Brain Open Chromatin Atlas The raw count matrix from Fullard et al. 29  is obtained as described in the Datasets section. We do not filter out any peaks. PCA and hierarchical trees are based on the log-normalized counts, adding an o↵set of 1 to avoid taking the log of zero. Hierarchical trees use the Euclidean distance between samples and are constructed using complete linkage. DA analysis is performed as described in the Mouse Tissue Atlas case study. The contrast matrix is defined for comparing the average expression of neuronal vs. non-neuronal samples.
Acknowledgements The authors thank  Michael Love  for his input on scaling normalization methods edgeR and  DESeq2 , and  Ameek Bindra , who contributed to this project through the  Undergraduate Research Apprenticeship Program (URAP)  of the  Department of Statistics at the University of California, Berkeley. KVdB is a postdoctoral fellow of the  Belgian American Educational Foundation (BAEF)  and is supported by the  Research Foundation Flanders (FWO)  grants  1246220N  and  G062219N. H-J C was supported by a postdoctoral research grant from the  UC Berkeley Siebel Stem Cell Center. SD and JN are supported by the  National Institutes of Health , grant  R01 DC007235.
17 17 Figure 2: GC-aware normalization methods cqn, FQ-FQ, and GC-FQ are successful in eliminating GC-content e↵ects on the di↵erential accessibility fold-change estimates. (a) Accessibility distributions for three replicates from the dataset of Philip et al. 17. The peaks are grouped into 10 equally-sized bins according to their GC-content (rows) and the accessibility distribution (kernel density estimate) is plotted for each bin. The distributional shapes are more comparable across samples for a particular GC-content bin, than they are across GC-content bins for a particular replicate. (b-d) There is no visible GC-content e↵ect on fold-changes estimated using edgeR following normalization with GC-aware methods cqn, EDASeq, and GC-FQ, in the mock comparisons for the dataset from Calderon et al.
e. , first-order kinetics) leading to an exponential increase in the biomass of individual cells over time  (Collins and Richmond, 1962; Koch and Schaechter, 1962; Koch, 1966). While the linear model of substrate assimilation is most commonly used to calculate assimilation rates in cells from nanoSIMS data (e. g. ,  Foster et al. , 2013; Klawonn et al. , 2016; Stryhanyuk et al. , 2018; Mills et al. , 2020) , the exponential model is used rather rarely (e. g. ,  Martínez-Pérez et al. , 2016; Berthelot et al. , 2019; Geerlings et al. , 2020 Geerlings et al. , , 2021;; Polerecky et al. , 2021). The impact of the substrate assimilation model (linear vs. exponential) on the calculated assimilation rate has previously not been studied in detail. Assimilation rates determined at single-cell resolution often reveal intercellular heterogeneity in populations.
In the latter case, θ(m) in Eq. 18 and Eq. 16 become equivalent as required (Supplementary Figure  4C ). Note that if γ is large (i. e. , C C max ) and each cell in the population grows at the same cell-specific rate r, Eq. 18 implies that 95% of cells divide within a time interval t 95 = 4 • C/r (Supplementary Figure  3 ). Thus, for a population with cell cycles that are highly synchronized, the parameter γ is equal to the ratio between the population doubling time, τ = C max /(2 • r) (Table  1 ), and the interval during which 95% of cells undergo division (i. e. , γ = τ/ t 95 ; Supplementary Figure  3 ). We performed Monte-Carlo simulations in Matlab to evaluate ζ from θ. Specifically, a cell j with a random initial C content C ij is selected based on the PDF θ. Then, for a given value of r and t, the isotopic composition of daughter cells originating from cell j is calculated as x E Sj = Z t, C ij , C max , r.
Calculating Rates of Cellular C Assimilation Using the findings presented above, we suggest a three-step procedure for calculating the cellular rates of C assimilation based on 13 C atom fractions measured by nanoSIMS (Figure  5 ). We emphasize the underlying assumptions of each step to provide the foundations of the formulas. Numerical implementation of the calculation steps, including the calculation of the best estimate and uncertainty of the rates, is provided as a Matlab script that can be interfaced with experimental data organized in a spreadsheet. The script is available at https://github. com/ lpolerecky/LARS. Previous studies documented that isotopic composition of cells can be significantly affected by label loss or dilution associated with common techniques of sample preparation for nanoSIMS analysis including chemical fixation, dehydration and resin embedding  (Musat et al.
Application to Hypothetical Data To illustrate the utility of the procedure described above, we calculated C assimilation rates based on hypothetical data from FIGURE 5 | Flowchart for quantifying C assimilation rates in single cells. First, cells are incubated with a 13 C-labeled substrate and their 13 C atom fractions, x, are measured by nanoSIMS. Second, the C content of the cells is constrained. Rate calculation then proceeds in three steps: (1) convert x to x E S to account for the isotopic composition of the substrate, (2) calculate the carbon-specific C assimilation rate, k (this rate is independent of the C content of the cell), and (3) calculate the cell-specific C assimilation rate, r. Depending on how well the C content of the measured cell can be constrained, the third step proceeds either by using the value of k from step 2 (Approach A and B) or using the value of x E S from step 1 and the inverse of the zig-zag function, Z -1 (Approach C).
The calculated r then represents the average rate across the cell cycle (Approach A), the instantaneous rate at the end of the SIP experiment (Approach B), or the average rate during the SIP experiment (Approach C). Steps 2 and 3 consider that the measured cell may be a product of a cell that divided during the incubation. The calculation approaches are implemented in a Matlab script available via GitHub (https://github. com/lpolerecky/LARS). cells incubated with a 13 C-labeled C source (x S,eff = 0. 1) for 2 hours. Measured cells had low, intermediate, and high levels of 13 C labeling (x 1 = 0. 02, x 2 = 0. 033, x 3 = 0. 09). Using Steps 1 and 2, these values yield carbon-specific assimilation rates of k 1 = 0. 053 h -1 , k 2 = 0. 14 h -1 , and k 3 = 1. 1 h -1 , which correspond to doubling times of about τ 1 = 13 h, τ 2 = 4.
ACKNOWLEDGMENTS This work was inspired by discussions and observations that emerged from the  Group for Aquatic Photosynthesis (GAP)  workshop held in Třeboň, Czech Republic, in August 2017. We are particularly grateful to  Ondřej Prášil  for organizing the workshop. We also thank  Celine Loussert-Fonta ,  Marcel Kuypers ,  Peter K. Weber , and  Martin Von Bergen  whose constructive comments helped improve the manuscript.