Apply DESeq2 for comparisons defined in a contrast matrix, requires the DESeq2 package.
Usage
PairwiseDESeq2(
data,
name.prefix = mode,
contrasts,
separate = FALSE,
mode = "total",
slot = "count",
normalization = NULL,
logFC = FALSE,
genes = NULL,
verbose = FALSE
)Arguments
- data
the grandR object
- name.prefix
the prefix for the new analysis name; a dot and the column names of the contrast matrix are appended; can be NULL (then only the contrast matrix names are used)
- contrasts
contrast matrix that defines all pairwise comparisons, generated using GetContrasts
- separate
model overdispersion separately for all pairwise comparison (TRUE), or fit a single model per gene, and extract contrasts (FALSE)
- mode
compute LFCs for "total", "new", or "old" RNA
- slot
which slot to use (should be a count slot, not normalized values)
- normalization
normalize on "total", "new", or "old" (see details)
- logFC
compute and add the log2 fold change as well
- genes
restrict analysis to these genes; NULL means all genes
- verbose
print status messages?
Value
a new grandR object including a new analysis table. The columns of the new analysis table are
- "M"
the base mean
- "S"
the log2FoldChange divided by lfcSE
- "P"
the Wald test P value
- "Q"
same as P but Benjamini-Hochberg multiple testing corrected
- "LFC"
the log2 fold change (only with the logFC parameter set to TRUE)
Details
DESeq2 by default performs size factor normalization. When computing differential expression of new RNA, it might be sensible to normalize w.r.t. to total RNA, i.e. use the size factors computed from total RNA instead of computed from new RNA. This can be accomplished by setting mode to "new", and normalization to "total"!
Normalization can also be a mode.slot! Importantly, do not specify a slot containing normalized values, but specify a slot of unnormalized values (which are used to compute the size factors for normalization!) Can also be a numeric vector of size factors with the same length as the data as columns. Then each value is divided by the corresponding size factor entry.
Examples
# \donttest{
sars <- ReadGRAND(system.file("extdata", "sars.tsv.gz", package = "grandR"),
design=c(Design$Condition,Design$dur.4sU,Design$Replicate))
#> Warning: Duplicate gene symbols (n=1, e.g. MATR3) present, making unique!
sars <- subset(sars,Coldata(sars,Design$dur.4sU)==2)
sars<-PairwiseDESeq2(sars,mode="total",
contrasts=GetContrasts(sars,contrast=c("Condition","Mock")))
sars<-PairwiseDESeq2(sars,mode="new",normalization="total",
contrasts=GetContrasts(sars,contrast=c("Condition","Mock")))
#> -- 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.
head(GetAnalysisTable(sars,column="Q"))
#> Gene Symbol Length Type total.SARS vs Mock.Q
#> UHMK1 ENSG00000152332 UHMK1 8478 Cellular 1.915748e-01
#> ATF3 ENSG00000162772 ATF3 2103 Cellular 5.615492e-10
#> PABPC4 ENSG00000090621 PABPC4 3592 Cellular 6.620783e-01
#> ROR1 ENSG00000185483 ROR1 5832 Cellular 5.054219e-01
#> ZC3H11A ENSG00000058673 ZC3H11A 11825 Cellular 2.141239e-02
#> ZBED6 ENSG00000257315 ZBED6 12481 Cellular 2.402991e-02
#> new.SARS vs Mock.Q
#> UHMK1 9.572768e-04
#> ATF3 1.178977e-27
#> PABPC4 5.475308e-08
#> ROR1 2.997945e-03
#> ZC3H11A 4.686341e-07
#> ZBED6 3.173834e-07
# }