Overview of TDEseqStatistical modeling

TDEseq is a temporal gene expression analysis approach that is primarily built upon the linear additive mixed models (LAMM) [39] framework to characterize the temporal gene expression changes for time-course (or longitudinal) scRNA-seq datasets (the “Materials and methods” section; Supplementary Text). Typically, we aim to detect one of four possible temporal gene expression patterns (i.e., growth, recession, peak, or trough) over multiple time points using both I-splines [35] and C-splines [36] basis functions (Fig. 1A) and examine one gene at a time. Briefly, in LAMM, we assume the log-normalized gene expression level of raw counts (the “Materials and methods” section), i.e., \(_\left(t\right)\) for gene \(g\), individual \(j\) and cell \(i\) at time point \(t\) is,

$$y_\left(t\right)=\boldsymbol^__g+\sum\nolimits_^Ks_k\left(t\right)\beta_+u_+e_,i=1,2,\cdots,n_j;g=1,2,\cdots,G;t=1,2,\cdots,T,j=1,\dots,M.$$

Fig. 1

Schematic overview of TDEseq and the methods comparison in simulations. A TDEseq is designed to perform temporal expression gene analysis of time-course scRNA-seq data. With a given gene, TDEseq determines one of four temporal expression patterns, i.e., growth, recession, peak, and trough. TDEseq combines the four p-values using the Cauchy combination rule as a final p-value, facilitating the detection of temporal gene expression patterns. B The quantile–quantile (QQ) plot shows the type I error control under the baseline parameter settings. The well-calibrated p-values will be expected laid on the diagonal line. The p-values generated from Mixed TDEseq (plum) and DESeq2 (brown) are reasonably well-calibrated, while Linear TDEseq (orange), tradeSeq (green), ImpulseDE2 (blue), Wilcoxon test (yellow) and edgeR (dark green) produced the p-values that are not well-calibrated. C The average power of 10 simulation replicates for temporal expression gene detection across a range of FDR cutoffs under the baseline parameter settings. Both versions of TDEseq exhibit high detection power of temporal expression genes, followed by DESeq2, edgeR, tradeSeq, and ImpulseDE2. Wilcoxon test does not fare well, presumably due to bias towards highly expressed genes. The TDEseq methods were highlighted using solid lines, while other methods were represented by dashed lines in the plots. D The comparison of Linear TDEseq, Mixed TDEseq, and ImpuseDE2 in terms of the accuracy of temporal expression pattern detection under the baseline parameter settings, at an FDR of 5%. The temporal expression genes detected by TDEseq demonstrated a higher accuracy than those detected by ImpluseDE2. E The quantile–quantile (QQ) plot shows the type I error control under the large batch effect parameter settings. The p-values generated from Mixed TDEseq coupled with scMerge (purple) and DESeq2 (brown) are reasonably well-calibrated, while Linear TDEseq (orange), Mixed TDEseq (plum), tradeSeq (green), ImpulseDE2 (blue), Wilcoxon test (yellow), and edgeR (dark green) generated the inflated p-values. F The average power of 10 simulation replicates the comparison of temporal expression gene detection across a range of FDR cutoffs under the large batch effect parameter settings. G The comparison of Linear TDEseq, Mixed TDEseq, and Mixed TDEseq coupled with scMerge and ImpuseDE2 in terms of the accuracy of temporal expression pattern detection under the large batch effect parameter settings, at an FDR of 5%. Since DESeq2, edgeR, tradeSeq, and Wilcoxon tests were not originally designed for pattern-specific detection we excluded them in the comparison. FDR denotes the false discovery rate

Where \(}}_\) is the cell-level or time-level covariate (e.g., cell size, or sequencing read depth), \(}_\) is its corresponding coefficient; \(}}_\) is a random vector to account for the variations from heterogeneous samples, i.e.,

$$}}_\sim MVN\left(0,_^}}_\right).$$

where \(}}_\) is a block diagonal matrix with a total of \(M\) block matrices, in which all elements of \(}}__\times _}\) are ones; \(_\) is the number of cells for the individual or replicate \(j\), and \(\sum_^_=N\); \(_\) is a random effect, which is an independent and identically distributed variable that follows a normal distribution with mean zero and variance \(_^\) to account for independent noise, i.e.,

$$_g=\left(e_,\cdots,e_\right)^\sim MVN(0,\sigma_g^2_).$$

