CQRM: A Composite Quantile Regression Model for Differential Expression Analysis of Non-Coding RNA Count Data

1. Introduction

Non-coding RNAs (ncRNAs) comprise a diverse class of RNA molecules that lack protein-coding potential but play essential roles in gene regulation and cellular function. These include long non-coding RNAs (lncRNAs), microRNAs (miRNAs), small nuclear RNAs (snRNAs) and circular RNAs (circRNAs), among others. Increasing evidence demonstrates that ncRNAs are involved in chromatin organization, transcriptional regulation, RNA processing, translation, and intercellular signaling, and contribute to numerous physiological and pathological processes across plants and animals [1,2,3,4,5]. Dysregulation of ncRNA expression has been implicated in a wide range of diseases, including cancer, neurodegenerative disorders, cardiovascular diseases, and psychiatric conditions, making ncRNAs attractive biomarkers and therapeutic targets [1,2,3,4,5].

Compared with coding genes, many ncRNAs exhibit lower expression levels, stronger tissue- and cell-type specificity, reduced evolutionary conservation, and greater expression heterogeneity [6,7]. These features suggest specialized regulatory functions but also create substantial analytical challenges for transcriptomic studies. The widespread adoption of bulk and single-cell RNA sequencing (scRNA-seq) technologies in recent decades has enabled transcriptome-wide profiling of ncRNA expression and stimulated growing interest in identifying differentially expressed (DE) ncRNAs associated with biological processes and disease phenotypes. However, DE analysis of ncRNAs is considerably more challenging than that of protein-coding genes, as many ncRNA species are characterized by low abundance, high sparsity, substantial biological variability, and multimodal expression distributions [8,9]. These characteristics are often biologically meaningful and may reflect heterogeneous transcriptional regulation, cell-state specificity, or context-dependent expression patterns. Consequently, DE analysis may manifest not only shifts in mean expression levels but also changes in variability, skewness, sparsity, or other distributional properties. Representative examples of such ncRNA expression distributions are shown in Figure S1.

For count data from RNA-seq or scRNA-seq, existing DE methods can be broadly classified into parametric count-based methods, nonparametric methods, and zero-inflation-adjusted methods (Table 1). Parametric methods such as edgeR [10] and DESeq2 [11] model gene expression counts using the negative binomial (NB) distribution and primarily focus on detecting differences in mean expression levels between conditions. Limma-voom [12] transforms count data and employs empirical Bayes moderated linear models to improve variance estimation and inference, particularly in small-sample settings. Although these methods have demonstrated excellent performance for mRNA analyses, the underlying assumptions may be violated for many ncRNAs whose expression distributions exhibit excessive zeros, heavy tails, or multimodality.

Nonparametric DE methods such as the Wilcoxon rank-sum test [13], SAMSeq [14], and NOISeq [15] relax distributional assumptions and often provide improved robustness in heterogeneous settings. The Wilcoxon test has been highlighted in recent comparative studies as a more robust approach for DE analysis with better false discovery rate (FDR) control and power performance [16]. However, rank-based methods generally focus on differences in central tendency and may have limited sensitivity when DE arises through more complex distributional changes in variability, skewness, sparsity, or multimodality. Previous benchmarking studies have reported that many conventional RNA-seq DE pipelines show reduced sensitivity and unstable FDR control for low-abundance transcripts, including lncRNAs and other ncRNAs [9]. For scRNA-seq data, several zero-inflated methods have been proposed, including zero-inflated negative binomial (ZINB) models [17] and a zero-inflated quantile regression model [18]. While these methods improve modeling of sparse expression data, they primarily focus on dropout mechanisms. Dedicated DE methods specifically designed for ncRNAs remain relatively limited. Existing methods such as lncDIFF [19] target count-transformed continuous data and have demonstrated improved performance for lowly expressed lncRNAs, but generally remain tied to specific distributional assumptions or particular ncRNA classes. Given the growing recognition of ncRNAs as critical regulators in health and disease, there remains a pressing need for flexible statistical frameworks capable of accommodating the diverse and heterogeneous expression patterns observed across ncRNA species.

