202110 Revisiting some queries about stimulation responses at 4/8/12 hours

1 Introduction

Dr. Mosser approached us with a query about what genes are changing in human macrophages shortly after stimulation. We have a couple/few experiments which may help address this. Therefore, I think I will copy my master sample sheet and extract only the ~20 samples which are relevant rather than read in all ~400 samples and subset down to the 20 of interest.

However, the sheet he gave me has a subset of the samples from these experiments, I think I will keep them all and subset down in case it turns out we want to consider some of them.

sample_sheet <- "sample_sheets/mosser_samples_202110.xlsx"
file_column <- "hg3891salmon"

2 Gather annotation data

meta_info <- check_metadata_year(metadata = sample_sheet, column = file_column)

## Checking the creation time on the first count table.

## hs_annot <- load_biomart_annotations(year = meta_info[["year"]],
##                                      overwrite=TRUE)$annotation
hs_annot <- load_biomart_annotations(year = meta_info[["year"]], month = 10)$annotation

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

rownames(hs_annot) <- make.names(
  paste0(hs_annot[["ensembl_transcript_id"]], ".",
         hs_annot[["transcript_version"]]),
  unique = TRUE)
hs_tx_gene <- hs_annot[, c("ensembl_gene_id", "ensembl_transcript_id")]
hs_tx_gene[["id"]] <- rownames(hs_tx_gene)
hs_tx_gene <- hs_tx_gene[, c("id", "ensembl_gene_id")]
new_hs_annot <- hs_annot
rownames(new_hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique = TRUE)

3 Create expressionset

Including the samples I decided to keep from the two experiments in question, there are 94 remaining in these experiments.

The actual sample IDs which are of interest are:

