Confident fold change
2021-10-26
Contents
This document shows typical Topconfects usage with limma, edgeR, or DESeq2.
The first step is to load a dataset. Here, we’re looking at RNA-seq data that
investigates the response of Arabodopsis thaliana to a bacterial pathogen.
Besides the experimental and control conditions, there is also a batch effect.
This dataset is also examined in section 4.2 of the edgeR user manual, and
I’ve followed the initial filtering steps in the edgeR manual.
library(topconfects)
library(NBPSeq)
library(edgeR)
library(limma)
library(dplyr)
library(ggplot2)
data(arab)
# Retrieve symbol for each gene
info <-
AnnotationDbi::select(
org.At.tair.db::org.At.tair.db,
rownames(arab), "SYMBOL") %>%
group_by(TAIR) %>%
summarize(
SYMBOL=paste(unique(na.omit(SYMBOL)),collapse="/"))
arab_info <-
info[match(rownames(arab),info$TAIR),] %>%
select(-TAIR) %>%
as.data.frame
rownames(arab_info) <- rownames(arab)
# Extract experimental design from sample names
Treat <- factor(substring(colnames(arab),1,4)) %>% relevel(ref="mock")
Time <- factor(substring(colnames(arab),5,5))
y <- DGEList(arab, genes=as.data.frame(arab_info))
# Keep genes with at least 3 samples having an RPM of more than 2
keep <- rowSums(cpm(y)>2) >= 3
y <- y[keep,,keep.lib.sizes=FALSE]
y <- calcNormFactors(y)1 limma analysis
1.1 Standard limma analysis steps
design <- model.matrix(~Time+Treat)
design[,]## (Intercept) Time2 Time3 Treathrcc
## 1 1 0 0 0
## 2 1 1 0 0
## 3 1 0 1 0
## 4 1 0 0 1
## 5 1 1 0 1
## 6 1 0 1 1fit <-
voom(y, design) %>%
lmFit(design)1.2 Apply topconfects
Find largest fold changes that we can be confident of at FDR 0.05.
confects <- limma_confects(fit, coef="Treathrcc", fdr=0.05)
confects## $table
## confect effect AveExpr name SYMBOL
## 1 3.090 4.849 6.567 AT3G46280
## 2 2.895 4.331 9.155 AT4G12500
## 3 2.895 4.493 6.053 AT2G19190 FRK1/SIRK
## 4 2.608 3.839 9.166 AT4G12490 AZI3
## 5 2.608 3.952 6.615 AT1G51800 IOS1
## 6 2.608 4.299 5.440 AT2G39530 CASPL4D1
## 7 2.606 3.908 7.899 AT2G39200 ATMLO12/MLO12
## 8 2.606 3.710 8.729 AT5G64120 AtPRX71/PRX71
## 9 2.595 5.346 1.807 AT5G31702
## 10 2.563 4.888 2.383 AT2G08986
## ...
## 2184 of 16526 non-zero log2 fold change at FDR 0.05
## Prior df 18.21.3 Looking at the result
Here the usual logFC values estimated by limma are shown as dots, with lines
to the confect value.
confects_plot(confects)
confects_plot_me overlays the confects (black) on a Mean-Difference Plot
(grey) (as might be produced by limma::plotMD). As we should expect, the very
noisy differences with low mean expression are removed if we look at the
confects.
confects_plot_me(confects)
Let’s compare this to the ranking we obtain from topTable.
fit_eb <- eBayes(fit)
top <- topTable(fit_eb, coef="Treathrcc", n=Inf)
rank_rank_plot(confects$table$name, rownames(top), "limma_confects", "topTable")
You can see that the top 19 genes from topTable are all within the top 40 for topconfects ranking, but topconfects has also highly ranked some other genes. These have a large effect size, and sufficient if not overwhelming evidence of this.
An MD-plot highlighting the positions of the top 40 genes in both rankings also illustrates the differences between these two ways of ranking genes.
plotMD(fit, legend="bottomleft", status=paste0(
ifelse(rownames(fit) %in% rownames(top)[1:40], "topTable ",""),
ifelse(rownames(fit) %in% confects$table$name[1:40], "confects ","")))
2 edgeR analysis
An analysis in edgeR produces similar results. Note that only quasi-likelihood testing from edgeR is supported.
2.1 Standard edgeR analysis
y <- estimateDisp(y, design, robust=TRUE)
efit <- glmQLFit(y, design, robust=TRUE)2.2 Apply topconfects
A step of 0.05 is used here merely so that the vignette will build quickly.
edger_confects calls edgeR::glmTreat repeatedly, which is necessarily slow.
In practice a smaller value such as 0.01 should be used.
econfects <- edger_confects(efit, coef="Treathrcc", fdr=0.05, step=0.05)
econfects## $table
## confect effect logCPM name SYMBOL
## 1 3.05 6.319 6.729 AT5G48430
## 2 2.95 4.785 8.097 AT3G46280
## 3 2.95 5.779 4.903 AT3G55150 ATEXO70H1/EXO70H1
## 4 2.95 4.497 7.374 AT2G19190 FRK1/SIRK
## 5 2.95 4.943 5.771 AT2G39380 ATEXO70H2/EXO70H2
## 6 2.95 5.321 5.419 AT1G51850 SIF2
## 7 2.95 5.414 5.199 AT2G44370
## 8 2.85 4.341 6.709 AT2G39530 CASPL4D1
## 9 2.55 4.337 6.373 AT1G51820 SIF4
## 10 2.55 4.527 5.587 AT2G30750 CYP71A12
## ...
## 1729 of 16526 non-zero log2 fold change at FDR 0.05
## Dispersion 0.045 to 0.31
## Biological CV 21.1% to 55.8%2.3 Looking at the result
confects_plot(econfects)
confects_plot_me(econfects)
etop <-
glmQLFTest(efit, coef="Treathrcc") %>%
topTags(n=Inf)
plotMD(efit, legend="bottomleft", status=paste0(
ifelse(rownames(efit) %in% econfects$table$name[1:40], "confects ", ""),
ifelse(rownames(efit) %in% rownames(etop)[1:40], "topTags ","")))
3 DESeq2 analysis
DESeq2 does its own filtering of lowly expressed genes, so we start from the original count matrix. The initial steps are as for a normal DESeq2 analysis.
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
countData = arab,
colData = data.frame(Time, Treat),
rowData = arab_info,
design = ~Time+Treat)
dds <- DESeq(dds)3.1 Apply topconfects
The contrast or coefficient to test is specified as in the DESeq2::results
function. The step of 0.05 is merely so that this vignette will build quickly,
in practice a smaller value such as 0.01 should be used. deseq2_confects
calls results repeatedly, and in fairness results
has not been optimized for this.
dconfects <- deseq2_confects(dds, name="Treat_hrcc_vs_mock", step=0.05)DESeq2 offers shrunken estimates of LFC. This is another sensible way of ranking genes. Let’s compare them to the confect values.
shrunk <- lfcShrink(dds, coef="Treat_hrcc_vs_mock", type="ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041dconfects$table$shrunk <- shrunk$log2FoldChange[dconfects$table$index]
dconfects## $table
## confect effect baseMean name filtered shrunk
## 1 3.40 4.785 586.43 AT3G46280 FALSE 4.637
## 2 3.15 4.501 354.89 AT2G19190 FALSE 4.363
## 3 3.10 4.320 2997.60 AT4G12500 FALSE 4.211
## 4 2.90 4.976 115.20 AT2G39380 FALSE 4.595
## 5 2.90 5.433 89.39 AT1G51850 FALSE 4.829
## 6 2.90 4.350 223.26 AT2G39530 FALSE 4.176
## 7 2.85 5.453 77.13 AT2G44370 FALSE 4.785
## 8 2.85 5.875 62.88 AT3G55150 FALSE 4.930
## 9 2.85 3.838 2532.56 AT4G12490 FALSE 3.763
## 10 2.80 3.969 448.55 AT1G51800 FALSE 3.864
## ...
## 1219 of 26222 non-zero effect size at FDR 0.05DESeq2 filters some genes, these are placed last in the table. If your intention
is to obtain a ranking of all genes, you should disable this with
deseq2_confects(..., cooksCutoff=Inf, independentFiltering=FALSE).
table(dconfects$table$filtered)##
## FALSE TRUE
## 11479 14743tail(dconfects$table)## rank index confect effect baseMean name filtered shrunk
## 26217 26217 26203 NA 0.8487098 6.324289 AT5G67460 TRUE 0.08464016
## 26218 26218 26206 NA -0.6852089 1.220362 AT5G67488 TRUE -0.01672116
## 26219 26219 26207 NA 1.4505333 4.657085 AT5G67490 TRUE 0.13703515
## 26220 26220 26210 NA -0.5940465 6.699483 AT5G67520 TRUE -0.06727656
## 26221 26221 26213 NA 0.6767779 1.542612 AT5G67550 TRUE 0.02225821
## 26222 26222 26219 NA 0.1805069 12.111601 AT5G67610 TRUE 0.030217573.2 Looking at the result
Shrunk LFC estimates are shown in red.
confects_plot(dconfects) +
geom_point(aes(x=shrunk, size=baseMean, color="lfcShrink"), alpha=0.75)
lfcShrink aims for a best estimate of the LFC, whereas confect is a
conservative estimate. lfcShrink can produce non-zero values for genes which
can’t be said to significantly differ from zero – it doesn’t do double duty as
an indication of significance – whereas the confect value will be NA in this
case. The plot below compares these two quantities. Only un-filtered genes are
shown (see above).
filter(dconfects$table, !filtered) %>%
ggplot(aes(
x=ifelse(is.na(confect),0,confect), y=shrunk, color=!is.na(confect))) +
geom_point() + geom_abline() + coord_fixed() + theme_bw() +
labs(color="Significantly\nnon-zero at\nFDR 0.05",
x="confect", y="lfcShrink using ashr")
4 Comparing results
rank_rank_plot(confects$table$name, econfects$table$name,
"limma confects", "edgeR confects")
rank_rank_plot(confects$table$name, dconfects$table$name,
"limma confects", "DESeq2 confects")
rank_rank_plot(econfects$table$name, dconfects$table$name,
"edgeR confects", "DESeq2 confects")
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DESeq2_1.34.0 SummarizedExperiment_1.24.0
## [3] Biobase_2.54.0 MatrixGenerics_1.6.0
## [5] matrixStats_0.61.0 GenomicRanges_1.46.0
## [7] GenomeInfoDb_1.30.0 IRanges_2.28.0
## [9] S4Vectors_0.32.0 BiocGenerics_0.40.0
## [11] ggplot2_3.3.5 dplyr_1.0.7
## [13] edgeR_3.36.0 limma_3.50.0
## [15] NBPSeq_0.3.0 topconfects_1.10.0
## [17] BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 bit64_4.0.5 RColorBrewer_1.1-2
## [4] httr_1.4.2 tools_4.1.1 bslib_0.3.1
## [7] irlba_2.3.3 utf8_1.2.2 R6_2.5.1
## [10] DBI_1.1.1 colorspace_2.0-2 withr_2.4.2
## [13] tidyselect_1.1.1 bit_4.0.4 compiler_4.1.1
## [16] DelayedArray_0.20.0 labeling_0.4.2 bookdown_0.24
## [19] sass_0.4.0 scales_1.1.1 SQUAREM_2021.1
## [22] genefilter_1.76.0 mixsqp_0.3-43 stringr_1.4.0
## [25] digest_0.6.28 rmarkdown_2.11 XVector_0.34.0
## [28] pkgconfig_2.0.3 htmltools_0.5.2 invgamma_1.1
## [31] fastmap_1.1.0 highr_0.9 org.At.tair.db_3.14.0
## [34] rlang_0.4.12 RSQLite_2.2.8 jquerylib_0.1.4
## [37] farver_2.1.0 generics_0.1.1 jsonlite_1.7.2
## [40] BiocParallel_1.28.0 RCurl_1.98-1.5 magrittr_2.0.1
## [43] GenomeInfoDbData_1.2.7 Matrix_1.3-4 Rcpp_1.0.7
## [46] munsell_0.5.0 fansi_0.5.0 lifecycle_1.0.1
## [49] stringi_1.7.5 yaml_2.2.1 zlibbioc_1.40.0
## [52] plyr_1.8.6 qvalue_2.26.0 grid_4.1.1
## [55] blob_1.2.2 parallel_4.1.1 crayon_1.4.1
## [58] lattice_0.20-45 Biostrings_2.62.0 splines_4.1.1
## [61] annotate_1.72.0 KEGGREST_1.34.0 locfit_1.5-9.4
## [64] magick_2.7.3 knitr_1.36 pillar_1.6.4
## [67] geneplotter_1.72.0 reshape2_1.4.4 XML_3.99-0.8
## [70] glue_1.4.2 evaluate_0.14 BiocManager_1.30.16
## [73] png_0.1-7 vctrs_0.3.8 gtable_0.3.0
## [76] purrr_0.3.4 assertthat_0.2.1 ashr_2.2-47
## [79] cachem_1.0.6 xfun_0.27 xtable_1.8-4
## [82] survival_3.2-13 truncnorm_1.0-8 tibble_3.1.5
## [85] AnnotationDbi_1.56.0 memoise_2.0.0 statmod_1.4.36
## [88] ellipsis_0.3.2