In this paper, we propose a Composite Quantile Regression Model (CQRM) for DE analysis of ncRNA sequencing data. Unlike conventional DE methods that primarily target differences in mean expression, CQRM integrates information across multiple quantiles of the expression distribution, enabling the detection of a wider range of DE patterns. By avoiding restrictive parametric assumptions, CQRM remains robust to the heterogeneous, sparse, and often non-standard expression distributions commonly observed in ncRNAs. In addition, its regression-based framework naturally accommodates complex experimental designs and covariate adjustment. Although quantile regression has previously been applied to DE analysis, such as in ZIAQ for single-cell RNA-seq data [18], existing approaches typically evaluate differential expression at individual quantiles. In contrast, CQRM aggregates evidence across multiple combinations of quantiles to capture global distributional differences between conditions, providing a flexible and robust framework for DE analysis across diverse non-coding RNA species.

2. Method

We propose a composite quantile regression-based framework for differential expression (DE) analysis of ncRNA count data from RNA-seq or scRNA-seq experiments. The method first models the association between expression counts and covariates at individual quantiles using quantile regression with jittering to accommodate the discrete nature of count data. It then aggregates evidence across multiple quantiles through a composite testing procedure, enabling the detection of differential expression patterns that may occur throughout the expression distribution rather than at a single quantile or solely in the mean.

2.1. Modeling ncRNA Count Data Using Quantile Regression with Jittering

Denote by Y the $n\times G$ matrix of count data from RNA-seq or scRNA-seq where $Y_{ig}$ indicates the count value for gth ncRNA ($1\le g\le G$) in sample or cell i ($1\le i\le n$). Denote by X the $n\times K$ matrix of covariates. In a typical DE analysis, we are interested in comparing the RNA expression profiles between a case and a control group. Without loss of generality, we assume one covariate ($K=1$), where $X_i=1$ when a sample is from the case group (e.g. cancer patients) and $X_i=0$ when a sample is from the control group (e.g. normal controls). Standard quantile regression (QR) assumes a continuous response variable. Directly applying QR to count data is challenging because the discreteness of the response leads to non-smoothness and non-standard asymptotic behavior of the conditional quantile estimator [20]. To address this issue, Machado and Santos Silva [21] proposed adding independent uniform noise to the count outcome, a procedure commonly referred to as jittering [22]. Specifically, let $U_{ig}$ denote the uniform random variable independent of $Y_{ig}$, and define the jittered response $Z_{ig}=Y_{ig}+U_{ig}$. The resulting variable $Z_{ig}$ is continuous, and its conditional quantiles have a one-to-one relationship with those of the original count variable under the jittering construction. After applying a monotonic transformation to the jittered count variable $T\left(Z_{ig}\right)={log}_2(Z_{ig}-\tau )$ if $Z_{ig}>\tau$ and ${log}_2\left(\zeta \right)$ otherwise, where $\zeta$ is a suitably small positive number, the conditional quantile function of the transformed response can be modeled using quantile regression as:

$$Q_{T\left(Z_{ig}\right)}\left(\tau \right|{\mathbf{x}}_i)={\alpha }_g\left(\tau \right)+{\mathbf{x}}_i^T{\beta }_g\left(\tau \right),$$

(1)

where $Q_{T\left(Z_{ig}\right)}\left(\tau \right|{\mathbf{x}}_i)$ denotes the conditional $\tau$-th quantile of $T\left(Z_{ig}\right)$ given covariates ${\mathbf{x}}_i$, ${\alpha }_g\left(\tau \right)$ is the quantile-specific intercept, and ${\beta }_g\left(\tau \right)$ is the vector of quantile-specific regression coefficients. The transformation $T\left(Z_{ig}\right)$ we employed here is the same as in Machado and Santos Silva [21], which works as a reasonable approximation in many situations [21].

The coefficient estimates $(\widehat{\alpha }\left(\tau \right),\widehat{\beta }\left(\tau \right))$ for any fixed $\tau$ can be obtained by

