S. cerevisiae differential expression.

1 S. cerevisiae differential expression

2 Filter the raw data

In sample_estimation, I created sc_filt which is precisely what I want.

3 Start with batch in the model

keepers <- list(
    "mcwu_vs_wcwu" = c("mtc_wtu", "wtc_wtu"),
    "wcmu_vs_wcwu" = c("wtc_mtu", "wtc_wtu"))

batch <- all_pairwise(input=sc_filt)

## Finished running DE analyses, collecting outputs.

## Comparing analyses 1/3: wtc_mtu_vs_mtc_wtu

## Comparing analyses 2/3: wtc_wtu_vs_mtc_wtu

## Comparing analyses 3/3: wtc_wtu_vs_wtc_mtu

fsva <- all_pairwise(input=sc_filt, model_batch="sva")

## The be method chose 1 surrogate variable(s).

## Estimate type 'sva' is shorthand for 'sva_unsupervised'.

## Other sva options include: sva_supervised and svaseq.

## Attempting sva unsupervised surrogate estimation.

## Finished running DE analyses, collecting outputs.

## Comparing analyses 1/3: wtc_mtu_vs_mtc_wtu

## Comparing analyses 2/3: wtc_wtu_vs_mtc_wtu

## Comparing analyses 3/3: wtc_wtu_vs_wtc_mtu

batch_write <- combine_de_tables(all_pairwise_result=batch,
                                 keepers=keepers,
                                 excel=paste0("excel/batch_in_model_differential-v", ver, ".xlsx"))

## Deleting the file excel/batch_in_model_differential-v20170515.xlsx before writing the tables.

## Writing a legend of columns.

## Working on 1/2: mcwu_vs_wcwu

## Found inverse table with wtc_wtu_vs_mtc_wtu

## Working on 2/2: wcmu_vs_wcwu

## Found inverse table with wtc_wtu_vs_wtc_mtu

## Adding venn plots for mcwu_vs_wcwu.

## Loading required package: methods

## Limma expression coefficients for mcwu_vs_wcwu; R^2: 0.928; equation: y = 0.968x + 0.298

## Edger expression coefficients for mcwu_vs_wcwu; R^2: 0.927; equation: y = 0.971x + 0.308

## DESeq2 expression coefficients for mcwu_vs_wcwu; R^2: 0.926; equation: y = 0.988x + 0.0829

## Adding venn plots for wcmu_vs_wcwu.

## Limma expression coefficients for wcmu_vs_wcwu; R^2: 0.95; equation: y = 0.899x + 0.839

## Edger expression coefficients for wcmu_vs_wcwu; R^2: 0.95; equation: y = 0.899x + 1.06

## DESeq2 expression coefficients for wcmu_vs_wcwu; R^2: 0.955; equation: y = 0.882x + 0.999

## Writing summary information.

## The sheet: pairwise_summary is in legend, mcwu_vs_wcwu, wcmu_vs_wcwu, pairwise_summary.

## Attempting to add the comparison plot to pairwise_summary at row: 20 and column: 1

## Performing save of the workbook.

batch_sig <- sm(extract_significant_genes(combined=batch_write,
                                          excel=paste0("excel/batch_in_model_significant-v", ver, ".xlsx")))

fsva_write <- combine_de_tables(all_pairwise_result=fsva,
                                keepers=keepers,
                                excel=paste0("excel/sva_in_model_differential-v", ver, ".xlsx"))

## Deleting the file excel/sva_in_model_differential-v20170515.xlsx before writing the tables.

## Writing a legend of columns.

## Working on 1/2: mcwu_vs_wcwu

## Found inverse table with wtc_wtu_vs_mtc_wtu

## Working on 2/2: wcmu_vs_wcwu

## Found inverse table with wtc_wtu_vs_wtc_mtu

## Adding venn plots for mcwu_vs_wcwu.

## Limma expression coefficients for mcwu_vs_wcwu; R^2: 0.936; equation: y = 0.996x + 0.102

