20190116: Differential Expression for TNSeq of Spyogenes. (2017)
atb
2019-01-15
1 Differential Expression, Spyogenes: 20190115
This worksheet performs a set of pairwise comparisons of the counts for this tnseq experiment given the sva surrogate adjustments.
This is being performed via limma/edger/deseq2, which one may rightly argue is not entirely appropriate for TNSeq data. On the other hand, as a general metric of changing fitness, it provides some interesting information.
1.1 Perform Pairwise analyses
all_pairwise() is responsible for following:
- Creating an experimental model which includes surrogate estimates from svaseq.
- Calling limma, EdgeR, and DESEq2 with the data and that model for every possible pairwise contrast in the data. (Optionally it will now do EBSeq, too.)
- Calling a basic method which simply subtracts log2 scaled values as a form of negative control.
- Collecting those results.
Among the summaries it provides is a heatmap of how similar each method was for every contrast.
## All contrasts are quite similar.
## The lowest correlation coefficient among the methods/contrasts is ~ 0.6 and is
## between the 'real' methods and my stupid basic method, which mostly tells me
## that sva actually did something.
initial_diff$comparison[[7]]
## Here we see that the log2FCs from limma and DESeq2 are very similar.
## The blue line is the identity; so DESeq is seeing slightly higher
## values than limma almost across the board.
initial_diff$comparison[[9]]
1.2 Combine the tables
combine_de_tables() does what it says on the tin. When provided an excel filename, it will also put the combined tables into the excel file and try to add in some pretty plots to the resulting excel file.
keepers <- list("t2t1" = c("rpmi_t2", "thy_t1"),
"t3t1" = c("rpmi_t3", "thy_t1"),
"t2zn" = c("rpmi_t2_highzn", "rpmi_t2_lowzn"),
"t3zn" = c("rpmi_t3_highzn", "rpmi_t3_lowzn"),
"time_highzn" = c("rpmi_t3_highzn", "rpmi_t2_highzn"),
"time_lowzn" = c("rpmi_t3_lowzn", "rpmi_t2_lowzn"),
"t2_rpmi_highzn" = c("rpmi_t2_highzn", "rpmi_t2"),
"t3_rpmi_highzn" = c("rpmi_t3_highzn", "rpmi_t3"),
"t2_rpmi_lowzn" = c("rpmi_t2_lowzn", "rpmi_t2"),
"t3_rpmi_lowzn" = c("rpmi_t3_lowzn", "rpmi_t3")
)
## Each element of keepers it a name of the contrast to perform followed by a 2 element list
## with the numerator and denominator.
initial_write <- sm(combine_de_tables(initial_diff,
excel=paste0("excel/", rundate, "_initial_diff-v", ver, ".xlsx"),
keepers=keepers))
The data generated by combine_de_tables() is supposed to give me some plots which we may use to evaluate how well the data came together as well as how similar the final results are for limma/deseq/edger.
1.3 Find significant genes
sig_write <- sm(extract_significant_genes(
initial_write,
excel=paste0("excel/", rundate, "_significant-v", ver, ".xlsx")))## If you wish to reproduce this exact build of hpgltools, invoke the following:## > git clone http://github.com/abelew/hpgltools.git## > git reset c89fbb1a6badaccdac9bfbc109bd1fe2d673b639## This is hpgltools commit: Sun Jan 13 20:46:37 2019 -0500: c89fbb1a6badaccdac9bfbc109bd1fe2d673b639R version 3.5.2 (2018-12-20)
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: parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: bindrcpp(v.0.2.2), ruv(v.0.9.7), hpgltools(v.2018.11), Biobase(v.2.42.0) and BiocGenerics(v.0.28.0)
loaded via a namespace (and not attached): tidyselect(v.0.2.5), lme4(v.1.1-19), RSQLite(v.2.1.1), AnnotationDbi(v.1.44.0), htmlwidgets(v.1.3), grid(v.3.5.2), BiocParallel(v.1.16.5), devtools(v.2.0.1), munsell(v.0.5.0), codetools(v.0.2-16), preprocessCore(v.1.45.0), units(v.0.6-2), withr(v.2.1.2), colorspace(v.1.3-2), GOSemSim(v.2.8.0), knitr(v.1.21), rstudioapi(v.0.9.0), stats4(v.3.5.2), Vennerable(v.3.1.0.9000), robustbase(v.0.93-3), DOSE(v.3.8.0), labeling(v.0.3), urltools(v.1.7.1), GenomeInfoDbData(v.1.2.0), bit64(v.0.9-7), farver(v.1.1.0), rprojroot(v.1.3-2), xfun(v.0.4), R6(v.2.3.0), doParallel(v.1.0.14), GenomeInfoDb(v.1.18.1), locfit(v.1.5-9.1), bitops(v.1.0-6), fgsea(v.1.8.0), gridGraphics(v.0.3-0), DelayedArray(v.0.8.0), assertthat(v.0.2.0), scales(v.1.0.0), ggraph(v.1.0.2), nnet(v.7.3-12), enrichplot(v.1.2.0), gtable(v.0.2.0), sva(v.3.30.1), processx(v.3.2.1), rlang(v.0.3.1), genefilter(v.1.64.0), splines(v.3.5.2), rtracklayer(v.1.42.1), lazyeval(v.0.2.1), acepack(v.1.4.1), europepmc(v.0.3), checkmate(v.1.9.0), yaml(v.2.2.0), reshape2(v.1.4.3), GenomicFeatures(v.1.34.1), backports(v.1.1.3), qvalue(v.2.14.1), Hmisc(v.4.1-1), RBGL(v.1.58.1), clusterProfiler(v.3.10.1), tools(v.3.5.2), usethis(v.1.4.0), ggplotify(v.0.0.3), ggplot2(v.3.1.0), gplots(v.3.0.1), RColorBrewer(v.1.1-2), sessioninfo(v.1.1.1), ggridges(v.0.5.1), Rcpp(v.1.0.0), plyr(v.1.8.4), base64enc(v.0.1-3), progress(v.1.2.0), zlibbioc(v.1.28.0), purrr(v.0.2.5), RCurl(v.1.95-4.11), ps(v.1.3.0), prettyunits(v.1.0.2), rpart(v.4.1-13), viridis(v.0.5.1), cowplot(v.0.9.4), S4Vectors(v.0.20.1), SummarizedExperiment(v.1.12.0), ggrepel(v.0.8.0), cluster(v.2.0.7-1), colorRamps(v.2.3), fs(v.1.2.6), variancePartition(v.1.12.1), magrittr(v.1.5), data.table(v.1.11.8), openxlsx(v.4.1.0), DO.db(v.2.9), triebeard(v.0.3.0), packrat(v.0.5.0), matrixStats(v.0.54.0), pkgload(v.1.0.2), hms(v.0.4.2), evaluate(v.0.12), xtable(v.1.8-3), pbkrtest(v.0.4-7), XML(v.3.98-1.16), IRanges(v.2.16.0), gridExtra(v.2.3), testthat(v.2.0.1), compiler(v.3.5.2), biomaRt(v.2.38.0), tibble(v.2.0.0), KernSmooth(v.2.23-15), crayon(v.1.3.4), minqa(v.1.2.4), htmltools(v.0.3.6), mgcv(v.1.8-26), corpcor(v.1.6.9), Formula(v.1.2-3), tidyr(v.0.8.2), geneplotter(v.1.60.0), DBI(v.1.0.0), tweenr(v.1.0.1), MASS(v.7.3-51.1), Matrix(v.1.2-15), cli(v.1.0.1), quadprog(v.1.5-5), gdata(v.2.18.0), bindr(v.0.1.1), igraph(v.1.2.2), GenomicRanges(v.1.34.0), pkgconfig(v.2.0.2), rvcheck(v.0.1.3), GenomicAlignments(v.1.18.1), foreign(v.0.8-71), xml2(v.1.2.0), foreach(v.1.4.4), annotate(v.1.60.0), XVector(v.0.22.0), stringr(v.1.3.1), callr(v.3.1.1), digest(v.0.6.18), graph(v.1.60.0), Biostrings(v.2.50.2), rmarkdown(v.1.11), fastmatch(v.1.1-0), htmlTable(v.1.13.1), edgeR(v.3.24.3), directlabels(v.2018.05.22), Rsamtools(v.1.34.0), gtools(v.3.8.1), nloptr(v.1.2.1), nlme(v.3.1-137), jsonlite(v.1.6), desc(v.1.2.0), viridisLite(v.0.3.0), limma(v.3.38.3), pillar(v.1.3.1), lattice(v.0.20-38), DEoptimR(v.1.0-8), httr(v.1.4.0), pkgbuild(v.1.0.2), survival(v.2.43-3), GO.db(v.3.7.0), glue(v.1.3.0), remotes(v.2.0.2), zip(v.1.0.0), UpSetR(v.1.3.3), iterators(v.1.0.10), pander(v.0.6.3), bit(v.1.1-14), ggforce(v.0.1.3), stringi(v.1.2.4), blob(v.1.1.1), DESeq2(v.1.22.2), latticeExtra(v.0.6-28), caTools(v.1.17.1.1), memoise(v.1.1.0) and dplyr(v.0.7.8)
this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))## Saving to 20190115_03_differential_expression-v20190115.rda.xz