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 and description for each gene
info <-
AnnotationDbi::select(
org.At.tair.db::org.At.tair.db,
rownames(arab), c("SYMBOL","GENENAME")) %>%
group_by(TAIR) %>%
summarize(
SYMBOL=paste(unique(na.omit(SYMBOL)),collapse="/"),
GENENAME=paste(unique(na.omit(GENENAME)),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)
## (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 1
Find largest fold changes that we can be confident of at FDR 0.05.
## $table
## rank index confect effect AveExpr name SYMBOL
## 1 9023 3.090 4.849117 6.566691 AT3G46280
## 2 10774 2.895 4.330654 9.155218 AT4G12500
## 3 4945 2.895 4.493384 6.052661 AT2G19190 FRK1
## 4 10773 2.608 3.838721 9.165703 AT4G12490
## 5 2562 2.608 3.951678 6.615271 AT1G51800 IOS1
## 6 6264 2.608 4.299425 5.440243 AT2G39530
## 7 6239 2.606 3.907806 7.898640 AT2G39200 ATMLO12/MLO12
## 8 16247 2.606 3.710435 8.728518 AT5G64120 AtPRX71/PRX71
## 9 14449 2.595 5.346437 1.806933 AT5G31702
## 10 4642 2.563 4.888477 2.382601 AT2G08986
## GENENAME
##
##
## Receptor-like protein kinase. Involved in early defense signaling.
##
##
##
## A member of a large family of seven-transmembrane domain proteins specific to plants, homologs of the barley mildew resistance locus o (MLO) protein. The Arabidopsis genome contains 15 genes encoding MLO proteins, with localization in plasma membrane. Phylogenetic analysis revealed four clades of closely-related AtMLO genes. ATMLO6 belongs to the clade IV, with AtMLO2, AtMLO3 and AtMLO12. The gene is expressed during early seedling growth, in root tips and cotyledon vascular system, in floral organs (anthers and stigma), and in fruit abscission zone, as shown by GUS activity patterns. The expression of several phylogenetically closely-related AtMLO genes showed similar or overlapping tissue specificity and analogous responsiveness to external stimuli, suggesting functional redundancy, co-function, or antagonistic function(s).
## encodes a cell wall bound peroxidase that is induced by hypo-osmolarity
##
##
## ...
## 2184 of 16526 non-zero log2 fold change at FDR 0.05
## Prior df 18.2
Here the usual logFC
values estimated by limma
are shown as dots, with lines to the confect
value.
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.
Let’s compare this to the ranking we obtain from topTable
.
fit_eb <- eBayes(fit)
top <- topTable(fit_eb, coef=4, 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.
An analysis in edgeR produces similar results. Note that only quasi-likelihood testing from edgeR is supported.
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.
## $table
## rank index confect effect logCPM name SYMBOL
## 1 15207 3.05 6.319297 6.729025 AT5G48430
## 2 9023 2.95 4.784830 8.096900 AT3G46280
## 3 9617 2.95 5.778590 4.902951 AT3G55150 ATEXO70H1/EXO70H1
## 4 4945 2.95 4.497055 7.373501 AT2G19190 FRK1
## 5 6253 2.95 4.942601 5.771324 AT2G39380 ATEXO70H2/EXO70H2
## 6 2565 2.95 5.321291 5.418659 AT1G51850
## 7 6619 2.95 5.414120 5.199481 AT2G44370
## 8 6264 2.85 4.340932 6.709049 AT2G39530
## 9 2564 2.55 4.336766 6.373316 AT1G51820
## 10 5662 2.55 4.526719 5.587339 AT2G30750 CYP71A12
## GENENAME
##
##
## A member of EXO70 gene family, putative exocyst subunits, conserved in land plants. Arabidopsis thaliana contains 23 putative EXO70 genes, which can be classified into eight clusters on the phylogenetic tree.
## Receptor-like protein kinase. Involved in early defense signaling.
## A member of EXO70 gene family, putative exocyst subunits, conserved in land plants. Arabidopsis thaliana contains 23 putative EXO70 genes, which can be classified into eight clusters on the phylogenetic tree.
##
##
##
##
## putative cytochrome P450
## ...
## 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%
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)
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.
DESeq2 offers shrunken estimates of LFC. This is another sensible way of ranking genes. Let’s compare them to the confect values.
## 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/kxw041
## $table
## rank index confect effect baseMean name filtered shrunk
## 1 14551 3.40 4.785076 586.43055 AT3G46280 FALSE 4.636524
## 2 8130 3.15 4.500949 354.89430 AT2G19190 FALSE 4.362663
## 3 17495 3.10 4.320279 2997.59547 AT4G12500 FALSE 4.210776
## 4 10088 2.90 4.976427 115.19576 AT2G39380 FALSE 4.595044
## 5 4063 2.90 5.433121 89.38966 AT1G51850 FALSE 4.829158
## 6 10106 2.90 4.350196 223.26235 AT2G39530 FALSE 4.176090
## 7 10602 2.85 5.453274 77.12549 AT2G44370 FALSE 4.784621
## 8 15411 2.85 5.875086 62.87522 AT3G55150 FALSE 4.930458
## 9 17494 2.85 3.838114 2532.56056 AT4G12490 FALSE 3.762826
## 10 4058 2.80 3.969472 448.54556 AT1G51800 FALSE 3.863933
## ...
## 1219 of 26222 non-zero effect size at FDR 0.05
DESeq2 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)
.
##
## FALSE TRUE
## 11479 14743
## rank index confect effect baseMean name filtered shrunk
## 26217 26217 26203 NA 0.8487098 6.324289 AT5G67460 TRUE 0.08464004
## 26218 26218 26206 NA -0.6852089 1.220362 AT5G67488 TRUE -0.01672112
## 26219 26219 26207 NA 1.4505333 4.657085 AT5G67490 TRUE 0.13703492
## 26220 26220 26210 NA -0.5940465 6.699483 AT5G67520 TRUE -0.06727648
## 26221 26221 26213 NA 0.6767779 1.542612 AT5G67550 TRUE 0.02225816
## 26222 26222 26219 NA 0.1805069 12.111601 AT5G67610 TRUE 0.03021755
Shrunk LFC estimates are shown in red.
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).
## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.11-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] DESeq2_1.28.0 SummarizedExperiment_1.18.0
## [3] DelayedArray_0.14.0 matrixStats_0.56.0
## [5] Biobase_2.48.0 GenomicRanges_1.40.0
## [7] GenomeInfoDb_1.24.0 IRanges_2.22.0
## [9] S4Vectors_0.26.0 BiocGenerics_0.34.0
## [11] ggplot2_3.3.0 dplyr_0.8.5
## [13] edgeR_3.30.0 limma_3.44.0
## [15] NBPSeq_0.3.0 topconfects_1.4.0
##
## loaded via a namespace (and not attached):
## [1] bit64_0.9-7 splines_4.0.0 assertthat_0.2.1
## [4] statmod_1.4.34 mixsqp_0.3-17 blob_1.2.1
## [7] GenomeInfoDbData_1.2.3 yaml_2.2.1 pillar_1.4.3
## [10] RSQLite_2.2.0 lattice_0.20-41 glue_1.4.0
## [13] digest_0.6.25 RColorBrewer_1.1-2 XVector_0.28.0
## [16] qvalue_2.20.0 colorspace_1.4-1 htmltools_0.4.0
## [19] Matrix_1.2-18 plyr_1.8.6 XML_3.99-0.3
## [22] pkgconfig_2.0.3 invgamma_1.1 genefilter_1.70.0
## [25] zlibbioc_1.34.0 purrr_0.3.4 xtable_1.8-4
## [28] scales_1.1.0 BiocParallel_1.22.0 tibble_3.0.1
## [31] annotate_1.66.0 farver_2.0.3 ellipsis_0.3.0
## [34] withr_2.2.0 ashr_2.2-47 survival_3.1-12
## [37] magrittr_1.5 crayon_1.3.4 memoise_1.1.0
## [40] evaluate_0.14 truncnorm_1.0-8 tools_4.0.0
## [43] lifecycle_0.2.0 stringr_1.4.0 munsell_0.5.0
## [46] locfit_1.5-9.4 irlba_2.3.3 AnnotationDbi_1.50.0
## [49] compiler_4.0.0 rlang_0.4.5 grid_4.0.0
## [52] RCurl_1.98-1.2 bitops_1.0-6 labeling_0.3
## [55] rmarkdown_2.1 gtable_0.3.0 DBI_1.1.0
## [58] reshape2_1.4.4 R6_2.4.1 knitr_1.28
## [61] bit_1.1-15.2 org.At.tair.db_3.10.1 stringi_1.4.6
## [64] SQUAREM_2020.2 Rcpp_1.0.4.6 vctrs_0.2.4
## [67] geneplotter_1.66.0 tidyselect_1.0.0 xfun_0.13