• Kajal’s 4hr M-CSF LPS stimulated: HPGL0917, HPGL0918, HPGL0919, HPGL0920, HPGL0921 (these have the ‘state’ column set to ‘lps’ and ‘expt_time’ set to ‘t4h’)
• M-CSF LPS 12hr stimulated: HPGL0932, HPGL0933, HPGL0934, HPGL0935, HPGL0936 (these have state set to ‘lps’ and time to ‘t12h’.
• Medimmune grown ex-vivo: HPGL0732, HPGL0741, HPGL0733, HPGL0742, HPGL0735, HPGL0744, HPGL0737, HPGL0746, HPGL0739, and HPGL0748 (these have state set to ‘uninfected’ and time across a range from t4h->t24h).

stimulators <- create_expt(sample_sheet, file_column="hg3891salmon",
                           gene_info=new_hs_annot, tx_gene_map=hs_tx_gene)

## Reading the sample metadata.

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

## Matched 19803 annotations and counts.

## Bringing together the count matrix and gene information.

## The mapped IDs are not the rownames of your gene information, changing them now.

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

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

## The final expressionset has 19803 rows and 94 columns.

## Quick check to see if my check_metadata_year worked
## I know a priori that ENSG00000169245 has changed version multiple
## times since 2018 and that it was not being included if I use a
## newer set of annotations.
"ENSG00000169245" %in% rownames(exprs(stimulators))

## [1] TRUE

## YAY

time_state <- paste0(pData(stimulators)[["expttime"]], "_",
                     pData(stimulators)[["state"]])
stimv2 <- set_expt_conditions(stimulators, fact=time_state)

4 Before subsetting

I want to look at these two groups of samples together before I focus in on the times/samples that Dr. Mosser is interested in. My reasoning: I should be able to tell if they are pathologically batched more easily now than after taking a bunch of samples away.

time_factor <- factor(pData(stimv2)[["expttime"]],
                      levels=c("t0h", "t4h", "t6h", "t8h", "t12h", "t24h", "t48h"))
pData(stimv2)[["expttime"]] <- time_factor

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

## Removing 7218 low-count genes (12585 remaining).

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

plot_pca(stim_norm, size_column="expttime", plot_labels=FALSE)$plot

stim_nb <- normalize_expt(stimv2, norm="quant", convert="cpm", transform="log2",
                            batch="svaseq", filter=TRUE)

## Warning in normalize_expt(stimv2, norm = "quant", convert = "cpm", transform =
## "log2", : Quantile normalization and sva do not always play well together.

## Removing 7218 low-count genes (12585 remaining).

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

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

## Setting 5709 low elements to zero.

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

plot_pca(stim_nb, size_column="expttime", plot_labels=FALSE)$plot

Ok, so it takes a while of staring at it to get a sense of what is going on, but I think it is eventually pretty clear that the various other stimulation types are confusing things pretty significantly.

stimv3 <- stimv2 %>%
  subset_expt(subset="expttime!='t0h'") %>%
  subset_expt(subset="expttime!='t6h'") %>%
  subset_expt(subset="expttime!='t48h'") %>%
  subset_expt(subset="state!='lpsado'") %>%
  subset_expt(subset="state!='lpspge'") %>%
  subset_expt(subset="state!='il4'")

## subset_expt(): There were 94, now there are 91 samples.

## subset_expt(): There were 91, now there are 87 samples.

## subset_expt(): There were 87, now there are 85 samples.

## subset_expt(): There were 85, now there are 74 samples.

## subset_expt(): There were 74, now there are 63 samples.

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

pData(stimv3)[, c("sampleid", "expttime", "condition")]

##          sampleid expttime       condition
## hpgl0732 hpgl0732      t4h  t4h_uninfected
## hpgl0733 hpgl0733      t4h         t4h_lps
## hpgl0734 hpgl0734      t8h  t8h_uninfected
## hpgl0735 hpgl0735      t8h         t8h_lps
## hpgl0736 hpgl0736     t12h t12h_uninfected
## hpgl0737 hpgl0737     t12h        t12h_lps
## hpgl0738 hpgl0738     t24h t24h_uninfected
## hpgl0739 hpgl0739     t24h        t24h_lps
## hpgl0741 hpgl0741      t4h  t4h_uninfected
## hpgl0742 hpgl0742      t4h         t4h_lps
## hpgl0743 hpgl0743      t8h  t8h_uninfected
## hpgl0744 hpgl0744      t8h         t8h_lps
## hpgl0745 hpgl0745     t12h t12h_uninfected
## hpgl0746 hpgl0746     t12h        t12h_lps
## hpgl0747 hpgl0747     t24h t24h_uninfected
## hpgl0748 hpgl0748     t24h        t24h_lps
## hpgl0754 hpgl0754      t4h  t4h_uninfected
## hpgl0755 hpgl0755      t4h         t4h_lps
## hpgl0756 hpgl0756      t8h  t8h_uninfected
## hpgl0757 hpgl0757      t8h         t8h_lps
## hpgl0758 hpgl0758     t12h t12h_uninfected
## hpgl0759 hpgl0759     t12h        t12h_lps
## hpgl0760 hpgl0760     t24h t24h_uninfected
## hpgl0761 hpgl0761     t24h        t24h_lps
## hpgl0912 hpgl0912      t4h  t4h_uninfected
## hpgl0913 hpgl0913      t4h  t4h_uninfected
## hpgl0914 hpgl0914      t4h  t4h_uninfected
## hpgl0915 hpgl0915      t4h  t4h_uninfected
## hpgl0916 hpgl0916      t4h  t4h_uninfected
## hpgl0917 hpgl0917      t4h         t4h_lps
## hpgl0918 hpgl0918      t4h         t4h_lps
## hpgl0919 hpgl0919      t4h         t4h_lps
## hpgl0920 hpgl0920      t4h         t4h_lps
## hpgl0921 hpgl0921      t4h         t4h_lps
## hpgl0932 hpgl0932     t12h        t12h_lps
## hpgl0933 hpgl0933     t12h        t12h_lps
## hpgl0934 hpgl0934     t12h        t12h_lps
## hpgl0935 hpgl0935     t12h        t12h_lps
## hpgl0936 hpgl0936     t12h        t12h_lps
## hpgl0937 hpgl0937     t12h t12h_uninfected
## hpgl0938 hpgl0938     t12h t12h_uninfected
## hpgl0939 hpgl0939     t12h t12h_uninfected
## hpgl0940 hpgl0940     t12h t12h_uninfected
## hpgl0941 hpgl0941     t12h t12h_uninfected
## hpgl0942 hpgl0942      t4h  t4h_uninfected
## hpgl0943 hpgl0943      t4h  t4h_uninfected
## hpgl0944 hpgl0944      t4h  t4h_uninfected
## hpgl0945 hpgl0945      t4h  t4h_uninfected
## hpgl0946 hpgl0946      t4h         t4h_lps
## hpgl0947 hpgl0947      t4h         t4h_lps
## hpgl0948 hpgl0948      t4h         t4h_lps
## hpgl0949 hpgl0949      t4h         t4h_lps
## hpgl0950 hpgl0950      t4h         t4h_lps
## hpgl0961 hpgl0961     t12h        t12h_lps
## hpgl0962 hpgl0962     t12h        t12h_lps
## hpgl0963 hpgl0963     t12h        t12h_lps
## hpgl0964 hpgl0964     t12h        t12h_lps
## hpgl0965 hpgl0965     t12h        t12h_lps
## hpgl0966 hpgl0966     t12h t12h_uninfected
## hpgl0967 hpgl0967     t12h t12h_uninfected
## hpgl0968 hpgl0968     t12h t12h_uninfected
## hpgl0969 hpgl0969     t12h t12h_uninfected
## hpgl0970 hpgl0970     t12h t12h_uninfected

time_state <- paste0(pData(stimv3)[["expttime"]], "_",
                     pData(stimv3)[["state"]])
stimv3 <- set_expt_conditions(stimv3, fact=time_state)

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

## Removing 7366 low-count genes (12437 remaining).

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

plot_pca(stimv3_norm, plot_labels=FALSE)$plot

stimv3_nb <- normalize_expt(stimv3, transform="log2", convert="cpm", norm="quant", filter=TRUE, batch="svaseq")

## Warning in normalize_expt(stimv3, transform = "log2", convert = "cpm", norm =
## "quant", : Quantile normalization and sva do not always play well together.

## Removing 7366 low-count genes (12437 remaining).

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

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

## Setting 4581 low elements to zero.

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

plot_pca(stimv3_nb)$plot

plot_pca(stimv3_nb, size_column="expttime")$plot

5 Smallest group

After a conversation with Najib, he suggested working with a very limited set of 4 and 8 hour samples, probably exclusively from the medimmune experiment.

stimv4 <- stimv3 %>%
  subset_expt(subset="expttime=='t4h'|expttime=='t8h'")

## subset_expt(): There were 63, now there are 31 samples.

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

## Removing 7745 low-count genes (12058 remaining).

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

plot_pca(stimv4_norm)$plot

stimv4_nb <- normalize_expt(stimv4, transform="log2", convert="cpm", norm="quant",
                            batch="svaseq", filter=TRUE)

## Warning in normalize_expt(stimv4, transform = "log2", convert = "cpm", norm =
## "quant", : Quantile normalization and sva do not always play well together.

## Removing 7745 low-count genes (12058 remaining).

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

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

## Setting 1585 low elements to zero.

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

plot_pca(stimv4_nb)$plot

5.1 Consider only 4 and 12 hours

stimv412 <- stimv3 %>%
  subset_expt(subset="expttime=='t4h'|expttime=='t12h'")

## subset_expt(): There were 63, now there are 51 samples.

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

## Removing 7497 low-count genes (12306 remaining).

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

plot_pca(stimv412_norm)$plot

stimv412_nb <- normalize_expt(stimv412, transform="log2", convert="cpm", norm="quant",
                              batch="svaseq", filter=TRUE)

## Warning in normalize_expt(stimv412, transform = "log2", convert = "cpm", :
## Quantile normalization and sva do not always play well together.

## Removing 7497 low-count genes (12306 remaining).

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

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

## Setting 3045 low elements to zero.

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

plot_pca(stimv412_nb)$plot

st_filt <- normalize_expt(stimv412, filter="simple")

## Removing 2617 low-count genes (17186 remaining).

initial <- all_pairwise(
    st_filt, model_batch=TRUE, force=TRUE,
    extra_contrasts="t12v4_lpsvuninf=(t12h_lps-t12h_uninfected)-(t4h_lps-t4h_uninfected)")

## Plotting a PCA before surrogate/batch inclusion.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
    "12v4_uninfected" = c("t12h_uninfected", "t4h_uninfected"),
    "12v4_lps" = c("t12h_lps", "t4h_lps"),
    "lpsvuninf_4hr" = c("t4h_lps", "t4h_uninfected"),
    "lpsvuninf_12hr" = c("t12h_lps", "t12h_uninfected"),
    "t12v4_lpsvuninf" = NULL)