$$({\widehat{\alpha }}_g\left(\tau \right),{\widehat{\beta }}_g\left(\tau \right))=\underset{\alpha \left(\tau \right),\beta \left(\tau \right)}{\mathrm{argmin}}{\displaystyle \sum _{i=1}^n}{\rho }_{\tau }(T\left(Z_{ig}\right)-{\alpha }_g-{\mathbf{x}}_i^T{\beta }_g\left(\tau \right)),$$

(2)

where ${\rho }_{\tau }\left(u\right)=u\{\tau -I(u<0)\}$.

To decrease the estimate’s dependence on the random draw of $U_{ig}$, we averaged the QR estimates from M repeatedly and independently “jittered” samples and used the “average-jittering” estimator ${\widehat{\beta }}_g^A\left(\tau \right)$ as our final estimate:

$${\widehat{\beta }}_g^A\left(\tau \right)={\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M}{\widehat{\beta }}_g^{\left(m\right)}\left(\tau \right),$$

(3)

where ${\widehat{\beta }}_g^{\left(m\right)}\left(\tau \right)$ is the estimate derived from m-th jitter sample. The “average-jittering” estimator is more efficient than the estimate obtained with a single sample. Following Machado and Santos Silva [21], we also showed the asymptotic behavior of ${\widehat{\beta }}_g^A\left(\tau \right)$ (see Appendix A).

Similarly, we derived $T^A\left(Z_{ig}\right)$ and ${\widehat{\alpha }}_g^A\left(\tau \right)$ as:

$$\begin{array}{cc}\hfill T^A\left(Z_{ig}\right)& ={\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M}T\left(Z_{ig}^{\left(m\right)}\right);\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\widehat{\alpha }}_g^A\left(\tau \right)& ={\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M}{\widehat{\alpha }}_g^{\left(m\right)}\left(\tau \right).\hfill \end{array}$$

(5)

2.2. Composite Quantile Regression for Differential Expression Testing

A key limitation of existing QR-based differential expression methods is that inference is typically conducted at a single quantile level. Such approaches are most sensitive to differential expression occurring at that particular portion of the distribution, but may lose power when expression differences are distributed across multiple quantiles or when the affected quantiles vary across ncRNA species. Given the heterogeneous and often multimodal expression distributions observed in ncRNA sequencing data, differential expression may manifest through changes in location, variability, sparsity, or other distributional characteristics. Motivated by this observation, we propose a composite QR framework that aggregates evidence across multiple quantiles to capture a broader spectrum of DE patterns.

For a fixed quantile level $\tau$, inference can be performed using a score test. Compared with asymptotic Wald-type inference, score tests avoid estimation of the conditional density function and generally exhibit improved finite-sample performance, particularly at extreme quantiles [23]. In contrast to bootstrap-based procedures, score tests are also computationally efficient, making them well suited for large-scale RNA-seq studies involving thousands of ncRNAs. For the quantile regression model introduced in Section 2.1, the hypothesis of interest at quantile $\tau$ (without loss of generality, we will drop subscript g from here on) is

$$H_0:\beta \left(\tau \right)=0\mathrm{vs}{H}_1:\beta \left(\tau \right)\ne 0.$$

(6)

The corresponding score statistic is

$$S\left(\tau \right)=n^{-1/2}{\displaystyle \sum _{i=1}^n}{\varphi }_{\tau }\{T^A\left(Z_i\right)-{\tilde{\alpha }}^A\left(\tau \right)\}{x}_i^{*},$$

(7)

