"predFC"(y, design=NULL, prior.count=0.125, offset=NULL, dispersion=NULL, weights=NULL, ...) "predFC"(y, design=NULL, prior.count=0.125, offset=NULL, dispersion=0, weights=NULL, ...)
glmFit. Usually equal to log library sizes.
designis given) or logCPM (if
design=NULL) on the log2 scale.
Specifically the function adds a small prior count to each observation before estimating the glm. The actual prior count that is added is proportion to the library size. This has the effect that any log-fold-change that was zero prior to augmentation remains zero and non-zero log-fold-changes are shrunk towards zero.
The prior counts can be viewed as equivalent to a prior belief that the log-fold changes are small, and the output can be viewed as posterior log-fold-changes from this Bayesian viewpoint. The output coefficients are called predictive log fold-changes because, depending on the prior, they may be a better prediction of the true log fold-changes than the raw estimates.
Log-fold changes for genes with low counts are shrunk more than those for genes with high counts. In particular, infinite log-fold-changes arising from zero counts are avoided. The exact degree to which this is done depends on the negative binomail dispersion.
design=NULL, then the function returns a matrix of the same size as
y containing log2 counts-per-million, with zero values for the counts avoided.
This equivalent to choosing
design to be the identity matrix with the same number of columns as
# generate counts for a two group experiment with n=2 in each group and 100 genes dispersion <- 0.1 y <- matrix(rnbinom(400,size=1/dispersion,mu=4),nrow=100) y <- DGEList(y,group=c(1,1,2,2)) design <- model.matrix(~group, data=y$samples) #estimate the predictive log fold changes predlfc<-predFC(y,design,dispersion=dispersion,prior.count=1) logfc <- predFC(y,design,dispersion=dispersion,prior.count=0) logfc.truncated <- pmax(pmin(logfc,100),-100) #plot predFC's vs logFC's plot(predlfc[,2],logfc.truncated[,2],xlab="Predictive log fold changes",ylab="Raw log fold changes") abline(a=0,b=1)