Particularly, the variable \(_\left(t\right)\) is a smoothing spline basis function, which involves either I-splines or C-splines to model monotone patterns (i.e., growth and recession) and quadratic patterns (i.e., peak and trough), respectively [40] (Fig. 1A). The I-splines are defined as \(_^\left(x\right)=__}^_^\left(v\right)dv\) [35], while C-splines are defined as \(_^\left(x\right)=__}^_^\left(v\right)d(v)\) [36] based on I-splines, and \(_\) is its corresponding coefficient that is restricted to \(_\ge 0\), and \(l (l=\mathrm,\cdots ,L)\) is the number of grid points. We set \(L\) to be a total number of grid points, which is equal to the number of time points in scRNA-seq studies; \(k\) denotes the order of the spline function; \(MVN\) denotes the multivariate normal distribution.

Hypothesis testing

In the LAMM model mentioned above, we are interested in examining whether a gene shows one of four temporal expression patterns, i.e., growth, recession, peak, or trough (Fig. 1A). Testing whether a gene expression displays temporal gene expression patterns can be translated into testing the null hypothesis \(}}_:}}_=0\). Parameter estimates and hypothesis testing in LAMM are notoriously difficult, as the LAMM likelihood involves M-splines [35] (non-linear) subject to nonnegative constraints that cannot be solved analytically. To make the LAMM model scalable estimation and inference, we developed an approximate inference algorithm based on a cone programming projection algorithm [41, 42]. With parameter estimates, we computed a p-value for each of the four patterns using the test statistics [43], which follow a mixture of beta distributions [44]. Afterward, we combined these four p-values through the Cauchy combination rule [45]. The Cauchy combination rule allows us to combine multiple potentially correlated p-values into a single p-value to determine whether a gene exhibits the temporal expression pattern or not (the “Materials and methods” section; Additional file 1: Supplementary Text).

We refer to the above method as the mixed version of TDEseq (we denoted as Mixed TDEseq). Besides the mixed version, we have also developed a linear version of TDEseq (to distinguish Mixed TDEseq, we denoted this as Linear TDEseq) for modeling the small or no sample heterogeneity inherited in time-course scRNA-seq data (Additional file 1: Supplementary Text). Both versions of TDEseq were implemented in the same R package with multiple threads computing capability. The software TDEseq, together with the reproducibility analysis code, is freely available at https://sqsun.github.io/software.html.

TDEseq generates well-calibrated p-values and exhibits powerful gene detection of temporal expression changes in simulations

To benchmark the robustness and performance of TDEseq, we simulated extensive scRNA-seq datasets using the Splatter R package [46] and compared two versions of TDEseq with other five existing approaches but not specific designs for time-course scRNA-seq data analysis, which are the two-sided Wilcoxon rank-sum test (Wilcoxon test), tradeSeq [47], ImpulseDE2 [25], edgeR [27, 48], and DESeq2 [26] (the “Materials and methods” section). The simulations were typically designed to assess the ability of TDEseq in terms of type I error control and temporal gene detection power with varying various parameter settings, including the number of time points (i.e., 4, 5, or 6), the number of cells for each sample in each time point (i.e., 50, 100, or 200; three replicates/samples for each time point), the expected UMI counts for each cell of scRNA-seq data (i.e., 7.0 as low, 9.7 as medium, and 13.8 as high), the effect size of temporal expression changes (i.e., 0.1 as low, 0.4 as medium, and 0.7 as high), and the sample-level unwanted technical variations (i.e., batch effects; 0 as no batch effects, 0.04 as medium, and 0.12 as high).

To do so, we considered a baseline simulation scenario: the number of time points as 5; the number of cells in each sample as 100; the expected UMI counts for each cell as 9.7; the batch effect size as 0.04; the time point-specific effect size as 0.4; all cells were measured by 10,000 genes, in which 1,000 genes were randomly assigned one of four possible temporal patterns (i.e., growth, recession, peak, and trough; Additional file 2: Fig. S1A-S1D) in power simulations. With the baseline parameter settings, we varied one parameter at a time to examine whether the gene was temporally expressed over multiple time points. Notably, the expected UMI counts under baseline settings were estimated from the lung adenocarcinoma progression scRNA-seq data [14] (the “Materials and methods” section).