funkytown <- combine_de_tables(initial, keepers=keepers, excel="lps_vs_infected.xlsx")

## Deleting the file lps_vs_infected.xlsx before writing the tables.

## Assuming t12v4_lpsvuninf is a table from extra_contrasts.

## Unable to find the table in the set of possible tables.

## The possible tables are: t12h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_lps, t4h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_uninfected, t4h_uninfected_vs_t12h_uninfected, t4h_uninfected_vs_t4h_lps

## Unable to find the table in the set of possible tables.

## The possible tables are: t12h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_lps, t4h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_uninfected, t4h_uninfected_vs_t12h_uninfected, t4h_uninfected_vs_t4h_lps, t12v4_lpsvuninf

## Unable to find the table in the set of possible tables.

## The possible tables are: t12h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_lps, t4h_uninfected_vs_t12h_lps, t4h_lps_vs_t12h_uninfected, t4h_uninfected_vs_t12h_uninfected, t4h_uninfected_vs_t4h_lps

Take top n genes from LPSvNS at 4 hours, plot their abundance at 4 and 12 hours as a scatter plot.

topn <- 50

t4hr_start <- funkytown[["data"]][["lpsvuninf_4hr"]]
descending_idx <- order(t4hr_start[["edger_logfc"]], decreasing = TRUE)
t4hr_start <- t4hr_start[descending_idx, ]

