TMRC2 Infecting macrophages: 202203

1 TODO

1. Compare these samples to previous in-vitro L.major infected macrophages.
2. Look at shared/unique genes with respect to antimonial treatment.
3. Send along rpkm counts/tables.
4. Perform enrichment analyses of high/low genes vs. DE genes.
5. Compare to MSigDB and/or GSVA

2 Introduction

3 Annotation

We take the annotation data from ensembl’s biomart instance. The genome which was used to map the data was hg38 revision 91. My default when using biomart is to load the data from 1 year before the current date, which provides annotations which match hg38 91 almost perfectly.

hs_annot <- load_biomart_annotations(year="2020", month="jan")

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

hs_annot <- hs_annot[["annotation"]]
hs_annot[["transcript"]] <- paste0(rownames(hs_annot), ".", hs_annot[["version"]])
rownames(hs_annot) <- make.names(hs_annot[["ensembl_gene_id"]], unique=TRUE)
tx_gene_map <- hs_annot[, c("transcript", "ensembl_gene_id")]

4 Sample Estimation

I used two mapping methods for this data, hisat2 and salmon. Most analyses use hisat2, which is a more traditional map-and-count method. In contrast, salmon uses what may be thought of as a indexed voting method (so that multi-matches are discounted and the votes split among all matches). Salmon also required a pre-existing database of known transcripts (though later versions may actually use mapping from things like hisat), while hisat uses the genome and a database of known transcripts (and optionally can search for splicing junctions to find new transcripts).

4.1 Generate expressionsets

The first thing to note is the large range in coverage. There are multiple samples with coverage which is too low to use. These will be removed shortly.

sanitize_columns <- c("donor", "drug", "condition", "batch", "macrophagetreatment",
                      "macrophageinfectionzymodeme")
macr <- create_expt("sample_sheets/tmrc2_macrophage_samples_202203.xlsx",
                    file_column = "hg38100hisatfile",
                    savefile = "tmrc2_macrophage.rda",
                    gene_info = hs_annot) %>%
  exclude_genes_expt(column = "gene_biotype", method = "keep",
                     pattern = "protein_coding") %>%
  subset_expt(nonzero = 11000) %>%
  set_expt_factors(columns = sanitize_columns) %>%
  set_expt_factors(columns=sanitize_columns, class="factor") %>%
  set_expt_batches(fact="macrophagetreatment") %>%
  set_expt_conditions(fact="macrophageinfectionzymodeme")

## Reading the sample metadata.

## The sample definitions comprises: 28 rows(samples) and 67 columns(metadata fields).

## Matched 21447 annotations and counts.

## Bringing together the count matrix and gene information.

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

## The final expressionset has 21481 rows and 28 columns.

## Before removal, there were 21481 genes, now there are 19923.

## The samples (and read coverage) removed when filtering 11000 non-zero genes are:

## named numeric(0)

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

5 Macrophages

These samples are rather different from all of the others. The following section is therefore written primarily in response to a separate set of emails from Olga and Maria Adelaida; here is a snippet:

Dear all, about the samples corresponding to infected macrophages with three sensitive (2.2) and three resistant (2.3) clinical strains of L. (V.) panamensis, I send you the results of parasite burden by detection of 7SLRNA. I think these results are interesting, but the sample size is very small. Doctor Najib or Trey could you please send me the quality data and PCA analysis of these samples?

and

Hi Doctor, thank you. These samples corresponding to primary macrophages infected with clinical strains 2.2 (n=3) and 2.3 (n = 3). These information is in the file: TMRC project 3: excel Host TMRC3 v1.1, rows 137 to 150.

Thus I added 3 columns to the end of the sample sheet which seek to include this information. The first is ‘drug’ and encodes both the infection state (no for the two controls and yes for everything else), the second is zymodeme which I took from the tmrc2 sample sheet, the last is drug, which is either no or sb.

table(pData(macr)[["condition"]])

## 
## none z2.2 z2.3 
##    4   12   12

table(pData(macr)[["batch"]])

## 
##      inf   inf_sb    uninf uninf_sb 
##       12       12        2        2

table(pData(macr)[["macrophagetreatment"]])

## 
##      inf   inf_sb    uninf uninf_sb 
##       12       12        2        2

table(pData(macr)[["macrophageinfectionzymodeme"]])

## 
## none z2.2 z2.3 
##    4   12   12

macr_libsize <- plot_libsize(macr, yscale=TRUE)
dev <- pp(file="macrophage_side_experiment/macr_libsize.png")
macr_libsize$plot
closed <- dev.off()
macr_libsize$plot