With the baseline parameter setup, we found that only Mixed TDEseq and DESeq2 generated the well-calibrated p-values under the null simulations, whereas all other methods produced the inflated or conserved p-values (Fig. 1B). Besides, for the power simulations, Linear TDEseq and Mixed TDEseq can produce a more powerful temporal expression pattern detection rate across a range of FDR cutoffs (Fig. 1C). Specifically, with a false discovery rate (FDR) of 5%, the power detection rate of both Linear and Mixed TDEseq was 43.3% and 40.8%, respectively, followed by DESeq2 was 38.7%, edgeR was 36.4%, tradeSeq was 25.9%, and ImpulseDE2 was 18.4%. Furthermore, we also examined the accuracy of pattern detection, finding both Linear and Mixed versions of TDEseq outperformed ImpluseDE2 (the sole method capable of identifying pattern-specific genes). Specifically, with an FDR of 5%, the averaged accuracy of pattern examination (with 10-time repeats) for Mixed TDEseq achieved 99.0% for growth, 100% for recession, 80.4% for peak, and 43.6% for trough. In contrast, ImpulseDE2 achieved 73.1% for growth, 83.5% for recession, 39.3% for peak, and 21.3% for trough (Fig. 1D).

In addition, we systematically examined the performance of the type I error control rate under other parameter settings. Our findings indicate that Mixed TDEseq consistently produces well-calibrated p-values (Additional file 2: Fig. S2A, S2B, S3A, S4A, and S4B) except when dealing with high UMI counts (Additional file 2: Fig. S3B). These observations are presumably due to the presence of sample-level variations or batch effects associated with high UMI counts. On the other hand, in terms of temporal expression gene detection and averaged accuracy of pattern examination, either Mixed or Linear TDEseq displayed more powerful performance across a range of parameter settings regardless of the number of time points (Additional file 2: Fig. S2C and S2D), the low expected UMI counts for each cell (Additional file 2: Fig. S3C), the large number of cells per sample (Additional file 2: Fig. S4D), and the small effect size setups (Additional file 2: Fig. S5), as well as the accuracy of pattern examination (Additional file 2: Fig. S6). Meanwhile, we found both pseudobulk-based methods, i.e., either edgeR or DESeq2, performed well with high UMI counts (Additional file 2: Fig. S2D) and small number of cells per sample (Additional file 2: Fig. S3C). These observations were probably consistent with the previous studies that DESeq2/edgeR performed well on log-normal distributed small sample size RNA-seq data [49, 50]. Taken together, we summarized the findings on detection power at an FDR of 5% across diverse parameter settings. The results demonstrated that the power of temporal gene detection increases with a rise in the number of time points, effect size, and UMI counts. Conversely, it diminishes as the number of cells within each individual increase, along with sample-level variations (i.e., batch effects) (Additional file 2: Fig. S7). Notably, the Wilcoxon test did not fare well in all power simulations, presumably due to failure to properly control the type I error rate.

In addition, we further examined the performance of TDEseq in other two temporal expression patterns: (1) a plateau in the first few time points, then another plateau in the last few time points (we referred to this pattern as a bi-plateau pattern; Additional file 2: Fig. S1E), and (2) a multi-mode pattern at begin time points then stable in last time points (we referred this pattern as a multi-modal pattern; Additional file 2: Fig. S1F). Under the bi-plateau pattern, Mixed TDEseq still displayed more powerful performance than other methods (Additional file 2: Fig. S8A), suggesting the shape-restricted spline function is flexible to capture bi-plateau patterns. In contrast, under the multi-modal pattern, all methods achieved low detection power, but edgeR and DESeq2 showed a higher performance than other methods (Additional file 2: Fig. S8B). However, this may not be a great issue since the multi-modal pattern may be a rare scenario in real data applications [25].

TDEseq coupled with batch removal strategy exhibits excellent performance in analyzing large heterogeneous scRNA-seq data

Intuitively, in situations with minimal or no sample-level variations (i.e., batch effects), it is reasonable to expect that trajectory-based differential expression methods (e.g., tradeSeq) would yield comparable results to temporal-based differential expression methods. To do so, we reduced the batch effect size to zero. As a result, we observed that both versions of TDEseq and tradeSeq generated well-calibrated p-values under the null simulations, whereas ImpulseDE2 demonstrated overly conservative p-values and the Wilcoxon test displays inflated p-values (Additional file 2: Fig. S9A). Again, both versions of TDEseq and all other approaches generated comparable results of temporal expression pattern (Additional file 2: Fig. S9B). As we know, the presence of unwanted batch effects poses substantial obstacles in detecting temporal expression changes. We therefore increased the batch effect size to 0.12. As a result, we found Mixed TDEseq outperformed other methods in terms of temporal expression pattern detection power (Fig. 1F). However, the p-values generated by Mixed TDEseq were not well-calibrated (Fig. 1E).

