Avoiding false discoveries in single-cell RNA-seq by revisiting the first Alzheimer’s disease dataset

  1. Alan E Murphy  Is a corresponding author
  2. Nurun Fancy
  3. Nathan Skene  Is a corresponding author
  1. UK Dementia Research Institute at Imperial College London, United Kingdom
  2. Department of Brain Sciences, Imperial College London, United Kingdom

Abstract

Mathys et al. conducted the first single-nucleus RNA-seq (snRNA-seq) study of Alzheimer’s disease (AD) (Mathys et al., 2019). With bulk RNA-seq, changes in gene expression across cell types can be lost, potentially masking the differentially expressed genes (DEGs) across different cell types. Through the use of single-cell techniques, the authors benefitted from increased resolution with the potential to uncover cell type-specific DEGs in AD for the first time. However, there were limitations in both their data processing and quality control and their differential expression analysis. Here, we correct these issues and use best-practice approaches to snRNA-seq differential expression, resulting in 549 times fewer DEGs at a false discovery rate of 0.05. Thus, this study highlights the impact of quality control and differential analysis methods on the discovery of disease-associated genes and aims to refocus the AD research field away from spuriously identified genes.

eLife assessment

This paper reports a useful finding on the impact of choices of quality control and differential analysis methods on the discovery of disease-associated gene expression signatures. The study provides a solid comparison of the data process by re-analysis of a large-scale snRNA-seq dataset for Alzheimer's disease. This paper would be of interest to the community as to rigorous analyses for large-scale single-cell datasets.

https://doi.org/10.7554/eLife.90214.3.sa0

Introduction

Mathys et al., 2019 undertook the first single-nucleus RNA-seq (snRNA-seq) study of Alzheimer’s disease (AD). The authors profiled the transcriptomes of approximately 80,000 cells from the prefrontal cortex, collected from 48 individuals – 24 of which presented with varying degrees of AD pathology. (Mathys et al., 2019) data processing and quality control (QC) strategy for their snRNA-seq data was state of the art at this time. Furthermore, the authors took extra measures in an attempt to ensure the reliability of their results. Here, we reanalyse this data as not a criticism of the study, but as an endeavour to raise awareness and provide recommendations for rigorous analysis of single-cell and single-nucleus RNA-seq data (sc/snRNA-seq) for future studies. Most importantly, we aim to ensure that the AD research field does not focus on spuriously identified genes.

Results and discussion

Our questions of Mathys et al., 2019 focus around their data processing and their differential expression (DE) analysis (Figure 1). Firstly, in relation to their processing approach, the authors discussed the high percentages of mitochondrial reads and low number of reads per cell present in their data. This is indicative of low cell quality (Ilicic et al., 2016); however, we believe the authors’ QC approach may not capture all of these low-quality cells. Moreover, the authors did not integrate the cells from different individuals to account for batch effects. As the field has matured since the authors’ work was published, dataset integration has become a common step in sc-RNA-seq protocols and is recommended by some to remove confounding sources of variation (Heumos et al., 2023; Amezquita et al., 2020; Tran et al., 2020). To gain advantage of these recent approaches, we used scFlow (Khozoie et al., 2021) to reprocess the authors’ data. This pipeline included the removal of empty droplets, nuclei with low read counts and doublets, followed by embedding and integration of cells from separate samples and cell typing. scFlow combines best-practice approaches for processing sc/snRNA-seq datasets; see ‘Materials and methods’ for a detailed explanation of these steps. Reprocessing resulted in 50,831 cells passing QC, approximately 20,000 less than the authors’ postprocessing set with differing cell-type proportions (Figures 2 and 3).

Pseudobulk differential expression results in far less dubious disease-related genes.

(a, b) The log2 fold change and -log10 false discovery rate (FDR) of the differentially expressed genes (DEGs) from the authors’ original work (Mathys et al.) and our reanalysis (Our analysis). In (b), we have marked an FDR of 5 × 10–7, dashed grey line, to highlight how small the p-values from Mathys et al.’s analysis are. For (a, b), n is based on the number of DEGs: 26 for our analysis and 23,923 for Mathys et al. (c–g) show the Pearson correlation between the cell counts after quality control (QC) and the number of DEGs identified - n is the 6 cell types tested. For (f, g) analysis, the samples have been randomly mixed between case and control patients - n = 100 random permutations. The different cell types are astrocytes (Astro), excitatory neurons (Exc), inhibitory neurons (Inh), microglia (Micro), oligodendrocytes (Oligo), and oligodendrocyte precursor cells (OPC).

The nuclei that were removed from our quality control approach as their proportion of mitochondrial reads were ≥10%, but kept in the authors’.

(a) shows the proportion of mitochondrial reads across the different cell types. (b) gives the number of removed nuclei which were kept by the authors. The different cell types are astrocytes (Ast), excitatory neurons (Ex), inhibitory neurons (In), microglia (Mic), oligodendrocytes (Oli), and oligodendrocyte precursor cells (Opc).

The proportion of cells left after quality control (QC) from the authors’ processing approach (Mathys et al.) and our standardised pipeline approach – scFlow (Our analysis).

