statistic tests for HTSet

This is a wrapper of several widely used statistical method for high through-put experimental data such as RNAseq. The limma-package performs linear model on continous data, or cooperate with the voom function to handle count data. The edgeR-package and DESeq2-package performs negative binomial generalized linear model on count data.

model_fit(
  object,
  design,
  coef,
  engine,
  args = list(),
  transform,
  adjust.method = "BH"
)

Arguments

object	HTSet
design	matrix. Number of rows must batch number of samples of object. Usually the output of model.matrix.
coef	character. The coefficient to perform statistical tests. Must be in the column names of design.
engine	character. The engine to perform statistical analysis. Supported are limma, edgeR, and DESeq2.
args	list. A list of argumnets to be parsed to the backend statistical engine.
transform	function. The transform to be passed to transform_by_sample. No transform is performed by default.
adjust.method	character. Method used to adjust the p-values for multiple testing. See p.adjust for the complete list of supported methods. Default is "BH" for Benjamini-Hochberg test.

Value

A list-like S3 class ModelFit object is returned with the elements as following.

results

A data.frame of the statistical test results for each gene/ feature.

logFC

estimate of the log2-fold-change corresponding to the effect.

mean

average log2-expression for the gene/feature accross all samples. Same as the AveExpr in limma's topTable, and the baseMean in DESeq2's results.

stat

the statistic value of the corresponding test. When using limma, this is the t-statistic value, same as the t column in the result of topTable. For edgeR, this is the f-statistic value for Quansi-likelihood test or the likelihood ratio statistic value for likelihood ratio test. Same as the column F or RT in the result of edgeR's topTags. As for DESeq2, this is the Wald statistic value for Wald test or the difference in deviance between the reduced model and the full model for likelihood ratio test. same as the stat column in the output of DESeq2's results.

pval

raw p-value

padj

adjusted p-value

adjust.method

Method used to correct for multiple testing

design

design matrix

df

degree of freedoms

distribution

The distribution that p values were calculated

engine

package used for statistical test

coef

coefficient tested

params

additional parameters parsed

See also

HTSet-class lmFit voom glmQLFit glmLRT DESeq

Examples

data(exrna)
design = model.matrix(~ Condition, data = exrna$pdata)
coef = "ConditionSystemic Lupus Erythematosus"
fit1 = model_fit(object = exrna, design = design, coef = coef, engine = "limma", args = list(voom = TRUE))
fit2 = model_fit(object = exrna, design = design, coef = coef, engine = "edgeR")
fit3 = model_fit(object = exrna, design = design, coef = coef, engine = "edgeR", args = list(model = "lrt"))
fit4 = model_fit(object = exrna, design = design, coef = coef, engine = "DESeq2")
#> Loading required namespace: DESeq2
#> using supplied model matrix
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> -- note: fitType='parametric', but the dispersion trend was not well captured by the
#>    function: y = a/x + b, and a local regression fit was automatically substituted.
#>    specify fitType='local' or 'mean' to avoid this message next time.
#> final dispersion estimates
#> fitting model and testing