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"))))
ver <- "202005"

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_one <- subset_expt(all_minus_one, subset="stage!='A96'")

## Using a subset expression.

## There were 11, now there are 8 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_minus_one, filter=TRUE, convert="cpm", norm="quant", transform="log2"))
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 Write count tables

written <- write_expt(all_minus_one,
                      excel=glue::glue("excel/fernanda_written-v{ver}.xlsx"))

## Writing the first sheet, containing a legend and some summary data.

## Writing the raw reads.

## Graphing the raw reads.

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Attempting mixed linear model with: ~  (1|condition)

## Fitting the expressionset to the model, this is slow.

## Dividing work into 100 chunks...

## 
## Total:60 s

## Error in `.rowNamesDF<-`(x, value = value) : invalid 'row.names' length

## A couple of common errors:
## An error like 'vtv downdated' may be because there are too many 0s, filter the data and rerun.
## An error like 'number of levels of each grouping factor must be < number of observations' means
## that the factor used is not appropriate for the analysis - it really only works for factors
## which are shared among multiple samples.

## Retrying with only condition in the model.

## Error in if (grep("object '.*' not found", possibleError$message) == 1) { : 
##   argument is of length zero

## Attempting again with only condition failed.

## Error in simple_varpart(expt, predictor = NULL, factors = c("condition",  : 
## 

## Writing the normalized reads.

## Graphing the normalized reads.

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## Warning in doTryCatch(return(expr), name, parentenv, handler): display list
## redraw incomplete

## varpart sees only 1 batch, adjusting the model accordingly.

## Attempting mixed linear model with: ~  (1|condition)

## Fitting the expressionset to the model, this is slow.

## Dividing work into 100 chunks...

## 
## Total:87 s

## Placing factor: condition at the beginning of the model.

## Writing the median reads by factor.