With regards to data quality, it is worth noting that over 99% of nuclei had less than 200 genes expressed (Table 1). While this QC step was not unique to our reprocessing, the authors made the same exclusion in their analysis (Mathys et al., 2019), it highlights the relatively low quality of the data which may be attributable to the early stage of snRNA-seq technology of the time. For example, Brase et al.’s recent study of snRNA-seq of autosomal-dominant AD (Brase et al., 2023) used a more stringent cut-off for the minimum number of expressed genes and still kept 27% (122 times more) of the assayed cells after all QC steps. Moreover, the authors discussed the high percentages of mitochondrial reads in their data. The differences in approaches to filtering based on the proportion mitochondrial reads accounts for the notable discrepancy in the number of nuclei after QC between our approach and the authors’. Our approach used a 10% cut-off for the proportion of mitochondrial reads in a nuclei, as set out in Amezquita et al.’s best-practice guidelines (Amezquita et al., 2020), which is less stringent than Seurat’s guidelines (5%) (Hao et al., 2021) or that from Heumos et al., 2023 (8% from a median absolute deviations [MAD]-based cut-off selection). Conversely, the authors filtered out high mitochondrial read nuclei based on clusters from their t-SNE projection of the data (Mathys et al., 2019). Even at our lenient cut-off, over 16,000 nuclei that were removed in our QC pipeline were kept by the authors’ Figure 2, explaining the discrepancy in the number of nuclei after QC. Based on Figure 2, it is clear that the authors’ approach was ineffective at removing nuclei with high proportions of mitochondrial reads which is indicative of cell death (Heumos et al., 2023; Ilicic et al., 2016) – both excitatory and inhibitory nuclei with higher than 75% reads from the mitochondria were kept in the final processed dataset by the authors. We have made the data from our alternative processing approach publicly available (through Synapse: https://doi.org/10.7303/syn51758062.1) so that researchers can utilise this resource free of low-quality nuclei.

Table 1
Overview of the aggregated number of cells across samples removed at each step of the quality control (QC) as part of scFlow.

Note that cells can fail QC for more than one check, so only the total failed and total passed rows will sum to 100%.

QC stepsTotal cellsPercentage
Pre-QC35,389,440
Total failed35,337,87499.85
 Minimum library size (n < 200)35,307,28199.77
 Maximum library size47420.01
 Minimum expressed genes (n < 200)35,312,43499.78
 Maximum library size/expressed genes (MAD> 4)21490.01
 Proportion of mitochondrial genes (≥ 0.1)1,097,7383.10
 Multiplets (pK = 0.0054)5810.00
Total passed51,5660.15
  1. MAD, median absolute deviation.

Our second question of Mathys et al., 2019 is their DE approach. The authors conducted a DE analysis between the controls and the patients with AD pathology, concentrating on six neuronal and glial cell types; excitatory neurons, inhibitory neurons, astrocytes, microglia, oligodendrocytes, and oligodendrocyte precursor cells, derived from the Allen Brain Atlas (Tasic et al., 2018). They performed downstream analysis on their identified differentially expressed genes (DEGs) and investigated some of the most compelling genes in more detail. Therefore, all findings put forward by their paper were based upon the validity of their DE approach. However, for this approach, the authors conducted a two-part, cell- and patient-level analysis. The cell-level analysis took each cell as an independent replicate, and the results of which were compared for consistency in directionality and rank of their DEGs against their patient-level analysis, a Poisson mixed model. The authors identified 1031 DEGs using this combinatorial approach – DEGs requiring a false discovery rate (FDR) < 0.01 in the cell-level and an FDR < 0.05 in the patient-level analysis. It is important to note that this cell-level DE approach, also known as pseudoreplication, overestimates the confidence in DEGs due to the statistical dependence between cells from the same patient not being considered (Murphy and Skene, 2022; Squair et al., 2021; Zimmerman et al., 2021; Lazic, 2010). When we inspect all DEGs identified at an FDR of 0.05 from the authors’ cell-level analysis, this number increases to 14,274. Pseudobulk DE analysis has recently been proven to give optimal performance compared to both mixed models and pseudoreplication approaches (Murphy and Skene, 2022; Squair et al., 2021; Crowell et al., 2020; Soneson and Robinson, 2018). It aggregates counts to individuals, thus accounting for the dependence between an individual’s cells.

Here, to compare the effect of the different DE approaches in isolation, we apply a pseudobulk DE approach (Chen et al., 2016) to the authors’ original processed data. We found 26 unique DEGs when considering the six cell types used by the authors (Table 2). This was 549 times fewer DEGs than that reported originally at an FDR of 0.05. When we compare these DEGs, we can see that the absolute log2 fold change (LFC) of our DEGs is 15 times larger than the authors’; median LFC of 2.34 and 0.16, despite the authors’ DEGs having an FDR score 8000 times smaller; median FDR of 2.89 × 10–7and 0.002 (Figure 1a and b). Although we examined a high correlation in the genes’ fold change values across our pseudobulk analysis and the authors’ pseudoreplication analysis (Pearson R of 0.87 for an adjusted p-value of 0.05, Table 3), the p-values and resulting DEGs vary considerably. The correspondence in fold change values is expected given the approaches are applied to the same dataset, whereas the probabilities, which pertain to the likelihood that a gene’s expressional changes is related to the case/control differences in AD, importantly do not align. We can show that this stark contrast is just an artefact of the authors taking cells as independent replicates and thus artificially inflating confidence by considering the Pearson correlation between the number of DEGs found and the cell counts (Figure 1c–e). There is a near perfect, positive correlation between DEG and cell counts for the authors’ pseudoreplication analysis (Figure 1c) and for the 1031 genes from the authors’ combinatorial approach (Figure 1d) which is not present in our pseudobulk reanalysis (Figure 1e).

Table 2
The differentially expressed genes from our reanalysis using the same processed data the authors used and pseudobulk differential expression approach.
CelllogFClogCPMLRp-Valueadj_pvalHGNC
Mic2.701789136.9979461926.14184153.17E-070.00061349ACRBP
Mic1.489300718.0624087728.63612178.73E-080.00019303APOC1
Mic1.093276698.6419976921.53230143.48E-060.00336416CD81
Mic–1.41576817.9388487523.99554679.66E-070.00135806CD83
Mic3.37827276.8618354832.08044011.48E-084.58E-05CLEC1B
Mic2.840724526.7437054221.77455093.07E-060.00316269EGF
Mic2.557696586.7834508718.04688722.16E-050.01699007ELOVL7
Mic–1.20560988.3319749922.66440451.93E-060.00229576IFI44L
Mic–1.66160697.1536663916.48012744.92E-050.03306938IFI6
Mic–1.98094257.0039628917.91808232.31E-050.01699007IFIT3
Mic2.765026726.7297880520.65436375.50E-060.00472825ITGA2B
Mic1.909634037.0155223316.32001895.35E-050.03448474MAP1A
Mic–1.81945088.2620888745.22210081.76E-111.36E-07NAMPT
Mic2.09450447.1104845620.80685245.08E-060.00462318NEXN
Mic–2.37897626.9389698522.39124412.22E-060.00245752NR4A2
Mic–2.85534626.7371386222.80298681.79E-060.00229576NR4A3
Mic3.328738296.8494272130.9553272.64E-086.81E-05PF4
Mic3.42139866.8732638333.26216578.05E-093.11E-05PKHD1L1
Mic3.645256776.9342217438.6612725.04E-102.60E-06PPBP
Mic2.304826798.1057044360.79326976.34E-159.81E-11PTPRG
Mic–1.03824688.1145026615.59682737.84E-050.04850839RORA
Mic2.546366496.6920298117.25326063.27E-050.02300507SDPR
Mic–0.96296178.843433417.93191312.29E-050.01699007SYTL3
Mic–1.42153747.9962980625.47362724.48E-070.00077092TMEM2
Mic2.989015966.7727664124.21008198.64E-070.00133637TUBB1
Opc–2.82747185.0337129222.13345812.54E-060.04176231EGR1
  1. CPM - Counts per Million, LR - fold change ratio, HGNC - HUGO Gene Nomenclature Committee.

Table 3
Pearson correlation between our pseudobulk differential expression analysis and the authors’ pseudoreplication analysis on all genes found to be significant at different adjusted p-value cut-offs from the authors’ pseudoreplication analysis.
Pseudoreplication adjusted p-value cut-offNumber of genes comparedPearson correlation
0.0120,1520.8646269
0.0523,9030.8708275
0.126,3820.8721126
0.2532,1170.8764692
0.542,0220.8751554
184,4670.826248

A further point which questions the authors’ DE approach is that they identified the vast majority of DEGs in the more abundant, neuronal cell types (Mathys et al., 2019). However, an increase in the number of cells is not the same as an increase in sample size since these cells are not independent from one another – they come from the same sample. Therefore, an increase in the number of cells should not necessarily result in an increase in the number of DEGs, whereas an increase in the number of samples would. This point is the major issue with pseudoreplication approaches which overestimate confidence when performing DE due to the statistical dependence between cells from the same patient not being considered (Squair et al., 2021; Lazic, 2010). In our opinion, it makes more sense to identify the majority of large effect size DEGs in microglia which recent work has established is the primary cell type by which the genetic risk for AD acts (Skene and Grant, 2016; McQuade and Blurton-Jones, 2019). This is what we found with our pseudobulk DE approach – 96% of all DEGs were in microglia (Table 2), whereas only 3% of the authors’ DEGs were in microglia.

Although it has been proven that pseudoreplication approaches result in false positives by artificially inflating the confidence from non-independent samples, we wanted to investigate the effect of the approach on the authors’ dataset. We ran the same cell-level analysis approach – a Wilcoxon rank-sum test and FDR multiple-testing correction – 100 times whilst randomly permuting the patient identifiers (Figure 1f). We would expect to find minimal DEGs with this approach given the random mixing of case and control patients. However, this pseudoreplication approach consistently found high numbers of DEGs, and we observe the same correlation between the number of cells and the number of DEGs as with the authors’ results. We did not observe the same pattern when running the same analysis with pseudobulk DE (Figure 1g). As a result, we conclude that integrating this pseudoreplication approach with a mixed model like the authors proposed just artificially inflates the test confidence for a random sample of the genes resulting in more false discoveries in cell types with bigger counts.

Up to this point, to compare the effect of the DE approaches in isolation, we analysed the same processed data from the authors as opposed to our reprocessed data. We also performed pseudobulk DE on our reprocessed data and found 16 unique DEGs (Table 4). It is worth noting that the fold change correlation between our two DE analyses (reprocessed data vs authors’ processed data) on the identified DEGs is only moderate (Pearson R of 0.57) and is lower than that of the correlation between pseudoreplication and pseudobulk on the same dataset (Table 3). This highlights the effect that the low quality high mitochondrial read cells have on DE analysis.

Table 4
The differentially expressed genes from our reanalysis using the reprocessed data and pseudobulk differential expression approach.
CelllogFClogCPMLRp-Valueadj_pvalensembl_idHGNC
OPC–4.15446634.9210080321.69114453.20E-060.04985906ENSG00000166573GALR1
Astro–4.58452764.796514322.23678472.41E-060.037634ENSG00000137959IFI44L
Micro–3.76166197.3287531626.81496882.24E-070.00077905ENSG00000077238IL4R
Micro–2.06814467.8873644117.59290952.74E-050.0346187ENSG00000105835NAMPT
Micro–1.67575567.5847250619.17368291.19E-050.02076348ENSG00000118257NRP2
Micro–3.15564036.8523265319.20646271.17E-050.02076348ENSG00000135363LMO2
Micro–3.43392656.929047219.59755899.56E-060.02076348ENSG00000138135CH25H
Micro–2.81831096.7750067616.9079593.92E-050.04550806ENSG00000142408CACNG8
Micro2.900766478.3456061745.51442661.52E-112.11E-07ENSG00000144724PTPRG
Micro3.258675896.9167101316.55191474.73E-050.0490155ENSG00000163106HPGDS
Micro–2.02909057.1232116616.47467464.93E-050.0490155ENSG00000171612SLC25A33
Micro–3.46573016.9330722119.78833018.65E-060.02076348ENSG00000172243CLEC7A
Micro–4.1728077.1681358334.35158074.60E-093.20E-05ENSG00000174600CMKLR1
Micro–3.19845886.8731055518.53358891.67E-050.0232342ENSG00000227531RP11-202G18.1
Micro3.405628876.938170318.55265021.65E-050.0232342ENSG00000228058RP11-552D4.1
Micro4.460733017.6655916329.77166794.86E-080.00022549ENSG00000253496RP11-13N12.1

In conclusion, the authors’ analysis has been highly influential in the field with numerous studies undertaken based on their results, something we show has uncertain foundations. However, we would like to highlight that the use of pseudoreplication in neuroscience research is not isolated to the authors’ work; others have used this approach (Fernandes et al., 2020; Lui et al., 2021; Wakhloo et al., 2020), and their results should be similarly scrutinised. Here, we provide our processed count matrix with metadata and also the DEGs identified using an independently validated, DE approach so that other researchers can use this rich dataset free from spurious nuclei or DEGs. While the number of DEGs found here is significantly lower, much greater confidence can be had that these are AD-relevant genes. The low number of DEGs found may also cause concern given the sample size and cost of collection and sequencing of such datasets. However, the increasing number of snRNA-seq studies being conducted for AD creates the opportunity to conduct differential meta-analyses to increase power. Further work is required in the field to develop methods to conduct such analysis, integrating studies and accounting for their heterogeneity, similar to that which has been done for bulk RNA-seq (Rau et al., 2014). Some such approaches have already been made in COVID-19 research which could be leveraged for neurodegenerative disease (Garg et al., 2021).

Materials and methods

Processing of sc/snRNA-seq dataset

The data reprocessing was conducted with scFlow (Khozoie et al., 2021), the steps of which are discussed in the following two sections.

Quality control of snRNAseq data

Request a detailed protocol

The raw snRNA-seq data (10.7303/syn18485175) and the ROSMAP metadata (10.7303/syn3157322) were downloaded from https://www.synapse.org/ upon acquiring appropriate approval. Downstream primary analyses of gene–cell matrices were performed using our scFlow pipeline (Khozoie et al., 2021). To determine ambient RNA profile and distinguish true nuclei from empty droplets, emptyDrops was used with a lower parameter of <100 counts, an alpha cut-off of ≤0.001, and with 10,000 Monte Carlo iterations (Lun et al., 2019). This approach has been recommended as best practice in the literature (Amezquita et al., 2020). Nuclei were then filtered for ≥200 total counts and ≥200 total expressed genes, which was defined as a minimum of 2 counts in at least three cells. We excluded any nuclei with total counts or total expressed genes with more than 4 MAD defined by an adaptive thresholding method. Nuclei were excluded if the proportion of counts mapping to mitochondrial genes was more than 10%, as set out in best-practice guidelines (Amezquita et al., 2020). Doublets were identified using the DoubletFinder algorithm, with a doublets-per-thousand-cells increment of eight cells (recommended by 10X Genomics), and a pK value of 0.005 (McGinnis et al., 2019). DoubletFinder was shown to be the best overall performing method in a recent benchmark (Xi and Li, 2021). The aggregated number of cells and proportions dropped at each step is given in Table 1 while a comparison of the proportion of cells in each cell type after reprocessing compared to the authors’ processed data is given in Figure 2. All files from the scFlow run, including QC statistics, are available in the GitHub repository in the scFlow_files folder (copy archived at Murphy, 2023). This includes sample-level genes and cells’ QC numbers.

Integration and clustering

Request a detailed protocol

The linked inference of genomic experimental relationships (LIGER) package was used to calculate integrative factors across samples (Welch et al., 2019). LIGER was recently found to be one of the top performing methods for batch-effect correction (Tran et al., 2020). LIGER parameters used included k: 30; lambda: 5.0; thresh: 0.0001; max_iters: 100; knn_k: 20; min_cells: 2; quantiles: 50; nstart: 10; resolution: 1; num_genes: 3000; and centre: false. Two-dimensional embeddings of the LIGER integrated factors were calculated using the Uniform Manifold Approximation and Projection (UMAP) algorithm with the following parameters: pca_dims: 50; n_neighbours: 35; init: spectral; metric: euclidean; n_epochs: 200; learning_rate: 1; min_dist: 0.4; spread: 0.85; set_op_mix_ratio: 1; local connectivity: 1; repulsion_strength: 1; negative_sample_rate: 5; and fast_sgd: false (McInnes et al., 2020). The Leiden community detection algorithm was used to detect clusters of cells from the 2D UMAP (LIGER) embeddings; a resolution parameter of 0.001 and a k value of 50 were used (Traag et al., 2019). This approach has been noted as best practice by a recent review (Heumos et al., 2023). Automated cell typing of the detected clusters was performed as previously described using the Expression Weighted Celltype Enrichment algorithm in scFlow against a previously generated cell-type data reference from the Allen Human Brain Atlas (Hodge et al., 2019; Skene and Grant, 2016). The top five marker genes for each automatically annotated cell type were determined using Monocle 3 and validated against canonical cell-type markers (Trapnell et al., 2014).

DE analysis

Request a detailed protocol

All DE analyses were conducted using pseudobulk DE approach with sum aggregation and edgeR LRT (Chen et al., 2016). Pseudobulk aggregates nuclei within a biological replicate (an individual) for each cell type, reducing the dropout issue in single-cell data and avoiding the false inflation of confidence from non-independent samples of pseudoreplication approaches (Squair et al., 2021; Murphy and Skene, 2022). The DE analysis pipeline is available at GitHub repository (copy archived at Murphy, 2023). This is a general use pipeline which can be run for any single-nucleus or single-cell transcriptomic dataset. Note that we report DEGs across AD and controls using the same processed data the authors used (Table 2) and using our reprocessed data (Table 4).

Code availability

Request a detailed protocol

The DE analysis pipeline is available at GitHub repository (copy archived at Murphy, 2023). This is a general use pipeline which can be run for any single-nucleus or single-cell transcriptomic dataset. The config file containing all the parameters used and QC overview file for the scFlow run is also available in this repository.

Data availability

The differentially expressed genes and processed count matrix from the original study are available with their manuscript. The count matrix and metadata from our reprocessing approach are available via the AD Knowledge Portal (https://adknowledgeportal.org). The AD Knowledge Portal is a platform for accessing data, analyses, and tools generated by the Accelerating Medicines Partnership (AMP-AD) Target Discovery Program and other National Institute on Aging (NIA)-supported programs to enable open-science practices and accelerate translational learning. The data, analyses, and tools are shared early in the research cycle without a publication embargo on secondary use. Data is available for general research use according to the following requirements for data access and data attribution (https://adknowledgeportal.org/DataAccess/Instructions). For access to content described in this article, see https://doi.org/10.7303/syn51758062.1. All other relevant scripts and data for working with this dataset and supporting the key findings of this study are available within the article or from our GitHub repository (copy archived at Murphy, 2023).

The following previously published data sets were used

References

Peer review

Joint Public Review:

Murphy, Fancy and Skene performed a reanalysis of snRNA-seq data from Alzheimer Disease (AD) patients and healthy controls published previously by Mathys et al. (2019), arriving at the conclusion that many of the transcriptional differences described in the original publication were false positives. This was achieved by revising the strategy for both quality control and differential expression analysis. With this re-analysis, the authors aim to raise awareness of the impact of data analysis choices for scRNA-seq data and to caution focus on putatively wrongly identified genes in the AD research community. The revised manuscript has been improved by separating QC and DE analysis, which makes interpretation of both steps more straightforward.

STRENGTHS:

The authors demonstrate that the choice of data analysis strategy can have a vast impact on the results of a study, which in itself may not be obvious to many researchers.

The authors apply a pseudobulk-based differential expression analysis strategy (essentially, adding up counts from all cells per individual and comparing those counts with standard RNA-seq differential expression tests), which is (a) in line with latest community recommendations, (b) different from the "default options" in most popular scRNA-seq analysis suites, and (c) explains the vastly different number of DEGs identified by the authors and the original publication. The recommendation of this approach together with a detailed assessment of the DEGs found by both methodologies could potentially be a useful finding for the research community. Unfortunately, it is currently not sufficiently substantiated.

All code and data used in this study are publicly available to the readers.

WEAKNESSES:

The authors interpret the fact that they found fewer DEGs with their method than the original paper as a good thing by making the assumption that all genes that were not found were false positives. However, they do not prove this, and it is likely that at least some genes were not found due to a lack of statistical power and not because they were actually "incorrect". The original paper also had performed independent validations of some genes that were not found here. I had raised this weakness in my first review, but it was not explicitly addressed and still pertains to the revised manuscript. The authors have added an analysis that shows that "pseudoreplication" is prone to false positive (FP) discoveries for high cell numbers (Fig. 1f), but this does not prove that all of Mathys' DEGs were wrong.

I am concerned that almost all DEGs found by the authors are in the rare cell types, foremost the rare microglia (see Fig. 1e). Indeed, there is a weak negative correlation between cell counts and numbers of DEGs (Fig. 1e), if the correlation analysis is to be believed (see next point). It is unclear to me how many cells the pseudo-bulk counts were based on for these cell types, but it seems that (a) there were few and (b) there were quite few reads per cells. If both are the case, the pseudobulk counts for these cell populations might be rather noisy and the DEG results liable to outliers with extreme fold changes. Supp. Fig. 3b now shows three examples of DEGs, of which one (EGR1) looks like the DE call is indeed largely driven by four outliers, while Supp. Fig 3a shows at least one gene (BEX1) that could be FP of the pseudobulk approach due to insufficient statistical power. The authors go on to cite two papers (one is their own, published in a journal with suspected lack of appropriate quality assurance measures https://predatoryreports.org/the-predatory-journals-1), to support that the finding of DEGs in microglia "makes more sense" (l. 127). In summary, neither the presented examples nor the supporting literature are convincing. Lastly, the authors even show themselves that their approach is liable to FPs if applied with very low cell numbers in the range of those for microglia and OPCs (Fig. 1g).

The correlation analysis between cell counts and number of DEGs found is weak. In all three cases (Fig. 1c, d, e) the correlation is largely driven by a single outlier data point.

The authors claim they improved the quality control of the dataset but offer no objective metric to assess this putative improvement. The authors' QC procedure removes some 20k cells that had not been filtered out by Mathys' et al. As the authors state themselves, this difference is mostly due to the removal of cells with a high mitochondrial read content. Murphy et al use a fixed threshold for the mitochondrial percentage of reads, while the original paper had removed cell clusters with an "abnormally high" mitochondrial read fraction. That also seems reasonable, given that some cells might have a higher mitochondrial read content for reasons other than being "low quality". Simply stating that Mathys' approach was ineffective at removing cells with high mitochondrial read content is a self-fulfilling prophecy given the difference in approach, and itself not proof that the original QC procedure was inferior.

Batch correction: "Dataset integration has become a common step in single-cell RNA-Seq protocols and is recommended to remove confounding sources of variation" (l. 38). While it is true that many authors now choose to perform an integration step as part of their analysis workflow, this is by no means uncontroversial as there is a risk of "over-integration" and loss of true biological differences. I had raised this point previously, but the authors chose not to address it (quoted text and line numbers updated). Given that there is controversy in the literature and "community opinion" on the topic of data integration, this is another example of the authors claiming superiority in analysis without showing proof.

Due to a lack of comparison with other methods and due to the fact that the author's methodology was only applied to a single dataset, the paper presents merely a case study, which could be useful but falls short of providing a general recommendation for a best practice workflow.

APPRAISAL:

The manuscript could help to increase awareness of data analysis choices in the community, but only if the superiority of the methodology was clearly demonstrated. However, the authors only show that there are differences but have no convincing (orthogonal) evidence that their methodology was indeed better. This applies to both QC and DE analysis.

https://doi.org/10.7554/eLife.90214.3.sa1

Author response

The following is the authors’ response to the original reviews.

Response to Reviewers

To whom it may concern,Thank you for your constructive feedback on our manuscript. I appreciate the time and effort that you and the reviewers have dedicated to providing your valuable feedback. We are grateful to the reviewers for their insightful comments and suggestions for our paper. I have been able to incorporate changes to reflect the majority of these suggestions provided. I have updated the analysis scripts (at https://github.com/neurogenomics/reanalysis_Mathys_2019) and have listed these changes in blue below:

eLife assessment:

This work is useful as it highlights the importance of data analysis strategies in influencing outcomes during differential gene expression testing. While the manuscript has the potential to enhance awareness regarding data analysis choices in the community, its value could be further enhanced by providing a more comprehensive comparison of alternative methods and discussing the potential differences in preprocessing, such as scFLOW. The current analysis, although insightful, appears incomplete in addressing these aspects.

We thank the reviewing editors for this note. We agree that the differences in preprocessing will affect the results and conceal which step in our reanalysis resulted in the discrepancies we noted. To address this, we have split out our reanalysis into two separate parts - In the main body of the text we discuss the differences resulting from just changing the differential expression approach where we use the same processed data as the authors to enable a fair comparison. Secondly, we still provide the reprocessed data and perform differential expression analysis on it and discuss the cause and impact the differences in the processing steps made to the results.

Reviewer 1:

I think readers would be interested to learn more about the genes that were found "significant" by the original paper but sorted out by the authors. Did they just fall short of the cutoffs? If so, how many more samples would have been required to ascertain significance? This would yield a recommendation for future studies and an overall more positive/productive spirit to the manuscript. On the other hand, I suspect a fraction of DEGs were false positives due to differences in the proportions of cells from different individuals compared to the original analysis. Which percentage of DEGs does this apply to? Again, this would raise awareness of the issue and support the use of pseudobulk approaches.

To investigate the relationship between the genes and how they differ across our analysis we have added a correlation analysis between our different DE approaches (using the same processed data), see paragraph 5 in the manuscript and supplementary table 3. In short, we find that there is a high correlation in the genes’ fold change values across our pseudobulk analysis and the author’s pseudoreplication analysis on the same dataset (pearson R of 0.87 for an adjusted p-value of 0.05) which is somewhat expected given the DE approaches are applied to the same dataset. However, the p-values, which pertain to the likelihood that a gene’s expressional changes is related to the case/control differences in AD, and resulting DEGs vary considerably due to the artificially inflated confidence of the author’s approach (Fig. 1c-e).Despite there being a correlation between the pseudoreplciation and pseudobulk approaches here, we do not think it makes sense to consider how many more samples would have been required to ascertain significance. The differences in results between the two approaches is not negatable with sample size as many DEGs identified by pseudoreplication will be false positives as highlighted in previous work1,2,3,4.However, perhaps we are misinterpreting the reviewer, who may have meant a power analysis which we have not conducted. Such an undertaking would require analysing a multitude of snRNA-Seq of large sample sizes to garner a confident estimate for power calculations based on pseudobulk approaches. Although we agree with the reviewer that this would be beneficial to the field, we do not believe it is in scope for this work.On the reviewer’s note regarding a fraction of DEGs being false positives due to differences in the proportions of cells from different individuals compared to the original analysis - We have analysed the same processed data the authors used to negate the differences caused by the differing processing steps. We thank the reviewer for this suggestion. We also give more insight into the cause of these differences, namely on filtering our nuclei with large proportions of mitochondrial reads and discuss their effect in paragraph 3 (also see Supplementary Figure 2).

Given there are only a few DEGs, it would be good to show more data about these genes to allow better assessment of the robustness of the results, i.e., boxplots of the pseudobulk counts in the compared groups and perhaps heatmaps of the raw counts prior to aggregation. This could rule out concerns about outliers affecting the results.

In Supplementary Figure 3, we have added boxplots of the sum pseudobulked, trimmed mean of M-values (TMM) normalised counts for three of our identified DEGs (b) and three of the authors’ DEGs which they discuss in their manuscript (a) to show the differences in counts across AD pathology and controls for these genes. We hope this gives some insight into the transcriptional changes highlighted by the differing approaches. In our opinion, there is a clear difference in the transcriptional signal in the genes identified from pseudobulk which is not present for the genes identified from the authors approach.

Overall, I believe the paper would deliver a clearer message by mainlining the QC from the original study and only changing the DE analysis. However, if keeping the part about QC/batch correction:

  • Assess to which degree changes in cell type proportion are indeed due to batch correction (as suggested in the text) and not filtering by looking at the annotated cell types in the original publication and those in your analysis.

  • Also perform the analysis without changing QC and state the # of DEGs in both cases, to at least allow some disentanglement of the effect of different steps of the analysis.

  • Please state the number of cells removed by each QC step in the supplementary note.

We thank the reviewer for this suggestion. We agree with performing the DE analysis on the same processed data as the original authors and have split out our reanalysis into two separate parts, primarily focussing on the discrepancies caused by the choice of differential expression (DE) approach. By splitting our analysis in this manner, we can identify the substantial differences in results caused by differing the DE approach in the study. Secondly, we can see how differences in preprocessing affects the DE results in isolation too – see paragraph 8 but in short, the fold change correlation between pseudobulk DE analyses on the reprocessed data vs authors processed data only had a moderate correlation (Pearson R of 0.57).

In regards to the number of cells removed by each QC step, we have added an aggregated view for all samples in supplementary table 3 and also give the full statistics per sample in our Github repository: https://github.com/neurogenomics/reanalysis_Mathys_2019. Moreover, we investigated the root cause in the differences in nuclei numbers, uncovering filtering down to mitochondrial read proportions as the main culprit (Supplementary Figure 2).

I recommend the authors read the following papers, assess whether their methodology agrees with them, and add citations as appropriate to support statements made in the manuscript.

We thank the reviewer for this comprehensive list. We have updated our manuscript and supplementary file and main text throughout to cite many of these where appropriate. We believe this helps add context to our decisions for the differing tools and approaches used as part of the processing pipeline with scFlow and the differential expression approach.

I believe the authors' intention was to show the results of their reanalysis not as a criticism of the original paper (which can hardly be faulted for their strategy which was state-of-the-art at the time and indeed they took extra measures attempting to ensure the reliability of their results), but primarily to raise awareness and provide recommendations for rigorous analysis of sc/snRNA-seq data for future studies.

We thank the reviewer for this note, this was exactly our intent. Furthermore, we are based in a dementia research institute and our aim is to ensure that ensure that the Alzheimer’s disease research field does not focus on spuriously identified genes.We have updated the text of the manuscript (start paragraph 2) to explicitly state this so our message is not misconstrued.

In my opinion, the purpose of the paper might be better served by focusing on the DE strategy without changing QC and instead detailing where/how DEGs were gained/lost and supporting whether these were false positives.

We agree that the differences in preprocessing will affect the results and conceal which step in our reanalysis resulted in the discrepancies we noted. To address this, we have split out our reanalysis into two separate parts - In the main body of the text we discuss the differences resulting from just changing the differential expression approach where we use the same processed data as the authors to enable a fair comparison. Secondly, we still provide the reprocessed data and perform differential expression analysis on it and discuss the impact the differences in the processing steps made to the results. As previously mentioned, we have also added further investigation into the DEGs identified, looking at the correlation across the differing approaches and plotting the counts for selected genes.

For instance, removal with a mitochondrial count of <5% seems harsh and might account for a large proportion of additional cells filtered out in comparison to the original analysis. There is no blanket "correct cutoff" for this percentage. For instance, the "classic" Seurat tutorial https://satijalab.org/seurat/articles/pbmc3k_tutorial.html uses the 5% threshold chosen by the authors, an MAD-based selection of cutoff arrived at 8% here https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html, another "best practices" guide choses by default 10% https://bioconductor.org/books/3.17/OSCA.basic/quality-control.html#quality-control-discarded, etc. Generally, the % of mitochondrial reads varies a lot between datasets.

Apologies, the 5% cut-off was a misprint – the actual cut-off used was 10% which, as the reviewer notes, is on the higher side of what is recommended. We have updated our manuscript to rectify this mistake and discuss the differences in the number of cells caused by the two approaches to mitochondrial filtering in the manuscript (paragraph 3). We found that over 16,000 nuclei that were removed in our QC pipeline were kept by the author’s (Supplementary Fig. 2), explaining the discrepancy in the number of nuclei after QC. Based on Supplementary Fig. 2, it is clear the author’s approach was ineffective at removing nuclei with high proportions of mitochondrial reads which is indicative of cell death5,6. We hope this alleviates the reviewer’s concerns around our alternative processing approach. Moreover, as mentioned, we swapped to compare the differences by DE approaches on the same data to avoid any effect by this.

Reviewer 2:

The paper would be better if the authors merged this work with the scFLOW paper so that they can justify their analysis pipeline and show it in an influential dataset.

We thank the reviewer for this note. We would like to clarify that the purpose of our work was not to show the scFlow analysis pipeline on an influential dataset but rather to raise awareness and provide recommendations for rigorous analysis of single-cell and single-nucleus RNA-Seq data (sc/snRNA-Seq) for future studies and to help redirect the focus of the Alzheimer’s disease research field away from possible spuriously identified genes. We have updated our manuscript text to highlight this (see start paragraph 2). Furthermore, we are aware our original approach reprocessing the data with scFlow will affect the results and conceal which step in our reanalysis resulted in the discrepancies we noted. Thus, we have split out our reanalysis into two separate parts - In the main body of the text we discuss the differences resulting from just changing the differential expression approach where we use the same processed data as the authors to enable a fair comparison. Secondly, we still provide the reprocessed data so that the community can benefit from it and perform differential expression analysis on it and discuss the impact the differences in the processing steps made to the results. We have also added further references supporting the choice of steps and tools used in scFlow in the supplementary text which should address the reviewer’s concerns about justifying the analysis pipeline. Moreover, we identified the cause of the nuclei count differences caused by the two processing approaches, namely on filtering our nuclei with large proportions of mitochondrial reads and discuss their effect in paragraph 3 (also see Supplementary Figure 2).

A major contribution is the use of the authors' own inhouse pipeline for data preparation (scFLOW), but this software is unpublished since 2021 and consequently not yet refereed. It isn't reasonable to take this pipeline as being validated in the field.

We believe our answer to the previous point addresses these concerns - We have added references supporting the choice of steps and tools used in scFlow in the supplementary text which should address the reviewer’s concerns about justifying the analysis pipeline. Moreover, as a result of the pipeline we identified that 16,000 of the nuclei kept by the authors are likely of low quality and indicative of cell death with high mitochondrial read proportions5,6.

They also worry that the significant findings in Mathys' paper are influenced by the number of cells of each type. I'm sure it is since power is a function of sample size, but is this a bad thing? It seems odd that their approach is not influenced by sample size.

We thank the reviewer for highlighting this point. As they noted, we conclude that the original authors number of DEGs is just a product of the number of cells. However, the reviewer states that ‘It seems odd that their approach is not influenced by sample size’. An increase in the number of cells is not an increase in sample size since these cells are not independent from one another - they come from the same sample. Therefore, an increase in the number of cells should not result in an increase in the number of DEGs whereas an increase in the number of samples would. This point is the major issue with pseudoreplication approaches which over-estimate the confidence when performing differential expression due to the statistical dependence between cells from the same patient not being considered. See these references for more information on this point1,2,7,8. We have added a discussion of this point to our manuscript in paragraph 6.

Moreover, recent work has established that the genetic risk for Alzheimer’s disease acts primarily via microglia9,10. Thus, it would be reasonable to expect that the majority of large effect size DEGs identified would be found in this cell type. This is what we found with our pseudobulk differential expression approach – 96% of all DEGs were in microglia. We have updated the text of our manuscript (paragraph 5) to highlight this last point.

References

1. Murphy, A. E. & Skene, N. G. A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. Nat. Commun. 13, 7851 (2022).

2. Squair, J. W. et al. Confronting false discoveries in single-cell differential expression. Nat. Commun. 12, 5692 (2021).

3. Crowell, H. L. et al. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat. Commun. 11, 6077 (2020).

4. Soneson, C. & Robinson, M. D. Bias, robustness and scalability in single-cell differential expression analysis. Nat. Methods 15, 255–261 (2018).

5. Ilicic, T. et al. Classification of low quality cells from single-cell RNA-seq data. Genome Biol. 17, 29 (2016).

6. Heumos, L. et al. Best practices for single-cell analysis across modalities. Nat. Rev. Genet. 24, 550–572 (2023).

7. Zimmerman, K. D., Espeland, M. A. & Langefeld, C. D. A practical solution to pseudoreplication bias in single-cell studies. Nat. Commun. 12, 738 (2021).

8. Lazic, S. E. The problem of pseudoreplication in neuroscientific studies: is it affecting your analysis? BMC Neurosci. 11, 5 (2010).

9. Skene, N. G. & Grant, S. G. N. Identification of Vulnerable Cell Types in Major Brain Disorders Using Single Cell Transcriptomes and Expression Weighted Cell Type Enrichment. Front. Neurosci. 0, (2016).

10. McQuade, A. & Blurton-Jones, M. Microglia in Alzheimer’s disease: Exploring how genetics and phenotype influence risk. J. Mol. Biol. 431, 1805–1817 (2019).

https://doi.org/10.7554/eLife.90214.3.sa2

Article and author information

Author details

  1. Alan E Murphy

    1. UK Dementia Research Institute at Imperial College London, London, United Kingdom
    2. Department of Brain Sciences, Imperial College London, London, United Kingdom
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft
    For correspondence
    alanmurph94@hotmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2487-8753
  2. Nurun Fancy

    1. UK Dementia Research Institute at Imperial College London, London, United Kingdom
    2. Department of Brain Sciences, Imperial College London, London, United Kingdom
    Contribution
    Data curation, Software, Formal analysis, Validation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6481-6266
  3. Nathan Skene

    1. UK Dementia Research Institute at Imperial College London, London, United Kingdom
    2. Department of Brain Sciences, Imperial College London, London, United Kingdom
    Contribution
    Conceptualization, Supervision, Methodology, Project administration, Writing - review and editing
    For correspondence
    n.skene@imperial.ac.uk
    Competing interests
    No competing interests declared

Funding

UK Research and Innovation (Future Leaders Fellowship (MR/T04327X/1))

  • Nathan Skene

UK Dementia Research Institute

  • Nathan Skene

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by a UKDRI Future Leaders Fellowship (grant number MR/T04327X/1) and the UK Dementia Research Institute, which receives its funding from UK DRI Ltd, funded by the UK Medical Research Council, Alzheimer’s Society and Alzheimer’s Research UK. The results published here are in whole or in part based on data obtained from the AD Knowledge Portal (https://adknowledgeportal.org). The data available in the AD Knowledge Portal would not be possible without the participation of research volunteers and the contribution of data by collaborating researchers. Study data were provided by the Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago. Data collection was supported through funding by NIA grants P30AG10161 (ROS), R01AG15819 (ROSMAP; genomics and RNAseq), R01AG17917 (MAP), R01AG30146, R01AG36836 (RNA-seq), U01AG32984 (genomic and whole-exome sequencing), U01AG46152, U01AG61356 (ROSMAP AMP-AD, targeted proteomics), U01AG46161 (TMT proteomics), U01AG61356 (whole genome sequencing, targeted proteomics, ROSMAP AMP-AD), the Illinois Department of Public Health (ROSMAP), and the Translational Genomics Research Institute (genomic). Additional phenotypic data can be requested at https://www.radc.rush.edu/.

Senior Editor

  1. Murim Choi, Seoul National University, Republic of Korea

Reviewing Editor

  1. Joon-Yong An, Korea University, Republic of Korea

Version history

  1. Preprint posted: June 14, 2023 (view preprint)
  2. Sent for peer review: July 11, 2023
  3. Preprint posted: September 1, 2023 (view preprint)
  4. Preprint posted: November 16, 2023 (view preprint)
  5. Version of Record published: December 4, 2023 (version 1)

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.90214. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2023, Murphy et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 909
    Page views
  • 83
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Alan E Murphy
  2. Nurun Fancy
  3. Nathan Skene
(2023)
Avoiding false discoveries in single-cell RNA-seq by revisiting the first Alzheimer’s disease dataset
eLife 12:RP90214.
https://doi.org/10.7554/eLife.90214.3

Share this article

https://doi.org/10.7554/eLife.90214

Further reading

    1. Cell Biology
    2. Chromosomes and Gene Expression
    Carolline Ascenção, Jennie R Sims ... Marcus B Smolka
    Research Article

    Meiotic sex chromosome inactivation (MSCI) is a critical feature of meiotic prophase I progression in males. While the ATR kinase and its activator TOPBP1 are key drivers of MSCI within the specialized sex body (SB) domain of the nucleus, how they promote silencing remains unclear given their multifaceted meiotic functions that also include DNA repair, chromosome synapsis, and SB formation. Here we report a novel mutant mouse harboring mutations in the TOPBP1-BRCT5 domain. Topbp1B5/B5 males are infertile, with impaired MSCI despite displaying grossly normal events of early prophase I, including synapsis and SB formation. Specific ATR-dependent events are disrupted, including phosphorylation and localization of the RNA:DNA helicase Senataxin. Topbp1B5/B5 spermatocytes initiate, but cannot maintain ongoing, MSCI. These findings reveal a non-canonical role for the ATR-TOPBP1 signaling axis in MSCI dynamics at advanced stages in pachynema and establish the first mouse mutant that separates ATR signaling and MSCI from SB formation.

    1. Chromosomes and Gene Expression
    Masaaki Sokabe, Christopher S Fraser
    Insight

    A new in vitro system called Rec-Seq sheds light on how mRNA molecules compete for the machinery that translates their genetic sequence into proteins.