countTest2 {MethylIT}R Documentation

Regression Test for Count

Description

Perform Poisson and Negative Binomial regression analysis to compare the counts from different groups, treatment and control. The difference between functions 'countTest2' and 'countTest' resides in the estimation of the prior weights used in Negative Binomial generalized linear model.

Usage

countTest2(DS, num.cores = 1, countFilter = TRUE, CountPerBp = NULL,
  minCountPerIndv = 3, maxGrpCV = NULL, FilterLog2FC = TRUE,
  pAdjustMethod = "BH", pvalCutOff = 0.05, MVrate = 0.98,
  Minlog2FC = 0.5, test = c("Wald", "LRT"), scaling = 1L,
  tasks = 0L, saveAll = FALSE, verbose = TRUE)

Arguments

DS

A 'glmDataSet' object, which is created with function link{glmDataSet}.

num.cores, tasks

Paramaters for parallele computation using package BiocParallel-package: the number of cores to use, i.e. at most how many child processes will be run simultaneously (see bplapply and the number of tasks per job (only for Linux OS).

countFilter

whether or not to filter the counts according to the minimum count per region per each individual/sample, which is setting by "minCountPerIndv".

CountPerBp

for each group the count per bp must be equal or greater than CountPerBp. The filter is applied if 'CountPerBp' is given and if 'x' DESeqDataSet object has the rowRanges as a GRanges object on it

minCountPerIndv

each gene or region must have more than 'minCountPerIndv' counts (on average) per individual in at least one group.

maxGrpCV

A numerical vector. Maximum coefficient of variance for each group. Defaul maxGrpCV = NULL. The numbers maxGrpCV[1] and maxGrpCV[2] will be taken as the maximun variances values permitted in control and in treatment groups, repectively. If only maxGrpCV[1] is provided, then maxGrpCV = c(maxGrpCV[1], maxGrpCV[1]). This parameter is addressed to prevent testing regions where intra-group variations are very large, e.g.: control = c(1,0,1,1) and traatment = c(1, 0, 1, 40). The coefficient of variance for the treatment group is 1.87, very high. The generalized linear regression analysis would yield statistical significant group differences, but evidently there is something wrong in one of the treatment samples. We would try the application of further statistical smoothing approach, but we prefer to leave the user decide which regions to test.

FilterLog2FC

if TRUE, the results are filtered using the minimun absolute value of log2FoldChanges observed to accept that a gene in the treatment is differentially expressed in respect to the control

pAdjustMethod

method used to adjust the results; default: BH

pvalCutOff

cutoff used then a p-value adjustment is performed

MVrate

Minimum Mean/Variance rate.

Minlog2FC

minimum logarithm base 2 of fold changes.

test

A character string matching one of "Wald" or "LRT". If test = "Wald", then the p-value of the Wald test for the coefficient of the independent variable (treatment group) will be reported. If test = "LRT", then the p-value from a likelihood ratio test given by anova function from stats packages will be the reported p-value for the group comparison when the best fitted model is the negative binomial. As suggested for glm, if best fitted model is Poisson or quasi-Poisson, then the best test is 'Chi-squared' or 'F-test', respectively. So, for the sake of simplicity, the corresponding suitable test will be applied when test = "LRT".

scaling

integer (default 1). Scaling factor estimate the signal density as: scaling x "DMP-Count-Per-Bp". For example, if scaling = 1000, then signal density denotes the number of DMPs in 1000 bp.

saveAll

if TRUE all the temporal results are returned

verbose

if TRUE, prints the function log to stdout

Details

A pairwise group comparison, control versus treatment, is performed. The experimental design settings must be introduced using function link{glmDataSet} to provide dataset (DS) object.

Value

a data frame or GRanges object (if the DS contain the GRanges information for each gene) with the test results and original count matrix, plus control and treatment signal densities and their variation.

Examples

set.seed(133) # Set a seed
## A GRanges object with the count matrix in the metacolumns is created
countData <- matrix(sample.int(200, 500, replace = TRUE), ncol = 4)
colnames(countData) <- c("A1","A2","B1","B2")
start <- seq(1, 25e4, 2000)
end <- start + 1000
chr <- c(rep("chr1", 70), rep("chr2", 55))
GR <- GRanges(seqnames = chr, IRanges(start = start, end = end))
mcols(GR) <- countData
## Gene IDs
names(GR) <- paste0("gene", 1:length(GR))

## An experiment design is set.
colData <- data.frame(condition = factor(c("A","A","B","B")),
                      c("A1","A2","B1","B2"),
                      row.names = 2)
## A RangedGlmDataSet is created
ds <- glmDataSet(GR = GR, colData = colData)

## The gneralized linear model pairwise group comparison, group 'A'
## ('control') versus 'B' (treatment) is performed.
countTest2(ds, num.cores = 1L, maxGrpCV = c(0.4, 0.4), verbose = FALSE)


[Package MethylIT version 0.3.1 ]