4 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 1352 low-count genes (21953 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.

keepers <- list(
    "epi_over_t60h" = c("CLBrEpi", "CLBrA60"),
    "epi_over_tryp" = c("CLBrEpi", "CLBrTryp"))

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

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

## batch_counts: Before batch/surrogate estimation, 350 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 (21953 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 350 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, 171845 entries are x>1: 98%.

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

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

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

## Attempting svaseq estimation with 2 surrogates.

## There are 28 (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.

all_minus_tables_sva <- combine_de_tables(
    all_minus_de_sva, keepers=keepers,
    excel=glue::glue("excel/clbr_analyses-sva-v{ver}.xlsx"))

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

## Writing a legend of columns.

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

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

## Found table with CLBrEpi_vs_CLBrA60

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

## Found inverse table with CLBrTryp_vs_CLBrEpi

## Adding venn plots for epi_over_t60h.

## Limma expression coefficients for epi_over_t60h; R^2: 0.884; equation: y = 0.919x + 0.296

## Deseq expression coefficients for epi_over_t60h; R^2: 0.883; equation: y = 0.912x + 0.743

## Edger expression coefficients for epi_over_t60h; R^2: 0.883; equation: y = 0.914x + 0.701

## Adding venn plots for epi_over_tryp.

## Limma expression coefficients for epi_over_tryp; R^2: 0.884; equation: y = 0.919x + 0.296

## Deseq expression coefficients for epi_over_tryp; R^2: 0.883; equation: y = 0.912x + 0.743

## Edger expression coefficients for epi_over_tryp; R^2: 0.883; equation: y = 0.914x + 0.701

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/clbr_analyses-sva-v202005.xlsx.

all_minus_sig_sva <- extract_significant_genes(
    all_minus_tables_sva,
    excel=glue::glue("excel/clbr_analyses_sig-sva-v{ver}.xlsx"))

## Writing a legend of columns.

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

## 1/2: Creating significant table up_limma_epi_over_t60h

## 2/2: Creating significant table up_limma_epi_over_tryp

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

## 1/2: Creating significant table up_edger_epi_over_t60h

## 2/2: Creating significant table up_edger_epi_over_tryp

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

## 1/2: Creating significant table up_deseq_epi_over_t60h

## 2/2: Creating significant table up_deseq_epi_over_tryp

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

## 1/2: Creating significant table up_ebseq_epi_over_t60h

## 2/2: Creating significant table up_ebseq_epi_over_tryp

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

## 1/2: Creating significant table up_basic_epi_over_t60h

## 2/2: Creating significant table up_basic_epi_over_tryp

## Adding significance bar plots.

all_minus_de_nob <- all_pairwise(all_minus_norm, model_batch=FALSE)

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Assuming no batch in model for testing pca.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

all_minus_tables_nob <- combine_de_tables(
    all_minus_de_nob, keepers=keepers,
    excel=glue::glue("excel/clbr_analyses-nob-v{ver}.xlsx"))

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

## Writing a legend of columns.

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

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

## Found table with CLBrEpi_vs_CLBrA60

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

## Found inverse table with CLBrTryp_vs_CLBrEpi

## Adding venn plots for epi_over_t60h.

## Limma expression coefficients for epi_over_t60h; R^2: 0.891; equation: y = 0.923x + 0.275

## Deseq expression coefficients for epi_over_t60h; R^2: 0.891; equation: y = 0.921x + 0.672

## Edger expression coefficients for epi_over_t60h; R^2: 0.891; equation: y = 0.921x + 0.643

## Adding venn plots for epi_over_tryp.

## Limma expression coefficients for epi_over_tryp; R^2: 0.891; equation: y = 0.923x + 0.275

## Deseq expression coefficients for epi_over_tryp; R^2: 0.891; equation: y = 0.921x + 0.672

## Edger expression coefficients for epi_over_tryp; R^2: 0.891; equation: y = 0.921x + 0.643

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/clbr_analyses-nob-v202005.xlsx.

all_minus_sig_nob <- extract_significant_genes(
    all_minus_tables_nob,
    excel=glue::glue("excel/clbr_analyses_nob-v{ver}.xlsx"))

## Writing a legend of columns.

## Printing significant genes to the file: excel/clbr_analyses_nob-v202005.xlsx

## 1/2: Creating significant table up_limma_epi_over_t60h

## 2/2: Creating significant table up_limma_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_nob-v202005.xlsx

## 1/2: Creating significant table up_edger_epi_over_t60h

## 2/2: Creating significant table up_edger_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_nob-v202005.xlsx

## 1/2: Creating significant table up_deseq_epi_over_t60h

## 2/2: Creating significant table up_deseq_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_nob-v202005.xlsx

## 1/2: Creating significant table up_ebseq_epi_over_t60h

## 2/2: Creating significant table up_ebseq_epi_over_tryp

## Printing significant genes to the file: excel/clbr_analyses_nob-v202005.xlsx

## 1/2: Creating significant table up_basic_epi_over_t60h

## 2/2: Creating significant table up_basic_epi_over_tryp

## Adding significance bar plots.

first_sig <- all_minus_sig_sva[["deseq"]][["ups"]][["epi_over_t60h"]]
second_sig <- all_minus_sig_sva[["deseq"]][["ups"]][["epi_over_tryp"]]
both <- c(rownames(first_sig), rownames(second_sig))
length(both)

## [1] 4636

union <- unique(both)
length(union)

## [1] 3551

first_sig <- all_minus_sig_sva[["deseq"]][["downs"]][["epi_over_t60h"]]
second_sig <- all_minus_sig_sva[["deseq"]][["downs"]][["epi_over_tryp"]]
both <- c(rownames(first_sig), rownames(second_sig))
length(both)

## [1] 4179

union <- unique(both)
length(union)

## [1] 3919

first_sig <- all_minus_sig_nob[["deseq"]][["ups"]][["epi_over_t60h"]]
second_sig <- all_minus_sig_nob[["deseq"]][["ups"]][["epi_over_tryp"]]
both <- c(rownames(first_sig), rownames(second_sig))
length(both)

## [1] 4790

union <- unique(both)
length(union)

## [1] 3690

first_sig <- all_minus_sig_nob[["deseq"]][["downs"]][["epi_over_t60h"]]
second_sig <- all_minus_sig_nob[["deseq"]][["downs"]][["epi_over_tryp"]]
both <- c(rownames(first_sig), rownames(second_sig))
length(both)

## [1] 4158

union <- unique(both)
length(union)

## [1] 3897

5 Epimastigote vs. not Epimastigote

epi_vs <- all_minus_one

## Hey doofus, you know that set_expt_conditions can do
## These next 5 lines without having to do this
## I mean damn, you wrote it one would expect you to know...
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_sva <- all_pairwise(epi_vs, model_batch="svaseq", filter=TRUE)

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

## batch_counts: Before batch/surrogate estimation, 350 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 (21953 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 350 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, 171845 entries are x>1: 98%.

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

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

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

## Attempting svaseq estimation with 1 surrogates.

## There are 5 (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_sva <- combine_de_tables(epi_vs_de_sva, excel="excel/epi_vs_all_sva-de.xlsx")

## Deleting the file excel/epi_vs_all_sva-de.xlsx before writing the tables.

## Writing a legend of columns.

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

## Working on table 1/1: yes_vs_no

## Adding venn plots for yes_vs_no.

## Limma expression coefficients for yes_vs_no; R^2: 0.876; equation: y = 0.972x + 0.102

## Deseq expression coefficients for yes_vs_no; R^2: 0.893; equation: y = 1.07x - 0.587

## Edger expression coefficients for yes_vs_no; R^2: 0.893; equation: y = 1.07x - 0.597

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/epi_vs_all_sva-de.xlsx.

epi_vs_sig_sva <- extract_significant_genes(epi_vs_table_sva, excel="excel/epi_vs_all_sva-sig.xlsx")

## Writing a legend of columns.

## Printing significant genes to the file: excel/epi_vs_all_sva-sig.xlsx

## 1/1: Creating significant table up_limma_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_sva-sig.xlsx

## 1/1: Creating significant table up_edger_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_sva-sig.xlsx

## 1/1: Creating significant table up_deseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_sva-sig.xlsx

## 1/1: Creating significant table up_ebseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_sva-sig.xlsx

## 1/1: Creating significant table up_basic_yes_vs_no

## Adding significance bar plots.

epi_vs_sig_sva$summary_df

##           limma_change_counts_up limma_change_counts_down
## yes_vs_no                   1749                     1516
##           edger_change_counts_up edger_change_counts_down
## yes_vs_no                   1733                     1658
##           deseq_change_counts_up deseq_change_counts_down
## yes_vs_no                   1764                     1625
##           ebseq_change_counts_up ebseq_change_counts_down
## yes_vs_no                   1125                     2019
##           basic_change_counts_up basic_change_counts_down
## yes_vs_no                   1225                      294

epi_vs_de_batch <- all_pairwise(epi_vs, model_batch=TRUE, filter=TRUE)

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Using limma's removeBatchEffect to visualize with(out) batch inclusion.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

epi_vs_table_batch <- combine_de_tables(epi_vs_de_batch, excel="excel/epi_vs_all_batch-de.xlsx")

## Deleting the file excel/epi_vs_all_batch-de.xlsx before writing the tables.

## Writing a legend of columns.

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

## Working on table 1/1: yes_vs_no

## Adding venn plots for yes_vs_no.

## Limma expression coefficients for yes_vs_no; R^2: 0.819; equation: y = 0.899x + 0.491

## Deseq expression coefficients for yes_vs_no; R^2: 0.865; equation: y = 1.08x - 0.689

## Edger expression coefficients for yes_vs_no; R^2: 0.864; equation: y = 1.08x - 0.677

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/epi_vs_all_batch-de.xlsx.

epi_vs_sig_batch <- extract_significant_genes(epi_vs_table_batch, excel="excel/epi_vs_all_batch-sig.xlsx")

## Writing a legend of columns.

## Printing significant genes to the file: excel/epi_vs_all_batch-sig.xlsx

## 1/1: Creating significant table up_limma_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_batch-sig.xlsx

## 1/1: Creating significant table up_edger_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_batch-sig.xlsx

## 1/1: Creating significant table up_deseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_batch-sig.xlsx

## 1/1: Creating significant table up_ebseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_batch-sig.xlsx

## 1/1: Creating significant table up_basic_yes_vs_no

## Adding significance bar plots.

epi_vs_sig_batch$summary_df

##           limma_change_counts_up limma_change_counts_down
## yes_vs_no                   2001                      901
##           edger_change_counts_up edger_change_counts_down
## yes_vs_no                   1397                     1154
##           deseq_change_counts_up deseq_change_counts_down
## yes_vs_no                   1496                     1096
##           ebseq_change_counts_up ebseq_change_counts_down
## yes_vs_no                   1125                     2019
##           basic_change_counts_up basic_change_counts_down
## yes_vs_no                   1225                      294

epi_vs_de_nobatch <- all_pairwise(epi_vs, model_batch=FALSE, filter=TRUE)

## Plotting a PCA before surrogates/batch inclusion.

## Not putting labels on the plot.

## Assuming no batch in model for testing pca.

## Not putting labels on the plot.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

epi_vs_table_nobatch <- combine_de_tables(epi_vs_de_nobatch, excel="excel/epi_vs_all_nobatch-de.xlsx")

## Deleting the file excel/epi_vs_all_nobatch-de.xlsx before writing the tables.

## Writing a legend of columns.

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

## Working on table 1/1: yes_vs_no

## Adding venn plots for yes_vs_no.

## Limma expression coefficients for yes_vs_no; R^2: 0.819; equation: y = 0.899x + 0.491

## Deseq expression coefficients for yes_vs_no; R^2: 0.865; equation: y = 1.08x - 0.689

## Edger expression coefficients for yes_vs_no; R^2: 0.864; equation: y = 1.08x - 0.677

## Writing summary information, compare_plot is: TRUE.

## Performing save of excel/epi_vs_all_nobatch-de.xlsx.

epi_vs_sig_nobatch <- extract_significant_genes(epi_vs_table_nobatch, excel="excel/epi_vs_all_nobatch-sig.xlsx")

## Writing a legend of columns.

## Printing significant genes to the file: excel/epi_vs_all_nobatch-sig.xlsx

## 1/1: Creating significant table up_limma_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_nobatch-sig.xlsx

## 1/1: Creating significant table up_edger_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_nobatch-sig.xlsx

## 1/1: Creating significant table up_deseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_nobatch-sig.xlsx

## 1/1: Creating significant table up_ebseq_yes_vs_no

## Printing significant genes to the file: excel/epi_vs_all_nobatch-sig.xlsx

## 1/1: Creating significant table up_basic_yes_vs_no

## Adding significance bar plots.

epi_vs_sig_nobatch$summary_df

##           limma_change_counts_up limma_change_counts_down
## yes_vs_no                   2001                      901
##           edger_change_counts_up edger_change_counts_down
## yes_vs_no                   1397                     1154
##           deseq_change_counts_up deseq_change_counts_down
## yes_vs_no                   1496                     1096
##           ebseq_change_counts_up ebseq_change_counts_down
## yes_vs_no                   1125                     2019
##           basic_change_counts_up basic_change_counts_down
## yes_vs_no                   1225                      294

6 Ontology searching

For these contrasts, we are doing ontology searches of a couple contrasts:

1. epimastigote vs. trypsomastigote
2. 60 hours vs. epimastigote

But first we need to get the correct ontology data frame and gene lengths. There is one aspect of this worth considering, in response to a query on the bioconductor email lists, I made a change to how I collect ontology information for the species of the eupathdb. Previously it was only the GOSLIM table, but now I include both that and the GO.db table. I do not know if this will affect these results.

Also, I am going to assume the default metric of ‘significant’ from the perspective of deseq2.

I might try something new, too, I also now have pretty good tables from reactome for the tryps.

all_go <- all_go[, c("GID", "GO")]
colnames(all_go) <- c("ID", "GO")
all_go <- as.data.frame(all_go)

colnames(all_lengths) <- c("ID", "length")

6.1 Epimastigote vs. Trypomastigote

sig_up <- all_minus_sig[["deseq"]][["ups"]][["epi_over_tryp"]]
sig_down <- all_minus_sig[["deseq"]][["downs"]][["epi_over_tryp"]]
epi_vs_tryp_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_tryp_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 1132 genes out of 2404 from the sig_genes in the go_db.

## Found 2404 genes out of 2404 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_tryp_goseq-v202005.xlsx.

epi_vs_tryp_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_tryp_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 385 genes out of 3525 from the sig_genes in the go_db.

## Found 3525 genes out of 3525 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_tryp_goseq-v202005.xlsx.

epi_vs_tryp_up_topgo <- simple_topgo(
    sig_genes=sig_up,
    go_db=all_go,
    excel=glue::glue("excel/clbr_up_epi_vs_tryp_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_up_epi_vs_tryp_topgo-v202005.xlsx.

epi_vs_tryp_down_topgo <- simple_topgo(
    sig_genes=sig_down,
    go_db=all_go,
    excel=glue::glue("excel/clbr_down_epi_vs_tryp_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_down_epi_vs_tryp_topgo-v202005.xlsx.

epi_vs_tryp_up_gostats <- simple_gostats(
    sig_genes=sig_up,
    gff="~/scratch/libraries/genome/tcruzi_all.gff",
    gff_type="gene", gff_id="ID",
    go_db=all_go,
    excel=glue::glue("excel/clbr_up_epi_vs_tryp_gostats-v{ver}.xlsx"))

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=TRUE)

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Returning a df with 24 columns and 25100 rows.

## simple_gostats(): gff_type is: gene. Change that if there are bad merges.

## simple_gostats(): type cds has 0 annotations.

## simple_gostats(): type gene has 25100 annotations.

## simple_gostats(): type exon has 0 annotations.

## simple_gostats(): type protein_coding has 0 annotations.

## simple_gostats(): the current annotations has: 25100 rows and 24 columns.

## simple_gostats(): Creating the gene set collection.  This is slow.

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:Rgraphviz':
## 
##     from, to

## The following object is masked from 'package:base':
## 
##     expand.grid

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:topGO':
## 
##     members

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:GOstats':
## 
##     makeGOGraph

## simple_gostats(): Performing MF GSEA.

## Found 229 over MF categories.

## simple_gostats(): Performing BP GSEA.

## Found 262 over BP categories.

## simple_gostats(): Performing CC GSEA.

## Found 80 over CC categories.

## simple_gostats(): Performing under MF GSEA.

## Found 75 under MF categories.

## simple_gostats(): Performing under BP GSEA.

## Found 186 under BP categories.

## simple_gostats(): Performing under CC GSEA.

## Found 30 under CC categories.

## Writing data to: excel/clbr_up_epi_vs_tryp_gostats-v202005.xlsx.

## Error in topGO::showSigOfNodes(mf_GOdata, mf_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in topGO::showSigOfNodes(bp_GOdata, bp_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in topGO::showSigOfNodes(cc_GOdata, cc_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in `[[<-.data.frame`(`*tmp*`, "ontology", value = "BP") : 
##   replacement has 1 row, data has 0

epi_vs_tryp_down_gostats <- simple_gostats(
    sig_genes=sig_down,
    gff="~/scratch/libraries/genome/tcruzi_all.gff",
    gff_type="gene", gff_id="ID",
    go_db=all_go,
    excel=glue::glue("excel/clbr_down_epi_vs_tryp_gostats-v{ver}.xlsx"))

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=TRUE)

## Trying attempt: rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Had a successful gff import with rtracklayer::import.gff3(gff, sequenceRegionsAsSeqinfo=FALSE)

## Returning a df with 24 columns and 25100 rows.

## simple_gostats(): gff_type is: gene. Change that if there are bad merges.

## simple_gostats(): type cds has 0 annotations.

## simple_gostats(): type gene has 25100 annotations.

## simple_gostats(): type exon has 0 annotations.

## simple_gostats(): type protein_coding has 0 annotations.

## simple_gostats(): the current annotations has: 25100 rows and 24 columns.

## simple_gostats(): Creating the gene set collection.  This is slow.

## simple_gostats(): Performing MF GSEA.

## Found 66 over MF categories.

## simple_gostats(): Performing BP GSEA.

## Found 94 over BP categories.

## simple_gostats(): Performing CC GSEA.

## Found 11 over CC categories.

## simple_gostats(): Performing under MF GSEA.

## Found 51 under MF categories.

## simple_gostats(): Performing under BP GSEA.

## Found 88 under BP categories.

## simple_gostats(): Performing under CC GSEA.

## Found 49 under CC categories.

## Writing data to: excel/clbr_down_epi_vs_tryp_gostats-v202005.xlsx.

## Error in topGO::showSigOfNodes(mf_GOdata, mf_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in topGO::showSigOfNodes(bp_GOdata, bp_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in topGO::showSigOfNodes(cc_GOdata, cc_over_nodes, useInfo = "all",  : 
##   No nodes were selected
## Error in `[[<-.data.frame`(`*tmp*`, "ontology", value = "BP") : 
##   replacement has 1 row, data has 0

6.1.1 Plots

## up
epi_vs_tryp_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_tryp_up_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_up_topgo[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_tryp_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_tryp_down_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_tryp_down_topgo[["pvalue_plots"]][["mfp_plot_over"]]

6.2 Epimastigote vs. 60 hours

sig_up <- all_minus_sig[["deseq"]][["ups"]][["epi_over_t60h"]]
sig_down <- all_minus_sig[["deseq"]][["downs"]][["epi_over_t60h"]]
epi_vs_t60_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_t60_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 854 genes out of 2489 from the sig_genes in the go_db.

## Found 2489 genes out of 2489 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_t60_goseq-v202005.xlsx.

epi_vs_t60_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_t60_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 210 genes out of 624 from the sig_genes in the go_db.

## Found 624 genes out of 624 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_t60_goseq-v202005.xlsx.

epi_vs_t60_up_topgo <- simple_topgo(
    sig_genes=sig_up,
    go_db=all_go,
    excel=glue::glue("excel/clbr_up_epi_vs_t60_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_up_epi_vs_t60_topgo-v202005.xlsx.

epi_vs_t60_down_topgo <- simple_topgo(
    sig_genes=sig_down,
    go_db=all_go,
    excel=glue::glue("excel/clbr_down_epi_vs_t60_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_down_epi_vs_t60_topgo-v202005.xlsx.

6.2.1 Plots

## up
epi_vs_t60_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_t60_up_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_up_topgo[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_t60_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_t60_down_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_t60_down_topgo[["pvalue_plots"]][["mfp_plot_over"]]

7 Epi vs everything else

sig_up <- epi_vs_sig[["deseq"]][["ups"]][["yes_vs_no"]]
sig_down <- epi_vs_sig[["deseq"]][["downs"]][["yes_vs_no"]]
epi_vs_all_up_goseq <- simple_goseq(
    sig_genes=sig_up,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_up_epi_vs_all_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 744 genes out of 1701 from the sig_genes in the go_db.

## Found 1701 genes out of 1701 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_up_epi_vs_all_goseq-v202005.xlsx.

epi_vs_all_down_goseq <- simple_goseq(
    sig_genes=sig_down,
    go_db=all_go,
    length_db=all_lengths,
    excel=glue::glue("excel/clbr_down_epi_vs_all_goseq-v{ver}.xlsx"))

## Using the row names of your table.

## Found 222 genes out of 1852 from the sig_genes in the go_db.

## Found 1852 genes out of 1852 from the sig_genes in the length_db.

## Using manually entered categories.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

## simple_goseq(): Calculating q-values

## simple_goseq(): Filling godata with terms, this is slow.

## Testing that go categories are defined.

## Removing undefined categories.

## Gathering synonyms.

## Gathering category definitions.

## simple_goseq(): Making pvalue plots for the ontologies.

## Writing data to: excel/clbr_down_epi_vs_all_goseq-v202005.xlsx.

epi_vs_all_up_topgo <- simple_topgo(
    sig_genes=sig_up,
    go_db=all_go,
    excel=glue::glue("excel/clbr_up_epi_vs_all_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_up_epi_vs_all_topgo-v202005.xlsx.

epi_vs_all_down_topgo <- simple_topgo(
    sig_genes=sig_down,
    go_db=all_go,
    excel=glue::glue("excel/clbr_down_epi_vs_all_topgo-v{ver}.xlsx"))

## simple_topgo(): Set densities=TRUE for ontology density plots.

## Writing data to: excel/clbr_down_epi_vs_all_topgo-v202005.xlsx.

7.0.1 Plots

## up
epi_vs_all_up_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_up_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_all_up_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_up_topgo[["pvalue_plots"]][["mfp_plot_over"]]

## down
epi_vs_all_down_goseq[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_down_goseq[["pvalue_plots"]][["mfp_plot_over"]]

epi_vs_all_down_topgo[["pvalue_plots"]][["bpp_plot_over"]]

epi_vs_all_down_topgo[["pvalue_plots"]][["mfp_plot_over"]]

8 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"))

9 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 4.0.0 (2020-04-24)

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

other attached packages: GOstats(v.2.54.0), GO.db(v.3.11.4), AnnotationDbi(v.1.50.0), IRanges(v.2.22.2), S4Vectors(v.0.26.1), Rgraphviz(v.2.32.0), graph(v.1.66.0), SparseM(v.1.78), topGO(v.2.40.0), ruv(v.0.9.7.1), BiocParallel(v.1.22.0), variancePartition(v.1.18.0), hpgltools(v.1.0), Biobase(v.2.48.0) and BiocGenerics(v.0.34.0)

loaded via a namespace (and not attached): R.utils(v.2.9.2), tidyselect(v.1.1.0), lme4(v.1.1-23), RSQLite(v.2.2.0), Rtsne(v.0.15), devtools(v.2.3.0), scatterpie(v.0.1.4), munsell(v.0.5.0), preprocessCore(v.1.50.0), codetools(v.0.2-16), statmod(v.1.4.34), withr(v.2.2.0), colorspace(v.1.4-1), GOSemSim(v.2.14.0), Category(v.2.54.0), knitr(v.1.28), rstudioapi(v.0.11), Vennerable(v.3.1.0.9000), robustbase(v.0.93-6), DOSE(v.3.14.0), labeling(v.0.3), urltools(v.1.7.3), GenomeInfoDbData(v.1.2.3), polyclip(v.1.10-0), bit64(v.0.9-7), farver(v.2.0.3), rprojroot(v.1.3-2), downloader(v.0.4), vctrs(v.0.3.1), generics(v.0.0.2), xfun(v.0.14), BiocFileCache(v.1.12.0), fastcluster(v.1.1.25), R6(v.2.4.1), doParallel(v.1.0.15), GenomeInfoDb(v.1.24.0), graphlayouts(v.0.7.0), locfit(v.1.5-9.4), bitops(v.1.0-6), fgsea(v.1.14.0), gridGraphics(v.0.5-0), DelayedArray(v.0.14.0), assertthat(v.0.2.1), scales(v.1.1.1), ggraph(v.2.0.3), enrichplot(v.1.8.1), gtable(v.0.3.0), sva(v.3.36.0), processx(v.3.4.2), tidygraph(v.1.2.0), rlang(v.0.4.6), genefilter(v.1.70.0), splines(v.4.0.0), rtracklayer(v.1.48.0), europepmc(v.0.4), BiocManager(v.1.30.10), yaml(v.2.2.1), reshape2(v.1.4.4), GenomicFeatures(v.1.40.0), backports(v.1.1.7), qvalue(v.2.20.0), RBGL(v.1.64.0), clusterProfiler(v.3.16.0), tools(v.4.0.0), usethis(v.1.6.1), ggplotify(v.0.0.5), ggplot2(v.3.3.1), ellipsis(v.0.3.1), 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.34.0), BiasedUrn(v.1.07), purrr(v.0.3.4), RCurl(v.1.98-1.2), ps(v.1.3.3), prettyunits(v.1.1.1), openssl(v.1.4.1), viridis(v.0.5.1), cowplot(v.1.0.0), SummarizedExperiment(v.1.18.1), ggrepel(v.0.8.2), colorRamps(v.2.3), fs(v.1.4.1), magrittr(v.1.5), data.table(v.1.12.8), openxlsx(v.4.1.5), DO.db(v.2.9), triebeard(v.0.3.0), goseq(v.1.40.0), matrixStats(v.0.56.0), pkgload(v.1.1.0), 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), gridExtra(v.2.3), testthat(v.2.3.2), compiler(v.4.0.0), biomaRt(v.2.44.0), tibble(v.3.0.1), KernSmooth(v.2.23-16), crayon(v.1.3.4), R.oo(v.1.23.0), minqa(v.1.2.4), htmltools(v.0.4.0), corpcor(v.1.6.9), mgcv(v.1.8-31), tidyr(v.1.1.0), geneplotter(v.1.66.0), DBI(v.1.1.0), tweenr(v.1.0.1), geneLenDataBase(v.1.24.0), dbplyr(v.1.4.4), MASS(v.7.3-51.5), rappdirs(v.0.3.1), boot(v.1.3-24), Matrix(v.1.2-18), cli(v.2.0.2), R.methodsS3(v.1.8.0), quadprog(v.1.5-8), gdata(v.2.18.0), igraph(v.1.2.5), GenomicRanges(v.1.40.0), pkgconfig(v.2.0.3), rvcheck(v.0.1.8), GenomicAlignments(v.1.24.0), xml2(v.1.3.2), foreach(v.1.5.0), annotate(v.1.66.0), XVector(v.0.28.0), AnnotationForge(v.1.30.1), stringr(v.1.4.0), callr(v.3.4.3), digest(v.0.6.25), Biostrings(v.2.56.0), rmarkdown(v.2.2), fastmatch(v.1.1-0), edgeR(v.3.30.3), GSEABase(v.1.50.1), directlabels(v.2020.1.31), curl(v.4.3), Rsamtools(v.2.4.0), 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.44.3), fansi(v.0.4.1), pillar(v.1.4.4), lattice(v.0.20-41), httr(v.1.4.1), DEoptimR(v.1.0-8), pkgbuild(v.1.0.8), survival(v.3.1-12), glue(v.1.4.1), remotes(v.2.1.1), zip(v.2.0.4), 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.28.1), caTools(v.1.18.0), memoise(v.1.1.0) and dplyr(v.1.0.0)

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 71e9e0510e6a865fc024d8cf57090f2f02602a35

## This is hpgltools commit: Fri Jun 26 13:06:16 2020 -0400: 71e9e0510e6a865fc024d8cf57090f2f02602a35

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

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

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