## Edger expression coefficients for mcwu_vs_wcwu; R^2: 0.937; equation: y = 0.996x + 0.0632

## DESeq2 expression coefficients for mcwu_vs_wcwu; R^2: 0.935; equation: y = 1.01x - 0.247

## Adding venn plots for wcmu_vs_wcwu.

## Limma expression coefficients for wcmu_vs_wcwu; R^2: 0.667; equation: y = 0.684x + 2.46

## Edger expression coefficients for wcmu_vs_wcwu; R^2: 0.673; equation: y = 0.688x + 3.12

## DESeq2 expression coefficients for wcmu_vs_wcwu; R^2: 0.731; equation: y = 0.725x + 2.48

## Writing summary information.

## The sheet: pairwise_summary is in legend, mcwu_vs_wcwu, wcmu_vs_wcwu, pairwise_summary.

## Attempting to add the comparison plot to pairwise_summary at row: 20 and column: 1

## Performing save of the workbook.

fsva_sig <- sm(extract_significant_genes(combined=fsva_write,
                                         excel=paste0("excel/sva_in_model_significant-v", ver, ".xlsx")))

sc_nobatch <- all_pairwise(input=sc_norm, model_batch=FALSE)

## Finished running DE analyses, collecting outputs.

## Comparing analyses 1/1: wtc_wtu_vs_mtc_wtu

nobatch_write <- combine_de_tables(all_pairwise_result=sc_nobatch,
                                   excel=paste0("excel/nobatch_differential-v", ver, ".xlsx"))

## Deleting the file excel/nobatch_differential-v20170515.xlsx before writing the tables.

## Writing a legend of columns.

## Working on table 1/1: wtc_wtu_vs_mtc_wtu

## This can do comparisons among the following columns in the pairwise result:

## mtc_wtu, wtc_wtu, wtc_wtu_vs_mtc_wtu

## Actually comparing wtc_wtu and mtc_wtu.

## This can do comparisons among the following columns in the pairwise result:

## mtc_wtu, wtc_wtu

## Actually comparing wtc_wtu and mtc_wtu.

## This can do comparisons among the following columns in the pairwise result:

## mtc_wtu, wtc_wtu

## Actually comparing wtc_wtu and mtc_wtu.

## Adding venn plots for wtc_wtu_vs_mtc_wtu.

## Limma expression coefficients for wtc_wtu_vs_mtc_wtu; R^2: 0.915; equation: y = 0.957x + 0.084

## Edger expression coefficients for wtc_wtu_vs_mtc_wtu; R^2: 0.969; equation: y = 0.987x + 0.165

## DESeq2 expression coefficients for wtc_wtu_vs_mtc_wtu; R^2: 0.973; equation: y = 0.988x + 0.0618

## Writing summary information.

## The sheet: pairwise_summary is in legend, wtc_wtu_vs_mtc_wtu, pairwise_summary.

## Attempting to add the comparison plot to pairwise_summary at row: 19 and column: 1

## Performing save of the workbook.

3.1 Compare sva to batch in model