To this end, to properly control the unwanted variables in the large batch effects scenario, we additionally carried out the batch effects correction procedure prior to performing temporal gene expression analysis. To do so, we benchmarked five existing batch removal methods that can return the corrected gene expression matrix, including MNN [51], scMerge [52], ZINB-WaVE [53], ComBat [54], and Limma [55]. As a result, with evaluation criterion iLISI score [56] (the “Materials and methods” section) for batch correction approaches, we found scMerge (0.49) achieved a higher alignment score than Limma (0.48), ComBat (0.48), MNN (0.11), and ZINB-WaVE (0.21; Additional file 2: Fig. S10). Moreover, we found Mixed TDEseq coupled with scMerge (Mixed TDEseq + scMerge) performed reasonably well in terms of the type I error control (Fig. 1E and Additional file 2: Fig. S11A) and was more powerful in detecting temporal expression genes (Fig. 1F and Additional file 2: Fig. S11B), suggesting this combination is suitable for time-course scRNA-seq data with strong sample-level variations (i.e., batch effects). Taken together, TDEseq coupled with scMerge may be an ideal combination for the identification of temporal gene expression patterns when time-course scRNA-seq data involves large heterogeneous samples.

TDEseq performs well in the intertwined cells among time points

The simulations above all display the time point-specific expression. To mimic the cell differentiation scenario where the same type of cells were intertwined among time points (Additional file 2: Fig. S12A), we simulated additional scRNA-seq datasets (denoted as smudged data) using the Symsim R package [57] (the “Materials and methods” section). Consequently, the pseudotime was inferred using Slingshot [6] according to the recommendation from the previous studies [31, 32]. In this simulation, we first took the inferred pseudotime as inputs in tradeSeq and ImpulseDE2, while the time points as inputs in both versions of TDEseq, edgeR, and DESeq2. As a result, we observed that the performance of Linear TDEseq was comparable with ImpluseDE2 in a small proportion of intertwined cells between time points (Additional file 2: Fig. S12B). With a medium proportion of intertwined cells (Additional file 2: Fig. S12C) and a large proportion of intertwined cells (Additional file 2: Fig. S12D), the pseudotime-based methods tradeSeq and ImpulseDE2 outperformed the time points-based methods, both versions of TDEseq, edgeR, and DESeq2. Furthermore, we took the inferred pseudotime as inputs in both versions of TDEseq, Linear TDEseq outperformed Mixed TDEseq, and ImpulseDE2, but not tradeSeq (Additional file 2: Fig. S12E).

In addition, we further examined the temporal expression patterns that were detected by Linear TDEseq. As a result, we found even though the pseudotime as inputs, TDEseq displayed four distinct temporal expression patterns (Additional file 2: Fig. S12F). Therefore, TDEseq was also useful for detecting temporal expression patterns with the pseudotime as inputs.

TDEseq detects drug-associated temporal expression changes of time-course scRNA-seq data

We first applied TDEseq on a drug-treatment time-resolved scRNA-seq dataset (Additional file 3: Table S1). The data were assayed by Well-TEMP-seq protocols to profile the transcriptional dynamics of colorectal cancer cells exposed to 5-AZA-CdR [23] (the “Materials and methods” section). This scRNA-seq dataset consists of D0, D1, D2, and D3 four time points (Fig. 2A), and each time point contains 4,000 cells. We expected these scRNA-seq datasets to exhibit minimal individual heterogeneity across multiple time points since Well-TEMP-seq addressed the cell lines within one chip [23] (Additional file 2: Fig. S13A). Therefore, we performed the temporal gene detection methods without batch effects correction. Since only one sample was involved in each time point, we excluded Mixed TDEseq, edgeR, and DESeq2 in this application.

Fig. 2

