Experimental function for performing 2x3 factor DESeq2 analyses. Output can
be passed to deseq_2x3_polar() and subsequently plotted. Example usage
would include comparing gene expression against a binary outcome e.g.
response vs non-response, across 3 drugs: the design would be ~ response
and group would refer to the medication column in the metadata.
deseq_2x3(object, design, group, ...)
| object | An object of class 'DESeqDataSet' containing full dataset |
|---|---|
| design | Design formula. The main contrast is taken from the last term of the formula and must be a binary factor. |
| group | Character value for the column with the 3-way grouping factor
within the sample information data |
| ... | Optional arguments passed to |
Returns a list of 3 DESeq2 results objects which can be passed onto
deseq_2x3_polar()
# \donttest{ # Basic DESeq2 set up library(DESeq2) #> Loading required package: S4Vectors #> Loading required package: stats4 #> Loading required package: BiocGenerics #> #> Attaching package: ‘BiocGenerics’ #> The following objects are masked from ‘package:stats’: #> #> IQR, mad, sd, var, xtabs #> The following objects are masked from ‘package:base’: #> #> anyDuplicated, aperm, append, as.data.frame, basename, cbind, #> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, #> get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, #> match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, #> Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort, #> table, tapply, union, unique, unsplit, which.max, which.min #> #> Attaching package: ‘S4Vectors’ #> The following objects are masked from ‘package:base’: #> #> expand.grid, I, unname #> Loading required package: IRanges #> Loading required package: GenomicRanges #> Loading required package: GenomeInfoDb #> Loading required package: SummarizedExperiment #> Loading required package: MatrixGenerics #> Loading required package: matrixStats #> #> Attaching package: ‘MatrixGenerics’ #> The following objects are masked from ‘package:matrixStats’: #> #> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, #> colCounts, colCummaxs, colCummins, colCumprods, colCumsums, #> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, #> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, #> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, #> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, #> colWeightedMeans, colWeightedMedians, colWeightedSds, #> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, #> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, #> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, #> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, #> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, #> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, #> rowWeightedMads, rowWeightedMeans, rowWeightedMedians, #> rowWeightedSds, rowWeightedVars #> Loading required package: Biobase #> Welcome to Bioconductor #> #> Vignettes contain introductory material; view with #> 'browseVignettes()'. To cite Bioconductor, see #> 'citation("Biobase")', and for packages 'citation("pkgname")'. #> #> Attaching package: ‘Biobase’ #> The following object is masked from ‘package:MatrixGenerics’: #> #> rowMedians #> The following objects are masked from ‘package:matrixStats’: #> #> anyMissing, rowMedians counts <- matrix(rnbinom(n=3000, mu=100, size=1/0.5), ncol=30) rownames(counts) <- paste0("gene", 1:100) cond <- rep(factor(rep(1:3, each=5), labels = c('A', 'B', 'C')), 2) resp <- factor(rep(1:2, each=15), labels = c('non.responder', 'responder')) metadata <- data.frame(drug = cond, response = resp) # Full dataset object construction dds <- DESeqDataSetFromMatrix(counts, metadata, ~response) #> converting counts to integer mode # Perform 3x DESeq2 analyses comparing binary response for each drug res <- deseq_2x3(dds, ~response, "drug") #> drug = A #> 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 #> drug = B #> 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 #> drug = C #> 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 # Generate polar object obj <- deseq_2x3_polar(res) # 2d plot radial_plotly(obj) # 3d plot volcano3D(obj) # }