CL-Brener specific material with Fernanda

Pulling together the various strands of material used in Fernanda’s manuscript. I think this should prove to just be a repetition of work already performed by Edson.

1 Collect annotation data

Actually, since this was performed with a previous release of the T.cruzi genome, it is probably smarter for me to load the annotations from my earlier savefile. Because the tritrypdb does not save the older annotation data with new releases, there is no guarantee that I will match the annotations completely if I make my own (they do provide a mapping of changed IDs which seems to work, but why waste the time messing with it when I still have the saved annotations from when I did it originally).

ver <- "20170810"
previous_file <- "01_annotation.Rmd"
tmp <- try(sm(loadme(filename=paste0(gsub(pattern="\\.Rmd", replace="", x=previous_file), "-v", ver, ".rda.xz"))))

2 Create a CL-Brener exclusive data set excluding the problematic sample

all_expt <- sm(create_expt(metadata="sample_sheets/cl14clbr_samples_combined_all.xlsx",
                           file_column="clbrenerfile", gene_info=all_genes))
all_expt <- exclude_genes_expt(all_expt)

## Before removal, there were 25100 entries.

## Now there are 23305 entries.

## Percent kept: 99.781, 99.815, 99.814, 99.818, 97.848, 98.584, 99.776, 99.449, 99.331, 99.457, 99.321, 99.395, 99.157, 98.435, 98.543, 99.306, 99.396, 99.653, 99.663, 99.249, 99.436, 98.720, 99.886, 99.761

## Percent removed: 0.219, 0.185, 0.186, 0.182, 2.152, 1.416, 0.224, 0.551, 0.669, 0.543, 0.679, 0.605, 0.843, 1.565, 1.457, 0.694, 0.604, 0.347, 0.337, 0.751, 0.564, 1.280, 0.114, 0.239

outlier_subset <- "sampleid!='HPGL0490'"
all_minus_one <- subset_expt(all_expt, subset=outlier_subset)

## Using a subset expression.

## There were 24, now there are 23 samples.

all_minus_one <- subset_expt(all_minus_one, subset="type!='CL14'")

## Using a subset expression.

## There were 23, now there are 11 samples.

all_minus_cl14 <- subset_expt(all_expt, subset="type!='CL14'")

## Using a subset expression.

## There were 24, now there are 12 samples.

tmp_expt <- sm(normalize_expt(all_minus_cl14, filter=TRUE, convert="cpm", norm="quant", transform="log2"))
pp(file="images/all_clbr_pca.pdf", image=plot_pca(tmp_expt)$plot)

## Writing the image to: images/all_clbr_pca.pdf and calling dev.off().

tmp_expt <- sm(normalize_expt(all_miOAnus_one, filter=TRUE, convert="cpm", norm="quant", transform="log2"))

## Error in normalize_expt(all_miOAnus_one, filter = TRUE, convert = "cpm", : object 'all_miOAnus_one' not found

pp(file="images/all_clbr_minus_one_pca.pdf", image=plot_pca(tmp_expt)$plot)

## Writing the image to: images/all_clbr_minus_one_pca.pdf and calling dev.off().

3 Perform a DE analysis

all_minus_norm <- normalize_expt(all_minus_one, filter=TRUE)

## This function will replace the expt$expressionset slot with:

## cbcb(data)

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unconverted.  It is often advisable to cpm/rpkm
##  the data to normalize for sampling differences, keep in mind though that rpkm
##  has some annoying biases, and voom() by default does a cpm (though hpgl_voom()
##  will try to detect this).

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 1245 low-count genes (22060 remaining).

## Step 2: not normalizing the data.

## Step 3: not converting the data.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

all_minus_de <- all_pairwise(all_minus_norm, model_batch="svaseq")

## batch_counts: Before batch/surrogate estimation, 241564 entries are x>1: 100%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

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

## Attempting svaseq estimation with 1 surrogates.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (22060 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 396 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 237917 entries are x>1: 98%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 4347 entries are 0<x<1: 2%.

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

## Attempting svaseq estimation with 2 surrogates.

## There are 44 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
    "epi_over_t96h" = c("CLBrEpi", "CLBrA96"),
    "epi_over_t60h" = c("CLBrEpi", "CLBrA60"),
    "epi_over_tryp" = c("CLBrEpi", "CLBrTryp"))