The time-resolved scRNA-seq data analysis for the HCT116 cell lines after 5-AZA-CdR treatment. A The experimental design of HCT116 cell lines treated with 5-AZA-CdR. The scRNA-seq data were assayed by Well-TEMP-seq protocols, consisting of four time points, i.e., D0, D1, D2, or D3 after treatment. B The quantile–quantile (QQ) plot shows the type I error control under the null simulations with permutation strategy. The well-calibrated p-values will be expected laid on the diagonal line. The p-values produced by Linear TDEseq (orange) and tradeSeq (green) are reasonably well-calibrated, while those from ImpulseDE2 (blue) are overly conservative. C The power comparison of temporal expression gene detection across a range of FDR cutoffs. Linear TDEseq was highlighted using solid lines, while other methods were represented by dashed lines in the plots. Linear TDEseq displays the powerful performance of temporal expression gene detection. D The heatmap demonstrates the pattern-specific temporal expression genes that were identified by Linear TDEseq. Gene expression levels were log-transformed and were standardized using z-scores for visualization. The top-ranked temporal expression genes identified by Linear TDEseq show distinct four patterns. E The Venn diagram shows the overlapping of the temporally expressed genes (FDR ≤ 0.05) identified by Linear TDEseq, tradeSeq, or ImpulseDE2. Those method-specific unique genes were enriched in the number of GO terms (NGO, BH-adjusted p-value < 0.05). The temporal expression genes detected by Linear TDEseq were enriched with a greater number of GO terms. F The UMAP shows two temporal expression genes, i.e., DKK1 and IFITM3, which were identified by Linear TDEseq but not by tradeSeq. G The bubble plot demonstrates the significant GO terms enriched by pattern-specific temporal expression genes, which were identified by Linear TDEseq. The peak-specific temporal expression genes enriched more significant GO terms. The Wilcoxon test was excluded from this comparison due to its poor performance in simulations. DESeq2 and edgeR were excluded from this comparison due to only one sample at each time point. FDR denotes the false discovery rate

Next, we examined the ability of Linear TDEseq in terms of type I error control. To do so, we utilized a permutation strategy (repeated 5 times) to construct a null distribution (the “Materials and methods” section). Consistent with simulation studies, we found Linear TDEseq can produce well-calibrated p-values while tradeSeq produced inflated p-values and ImpulseDE2 produced overly conserved p-values (Fig. 2B). Besides, in terms of temporally expressed gene detection, Linear TDEseq outperformed other methods across a range of FDR cutoffs (Fig. 2C), even in pattern-specific temporal expression gene detection (Additional file 2: Fig. S13B). For example, Linear TDEseq identified a total of 5,596 temporally expressed genes at an FDR of 5%, including 1,341 growth genes, 1,177 recession genes, 225 trough genes, and 2,853 peak genes, which displayed four distinct temporal expression patterns (Fig. 2D). In contrast, ImpulseDE2 identified a total of 4,792 temporally expressed genes and tradeSeq detected a total of 2,672 temporally expressed genes (Additional file 3: Table S3). Overall, besides the 2,427 common shared temporally expressed genes detected by Linear TDEseq, tradeSeq, and ImpulseDE2 methods (Fig. 2E), a total of 559 temporally expressed genes were uniquely detected by TDEseq, which were also significantly enriched in cell cycle DNA replication (GO:0044786; BH adjusted p-value = 7.89e − e) and response to interleukin1 (GO:0070555; BH adjusted p-value = 0.031). In contrast, tradeSeq or ImpulseDE2 unique genes were not enriched in 5-AZA-CdR treatment response associated GO terms (Fig. 2E). Specifically, we found tumor suppressor genes which were a target of 5-AZA-CdR, i.e., DKK1 [23, 58] was identified by TDEseq as top-ranked significant temporal expression genes (p-value < 1e − 300, FDR = 0), while not being detected by tradeSeq (p-value = 0.82, FDR = 0.90), probably due to though this gene had clearly peak pattern, the log fold change was small enough, and it was difficult to detect with penalized splines; besides, a 5-AZA-CdR response gene IFITM3 [58, 59] was also identified by TDEseq as top-ranked significant genes (p-value < 1e − 300, FDR = 0), but not detected by tradeSeq (p-value = 0.19, FDR = 0.36, Fig. 2F).