where ${\varphi }_{\tau }\left(u\right)=\tau -I\{u<0\}$ is the quantile score function, ${\tilde{\alpha }}^A\left(\tau \right)$ denotes the estimated intercept under the null, and $x_i^{*}$ is the ith element of the residualized covariate vector $X^{*}={(x_1^{*},\dots ,x_n^{*})}^T$. Specifically, $X^{*}=(I-P_{\mathbf{1}})X$, where $P_{\mathbf{1}}$ is the projection matrix onto the space spanned by the intercept. Thus, under the one-covariate setting, $X^{*}$ is an $n\times 1$ vector representing the covariate after removing its projection onto the intercept space. $S\left(\tau \right)$ measures the association between the transformed expression $T^A\left(Z_i\right)$ and X at the $\tau$-th conditional quantile, i.e., after removing the collinearity between X and the intercept via $X^{*}$. Rather than testing each quantile separately, our goal is to assess differential expression across multiple quantiles simultaneously. Let ${\tau }_1,{\tau }_2,\dots ,{\tau }_K$ denote a prespecified set of quantile levels and let ${\tau }^{\left(w\right)}$ denote a specific non-empty subset of $\{{\tau }_1,{\tau }_2,\dots ,{\tau }_K\}$, with up to $2^K-1$ possible subsets. The choice of quantiles involves a trade-off between capturing stronger signals and introducing additional noise. In practice, a tuning procedure is recommended to determine an appropriate set of quantiles that balances statistical power and Type I error control. In this study, we use the 0.2, 0.4, 0.6, and 0.8 quantiles for all CQRM-based methods. For any non-empty subset w, let ${\tau }^{\left(w\right)}={({\tau }_{\left(1\right)},{\tau }_{\left(2\right)},\dots ,{\tau }_{\left(h\right)})}^T$ denote the corresponding quantile vector of size h, and define the vector of rank-score statistics as $S^{\left(w\right)}={(S\left({\tau }_{\left(1\right)}\right),S\left({\tau }_{\left(2\right)}\right),\dots ,S\left({\tau }_{\left(h\right)}\right))}^T$. We can then conduct a composite score test aggregating information across these quantiles. Under the null, $S^{\left(w\right)}$ follows an asymptotically multivariate normal distribution:

$$S^{\left(w\right)}\sim N(0,{\Sigma }^{\left(w\right)}),$$

(8)

where ${\Sigma }^{\left(w\right)}$ is the $h\times h$ covariance matrix whose $(j,j^{\prime })$-th element is

$${\sigma }_{j{j}^{\prime }}=n^{-1}\{min({\tau }_j,{\tau }_{j^{\prime }})-{\tau }_j{\tau }_{j^{\prime }}\}{X}^{*T}{X}^{*},j,j^{\prime }=1,\dots ,h.$$

(9)

In particular, when $j=j^{\prime }$, this reduces to the diagonal elements

$${\sigma }_{jj}=n^{-1}{\tau }_j(1-{\tau }_j)X^{*T}{X}^{*},j=1,\dots ,h,$$

(10)

whereas when $j\ne j^{\prime }$, the off-diagonal elements are

$${\sigma }_{j{j}^{\prime }}=n^{-1}\{min({\tau }_j,{\tau }_{j^{\prime }})-{\tau }_j{\tau }_{j^{\prime }}\}{X}^{*T}{X}^{*},j,j^{\prime }=1,\dots ,h,j\ne j^{\prime }.$$

(11)

The composite rank score test statistic is then constructed based on $S^{\left(w\right)}$:

$$R^{\left(w\right)}={S^{\left(w\right)}}^T{{\Sigma }^{\left(w\right)}}^{-1}{S}^{\left(w\right)},$$

(12)

which is shown to asymptotically follow a ${\chi }^2$ distribution with degrees of freedom equal to h under the null hypothesis. This statistic jointly evaluates differential expression across all quantiles contained in the combination and therefore captures a wider range of distributional changes than single-quantile tests.

To adaptively identify the most informative quantile combination for each ncRNA, we evaluate all $2^K-1$ non-empty quantile combinations and compute the corresponding p-values $p^{\left(w\right)}$. We then employ a minimum-p-value (minP) procedure [24] and treat the quantile combination with the minimum p-value as the one reflecting the significance of the coefficient in the hypothesis test for differential expression, i.e. $w_{best}=\underset{w}{\mathrm{argmin}}{p}^{\left(w\right)}$.

Since the minimum p-value MinP $p^{\left(w_{best}\right)}$ is obtained by searching across multiple quantile combinations, it does not follow a uniform distribution under the null hypothesis. A straightforward approach is to estimate its null distribution using permutation. However, such a procedure is computationally prohibitive for large-scale RNA-seq studies because the quantile regression model must be repeatedly refitted for every gene and every permutation.