all_minus_tables <- combine_de_tables(
    all_minus_de, keepers=keepers,
    excel=glue::glue("excel/clbr_analyses-v{ver}.xlsx"))

## Deleting the file excel/clbr_analyses-v20170810.xlsx before writing the tables.

## Writing a legend of columns.

## Printing a pca plot before/after surrogates/batch estimation.

## Working on 1/3: epi_over_t96h which is: CLBrEpi/CLBrA96.

## Found table with CLBrEpi_vs_CLBrA96

## Working on 2/3: epi_over_t60h which is: CLBrEpi/CLBrA60.

## Found table with CLBrEpi_vs_CLBrA60

## Working on 3/3: epi_over_tryp which is: CLBrEpi/CLBrTryp.

## Found inverse table with CLBrTryp_vs_CLBrEpi

## Adding venn plots for epi_over_t96h.

## Limma expression coefficients for epi_over_t96h; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_t96h; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_t96h; R^2: 0.894; equation: y = 0.839x + 1.31

## Adding venn plots for epi_over_t60h.

## Limma expression coefficients for epi_over_t60h; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_t60h; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_t60h; R^2: 0.894; equation: y = 0.839x + 1.31

## Adding venn plots for epi_over_tryp.

## Limma expression coefficients for epi_over_tryp; R^2: 0.895; equation: y = 0.971x + 0.0236

## Deseq expression coefficients for epi_over_tryp; R^2: 0.894; equation: y = 0.838x + 1.33

## Edger expression coefficients for epi_over_tryp; R^2: 0.894; equation: y = 0.839x + 1.31

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/clbr_analyses-v20170810.xlsx.

all_minus_sig <- extract_significant_genes(
    all_minus_tables,
    excel=glue::glue("excel/clbr_analyses_sig-v{ver}.xlsx"))

## Writing a legend of columns.

## Writing excel data according to limma for epi_over_t96h: 1/15.

## After (adj)p filter, the up genes table has 6602 genes.

## After (adj)p filter, the down genes table has 6922 genes.

## After fold change filter, the up genes table has 2081 genes.

## After fold change filter, the down genes table has 2132 genes.

## Writing excel data according to limma for epi_over_t60h: 2/15.

## After (adj)p filter, the up genes table has 4542 genes.

## After (adj)p filter, the down genes table has 6896 genes.

## After fold change filter, the up genes table has 1756 genes.

## After fold change filter, the down genes table has 838 genes.

## Writing excel data according to limma for epi_over_tryp: 3/15.

## After (adj)p filter, the up genes table has 9459 genes.

## After (adj)p filter, the down genes table has 7640 genes.

## After fold change filter, the up genes table has 3339 genes.

## After fold change filter, the down genes table has 2836 genes.

## Printing significant genes to the file: excel/clbr_analyses_sig-v20170810.xlsx

## 1/3: Creating significant table up_limma_epi_over_t96h

## 2/3: Creating significant table up_limma_epi_over_t60h

## 3/3: Creating significant table up_limma_epi_over_tryp

## Writing excel data according to edger for epi_over_t96h: 1/15.

## After (adj)p filter, the up genes table has 6168 genes.

## After (adj)p filter, the down genes table has 7462 genes.

## After fold change filter, the up genes table has 1776 genes.

## After fold change filter, the down genes table has 2617 genes.

## Writing excel data according to edger for epi_over_t60h: 2/15.

## After (adj)p filter, the up genes table has 6286 genes.

## After (adj)p filter, the down genes table has 5885 genes.

## After fold change filter, the up genes table has 2484 genes.

## After fold change filter, the down genes table has 664 genes.

## Writing excel data according to edger for epi_over_tryp: 3/15.

## After (adj)p filter, the up genes table has 7790 genes.

## After (adj)p filter, the down genes table has 8071 genes.

## After fold change filter, the up genes table has 2308 genes.

## After fold change filter, the down genes table has 3559 genes.

## Printing significant genes to the file: excel/clbr_analyses_sig-v20170810.xlsx

## 1/3: Creating significant table up_edger_epi_over_t96h

## 2/3: Creating significant table up_edger_epi_over_t60h

## 3/3: Creating significant table up_edger_epi_over_tryp

## Writing excel data according to deseq for epi_over_t96h: 1/15.

## After (adj)p filter, the up genes table has 6368 genes.

