Description Usage Arguments Details Value Note Author(s) References See Also Examples
Fit a quasi-likelihood negative binomial generalized log-linear model to count data. Conduct genewise statistical tests for a given coefficient or contrast.
1 2 3 4 5 6 7 8 9 10 11 | ## S3 method for class 'DGEList'
glmQLFit(y, design=NULL, dispersion=NULL, abundance.trend=TRUE,
robust=FALSE, winsor.tail.p=c(0.05, 0.1), ...)
## S3 method for class 'SummarizedExperiment'
glmQLFit(y, design=NULL, dispersion=NULL, abundance.trend=TRUE,
robust=FALSE, winsor.tail.p=c(0.05, 0.1), ...)
## Default S3 method:
glmQLFit(y, design=NULL, dispersion=NULL, offset=NULL, lib.size=NULL,
weights=NULL, abundance.trend=TRUE, AveLogCPM=NULL, robust=FALSE,
winsor.tail.p=c(0.05, 0.1), ...)
glmQLFTest(glmfit, coef=ncol(glmfit$design), contrast=NULL, poisson.bound=TRUE)
|
y |
a matrix of counts, or a |
design |
numeric matrix giving the design matrix for the genewise linear models. |
dispersion |
numeric scalar, vector or matrix of negative binomial dispersions. If |
offset |
numeric matrix of same size as |
lib.size |
numeric vector of length |
weights |
numeric matrix of same size as |
abundance.trend |
logical, whether to allow an abundance-dependent trend when estimating the prior values for the quasi-likelihood multiplicative dispersion parameter. |
AveLogCPM |
average log2-counts per million, the average taken over all libraries in |
robust |
logical, whether to estimate the prior QL dispersion distribution robustly. |
winsor.tail.p |
numeric vector of length 2 giving proportion to trim (Winsorize) from lower and upper tail of the distribution of genewise deviances when estimating the hyperparameters. Positive values produce robust empirical Bayes ignoring outlier small or large deviances. Only used when |
... |
other arguments are passed to |
glmfit |
a |
coef |
integer or character index vector indicating which coefficients of the linear model are to be tested equal to zero. Ignored if |
contrast |
numeric vector or matrix specifying one or more contrasts of the linear model coefficients to be tested equal to zero. |
poisson.bound |
logical, if |
glmQLFit
and glmQLFTest
implement the quasi-likelihood (QL) methods of Lund et al (2012), with some enhancements and with slightly different glm, trend and FDR methods.
See Lun et al (2016) or Chen et al (2016) for tutorials describing the use of glmQLFit
and glmQLFit
as part of a complete analysis pipeline.
Another case study using glmQLFit
and glmQLFTest
is given in Section 4.7 of the edgeR User's Guide.
glmQLFit
is similar to glmFit
except that it also estimates QL dispersion values.
It calls the limma function squeezeVar
to conduct empirical Bayes moderation of the genewise QL dispersions.
If robust=TRUE
, then the robust hyperparameter estimation features of squeezeVar
are used (Phipson et al, 2013).
If abundance.trend=TRUE
, then a prior trend is estimated based on the average logCPMs.
glmQLFit
gives special attention to handling of zero counts, and in particular to situations when fitted values of zero provide no useful residual degrees of freedom for estimating the QL dispersion (Lun and Smyth, 2017).
The usual residual degrees of freedom are returned as df.residual
while the adjusted residual degrees of freedom are returned as df.residuals.zeros
.
glmQLFTest
is similar to glmLRT
except that it replaces likelihood ratio tests with empirical Bayes quasi-likelihood F-tests.
The p-values from glmQLFTest
are always greater than or equal to those that would be obtained from glmLRT
using the same negative binomial dispersions.
glmQLFit
produces an object of class DGEGLM
with the same components as produced by glmFit
, plus:
df.residual.zeros |
a numeric vector containing the number of effective residual degrees of freedom for each gene, taking into account any treatment groups with all zero counts. |
df.prior |
a numeric vector or scalar, giving the prior degrees of freedom for the QL dispersions. |
var.prior |
a numeric vector of scalar, giving the location of the prior distribution for the QL dispersions. |
var.post |
a numeric vector containing the posterior empirical Bayes QL dispersions. |
df.prior
is a vector of length nrow(y)
if robust=TRUE
, otherwise it has length 1.
var.prior
is a vector of length nrow(y)
if abundance.trend=TRUE
, otherwise it has length 1.
glmQFTest
produce an object of class DGELRT
with the same components as produced by glmLRT
, except that the table$LR
column becomes table$F
and contains quasi-likelihood F-statistics.
It also stores df.total
, a numeric vector containing the denominator degrees of freedom for the F-test, equal to df.prior + df.residual.zeros
.
The negative binomial dispersions dispersion
supplied to glmQLFit
and glmQLFTest
must be based on a global model, that is, they must be either trended or common dispersions.
It is not correct to supply genewise dispersions because glmQLFTest
estimates genewise variability using the QL dispersion.
Yunshun Chen, Aaron Lun, Davis McCarthy and Gordon Smyth
Chen Y, Lun ATL, and Smyth, GK (2016). From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5, 1438. http://f1000research.com/articles/5-1438
Lun, ATL, Chen, Y, and Smyth, GK (2016). It's DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods in Molecular Biology 1418, 391-416. http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf (Preprint 8 April 2015)
Lund, SP, Nettleton, D, McCarthy, DJ, and Smyth, GK (2012). Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Statistical Applications in Genetics and Molecular Biology Volume 11, Issue 5, Article 8. http://www.statsci.org/smyth/pubs/QuasiSeqPreprint.pdf
Lun, ATL, and Smyth, GK (2017). No counts, no variance: allowing for loss of degrees of freedom when assessing biological variability from RNA-seq data. Statistical Applications in Genetics and Molecular Biology 16(2), 83-93. https://doi.org/10.1515/sagmb-2017-0010
Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. https://projecteuclid.org/euclid.aoas/1469199900
topTags
displays results from glmQLFTest
.
plotQLDisp
can be used to visualize the distribution of QL dispersions after EB shrinkage from glmQLFit
.
The QuasiSeq
package gives an alternative implementation of the Lund et al (2012) methods.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | nlibs <- 4
ngenes <- 1000
dispersion.true <- 1/rchisq(ngenes, df=10)
design <- model.matrix(~factor(c(1,1,2,2)))
# Generate count data
y <- rnbinom(ngenes*nlibs,mu=20,size=1/dispersion.true)
y <- matrix(y,ngenes,nlibs)
d <- DGEList(y)
d <- calcNormFactors(d)
# Fit the NB GLMs with QL methods
d <- estimateDisp(d, design)
fit <- glmQLFit(d, design)
results <- glmQLFTest(fit)
topTags(results)
fit <- glmQLFit(d, design, robust=TRUE)
results <- glmQLFTest(fit)
topTags(results)
fit <- glmQLFit(d, design, abundance.trend=FALSE)
results <- glmQLFTest(fit)
topTags(results)
|
Loading required package: limma
Coefficient: factor(c(1, 1, 2, 2))2
logFC logCPM F PValue FDR
118 -2.611611 9.962561 14.523192 0.001769154 0.7560927
728 2.186866 9.885578 13.775226 0.002163362 0.7560927
31 1.993878 10.029548 13.602231 0.002268278 0.7560927
914 -1.707966 10.370731 9.607287 0.007488196 0.9416704
315 1.680880 10.365811 9.371487 0.008091077 0.9416704
192 1.561621 10.096980 8.859932 0.009601046 0.9416704
300 -1.608652 10.411703 8.176986 0.012150275 0.9416704
673 1.586165 9.904900 8.084169 0.012553623 0.9416704
1000 1.471175 10.161560 7.925976 0.013277154 0.9416704
217 -1.350659 10.068246 7.350414 0.016347195 0.9416704
Coefficient: factor(c(1, 1, 2, 2))2
logFC logCPM F PValue FDR
118 -2.611611 9.962561 15.172781 0.0009602732 0.7257361
31 1.993878 10.029548 13.643294 0.0015025209 0.7257361
728 2.186866 9.885578 12.474696 0.0021772083 0.7257361
914 -1.707966 10.370731 9.283936 0.0065324780 0.9188915
315 1.680880 10.365811 9.025923 0.0071846788 0.9188915
192 1.561621 10.096980 8.701079 0.0081115075 0.9188915
300 -1.608652 10.411703 8.089675 0.0102415132 0.9188915
673 1.586165 9.904900 7.424790 0.0132975830 0.9188915
1000 1.471175 10.161560 7.214981 0.0144648892 0.9188915
217 -1.350659 10.068246 7.184654 0.0146429392 0.9188915
Coefficient: factor(c(1, 1, 2, 2))2
logFC logCPM F PValue FDR
118 -2.611611 9.962561 14.472546 0.001775936 0.8180217
728 2.186866 9.885578 13.996117 0.002018623 0.8180217
31 1.993878 10.029548 13.286001 0.002454065 0.8180217
914 -1.707966 10.370731 9.215992 0.008472571 0.9407065
315 1.680880 10.365811 8.981174 0.009167113 0.9407065
421 -2.335122 9.314667 8.979909 0.009171027 0.9407065
192 1.561621 10.096980 8.546197 0.010633920 0.9407065
673 1.586165 9.904900 8.163605 0.012150556 0.9407065
300 -1.608652 10.411703 7.910448 0.013290996 0.9407065
213 1.760680 9.682400 7.866366 0.013501941 0.9407065
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.