batch_fsva <- batch_write$data$mcwu_vs_wcwu[, c("limma_logfc", "limma_adjp")]
batch_fsva <- merge(batch_fsva, fsva_write$data$mcwu_vs_wcwu[, c("limma_logfc", "limma_adjp")], by="row.names")
batch_fsva[, 1] <- as.numeric(batch_fsva[, 2])
batch_fsva[, 3] <- as.numeric(batch_fsva[, 4])
test_batch_fsva <- plot_linear_scatter(df=batch_fsva[, c(2, 4)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

test_batch_fsva$scatter

4 Perform a cross reference

Carol has a sheet of data from the Jacobson lab against which she wishes to cross reference this data.

jac <- read.csv("external_data/jacobson.csv")
batch_table <- batch_write$data$wcmu_vs_wcwu
merged <- merge(jac, batch_table, by.x="Gene", by.y="geneid")
dim(merged)

## [1] 1549   46

merged <- merged[, c("B.upf1", "limma_logfc", "deseq_logfc", "edger_logfc")]
merged[[1]] <- as.numeric(merged[[1]])
merged[[1]] <- log2(x=merged[[1]] + 0.5)
merged[[2]] <- as.numeric(merged[[2]])
merged[[3]] <- as.numeric(merged[[3]])
merged[[4]] <- as.numeric(merged[[4]])

limma_plot <- plot_linear_scatter(df=merged[, c(1,2)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

pdf(file="images/linear_scatter_limma_jacobson.pdf")
limma_plot$scatter
dev.off()

## png 
##   2

limma_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 32, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6056 0.6650
## sample estimates:
##    cor 
## 0.6362

deseq_plot <- plot_linear_scatter(df=merged[, c(1,3)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

deseq_plot$scatter

deseq_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 30, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5686 0.6323
## sample estimates:
##    cor 
## 0.6014

edger_plot <- plot_linear_scatter(df=merged[, c(1,4)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

edger_plot$scatter

edger_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 31, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5931 0.6539
## sample estimates:
##    cor 
## 0.6244

jac <- read.csv("external_data/jacobson.csv")
sva_table <- fsva_write$data$wcmu_vs_wcwu
merged <- merge(jac, sva_table, by.x="Gene", by.y="geneid")
dim(merged)

## [1] 1549   46

merged <- merged[, c("B.upf1", "limma_logfc", "deseq_logfc", "edger_logfc")]
merged[[1]] <- as.numeric(merged[[1]])
merged[[1]] <- log2(x=merged[[1]] + 0.5)
merged[[2]] <- as.numeric(merged[[2]])
merged[[3]] <- as.numeric(merged[[3]])
merged[[4]] <- as.numeric(merged[[4]])
limma_plot <- plot_linear_scatter(df=merged[, c(1,2)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

limma_plot$scatter

limma_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 7.3, df = 1500, p-value = 4e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1346 0.2308
## sample estimates:
##    cor 
## 0.1831

cor.test(merged[, 1], merged[, 2], method="spearman")

## Warning in cor.test.default(merged[, 1], merged[, 2], method = "spearman"):
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  merged[, 1] and merged[, 2]
## S = 4.3e+08, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.3013

deseq_plot <- plot_linear_scatter(df=merged[, c(1,3)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

deseq_plot$scatter

deseq_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 11, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2123 0.3052
## sample estimates:
##    cor 
## 0.2593

edger_plot <- plot_linear_scatter(df=merged[, c(1,4)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

edger_plot$scatter

edger_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 6.1, df = 1500, p-value = 1e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1046 0.2018
## sample estimates:
##    cor 
## 0.1536

jac <- read.csv("external_data/jacobson.csv")
nobatch_table <- nobatch_write$data$wtc_wtu_vs_mtc_wtu
merged <- merge(jac, nobatch_table, by.x="Gene", by.y="geneid")
dim(merged)

## [1] 1529   46

merged <- merged[, c("B.upf1", "limma_logfc", "deseq_logfc", "edger_logfc")]
merged[[1]] <- as.numeric(merged[[1]])
merged[[1]] <- log2(x=merged[[1]] + 0.5)
merged[[2]] <- as.numeric(merged[[2]])
merged[[3]] <- as.numeric(merged[[3]])
merged[[4]] <- as.numeric(merged[[4]])
limma_plot <- plot_linear_scatter(df=merged[, c(1,2)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

limma_plot$scatter

limma_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 7.6, df = 1500, p-value = 5e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1425 0.2391
## sample estimates:
##    cor 
## 0.1913

cor.test(merged[, 1], merged[, 2], method="spearman")

## Warning in cor.test.default(merged[, 1], merged[, 2], method = "spearman"):
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  merged[, 1] and merged[, 2]
## S = 4.8e+08, p-value = 6e-14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.1904

deseq_plot <- plot_linear_scatter(df=merged[, c(1,3)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

deseq_plot$scatter

deseq_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 8.5, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1635 0.2593
## sample estimates:
##    cor 
## 0.2119

edger_plot <- plot_linear_scatter(df=merged[, c(1,4)], loess=TRUE)

## Used Bon Ferroni corrected t test(s) between columns.

edger_plot$scatter

edger_plot$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 8.5, df = 1500, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1634 0.2591
## sample estimates:
##    cor 
## 0.2118

library("pander")
pander(sessionInfo())

R version 3.3.3 (2017-03-06)

**Platform:** x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: methods, stats, graphics, grDevices, utils, datasets and base

other attached packages: pander(v.0.6.0), Vennerable(v.3.1.0.9000), ruv(v.0.9.6), hpgltools(v.2017.01) and rmarkdown(v.1.5)

loaded via a namespace (and not attached): nlme(v.3.1-131), bitops(v.1.0-6), devtools(v.1.13.1), doParallel(v.1.0.10), RColorBrewer(v.1.1-2), rprojroot(v.1.2), GenomeInfoDb(v.1.10.3), tools(v.3.3.3), backports(v.1.0.5), R6(v.2.2.1), rpart(v.4.1-11), Hmisc(v.4.0-3), DBI(v.0.6-1), lazyeval(v.0.2.0), BiocGenerics(v.0.20.0), mgcv(v.1.8-17), colorspace(v.1.3-2), nnet(v.7.3-12), withr(v.1.0.2), gridExtra(v.2.2.1), DESeq2(v.1.14.1), compiler(v.3.3.3), preprocessCore(v.1.36.0), graph(v.1.52.0), Biobase(v.2.34.0), htmlTable(v.1.9), xml2(v.1.1.1), rtracklayer(v.1.34.2), labeling(v.0.3), scales(v.0.4.1), checkmate(v.1.8.2), DEoptimR(v.1.0-8), robustbase(v.0.92-7), genefilter(v.1.56.0), RBGL(v.1.50.0), commonmark(v.1.2), stringr(v.1.2.0), digest(v.0.6.12), Rsamtools(v.1.26.2), foreign(v.0.8-68), XVector(v.0.14.1), base64enc(v.0.1-3), htmltools(v.0.3.6), limma(v.3.30.13), htmlwidgets(v.0.8), rlang(v.0.1.1), RSQLite(v.1.1-2), BiocParallel(v.1.8.2), gtools(v.3.5.0), acepack(v.1.4.1), RCurl(v.1.95-4.8), magrittr(v.1.5), Formula(v.1.2-1), Matrix(v.1.2-10), Rcpp(v.0.12.11), munsell(v.0.4.3), S4Vectors(v.0.12.2), stringi(v.1.1.5), yaml(v.2.1.14), edgeR(v.3.16.5), SummarizedExperiment(v.1.4.0), zlibbioc(v.1.20.0), plyr(v.1.8.4), grid(v.3.3.3), parallel(v.3.3.3), ggrepel(v.0.6.5), crayon(v.1.3.2), lattice(v.0.20-35), Biostrings(v.2.42.1), splines(v.3.3.3), GenomicFeatures(v.1.26.4), annotate(v.1.52.1), locfit(v.1.5-9.1), knitr(v.1.16), GenomicRanges(v.1.26.4), geneplotter(v.1.52.0), reshape2(v.1.4.2), codetools(v.0.2-15), biomaRt(v.2.30.0), stats4(v.3.3.3), XML(v.3.98-1.7), evaluate(v.0.10), latticeExtra(v.0.6-28), data.table(v.1.10.4), foreach(v.1.4.3), testthat(v.1.0.2), gtable(v.0.2.0), ggplot2(v.2.2.1), openxlsx(v.4.0.17), xtable(v.1.8-2), roxygen2(v.6.0.1), survival(v.2.41-3), tibble(v.1.3.1), iterators(v.1.0.8), GenomicAlignments(v.1.10.1), AnnotationDbi(v.1.36.2), memoise(v.1.1.0), IRanges(v.2.8.2), cluster(v.2.0.6) and sva(v.3.22.0)