macr_nonzero <- plot_nonzero(macr)
dev <- pp(file="macrophage_side_experiment/macr_nonzero.png")
macr_nonzero$plot
closed <- dev.off()
macr_nonzero$plot

macr_norm <- sm(normalize_expt(macr, transform="log2", norm="quant", convert="cpm", filter=TRUE))
norm_pca <- plot_pca(macr_norm, plot_labels=FALSE)
norm_pca$plot

dev <- pp(file="macrophage_side_experiment/pca_colors_zymodeme.png")
norm_pca$plot
closed <- dev.off()
norm_pca$plot

treatment_macr <- set_expt_conditions(macr, fact="macrophagetreatment") %>%
  set_expt_batches(fact="macrophageinfectionzymodeme")
treatment_norm <- sm(normalize_expt(treatment_macr, transform="log2", norm="quant", convert="cpm", filter=TRUE))
treatment_pca <- plot_pca(treatment_norm, plot_labels=FALSE)
treatment_pca$plot

dev <- pp(file="macrophage_side_experiment/pca_colors_treatment.png")
treatment_pca$plot
closed <- dev.off()
treatment_pca$plot

plot_3d_pca(norm_pca)$plot

macr_nb <- sm(normalize_expt(macr, norm = "quant", convert = "cpm",
                             transform = "log2", batch = "svaseq", filter = TRUE))
nb_pca <- plot_pca(macr_nb, plot_labels=FALSE)
dev <- pp(file="macrophage_side_experiment/normbatch_pca.png")
nb_pca$plot
closed <- dev.off()
nb_pca$plot

5.1 Write the data

macr_written <- sm(write_expt(macr, excel="macrophage_side_experiment/macrophage_expt.xlsx"))

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

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

5.2 Perform DE

There were a couple of contrasts explicitly requested for this data:

• Macrophages infected with 2.3 without SbV vs macrophagues infected with 2.2 without SbV
• Macrophages infected with 2.3 without SbV vs macrophages uninfected (M0)

Unfortunately, I think to really answer these questions we will require more replicates of a few of these conditions, most notably the uninfected samples.

treatment_macr <- set_expt_conditions(macr, fact="macrophagetreatment") %>%
  set_expt_batches(fact="macrophageinfectionzymodeme")
##tmp <- normalize_expt(macr, filter=TRUE)
##zy_de <- deseq_pairwise(tmp, model_batch="svaseq")
##zy_ed <- edger_pairwise(tmp, model_batch="svaseq")
combined_condition <- paste0(pData(macr)[["macrophagetreatment"]], "_",
                             pData(macr)[["macrophageinfectionzymodeme"]])
pData(macr[["expressionset"]])[["combined"]] <- combined_condition
macr <- set_expt_conditions(macr, fact = "combined")
zymo_de <- all_pairwise(macr, model_batch = "svaseq", filter = TRUE,
                        do_ebseq = FALSE)

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

## Setting 801 low elements to zero.

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

## Finished running DE analyses, collecting outputs.

## Comparing analyses.

keepers <- list(
    "z22_uninf" = c("infz22", "uninfnone"),
    "z23_uninf" = c("infz23", "uninfnone"),
    "z22_sbuninf" = c("infsbz22", "uninfsbnone"),
    "z23_sbuninf" = c("infsbz23", "uninfsbnone"),
    "zymo_nodrug" = c("infz23", "infz22"),
    "zymo_sb" = c("infsbz23", "infsbz22"))

zymo_table <- combine_de_tables(
    zymo_de, keepers=keepers,
    excel=glue::glue("excel/macrophage_de-v{ver}.xlsx"))

dev <- pp(file="images/z22_vs_uninfected_ma.png")
zymo_table$plots$z22_uninf$deseq_ma_plots$plot
closed <- dev.off()
zymo_table$plots$z22_uninf$deseq_ma_plots$plot

dev <- pp(file="images/z22_vs_sb_uninfected_ma.png")
zymo_table$plots$z22_sbuninf$deseq_ma_plots$plot
closed <- dev.off()
zymo_table$plots$z22_sbuninf$deseq_ma_plots$plot

dev <- pp(file="images/z23_vs_uninfected_ma.png")
zymo_table$plots$z23_uninf$deseq_ma_plots$plot
closed <- dev.off()
zymo_table$plots$z23_uninf$deseq_ma_plots$plot