To improve computational efficiency, we developed an approximate inference procedure (CQRM-Approx) that avoids gene-level permutations by directly approximating the null distribution of the MinP statistic. Under the null hypothesis, the composite score vector asymptotically follows a multivariate normal distribution $\mathbf{S}\sim N(\mathbf{0},\Sigma )$ (), where $\Sigma$ is estimated from the sample covariance matrix and is shared across genes. We first generate B independent realizations ${\mathbf{S}}^{\left(1\right)},\dots ,{\mathbf{S}}^{\left(B\right)}$ from this null distribution. For each realization, p-values are computed for all $2^K-1$ candidate quantile combinations using the corresponding covariance submatrices. Then we record the minimum p-value among all quantile combinations with the same degree of freedom as $w_{best}$, denoted by $p^{\left(b\right)}$. Repeating this procedure across all B realizations yields an empirical null distribution of the MinP statistic corresponding to the same degree of freedom as $w_{best}$.

For a given gene, the final adjusted p-value is computed as: $p_{final}={\displaystyle \frac{{\displaystyle \sum _{b=1}^B}I\{p^{\left(b\right)}<p^{\left(w_{best}\right)}\}}{B}}$, which represents the proportion of null MinP statistics that are at least as significant as the observed MinP statistic. By generating the null distribution directly from the asymptotic score distribution, CQRM-Approx eliminates the need for repeated permutation and model refitting while substantially reducing computational cost. Theoretical complexity analysis (Table S1) and empirical runtime benchmarking (Figure S2, Table S2) demonstrate that CQRM-Approx achieves considerable computational savings relative to the permutation-based implementation with minimal loss of accuracy.

Unlike conventional differential expression methods that primarily target mean shifts, the proposed CQRM method aggregates evidence across multiple quantiles and adaptively identifies the quantile combination that best characterizes differential expression. As a result, CQRM is capable of detecting a broader range of distributional changes, including shifts affecting only the lower or upper tails, heterogeneous expression variability, and multimodal expression patterns commonly observed in ncRNA sequencing studies. To investigate alternative patterns of differential expression, we also developed two related procedures: CQRM-Region and CQRM-Best. CQRM-Region restricts the candidate quantile combinations to contiguous quantile intervals, thereby favoring differential expression signals that extend smoothly across neighboring portions of the expression distribution. In contrast, CQRM-Best considers all individual quantiles and selects the most significant quantile-specific test. Next, we evaluate the performance of the proposed CQRM variants and compare them with existing state-of-the-art DE methods using semi-synthetic data simulations.

3. Results

3.1. Simulation

3.1.1. Simulation Settings

We conducted semi-synthetic simulations based on the kidney renal clear cell carcinoma (KIRC) RNA-seq dataset from The Cancer Genome Atlas (TCGA) [25]. The dataset contains expression measurements for 20,531 genes across 72 paired tumor and adjacent normal tissue samples. Using a real RNA-seq dataset allows the simulations to preserve realistic count distributions, mean–variance relationships, sparsity patterns, and other sources of biological and technical variation commonly observed in ncRNA sequencing studies.

Genes with zero counts in more than 60% of samples were first removed. To construct a set of putative true differentially expressed genes (DEGs), we applied three representative DE methods from different methodological classes: DESeq2 (negative-binomial based), Wilcoxon rank-sum test (nonparametric), and CQRM (quantile-regression based). Genes identified as significant by all three methods at a highly stringent false discovery rate (FDR) threshold of $3e^{-6}$ were designated as true DEGs. The remaining genes were treated as non-DE genes. A semi-synthetic dataset was then generated by preserving the original counts of the true DEGs while randomly permuting sample labels for the non-DE genes between the tumor and normal groups. This strategy preserves realistic expression distributions while establishing a known ground truth for performance evaluation. To better mimic the heterogeneous and noisy expression characteristics frequently observed in ncRNAs, genes with expression variance below the 10th percentile were excluded from subsequent analyses.