Finally, we performed gene set enrichment analysis (GSEA) on the pattern-specific temporal expression genes to examine top GO terms enriched by the given gene lists (the “Materials and methods” section). Specifically, with an FDR of 5%, a total of 1,341 growth-specific temporal expression genes were detected by Linear TDEseq. These genes were enriched in a total of 179 GO terms. Because the 5-AZA-CdR treatment leads HCT116 cells to a viral mimicry state, and triggers the antiviral response [60], we expected a result of an immune response that drives the immune-associated genes is upregulated. Indeed, the GO terms contain many immune response terms, such as the cell activation involved in the immune response process (GO: 0002263; BH-adjusted p-value = 8.90e − 7), indicating immune response was activated by the 5-AZA-CdR treatment; a total of 1,177 recession-specific temporal expression genes were enriched in a total of 244 GO terms, e.g., many regulations of histone methylation terms such as positive regulation of histone H3-K4 methylation (GO: 0051571; BH-adjusted p-value = 0.047), implying DNA methylation inhibitions and gene expression regulation were occurred after 5-AZA-CdR treatment, due to the global DNA demethylation effects of 5-AZA-CdR [61]; a total of 2,853 peak-specific temporal expression genes were enriched in a total of 249 GO terms. For example, the ATP metabolic process pathways, particularly oxidative phosphorylation (GO: 0006119; BH-adjusted p-value = 3.30e − 7), are impacted by the increase of intracellular ROS and mitochondrial superoxide induced by 5-AZA-CdR. However, this effect diminishes over time [62]; similarly, a total of 225 trough-specific temporal expression genes were enriched in a total of 89 GO terms, with a significant portion belonging to cell cycle pathways, including mitotic nuclear division (GO: 0140014; BH-adjusted p-value = 1.06e − 6). These findings suggest that 5-AZA-CdR treatment may lead to the suppression of tumor cell proliferation and division [60] (Fig. 2G).

TDEseq detects hepatic cell differentiation-associated temporal expression genes of time-course scRNA-seq data

We next applied TDEseq to a hepatoblast-to-hepatocyte transition study from the C57BL/6 and C3H embryo mice livers [38] (the “Materials and methods” section; Additional file 3: Table S1). This scRNA-seq dataset consists of 7 developmental stages from 13 samples, including E10.5 (54 cells from 1 sample), E11.5 (70 cells from 2 samples), E12.5 (41 cells from 2 samples), E13.5 (65 cells from 2 samples), E14.5 (70 cells from 2 samples), E15.5 (77 cells from 2 samples), and E17.5 (70 cells from 2 samples) [63] (Fig. 3A). Compared with the above time-resolved scRNA-seq data, this time-course scRNA-seq dataset contains multiple samples at each stage, exhibiting small individual heterogeneity across all developmental stages (Additional file 2: Fig. S14A). Therefore, we carried out both versions of TDEseq that would be expected to be comparable in such a scenario and excluded edgeR and DESeq2 in this application due to one or two samples involved in each time point.

Fig. 3

The time-course scRNA-seq data analysis for mouse fetal liver development. A The experimental design of mouse fetal liver sample collection. The scRNA-seq data were assayed on the FACS isolated cell populations, consisting of seven liver developmental stages, i.e., E10.5, E11.5, E12.5, E13.5, E14.5, E15.5, and E17.5. B The quantile–quantile (QQ) plot shows the type I error control under the permutation strategy. The well-calibrated p-values will be expected laid on the diagonal line. The p-values produced by Linear TDEseq (orange), Mixed TDEseq (plum), and tradeSeq (green) are reasonably well-calibrated, while those from ImpulseDE2 (blue) are overly conservative. C The power comparison of temporal expression gene detection across a range of FDR cutoffs. The TDEseq methods were highlighted using solid lines, while other methods were represented by dashed lines in the plots. Both versions of TDEseq display the powerful performance of temporal expression gene detection. D The heatmap demonstrates the pattern-specific temporal expression genes that were identified by Linear TDEseq. Gene expression levels were log-transformed and were standardized using z-scores for visualization. The top-ranked temporal expression genes identified by Linear TDEseq show distinct four patterns. E The Venn diagram shows the overlapping of the temporally expressed genes (FDR ≤ 0.05) identified by Linear TDEseq, tradeSeq, or ImpulseDE2. Those method-specific unique genes were enriched in the number of GO terms (NGO, BH-adjusted p-value < 0.05). The temporal expression genes detected by Linear TDEseq were enriched more GO terms. F The UMAP shows two temporal expression genes, i.e., Atf4 and Itgb1, which were uniquely identified by Linear TDEseq. G The bubble plot demonstrates the significant GO terms enriched by pattern-specific temporal expression genes, which were identified by Linear TDEseq. The recession-specific temporal expression genes enriched more significant GO terms, whereas trough-specific temporal expression genes were not enriched in any GO terms. The Wilcoxon test was excluded from this comparison due to its poor performance in simulations. DESeq2 and edgeR were excluded from this comparison due to only one or two samples at each time point. FDR denotes the false discovery rate

To do so, we first examined the ability of TDEseq in terms of type I error control using permutation strategies (the “Materials and methods” section). Consistent with the simulation results, Linear TDEseq, Mixed TDEseq, and tradeSeq could produce the well-calibrated p-values whereas ImpulseDE2 generated overly conservative p-values (Fig. 3B). Besides, in terms of temporally expressed gene dete