## After (adj)p filter, the down genes table has 6246 genes.

## After fold change filter, the up genes table has 1921 genes.

## After fold change filter, the down genes table has 2471 genes.

## Writing excel data according to deseq for epi_over_t60h: 2/15.

## After (adj)p filter, the up genes table has 6105 genes.

## After (adj)p filter, the down genes table has 5255 genes.

## After fold change filter, the up genes table has 2489 genes.

## After fold change filter, the down genes table has 624 genes.

## Writing excel data according to deseq for epi_over_tryp: 3/15.

## After (adj)p filter, the up genes table has 7638 genes.

## After (adj)p filter, the down genes table has 7387 genes.

## After fold change filter, the up genes table has 2404 genes.

## After fold change filter, the down genes table has 3525 genes.

## Printing significant genes to the file: excel/clbr_analyses_sig-v20170810.xlsx

## 1/3: Creating significant table up_deseq_epi_over_t96h

## 2/3: Creating significant table up_deseq_epi_over_t60h

## 3/3: Creating significant table up_deseq_epi_over_tryp

## Writing excel data according to ebseq for epi_over_t96h: 1/15.

## After (adj)p filter, the up genes table has 4710 genes.

## After (adj)p filter, the down genes table has 6087 genes.

## After fold change filter, the up genes table has 1534 genes.

## After fold change filter, the down genes table has 2442 genes.

## Writing excel data according to ebseq for epi_over_t60h: 2/15.

## After (adj)p filter, the up genes table has 1918 genes.

## After (adj)p filter, the down genes table has 8088 genes.

## After fold change filter, the up genes table has 985 genes.

## After fold change filter, the down genes table has 2494 genes.

## Writing excel data according to ebseq for epi_over_tryp: 3/15.

## After (adj)p filter, the up genes table has 7550 genes.

## After (adj)p filter, the down genes table has 5044 genes.

## After fold change filter, the up genes table has 2854 genes.

## After fold change filter, the down genes table has 3053 genes.

## Printing significant genes to the file: excel/clbr_analyses_sig-v20170810.xlsx

## 1/3: Creating significant table up_ebseq_epi_over_t96h

## 2/3: Creating significant table up_ebseq_epi_over_t60h

## 3/3: Creating significant table up_ebseq_epi_over_tryp

## Writing excel data according to basic for epi_over_t96h: 1/15.

## After (adj)p filter, the up genes table has 6147 genes.

## After (adj)p filter, the down genes table has 6529 genes.

## After fold change filter, the up genes table has 1903 genes.

## After fold change filter, the down genes table has 1999 genes.

## Writing excel data according to basic for epi_over_t60h: 2/15.

## After (adj)p filter, the up genes table has 3660 genes.

## After (adj)p filter, the down genes table has 5316 genes.

## After fold change filter, the up genes table has 1453 genes.

## After fold change filter, the down genes table has 708 genes.

## Writing excel data according to basic for epi_over_tryp: 3/15.

## After (adj)p filter, the up genes table has 2853 genes.

## After (adj)p filter, the down genes table has 2402 genes.

## After fold change filter, the up genes table has 1261 genes.

## After fold change filter, the down genes table has 1315 genes.

## Printing significant genes to the file: excel/clbr_analyses_sig-v20170810.xlsx

## 1/3: Creating significant table up_basic_epi_over_t96h

## 2/3: Creating significant table up_basic_epi_over_t60h

## 3/3: Creating significant table up_basic_epi_over_tryp

## Adding significance bar plots.

4 Epimastigote vs. not Epimastigote

epi_vs <- all_minus_one
epi_vs_meta <- pData(epi_vs)
epi_vs_meta[["epimastigote"]] <- "no"
epi_samples <- epi_vs_meta[["stage"]] == "Epi"
epi_vs_meta[epi_samples, "epimastigote"] <- "yes"
pData(epi_vs[["expressionset"]]) <- epi_vs_meta
epi_vs <- set_expt_conditions(epi_vs, fact="epimastigote")

epi_vs_de <- all_pairwise(epi_vs, model_batch="svaseq", filter=TRUE)

## batch_counts: Before batch/surrogate estimation, 241564 entries are x>1: 100%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

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

## Attempting svaseq estimation with 1 surrogates.

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Using svaseq to visualize before/after batch inclusion.

