For all possible absolute log2 fold changes, which genes have at least this fold change at a specified False Discovery Rate? This is built by repeatedly calling DESeq2::results with the "greaterAbs" alternative hypothesis.

deseq2_confects(object, ..., fdr = 0.05, step = 0.01)

Arguments

object

Object produced by the DESeq2::DESeq function.

...

Further arguments to DESeq2::results. At a minimum you should specify either contrast= or name=.

fdr

False Discovery Rate to control for.

step

Granularity of log2 fold changes to test.

Value

See nest_confects for details of how to interpret the result.

The filtered column in the result indicates whether DESeq2 filtered the gene. Such genes do not count toward the total number of genes when controlling FDR. If your intention is to obtain a ranking of all genes, you should disable this with deseq2_confects(..., cooksCutoff=Inf, independentFiltering=FALSE).

Details

Results are presented in a table such that for any given LFC, if the reader chooses the genes with abs(confect) less than this they are assured that this set of genes has at least this LFC (with the specified FDR). The confect column may also be viewed as a confidence bound on the LFC of each gene, with a dynamic correction for multiple testing.

Examples

# Generate some random data n <- 20 folds <- seq(-8,8,length.out=n) row_means <- runif(n, min=0, max=5) lib_scale <- c(1,2,3,4) means <- 2^(outer(folds, c(-0.5,-0.5,0.5,0.5))) * row_means * rep(lib_scale,each=n) counts <- rnbinom(length(means), mu=means, size=1/0.1) dim(counts) <- dim(means) group <- factor(c("A","A","B","B")) # Apply DESeq2 library(DESeq2)
#> Loading required package: S4Vectors
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: parallel
#> #> Attaching package: ‘BiocGenerics’
#> The following objects are masked from ‘package:parallel’: #> #> clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, #> clusterExport, clusterMap, parApply, parCapply, parLapply, #> parLapplyLB, parRapply, parSapply, parSapplyLB
#> The following object is masked from ‘package:limma’: #> #> plotMA
#> The following objects are masked from ‘package:stats’: #> #> IQR, mad, sd, var, xtabs
#> The following objects are masked from ‘package:base’: #> #> anyDuplicated, append, as.data.frame, basename, cbind, colMeans, #> colnames, colSums, 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, rowMeans, rownames, #> rowSums, sapply, setdiff, sort, table, tapply, union, unique, #> unsplit, which, which.max, which.min
#> #> Attaching package: ‘S4Vectors’
#> The following object is masked from ‘package:base’: #> #> expand.grid
#> Loading required package: IRanges
#> Loading required package: GenomicRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: SummarizedExperiment
#> 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")'.
#> Loading required package: DelayedArray
#> Loading required package: matrixStats
#> #> Attaching package: ‘matrixStats’
#> The following objects are masked from ‘package:Biobase’: #> #> anyMissing, rowMedians
#> Loading required package: BiocParallel
#> #> Attaching package: ‘DelayedArray’
#> The following objects are masked from ‘package:matrixStats’: #> #> colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
#> The following objects are masked from ‘package:base’: #> #> aperm, apply
dds <- DESeqDataSetFromMatrix( countData = counts, colData = data.frame(group=group), design = ~group)
#> converting counts to integer mode
dds <- DESeq(dds)
#> 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.
#> Warning: Estimated rdf < 1.0; not estimating variance
#> final dispersion estimates
#> fitting model and testing
# Find top confident effect sizes deseq2_confects(dds, name="group_B_vs_A", step=0.1)
#> $table #> rank index confect effect baseMean filtered #> 1 20 2.5 8.788926 114.938015 FALSE #> 2 16 1.9 7.614297 54.555396 FALSE #> 3 18 1.9 6.611207 52.719543 FALSE #> 4 19 1.6 6.867433 30.503372 FALSE #> 5 2 -1.5 -6.570242 21.891665 FALSE #> 6 17 1.5 5.830777 30.928070 FALSE #> 7 4 -0.6 -5.572645 11.088803 FALSE #> 8 3 0.0 -3.733286 18.576539 FALSE #> 9 13 0.0 4.031585 9.281007 FALSE #> 10 6 NA -4.744159 2.298399 FALSE #> ... #> 9 of 20 non-zero effect size at FDR 0.05