Two factors were considered in the simulation:

(i)

Expression level (low vs. moderate counts). To investigate performance under different count regimes, genes were stratified according to their mean expression levels. Genes with mean counts below the 20th percentile were classified as the low-count group, representing the sparse and low-abundance expression patterns commonly observed in ncRNAs. Genes with mean counts between the 20th and 50th percentiles were classified as the moderate-count group.

(ii)

Signal strength (weak to strong): True DEGs were further categorized into three groups of equal size according to increasingly stringent significance thresholds, representing weak, moderate, and strong signals, respectively.

For each simulation replicate, a dataset containing 2,000 genes was generated, with 40% designated as true DEGs and 60% as non-DE genes. A relatively high DEG proportion (40%) was chosen to reflect transcriptome-wide perturbations commonly observed in cancer studies and to ensure stable evaluation across a broad range of signal strengths. Each group (tumor and normal) contained 50 samples subsampled from the original data. Results are based on 50 simulation replicates. We compared nine DE methods: DESeq2, edgeR, limma-voom, Wilcoxon rank-sum test, NOISeq, ZIAQ, CQRM, CQRM-Best, and CQRM-Region. We evaluated both detection power and false discovery rate (FDR) control. Detection power was assessed as the proportion of true DE genes among the top-ranked genes identified by each method. To evaluate error control, we computed the empirical FDR among genes declared significant at nominal FDR thresholds of 0.1, 0.05, 0.01, and 0.001. All results were summarized over 50 independent simulation replicates.

3.1.2. Simulation Results

The simulation results are summarized in Figure 1 and Figure 2. In the low-count setting (Figure 1(d)–(f)), the negative-binomial-based methods DESeq2 and edgeR exhibit substantial FDR inflation, particularly at moderate and large nominal FDR thresholds. Although both methods are highly sensitive and frequently rank among the most powerful approaches, this gain in power comes at the expense of reduced error control. As signal strength decreases, the FDR inflation becomes less pronounced, although both methods remain above the nominal level in most low-count scenarios. In contrast, NOISeq consistently exhibits conservative FDR control across all settings, with empirical FDR well below the nominal thresholds. This conservative behavior is accompanied by substantially lower power than competing methods, particularly under weak-signal conditions. Limma-voom generally maintains reasonable FDR control and is less sensitive to count sparsity than the negative-binomial-based approaches, although its power is typically lower than that of the leading methods. The Wilcoxon rank-sum test demonstrates strong power across nearly all scenarios and often performs comparably to the most powerful methods. However, under low-count settings it exhibits noticeable FDR inflation, particularly at weaker signal strength. ZIAQ demonstrates overall strong performance with both good FDR control and high power in low-count settings. However, its performance deteriorates under weak-signal scenarios, particularly in the moderate-count setting, where it exhibits a noticeable loss of power relative to CQRM methods.

Quantile regression-based methods generally achieve good power with well-controlled FDR, especially in low-count settings. Among them, CQRM consistently provides the best FDR calibration, with empirical FDR values generally closest to the nominal levels across a wide range of settings. CQRM-Best and CQRM-Region tend to be slightly more conservative, often yielding empirical FDR values below the target thresholds. Importantly, all three CQRM variants maintain competitive power while providing substantially improved error control relative to DESeq2, edgeR, and Wilcoxon in challenging low-count scenarios. In moderate-count data with very weak signals, CQRM-based methods demonstrate slightly underpowered behavior compared to DESeq2 and Wilcoxon, which may be because the data are closer to their assumptions in this setting. Another simulation with different sample sizes was conducted, and the power of all DE methods increased with sample size, as expected (Figure S3 and Figure S4).

Taken together, the results suggest that CQRM offers a robust compromise between statistical power and FDR control, particularly in low-count settings that resemble the expression characteristics of many ncRNAs. While DESeq2, edgeR, and Wilcoxon can achieve higher power in some scenarios, this advantage is often accompanied by substantial inflation of false discoveries. In contrast, CQRM maintains competitive sensitivity while providing more reliable error control across a broad range of count levels and signal strengths.