## Performing a test normalization with: raw

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (22060 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 396 values equal to 0, adding 1 to the matrix.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 237917 entries are x>1: 98%.

## batch_counts: Before batch/surrogate estimation, 396 entries are x==0: 0%.

## batch_counts: Before batch/surrogate estimation, 4347 entries are 0<x<1: 2%.

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

## Attempting svaseq estimation with 2 surrogates.

## There are 7 (0%) elements which are < 0 after batch correction.

## Setting low elements to zero.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

epi_vs_table <- combine_de_tables(epi_vs_de)

## Writing a legend of columns.

## Working on table 1/1: yes_vs_no

epi_vs_sig <- extract_significant_genes(epi_vs_table)

## Writing a legend of columns.

## Writing excel data according to limma for yes_vs_no: 1/5.

## After (adj)p filter, the up genes table has 5295 genes.

## After (adj)p filter, the down genes table has 5661 genes.

## After fold change filter, the up genes table has 1639 genes.

## After fold change filter, the down genes table has 1301 genes.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_limma_yes_vs_no

## Writing excel data according to edger for yes_vs_no: 1/5.

## After (adj)p filter, the up genes table has 4794 genes.

## After (adj)p filter, the down genes table has 5376 genes.

## After fold change filter, the up genes table has 1615 genes.

## After fold change filter, the down genes table has 1962 genes.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_edger_yes_vs_no

## Writing excel data according to deseq for yes_vs_no: 1/5.

## After (adj)p filter, the up genes table has 5754 genes.

## After (adj)p filter, the down genes table has 5871 genes.

## After fold change filter, the up genes table has 1701 genes.

## After fold change filter, the down genes table has 1852 genes.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_deseq_yes_vs_no

## Writing excel data according to ebseq for yes_vs_no: 1/5.

## After (adj)p filter, the up genes table has 4001 genes.

## After (adj)p filter, the down genes table has 4474 genes.

## After fold change filter, the up genes table has 1268 genes.

## After fold change filter, the down genes table has 2272 genes.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_ebseq_yes_vs_no

## Writing excel data according to basic for yes_vs_no: 1/5.

## After (adj)p filter, the up genes table has 4456 genes.

## After (adj)p filter, the down genes table has 4723 genes.

## After fold change filter, the up genes table has 1512 genes.

## After fold change filter, the down genes table has 929 genes.

## Printing significant genes to the file: excel/significant_genes.xlsx

## 1/1: Creating significant table up_basic_yes_vs_no

## Adding significance bar plots.

5 A new request

The following, slightly mysterious request has appeared in my TODO list:

“NEW Stuff with Santuza” 1. DE analyses of CL Brener 1) Epi vs Trypo; 2) Epi vs ama60 2. GO of above

This is peculiar because of how unsatisfactory everyone seemed to find the epimastigote samples previously. I will therefore resurrect the expressionset containing all of the data and seek first to show why this is either a good or bad idea. Maybe it is fine, I do not recall.

The data structure which still includes the epimastigote samples was called ‘all_expt_all’ because it contained all samples for all haplotypes.

all_norm <- sm(normalize_expt(all_expt_all, filter=TRUE, transform="log2", convert="cpm",
                              norm="quant"))
all_plots <- sm(graph_metrics(all_norm))
all_plots$pcaplot

Well, the above PCA plot shows pretty clearly why we might be concerned about the CLBr Tryp samples, one of them looks very much like an A96. The CLBr Epis look ok, I think, but the plot is pretty crowded, I am going to pull only the 3 types of interest and see if it is clearer.

tmp <- subset_expt(all_expt_all, subset="type=='CLBr'&stage!='A96'")
tmp_norm <- sm(normalize_expt(tmp, filter=TRUE, transform="log2", convert="cpm", norm="quant"))
tt <- plot_pca(tmp_norm)
tt$plot

Ok, so I am pretty sure that the only screwball is hpgl0490. I am going to leave it in place for the moment even though it looks funky. I am also a bit torn about what is more appropriate: perform the differential expression analysis with all samples or just the 9 of interest to Santuza for this question. I think it is more appropriate to do all of them.

## This will take a while...
all_norm_all <- sm(normalize_expt(all_expt_all, filter=TRUE))
all_de <- all_pairwise(all_norm_all, model_batch="svaseq")

