.Rmd file containing code to perform differential expression analysis with DESeq2 with custom model matrixR/generateRmdCodeDiffExpPhylo.R
DESeq2.length.createRmd.RdA function to generate code that can be run to perform differential expression analysis of RNAseq data (comparing two conditions) using the DESeq2 package. The code is written to a .Rmd file. This function is generally not called by the user, the main interface for performing differential expression analysis is the runDiffExp function.
DESeq2.length.createRmd(
data.path,
result.path,
codefile,
fit.type,
test,
beta.prior = TRUE,
independent.filtering = TRUE,
cooks.cutoff = TRUE,
impute.outliers = TRUE,
extra.design.covariates = NULL,
nas.as.ones = FALSE
)The path to a .rds file containing the phyloCompData object that will be used for the differential expression analysis.
The path to the file where the result object will be saved.
The path to the file where the code will be written.
The fitting method used to get the dispersion-mean relationship. Possible values are "parametric", "local" and "mean".
The test to use. Possible values are "Wald" and "LRT".
Whether or not to put a zero-mean normal prior on the non-intercept coefficients. Default is TRUE.
Whether or not to perform independent filtering of the data. With independent filtering=TRUE, the adjusted p-values for genes not passing the filter threshold are set to NA.
The cutoff value for the Cook's distance to consider a value to be an outlier. Set to Inf or FALSE to disable outlier detection. For genes with detected outliers, the p-value and adjusted p-value will be set to NA.
Whether or not the outliers should be replaced by a trimmed mean and the analysis rerun.
A vector containing the names of extra control variables to be passed to the design matrix of DESeq2. All the covariates need to be a column of the sample.annotations data frame from the phyloCompData object, with a matching column name. The covariates can be a numeric vector, or a factor. Note that "condition" factor column is always included, and should not be added here. See Details.
Whether or not adjusted p values that are returned as NA by DESeq2 should be set to 1. This option is useful for comparisons with other methods. For more details, see section "I want to benchmark DESeq2 comparing to other DE tools" from the DESeq2 vignette (available by running vignette("DESeq2", package = "DESeq2")). Default to FALSE.
The function generates a .Rmd file containing the code for performing the differential expression analysis. This file can be executed using e.g. the knitr package.
For more information about the methods and the interpretation of the parameters, see the DESeq2 package and the corresponding publications.
The lengths matrix is used as a normalization factor and applied to the DESeq2
model in the way explained in normalizationFactors
(see examples of this function).
The provided matrix will be multiplied by the default normalization factor
obtained through the estimateSizeFactors function.
The design model used in the DESeqDataSetFromMatrix
uses the "condition" column of the sample.annotations data frame from the phyloCompData object
as well as all the covariates named in extra.design.covariates.
For example, if extra.design.covariates = c("var1", "var2"), then
sample.annotations must have two columns named "var1" and "var2", and the design formula
in the DESeqDataSetFromMatrix function will be:
~ condition + var1 + var2.
Anders S and Huber W (2010): Differential expression analysis for sequence count data. Genome Biology 11:R106
Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8.
try(
if (require(DESeq2)) {
tmpdir <- normalizePath(tempdir(), winslash = "/")
## Simulate data
mydata.obj <- generateSyntheticData(dataset = "mydata", n.vars = 1000,
samples.per.cond = 5, n.diffexp = 100,
id.species = 1:10,
lengths.relmeans = rpois(1000, 1000),
lengths.dispersions = rgamma(1000, 1, 1),
output.file = file.path(tmpdir, "mydata.rds"))
## Add covariates
## Model fitted is count.matrix ~ condition + test_factor + test_reg
sample.annotations(mydata.obj)$test_factor <- factor(rep(1:2, each = 5))
sample.annotations(mydata.obj)$test_reg <- rnorm(10, 0, 1)
saveRDS(mydata.obj, file.path(tmpdir, "mydata.rds"))
## Diff Exp
runDiffExp(data.file = file.path(tmpdir, "mydata.rds"), result.extent = "DESeq2",
Rmdfunction = "DESeq2.length.createRmd",
output.directory = tmpdir, fit.type = "parametric",
test = "Wald",
extra.design.covariates = c("test_factor", "test_reg"))
})
#> Warning: Vector 'id.species' must be a factor. Transforming.
#> Error in sample.annotations(mydata.obj)$test_factor <- factor(rep(1:2, :
#> could not find function "sample.annotations"