3.2. Case Study 1: DE Analysis of Different ncRNA Species Across Human Cell Lines

To evaluate the performance of CQRM on diverse non-coding RNA species, we analyzed the broad-spectrum RNA sequencing dataset GSE151334 generated using the scRNA-seq tool Smart-seq-total [26]. Unlike conventional RNA-seq protocols that primarily capture mRNAs, Smart-seq-total profiles the entire cellular RNA content, including mRNAs, lncRNAs, miRNAs, snoRNAs, snRNAs, and other non-coding RNA species. We analyzed the total RNAomes of human primary dermal fibroblasts, HEK293T cells (human embryonic kidney cell line), and MCF7 cells (human breast cancer cell line) from the original study. Cells with fewer than 2,000 detected genes or more than two million sequencing reads were excluded. RNAs lacking annotation were removed prior to analysis. Overall, expression levels are low in this dataset, with median counts for both ncRNA and total RNA close to zero and mean counts below 50. Pairwise differential expression analyses were then performed among fibroblasts, HEK293T, and MCF7 cells using edgeR, Wilcoxon, and CQRM. The ZIAQ method encountered errors on this dataset and was therefore excluded from the analysis.

The results are summarized in Figure 3 and Figure 4 (ncRNA includes lncRNAs, miRNAs, snoRNAs, snRNAs and other rare ncRNA species; all RNAs include both mRNAs and ncRNAs). Because the true DE status is unknown in real data, genes identified as significant (Benjamini-Hochberg adjusted $p<0.05$) by multiple methods can be regarded as higher-confidence DE candidates. Therefore, genes consistently detected across methods were used to assess the extent to which CQRM recovers robust DE signals. The Venn diagrams in Figure 3 compare the overlap among CQRM, edgeR, and Wilcoxon. Across all three pairwise cell line comparisons, CQRM shared a substantial proportion of DE genes with both edgeR and Wilcoxon, indicating its ability to recover strong and reproducible DE signals. CQRM also exhibits conservative behavior in DE analysis, with comparatively fewer uniquely detected genes than the other methods, consistent with the stable FDR control observed in simulation.

Inspection of the expression distributions further illustrates the advantage of CQRM. As shown in the box plots (Figure 4), the top 100 DEGs identified by CQRM tended to have lower average expression levels (mean counts) than those identified by edgeR and Wilcoxon for both ncRNAs and all RNAs, suggesting improved sensitivity for low-abundance transcripts commonly encountered among ncRNAs.

3.3. Case Study 2: DE Analysis of lncRNA During Human Organ Development