epi_tryp_60_keepers <- list(
  "clbrepi_vs_tryp" = c("CLBr.Epi", "CLBr.Tryp"),
  "clbrepi_vs_a60" = c("CLBr.Epi", "CLBr.A60"),
  "cl14epi_vs_tryp" = c("CL14.Epi", "CL14.Tryp"),
  "cl14epi_vs_a60" = c("CL14.Epi", "CL14.A60"))
peculiar_data <- combine_de_tables(
  all_de,
  keepers=epi_tryp_60_keepers,
  excel=paste0("excel/epi_tryp_60-v", ver, ".xlsx"))
peculiar_sig <- extract_significant_genes(
  peculiar_data,
  excel=paste0("excel/epi_tryp_60_sig-v", ver, ".xlsx"))

6 Perform goseq with this

all_go <- all_go[, c("GID", "GO")]
colnames(all_go) <- c("ID", "GO")
all_go <- as.data.frame(all_go)
sig <- peculiar_sig$limma$ups$CLBr.Epi_vs_CLBr.Tryp
colnames(all_lengths) <- c("ID", "length")

clbr_up_epi_vs_tryp_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_up_epi_vs_tryp_goseq-v", ver, ".xlsx"))
clbr_up_epi_vs_tryp_goseq$pvalue_plots$bpp_plot_over

clbr_down_epi_vs_tryp_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_down_epi_vs_tryp_goseq-v", ver, ".xlsx"))
clbr_down_epi_vs_tryp_goseq$pvalue_plots$bpp_plot_over

clbr_up_epi_vs_a60_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_up_epi_vs_a60_goseq-v", ver, ".xlsx"))
clbr_up_epi_vs_a60_goseq$pvalue_plots$bpp_plot_over

clbr_down_epi_vs_a60_goseq <- simple_goseq(
  sig_genes=sig,
  go_db=all_go,
  length_db=all_lengths,
  excel=paste0("excel/clbr_down_epi_vs_a60_goseq-v", ver, ".xlsx"))
clbr_down_epi_vs_a60_goseq$pvalue_plots$bpp_plot_over

pander::pander(sessionInfo())

R version 3.6.1 (2019-07-05)

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: parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ruv(v.0.9.7.1), hpgltools(v.1.0), Biobase(v.2.46.0) and BiocGenerics(v.0.32.0)