dev <- pp(file="images/z23_vs_sb_uninfected_ma.png")
zymo_table$plots$z23_sbuninf$deseq_ma_plots$plot
closed <- dev.off()
zymo_table$plots$z23_sbuninf$deseq_ma_plots$plot

zymo_sig <- extract_significant_genes(zymo_table, according_to="deseq")

z22_ma <- plot_ma_de(zymo_table[["data"]][["z22_uninf"]], expr_col="deseq_basemean",
                 fc_col="deseq_logfc", p_col="deseq_adjp")
z22_ma$plot

z22_uninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z22_uninf"]]
dim(z22_uninf_up_sig)

## [1] 436  47

z22_uninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z22_uninf"]]
dim(z22_uninf_down_sig)

## [1]  6 47

z22_sbuninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z22_sbuninf"]]
dim(z22_sbuninf_up_sig)

## [1] 20 47

z22_sbuninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z22_sbuninf"]]
dim(z22_sbuninf_down_sig)

## [1]  0 47

z23_ma <- plot_ma_de(zymo_table[["data"]][["z23_uninf"]], expr_col="deseq_basemean",
                     fc_col="deseq_logfc", p_col="deseq_adjp")
z23_ma$plot

z23_uninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z23_uninf"]]
dim(z23_uninf_up_sig)

## [1] 643  47

z23_uninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z23_uninf"]]
dim(z23_uninf_down_sig)

## [1] 372  47

z23_sbuninf_up_sig <- zymo_sig[["deseq"]][["ups"]][["z23_sbuninf"]]
dim(z23_sbuninf_up_sig)

## [1] 530  47

z23_sbuninf_down_sig <- zymo_sig[["deseq"]][["downs"]][["z23_sbuninf"]]
dim(z23_sbuninf_down_sig)

## [1] 323  47

5.3 Compare groups of genes

z22_z23_nosb_groups <- list(
    "z22_up" = rownames(z22_uninf_up_sig),
    "z23_up" = rownames(z23_uninf_up_sig))
z22_z23_nosb_venn <- Vennerable::Venn(Sets = z22_z23_nosb_groups)
dev <- pp(file="images/z22_z23_nosb_venn.png")
z22_z23_nosb_venn_plot <- Vennerable::plot(z22_z23_nosb_venn, doWeights = FALSE)
closed <- dev.off()
Vennerable::plot(z22_z23_nosb_venn, doWeights = FALSE)

annot <- fData(macr)
z23_nosb_unique <- z22_z23_nosb_venn@IntersectionSets[["01"]]
wanted <- rownames(annot) %in% z23_nosb_unique
z23_wanted <- annot[wanted, ]
written <- write_xlsx(z23_wanted, excel="excel/z23_nosb_unique_genes.xlsx")

z22_nosb_unique <- z22_z23_nosb_venn@IntersectionSets[["10"]]
wanted <- rownames(annot) %in% z22_nosb_unique
z22_wanted <- annot[wanted, ]
written <- write_xlsx(z22_wanted, excel="excel/z22_nosb_unique_genes.xlsx")

z22_nosb_shared <- z22_z23_nosb_venn@IntersectionSets[["11"]]
wanted <- rownames(annot) %in% z22_nosb_shared
shared_wanted <- annot[wanted, ]
written <- write_xlsx(shared_wanted, excel="excel/z22_z23_nosb_shared_genes.xlsx")

z22_unique_gp <- simple_gprofiler(rownames(z22_wanted))

## Performing gProfiler GO search of 333 genes against hsapiens.

## GO search found 74 hits.

## Performing gProfiler KEGG search of 333 genes against hsapiens.

## KEGG search found 3 hits.

## Performing gProfiler REAC search of 333 genes against hsapiens.

## REAC search found 13 hits.

## Performing gProfiler MI search of 333 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 333 genes against hsapiens.

## TF search found 43 hits.

## Performing gProfiler CORUM search of 333 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 333 genes against hsapiens.

## HP search found 0 hits.

dev <- pp(file="images/z22_unique_bp.png", height=9, width=6)
z22_unique_gp[["pvalue_plots"]][["bpp_plot_over"]]
closed <- dev.off()
z22_unique_gp[["pvalue_plots"]][["bpp_plot_over"]]

dev <- pp(file="images/z22_unique_kegg.png")
z22_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]
closed <- dev.off()
z22_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]

dev <- pp(file="images/z22_unique_reac.png")
z22_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]
closed <- dev.off()
z22_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]

z23_unique_gp <- simple_gprofiler(rownames(z23_wanted))