t12hr_start <- funkytown[["data"]][["lpsvuninf_12hr"]]
t12hr_start <- t12hr_start[descending_idx, ]

topn_4hr <- head(t4hr_start, n=topn)
topn_12hr <- head(t12hr_start, n=topn)

combined <- as.data.frame(cbind(topn_4hr[["edger_logfc"]], topn_12hr[["edger_logfc"]]))
rownames(combined) <- make.names(head(t12hr_start[["hgncsymbol"]],
                                      n=topn), unique=TRUE)

plotted <- plot_scatter(combined) +
  scale_x_continuous(limits = c(0, 12)) +
  scale_y_continuous(limits = c(0, 12)) +
  ggplot2::geom_abline(slope=1, intercept=0)

## ggplotly(plotted)

all <- ggplt(plotted, filename="fifty_genes.html")

5.2 Comparing to Dr. Mosser’s previous table

Dr. Mosser provided me a table of previous results, the first two genes listed are CXCL10 and CXCL11 (ENSG00000169245 and 169248). I am not finding these two genes in any of my result tables. I think they got filtered out, even when I use a less restrictive filter. Let us try to understand what is going on.

"ENSG00000169245" %in% rownames(exprs(stimv412))

## [1] TRUE

## Hmm ok...
names <- fData(stimv412)[["hgnc_symbol"]]
cxcl_names <- grep(pattern="CXCL", x=names)
names[cxcl_names]

##  [1] "CXCL2"  "CXCL12" "CXCL6"  "CXCL9"  "CXCL14" "CXCL13" "CXCL16" "CXCL3" 
##  [9] "CXCL5"  "CXCL1"  "CXCL10" "CXCL11" "CXCL8"  "CXCL17"

