dataesr/scientific-paragraphs-categorization huggingface dataset:
This paper develops TWO-SIGMA-G, a set-level framework for DE testing in scRNA-seq data using the competitive null hypothesis. To test for DE at the gene-level, TWO-SIGMA-G uses our recently developed TWO-SIGMA  [38] , providing a flexible mixed-effects zero-inflated negative binomial regression model for good fit to the data at the gene-level. The use of a regression-based framework allows complex experiments including multiple covariates to be analyzed and provides many choices for gene-level statistics depending on the priority of the analysis. To avoid the inflated type-I error often caused by commonly existing positive IGC in a biological gene regulation pathway, IGC is estimated using an innovative residual-based approach and explicitly adjusted for at the set-level. We demonstrate TWO-SIGMA-G outperforms existing competitive gene set testing methods using extensive simulation scenarios.
We further use individual-level calculations to help mitigate the impacts of the large individual heterogeneity often seen in scRNA-seq datasets. In simulations, we found that this IGC estimate preserves type-I error in a conservative manner while still producing improved power in a variety of realistic scenarios. The estimate of the IGC is virtually free computationally in that the model is not refit via permutation or bootstrapping. Figure  2  shows the average set-level IGC estimates from TWO-SIGMA-G for each of the two comparisons in the Alzheimer's dataset described more in the real data analysis section. A non-negligible correlation exists in the residual space for both comparisons, with half of the sets having a correlation larger than 0. 02. For each comparison, over 98% of the estimated pairwise correlations are positive.
TWO-SIGMA-G for Set-Level Testing We extend our TWO-SIGMA method  [38]  to competitive gene set testing via TWO-SIGMA-Geneset (TWO-SIGMA-G), on overview of which is shown in Figure  1. Briefly, TWO-SIGMA uses a zero-inflated negative binomial regression model to test for DE at the gene level in scRNA-seq data. It is flexible and can be customized in several different ways. First, the zero-inflation component can be removed from the model entirely (as was done in the real data analysis section), leaving a standard negative binomial regression model. Second, the model can additionally include random effect terms to account for cell-cell correlation within the same sample and limit type-I error inflation in gene-level DE inference. Finally, many gene-level statistics measuring the evidence of DE can be used for set-level testing.
For example, the likelihood ratio or z statistic corresponding to a treatment effect are commonly used gene-level statistics. We discuss uses for other, more complex, gene-level statistics based on custom contrasts of regression parameters in the real data analysis section. TWO-SIGMA-G employs the Wilcoxon rank-sum test to compare the statistics of genes in the test set to the statistics of genes in the reference set, and therefore uses the sum of the ranks in the test set as the set-level summary statistic. In using the ranks, TWO-SIGMA-G provides robustness against the influence of very large gene-level statistics. Traditionally, the Wilcoxon rank-sum test assumes that observations within a group are independent. However, as mentioned, IGC is expected given the construction of gene sets as harmonious biological pathways, and can inflate type-I error if ignored  [43].
Using this modified variance, and the known mean of rank-sum statistics under the null, set-level p-values are computed analytically using a standard normal approximation  [43]. The reference set used in TWO-SIGMA-G can be chosen in one of two ways: either using a random sample of other genes of size m 1 or as the collection of all genes not in the test set under consideration. In addition to producing set-level significance, TWO-SIGMA-G also identifies the directionality of sets as up or down-regulated. Whether or not a zero-inflation component is included in gene-level models, directionality is produced by averaging gene-level log fold-change estimates in the test set to produce a set-level effect size and taking the sign of the result. These effect sizes are demonstrated further in the real data analysis section.
Larger gene sets are often more significant, even if additional genes represent noise, and other gene sets not being tested can influence results in counterintuitive ways  [6, 36]. Although GSEA is not strictly a competitive test, we still include it for comparison given its popularity for performing gene set testing in both bulk and scRNA-seq data  [31, 7]. GSEA was accessed via the GSEA R package (version 1. 2, available at https://github. com/GSEA-MSigDB/GSEA R). Other than the choice of 200 permutations, default options from the GSEA R package were used. fGSEA was developed as an alternative to GSEA which relies on gene-level summary statistics to improve the computational performance of GSEA. Rather than estimating the null distribution of the enrichment statistic for each set separately, as is done in GSEA, fGSEA employs a custom procedure to simultaneously estimate the null distribution for all sets.
was accessed via fgsea R package (version 1. 10. 1). Ten thousand permutations were used, and otherwise default options from the fgsea R package were used. PAGE compares the log fold-change values in the test set to those from the complement set of genes using a one-sample z-test, where the sample mean and variance are estimated using all genes. PAGE assumes genes are independent and does not adjust for IGC, but is computationally efficient as compared to GSEA. It may fail to preserve type-I error in some cases, however. PAGE was accessed via the PGSEA function of the PGSEA R package (version 1. 60). CAMERA was developed as a competitive gene set test for microarray or RNA-seq data. Gene-level statistics are first constructed using a linear model, meaning that CAMERA can accommodate complex experimental designs beyond a two-group comparison.
Set-level p-values are then computed using modifications of the t-test or Wilcoxon rank-sum test that allow for a common pairwise correlation in the test set. Rather than using the raw data to estimate the IGC, CAMERA uses the residuals from the linear model. As discussed above, use of the residuals means that the variation in gene expression explained by the covariates is removed, giving the most reliable estimate of the correlation between the gene-level statistics in the test set. By avoiding permutation, and unlike some early approaches, CAMERA provides a statistically valid, computationally efficient test of a precisely defined and fully specified null hypothesis  [13]. The hypothesis corresponds to a test that the average absolute value of each coefficient in the test set is larger in magnitude as compared to the reference set.
Two scRNA-seq data sets We performed gene set testing on two different datasets to illustrate the usefulness of TWO-SIGMA-G. Gene sets were taken from the Molecular Signatures Database (mSigDB)  [34, 20]  version 7, c2 collection, accessed via the msigdf R package (https://github. com/ToledoEM/msigdf). All gene sets with at least two genes present after filtering were analyzed in each of two datasets described below. We used Fisher's method to combine p-values over all cell types and create a consensus ranking of pathways (see, for example, Figure  6 (A) ). Both datasets were from 10X UMI-based single cell RNA-seq platforms. As with other UMI-based scRNA-seq count data, we found that this data was not consistent with zero-inflation  [35] , and thus we fit the TWO-SIGMA model without the zero-inflation component at the gene-level.
HIV Dataset Our first dataset consists of 11,630 single cells collected from humanized mice  [5]. We kept only genes with a proportion of zeros no higher than the mean percentage, leaving the most relevant and highly expressed 3,549 genes. A total of nine cell types are present in the data: 4,249 natural killer (NK) cells, 2,085 erythroid cells, 1,421 innate lymphoid cells (ILC), 1,205 B cells, 1,088 myeloid dendritic cells (mDC), 821 plasmacytoid dendritic cells (pDC), 555 progenitor cells, 126 macrophage cells, and 80 mast cells. The read counts are then treated as the outcome of interest. Because the primary interest is in comparisons between HIV and mock cells within cell-type, we categorize cells into one of 2*9 = 18 mutually exclusive groups. An ANCOVA model additionally adjusting for CDR was fit as a way to test for cell-type specific differences in expression levels comparing HIV to mock.
Alzheimer's Dataset Our second dataset consists of 70,634 single cells from human donors  [25]. A total of 48 individual donors are present, categorized into three pathology groups: 24 individuals are control patients free of a diagnosis of AD, 12 were diagnosed with early stage AD, and 12 were diagnosed with late stage AD. Cells from the six most common cell-types were analyzed: 34,976 excitatory neurons (Ex), 18,235 oligodendrocytes (Oli), 9,196 inhibitory neurons (In), 3,392 astrocytes (Ast), 2,627 oligodendrocyte progenitor cells (Opc) and 1,920 microglia cells. We did not remove cells beyond what was done in the original manuscript because extensive quality control was performed on the dataset we used. We chose to filter the original 17,926 genes to the 6,048 most highly expressed genes by removing genes unexpressed in at least 90% of cells.
Overview of TWO-SIGMA-G and Summary of Simulation Design TWO-SIGMA-G is a gene set testing method based on differential expression specifically designed for scRNAseq data. Because proper gene-level DE testing is key for a gene set test, TWO-SIGMA-G employs a mixedeffects zero-inflated negative binomial (ZINB) regression model we previously developed  [38]  to test for DE at the gene-level. This model can simultaneously capture the distribution of scRNA-seq datasets, account for within-sample correlation, analyze complex experiments by controlling for additional covariates, and provide options for gene-level test statistics. If zero-inflation exists at the gene-level, a gene-level model including both logistic and negative binomial components is fit. Otherwise, a one-component negative binomial model is used.
The standard output of TWO-SIGMA-G includes gene-level DE summary statistics and associated p-values, set-level p-values and IGC estimates, and set-level effect sizes to characterize pathways as 'Up' or 'Down' regulated. When multiple cell types are included as a contrast, all information except IGC estimates is cell-type-specific, allowing users to investigate cellular heterogeneity. We applied TWO-SIGMA-G to two scRNA-seq datasets, one experimental and one observational. TWO-SIGMA-G is freely available in the twosigma R package at https://github. com/edvanburen/twosigma. Full details of our proposed TWO-SIGMA-G is described in the Methods section. We also performed extensive simulations to compare TWO-SIGMA-G to six other methods designed for gene set testing, the full details of which are given in Supplementary Section S1.
Analysis of HIV data reveals biologically expected findings We analyzed a dataset of 11,630 single-cells collected from of 4 humanized donor mice, two of which were infected with HIV and two were given a mock treatment  [5]  (also described in the Methods section). A total of 3,549 genes and 4,772 gene sets from the Molecular Signatures Database (MSigDB)  [34, 20]  were analyzed. The number of differentially expressed genes and sets in all cell types comparing HIV to a mock treatment is shown in Table  1. An increase in the number of significantly DE genes does not always correspond to more significantly DE gene sets. For example, erythroid cells have the second largest number of DE genes, but rank eighth in terms of the number of DE gene sets. This observation is expected in a competitive gene set test like TWO-SIGMA-G, because competitive tests focus on the relative signal of gene sets as compared to a background reference set of genes.
Analysis of Alzheimer's data reveals cell-type-specific heterogeneity in setlevel expression We chose to use a second scRNA-seq dataset to demonstrate how TWO-SIGMA-G can handle a more complex application. Specifically, we use the scRNA-seq data of  [25]  (see the Methods section for more details) to analyze changes in gene expression as Alzheimer's Disease (AD) progresses. The data provides gene expression across three distinct pathology groups: control (AD free), early stage AD progression, and late stage AD progression. We focus on two relevant comparisons: late vs. early stage AD, and early stage AD vs. control. A total of 6,048 genes and 5,074 gene sets from the Molecular Signatures Database (mSigDB)  [34, 20]  were analyzed. We used TWO-SIGMA-G to compare early stage AD patients to control in each of the cell types and provide cell-type-specific results (Figure  8 ).
For a general application, however, TWO-SIGMA-G can help researchers to investigate cell-type-specific heterogeneity when cellular functions are not understood. The ability to test complex gene-level hypotheses as contrasts of regression parameters increases the diversity of cell-type-specific hypotheses that can be explored. Finally, we also compared the performance of TWO-SIGMA-G to MAST and CAMERA for both comparisons using the Alzheimer's data (Supplementary Figures  S24 to S27 ) in Ex cells (the most common cell type). Once again, we focused on comparisons to MAST and CAMERA because both can use regression models to control for additional covariates in this observational dataset and both estimate and adjust for IGC, which we found to be distinctly positive in this dataset (Figure  2 ).
However, we suggest that inference in real data analyses is likely not influenced greatly at the set-level by the presence or absence of gene-level random effects (Supplementary Figure  S29 ). TWO-SIGMA-G is implemented in the twosigma R package (https://github. com/edvanburen/twosigma), which is computationally efficient and allows for parallelization. To benchmark computational performance, we ran a modified version of our HIV data analysis, testing for a treatment effect of HIV pooled over all cell 13 types. This modification to a one degree of freedom hypothesis allows us to test identical hypotheses in TWO-SIGMA-G, MAST, and CAMERA to provide a fairer comparison of computation. Using three computing cores on a MacBook Pro laptop, the methods had the following respective runtimes: 33. 2 minutes for TWO-SIGMA-G, 33.
5 minutes for MAST (25 bootstrap replications), and 5 seconds for CAMERA. TWO-SIGMA-G shows slightly improved yet nearly identical computational performance to MAST in the presence of the other advantages for performing gene set testing in scRNA-seq data described throughout this paper. Unlike bulk RNA-seq data, many genes are often uncaptured or fail to survive filtering in scRNA-seq data. In gene set analysis, we must therefore assume that a gene set can be represented by the genes that exist in the dataset. In both the HIV and Alzheimer's datasets, we typically have around 40% representation regardless of set size after gene filtering (Supplementary Figure  S30 ). Given the biologically meaningful and interpretable results we presented, the absence of these genes does not seem to threaten the ability of scRNA-seq gene set analyses to contribute new biological insights.
Data Availability Both datasets analyzed in this manuscript are publicly available. The HIV dataset is available at the Gene Expression Omnibus under accession GSE148796. The Alzheimer's dataset is available upon completion of a data usage agreement at the Rush Alzheimer's Disease Center (RADC) Research Resource Sharing Hub (https://www. radc. rush. edu/docs/omics. htm) under "snRNA-seq PFC. " TWO-SIGMA-G is implemented in the function twosigmag in the twosigma R package, which is freely available on GitHub at https://github. com/edvanburen/twosigma.
Acknowledgements None at this time. Figure  3 : Type-I error performance for CAMERA, GSEA, MAST, and TWO-SIGMA-G using a reference set size of 30 genes. GSEA, CAMERA, and MAST were chosen for comparison here because all utilize the raw expression data, and the latter two use a regression modeling framework that can explicitly analyze complex experimental designs similar to the data analyzed in the real data analysis section. See Figure  4  for additional comparisons to summary-statistic based gene set testing methods fGSEA, iDEA, and PAGE. Each subfigure varies the existence of IGC between genes in the test set and the presence of gene-level random effect terms in the gene-level model (CAMERA and GSEA never include gene-level random effect terms). Within each subfigure, both unadjusted and adjusted set-level p-values are plotted, where available.
Scenarios along the x-axis of each subfigure vary the percentage of differentially expressed genes (with the same effect size) in the test and reference sets. For example, "T80,R50" corresponds to the configuration under the alternative hypothesis in which 80% of test set genes are DE and 50% of reference set genes are DE. Significance is determined using a p-value threshold of 0. 05, and 30 genes were included in the test and reference sets (Supplementary Figure  S7  presents results using a reference set of 100 genes). Note that GSEA did not output p-values for the "R0" scenarios. See the Methods section of the main text and Supplementary Section S1 for more details regarding the simulation procedure. Table  1 : Shows the number of differentially expressed genes (using TWO-SIGMA) and gene sets (using TWO-SIGMA-G) after FDR-adjustment for the HIV dataset.
-adjusted DE Genes in HIV Dataset Among 4 Most Common Cell Types Overlap Between FDR 477 63 (27. 6%) (3. 65%) 22 (1. 27%) 17 (0. 98%) count 41 27 (2. 37%) (1. 56%) 67 400 102 (3. 88%) 200 300 69 (5. 9%) 35 100 (3. 99%) 92 46 (2. 03%) (5. 32%) (2. 66%) 55 (3. 18%).
Thus, potential GABAergic currents appear outward at potentials more depolarized than ECl whereas cationic currents are inward. All electrophysiological experiments were performed with a MultiClamp700B (Axon Instruments) amplifier connected to an acquisition board (Digidata 1440; Axon Instruments) attached to a computer running pCLAMP software (Axon Instruments). Axonal fibers from the PFC expressing ChR2 were stimulated in the VLPO by pulses of 470 nm light of 5ms every 5s for 1 to 5 min using a diode laser (pE-2, Coolled). The currents evoked by the photo-stimulations were recorded in voltage-clamp configuration and characterized using the application of the following antagonists: kynurenate (1mM; Sigma); tetrodotoxin (TTX, 1µM; Sigma); 4 aminopyridine (4-AP, 100µM; Sigma). The currents recorded at the same imposed potential were averaged for each cell to perform the quantification.
In vivo Electrophysiological recordings Signals from all electrodes were recorded using an Intan Technologies amplifier chip (RHD2216, sampling rate 20 KHz). Local Fields Potentials (LFPs) were sampled and stored at 1250 Hz. Analyses were performed with custom-made Matlab programs, based on generic code that can be downloaded at http://www. battaglia. nl/computing/ and http://fmatoolbox. sourceforge. net/. Local field potentials were recorded using tungsten wires with PFA isolation (0. 002" bare). Tetrodes were plated to reduce impedances to around 100 kΩ using gold solution (Neuralynx).
Sleep scoring Once the data had been pre-processed, an automatic analysis of vigilance states developed in the team  (Bagur et al. , 2018)  was carried out using Matlab software. This analysis is based on the signals of the local field potentials (LFP) of the olfactory bulb (OB) and the hippocampus (HPC). Briefly, the gamma power in the OB is low during both REM and NREM sleep, allowing us to distinguish wakefulness from sleep states. NREM sleep is distinguished from REM sleep based on the theta (5-10Hz) / delta (2-4) power ratio in the HPC. Indeed, during NREM sleep there is a strong presence of slow waves (high power in the delta band), while theta oscillations are highly prominent during REM sleep.
Data analysis and statistics. Data analysis and statistics were performed in the Matlab environment (MathWorks) and R (version 4. 2. 2). All values are expressed as means ± s. e. m and statistical significance on all figures uses the following convention of p-values: * p <0. 05, ** p < 0. 01, *** p <0. 001. A normality test was performed on each dataset using the Shapiro-Wilk test. The parametric tests (paired or unpaired t-tests, one-way or two-way ANOVAs) were used if the datasets were normally distributed or if no non-parametric equivalent was available. For all optogenetic experiments, paired or unpaired t-tests were performed accordingly. Paired t-tests were also employed for the effect of the chemogenetic inhibition of the PFC-VLPO pathway under basal conditions (without SDS procedure).
Results We analyzed data from concurrent whole-head MEG and pupil recordings in 81 human participants, collected at three MEG laboratories (Glasgow, UK; Hamburg & Münster, Germany). All recordings were carried out under constant illumination and in sound-attenuated environments, thus removing fluctuations in pupil size due to changes in sensory input. Participants were instructed to keep their eyes open, while fixating and resting otherwise. The length of individual recordings varied between 5 and 10 min (see Materials and methods for further details). For each of the three MEG laboratories, we obtained time courses of MEG activity and of pupil diameter fluctuations. A first control analysis of the recordings (Figure  1 -figure supplements 1 and 2) established that they were highly comparable between laboratories, despite the differences in locations, setups, and participants.
For all three MEG laboratories, standardized-MI spectra showed similar patterns of widespread pupil-brain associations across frequencies with prominent peaks in the 8-32 Hz range (Figure  2A ). MI scalp distributions within this frequency range further underpinned the commonalities between recording sites despite different MEG systems and setups (Figure  2B ). Significant relationships with pupil size were also evident for the low-frequency (4-8 Hz) and high-frequency (64-128 Hz) components of MEG power (Figure  2A ) in all three datasets. A possible confound in identifying relationships between pupil and cortical dynamics are eye movements: Saccades change pupil size  (Mathôt, 2018; Mathôt et al. , 2015)  and oculomotor behavior is functionally linked to cortical alpha oscillations  (Popov et al.
Fluctuations in pupil, and pupil derivative, may primarily reflect cholinergic and noradrenergic activity, respectively  (Reimer et al. , 2016). Consistent for all three datasets, the cross-correlation analysis for pupil derivative revealed peaks closer to zero lag for low and high frequencies (2-4 Hz negative peak at 277 ms, 64-128 Hz positive peak at -90 ms; Figure  3 -figure supplement 1). In contrast, correlations with activity in the 8-16 Hz range peaked somewhat later (700 ms; negative peak; Figure  3 -figure supplement 1). The shorter lags for the pupil derivative are consistent with previous findings in rodents  (Reimer et al. , 2016). It should be noted that, due to the nature of the employed wavelet analysis, the temporal resolution of the spectral estimates varies across frequency bands, which poses a potential challenge when comparing peak latencies across frequencies.