loaded via a namespace (and not attached): R.utils(v.2.9.2), tidyselect(v.1.0.0), lme4(v.1.1-23), htmlwidgets(v.1.5.1), RSQLite(v.2.2.0), AnnotationDbi(v.1.48.0), grid(v.3.6.1), BiocParallel(v.1.20.1), devtools(v.2.3.0), munsell(v.0.5.0), codetools(v.0.2-16), preprocessCore(v.1.48.0), statmod(v.1.4.34), withr(v.2.2.0), colorspace(v.1.4-1), GOSemSim(v.2.12.1), knitr(v.1.28), rstudioapi(v.0.11), stats4(v.3.6.1), Vennerable(v.3.1.0.9000), robustbase(v.0.93-6), DOSE(v.3.12.0), labeling(v.0.3), urltools(v.1.7.3), GenomeInfoDbData(v.1.2.2), polyclip(v.1.10-0), bit64(v.0.9-7), farver(v.2.0.3), rprojroot(v.1.3-2), vctrs(v.0.2.4), xfun(v.0.13), BiocFileCache(v.1.10.2), R6(v.2.4.1), doParallel(v.1.0.15), GenomeInfoDb(v.1.22.1), graphlayouts(v.0.7.0), locfit(v.1.5-9.4), bitops(v.1.0-6), fgsea(v.1.12.0), gridGraphics(v.0.5-0), DelayedArray(v.0.12.3), assertthat(v.0.2.1), scales(v.1.1.0), nnet(v.7.3-14), ggraph(v.2.0.2), enrichplot(v.1.6.1), gtable(v.0.3.0), sva(v.3.34.0), processx(v.3.4.2), tidygraph(v.1.1.2), rlang(v.0.4.6), genefilter(v.1.68.0), splines(v.3.6.1), rtracklayer(v.1.46.0), acepack(v.1.4.1), checkmate(v.2.0.0), europepmc(v.0.3), BiocManager(v.1.30.10), yaml(v.2.2.1), reshape2(v.1.4.4), GenomicFeatures(v.1.38.2), backports(v.1.1.6), Hmisc(v.4.4-0), qvalue(v.2.18.0), RBGL(v.1.62.1), clusterProfiler(v.3.14.3), tools(v.3.6.1), usethis(v.1.6.1), ggplotify(v.0.0.5), ggplot2(v.3.3.0), ellipsis(v.0.3.0), gplots(v.3.0.3), RColorBrewer(v.1.1-2), sessioninfo(v.1.1.1), ggridges(v.0.5.2), Rcpp(v.1.0.4.6), plyr(v.1.8.6), base64enc(v.0.1-3), progress(v.1.2.2), zlibbioc(v.1.32.0), purrr(v.0.3.4), RCurl(v.1.98-1.2), ps(v.1.3.2), prettyunits(v.1.1.1), rpart(v.4.1-15), openssl(v.1.4.1), viridis(v.0.5.1), cowplot(v.1.0.0), S4Vectors(v.0.24.4), cluster(v.2.1.0), SummarizedExperiment(v.1.16.1), ggrepel(v.0.8.2), colorRamps(v.2.3), fs(v.1.4.1), variancePartition(v.1.16.1), magrittr(v.1.5), data.table(v.1.12.8), DO.db(v.2.9), openxlsx(v.4.1.4), triebeard(v.0.3.0), matrixStats(v.0.56.0), pkgload(v.1.0.2), hms(v.0.5.3), evaluate(v.0.14), xtable(v.1.8-4), pbkrtest(v.0.4-8.6), XML(v.3.99-0.3), jpeg(v.0.1-8.1), IRanges(v.2.20.2), gridExtra(v.2.3), testthat(v.2.3.2), compiler(v.3.6.1), biomaRt(v.2.42.1), tibble(v.3.0.1), KernSmooth(v.2.23-17), crayon(v.1.3.4), R.oo(v.1.23.0), minqa(v.1.2.4), htmltools(v.0.4.0), mgcv(v.1.8-31), corpcor(v.1.6.9), Formula(v.1.2-3), geneplotter(v.1.64.0), tidyr(v.1.0.2), DBI(v.1.1.0), tweenr(v.1.0.1), dbplyr(v.1.4.3), MASS(v.7.3-51.6), rappdirs(v.0.3.1), boot(v.1.3-25), Matrix(v.1.2-18), cli(v.2.0.2), R.methodsS3(v.1.8.0), gdata(v.2.18.0), igraph(v.1.2.5), GenomicRanges(v.1.38.0), pkgconfig(v.2.0.3), rvcheck(v.0.1.8), GenomicAlignments(v.1.22.1), foreign(v.0.8-76), xml2(v.1.3.2), foreach(v.1.5.0), annotate(v.1.64.0), XVector(v.0.26.0), stringr(v.1.4.0), callr(v.3.4.3), digest(v.0.6.25), graph(v.1.64.0), Biostrings(v.2.54.0), rmarkdown(v.2.1), fastmatch(v.1.1-0), htmlTable(v.1.13.3), edgeR(v.3.28.1), curl(v.4.3), Rsamtools(v.2.2.3), gtools(v.3.8.2), nloptr(v.1.2.2.1), lifecycle(v.0.2.0), nlme(v.3.1-147), jsonlite(v.1.6.1), desc(v.1.2.0), viridisLite(v.0.3.0), askpass(v.1.1), limma(v.3.42.2), fansi(v.0.4.1), pillar(v.1.4.3), lattice(v.0.20-41), DEoptimR(v.1.0-8), httr(v.1.4.1), pkgbuild(v.1.0.7), survival(v.3.1-12), GO.db(v.3.10.0), glue(v.1.4.0), remotes(v.2.1.1), zip(v.2.0.4), png(v.0.1-7), iterators(v.1.0.12), pander(v.0.6.3), bit(v.1.1-15.2), ggforce(v.0.3.1), stringi(v.1.4.6), blob(v.1.2.1), DESeq2(v.1.26.0), latticeExtra(v.0.6-29), caTools(v.1.18.0), memoise(v.1.1.0) and dplyr(v.0.8.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 5763684c1dc9aaca47da1768d30388741268c02c

## This is hpgltools commit: Sat Apr 25 14:42:57 2020 -0400: 5763684c1dc9aaca47da1768d30388741268c02c

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

## Saving to 01_annotation-v20170810.rda.xz

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