Droplet-based single-cell RNA-sequencing methods measure gene expression on tens or hundreds of thousands of cells at the single-cell level. Gene expression measurements in droplet technology are often in the form of low counts with a large fraction of zero values, and difficult to estimate the exact statistical distribution. We decomposed the gene expression distribution into two parts (Figure 1B): the status of genes (“on/off”) and the distribution shape of gene “on” parts (non-zero part). These two parts were closely associated with cell type, cell condition, or other biologic-driven factors. For the first part, the gene expression state ratio was defined as the times of the gene with the positive count in a group of cells. For the gene “On” part, we described the distribution shape by location parameter ( μ ) and scale parameter ( σ ) of the “On” parts (Figure 1B). Therefore, the whole gene expression pattern can be approximated by three parameters: the zero proportion of gene expression ( p ), the mean of the “On” parts ( μ ), and the variance of the “On” parts ( σ ).We assumed that non-DE genes have the same expression distribution shape in pre-defined groups (Figure 1A). We tested three parameters ( H 0 : p j 1 = p j 2 ,   μ j 1 = μ j 2 ,   σ j 1 2 = σ j 2 2 ) to identify whether a given gene is a DE gene (Figure 1C). Due to low counts, sparsity, and complexity of gene expression, it is challenging to estimate the exact distribution of every gene and construct a suitable statistic for the hypothesis H 0 when the theoretical distribution of genes is unknown. Instead of generating the test statistic based on the assumed distribution, we tested the complex null hypothesis H 0 using Fisher’s (Zappia et al., 2017) theory of combination test. H 0 :   θ = ( p j 1 ,   p j 2 ,   μ j 2 , μ j 2 . σ j 1 , σ j 2 ) ∈ Θ 0 ,   Θ 0 = { θ ∈ Θ : p j 1 = p j 2 , μ j 1 = μ j 2 ,     σ j 1 = σ j 2 } H A :   θ ∈ Θ 0 C { H 01 : p j 1 = p j 2 ;   H A 1 : p j 1 ≠ p j 2 H 02 : μ j 1 = μ j 2 ;   H A 2 : μ j 1 ≠ μ j 2 H 03 : σ j 1 2 = σ j 2 2 ;   H A 3 : σ j 1 2 ≠ σ j 2 2

We split the complex null hypothesis H 0 into three simple null hypotheses H 0 i and got a new statistic Q = − 2 ∑ log ( L i ) by combining three individual p-values L i . Each p-value L i was obtained by testing the simple null hypothesis H 0 i . The chi-square distribution was used to approximate the p-value of Q . Underlying this test framework, we easily captured various differences in gene expression and constructed a test for gene expression patterns without many assumptions. Moreover, we only calculated three simple observed test statistics and got the new statistic Q by combining three individual p-values L i . We could quickly identify differential expression (DE) genes in millions-scale scRNA data. The computation cost is almost negligible. If the new statistic Q is larger than the critical value, we reject the null hypothesis and identify the gene as a DE gene. We examined one gene at a time and implemented FDR correction for p-values of all genes.

HEART proposed a combination test to catch various sources of differences in gene expression patterns between two pre-defined groups. To validate the performance of HEART, we used two simulation data generation mechanisms to compare HEART and other six popular DE methods, including five model-based DE methods (DESeq2, MAST, Monocle3, NBID, and SC2P) and a default test in Seurat (Seurat-W). Simulation details were provided in the “Methods”. Briefly, the artificial simulation tool, Splatter package (Zappia et al., 2017), generated datasets in simulation1. Simulation2 datasets used a semi-simulation mechanism based on actual scRNA-seq data (PBMC68K)to create simulation datasets. In both simulations, we varied the number of samples and DE strength for DE genes. We evaluated the ability to identify DE genes, FDR control under the null hypothesis, and computational efficiency under various alternatives by a series of indexes: F 1 score, TPR, precision, computational time, etc.

In simulation 1, we evaluated the performances of each method on simulation datasets with the same simulation settings. HEART, Monocle3, and NBID perform better than other methods (Figure 2A; Supplementary Figure S2). They had higher F 1 scores than other methods and achieved a good balance between TPR and precision. Seurat had low precisions, because it was apt to identify the gene with mild signals. DESeq2 maintained high accuracy on medium-scale data (under 10000 cells), but it shows FDR inflation on the large-scale datasets (Supplementary Figure S2). Regarding running time, HEART and Seurat had incomparable advantages (Figure 2D, under 2 min on the datasets of 20000 cells). Although NBID and DESeq2 had good accuracy, they required a lot of running time (Figure 2D, more than 1 h on the datasets of 20000 cells with 11000 genes).

In Simulation2, we generated semi-simulation data from real scRNA-seq datasets instead of simulation datasets from artificial protocols (Figure 2B; Supplementary Figure S3) (Chen et al., 2018). We chose each cell subtype with various sample sizes from PBMC68K (Zheng et al., 2017) as source data to test the stability and scalability of each DE method. HEART, NBID, and Monocle3 have higher F 1 scores in different simulation datasets than other methods. When the sample size was adequate, HEART had good and stable performances, regardless of the statistical characteristics of the datasets. Seurat performed unstably on different datasets. DESeq2, MAST, and SC2P cannot detect DE genes in most scenarios. Importantly, HEART was much more computationally efficient than the other methods (Figure 2E). For the 20000-cells scale datasets, HEART completed computation in about 1–2 min, but NBID and DESeq2 needed 5–7 h for the same scale datasets. HEART was applicable to data with the sample size exceeding around millions of cells in theory. We generated null simulations without swapping genes to test the bias in p-value estimation for each method (Supplementary Figure S4). HEART controlled the type 1 error well.

Generally, HEART was an accurate, practical and scalable method for DE gene detection. In all semi-simulation scenarios, HEART and NBID performed better than other methods and had relatively stable performances on datasets with various characteristics. Other methods had poor performances on some semi-simulation datasets. As the sample size increases, the performances of HEART, NBID, and Monocle3 become better. However, HEART identified DE genes in the simulation scenarios with weak DE strength of differences, which means HEART was more sensitive than other competing DE methods (Supplementary Figure S3; Figure 3). The performance of NBID was slightly better than HEART in some scenarios, but it took a lot of time to run. (Simulation1 of 20000 cells: NBID: F 1 score = 0.871 running time = 6482 s; HEART: F 1 score = 0.84, running time = 52 s. Simulation2 of CD8⁺ cytotoxic T cells: NBID: F 1 score = 0.97, running time = 16205 s; HEART: F 1 score = 0.94, running time = 94 s)

Read count and unique molecular identifier (UMI) count are two main quantification schemes in single-cell RNA-sequencing technologies and have different statistical characterizations. Some literature (Zilionis et al., 2017; Chen et al., 2018; Kashima et al., 2020; Sarkar and Stephens, 2021) suggested that read count data have higher count levels, more sparsity and more variability than UMI counts data. To assess the accuracy and robustness of HEART on different quantification mechanisms, we applied HEART and other six DE methods (Seurat, DESeq2, MAST, Monocle3, NBID, and SC2P) on two real single-cell datasets from quantification schemes. A human brain dataset (Darmanis et al., 2015) (GSE67835) based on read count quantification schemes and a dataset of peripheral blood mononuclear cells (PBMC68K (Zheng et al., 2017)) quantified by UMI counts.

Colorectal cancer (CRC) is the most commonly diagnosed cancers in the world.20% of individuals with newly diagnosed colorectal cancer have metastatic disease upon presentation, and another 25% of those who initially have localized illness will eventually acquire metastases (Biller and Schrag, 2021). Distant metastasis was the main cause of death in patients with colorectal cancer, but the exact metastasis mechanism was still unknown. (Zhang et al., 2014). ScRNA-seq technology provided a new opportunity to investigate the association between genes and the mechanism of tumor initiation, progression, and metastasis (Lawson et al., 2018). Therefore, we applied HEART in a single-cell dataset (containing three sub-datasets: PBMC, normal tissue, and tumor tissue) of a stage III colorectal cancer patient. We used HEART to identify DE genes between tumor and normal fibroblasts and between tumor and normal epitheliums. Furthermore, we found two subpopulations of megakaryocytes (MKs) (Wang et al., 2021) in the PBMCs and utilized HEART to detect 207 DE genes on the2 MK subtype clusters to characterize functional differences and underlying molecular mechanisms. Highly expressed genes in the cluster MK3 (Satija et al., 2015; Stuart et al., 2019; Fa et al., 2021; Wang et al., 2021), such as CCL5, TUBB1, MYL9, HIST1H2AC, etc. (Figure 4A), were associated with early platelet production. Another subpopulation, MK5, with high CD74 and PLAC8 might be a less mature MK population. Moreover, we observed that many DE genes between MK3 and MK5 cells overlap with DE genes between tumor and normal epitheliums and DE genes between tumor and normal fibroblasts (Figure 4B , Figure 4C). They had similar expression patterns in the MK5 cells, tumor epitheliums, and fibroblasts (Figure 4D) and were related to colorectal cancer progression or metastasis. The Violin plot showed similar distribution shapes of CTTN in MK5 cells and epithelial tumor cells. The gene CTTN has been reported overexpressed in various cancers, including colorectal cancer, and had the function of promoting tumor cell migration (Luo et al., 2006; Jing et al., 2016; Zhang et al., 2017). Furthermore, S100A4 (Helfman et al., 2005; Nader et al., 2020), S100A6 (Komatsu et al., 2000), UBA52 (Zhou et al., 2019), FAU (Pickard et al., 2011), and VIM (Luque-Garcia et al., 2010; Xu et al., 2017), etc. Also had similar expression patterns between MK5 cells and tumor epitheliums and fibroblasts. S100A4 and S100A6 play an important role in tumor metastases, including colorectal tumor metastasis (Komatsu et al., 2000). Recent studies have proved that the bloodstream plays a crucial role in tumor metastasis and tumor immune escape (Lawson et al., 2018). The cooperation of hematopoiesis, megakaryocytes, and platelet production aided CTCs in escaping the immune system and disseminating within the bloodstream to establish distant organ metastasis. We also validated the expression pattern of these genes in the spatial transcriptome data of two other stage IV colorectal cancer patients (Supplementary Figure S4), which showed spatial patterns of high expression in cancer cells.

Consequently, we supposed that a series of genes, CTTN, S100A4, S100A6, etc., were potential colorectal cancer metastasis biomarkers. The MK5 subpopulation with highly-expressed above potential biomarkers might be a cluster related to colorectal cancer metastasis and have a circulating tumor cell (CTC). The exact mechanism between MK5 and colorectal tumor metastasis warranted further investigation.

We used three actual scRNA-seq datasets in applications. The PBMC68K is available from https://support.10xgenomics.com/single-cellgene-expression/datasets. The human brain dataset can be obtained by R package SC2P or the GEO database repository under accession code GSE67835. The scRNA-seq data of one stage III colorectal cancer patient have been deposited in the OMIX, China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac.cn/omix: accession no. OMIX002120). The spatial transcriptomic data of two colorectal cancer patients are available from (http://www.cancerdiversity.asia/scCRLM).

We used two simulation data generation mechanisms to generate scenarios with different settings. Each design had 20 replications. The popular artificial protocol, Splatter (Zappia et al., 2017), generated simulation datasets in simulation 1. Each scenario contained 10000 genes (1000 DE genes and 10000 non-DE genes) and two underlying subpopulations. We varied the number of samples (1,000, 2000, 5,000, 10000, 20000) and DE strength for DE genes (de.factor = 0.3, 0.5). De.factor is the differential expression factor produced from a log-normal distribution. A high de.factor can result in the strong DE strength of DE genes between groups (More details of parameters in Supplementary Material S1).

Simulation 2 adopted a semi-simulation mechanism based on actual scRNA datasets to recover the multimodality and biological characteristic complexity of actual scRNA-seq data (Figure 2C, Supplementary Material S1) (Chen et al., 2018). First, we randomly divided the real scRNA-seq dataset into two parts regarded as two groups of cells. The second step was to create differentially expressed genes. We ranked the mean counts of all genes of the second group of cells and chose 200 genes, starting with the one having a mean count just above s 1 . We selected another 200 genes beginning with the mean count just above s 2 = F C × s 1 . Then, we swapped the gene expression of these two equal numbers sets of selected genes in the second group of cells and got a simulation dataset with 2 cell groups with a known DE genes list. The parameter FC controlled the DE strength of DE genes between groups. We considered three DE strengths of DE genes: weak (FC = 1.5), moderate (FC = 2), and strong (FC = 2.5). In simulation 2, we chose PBMC68K (Zheng et al., 2017) as source data. PBMC68K consisted of transcription profiles of −68000 peripheral blood mononuclear cells and had 11 different cell subtypes with sample sizes ranging from −90 to −20000 (more details in Supplementary Material S1). We generated three simulation scenarios for each subtype of cells with three different levels (weak, moderate, and strong) of difference to test the sensitivity of detecting the DE genes.

All DE gene lists in simulation datasets were artificially set. We calculated all method performance indices according to known DE gene lists. Due to the unattainability of the whole accurate DE genes list of different cell groups in real single-cell data, we used different standards to set three potential DE gene lists and calculated all method performance indices.

Standard 1. Known DE genes from the literature.

Standard 2. The top 500 genes are ranked by the chosen number of times by all methods.

Standard 3. The top 1,000 genes are ranked by the chosen number of times by all methods.

For Standard 1, we collected dozens of known DE genes from various experiments based on bulk RNA-seq in the literature (Liu et al., 2001; Weng et al., 2012; Zhang et al., 2014; Darmanis et al., 2015). They were partial genes of the whole true DE genes between different cell clusters. For Standard 2 and 3, we ranked all genes’ chosen number of times by all methods and set the top 500 and 1,000 genes as potential DE genes between different cell clusters.

On the basis of the DE gene list in 2.3, we calculated a series of indices: F 1 scores, true positive rate (TPR, recall), false discovery rate (FDR), and time consumption to assess the performance of all methods. All indices were presented as the average value of 20 replications.

For this H 01 , we compared the positive expression ratio of the gene j in the two groups of cells. The total numbers of two groups of cells are n 1 , n 2 , respectively. And the numbers of positive expressions of the gene j in the two groups of cells are m j 1 = ∑ I ( y j 1 i ≠ 0 ) , m j 2 = ∑ I ( y j 2 i ≠ 0 ) , respectively. y j g i denotes the UMI count of the gene j of cell i in the group g = 1,2 . p j g is the gene j ’s positive expression proportion in the group g . p ^ j g is the estimator of p j g . Hence, the positive expression ratios of the gene j in group 1 and group 2 are p ^ j 1 = m j 1 n 1 , p ^ j 2 = m j 2 n 2 , respectively. H 01 : p j 1 = p j 2 ; H A 1 : p j 1 ≠ p j 2 z = p ^ j 1 − p ^ j 2 p ^ ∗ ( 1 − p ^ ∗ ) ( 1 n 1 + 1 n 2 ) ∼ N ( 0,1 ) w h e r e , p ^ ∗ = m j 1 + m j 2 n 1 + n 2 L 1 = 2 P ( Z > | z |   |   H 01   i s   t r u e )

In terms of hypotheses H 02 and H 03 , only the “On” state of each gene is involved in calculations. For hypothesis H 02 , we used the Student’s t-test to determine whether the two groups differ significantly on the central location of gene j ’s expression of the “On” state. H 02 : μ j 1 = μ j 2 ; H A 2 : μ j 1 ≠ μ j 2 t = x ¯ j 1 − x ¯ j 2 s j 1 2 m j 1 + s j 2 2 m j 2 ∼ t ( d f t ) L 2 = 2 P ( T > | t |   |   H 02   i s   t r u e )

μ j g is the mean of the gene j in the group g on the positive part (“on” state). x ¯ j g = ∑ x j g i m j g is the estimator of the μ j g . s j g is the estimator of the σ j g 2 , which is the variance of the gene j in the group g ′ s ‘on’ part.Where, x j g = { y j g i , w h i c h   y j g i > 0 }

For this H 03 , we used the Brown–Forsythe test to test the equality of scattering of gene j ’s positive expression. H 03 :   σ j 1 2 = σ j 2 2 ;   H A 3 :   σ j 1 2 ≠ σ j 2 2 W = ∑ g = 1 G ( m j g − 1 ) G − 1 ∑ g = 1 G m j g ( z ¯ j g − z ¯ j ) 2 ∑ g = 1 G ∑ i = 1 m j g ( z j g i − z ¯ j g ) 2 ∼ F ( G − 1 ,   ∑ g = 1 G ( m j g − 1 ) ) ,   z j g i = | x j g i − x ∼ j g | L 3 = P ( F ( 1 , m j 1 + m j 2 − 2 ) > W   | H 03   i s   t r u e ) Where, x ∼ j g in z j g i = | x j g i − x ∼ j g | is the median of the g-th subgroup. Then we performed statistical tests on each null hypothesis H 0 i , respectively. The p-value of each test is recorded as L i . We obtained a new statistic Q by combining three individual p-values L i of the statistics for each null hypothesis H 0 i . Q = − 2 ∑ i = 1 3 log ⁡ L i

Q follows the χ 2 distribution. If L i is independent, Q ∼ χ 2 ( 6 ) . The degree of freedom of Q is not equal to 6 in most scenarios because of the correlation of L i . To solve this problem, we obtained the freedom which is close to the real data distribution by sup d f   L ( d f | Q ) (more details inSupplementary Material S1).