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.
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
)A 'glmDataSet' object, which is created with function
glmDataSet.
Parameters for parallel 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).
whether or not to filter the counts according to the minimum count per region per each individual/sample, which is setting by 'minCountPerIndv'.
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
each gene or region must have more than 'minCountPerIndv' counts (on average) per individual in at least one group.
A numerical vector. Maximum coefficient of variance for each group. Default maxGrpCV = NULL. The numbers 'maxGrpCV1' and 'maxGrpCV2' will be taken as the maximun variances values permitted in control and in treatment groups, respectively. If only 'maxGrpCV1' is provided, then maxGrpCV = c('maxGrpCV1', 'maxGrpCV1'). This parameter is addressed to prevent testing regions where intra-group variations are very large, e.g.: control = c(1,0,1,1) and treatment = 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.
if TRUE, the results are filtered using the minimum absolute value of log2FoldChanges observed to accept that a gene in the treatment is differentially expressed in respect to the control
method used to adjust the results; default: BH
cutoff used then a p-value adjustment is performed.
Minimum Mean/Variance rate.
minimum logarithm base 2 of fold changes.
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'.
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.
if TRUE all the temporal results are returned
if TRUE, prints the function log to stdout
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.
A pairwise group comparison, control versus treatment, is performed.
The experimental design settings must be introduced using function
glmDataSet to provide dataset (DS) object.
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)
#> GRanges object with 4 ranges and 12 metadata columns:
#> seqnames ranges strand | A1 A2 B1
#> <Rle> <IRanges> <Rle> | <integer> <integer> <integer>
#> gene69 chr1 136001-137001 * | 57 67 117
#> gene98 chr2 194001-195001 * | 34 39 116
#> gene101 chr2 200001-201001 * | 166 168 117
#> gene114 chr2 226001-227001 * | 33 34 111
#> B2 log2FC scaled.deviance pvalue model
#> <integer> <numeric> <numeric> <numeric> <character>
#> gene69 102 0.568790 26.6591 4.17025e-07 Neg.Binomial
#> gene98 151 1.296789 82.6185 9.57148e-18 Neg.Binomial
#> gene101 71 -0.574699 13.0830 3.00075e-04 Neg.Binomial
#> gene114 158 1.390019 64.7083 7.92421e-15 Neg.Binomial
#> adj.pval CT.SignalDensity TT.SignalDensity SignalDensityVariation
#> <numeric> <numeric> <numeric> <numeric>
#> gene69 5.56034e-07 0.0619381 0.1093906 0.0474525
#> gene98 3.82859e-17 0.0364635 0.1333666 0.0969031
#> gene101 3.00075e-04 0.1668332 0.0939061 -0.0729271
#> gene114 1.58484e-14 0.0334665 0.1343656 0.1008991
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths