This document shows typical Topconfects usage.

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)

limma analysis

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         1
fit <-
    voom(y, design) %>%
    lmFit(design) %>%
    eBayes()

fit$df.prior
## [1] 18.24598

Apply topconfects

Find largest fold changes that we can be confident of at FDR 0.05.

confects <- limma_confects(fit, coef=4, fdr=0.05)

confects
## $table
##  rank index confect effect   AveExpr  name      SYMBOL       
##   1    9023 3.09    4.849117 6.566691 AT3G46280              
##   2   10774 2.89    4.330654 9.155218 AT4G12500              
##   3    4945 2.89    4.493384 6.052661 AT2G19190 FRK1         
##   4   10773 2.60    3.838721 9.165703 AT4G12490              
##   5    2562 2.60    3.951678 6.615271 AT1G51800 IOS1         
##   6    6264 2.60    4.299425 5.440243 AT2G39530              
##   7    6239 2.60    3.907806 7.898640 AT2G39200 ATMLO12/MLO12
##   8   16247 2.60    3.710435 8.728518 AT5G64120 AtPRX71/PRX71
##   9   14449 2.59    5.346437 1.806933 AT5G31702              
##  10    4642 2.56    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

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 (white) (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.

top <- topTable(fit, 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.

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 ","")))

edgeR analysis

An analysis in edgeR produces similar results.

Standard edgeR analysis

y <- estimateDisp(y, design, robust=TRUE)
efit <- glmQLFit(y, design, robust=TRUE)

Apply topconfects

econfects <- edger_confects(efit, coef=4, fdr=0.05)

econfects
## $table
##  rank index confect effect   logCPM   name      SYMBOL           
##   1   15207 3.06    6.319297 6.729025 AT5G48430                  
##   2    9023 2.98    4.784830 8.096900 AT3G46280                  
##   3    9617 2.98    5.778590 4.902951 AT3G55150 ATEXO70H1/EXO70H1
##   4    4945 2.98    4.497055 7.373501 AT2G19190 FRK1             
##   5    6253 2.98    4.942601 5.771324 AT2G39380 ATEXO70H2/EXO70H2
##   6    2565 2.98    5.321291 5.418659 AT1G51850                  
##   7    6619 2.98    5.414120 5.199481 AT2G44370                  
##   8    6264 2.85    4.340932 6.709049 AT2G39530                  
##   9    2564 2.58    4.336766 6.373316 AT1G51820                  
##  10    5662 2.57    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%

Looking at the result

confects_plot(econfects)

etop <-
    glmQLFTest(efit, coef=4) %>%
    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 ","")))

Comparing limma and edgeR results

rank_rank_plot(confects$table$name, econfects$table$name, "limma confects", "edgeR confects")


sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] bindrcpp_0.2      ggplot2_2.2.1     dplyr_0.7.4       edgeR_3.20.1     
## [5] limma_3.34.0      NBPSeq_0.3.0      topconfects_0.0.1
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.13         compiler_3.4.2       plyr_1.8.4          
##  [4] bindr_0.1            tools_3.4.2          statmod_1.4.30      
##  [7] bit_1.1-12           digest_0.6.12        RSQLite_2.0         
## [10] evaluate_0.10.1      memoise_1.1.0        tibble_1.3.4        
## [13] gtable_0.2.0         lattice_0.20-35      pkgconfig_2.0.1     
## [16] rlang_0.1.4          DBI_0.7              parallel_3.4.2      
## [19] yaml_2.1.14          org.At.tair.db_3.4.2 stringr_1.2.0       
## [22] knitr_1.17           IRanges_2.12.0       S4Vectors_0.16.0    
## [25] bit64_0.9-7          stats4_3.4.2         locfit_1.5-9.1      
## [28] rprojroot_1.2        grid_3.4.2           Biobase_2.38.0      
## [31] glue_1.2.0           qvalue_2.10.0        R6_2.2.2            
## [34] AnnotationDbi_1.40.0 rmarkdown_1.6        blob_1.1.0          
## [37] purrr_0.2.4          reshape2_1.4.2       magrittr_1.5        
## [40] BiocGenerics_0.24.0  backports_1.1.1      scales_0.5.0        
## [43] htmltools_0.3.6      splines_3.4.2        assertthat_0.2.0    
## [46] colorspace_1.3-2     labeling_0.3         stringi_1.1.5       
## [49] lazyeval_0.2.1       munsell_0.4.3