## Performing gProfiler GO search of 540 genes against hsapiens.

## GO search found 459 hits.

## Performing gProfiler KEGG search of 540 genes against hsapiens.

## KEGG search found 16 hits.

## Performing gProfiler REAC search of 540 genes against hsapiens.

## REAC search found 27 hits.

## Performing gProfiler MI search of 540 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 540 genes against hsapiens.

## TF search found 85 hits.

## Performing gProfiler CORUM search of 540 genes against hsapiens.

## CORUM search found 2 hits.

## Performing gProfiler HP search of 540 genes against hsapiens.

## HP search found 0 hits.

dev <- pp(file="images/z23_unique_bp.png", height=9, width=6)
z23_unique_gp[["pvalue_plots"]][["bpp_plot_over"]]
closed <- dev.off()
z23_unique_gp[["pvalue_plots"]][["bpp_plot_over"]]

dev <- pp(file="images/z23_unique_kegg.png")
z23_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]
closed <- dev.off()
z23_unique_gp[["pvalue_plots"]][["kegg_plot_over"]]

dev <- pp(file="images/z23_unique_reac.png")
z23_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]
closed <- dev.off()
z23_unique_gp[["pvalue_plots"]][["reactome_plot_over"]]

dev <- pp(file="images/z23_unique_tf.png")
z23_unique_gp[["pvalue_plots"]][["tf_plot_over"]]
closed <- dev.off()
z23_unique_gp[["pvalue_plots"]][["tf_plot_over"]]

z22_z23_shared_gp <- simple_gprofiler(rownames(shared_wanted))

## Performing gProfiler GO search of 103 genes against hsapiens.

## GO search found 92 hits.

## Performing gProfiler KEGG search of 103 genes against hsapiens.

## KEGG search found 15 hits.

## Performing gProfiler REAC search of 103 genes against hsapiens.

## REAC search found 6 hits.

## Performing gProfiler MI search of 103 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 103 genes against hsapiens.

## TF search found 14 hits.

## Performing gProfiler CORUM search of 103 genes against hsapiens.

## CORUM search found 4 hits.

## Performing gProfiler HP search of 103 genes against hsapiens.

## HP search found 0 hits.

dev <- pp(file="images/z23_shared_bp.png", height=9, width=6)
z22_z23_shared_gp[["pvalue_plots"]][["bpp_plot_over"]]
closed <- dev.off()
z22_z23_shared_gp[["pvalue_plots"]][["bpp_plot_over"]]

dev <- pp(file="images/z23_shared_kegg.png")
z22_z23_shared_gp[["pvalue_plots"]][["kegg_plot_over"]]
closed <- dev.off()
z22_z23_shared_gp[["pvalue_plots"]][["kegg_plot_over"]]

dev <- pp(file="images/z23_shared_reac.png")
z22_z23_shared_gp[["pvalue_plots"]][["reactome_plot_over"]]
closed <- dev.off()
z22_z23_shared_gp[["pvalue_plots"]][["reactome_plot_over"]]

dev <- pp(file="images/z23_shared_tf.png")
z22_z23_shared_gp[["pvalue_plots"]][["tf_plot_over"]]
closed <- dev.off()
z22_z23_shared_gp[["pvalue_plots"]][["tf_plot_over"]]

z22_uninf_up_gp <- simple_gprofiler(z22_uninf_up_sig)

## Performing gProfiler GO search of 436 genes against hsapiens.

## GO search found 146 hits.

## Performing gProfiler KEGG search of 436 genes against hsapiens.

## KEGG search found 8 hits.

## Performing gProfiler REAC search of 436 genes against hsapiens.

## REAC search found 18 hits.

## Performing gProfiler MI search of 436 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 436 genes against hsapiens.

## TF search found 63 hits.

## Performing gProfiler CORUM search of 436 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 436 genes against hsapiens.

## HP search found 3 hits.

dev <- pp(file="images/z22_uninf_up_gp_bp.png", height=20, width=10)
z22_uninf_up_gp$pvalue_plots$bpp_plot_over
closed <- dev.off()
z22_uninf_up_gp$pvalue_plots$bpp_plot_over

z22_uninf_down_gp <- simple_gprofiler(z22_uninf_down_sig)

## Performing gProfiler GO search of 6 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 6 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 6 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 6 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 6 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 6 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 6 genes against hsapiens.

## HP search found 0 hits.

z22_sbuninf_up_gp <- simple_gprofiler(z22_sbuninf_up_sig)