## So, does this mean that something changed with the annotations such
## that CXCL10 is no longer included, or was it not included in the
## salmon mapping?  Well, I still have my annotations loaded, lets
## check there first...
names <- hs_annot[["hgnc_symbol"]]
cxcl_names <- grep(pattern="CXCL", x=names)
names[cxcl_names]

##  [1] "CXCL6"   "CXCL9"   "CXCL13"  "CXCL16"  "CXCL3"   "CXCL5"   "CXCL2"  
##  [8] "CXCL10"  "CXCL11"  "CXCL8"   "CXCL14"  "CXCL17"  "CXCL12"  "CXCL12" 
## [15] "CXCL12"  "CXCL1"   "CXCL12"  "CXCL12"  "CXCL12"  "CXCL8"   "CXCL8"  
## [22] "CXCL12"  "CXCL12"  "CXCL1P1" "CXCL3"   "CXCL6"   "CXCL11"  "CXCL13" 
## [29] "CXCL2"   "CXCL1"   "CXCL2"   "CXCL3"   "CXCL3"   "CXCL14"  "CXCL6"  
## [36] "CXCL16"  "CXCL16"  "CXCL16"  "CXCL16"  "CXCL17"

## ok, so the salmon quantification did not include CXCL10.

## Lets look at the salmon tables...
## I see that CXCL10 has 1 transcript: ENST00000306602.3
## Looking at one of the relevant count tables (hpgl0968) for hg38_99,
## I see that this transcript was included (but with 0 counts).
## Looking at other samples, I see that it was observed a small number
## of times -- but it definitely _was_ included.

## Let us check that the tx_gene_map is ok...
"ENST00000306602.3" %in% rownames(hs_tx_gene)

## [1] FALSE

## ok, so it is definitely there, double-check the annotations.
"ENSG00000169245" %in% rownames(new_hs_annot)

## [1] TRUE

stimulators <- create_expt(sample_sheet, file_column="hg3891salmon")

## Reading the sample metadata.

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

## Matched 86451 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final expressionset has 86451 rows and 94 columns.

"ENST00000306602.2" %in% rownames(exprs(stimulators))

## [1] TRUE

5.3 Compare 12 hour subtraction to 4 hour

t12hr_df <- funkytown[["data"]][["lpsvuninf_12hr"]]
t4hr_df <- funkytown[["data"]][["lpsvuninf_4hr"]]

cpm_df <- merge(t12hr_df, t4hr_df, by="row.names")
rownames(cpm_df) <- make.names(cpm_df[["hgncsymbol.x"]], unique=TRUE)
cpm_df[["Row.names"]] <- NULL

test <- plot_scatter(cpm_df[, c("deseq_logfc.x", "deseq_logfc.y")])

6 Clarification

Take the top-n most increased genes at 4 hours after LPS stimulation. Given this set, ask the question: what happened to them?

One way to provide this answer is to pull out their average cpm/rpkm values at each of the 4 relevant times/conditions and make a table.

topn <- 100
lps_order <- order(t4hr_df[["deseq_logfc"]], decreasing=TRUE)
lps_reordered <- head(t4hr_df[lps_order, ], n=topn)

ids_of_interest <- rownames(lps_reordered)
annot_interest <- lps_reordered[, c("hgncsymbol", "description")]

## We can cheat a little and get the average cpm just to get an initial sense
t4h_lpsvcontrol_cpm <- lps_reordered[, c("basic_nummed", "basic_denmed")]
t12h_lpsvcontrol_cpm <- t12hr_df[ids_of_interest, c("basic_nummed", "basic_denmed")]

combined_lpsvcontrol <- merge(t4h_lpsvcontrol_cpm, t12h_lpsvcontrol_cpm, by="row.names")
combined_lpsvcontrol <- merge(annot_interest, combined_lpsvcontrol, by.x="row.names", by.y="Row.names")
colnames(combined_lpsvcontrol) <- c("rownames", "hgncsymbol", "description",
                                    "t4h_lps", "t4h_control", "t12h_lps", "t12h_control")
