It appears that it is possible though somewhat difficult to apply batch estimations generated by sva to the model given to DESeq/EdgeR/limma. In the case of limma it is fairly simple, but in the other two it is a bit more difficult. There is a nice discussion of this at: https://www.biostars.org/p/156186/ I am attempting to apply that logic to this data with limited success.
keepers <- list(
"wt_stex" = c("wt_st", "wt_ex"),
"mt_stex" = c("mt_st", "mt_ex"),
"st_wtmt" = c("mt_st", "wt_st"),
"ex_wtmt" = c("mt_ex", "wt_ex"))
pa_de <- all_pairwise(input=pa_expt,
model_batch=FALSE,
limma_method="robust",
keepers=keepers,
combined_excel=paste0("excel/pa_nobatch-v", ver, ".xlsx"))## Assuming no batch in model for testing pca.## Finished running DE analyses, collecting outputs.## Comparing analyses 1/15: mt_st_vs_mt_ex## Comparing analyses 2/15: mt_undef_vs_mt_ex## Comparing analyses 3/15: wt_ex_vs_mt_ex## Comparing analyses 4/15: wt_st_vs_mt_ex## Comparing analyses 5/15: wt_undef_vs_mt_ex## Comparing analyses 6/15: mt_undef_vs_mt_st## Comparing analyses 7/15: wt_ex_vs_mt_st## Comparing analyses 8/15: wt_st_vs_mt_st## Comparing analyses 9/15: wt_undef_vs_mt_st## Comparing analyses 10/15: wt_ex_vs_mt_undef## Comparing analyses 11/15: wt_st_vs_mt_undef## Comparing analyses 12/15: wt_undef_vs_mt_undef## Comparing analyses 13/15: wt_st_vs_wt_ex## Comparing analyses 14/15: wt_undef_vs_wt_ex## Comparing analyses 15/15: wt_undef_vs_wt_st## Invoking combine_de_tables().## Deleting the file excel/pa_nobatch-v20171019.xlsx before writing the tables.## Writing a legend of columns.## Printing a pca plot before/after surrogates/batch estimation.## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse## Working on 1/4: wt_stex## Found table with wt_st_vs_wt_ex## Working on 2/4: mt_stex## Found table with mt_st_vs_mt_ex## Working on 3/4: st_wtmt## Found inverse table with wt_st_vs_mt_st## Working on 4/4: ex_wtmt## Found inverse table with wt_ex_vs_mt_ex## Adding venn plots for wt_st_vs_wt_ex.## Limma expression coefficients for wt_st_vs_wt_ex; R^2: 0.67; equation: y = 0.811x + 0.815## Edger expression coefficients for wt_st_vs_wt_ex; R^2: 0.681; equation: y = 0.743x + 1.92## DESeq2 expression coefficients for wt_st_vs_wt_ex; R^2: 0.673; equation: y = 0.735x + 2.11## Adding venn plots for mt_st_vs_mt_ex.## Limma expression coefficients for mt_st_vs_mt_ex; R^2: 0.915; equation: y = 0.952x + 0.145## Edger expression coefficients for mt_st_vs_mt_ex; R^2: 0.911; equation: y = 0.849x + 1.15## DESeq2 expression coefficients for mt_st_vs_mt_ex; R^2: 0.91; equation: y = 0.845x + 1.23## Adding venn plots for mt_st_vs_wt_st.## Limma expression coefficients for mt_st_vs_wt_st; R^2: 0.796; equation: y = 0.903x + 0.55## Edger expression coefficients for mt_st_vs_wt_st; R^2: 0.795; equation: y = 0.968x + 0.23## DESeq2 expression coefficients for mt_st_vs_wt_st; R^2: 0.786; equation: y = 0.965x + 0.21## Adding venn plots for mt_ex_vs_wt_ex.## Limma expression coefficients for mt_ex_vs_wt_ex; R^2: 0.943; equation: y = 0.975x + 0.133## Edger expression coefficients for mt_ex_vs_wt_ex; R^2: 0.946; equation: y = 1.06x - 0.518## DESeq2 expression coefficients for mt_ex_vs_wt_ex; R^2: 0.944; equation: y = 1.07x - 0.565## Writing summary information.## Attempting to add the comparison plot to pairwise_summary at row: 22 and column: 1##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%## Performing save of the workbook.pa_sig <- extract_significant_genes(pa_de$combined, lfc=3,
excel=paste0("excel/pa_nobatch_significant-v", ver, ".xlsx"),
according_to="all")## Writing a legend of columns.## Writing excel data for wt_st_vs_wt_ex: 1/16.## After (adj)p filter, the up genes table has 1906 genes.## After (adj)p filter, the down genes table has 2509 genes.## After fold change filter, the up genes table has 332 genes.## After fold change filter, the down genes table has 124 genes.## Writing excel data for mt_st_vs_mt_ex: 2/16.## After (adj)p filter, the up genes table has 1242 genes.## After (adj)p filter, the down genes table has 1720 genes.## After fold change filter, the up genes table has 58 genes.## After fold change filter, the down genes table has 7 genes.## Writing excel data for mt_st_vs_wt_st: 3/16.## After (adj)p filter, the up genes table has 2092 genes.## After (adj)p filter, the down genes table has 1696 genes.## After fold change filter, the up genes table has 81 genes.## After fold change filter, the down genes table has 189 genes.## Writing excel data for mt_ex_vs_wt_ex: 4/16.## After (adj)p filter, the up genes table has 1085 genes.## After (adj)p filter, the down genes table has 1071 genes.## After fold change filter, the up genes table has 39 genes.## After fold change filter, the down genes table has 20 genes.## Printing significant genes to the file: excel/pa_nobatch_significant-v20171019.xlsx## 1/4: Creating significant table up_1limma_wt_st_vs_wt_ex## 2/4: Creating significant table up_2limma_mt_st_vs_mt_ex## 3/4: Creating significant table up_3limma_mt_st_vs_wt_st## 4/4: Creating significant table up_4limma_mt_ex_vs_wt_ex## Writing excel data for wt_st_vs_wt_ex: 5/16.## After (adj)p filter, the up genes table has 1637 genes.## After (adj)p filter, the down genes table has 1728 genes.## After fold change filter, the up genes table has 307 genes.## After fold change filter, the down genes table has 107 genes.## Writing excel data for mt_st_vs_mt_ex: 6/16.## After (adj)p filter, the up genes table has 931 genes.## After (adj)p filter, the down genes table has 969 genes.## After fold change filter, the up genes table has 71 genes.## After fold change filter, the down genes table has 6 genes.## Writing excel data for mt_st_vs_wt_st: 7/16.## After (adj)p filter, the up genes table has 1310 genes.## After (adj)p filter, the down genes table has 1265 genes.## After fold change filter, the up genes table has 66 genes.## After fold change filter, the down genes table has 174 genes.## Writing excel data for mt_ex_vs_wt_ex: 8/16.## After (adj)p filter, the up genes table has 386 genes.## After (adj)p filter, the down genes table has 489 genes.## After fold change filter, the up genes table has 41 genes.## After fold change filter, the down genes table has 25 genes.## Printing significant genes to the file: excel/pa_nobatch_significant-v20171019.xlsx## 1/4: Creating significant table up_1edger_wt_st_vs_wt_ex## 2/4: Creating significant table up_2edger_mt_st_vs_mt_ex## 3/4: Creating significant table up_3edger_mt_st_vs_wt_st## 4/4: Creating significant table up_4edger_mt_ex_vs_wt_ex## Writing excel data for wt_st_vs_wt_ex: 9/16.## After (adj)p filter, the up genes table has 1743 genes.## After (adj)p filter, the down genes table has 1710 genes.## After fold change filter, the up genes table has 326 genes.## After fold change filter, the down genes table has 97 genes.## Writing excel data for mt_st_vs_mt_ex: 10/16.## After (adj)p filter, the up genes table has 963 genes.## After (adj)p filter, the down genes table has 1161 genes.## After fold change filter, the up genes table has 61 genes.## After fold change filter, the down genes table has 6 genes.## Writing excel data for mt_st_vs_wt_st: 11/16.## After (adj)p filter, the up genes table has 1332 genes.## After (adj)p filter, the down genes table has 1391 genes.## After fold change filter, the up genes table has 61 genes.## After fold change filter, the down genes table has 184 genes.## Writing excel data for mt_ex_vs_wt_ex: 12/16.## After (adj)p filter, the up genes table has 522 genes.## After (adj)p filter, the down genes table has 533 genes.## After fold change filter, the up genes table has 40 genes.## After fold change filter, the down genes table has 22 genes.## Printing significant genes to the file: excel/pa_nobatch_significant-v20171019.xlsx## 1/4: Creating significant table up_1deseq_wt_st_vs_wt_ex## 2/4: Creating significant table up_2deseq_mt_st_vs_mt_ex## 3/4: Creating significant table up_3deseq_mt_st_vs_wt_st## 4/4: Creating significant table up_4deseq_mt_ex_vs_wt_ex## Writing excel data for wt_st_vs_wt_ex: 13/16.## After (adj)p filter, the up genes table has 1092 genes.## After (adj)p filter, the down genes table has 1542 genes.## After fold change filter, the up genes table has 292 genes.## After fold change filter, the down genes table has 122 genes.## Writing excel data for mt_st_vs_mt_ex: 14/16.## After (adj)p filter, the up genes table has 311 genes.## After (adj)p filter, the down genes table has 504 genes.## After fold change filter, the up genes table has 50 genes.## After fold change filter, the down genes table has 5 genes.## Writing excel data for mt_st_vs_wt_st: 15/16.## After (adj)p filter, the up genes table has 1330 genes.## After (adj)p filter, the down genes table has 1186 genes.## After fold change filter, the up genes table has 73 genes.## After fold change filter, the down genes table has 190 genes.## Writing excel data for mt_ex_vs_wt_ex: 16/16.## After (adj)p filter, the up genes table has 49 genes.## After (adj)p filter, the down genes table has 56 genes.## After fold change filter, the up genes table has 8 genes.## After fold change filter, the down genes table has 10 genes.## Printing significant genes to the file: excel/pa_nobatch_significant-v20171019.xlsx## 1/4: Creating significant table up_1basic_wt_st_vs_wt_ex## 2/4: Creating significant table up_2basic_mt_st_vs_mt_ex## 3/4: Creating significant table up_3basic_mt_st_vs_wt_st## 4/4: Creating significant table up_4basic_mt_ex_vs_wt_ex## Adding significance bar plots.tt <- pp(file="images/sig_barplot.png", image=pa_sig$sig_bar_plots$deseq)## Writing the image to: images/sig_barplot.png and calling dev.off().index.html sample_estimation.html
pander::pander(sessionInfo())R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
locale: LC_CTYPE=en_US.utf8, LC_NUMERIC=C, LC_TIME=en_US.utf8, LC_COLLATE=en_US.utf8, LC_MONETARY=en_US.utf8, LC_MESSAGES=en_US.utf8, LC_PAPER=en_US.utf8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.utf8 and LC_IDENTIFICATION=C
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: foreach(v.1.4.4), Vennerable(v.3.1.0.9000) and hpgltools(v.2018.03)
loaded via a namespace (and not attached): colorspace(v.1.3-2), rprojroot(v.1.3-2), htmlTable(v.1.12), corpcor(v.1.6.9), XVector(v.0.20.0), GenomicRanges(v.1.32.4), base64enc(v.0.1-3), rstudioapi(v.0.7), roxygen2(v.6.0.1), ggrepel(v.0.8.0), bit64(v.0.9-7), AnnotationDbi(v.1.42.1), xml2(v.1.2.0), codetools(v.0.2-15), splines(v.3.5.1), doParallel(v.1.0.11), robustbase(v.0.93-1.1), geneplotter(v.1.58.0), knitr(v.1.20), Formula(v.1.2-3), Rsamtools(v.1.32.2), annotate(v.1.58.0), cluster(v.2.0.7-1), graph(v.1.58.0), compiler(v.3.5.1), httr(v.1.3.1), backports(v.1.1.2), assertthat(v.0.2.0), Matrix(v.1.2-14), lazyeval(v.0.2.1), limma(v.3.36.2), acepack(v.1.4.1), htmltools(v.0.3.6), prettyunits(v.1.0.2), tools(v.3.5.1), bindrcpp(v.0.2.2), gtable(v.0.2.0), glue(v.1.3.0), GenomeInfoDbData(v.1.1.0), reshape2(v.1.4.3), dplyr(v.0.7.6), Rcpp(v.0.12.17), Biobase(v.2.40.0), Biostrings(v.2.48.0), preprocessCore(v.1.42.0), rtracklayer(v.1.40.3), iterators(v.1.0.10), stringr(v.1.3.1), openxlsx(v.4.1.0), gtools(v.3.8.1), devtools(v.1.13.6), XML(v.3.98-1.12), edgeR(v.3.22.3), DEoptimR(v.1.0-8), directlabels(v.2018.05.22), zlibbioc(v.1.26.0), MASS(v.7.3-50), scales(v.0.5.0), doSNOW(v.1.0.16), hms(v.0.4.2), RBGL(v.1.56.0), parallel(v.3.5.1), SummarizedExperiment(v.1.10.1), RColorBrewer(v.1.1-2), yaml(v.2.1.19), memoise(v.1.1.0), gridExtra(v.2.3), pander(v.0.6.2), ggplot2(v.3.0.0), biomaRt(v.2.36.1), rpart(v.4.1-13), latticeExtra(v.0.6-28), stringi(v.1.2.3), RSQLite(v.2.1.1), genefilter(v.1.62.0), S4Vectors(v.0.18.3), checkmate(v.1.8.5), GenomicFeatures(v.1.32.0), BiocGenerics(v.0.26.0), zip(v.1.0.0), BiocParallel(v.1.14.2), GenomeInfoDb(v.1.16.0), rlang(v.0.2.1), pkgconfig(v.2.0.1), commonmark(v.1.5), matrixStats(v.0.53.1), bitops(v.1.0-6), evaluate(v.0.11), lattice(v.0.20-35), purrr(v.0.2.5), bindr(v.0.1.1), labeling(v.0.3), GenomicAlignments(v.1.16.0), htmlwidgets(v.1.2), bit(v.1.1-14), tidyselect(v.0.2.4), plyr(v.1.8.4), magrittr(v.1.5), DESeq2(v.1.20.0), R6(v.2.2.2), snow(v.0.4-2), IRanges(v.2.14.10), Hmisc(v.4.1-1), DelayedArray(v.0.6.1), DBI(v.1.0.0), pillar(v.1.3.0), foreign(v.0.8-70), withr(v.2.1.2), survival(v.2.42-6), RCurl(v.1.95-4.11), nnet(v.7.3-12), tibble(v.1.4.2), crayon(v.1.3.4), rmarkdown(v.1.10), progress(v.1.2.0), locfit(v.1.5-9.1), grid(v.3.5.1), data.table(v.1.11.4), blob(v.1.1.1), digest(v.0.6.15), xtable(v.1.8-2), stats4(v.3.5.1), munsell(v.0.5.0) and quadprog(v.1.5-5)
message(paste0("This is hpgltools commit: ", get_git_commit()))## If you wish to reproduce this exact build of hpgltools, invoke the following:## > git clone http://github.com/abelew/hpgltools.git## > git reset c730ef178f8e57bbf3819e21cf5e6cfe879e6328## R> packrat::restore()## This is hpgltools commit: Fri Jul 13 17:21:39 2018 -0400: c730ef178f8e57bbf3819e21cf5e6cfe879e6328this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))## Saving to 02_sample_estimation-v20171019.rda.xztmp <- sm(saveme(filename=this_save))