Histone Deacetylase 5 under/over expression in R.norvegicus.

1 TODO:

Check what was done with the over expression vector sequence, did I definitely map against it? (I think so).

2 Introduction

I still need to do some reading to understand this series of experiments; but I believe the idea was to use a short-hairpin to lower expression of hdac5 in some cases, and an adenovirus overexpression vector in others.

3 Preprocessing

I passed these samples to my default preprocessing pipeline on the umiacs cluster with the following command:

cd preprocessing
start=$(pwd)
for i in $(/bin/ls -d *); do
    cd $i
    cyoa --task pipe --method prnaseq --species rnorvegicus_7.2_107 --gff_type gene \
         --stranded reverse --gff_tag ID \
         --input $(/bin/ls *.fastq.gz | tr '\n' ':' | sed 's/:$//g')
    cd $start
done

The above command iterates over every directory of raw data I created and runs my rnaseq (–method prnaseq) pipeline (–task pipe) using a new copy of the ensembl Rattus norvegicus genome (version 7.2) from ensembl’s release 107 (–species rnorvegicus_7.2_107). When I downloaded it, I saved a gene feature file and the fasta genome with the same prefix and therefore htseq will use the gene annotations (–gff_type gene) by identifier ‘ID’ (–gff_tag ID).

I did not provide any instructions for trimming nor mapping, so it will use my defaults. All of the scripts run on the cluster are therefore located in the preprocessing/*/scripts/ directories in files named lexically in the order they were performed.

4 Genome Annotation

I will extract annotations from ensembl, since it is my favorite. My load_biomart_annotation() function attempts to find an archive server from one year ago, this may be made explicit if one wishes.

One important caveat, when creating the expressionset data structure, the gene IDs must match and the gff usually uses IDs which are a specific combination of columns from the ensembl table. It usually takes a moment to figure out how the two data sources line up; as a result sometimes it is easier to just extract annotations directly from the gff file.

columns <- c("ensembl_gene_id", "version", "ensembl_transcript_id", "transcript_version",
             "description", "gene_biotype", "entrezgene_description", "rgd_symbol", "band")
rn_annot <- load_biomart_annotations(species = "rnorvegicus", gene_requests = columns,
                                     year = "2020", month = 8)

## The biomart annotations file already exists, loading from it.

5 Metadata collection

I often just copy/paste the relevant fields from the output files created by trimomatic/htseq/hisat/etc. However, I did make a fun function: gather_preprocessing_metadata() which in theory at least can fill out a metadata xlsx file.

written <- gather_preprocessing_metadata("sample_sheets/all_samples.xlsx")

I noticed a couple problems with my function:

1. It does not like dates in the MM/DD/YYYY format. There are a couple reasons for which this may have happened.
2. It did not save the condition/batch columns I created, so I copy/pasted them from the previous sheet.

6 Initial Data Structure

Let us create an expression set of the metadata and annotations. Note that the

annot <- rn_annot[["annotation"]]
rownames(annot) <- make.names(paste0("gene:", annot[["ensembl_gene_id"]]), unique = TRUE)
rownames(annot) <- gsub(pattern="\\.", replacement = ":", x = rownames(annot))

rn_expt <- create_expt(metadata="sample_sheets/all_samples_modified.xlsx",
                       gene_info=annot, file_column="hisatcounttable")

## Reading the sample metadata.

## The sample definitions comprises: 24 rows(samples) and 30 columns(metadata fields).

## Matched 18641 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 23139 features and 24 samples.

written <- write_expt(rn_expt, excel="excel/rn_expt.xlsx")

## Deleting the file excel/rn_expt.xlsx before writing the tables.

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

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

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff,  : 
##   The variables specified in this model are redundant,
## so the design matrix is not full rank

## Retrying with only condition in the model.
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following object is masked from 'package:S4Vectors':
## 
##     expand
## 
## 
## Total:82 s
## `geom_smooth()` using formula = 'y ~ x'

## Error in checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff,  : 
##   The variables specified in this model are redundant,
## so the design matrix is not full rank

## Retrying with only condition in the model.
## 
## Total:65 s

7 Initial plots

Let us look briefly at a couple of likely plots of the data.

plot_libsize(rn_expt)$plot

plot_nonzero(rn_expt)$plot

## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

plot_variance_coefficients(rn_expt)$plot

## Naively calculating coefficient of variation/dispersion with respect to condition.
## Finished calculating dispersion estimates.

This is worrisome. The trimming kept ~95% of the reads, which is awesome. Mapping appears to have been in the 90%s. I would therefore have assumed the numbers of reads mapped to genes to be higher than 1.1<x<3.7 million.

My first place to look into this is the fastqc output. My initial hypothesis is that the sequence duplication file will show a huge number of reads mapping to rRNA or some other ncRNA. It appears that guess was incorrect, there do not appear to be significantly over represented sequences in this data. Maybe lets look at the reads in IGV?

rn_norm <- normalize_expt(rn_expt, filter=TRUE, convert="cpm",
                          norm="quant", transform="log2")

## Removing 8998 low-count genes (14141 remaining).

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

plot_pca(rn_norm)$plot

## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

rn_nbsva <- normalize_expt(rn_expt, filter=TRUE, convert="cpm",
                           batch="svaseq", transform="log2")

## Removing 8998 low-count genes (14141 remaining).

## Setting 920 low elements to zero.

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

plot_pca(rn_nbsva)$plot

## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

7.1 Test other factors

The sample metadata is a little limited vis a vis the ways in which the samples are categorized. We have, in effect, 4 primary factors:

• The expression method used: either (I am assuming this from the sample name) a short hairpin which was presumably transfected, or an adenovirus overexpression transfection.
• The perturbation goal: There are three goals used in these transfections, either do nothing, over-express HDAC5, or lower its expression.
• The target of the transfection: This is either luc (luciferase I assume), GFP, or histone deacetylase 5 (which of course has its own tremendous set of changes it causes to the transcriptional landscape).
• The combination of the above: There are at least two which one may express this, I chose 4 categories: control_kd, hdac5_kd, control_over, and hdac5_over. One may instead choose to combine the controls into one category, though I would argue (and the data agrees with me) that doing so is wrong; but in the ideal world the perturbation of the luciferase/GFP targets would do nothing and these categories would be identical.

With the above in mind, let us consider some of the other categorizations of the data and see if any make sense. I know already from the first picture that knockdown and overexpression are the two primary sources of variance.

rn_var <- simple_varpart(
    rn_expt, chosen_factor = "insertion",
    factors = c("replicatenumber", "insertion", "perturbation", "worked"))

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable

## 
## Total:149 s

pp(file = "images/rn_varpart.png", image = rn_var$partition_plot)

## Warning in pp(file = "images/rn_varpart.png", image = rn_var$partition_plot):
## There is no device to shut down.

kd_over <- set_expt_conditions(rn_expt, fact = "insertion")

## 
##   adeno_assoc short_hairpin 
##            14            10

## This will, by definition, just be the exact plot as before colored differently.
kd_over_norm <- normalize_expt(kd_over, transform = "log2", convert = "cpm",
                               filter = TRUE, norm = "quant")

## Removing 8998 low-count genes (14141 remaining).

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

pp(file = "images/over_pca.png", image = plot_pca(kd_over_norm)$plot)

## Warning in pp(file = "images/over_pca.png", image =
## plot_pca(kd_over_norm)$plot): There is no device to shut down.

perturb <- set_expt_conditions(rn_expt, fact = "perturbation")

## 
##     control   knockdown overexpress 
##          12           5           7

## This will, by definition, just be the exact plot as before colored differently.
perturb_norm <- normalize_expt(perturb, transform = "log2", convert = "cpm",
                               filter = TRUE, norm = "quant")

## Removing 8998 low-count genes (14141 remaining).
## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.

plot_pca(perturb_norm)$plot

## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

perturb_nb <- normalize_expt(perturb, transform = "log2", convert = "cpm",
                             filter = TRUE, batch = "svaseq")

## Removing 8998 low-count genes (14141 remaining).

## Setting 662 low elements to zero.

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

plot_pca(perturb_nb)$plot

target <- set_expt_conditions(rn_expt, fact = "target")

## 
##        GFP      HDAC5 luciferase 
##          7         12          5

## This will, by definition, just be the exact plot as before colored differently.
target_norm <- normalize_expt(target, transform = "log2", convert = "cpm",
                              filter = TRUE, norm = "quant")

## Removing 8998 low-count genes (14141 remaining).

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

plot_pca(target_norm)$plot
target_nb <- normalize_expt(target, transform = "log2", convert = "cpm",
                            filter = TRUE, batch = "svaseq")

## Removing 8998 low-count genes (14141 remaining).

## Setting 560 low elements to zero.

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

plot_pca(target_nb)$plot

8 Separate sets

NOTE!: We made a column in the sample sheet called ‘worked’, it is a categorical defined by its neighboring column ‘eYFP ratio’, which in turn is defined as the percentage of the unmapped to mouse reads which successfully mapped to the eYFP vector. That ratio ranged from 0.0 to 0.03 (e.g. absolute 0 to 3% of the remaining reads). I therefore made categories: none, hundredths, tenths, ones.

We then made the decision that a successful sample should have more than 0.09% reads mapped to YFP and excluded the others.

Simultaneously, experimental batch c1e3 is off on its own with only 2 samples and so we decided to exclude it. (There were two dates on which samples were collected and 3 dates of RNA isolation; c1e3 is the first sample date, but RNA on a later date and comprised 1 control and 1 over expression sample).

kd <- subset_expt(rn_expt, subset = "insertion=='short_hairpin'") %>%
  subset_expt(subset = "worked!='hundredths'")

## subset_expt(): There were 24, now there are 10 samples.

## subset_expt(): There were 10, now there are 5 samples.

kd_norm <- normalize_expt(kd, transform = "log2", convert = "cpm",
                          norm = "quant", filter = TRUE)

## Removing 9775 low-count genes (13364 remaining).

pp(file = "images/knockdown_pca.png", image = plot_pca(kd_norm)$plot)

## Warning in pp(file = "images/knockdown_pca.png", image =
## plot_pca(kd_norm)$plot): There is no device to shut down.

kd_nb <- normalize_expt(kd, transform = "log2",
                        convert = "cpm", batch = "svaseq", filter = TRUE)

## Removing 9775 low-count genes (13364 remaining).

## Setting 3 low elements to zero.

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

plot_pca(kd_nb)$plot

oe <- subset_expt(rn_expt, subset = "insertion=='adeno_assoc'") %>%
  subset_expt(subset = "batch!='c1e3'")

## subset_expt(): There were 24, now there are 14 samples.

## subset_expt(): There were 14, now there are 12 samples.

oe_norm <- normalize_expt(oe, transform = "log2",
                          convert = "cpm", norm = "quant", filter = TRUE)

## Removing 10067 low-count genes (13072 remaining).

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

plot_pca(oe_norm)$plot
oe_nb <- normalize_expt(oe, transform = "log2",
                        convert = "cpm", batch = "svaseq", filter = TRUE)

## Removing 10067 low-count genes (13072 remaining).

## Setting 40 low elements to zero.

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

plot_pca(oe_nb)$plot

9 Attempt a perturbation differential expression.

kd_de <- all_pairwise(kd, model_batch = FALSE, filter = TRUE)

## This DE analysis will perform all pairwise comparisons among:

## 
## control_kd   hdac5_kd 
##          3          2

## This will pre-filter the input data using normalize_expt's: TRUE argument.

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

kd_keeper <- list("down_control" = c("hdac5kd", "controlkd"))
kd_table <- combine_de_tables(kd_de, keepers = kd_keeper, excel = "excel/kd_table.xlsx")

## Deleting the file excel/kd_table.xlsx before writing the tables.

## Getting the genes >= 1.5 stdevs away from the mean of all.
## Getting the genes >= 1.5 stdevs away from the mean of all.
## Getting the genes >= 1.5 stdevs away from the mean of all.

kd_sig <- extract_significant_genes(kd_table, according_to = "deseq", p = 0.1,
                                    lfc = 0.38, excel = "excel/kd_sig.xlsx")

## Deleting the file excel/kd_sig.xlsx before writing the tables.

kd_up <- kd_sig[["deseq"]][["ups"]][["down_control"]]
kd_down <- kd_sig[["deseq"]][["downs"]][["down_control"]]
kd_gp <- simple_gprofiler(kd_down, species = "rnorvegicus")

## Warning in simple_gprofiler2(...): Hey, it looks like you forgot to strip off
## the htseq prefix for the gene IDs.

## No results to show
## Please make sure that the organism is correct or set significant = FALSE

## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE

## Error in gprofiler2::gost(query = gene_ids, organism = species, evcodes = evcodes,  : 
##   Bad request, response code 400

## No results to show
## Please make sure that the organism is correct or set significant = FALSE
## No results to show
## Please make sure that the organism is correct or set significant = FALSE

oe_de <- all_pairwise(oe, model_batch = "svaseq", filter = TRUE)

## This DE analysis will perform all pairwise comparisons among:

## 
## control_over   hdac5_over 
##            6            6

## This analysis will include surrogate estimates from: svaseq.

## This will pre-filter the input data using normalize_expt's: TRUE argument.

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

## Setting 40 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

oe_keeper <- list("up_control" = c("hdac5over", "controlover"))
oe_table <- combine_de_tables(oe_de, keepers = oe_keeper, excel = "excel/oe_table.xlsx")

## Deleting the file excel/oe_table.xlsx before writing the tables.

## Getting the genes >= 1.5 stdevs away from the mean of all.
## Getting the genes >= 1.5 stdevs away from the mean of all.
## Getting the genes >= 1.5 stdevs away from the mean of all.

oe_sig <- extract_significant_genes(oe_table, according_to = "deseq", p = 0.1, lfc = 0.38,
                                    excel = "excel/oe_sig.xlsx")

## Deleting the file excel/oe_sig.xlsx before writing the tables.

kd_up <- kd_sig[["deseq"]][["ups"]][["down_control"]]
nrow(kd_up)

## [1] 539

kd_down <- kd_sig[["deseq"]][["downs"]][["down_control"]]
nrow(kd_down)

## [1] 588

oe_up <- oe_sig[["deseq"]][["ups"]][["up_control"]]
nrow(oe_up)

## [1] 359

oe_down <- oe_sig[["deseq"]][["downs"]][["up_control"]]
nrow(oe_down)

## [1] 106

up_lst <- list(
    "kd_up" = rownames(kd_up),
    "oe_down" = rownames(oe_down))
up_inter <- Vennerable::Venn(up_lst)

down_lst <- list(
    "kd_down" = rownames(kd_down),
    "oe_up" = rownames(oe_up))
down_inter <- Vennerable::Venn(down_lst)

up_up <- list(
    "kd_up" = rownames(kd_up),
    "oe_up" = rownames(oe_up))
Vennerable::Venn(up_up)

## A Venn object on 2 sets named
## kd_up,oe_up 
##  00  10  01  11 
##   0 528 348  11

down_down <- list(
    "kd_down" = rownames(kd_down),
    "oe_down" = rownames(oe_down))
Vennerable::Venn(down_down)

## A Venn object on 2 sets named
## kd_down,oe_down 
##  00  10  01  11 
##   0 563  81  25

pander::pander(sessionInfo())

R version 4.2.0 (2022-04-22)

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

other attached packages: ruv(v.0.9.7.1), lme4(v.1.1-31), Matrix(v.1.5-3), BiocParallel(v.1.32.5), variancePartition(v.1.28.4), hpgltools(v.1.0), testthat(v.3.1.6), reticulate(v.1.28), SummarizedExperiment(v.1.28.0), GenomicRanges(v.1.50.2), GenomeInfoDb(v.1.34.9), IRanges(v.2.32.0), S4Vectors(v.0.36.1), MatrixGenerics(v.1.10.0), matrixStats(v.0.63.0), Biobase(v.2.58.0) and BiocGenerics(v.0.44.0)

loaded via a namespace (and not attached): rappdirs(v.0.3.3), rtracklayer(v.1.58.0), R.methodsS3(v.1.8.2), tidyr(v.1.3.0), ggplot2(v.3.4.1), clusterGeneration(v.1.3.7), bit64(v.4.0.5), knitr(v.1.42), R.utils(v.2.12.2), DelayedArray(v.0.24.0), data.table(v.1.14.6), KEGGREST(v.1.38.0), RCurl(v.1.98-1.10), doParallel(v.1.0.17), generics(v.0.1.3), preprocessCore(v.1.60.2), GenomicFeatures(v.1.50.4), callr(v.3.7.3), RhpcBLASctl(v.0.23-42), cowplot(v.1.1.1), usethis(v.2.1.6), RSQLite(v.2.2.20), shadowtext(v.0.1.2), bit(v.4.0.5), enrichplot(v.1.18.3), xml2(v.1.3.3), httpuv(v.1.6.8), assertthat(v.0.2.1), viridis(v.0.6.2), xfun(v.0.37), hms(v.1.1.2), jquerylib(v.0.1.4), IHW(v.1.26.0), evaluate(v.0.20), promises(v.1.2.0.1), DEoptimR(v.1.0-11), fansi(v.1.0.4), restfulr(v.0.0.15), progress(v.1.2.2), caTools(v.1.18.2), dbplyr(v.2.3.0), geneplotter(v.1.76.0), igraph(v.1.4.0), DBI(v.1.1.3), htmlwidgets(v.1.6.1), purrr(v.1.0.1), ellipsis(v.0.3.2), crosstalk(v.1.2.0), dplyr(v.1.1.0), backports(v.1.4.1), gprofiler2(v.0.2.1), annotate(v.1.76.0), aod(v.1.3.2), biomaRt(v.2.54.0), vctrs(v.0.5.2), remotes(v.2.4.2), cachem(v.1.0.6), withr(v.2.5.0), ggforce(v.0.4.1), HDO.db(v.0.99.1), robustbase(v.0.95-0), GenomicAlignments(v.1.34.0), treeio(v.1.22.0), fdrtool(v.1.2.17), prettyunits(v.1.1.1), DOSE(v.3.24.2), ape(v.5.6-2), lazyeval(v.0.2.2), crayon(v.1.5.2), genefilter(v.1.80.3), slam(v.0.1-50), labeling(v.0.4.2), edgeR(v.3.40.2), pkgconfig(v.2.0.3), tweenr(v.2.0.2), nlme(v.3.1-162), pkgload(v.1.3.2), devtools(v.2.4.5), rlang(v.1.0.6), lifecycle(v.1.0.3), miniUI(v.0.1.1.1), downloader(v.0.4), filelock(v.1.0.2), BiocFileCache(v.2.6.0), directlabels(v.2021.1.13), rprojroot(v.2.0.3), polyclip(v.1.10-4), graph(v.1.76.0), aplot(v.0.1.9), lpsymphony(v.1.26.3), boot(v.1.3-28.1), processx(v.3.8.0), png(v.0.1-8), viridisLite(v.0.4.1), rjson(v.0.2.21), bitops(v.1.0-7), R.oo(v.1.25.0), gson(v.0.0.9), KernSmooth(v.2.23-20), pander(v.0.6.5), Biostrings(v.2.66.0), blob(v.1.2.3), stringr(v.1.5.0), qvalue(v.2.30.0), remaCor(v.0.0.11), gridGraphics(v.0.5-1), scales(v.1.2.1), memoise(v.2.0.1), GSEABase(v.1.60.0), magrittr(v.2.0.3), plyr(v.1.8.8), gplots(v.3.1.3), zlibbioc(v.1.44.0), compiler(v.4.2.0), scatterpie(v.0.1.8), BiocIO(v.1.8.0), RColorBrewer(v.1.1-3), DESeq2(v.1.38.3), Rsamtools(v.2.14.0), cli(v.3.6.0), XVector(v.0.38.0), urlchecker(v.1.0.1), patchwork(v.1.1.2), ps(v.1.7.2), MASS(v.7.3-58.2), mgcv(v.1.8-41), tidyselect(v.1.2.0), stringi(v.1.7.12), highr(v.0.10), yaml(v.2.3.7), GOSemSim(v.2.24.0), locfit(v.1.5-9.7), ggrepel(v.0.9.3), grid(v.4.2.0), sass(v.0.4.5), fastmatch(v.1.1-3), tools(v.4.2.0), parallel(v.4.2.0), rstudioapi(v.0.14), foreach(v.1.5.2), gridExtra(v.2.3), Rtsne(v.0.16), farver(v.2.1.1), ggraph(v.2.1.0), digest(v.0.6.31), shiny(v.1.7.4), quadprog(v.1.5-8), Rcpp(v.1.0.10), broom(v.1.0.3), later(v.1.3.0), httr(v.1.4.4), AnnotationDbi(v.1.60.0), Rdpack(v.2.4), colorspace(v.2.1-0), brio(v.1.1.3), XML(v.3.99-0.13), fs(v.1.6.1), RBGL(v.1.74.0), yulab.utils(v.0.0.6), PROPER(v.1.30.0), tidytree(v.0.4.2), graphlayouts(v.0.8.4), ggplotify(v.0.1.0), plotly(v.4.10.1), sessioninfo(v.1.2.2), xtable(v.1.8-4), jsonlite(v.1.8.4), nloptr(v.2.0.3), ggtree(v.3.6.2), tidygraph(v.1.2.3), corpcor(v.1.6.10), ggfun(v.0.0.9), Vennerable(v.3.1.0.9000), R6(v.2.5.1), RUnit(v.0.4.32), profvis(v.0.3.7), pillar(v.1.8.1), htmltools(v.0.5.4), mime(v.0.12), glue(v.1.6.2), fastmap(v.1.1.0), minqa(v.1.2.5), clusterProfiler(v.4.6.0), codetools(v.0.2-19), fgsea(v.1.24.0), pkgbuild(v.1.4.0), mvtnorm(v.1.1-3), utf8(v.1.2.3), lattice(v.0.20-45), bslib(v.0.4.2), tibble(v.3.1.8), sva(v.3.46.0), pbkrtest(v.0.5.2), curl(v.5.0.0), gtools(v.3.9.4), zip(v.2.2.2), GO.db(v.3.16.0), openxlsx(v.4.2.5.2), survival(v.3.5-3), limma(v.3.54.1), rmarkdown(v.2.20), desc(v.1.4.2), munsell(v.0.5.0), fastcluster(v.1.2.3), GenomeInfoDbData(v.1.2.9), iterators(v.1.0.14), reshape2(v.1.4.4), gtable(v.0.3.1) and rbibutils(v.2.2.13)

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 c100d5666f5032d24c933739015d267ef651c323

## This is hpgltools commit: Wed Mar 1 09:50:14 2023 -0500: c100d5666f5032d24c933739015d267ef651c323

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

## Saving to index-v202303.rda.xz

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