| Title: | Differential Expression Analysis with Long Read RNA-Seq Data |
|---|---|
| Description: | Provides hurdle negative binomial models for differential expression analysis with long-read RNA-Seq data. |
| Authors: | Ziyang Liu [aut, cre] (ORCID: <https://orcid.org/0009-0004-2098-434X>), Hongxu Ding [aut, fnd], Xiaoxiao Sun [aut, fnd], Ziyuan Wang [aut, fnd] |
| Maintainer: | Ziyang Liu <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.1.1 |
| Built: | 2026-06-25 21:34:23 UTC |
| Source: | https://github.com/bioc/LRDE |
The LRDE package provides statistical methods for differential expression analysis of long-read RNA sequencing (RNA-Seq) data using a hurdle negative binomial generalized linear model (hurdle-NB GLM).
It implements procedures for:
Estimation of sample-specific size factors for normalization
Modeling zero inflation via group-specific expression probabilities
Gene-wise (tag-wise) dispersion estimation
Statistical testing for differential expression
These methods are designed to address key challenges in long-read RNA-Seq data, including limited sample sizes and excess zero counts (dropout events).
The main functions in this package include:
prepareDGE: Prepare count data for analysis. Converts supported
input types (matrix, data.frame, DGEList,
DESeqDataSet, and SummarizedExperiment) into a standardized format.
sizeFactorsEst: Estimate sample-specific size factors for normalization.
tagwiseEst: Estimate gene-specific (tag-wise) dispersion parameters
for a hurdle negative binomial model using prior information from bin-level estimates.
hurdle.LRT: Perform gene-wise likelihood ratio tests (LRT) for differential expression.
hurdle.Wald.Test: Perform gene-wise Wald tests for differential expression.
Typical workflow:
Prepare data using prepareDGE
Normalize counts with sizeFactorsEst
Estimate tag-wise dispersions using tagwiseEst
Perform differential expression testing with
hurdle.LRT or hurdle.Wald.Test
Ziyang Liu [email protected]
prepareDGE,
sizeFactorsEst,
tagwiseEst,
hurdle.LRT,
hurdle.Wald.Test
# Load the package library(LRDE) # Simulate count data set.seed(123) mat <- matrix(rnbinom(300, size = 5, mu = 5), nrow = 50) grp <- factor(c("A", "A", "A", "B", "B", "B")) # Prepare data y <- prepareDGE(mat, grp) # Normalize counts y <- sizeFactorsEst(y) # Estimate dispersions y <- tagwiseEst(y) # Differential expression testing y <- hurdle.Wald.Test(y) y <- hurdle.LRT(y) # Access results head(y$lrt_stats) head(y$p.values)# Load the package library(LRDE) # Simulate count data set.seed(123) mat <- matrix(rnbinom(300, size = 5, mu = 5), nrow = 50) grp <- factor(c("A", "A", "A", "B", "B", "B")) # Prepare data y <- prepareDGE(mat, grp) # Normalize counts y <- sizeFactorsEst(y) # Estimate dispersions y <- tagwiseEst(y) # Differential expression testing y <- hurdle.Wald.Test(y) y <- hurdle.LRT(y) # Access results head(y$lrt_stats) head(y$p.values)
Performs gene-wise likelihood ratio tests (LRTs) for distributional differential expression using a hurdle negative binomial model. The test jointly evaluates differences in the zero probability and positive-count mean between groups while holding the tag-wise dispersion fixed.
hurdle_LRT.dist(y)hurdle_LRT.dist(y)
y |
A list-like object returned from
|
For each gene:
The null model assumes a shared nonzero probability and a shared positive-count mean across groups.
The alternative model estimates group-specific nonzero probabilities and positive-count means.
The dispersion is fixed at its tag-wise estimate obtained from
tagwiseEst.
The LRT uses two degrees of freedom when both groups contain positive counts and one degree of freedom when either group contains only zero counts.
Thus, rejection of the null hypothesis indicates a difference in the overall expression distribution arising from the zero probability, positive-count mean, or both.
The input object y with two additional elements:
Numeric vector of LRT statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_LRT.dist(y)set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_LRT.dist(y)
Performs gene-wise likelihood ratio tests (LRTs) for differential expression in the expression mean using a hurdle negative binomial model with fixed zero probabilities and tag-wise dispersions.
hurdle_LRT.mean(y)hurdle_LRT.mean(y)
y |
A list-like object returned from
|
For each gene:
The null model assumes a single shared mean across groups.
The alternative model estimates group-specific means.
Zero probabilities and dispersions are fixed at estimates obtained
from tagwiseEst.
When one group has all zero counts, a one-sided Z-test is applied instead.
The input object y with two additional elements:
Numeric vector of LRT statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_LRT.mean(y)set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_LRT.mean(y)
Performs gene-wise Wald tests for distributional differential expression using a hurdle negative binomial model. The test jointly evaluates differences in the zero probability and positive-count mean between groups while holding the tag-wise dispersion fixed.
hurdle_Wald_Test.dist(y)hurdle_Wald_Test.dist(y)
y |
A list-like object returned from
|
For each gene:
When both groups contain positive counts, the model estimates group-specific nonzero probabilities and mean parameters while fixing the dispersion at its tag-wise estimate.
A joint two-degree-of-freedom Wald test evaluates the nonzero probability and positive-count mean differences between groups.
When either group contains only zero counts, the test compares only the nonzero probabilities using corrected log odds, with 0.5 added to each positive- and zero-count frequency.
The resulting squared Z statistic is compared with a chi-squared distribution with one degree of freedom in the all-zero-group case.
Rejection of the null hypothesis indicates a difference in the overall expression distribution arising from the zero probability, positive-count mean, or both.
The input object y with two additional elements:
Numeric vector of Wald statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_Wald_Test.dist(y)set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_Wald_Test.dist(y)
Performs gene-wise Wald tests for differential expression in the expression mean using a hurdle negative binomial model with fixed zero probabilities and tag-wise dispersions.
hurdle_Wald_Test.mean(y)hurdle_Wald_Test.mean(y)
y |
A list-like object returned from
|
For each gene:
The model estimates group-specific mean parameters while fixing the
zero probabilities and dispersion at estimates obtained from
tagwiseEst.
When both groups contain positive counts, a two-sided Wald test is applied to the log mean difference between groups.
When either group contains only zero counts, a one-sided Z-test is applied to the estimated log mean of the nonzero group using the corresponding Hessian-based standard error.
The input object y with two additional elements:
Numeric vector of Wald or one-sided Z statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_Wald_Test.mean(y)set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) y <- hurdle_Wald_Test.mean(y)
Performs gene-wise differential expression testing using a hurdle negative binomial model. The function can test differences in either the expression mean or the overall expression distribution between groups.
hurdle.LRT(y, test = c("mean", "distribution"))hurdle.LRT(y, test = c("mean", "distribution"))
y |
A list-like object returned from |
test |
Character string specifying the hypothesis to test. Use
|
When test = "mean", the function calls
hurdle_LRT.mean. When test = "distribution", it calls
hurdle_LRT.dist.
The input object y with two additional elements:
Numeric vector of test statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) mean_test <- hurdle.LRT(y, test = "mean") distribution_test <- hurdle.LRT(y, test = "distribution")set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) mean_test <- hurdle.LRT(y, test = "mean") distribution_test <- hurdle.LRT(y, test = "distribution")
Performs gene-wise differential expression testing using a hurdle negative binomial model. The function can test differences in either the expression mean or the overall expression distribution between groups.
hurdle.Wald.Test(y, test = c("mean", "distribution"))hurdle.Wald.Test(y, test = c("mean", "distribution"))
y |
A list-like object returned from |
test |
Character string specifying the hypothesis to test. Use
|
When test = "mean", the function calls
hurdle_Wald_Test.mean. When
test = "distribution", it calls
hurdle_Wald_Test.dist.
The input object y with two additional elements:
Numeric vector of Wald statistics for each gene.
Numeric vector of corresponding p-values.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) mean_test <- hurdle.Wald.Test(y, test = "mean") distribution_test <- hurdle.Wald.Test(y, test = "distribution")set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) mean_test <- hurdle.Wald.Test(y, test = "mean") distribution_test <- hurdle.Wald.Test(y, test = "distribution")
Converts various supported input types to a standardized list format for downstream
differential expression analysis. Supports matrix, data.frame, DGEList,
DESeqDataSet, and SummarizedExperiment objects.
prepareDGE(data, group)prepareDGE(data, group)
data |
A numeric matrix, data.frame, or supported object containing counts. |
group |
A vector of group labels for the columns/samples of |
This function performs input validation.
Checks for non-negative numeric values and absence of NA.
Ensures group labels match the number of samples.
Automatically assigns column names if missing.
Returns a list suitable for use with hurdle model-based DE functions.
A list with two elements:
An integer matrix of counts.
A data.frame containing sample-level metadata: group, lib.size, and size.factor.
# Example with a matrix set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) # Example with a SummarizedExperiment if (requireNamespace("SummarizedExperiment", quietly = TRUE)) { se <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = mat)) y_se <- prepareDGE(se, grp) y }# Example with a matrix set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) # Example with a SummarizedExperiment if (requireNamespace("SummarizedExperiment", quietly = TRUE)) { se <- SummarizedExperiment::SummarizedExperiment(assays = list(counts = mat)) y_se <- prepareDGE(se, grp) y }
Computes sample-specific size factors for long-read RNA-Seq data, used to normalize counts for differential expression analysis.
sizeFactorsEst(y, type = c("poscounts", "ratio"), locfunc = stats::median)sizeFactorsEst(y, type = c("poscounts", "ratio"), locfunc = stats::median)
y |
A count matrix ( |
type |
Character string specifying the method for estimating size factors:
Default is |
locfunc |
Function to summarize log-ratios across genes. Defaults to |
This function implements two methods for size factor estimation:
poscounts: Computes a geometric mean of positive counts per gene, then calculates ratios for each sample. Normalizes so that the geometric mean of size factors equals 1.
ratio: Uses the mean of log-counts per gene across samples to compute ratios.
The function automatically normalizes counts using the estimated size factors
and stores gene-level normalized means in baseMean.
A list (same structure as prepareDGE() output) with:
Original count matrix (integer).
Data frame with sample information, updated size.factor.
Normalized mean of counts per gene.
# Using a count matrix #' set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y, type = "poscounts")# Using a count matrix #' set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y, type = "poscounts")
Estimate gene-specific (tag-wise) dispersion parameters for a hurdle negative binomial model using prior information derived from bin-level estimates.
tagwiseEst(y, small.sample = TRUE, prior = TRUE)tagwiseEst(y, small.sample = TRUE, prior = TRUE)
y |
A list object created by |
small.sample |
Logical. If |
prior |
Logical. Whether to use prior information for dispersion
shrinkage when |
When small.sample = TRUE, the function calls
tagwiseEst.smallSample, which performs the following steps:
Retrieves bin-level prior estimates of zero probabilities and log-dispersion
for each gene via priorEst.
Fixes the zero probabilities and optimizes only the mean parameters and dispersion for each gene individually.
Uses the internal function nll_hurdle_fixed_P to compute the negative
log-likelihood with fixed zero probabilities.
When small.sample = FALSE, the function calls
tagwiseEst.largeSample, which fits a hurdle negative binomial model
separately to each gene without borrowing information across genes.
The resulting tagwise.disp will be used for downstream differential
expression analysis.
The input y object augmented with:
Numeric vector of estimated gene-wise dispersions.
Numeric matrix of fixed zero probabilities for each gene and group.
set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) head(y$tagwise.disp) head(y$zero_prob_matrix)set.seed(123) mat <- matrix(rnbinom(30, size = 5, mu = 5), nrow = 5) grp <- c("A", "A", "A", "B", "B", "B") y <- prepareDGE(mat, grp) y <- sizeFactorsEst(y) y <- tagwiseEst(y) head(y$tagwise.disp) head(y$zero_prob_matrix)