2.1 No batch in the model

hs_macr_nobatch <- sm(all_pairwise(input=hs_macr_lowfilt, model_batch=FALSE, limma_method="robust"))

## wow, all tools including basic agree almost completely
medians_by_condition <- hs_macr_nobatch$basic$medians
hs_macr_nobatch_tables <- sm(combine_de_tables(hs_macr_nobatch,
                                               excel=paste0("excel/hs_macr_nobatch-v", ver, ".xlsx"),
                                               keepers=hs_contrasts,
                                               extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

hs_macr_nobatch_sig <- sm(extract_significant_genes(hs_macr_nobatch_tables,
                                                    excel=paste0("excel/hs_macr_nobatch_significant-v", ver, ".xlsx"),
                                                    according_to="all"))

## Error in extract_significant_genes(hs_macr_nobatch_tables, excel = paste0("excel/hs_macr_nobatch_significant-v", : object 'hs_macr_nobatch_tables' not found

2.2 Batch in the model

In this attempt, we add a batch factor in the experimental model and see how it does.

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
hs_macr_batch <- sm(all_pairwise(input=hs_macr_lowfilt, limma_method="robust", parallel=FALSE))

medians_by_condition <- hs_macr_batch$basic$medians
hs_macr_batch_tables <- sm(combine_de_tables(
  hs_macr_batch,
  excel=paste0("excel/hs_macr_batchmodel-v", ver, ".xlsx"),
  sig_excel=paste0("excel/hs_macro_batchmodel_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

4.1 Batch estimated with SVA

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
hs_macr_sva <- sm(all_pairwise(
  input=hs_macr_lowfilt,
  model_batch="svaseq",
  limma_method="robust"))

medians_by_condition <- hs_macr_sva$basic$medians
hs_macr_sva_tables <- sm(combine_de_tables(
  hs_macr_sva,
  excel=paste0("excel/hs_macr_sva-v", ver, ".xlsx"),
  sig_excel=paste0("excel/hs_macr_sva_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

hs_macr_sva_ma_limma <- extract_de_plots(
  pairwise=hs_macr_sva,
  type="limma",
  table="sh_vs_chr")

hs_macr_sva_ma_limma$ma$plot

4.2 Batch correction via ruv residuals

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
## Bizarrely, sometimes if one runs this, it gives an error "Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'RUVr' for signature '"matrix", "logical", "numeric", "NULL"'"  -- however, if one then simply runs it again it works fine.
## I am going to assume that it is because I do not explicitly invoke the library.
## library(ruv)  ## hopefully a small code change made this not needed.
testme <- get_model_adjust(input=hs_macr_lowfilt, estimate_type="ruv_residuals")

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

## Attempting ruvseq residual surrogate estimation with 2 surrogates.

## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'RUVr' for signature '"matrix", "logical", "numeric", "NULL"'

## hmm I got the RUVr error again, but when I ran it manually did not.
## Even more strangely, if I just run the same thing again, no error...
testme <- get_model_adjust(input=hs_macr_lowfilt, estimate_type="ruv_residuals")

## The be method chose 2 surrogate variable(s).
## Attempting ruvseq residual surrogate estimation with 2 surrogates.

hs_macr_ruvres <- sm(all_pairwise(
  input=hs_macr_lowfilt,
  model_batch="ruv_residuals",
  limma_method="robust"))

medians_by_condition <- hs_macr_ruvres$basic$medians
hs_macr_ruvres_tables <- sm(combine_de_tables(
  hs_macr_ruvres,
  excel=paste0("excel/hs_macr_ruvres-v", ver, ".xlsx"),
  sig_excel=paste0("excel/hs_macr_ruvres_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

4.3 Batch correction with pca

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
hs_macr_pca <- sm(all_pairwise(
  input=hs_macr_lowfilt,
  model_batch="pca",
  limma_method="robust"))

medians_by_condition <- hs_macr_pca$basic$medians
hs_macr_pca_tables <- sm(combine_de_tables(
  hs_macr_pca,
  excel=paste0("excel/hs_macr_pca-v", ver, ".xlsx"),
  sig_excel=paste0("excel/hs_macr_pca_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

4.4 Batch correction with ruv empirical

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
hs_macr_ruvemp <- sm(all_pairwise(
  input=hs_macr_lowfilt,
  model_batch="ruv_empirical",
  limma_method="robust"))

medians_by_condition <- hs_macr_ruvemp$basic$medians
hs_macr_ruvemp_tables <- sm(combine_de_tables(
  hs_macr_ruvemp,
  excel=paste0("excel/hs_macr_ruvemp-v", ver, ".xlsx"),
  sig_excel=paste0("excel/hs_macr_ruvemp_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

4.5 Batch correction with combat

Then repeat with the batch-corrected data and see the differences.

hs_macr_combat <- sm(all_pairwise(
  input=hs_macr_combat_norm,
  force=TRUE,
  limma_method="robust"))

## Error in cpm.default(count_table, ...): library sizes should be finite and non-negative

medians_by_condition <- hs_macr_combat$basic$medians

## Error in eval(expr, envir, enclos): object 'hs_macr_combat' not found

hs_macr_combat_tables <- sm(combine_de_tables(
  hs_macr_combat,
  excel=paste0("excel/hs_macr_combat-v", ver, ".xlsx"),
  sig_excep=paste0("excel/hs_macr_combat_significant-v", ver, ".xlsx"),
  keepers=hs_contrasts,
  extra_annot=medians_by_condition))

## Error in combine_de_tables(hs_macr_combat, excel = paste0("excel/hs_macr_combat-v", : object 'hs_macr_combat' not found

hs_macr_combat_ma_limma <- extract_de_plots(
  pairwise=hs_macr_combat,
  type="limma",
  table="sh_vs_chr")

## Error in extract_de_plots(pairwise = hs_macr_combat, type = "limma", table = "sh_vs_chr"): object 'hs_macr_combat' not found

hs_macr_combat_ma_limma$ma$plot

## Error in eval(expr, envir, enclos): object 'hs_macr_combat_ma_limma' not found

hs_macr_combat_ma_edger <- extract_de_plots(
  pairwise=hs_macr_combat,
  type="edger",
  table="sh_vs_chr")

## Error in extract_de_plots(pairwise = hs_macr_combat, type = "edger", table = "sh_vs_chr"): object 'hs_macr_combat' not found

hs_macr_combat_ma_edger$ma$plot

## Error in eval(expr, envir, enclos): object 'hs_macr_combat_ma_edger' not found

hs_macr_combat_ma_deseq <- extract_de_plots(
  pairwise=hs_macr_combat,
  type="deseq",
  table="sh_vs_chr")

## Error in extract_de_plots(pairwise = hs_macr_combat, type = "deseq", table = "sh_vs_chr"): object 'hs_macr_combat' not found

hs_macr_combat_ma_deseq$ma$plot

## Error in eval(expr, envir, enclos): object 'hs_macr_combat_ma_deseq' not found

6.1 Compare DeSeq / Basic without batch in model

hs_macr_nobatch_basic <- merge(
  hs_macr_nobatch$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_nobatch_basic) <- hs_macr_nobatch_basic[["Row.names"]]
hs_macr_nobatch_logfc <- hs_macr_nobatch_basic[, c("logFC.x", "logFC.y")]
colnames(hs_macr_nobatch_logfc) <- c("nobatch", "basic")
lfc_nb_b <- sm(plot_linear_scatter(hs_macr_nobatch_logfc, pretty_colors=FALSE))
lfc_nb_b$scatter

lfc_nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 120, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8112 0.8266
## sample estimates:
##    cor 
## 0.8191

hs_macr_nobatch_p <- hs_macr_nobatch_basic[, c("P.Value","p")]
hs_macr_nobatch_p[[2]] <- as.numeric(hs_macr_nobatch_p[[2]])
colnames(hs_macr_nobatch_p) <- c("nobatch","basic")
hs_macr_nobatch_p <- -1 * log(hs_macr_nobatch_p)
hs_macr_p_nb_b <- sm(plot_linear_scatter(hs_macr_nobatch_p, pretty_colors=FALSE))
hs_macr_p_nb_b$scatter

hs_macr_p_nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 72, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6379 0.6649
## sample estimates:
##    cor 
## 0.6516

6.2 Compare SVA to batch in model, DESeq

hs_macr_sva_batch <- merge(
  hs_macr_sva$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_sva_batch) <- hs_macr_sva_batch[["Row.names"]]
hs_macr_sva_logfc <- hs_macr_sva_batch[, c("logFC.x","logFC.y")]
colnames(hs_macr_sva_logfc) <- c("sva","batch")
hs_macr_lfc_b_s <- sm(plot_linear_scatter(hs_macr_sva_logfc, pretty_colors=FALSE))
hs_macr_lfc_b_s$scatter

hs_macr_lfc_b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 1000, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9967 0.9970
## sample estimates:
##    cor 
## 0.9968

lfc_b_s <- hs_macr_sva_batch[, c("logFC.x", "logFC.y", "P.Value.x", "P.Value.y")]
colnames(lfc_b_s) <- c("l2fcsva", "l2fcbatch", "psva", "pbatch")
hs_macr_lfc_b_s$scatter

cutoff <- 0.1
lfc_b_s$state <- ifelse(lfc_b_s$psva > cutoff & lfc_b_s$pbatch > cutoff, "bothinsig",
                 ifelse(lfc_b_s$psva <= cutoff & lfc_b_s$pbatch <= cutoff, "bothsig",
                 ifelse(lfc_b_s$psva <= cutoff, "svasig", "batchsig")))
##lfcp_b_s$lfcstate <- ifelse(lfcp_b_s$l2fcsva >= 0.75 & lfcp_b_s$l2fcbatch, "", "")
num_bothinsig <- sum(lfc_b_s$state == "bothinsig")
num_bothsig <- sum(lfc_b_s$state == "bothsig")
num_svasig <- sum(lfc_b_s$state == "svasig")
num_batchsig <- sum(lfc_b_s$state == "batchsig")

library(ggplot2)

## Need help? Try Stackoverflow:
## https://stackoverflow.com/tags/ggplot2.

aes_color = "(l2fcsva >= 0.75 | l2fcsva <= -0.75 | l2fcbatch >= 0.75 | l2fcbatch <= -0.75)"

plt <- ggplot2::ggplot(lfc_b_s, aes_string(x="l2fcsva", y="l2fcbatch")) +
    ## ggplot2::geom_point(stat="identity", size=2, alpha=0.2, aes_string(shape="as.factor(aes_color)", colour="as.factor(state)", fill="as.factor(state)")) +
    ggplot2::geom_abline(colour="blue", slope=1, intercept=0, size=0.5) +
    ggplot2::geom_hline(yintercept=c(-0.75, 0.75), color="red", size=0.5) +
    ggplot2::geom_vline(xintercept=c(-0.75, 0.75), color="red", size=0.5) +
    ggplot2::geom_point(stat="identity", size=2, alpha=0.2, aes_string(colour="as.factor(state)", fill="as.factor(state)")) +
    ggplot2::scale_color_manual(name="state", values=c("bothinsig"="grey", "bothsig"="forestgreen", "svasig"="darkred", "batchsig"="darkblue")) +
    ggplot2::scale_fill_manual(name="state", values=c("bothinsig"="grey", "bothsig"="forestgreen", "svasig"="darkred", "batchsig"="darkblue"),
                               labels=c(
                                   paste0("Both InSig.: ", num_bothinsig),
                                   paste0("Both Sig.: ", num_bothsig),
                                   paste0("Sva Sig.: ", num_svasig),
                                   paste0("Batch Sig.: ", num_batchsig)),
                               guide=ggplot2::guide_legend(override.aes=aes(size=3, fill="grey"))) +
    ggplot2::guides(fill=ggplot2::guide_legend(override.aes=list(size=3))) +
    ggplot2::theme_bw()
plt

6.3 Compare ruvresid to batch in model, DESeq

hs_macr_batch_ruvresid_deseq <- merge(
  hs_macr_ruvres$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_ruvresid_deseq) <- hs_macr_batch_ruvresid_deseq[["Row.names"]]
hs_macr_batch_ruvresid_logfc <- hs_macr_batch_ruvresid_deseq[, c("logFC.x","logFC.y")]
colnames(hs_macr_batch_ruvresid_logfc) <- c("nobatch","basic")
lfc_ruv_bat <- plot_linear_scatter(hs_macr_batch_ruvresid_logfc, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

lfc_ruv_bat$scatter

lfc_ruv_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 140, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8574 0.8693
## sample estimates:
##    cor 
## 0.8634

6.4 Compare pca to batch in model, DESeq

hs_macr_batch_pca_deseq <- merge(
  hs_macr_pca$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_pca_deseq) <- hs_macr_batch_pca_deseq[["Row.names"]]
hs_macr_batch_pca_logfc <- hs_macr_batch_pca_deseq[, c("logFC.x", "logFC.y")]
colnames(hs_macr_batch_pca_logfc) <- c("nobatch","basic")
lfc_pca_bat <- plot_linear_scatter(hs_macr_batch_pca_logfc, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

lfc_pca_bat$scatter

lfc_pca_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 50, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4989 0.5332
## sample estimates:
##    cor 
## 0.5162

6.5 Compare ruv empirical to batch in model, DESeq

hs_macr_batch_ruvemp_deseq <- merge(
  hs_macr_ruvemp$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_ruvemp_deseq) <- hs_macr_batch_ruvemp_deseq[["Row.names"]]
hs_macr_batch_ruvemp_logfc <- hs_macr_batch_ruvemp_deseq[, c("logFC.x","logFC.y")]
colnames(hs_macr_batch_ruvemp_logfc) <- c("nobatch","basic")
lfc_ruvemp_bat <- sm(plot_linear_scatter(hs_macr_batch_ruvemp_logfc, pretty_colors=FALSE))
lfc_ruvemp_bat$scatter

lfc_ruvemp_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 57, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5487 0.5806
## sample estimates:
##    cor 
## 0.5648

6.6 Compare combat to batch in model, DESeq

hs_macr_combat_batch <- merge(
  hs_macr_combat$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")

## Error in merge(hs_macr_combat$deseq$all_tables$sh_vs_chr, hs_macr_batch$deseq$all_tables$sh_vs_chr, : object 'hs_macr_combat' not found

rownames(hs_macr_combat_batch) <- hs_macr_combat_batch[["Row.names"]]

## Error in eval(expr, envir, enclos): object 'hs_macr_combat_batch' not found

hs_macr_combat_batch <- hs_macr_combat_batch[, c("logFC.x","logFC.y")]

## Error in eval(expr, envir, enclos): object 'hs_macr_combat_batch' not found

colnames(hs_macr_combat_batch) <- c("batch","combat")

## Error in colnames(hs_macr_combat_batch) <- c("batch", "combat"): object 'hs_macr_combat_batch' not found

b_c <- plot_linear_scatter(hs_macr_combat_batch, pretty_colors=FALSE)

## Error in data.frame(df[, c(1, 2)]): object 'hs_macr_combat_batch' not found

b_c$scatter

## Error in eval(expr, envir, enclos): object 'b_c' not found

b_c$correlation

## Error in eval(expr, envir, enclos): object 'b_c' not found

6.7 Compare no batch to batch in model, limma

hs_macr_nobatch_batch <- merge(
  hs_macr_nobatch$limma$all_tables$sh_vs_chr,
  hs_macr_batch$limma$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_nobatch_batch) <- hs_macr_nobatch_batch[["Row.names"]]
hs_macr_nobatch_batch <- hs_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(hs_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- plot_linear_scatter(hs_macr_nobatch_batch, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 79, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6713 0.6962
## sample estimates:
##    cor 
## 0.6839

6.8 Batch in model vs. SVA, limma

hs_macr_batch_sva <- merge(
  hs_macr_batch$limma$all_tables$sh_vs_chr,
  hs_macr_sva$limma$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_sva) <- hs_macr_batch_sva[["Row.names"]]
hs_macr_batch_sva <- hs_macr_batch_sva[, c("logFC.x","logFC.y")]
colnames(hs_macr_batch_sva) <- c("batch","sva")
b_s <- plot_linear_scatter(hs_macr_batch_sva, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 460, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9830 0.9845
## sample estimates:
##    cor 
## 0.9837

6.9 Batch in model vs. combat, limma

hs_macr_batch_combat <- merge(
  hs_macr_batch$limma$all_tables$sh_vs_chr,
  hs_macr_combat$limma$all_tables$sh_vs_chr,
  by="row.names")

## Error in as.data.frame(y): object 'hs_macr_combat' not found

rownames(hs_macr_batch_combat) <- hs_macr_batch_combat[["Row.names"]]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

hs_macr_batch_combat <- hs_macr_batch_combat[, c("logFC.x","logFC.y")]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

colnames(hs_macr_batch_combat) <- c("batch","combat")

## Error in colnames(hs_macr_batch_combat) <- c("batch", "combat"): object 'hs_macr_batch_combat' not found

b_c <- plot_linear_scatter(hs_macr_batch_combat, pretty_colors=FALSE)

## Error in data.frame(df[, c(1, 2)]): object 'hs_macr_batch_combat' not found

b_c$scatter

## Error in eval(expr, envir, enclos): object 'b_c' not found

b_c$correlation

## Error in eval(expr, envir, enclos): object 'b_c' not found

6.10 Nobatch vs. batch in model, edger

hs_macr_nobatch_batch <- merge(
  hs_macr_nobatch$edger$all_tables$sh_vs_chr,
  hs_macr_batch$edger$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_nobatch_batch) <- hs_macr_nobatch_batch[["Row.names"]]
hs_macr_nobatch_batch <- hs_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(hs_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- sm(plot_linear_scatter(hs_macr_nobatch_batch, pretty_colors=FALSE))
nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 100, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7671 0.7857
## sample estimates:
##    cor 
## 0.7765

6.11 Batch in model vs. SVA, edger

hs_macr_batch_sva <- merge(
  hs_macr_batch$edger$all_tables$sh_vs_chr,
  hs_macr_sva$edger$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_sva) <- hs_macr_batch_sva[["Row.names"]]
hs_macr_batch_sva <- hs_macr_batch_sva[, c("logFC.x","logFC.y")]
colnames(hs_macr_batch_sva) <- c("batch","sva")
b_s <- plot_linear_scatter(hs_macr_batch_sva, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 340, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9698 0.9725
## sample estimates:
##    cor 
## 0.9712

6.12 Batch in model vs. combat, edger

hs_macr_batch_combat <- merge(
  hs_macr_batch$edger$all_tables$sh_vs_chr,
  hs_macr_combat$edger$all_tables$sh_vs_chr,
  by="row.names")

## Error in as.data.frame(y): object 'hs_macr_combat' not found

rownames(hs_macr_batch_combat) <- hs_macr_batch_combat[["Row.names"]]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

hs_macr_batch_combat <- hs_macr_batch_combat[, c("logFC.x","logFC.y")]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

colnames(hs_macr_batch_combat) <- c("batch","combat")

## Error in colnames(hs_macr_batch_combat) <- c("batch", "combat"): object 'hs_macr_batch_combat' not found

b_c <- plot_linear_scatter(hs_macr_batch_combat, pretty_colors=FALSE)

## Error in data.frame(df[, c(1, 2)]): object 'hs_macr_batch_combat' not found

b_c$scatter

## Error in eval(expr, envir, enclos): object 'b_c' not found

b_c$correlation

## Error in eval(expr, envir, enclos): object 'b_c' not found

6.13 Compare nobatch vs. batch, deseq

hs_macr_nobatch_batch <- merge(
  hs_macr_nobatch$deseq$all_tables$sh_vs_chr,
  hs_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_nobatch_batch) <- hs_macr_nobatch_batch[["Row.names"]]
hs_macr_nobatch_batch <- hs_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(hs_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- sm(plot_linear_scatter(hs_macr_nobatch_batch, pretty_colors=FALSE))
nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 79, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6727 0.6975
## sample estimates:
##    cor 
## 0.6853

6.14 Compare batch vs. SVA, deseq

hs_macr_batch_sva <- merge(
  hs_macr_batch$deseq$all_tables$sh_vs_chr,
  hs_macr_sva$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(hs_macr_batch_sva) <- hs_macr_batch_sva[["Row.names"]]
hs_macr_batch_sva <- hs_macr_batch_sva[, c("logFC.x", "logFC.y")]
colnames(hs_macr_batch_sva) <- c("batch", "sva")
b_s <- sm(plot_linear_scatter(hs_macr_batch_sva, pretty_colors=FALSE))
b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 1000, df = 7000, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9967 0.9970
## sample estimates:
##    cor 
## 0.9968

6.15 Batch in model vs. combat, deseq

hs_macr_batch_combat <- merge(
  hs_macr_batch$deseq$all_tables$sh_vs_chr,
  hs_macr_combat$deseq$all_tables$sh_vs_chr,
  by="row.names")

## Error in as.data.frame(y): object 'hs_macr_combat' not found

rownames(hs_macr_batch_combat) <- hs_macr_batch_combat[["Row.names"]]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

hs_macr_batch_combat <- hs_macr_batch_combat[, c("logFC.x", "logFC.y")]

## Error in eval(expr, envir, enclos): object 'hs_macr_batch_combat' not found

colnames(hs_macr_batch_combat) <- c("batch", "combat")

## Error in colnames(hs_macr_batch_combat) <- c("batch", "combat"): object 'hs_macr_batch_combat' not found

b_c <- sm(plot_linear_scatter(hs_macr_batch_combat, pretty_colors=FALSE))

## Error in data.frame(df[, c(1, 2)]): object 'hs_macr_batch_combat' not found

b_c$scatter

## Error in eval(expr, envir, enclos): object 'b_c' not found

b_c$correlation

## Error in eval(expr, envir, enclos): object 'b_c' not found

7.1 No batch in the model

lp_macr_nobatch <- sm(all_pairwise(lp_macr_lowfilt, limma_method="robust", model_batch=FALSE))
## wow, all tools including basic agree almost completely
medians_by_condition <- lp_macr_nobatch$basic$medians
lp_macr_nobatch_tables <- sm(combine_de_tables(
  lp_macr_nobatch,
  excel=paste0("excel/lp_macr_nobatch-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_nobatch_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.2 Batch in the model

In this attempt, we add a batch factor in the experimental model and see how it does.

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
lp_macr_batch <- sm(all_pairwise(lp_macr_lowfilt, limma_method="robust"))
medians_by_condition <- lp_macr_batch$basic$medians
lp_macr_batch_tables <- sm(combine_de_tables(
  lp_macr_batch,
  excel=paste0("excel/lp_macr_batchmodel-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_batchmodel_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.3 Batch estimated with SVA

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
lp_macr_sva <- sm(all_pairwise(lp_macr_lowfilt, limma_method="robust", model_batch="sva"))
medians_by_condition <- lp_macr_sva$basic$medians
lp_macr_sva_tables <- sm(combine_de_tables(
  lp_macr_sva,
  excel=paste0("excel/lp_macr_sva-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_sva_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.4 Batch correction via ruv residuals

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
## Bizarrely, sometimes if one runs this, it gives an error "Error in (function
## (classes, fdef, mtable) : unable to find an inherited method for function
## 'RUVr' for signature '"matrix", "logical", "numeric", "NULL"'"  -- however,
## if one then simply runs it again it works fine.
## I changed model_surrogates.R to take care of this I think.
lp_macr_ruvres <- sm(all_pairwise(lp_macr_lowfilt, model_batch="ruv_residuals", limma_method="robust"))
medians_by_condition <- lp_macr_ruvres$basic$medians
lp_macr_ruvres_tables <- sm(combine_de_tables(
  lp_macr_ruvres,
  excel=paste0("excel/lp_macr_ruvres-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_ruvres_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.5 Batch correction with pca

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
lp_macr_pca <- sm(all_pairwise(lp_macr_lowfilt, model_batch="pca", limma_method="robust"))
medians_by_condition <- lp_macr_pca$basic$medians
lp_macr_pca_tables <- sm(combine_de_tables(
  lp_macr_pca,
  excel=paste0("excel/lp_macr_pca-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_pca_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.6 Batch correction with ruv empirical

## Here just let all_pairwise run on filtered data and do its normal ~ 0 + condition + batch analyses
lp_macr_ruvemp <- sm(all_pairwise(lp_macr_lowfilt, model_batch="ruv_empirical", limma_method="robust"))
medians_by_condition <- lp_macr_ruvemp$basic$medians
lp_macr_ruvemp_tables <- sm(combine_de_tables(
  lp_macr_ruvemp,
  excel=paste0("excel/lp_macr_ruvemp-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_ruvemp_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

7.7 Batch correction with combat

Then repeat with the batch-corrected data and see the differences.

lp_macr_combat <- all_pairwise(lp_macr_combat_norm, force=TRUE,
                               limma_method="robust", parallel=FALSE)

## Using limma's removeBatchEffect to visualize before/after batch inclusion.

## Starting basic pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating median and variance tables.

## Basic step 2/3: Performing 3 comparisons.

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## Starting DESeq2 pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you,
## like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Choosing the non-intercept containing model.

## DESeq2 step 1/5: Including batch and condition in the deseq model.

## converting counts to integer mode

## DESeq2 step 2/5: Estimate size factors.

## DESeq2 step 3/5: Estimate dispersions.

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Using a parametric fitting seems to have worked.

## DESeq2 step 4/5: nbinomWaldTest.

## Plotting dispersions.

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

## EdgeR/DESeq expect raw data as input, reverting to the count filtered data.

## Starting EBSeq pairwise subset.

## Starting EBTest of chr vs. sh.

## Starting edgeR pairwise comparisons.

## About to round the data, this is a pretty terrible thing to do. But if you,
## like me, want to see what happens when you put non-standard data into deseq, then here you go.

## Warning in choose_binom_dataset(input, force = force): This data was
## inappropriately forced into integers.

## Choosing the non-intercept containing model.

## EdgeR step 1/9: importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

## Starting limma pairwise comparison.

## Leaving the data alone, regardless of normalization state.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## Choosing the non-intercept containing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method=quantile for voom.

## Limma step 3/6: running lmFit with method: robust.

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Warning in rlm.default(x = X, y = y, weights = w, ...): 'rlm' failed to
## converge in 20 steps

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Limma step 5/6: Running eBayes with robust=FALSE and trend=FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/1: Creating table: sh_vs_chr.  Adjust=BH

## Limma step 6/6: 1/2: Creating table: chr.  Adjust=BH

## Limma step 6/6: 2/2: Creating table: sh.  Adjust=BH

## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr
## Comparing analyses 1/1: sh_vs_chr

medians_by_condition <- lp_macr_combat$basic$medians
lp_macr_combat_tables <- sm(combine_de_tables(
  lp_macr_combat,
  excel=paste0("excel/lp_macr_combat-v", ver, ".xlsx"),
  sig_excel=paste0("excel/lp_macr_combat_significant-v", ver, ".xlsx"),
  keepers=lp_contrasts,
  extra_annot=medians_by_condition))

## Error in .subset2(x, i, exact = exact): subscript out of bounds

9.1 Compare DeSeq / Basic without batch in model

lp_macr_nobatch_basic <- merge(
  lp_macr_nobatch$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_nobatch_basic) <- lp_macr_nobatch_basic[["Row.names"]]
lp_macr_nobatch_logfc <- lp_macr_nobatch_basic[, c("logFC.x", "logFC.y")]
colnames(lp_macr_nobatch_logfc) <- c("nobatch","basic")
lfc_nb_b <- sm(plot_linear_scatter(lp_macr_nobatch_logfc, pretty_colors=FALSE))
lfc_nb_b$scatter

lfc_nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 89, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6805 0.7025
## sample estimates:
##    cor 
## 0.6917

lp_macr_nobatch_p <- lp_macr_nobatch_basic[, c("P.Value","p")]
lp_macr_nobatch_p[[2]] <- as.numeric(lp_macr_nobatch_p[[2]])
colnames(lp_macr_nobatch_p) <- c("nobatch","basic")
lp_macr_nobatch_p <- -1 * log(lp_macr_nobatch_p)
p_nb_b <- sm(plot_linear_scatter(lp_macr_nobatch_p, pretty_colors=FALSE))
p_nb_b$scatter

p_nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 120, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7760 0.7922
## sample estimates:
##    cor 
## 0.7842

9.2 Compare SVA to batch in model, DESeq

lp_macr_sva_batch <- merge(
  lp_macr_sva$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_sva_batch) <- lp_macr_sva_batch[["Row.names"]]
lp_macr_sva_logfc <- lp_macr_sva_batch[, c("logFC.x","logFC.y")]
colnames(lp_macr_sva_logfc) <- c("sva","batch")
lfc_b_s <- sm(plot_linear_scatter(lp_macr_sva_logfc, pretty_colors=FALSE))
lfc_b_s$scatter

lfc_b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 55, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4939 0.5251
## sample estimates:
##    cor 
## 0.5097

lp_macr_sva_p <- lp_macr_sva_batch[, c("P.Value.x","P.Value.y")]
lp_macr_sva_p[[2]] <- as.numeric(lp_macr_sva_p[[2]])
colnames(lp_macr_sva_p) <- c("sva","batch")
lp_macr_sva_p <- -1 * log(lp_macr_sva_p)
p_b_s <- sm(plot_linear_scatter(lp_macr_sva_p, pretty_colors=FALSE))
p_b_s$scatter

p_b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 25, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2381 0.2774
## sample estimates:
##    cor 
## 0.2579

lfcp_b_s <- lp_macr_sva_batch[, c("logFC.x", "logFC.y", "P.Value.x", "P.Value.y")]
colnames(lfcp_b_s) <- c("l2fcsva", "l2fcbatch", "psva", "pbatch")
lfc_b_s$scatter

cutoff <- 0.1
lfcp_b_s$state <- ifelse(lfcp_b_s$psva > cutoff & lfcp_b_s$pbatch > cutoff, "bothinsig",
                  ifelse(lfcp_b_s$psva <= cutoff & lfcp_b_s$pbatch <= cutoff, "bothsig",
                  ifelse(lfcp_b_s$psva <= cutoff, "svasig", "batchsig")))
##lfcp_b_s$lfcstate <- ifelse(lfcp_b_s$l2fcsva >= 0.75 & lfcp_b_s$l2fcbatch, "", "")
num_bothinsig <- sum(lfcp_b_s$state == "bothinsig")
num_bothsig <- sum(lfcp_b_s$state == "bothsig")
num_svasig <- sum(lfcp_b_s$state == "svasig")
num_batchsig <- sum(lfcp_b_s$state == "batchsig")

library(ggplot2)
aes_color = "(l2fcsva >= 0.75 | l2fcsva <= -0.75 | l2fcbatch >= 0.75 | l2fcbatch <= -0.75)"

plt <- ggplot2::ggplot(lfcp_b_s, aes_string(x="l2fcsva", y="l2fcbatch")) +
    ## ggplot2::geom_point(stat="identity", size=2, alpha=0.2, aes_string(shape="as.factor(aes_color)", colour="as.factor(state)", fill="as.factor(state)")) +
    ggplot2::geom_abline(colour="blue", slope=1, intercept=0, size=0.5) +
    ggplot2::geom_hline(yintercept=c(-0.75, 0.75), color="red", size=0.5) +
    ggplot2::geom_vline(xintercept=c(-0.75, 0.75), color="red", size=0.5) +
    ggplot2::geom_point(stat="identity", size=2, alpha=0.2, aes_string(colour="as.factor(state)", fill="as.factor(state)")) +
    ggplot2::scale_color_manual(name="state", values=c("bothinsig"="grey", "bothsig"="forestgreen", "svasig"="darkred", "batchsig"="darkblue")) +
    ggplot2::scale_fill_manual(name="state", values=c("bothinsig"="grey", "bothsig"="forestgreen", "svasig"="darkred", "batchsig"="darkblue"),
                               labels=c(
                                   paste0("Both InSig.: ", num_bothinsig),
                                   paste0("Both Sig.: ", num_bothsig),
                                   paste0("Sva Sig.: ", num_svasig),
                                   paste0("Batch Sig.: ", num_batchsig)),
                               guide=ggplot2::guide_legend(override.aes=aes(size=3, fill="grey"))) +
    ggplot2::guides(fill=ggplot2::guide_legend(override.aes=list(size=3))) +
    ggplot2::theme_bw()
plt

9.3 Compare ruvresid to batch in model, DESeq

lp_macr_batch_ruvresid_deseq <- merge(
  lp_macr_ruvres$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_ruvresid_deseq) <- lp_macr_batch_ruvresid_deseq[["Row.names"]]
lp_macr_batch_ruvresid_logfc <- lp_macr_batch_ruvresid_deseq[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_ruvresid_logfc) <- c("nobatch","basic")
lfc_ruv_bat <- plot_linear_scatter(lp_macr_batch_ruvresid_logfc, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

lfc_ruv_bat$scatter

lfc_ruv_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 130, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7949 0.8099
## sample estimates:
##    cor 
## 0.8025

9.4 Compare pca to batch in model, DESeq

lp_macr_batch_pca_deseq <- merge(
  lp_macr_pca$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_pca_deseq) <- lp_macr_batch_pca_deseq[["Row.names"]]
lp_macr_batch_pca_logfc <- lp_macr_batch_pca_deseq[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_pca_logfc) <- c("nobatch","basic")
lfc_pca_bat <- plot_linear_scatter(lp_macr_batch_pca_logfc, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

lfc_pca_bat$scatter

lfc_pca_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 64, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5499 0.5786
## sample estimates:
##    cor 
## 0.5644

9.5 Compare ruv empirical to batch in model, DESeq

lp_macr_batch_ruvemp_deseq <- merge(
  lp_macr_ruvemp$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$basic$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_ruvemp_deseq) <- lp_macr_batch_ruvemp_deseq[["Row.names"]]
lp_macr_batch_ruvemp_logfc <- lp_macr_batch_ruvemp_deseq[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_ruvemp_logfc) <- c("nobatch","basic")
lfc_ruvemp_bat <- sm(plot_linear_scatter(lp_macr_batch_ruvemp_logfc, pretty_colors=FALSE))
lfc_ruvemp_bat$scatter

lfc_ruvemp_bat$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 100, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7289 0.7480
## sample estimates:
##    cor 
## 0.7386

9.6 Compare combat to batch in model, DESeq

lp_macr_combat_batch <- merge(
  lp_macr_combat$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_combat_batch) <- lp_macr_combat_batch[["Row.names"]]
lp_macr_combat_batch <- lp_macr_combat_batch[, c("logFC.x","logFC.y")]
colnames(lp_macr_combat_batch) <- c("batch","combat")
b_c <- plot_linear_scatter(lp_macr_combat_batch, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_c$scatter

b_c$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 150, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8500 0.8613
## sample estimates:
##    cor 
## 0.8558

9.7 Compare no batch to batch in model, limma

lp_macr_nobatch_batch <- merge(
  lp_macr_nobatch$limma$all_tables$sh_vs_chr,
  lp_macr_batch$limma$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_nobatch_batch) <- lp_macr_nobatch_batch[["Row.names"]]
lp_macr_nobatch_batch <- lp_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(lp_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- plot_linear_scatter(lp_macr_nobatch_batch, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 150, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8441 0.8558
## sample estimates:
##    cor 
## 0.8501

9.8 Batch in model vs. SVA, limma

lp_macr_batch_sva <- merge(
  lp_macr_batch$limma$all_tables$sh_vs_chr,
  lp_macr_sva$limma$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_sva) <- lp_macr_batch_sva[["Row.names"]]
lp_macr_batch_sva <- lp_macr_batch_sva[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_sva) <- c("batch","sva")
b_s <- plot_linear_scatter(lp_macr_batch_sva, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 45, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4216 0.4556
## sample estimates:
##    cor 
## 0.4387

9.9 Batch in model vs. combat, limma

lp_macr_batch_combat <- merge(
  lp_macr_batch$limma$all_tables$sh_vs_chr,
  lp_macr_combat$limma$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_combat) <- lp_macr_batch_combat[["Row.names"]]
lp_macr_batch_combat <- lp_macr_batch_combat[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_combat) <- c("batch","combat")
b_c <- plot_linear_scatter(lp_macr_batch_combat, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_c$scatter

b_c$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 360, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9663 0.9690
## sample estimates:
##    cor 
## 0.9677

9.10 Nobatch vs. batch in model, edger

lp_macr_nobatch_batch <- merge(
  lp_macr_nobatch$edger$all_tables$sh_vs_chr,
  lp_macr_batch$edger$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_nobatch_batch) <- lp_macr_nobatch_batch[["Row.names"]]
lp_macr_nobatch_batch <- lp_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(lp_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- sm(plot_linear_scatter(lp_macr_nobatch_batch, pretty_colors=FALSE))
nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 100, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7371 0.7558
## sample estimates:
##    cor 
## 0.7466

9.11 Batch in model vs. SVA, edger

lp_macr_batch_sva <- merge(
  lp_macr_batch$edger$all_tables$sh_vs_chr,
  lp_macr_sva$edger$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_sva) <- lp_macr_batch_sva[["Row.names"]]
lp_macr_batch_sva <- lp_macr_batch_sva[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_sva) <- c("batch","sva")
b_s <- plot_linear_scatter(lp_macr_batch_sva, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 56, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4966 0.5277
## sample estimates:
##    cor 
## 0.5123

9.12 Batch in model vs. combat, edger

lp_macr_batch_combat <- merge(
  lp_macr_batch$edger$all_tables$sh_vs_chr,
  lp_macr_combat$edger$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_combat) <- lp_macr_batch_combat[["Row.names"]]
lp_macr_batch_combat <- lp_macr_batch_combat[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_combat) <- c("batch","combat")
b_c <- plot_linear_scatter(lp_macr_batch_combat, pretty_colors=FALSE)

## Warning in plot_multihistogram(df): NAs introduced by coercion

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

b_c$scatter

b_c$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 150, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8500 0.8613
## sample estimates:
##    cor 
## 0.8558

9.13 Compare nobatch vs. batch, deseq

lp_macr_nobatch_batch <- merge(
  lp_macr_nobatch$deseq$all_tables$sh_vs_chr,
  lp_macr_batch$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_nobatch_batch) <- lp_macr_nobatch_batch[["Row.names"]]
lp_macr_nobatch_batch <- lp_macr_nobatch_batch[, c("logFC.x","logFC.y")]
colnames(lp_macr_nobatch_batch) <- c("nobatch","batch")
nb_b <- sm(plot_linear_scatter(lp_macr_nobatch_batch, pretty_colors=FALSE))
nb_b$scatter

nb_b$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 110, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7400 0.7585
## sample estimates:
##    cor 
## 0.7494

9.14 Compare batch vs. SVA, deseq

lp_macr_batch_sva <- merge(
  lp_macr_batch$deseq$all_tables$sh_vs_chr,
  lp_macr_sva$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_sva) <- lp_macr_batch_sva[["Row.names"]]
lp_macr_batch_sva <- lp_macr_batch_sva[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_sva) <- c("batch","sva")
b_s <- sm(plot_linear_scatter(lp_macr_batch_sva, pretty_colors=FALSE))
b_s$scatter

b_s$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 55, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4939 0.5251
## sample estimates:
##    cor 
## 0.5097

9.15 Batch in model vs. combat, deseq

lp_macr_batch_combat <- merge(
  lp_macr_batch$deseq$all_tables$sh_vs_chr,
  lp_macr_combat$deseq$all_tables$sh_vs_chr,
  by="row.names")
rownames(lp_macr_batch_combat) <- lp_macr_batch_combat[["Row.names"]]
lp_macr_batch_combat <- lp_macr_batch_combat[, c("logFC.x","logFC.y")]
colnames(lp_macr_batch_combat) <- c("batch","combat")
b_c <- sm(plot_linear_scatter(lp_macr_batch_combat, pretty_colors=FALSE))
b_c$scatter

b_c$correlation

## 
##  Pearson's product-moment correlation
## 
## data:  df[, 1] and df[, 2]
## t = 150, df = 8700, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8500 0.8613
## sample estimates:
##    cor 
## 0.8558

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: ggplot2(v.3.0.0), edgeR(v.3.22.3), ruv(v.0.9.7) and hpgltools(v.2018.03)

loaded via a namespace (and not attached): backports(v.1.1.2), Hmisc(v.4.1-1), aroma.light(v.3.10.0), plyr(v.1.8.4), lazyeval(v.0.2.1), splines(v.3.5.1), BiocParallel(v.1.14.2), GenomeInfoDb(v.1.16.0), sva(v.3.28.0), digest(v.0.6.15), foreach(v.1.4.4), BiocInstaller(v.1.30.0), htmltools(v.0.3.6), gdata(v.2.18.0), magrittr(v.1.5), checkmate(v.1.8.5), memoise(v.1.1.0), cluster(v.2.0.7-1), doParallel(v.1.0.11), openxlsx(v.4.1.0), limma(v.3.36.2), Biostrings(v.2.48.0), annotate(v.1.58.0), matrixStats(v.0.54.0), R.utils(v.2.6.0), blockmodeling(v.0.3.1), prettyunits(v.1.0.2), colorspace(v.1.3-2), blob(v.1.1.1), ggrepel(v.0.8.0), dplyr(v.0.7.6), crayon(v.1.3.4), RCurl(v.1.95-4.11), graph(v.1.58.0), roxygen2(v.6.1.0), genefilter(v.1.62.0), lme4(v.1.1-17), bindr(v.0.1.1), survival(v.2.42-6), iterators(v.1.0.10), glue(v.1.3.0), registry(v.0.5), gtable(v.0.2.0), zlibbioc(v.1.26.0), XVector(v.0.20.0), DelayedArray(v.0.6.4), DEoptimR(v.1.0-8), BiocGenerics(v.0.26.0), DESeq(v.1.32.0), scales(v.1.0.0), rngtools(v.1.3.1), DBI(v.1.0.0), bibtex(v.0.4.2), Rcpp(v.0.12.18), xtable(v.1.8-2), progress(v.1.2.0), htmlTable(v.1.12), foreign(v.0.8-71), bit(v.1.1-14), OrganismDbi(v.1.22.0), preprocessCore(v.1.42.0), Formula(v.1.2-3), EBSeq(v.1.20.0), stats4(v.3.5.1), htmlwidgets(v.1.2), httr(v.1.3.1), gplots(v.3.0.1), RColorBrewer(v.1.1-2), acepack(v.1.4.1), R.methodsS3(v.1.7.1), pkgconfig(v.2.0.1), XML(v.3.98-1.15), nnet(v.7.3-12), locfit(v.1.5-9.1), tidyselect(v.0.2.4), labeling(v.0.3), rlang(v.0.2.1), reshape2(v.1.4.3), AnnotationDbi(v.1.42.1), munsell(v.0.5.0), tools(v.3.5.1), RSQLite(v.2.1.1), devtools(v.1.13.6), evaluate(v.0.11), stringr(v.1.3.1), yaml(v.2.2.0), knitr(v.1.20), bit64(v.0.9-7), pander(v.0.6.2), robustbase(v.0.93-2), zip(v.1.0.0), caTools(v.1.17.1.1), purrr(v.0.2.5), bindrcpp(v.0.2.2), EDASeq(v.2.14.1), doRNG(v.1.7.1), RBGL(v.1.56.0), nlme(v.3.1-137), R.oo(v.1.22.0), xml2(v.1.2.0), biomaRt(v.2.36.1), compiler(v.3.5.1), pbkrtest(v.0.4-7), rstudioapi(v.0.7), testthat(v.2.0.0), variancePartition(v.1.10.0), statmod(v.1.4.30), tibble(v.1.4.2), geneplotter(v.1.58.0), stringi(v.1.2.4), GenomicFeatures(v.1.32.1), lattice(v.0.20-35), Matrix(v.1.2-14), commonmark(v.1.5), nloptr(v.1.0.4), pillar(v.1.3.0), data.table(v.1.11.4), bitops(v.1.0-6), corpcor(v.1.6.9), rtracklayer(v.1.40.4), GenomicRanges(v.1.32.6), colorRamps(v.2.3), hwriter(v.1.3.2), R6(v.2.2.2), latticeExtra(v.0.6-28), directlabels(v.2018.05.22), ShortRead(v.1.38.0), KernSmooth(v.2.23-15), gridExtra(v.2.3), IRanges(v.2.14.10), codetools(v.0.2-15), MASS(v.7.3-50), gtools(v.3.8.1), assertthat(v.0.2.0), SummarizedExperiment(v.1.10.1), pkgmaker(v.0.27), DESeq2(v.1.20.0), rprojroot(v.1.3-2), RUVSeq(v.1.14.0), withr(v.2.1.2), GenomicAlignments(v.1.16.0), Rsamtools(v.1.32.2), S4Vectors(v.0.18.3), GenomeInfoDbData(v.1.1.0), mgcv(v.1.8-24), parallel(v.3.5.1), hms(v.0.4.2), quadprog(v.1.5-5), grid(v.3.5.1), rpart(v.4.1-13), minqa(v.1.2.4), rmarkdown(v.1.10), Biobase(v.2.40.0) and base64enc(v.0.1-3)

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 8a4df418d36491bab7b987cd21768b3d7cab1269

## R> packrat::restore()

## This is hpgltools commit: Thu Aug 16 10:45:44 2018 -0400: 8a4df418d36491bab7b987cd21768b3d7cab1269

this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
message(paste0("Saving to ", this_save))

## Saving to 03_expression_macrophage-v20170820.rda.xz

tmp <- sm(saveme(filename=this_save))