## Performing gProfiler GO search of 20 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 20 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 20 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 20 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 20 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 20 genes against hsapiens.

## CORUM search found 4 hits.

## Performing gProfiler HP search of 20 genes against hsapiens.

## HP search found 0 hits.

z22_sbuninf_down_gp <- simple_gprofiler(z22_sbuninf_down_sig)

## Performing gProfiler GO search of 0 genes against hsapiens.

## GO search found 0 hits.

## Performing gProfiler KEGG search of 0 genes against hsapiens.

## KEGG search found 0 hits.

## Performing gProfiler REAC search of 0 genes against hsapiens.

## REAC search found 0 hits.

## Performing gProfiler MI search of 0 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 0 genes against hsapiens.

## TF search found 0 hits.

## Performing gProfiler CORUM search of 0 genes against hsapiens.

## CORUM search found 0 hits.

## Performing gProfiler HP search of 0 genes against hsapiens.

## HP search found 0 hits.

z23_uninf_up_gp <- simple_gprofiler(z23_uninf_up_sig)

## Performing gProfiler GO search of 643 genes against hsapiens.

## GO search found 498 hits.

## Performing gProfiler KEGG search of 643 genes against hsapiens.

## KEGG search found 27 hits.

## Performing gProfiler REAC search of 643 genes against hsapiens.

## REAC search found 46 hits.

## Performing gProfiler MI search of 643 genes against hsapiens.

## MI search found 0 hits.

## Performing gProfiler TF search of 643 genes against hsapiens.

## TF search found 102 hits.

## Performing gProfiler CORUM search of 643 genes against hsapiens.

## CORUM search found 1 hits.

## Performing gProfiler HP search of 643 genes against hsapiens.

## HP search found 7 hits.

dev <- pp(file="images/z23_uninf_up_gp_bp.png", height=20, width=10)
z23_uninf_up_gp$pvalue_plots$bpp_plot_over
closed <- dev.off()
z23_uninf_up_gp$pvalue_plots$bpp_plot_over

6 Parasite reads

sanitize_columns <- c("donor", "drug", "condition", "batch", "macrophagetreatment",
                      "macrophageinfectionzymodeme")
library(org.Lpanamensis.MHOMCOL81L13.v46.eg.db)

## Loading required package: AnnotationDbi

## 

pan_db <- org.Lpanamensis.MHOMCOL81L13.v46.eg.db
all_lp_annot <- sm(load_orgdb_annotations(
    pan_db,
    keytype = "gid",
    fields = c("annot_gene_entrez_id", "annot_gene_name",
               "annot_strand", "annot_chromosome", "annot_cds_length",
               "annot_gene_product")))$genes
macr_lp <- create_expt(metadata="sample_sheets/tmrc2_macrophage_samples_202203.xlsx",
                       gene_info=all_lp_annot, file_column = "lpanamensisv36hisatfile")

## Reading the sample metadata.

## The sample definitions comprises: 28 rows(samples) and 67 columns(metadata fields).

## Warning in file(file, "rt"): cannot open file '/mnt/cbcb/fs01_abelew/
## cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/
## preprocessing/TMRC30051/outputs/hisat2_lpanamensis_v36/
## r1_trimmed.count_lpanamensis_v36_sno_gene_ID.count.xz': No such file or
## directory

## Error in file(file, "rt"): cannot open the connection

  exclude_genes_expt(column = "gene_biotype", method = "keep",
                     pattern = "protein_coding") %>%
  subset_expt(subset="typeofcells=='Macrophages'") %>%
  subset_expt(nonzero = 11000) %>%
  set_expt_factors(columns = sanitize_columns) %>%
  set_expt_factors(columns=sanitize_columns, class="factor") %>%
  set_expt_batches(fact="macrophagetreatment") %>%
  set_expt_conditions(fact="macrophageinfectionzymodeme")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'pData': error in evaluating the argument 'object' in selecting a method for function 'fData': error in evaluating the argument 'object' in selecting a method for function 'fData': error in evaluating the argument 'object' in selecting a method for function 'sampleNames': error in evaluating the argument 'object' in selecting a method for function 'sampleNames': argument "expt" is missing, with no default

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename=savefile))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 64dc311322615b4a8bfef4517ab3c4ba0d0b0e59

## This is hpgltools commit: Mon Apr 4 16:43:54 2022 -0400: 64dc311322615b4a8bfef4517ab3c4ba0d0b0e59

## Saving to tmrc2_macrophage_202203.rda.xz

tmp <- loadme(filename=savefile)