In the second example, we used CQRM to analyze a large-scale developmental transcriptome dataset spanning multiple organs and developmental stages [27]. Raw bulk RNA-seq count data spanning developmental stages of multiple human organs (brain, cerebellum, heart, kidney, liver, ovary, and testis) were obtained from LnExpDB (https://ngdc.cncb.ac.cn/lncexpdb/downloads; E-MTAB-6814). For the present analysis, only annotated lncRNAs were retained. LncRNAs with detectable expression in fewer than 20% of samples were excluded to reduce the influence of extremely sparse transcripts. Developmental stages were dichotomized into early and late developmental periods. Samples labeled as weeks post-conception or days post-birth were classified as early-stage, whereas samples labeled as years post-birth were classified as late-stage. Organs with insufficient sample sizes or highly imbalanced numbers of early- and late-stage samples (heart, kidney, liver, and ovary) were excluded from downstream analyses. Consequently, differential expression analyses were conducted for brain, cerebellum, and testis. DE analyses were carried out between early and late developmental stages for each organ using edgeR, Wilcoxon, ZIAQ and CQRM methods.

As shown in Figure 5, CQRM identified a large core set of lncRNAs shared with edgeR and Wilcoxon across all three organs. Consistent with the simulation studies and Case Study 1, the lncRNAs prioritized by CQRM tended to exhibit lower average expression levels than those identified by edgeR or Wilcoxon (Figure 6). To evaluate the functional relevance of the identified lncRNAs, we performed enrichment analysis for lncRNAs using the curated lncRNA resource NcPath [28], the result of which is included in supplementary Tables S3–S5.

Top signaling pathways targeted by lncRNAs detected by CQRM in the brain include those involved in neural metabolism and membrane formation (e.g. cholesterol metabolism, fat digestion and absorption, phosphatidylinositol signaling), stem-cell and developmental regulation (e.g., pluripotency signaling and autophagy) and neuroimmune regulation (e.g., cytokine signaling, Th1/Th2 differentiation), all of which have established roles in brain development and maturation [29,30,31]. For cerebellum, genes targeted by lncRNAs detected by CQRM were mainly enriched in pathways involved in neural circuit formation and cerebellar maturation (e.g. axon guidance, Rap1 signaling pathway, circadian rhythm). For testis, genes targeted by lncRNAs detected by CQRM were mainly enriched in pathways involved in endocrine and reproductive signaling (e.g. insulin secretion, thyroid hormone signaling pathway, oxytocin signaling pathway). Notably, these enriched pathways are highly relevant to the biological functions and developmental processes of their respective organs, supporting the ability of CQRM to identify biologically meaningful lncRNAs and their target genes. The distinct enrichment patterns observed across organs further highlight the tissue-specific regulatory roles of lncRNAs during development.

4. Discussion

In this paper, we developed CQRM, a composite quantile regression framework for DE analysis of ncRNA count data. Unlike conventional DE methods that primarily target differences in mean expression, CQRM aggregates information across multiple quantiles and detects a broader range of distributional changes, including shifts in variability, skewness, sparsity, and multimodality. These features are commonly observed across diverse ncRNA species and often reflect biologically meaningful regulatory processes. In semi-synthetic simulations generated based on the TCGA KIRC RNA-seq data, CQRM achieved competitive power while maintaining well controlled FDR across a range of count levels and signal strengths, especially in low-count settings. The two real-data applications further support the practical utility of CQRM. In both the broad-spectrum Smart-seq-total scRNA-seq dataset and the human organ developmental lncRNA sequencing dataset, CQRM demonstrated increased sensitivity for low-abundance ncRNAs while retaining the majority of consensus differential expression signals identified across methods. Moreover, for the organ developmental example, target genes of lncRNAs prioritized by CQRM are enriched in pathways highly related to the developmental process of specific organs, suggesting that the signals detected by CQRM are biologically meaningful rather than artifacts of increased false-positive findings.

To date, the majority of DE methods have been developed for protein-coding genes and rely on distributional assumptions that are often well suited to mRNA expression profiles. Although several recent studies have focused on identifying critical ncRNA regulators and reconstructing ncRNA-mediated regulatory networks [32,33,34], less efforts were made to develop DE methods specifically tailored to the unique characteristics of ncRNA sequencing data. Many ncRNAs exhibit low abundance, high sparsity, high biological variability, and heterogeneous or multimodal expression distributions [9,35]. Consequently, DE may arise through changes in multiple aspects of the expression distribution rather than solely through shifts in mean expression. CQRM addresses this gap by providing a flexible distribution-aware framework that aggregates evidence across multiple quantiles, which can improve statistical efficiency and robustness relative to single-quantile analyses, particularly under heterogeneous and non-Gaussian error distributions [36]. Although quantile regression has previously been applied for DE analysis in scRNA-seq data such as ZIAQ [18], CQRM captures a broader spectrum of DE patterns without relying on explicit zero-inflation assumptions thus shows particular advantages for ncRNA expression data.

The current framework analyzes each ncRNA independently and therefore does not explicitly leverage regulatory relationships, co-expression patterns, or pathway structure that may further improve power and biological interpretability. Future work will focus on integrating prior biological knowledge, gene regulatory networks, and multi-omics data within the composite quantile framework to better prioritize critical ncRNA markers and characterize ncRNA-mediated regulatory mechanisms. An R package implementing CQRM, together with codes required to reproduce the analyses presented in this study, are publicly available at: https://github.com/iamverywell/CQRM.