written <- write_xlsx(combined_lpsvcontrol, excel="cpm_lpsvcontrol_t4hvt12h.xlsx")

wanted <- c("hpgl0732", "hpgl0741", "hpgl0733", "hpgl0742", "hpgl0734", "hpgl0743", "hpgl0735", "hpgl0744")
stimv5 <- stimv2 %>%
  subset_expt(ids=wanted)

## subset_expt(): There were 94, now there are 8 samples.

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

## Removing 8840 low-count genes (10963 remaining).

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

plot_pca(stimv5_norm)$plot

pander::pander(sessionInfo())

R version 4.1.2 (2021-11-01)

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

other attached packages: ruv(v.0.9.7.1), ggplot2(v.3.3.5), hpgltools(v.1.0), testthat(v.3.1.1), SummarizedExperiment(v.1.24.0), GenomicRanges(v.1.46.1), GenomeInfoDb(v.1.30.0), IRanges(v.2.28.0), S4Vectors(v.0.32.3), MatrixGenerics(v.1.6.0), matrixStats(v.0.61.0), Biobase(v.2.54.0) and BiocGenerics(v.0.40.0)

loaded via a namespace (and not attached): rappdirs(v.0.3.3), rtracklayer(v.1.54.0), R.methodsS3(v.1.8.1), tidyr(v.1.1.4), bit64(v.4.0.5), knitr(v.1.37), R.utils(v.2.11.0), DelayedArray(v.0.20.0), data.table(v.1.14.2), KEGGREST(v.1.34.0), RCurl(v.1.98-1.5), doParallel(v.1.0.16), generics(v.0.1.1), GenomicFeatures(v.1.46.1), preprocessCore(v.1.56.0), callr(v.3.7.0), usethis(v.2.1.5), RSQLite(v.2.2.9), shadowtext(v.0.1.0), bit(v.4.0.4), tzdb(v.0.2.0), enrichplot(v.1.14.1), xml2(v.1.3.3), httpuv(v.1.6.4), assertthat(v.0.2.1), viridis(v.0.6.2), xfun(v.0.29), tximport(v.1.22.0), hms(v.1.1.1), jquerylib(v.0.1.4), evaluate(v.0.14), promises(v.1.2.0.1), IHW(v.1.22.0), DEoptimR(v.1.0-9), fansi(v.0.5.0), restfulr(v.0.0.13), progress(v.1.2.2), caTools(v.1.18.2), dbplyr(v.2.1.1), igraph(v.1.2.10), DBI(v.1.1.2), geneplotter(v.1.72.0), htmlwidgets(v.1.5.4), purrr(v.0.3.4), ellipsis(v.0.3.2), crosstalk(v.1.2.0), dplyr(v.1.0.7), backports(v.1.4.1), annotate(v.1.72.0), biomaRt(v.2.50.1), vctrs(v.0.3.8), remotes(v.2.4.2), cachem(v.1.0.6), withr(v.2.4.3), ggforce(v.0.3.3), robustbase(v.0.93-9), vroom(v.1.5.7), GenomicAlignments(v.1.30.0), treeio(v.1.18.1), fdrtool(v.1.2.17), prettyunits(v.1.1.1), DOSE(v.3.20.1), ape(v.5.6), lazyeval(v.0.2.2), crayon(v.1.4.2), genefilter(v.1.76.0), edgeR(v.3.36.0), pkgconfig(v.2.0.3), slam(v.0.1-49), labeling(v.0.4.2), tweenr(v.1.0.2), nlme(v.3.1-153), pkgload(v.1.2.4), devtools(v.2.4.3), rlang(v.0.4.12), lifecycle(v.1.0.1), downloader(v.0.4), filelock(v.1.0.2), BiocFileCache(v.2.2.0), rprojroot(v.2.0.2), polyclip(v.1.10-0), graph(v.1.72.0), Matrix(v.1.3-4), aplot(v.0.1.1), lpsymphony(v.1.22.0), boot(v.1.3-28), processx(v.3.5.2), png(v.0.1-7), viridisLite(v.0.4.0), rjson(v.0.2.20), bitops(v.1.0-7), R.oo(v.1.24.0), KernSmooth(v.2.23-20), pander(v.0.6.4), Biostrings(v.2.62.0), blob(v.1.2.2), stringr(v.1.4.0), qvalue(v.2.26.0), readr(v.2.1.1), gridGraphics(v.0.5-1), scales(v.1.1.1), memoise(v.2.0.1), magrittr(v.2.0.1), plyr(v.1.8.6), gplots(v.3.1.1), zlibbioc(v.1.40.0), compiler(v.4.1.2), scatterpie(v.0.1.7), BiocIO(v.1.4.0), RColorBrewer(v.1.1-2), lme4(v.1.1-27.1), DESeq2(v.1.34.0), Rsamtools(v.2.10.0), cli(v.3.1.0), XVector(v.0.34.0), patchwork(v.1.1.1), ps(v.1.6.0), MASS(v.7.3-54), mgcv(v.1.8-38), tidyselect(v.1.1.1), stringi(v.1.7.6), highr(v.0.9), yaml(v.2.2.1), GOSemSim(v.2.20.0), locfit(v.1.5-9.4), ggrepel(v.0.9.1), grid(v.4.1.2), sass(v.0.4.0), fastmatch(v.1.1-3), tools(v.4.1.2), parallel(v.4.1.2), rstudioapi(v.0.13), foreach(v.1.5.1), gridExtra(v.2.3), farver(v.2.1.0), ggraph(v.2.0.5), digest(v.0.6.29), shiny(v.1.7.1), Rcpp(v.1.0.7), broom(v.0.7.10), later(v.1.3.0), httr(v.1.4.2), AnnotationDbi(v.1.56.2), colorspace(v.2.0-2), XML(v.3.99-0.8), fs(v.1.5.2), splines(v.4.1.2), RBGL(v.1.70.0), yulab.utils(v.0.0.4), PROPER(v.1.26.0), tidytree(v.0.3.6), graphlayouts(v.0.7.2), ggplotify(v.0.1.0), plotly(v.4.10.0), sessioninfo(v.1.2.2), xtable(v.1.8-4), jsonlite(v.1.7.2), nloptr(v.1.2.2.3), ggtree(v.3.2.1), tidygraph(v.1.2.0), corpcor(v.1.6.10), ggfun(v.0.0.4), R6(v.2.5.1), Vennerable(v.3.1.0.9000), pillar(v.1.6.4), htmltools(v.0.5.2), mime(v.0.12), glue(v.1.6.0), fastmap(v.1.1.0), minqa(v.1.2.4), clusterProfiler(v.4.2.1), BiocParallel(v.1.28.3), codetools(v.0.2-18), fgsea(v.1.20.0), pkgbuild(v.1.3.1), utf8(v.1.2.2), lattice(v.0.20-45), bslib(v.0.3.1), tibble(v.3.1.6), sva(v.3.42.0), pbkrtest(v.0.5.1), curl(v.4.3.2), gtools(v.3.9.2), zip(v.2.2.0), GO.db(v.3.14.0), openxlsx(v.4.2.5), survival(v.3.2-13), limma(v.3.50.0), rmarkdown(v.2.11), desc(v.1.4.0), munsell(v.0.5.0), DO.db(v.2.9), GenomeInfoDbData(v.1.2.7), iterators(v.1.0.13), variancePartition(v.1.24.0), reshape2(v.1.4.4) and gtable(v.0.3.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 1c414fc2f002f2eb5020f3753fa29439b55b8a33

## This is hpgltools commit: Tue Jan 18 14:27:47 2022 -0500: 1c414fc2f002f2eb5020f3753fa29439b55b8a33

message(paste0("Saving to ", savefile))

## Error in paste0("Saving to ", savefile): object 'savefile' not found

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

## Error in file.path(getwd(), directory, filename): object 'savefile' not found

## Use me to reload the data
tmp <- loadme(filename = savefile)