TMRC2 202404: Visualizing Macrophage samples.

1 Introduction

This document will visualize the TMRC2 samples before completing the various differential expression and variant analyses in the hopes of getting an understanding of how the various samples relate to each other.

2 Library sizes and nonzero genes

2.1 Human data

This is really comprised of two datasets, one which is purely the comparison of the various parasite strains, and another which is derived from an experiment performed at CIDEIM which compares human macrophages as well as U937 cells following a infection and/or treatment with the antimonial. With that in mind, the next few plots are of this experiment.

2.1.1 Library sizes

Showing this plot twice in order to make explicit that our range of coverage really is pretty extraordinary.

hs_lib <- plot_libsize(hs_macrophage, yscale = "log2")
pp(file = glue("images/hs_macrophage_libsize-v{ver}.svg"), width = 15, height = 9)
hs_lib$plot
dev.off()

## png 
##   2

hs_lib

## Library sizes of 68 samples, 
## ranging from 15,958,141 to 128,637,602.

hs_lib_log <- plot_libsize(hs_macrophage)
hs_lib_log

## Library sizes of 68 samples, 
## ranging from 15,958,141 to 128,637,602.

Potential start for a figure legend:

Library sizes of the protein coding gene counts observed per sample. The samples were mapped with hisat2 using the hg38 revision 100 human genome; the alignments were sorted, indexed, and counted via htseq using the gene features and non-protein coding features were excluded. The per-sample sums of the remaining matrix were plotted to check that the relative sample coverage is sufficient and not too divergent across samples.

hs_non <- plot_nonzero(hs_macrophage, cutoff = 0.65)

## 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.

pp(file = glue("images/hs_macrophage_nonzero-v{ver}.svg"))
hs_non$plot

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

dev.off()

## png 
##   2

hs_non

## A non-zero genes plot of 68 samples.
## These samples have an average 40.84 CPM coverage and 14870 genes observed, ranging from 13978 to
## 16421.

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

Differences in relative gene content with respect to sequencing coverage. The per-sample number of observed genes was plotted with respect to the relative CPM coverage in order to check that the samples are sufficiently and similarly diverse. Many samples were observed near the putative asymptote of likely gene content; no samples were observed with fewer than 65% of the human genes included.

hs_box <- plot_boxplot(hs_macrophage)

## 449491 entries are 0.  We are on a log scale, adding 1 to the data.

pp(file = "images/hs_macrophage_boxplot.svg")
hs_box
dev.off()

## png 
##   2

hs_box

The distribution of observed counts / gene for all samples was plotted as a boxplot on the log2 scale. No genes were observed as explicit outliers and the range of mean coverage spanned an order of magnitude from 20-200 reads/gene. Quartile boxes were colored according to infection status and drug treatment.

filter_plot <- plot_libsize_prepost(hs_macrophage)
pp(file = "images/hs_macrophage_lowgene.svg")
filter_plot$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

dev.off()

## png 
##   2

filter_plot

## A comparison of the counts before and after filtering.
## The number of genes with low coverage changes by 5999-8102 genes.

## Warning: Using alpha for a discrete variable is not advised.

pp(file = "images/hs_macrophage_lowcount.svg")
filter_plot$count_plot
dev.off()

## png 
##   2

filter_plot$count_plot

Numbers of low-count genes before and after filtering. The height of each bar represents the number of low-count genes (>= 2 counts) before performing a low-count filter. The lower bar represents the number of genes remaining after low-count filtering, the number inside the bar is the difference.

When the low-count filter is applied, some samples have significantly more than the cutoff number of counts/gene (2). As a result, when the number of total counts removed is plotted, that sum is ubiquitously more than twice the number of removed gene and only approaches that threshold for the samples with lowest coverage. This suggests that a coefficient of variance-based filter may be more appropriate under some circumstances.

2.2 Distribution Visualizations

The distribution of samples observed in the macrophage dataset is evocative and suggests that there are very clear distinctions between the two strains as well as the drug treatment.

2.2.1 PCA

There are a few ways we can express the sample distribution via PCA; in this instance we explicitly concatenate the infection and drug treatment status into one factor.

As of 202212, we now have a set of samples from two different cell types: macrophages and U937. Thus the first order of business is to observe their similarities/differences.

type_batch <- set_expt_batches(hs_macrophage, fact="typeofcells")

## The number of samples by batch are:

## 
## Macrophages        U937 
##          54          14

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

## Removing 9198 low-count genes (12283 remaining).

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

type_batch_pca <- plot_pca(type_batch_norm, plot_title = "PCA of macrophage expression values",
                           plot_alpha = 0.6, plot_labels = FALSE)

pp(file = glue("images/type_batch_macrophage_norm_pca-v{ver}.svg"))
type_batch_pca$plot
dev.off()

## png 
##   2

type_batch_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf, inf_sb, uninf, uninf_sb
## Shapes are defined by Macrophages, U937.

The differences, as expected are stark. Thus, let us consider the macrophage and U937 samples separately. In addition, consider the samples behind the lense of sample date, strain, infection status, etc.

hs_macr <- subset_expt(hs_macrophage, subset = "typeofcells=='Macrophages'")

## The samples excluded are: TMRC30309, TMRC30293, TMRC30294, TMRC30291, TMRC30292, TMRC30307, TMRC30308, TMRC30310, TMRC30331, TMRC30311, TMRC30332, TMRC30305, TMRC30306, TMRC30330.

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

written <- write_expt(
  hs_macr,
  excel = glue("analyses/macrophage_de/{ver}/read_counts/hs_macr-v{ver}.xlsx"))

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

## 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.
## 353608 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Changed 353608 zero count features.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Initial model failed:
## The variables specified in this model are redundant,
## so the design matrix is not full rank

## `geom_smooth()` using formula = 'y ~ x'
## Subsetting on features.
## 
## remove_genes_expt(), before removal, there were 11756 genes, now there are 11756.

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Initial model failed:
## The variables specified in this model are redundant,
## so the design matrix is not full rank

## Retrying with only condition in the model.

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

## Removing 9725 low-count genes (11756 remaining).
## transform_counts: Found 5 values equal to 0, adding 1 to the matrix.

macr_pca <- plot_pca(macr_norm, plot_title = "PCA of macrophage expression values",
                     plot_alpha = 0.6, plot_labels = FALSE)

pp(file = glue("images/macr_norm_pca_labeled-v{ver}.png"), width=20, height=12)
macr_pca$plot +
  geom_text_repel(aes(label=sampleid), force=1, label.padding=10, box.padding=10, max.overlaps=100)

## Warning in geom_text_repel(aes(label = sampleid), force = 1, label.padding =
## 10, : Ignoring unknown parameters: `label.padding`

dev.off()

## png 
##   2

pp(file = glue("images/macr_norm_pca-v{ver}.svg"))
macr_pca$plot
dev.off()

## png 
##   2

macr_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf, inf_sb, uninf, uninf_sb
## Shapes are defined by none, z2.2, z2.3.

hs_macr_datebatch <- set_expt_batches(hs_macr, fact = "oldnew")

## The number of samples by batch are:

## 
##  current previous 
##       27       27

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

## Removing 9725 low-count genes (11756 remaining).
## transform_counts: Found 5 values equal to 0, adding 1 to the matrix.

macr_date_pca <- plot_pca(
  macr_date_norm,
  plot_title = "PCA of macrophage expression values, date batch",
  plot_alpha = 0.6, plot_labels = FALSE)
pp(file = glue("images/macr_norm_date_pca-v{ver}.svg"))
macr_date_pca$plot
dev.off()

## png 
##   2

macr_date_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf, inf_sb, uninf, uninf_sb
## Shapes are defined by current, previous.

Some of the same questions/observations are interesting for the U937 samples.

Note to self, I already created this subset in the tmrc2_datasets file, this is redundant.

hs_u937 <- subset_expt(hs_macrophage, subset = "typeofcells!='Macrophages'")

## The samples excluded are: TMRC30051, TMRC30057, TMRC30059, TMRC30060, TMRC30061, TMRC30062, TMRC30063, TMRC30064, TMRC30065, TMRC30066, TMRC30067, TMRC30069, TMRC30117, TMRC30243, TMRC30244, TMRC30245, TMRC30246, TMRC30247, TMRC30248, TMRC30249, TMRC30250, TMRC30251, TMRC30252, TMRC30266, TMRC30267, TMRC30268, TMRC30286, TMRC30326, TMRC30316, TMRC30317, TMRC30322, TMRC30323, TMRC30328, TMRC30318, TMRC30319, TMRC30324, TMRC30325, TMRC30320, TMRC30321, TMRC30327, TMRC30312, TMRC30297, TMRC30298, TMRC30299, TMRC30300, TMRC30295, TMRC30296, TMRC30303, TMRC30304, TMRC30301, TMRC30302, TMRC30314, TMRC30315, TMRC30313.

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

u937_written <- write_expt(
  hs_u937, excel = glue("analyses/macrophage_de/{ver}/read_counts/hs_u937_reads-v{ver}.xlsx"))

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

## 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.
## 95883 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Changed 95883 zero count features.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Initial model failed:
## The variables specified in this model are redundant,
## so the design matrix is not full rank

## `geom_smooth()` using formula = 'y ~ x'
## Subsetting on features.
## 
## remove_genes_expt(), before removal, there were 10751 genes, now there are 10751.

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Initial model failed:
## The variables specified in this model are redundant,
## so the design matrix is not full rank

## Retrying with only condition in the model.

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

## Removing 10730 low-count genes (10751 remaining).

u937_pca <- plot_pca(u937_norm, plot_title = "PCA of U937 expression values",
                     plot_alpha = 0.6, plot_labels = FALSE)
pp(file = glue("images/u937_norm_pca-v{ver}.svg"))
u937_pca$plot
dev.off()

## png 
##   2

u937_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf, inf_sb, uninf, uninf_sb
## Shapes are defined by none, z2.2, z2.3.

Having taken some views of the samples, repeat with sva.

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

## Removing 9725 low-count genes (11756 remaining).

## Setting 1746 low elements to zero.

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

hs_nb_pca <- plot_pca(hs_nb, plot_title = "PCA of macrophage expression values post-sva",
                      plot_alpha = 0.6, plot_labels = FALSE)

pp(file = "images/hs_macr_nb_pca.svg")
hs_nb_pca$plot
dev.off()

## png 
##   2

hs_nb_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf, inf_sb, uninf, uninf_sb
## Shapes are defined by none, z2.2, z2.3.

Some likely text for a figure legend might include something like the following (paraphrased from Najib’s 2016 dual transcriptome profiling paper (10.1128/mBio.00027-16)):

Expression profiles of the human macrophages in response to drug treatment and/or parasite infection by two strains. Each glyph represents one sample, the circular samples were not infected; green circles were not treated while yellow circles were. Zymodeme 2.2 infected samples are squares with/out antimonial (purple and green respectively). Zymodeme 2.3 infected samples are diamonds with/out antimonial treatment (pink and orange respectively). This analysis was performed following a low-count filter, cpm conversion, quantile normalization, and a log2 transformation. The second plot is identical except surrogate estimates were derived by svaseq (which estimated 3 surrogates) after normalization and those estimates were used to modify the normalized counts before log2 transformation.

Some interpretation for this figure might include:

When PCA was performed on the human macrophage transcriptome data, the first principal component coincided with drug treatment and the second with differences between zymodemes. The uninfected samples were observed to have very similar transcriptional profiles to the corresponding zymodeme 2.2 samples, suggesting that infection with zymodeme 2.2 has a relatively minimal impact on the macrophage. In contrast, the zymodeme 2.3 samples, excepting one drug treated sample, displayed a significant shift from the uninfected and z2.2. This shift was particularly evident in the non-treated samples. When sva was applied to this data, it slightly improved the observed differences between the drug treatments and zymodeme infections, but did not significantly change the relationship between z2.2 and the uninfected samples.

2.3 Compare to drug + infection/zymodeme

The primary contrasts actually examined in 03differential_expression.Rmd are after combining the metadata factors associated with the drug treatment and infection strain (or uninfected). With that in mind, I should make a version of the above using the same logic.

Thus, for the moment I am going to copy the appropriate block straight out of the 03 file.

new_conditions <- paste0(pData(hs_macrophage)[["macrophagetreatment"]], "_",
  pData(hs_macrophage)[["macrophagezymodeme"]])
## Note the sanitize() call is redundant with the addition of sanitize() in the
## datastructures file, but I don't want to wait to rerun that.
hs_macr <- set_expt_conditions(hs_macrophage, fact = new_conditions) %>%
  sanitize_expt_pData(column = "drug") %>%
  subset_expt(subset = "typeofcells!='U937'")

## The numbers of samples by condition are:

## 
##      inf_z22      inf_z23    infsb_z22    infsb_z23   uninf_none uninfsb_none 
##           14           15           15           14            5            5

## The samples excluded are: TMRC30309, TMRC30293, TMRC30294, TMRC30291, TMRC30292, TMRC30307, TMRC30308, TMRC30310, TMRC30331, TMRC30311, TMRC30332, TMRC30305, TMRC30306, TMRC30330.

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

I am reasonably certain this is the dataset of interest.

2.3.1 Recreate PCA with labels

So, now let us repeat the pca with the extra labels and larger canvas.

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

## Removing 9725 low-count genes (11756 remaining).

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

de_pca <- plot_pca(hs_macr_norm)
## ok, that looks just like the previous image, now lets make it bigger with labels...

pp(file = "images/infection_drug_pca.png", width = 20, height = 12)
de_pca$plot +
  geom_text_repel(aes(label = sampleid), force = 1, label.padding = 10,
                  box.padding = 10, max.overlaps = 100)

## Warning in geom_text_repel(aes(label = sampleid), force = 1, label.padding =
## 10, : Ignoring unknown parameters: `label.padding`

dev.off()

## png 
##   2

de_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by inf_z22, inf_z23, infsb_z22, infsb_z23, uninf_none, uninfsb_none
## Shapes are defined by none, z2.2, z2.3.

2.3.2 Correlation heatmaps

corheat <- plot_corheat(macr_norm, plot_title = "Correlation heatmap of macrophage
                 expression values
")
corheat

## A heatmap of pairwise sample correlations ranging from: 
## 0.874321963018488 to 0.99652359119492.

corheat <- plot_corheat(hs_nb, plot_title = "Correlation heatmap of macrophage
                 expression values
")
corheat

## A heatmap of pairwise sample correlations ranging from: 
## 0.928918488834237 to 0.996693097072659.

disheat <- plot_disheat(macr_norm, plot_title = "Euclidean heatmap of macrophage
                 expression values
")
disheat

## A heatmap of pairwise sample distances ranging from: 
## 19.7826122756158 to 118.945131079151.

disheat_nb <- plot_disheat(hs_nb, plot_title = "Euclidean heatmap of macrophage
                 expression values
")
disheat_nb

## A heatmap of pairwise sample distances ranging from: 
## 19.0138353587251 to 90.5582333469978.

plot_sm(macr_norm)$plot

Potential start for a figure legend:

Global relationships among the human macrophage transcriptional profiles. Pairwise pearson correlations and Euclidean distances were calculated using the normalized expression matrices. Colors along the top row delineate the experimental conditions (same colors as the PCA) and the colors along the first column delineate infection status (purple: z2.3, yellow: z2.2, green: uninfected). Samples were clustered by nearest neighbor clustering and each colored tile describes one correlation value between two samples (red to white delineates pearson correlation values of the 11,460 normalized gene values between two samples ranging from <= 0.9 to >= 0.98) or the euclidean distance between two samples (dark blue to white delineates identical to a normalized euclidean distance of >= 90).

Some interpretation for this figure might include:

When the global relationships among the samples were distilled down to individual euclidean distances or pearson correlation coefficients between pairs of samples, the primary clustering among samples observed was according to drug treatment. Secondary clades intermingled the z2.2 and uninfected samples. The data before svaseq provides weak evidence for the hypothesis that sample TMRC30062 (z2.3 drug treated) is actually a z2.2 drug treated sample. This hypothesis was discounted by manually examining the relatively few parasite reads in IGV and comparing the observed variant positions to other known (not drug treated) z2.3 and z2.2 samples (in this case I compared TMRC30062, the sample in question, to TMRC30061 (the same strain without drug treatment), TMRC30063 (another z2.3 strain without drug treatment) and TMRC30286 (a z2.2 (strain ID 11075) sample which was not treated); because the drug treated sample has few reads, it is difficult to find z2.2-specific variants; but positive matches were readily identifiable to previously characterized z2.3-specific variants. Three such locations are shown in the following image (chromosome 31, 682.4Kb, 691.5kb, 673kb):

• Note, I commented out the picture because I am worried it will mess up rendering in the container.

I also wrote a scoring function which sums up all observed variant positions by putative zymodeme status and for this sample it found 11 positions which were z2.2 specific, 119 which were z2.3 specific, 2808 positions (out of 3541: 0.793) which did not have z2.2 specific variants and 71,944 positions (out of 81,556: 0.882) which did not have z2.3 specific variants). The proportions of strain specific (un)observed variant positions is interesting because it changed over time/celltype for some people. The observed changes might just be noise in the data (we only observed 130 specific positions out of ~80,000 in this sample, for example), but in at least some cases it seems evocative.

3 Sources of variance

One thing I did not realize until this very moment, weird combinations of donor/drug/treatment/zymodeme are confounded. Drug and treatment are confounded by definition, but it looks like donor and treatment are, too. This limits the ability to run variancePartition slightly.

3.1 Donor, drug, zymodeme

Let us first run variancePartition in the order donor,drug,zymodeme; which, as Theresa noted, should bias the results in that order.

Oddly, when I did this, it appears to not change the results at all, let us look a little more closely.

table(pData(hs_macr)[["donor"]])

## 
## d01 d02 d09 d81 
##  13  14  13  14

table(pData(hs_macr)[["drug"]])

## 
## antimony     none 
##       27       27

table(pData(hs_macr)[["macrophagezymodeme"]])

## 
## none  z22  z23 
##    8   23   23

table(pData(hs_macr)[["rnaextractiondate"]])

## 
## 20200113 20210831 20211118 20211221 20220829 20220919 
##       13       11        1        2       14       13

table(pData(hs_macr)[["oldnew"]])

## 
##  current previous 
##       27       27

donor_drug_zymo_varpart <- simple_varpart(
    hs_macr,
    factors = c("donor", "drug", "macrophagezymodeme"))

## Subsetting on features.

## remove_genes_expt(), before removal, there were 17804 genes, now there are 15325.

donor_drug_zymo_varpart

## The result of using variancePartition with the model:
## ~ donor + drug + macrophagezymodeme

by_donor <- replot_varpart_percent(donor_drug_zymo_varpart,
                                   column = "donor")
by_donor

by_zymo <- replot_varpart_percent(donor_drug_zymo_varpart,
                                  column = "macrophagezymodeme")
by_zymo

top_ids <- gsub(x = head(rownames(by_zymo[["resorted"]]), n = 100),
                pattern = "^gene:", replacement = "")
high_zymo_variance_gp <- simple_gprofiler(top_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

## Add a little logic here to use enrichplot::dotplot().

high_zymo_variance_gp

## A set of ontologies produced by gprofiler using 100
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are 122 GO hits, 6, KEGG hits, 12 reactome hits, 7 wikipathway hits, 80 transcription factor hits, 1 miRNA hits, 0 HPA hits, 0 HP hits, and 0 CORUM hits.

## Category BP is the most populated with 30 hits.

written <- write_gprofiler_data(
  high_zymo_variance_gp,
  excel = glue("analyses/macrophage_de/{ver}/gprofiler/high_zymodeme_variance-v{ver}.xlsx"))

## Loading required namespace: Vennerable

by_drug <- replot_varpart_percent(donor_drug_zymo_varpart, column = "drug")
by_drug

top_ids <- gsub(x = head(rownames(by_drug[["resorted"]]), n = 100),
                pattern = "^gene:", replacement = "")
high_drug_variance_gp <- simple_gprofiler(top_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
## No results to show
## Please make sure that the organism is correct or set significant = FALSE

## Add a little logic here to use enrichplot::dotplot().

high_drug_variance_gp

## A set of ontologies produced by gprofiler using 100
## genes against the hsapiens annotations and significance cutoff 0.05.
## There are 22 GO hits, 0, KEGG hits, 4 reactome hits, 1 wikipathway hits, 24 transcription factor hits, 0 miRNA hits, 0 HPA hits, 0 HP hits, and 0 CORUM hits.

## Category TF is the most populated with 24 hits.

high_drug_variance_gp$pvalue_plots$MF

written <- write_gprofiler_data(
  high_drug_variance_gp,
  excel = glue("analyses/macrophage_de/{ver}/gprofiler/high_drug_variance-v{ver}.xlsx"))

time_drug_zymo_varpart <- simple_varpart(
    hs_macr,
    factors = c("oldnew", "drug", "macrophagezymodeme"))

## Subsetting on features.

## remove_genes_expt(), before removal, there were 17804 genes, now there are 15325.

time_drug_zymo_varpart$partition_plot

3.1.1 Compare to PCA with drug as color and donor as shape.

drug_donor_test <- set_expt_conditions(hs_macr, fact = "drug") %>%
  set_expt_batches(fact = "donor")

## The numbers of samples by condition are:

## 
## antimony     none 
##       27       27

## The number of samples by batch are:

## 
## d01 d02 d09 d81 
##  13  14  13  14

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

## Removing 9725 low-count genes (11756 remaining).

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

plot_pca(drug_donor_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by antimony, none
## Shapes are defined by d01, d02, d09, d81.

drug_norm_cv <- plot_variance_coefficients(drug_donor_norm)

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

## Finished calculating dispersion estimates.

drug_norm_cv

## Plot describing the observed variance coefficients on a per-gene basis.

3.1.2 Same comparison, but this time with infected state as shape

drug_infectedp_test <- set_expt_conditions(hs_macr, fact = "drug") %>%
  set_expt_batches(fact = "batch")

## The numbers of samples by condition are:

## 
## antimony     none 
##       27       27

## The number of samples by batch are:

## 
## none z2.2 z2.3 
##    8   23   23

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

## Removing 9725 low-count genes (11756 remaining).

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

plot_pca(drug_infectedp_norm)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by antimony, none
## Shapes are defined by none, z2.2, z2.3.

3.2 Drug zymodeme donor

drug_zymo_donor_varpart <- simple_varpart(
    hs_macr,
    factors = c("drug", "macrophagezymodeme", "donor"))

## Subsetting on features.

## remove_genes_expt(), before removal, there were 17804 genes, now there are 15325.

donor_drug_zymo_varpart$partition_plot

donor_drug_zymo_varpart$model_used

## ~donor + drug + macrophagezymodeme
## <environment: 0x557e7f6b8140>

drug_zymo_donor_varpart$partition_plot

drug_zymo_donor_varpart$model_used

## ~drug + macrophagezymodeme + donor
## <environment: 0x557e971eba08>

I think the above shows that the order of elements used in the model do not affect the results.

4 Examine Donors

4.1 Donor comparison

combined_batch <- paste0(pData(hs_macr)[["drug"]], "_",
  pData(hs_macr)[["macrophagezymodeme"]])
hs_donors <- set_expt_conditions(hs_macr, fact = "donor") %>%
  set_expt_batches(fact = "macrophagezymodeme")

## The numbers of samples by condition are:

## 
## d01 d02 d09 d81 
##  13  14  13  14

## The number of samples by batch are:

## 
## none  z22  z23 
##    8   23   23

donor_norm <- normalize_expt(hs_donors, filter = TRUE,
                             convert = "cpm", transform = "log2")

## Removing 9725 low-count genes (11756 remaining).

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

donor_norm_cv <- plot_variance_coefficients(donor_norm)

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

## Finished calculating dispersion estimates.

donor_norm_cv

## Plot describing the observed variance coefficients on a per-gene basis.

pp(file = "images/tmrc2_macrophages_donors.svg")
plot_pca(donor_norm, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by d01, d02, d09, d81
## Shapes are defined by none, z22, z23.

dev.off()

## png 
##   2

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

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

## Removing 9725 low-count genes (11756 remaining).

## Setting 1456 low elements to zero.

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

pp(file = "images/tmrc2_macrophages_donors_sva.png")
plot_pca(donor_nb, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by d01, d02, d09, d81
## Shapes are defined by none, z22, z23.

dev.off()

## png 
##   2

4.2 Drug comparison

hs_drug <- set_expt_conditions(hs_macr, fact = "drug")

## The numbers of samples by condition are:

## 
## antimony     none 
##       27       27

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

## Removing 9725 low-count genes (11756 remaining).

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

pp(file = "images/tmrc2_macrophages_drugs.svg")
plot_pca(drug_norm, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by antimony, none
## Shapes are defined by none, z2.2, z2.3.

dev.off()

## png 
##   2

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

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

## Removing 9725 low-count genes (11756 remaining).

## Setting 1120 low elements to zero.

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

pp(file = "images/tmrc2_macrophages_drugs_sva.svg")
plot_pca(drug_nb, plot_labels = FALSE)

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by antimony, none
## Shapes are defined by none, z2.2, z2.3.

dev.off()

## png 
##   2

4.3 Parasite data

4.3.1 Library sizes, non-zero plot, and distribution

lp_lib <- plot_libsize(lp_macrophage)
pp(file = "images/lp_macrophage_libsize.svg")
lp_lib$plot
dev.off()

## png 
##   2

lp_lib

## Library sizes of 50 samples, 
## ranging from 5,742 to 7,646,644.

lp_non <- plot_nonzero(lp_macrophage)

## The following samples have less than 5661.5 genes.

## [1] "TMRC30249" "TMRC30300" "TMRC30302" "TMRC30292" "TMRC30331" "TMRC30332"

## 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.

pp(file = "images/lp_macrophage_nonzero.svg")
lp_non$plot

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

dev.off()

## png 
##   2

lp_non

## A non-zero genes plot of 50 samples.
## These samples have an average 1.266 CPM coverage and 7818 genes observed, ranging from 2849 to
## 8596.

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

lp_box <- plot_boxplot(lp_macrophage)

## 44583 entries are 0.  We are on a log scale, adding 1 to the data.

pp(file = "images/lp_macrophage_boxplot.svg")
lp_box
dev.off()

## png 
##   2

lp_box

filter_plot <- plot_libsize_prepost(lp_macrophage)
pp(file = "images/lp_macrophage_lowgene.svg")
filter_plot$lowgene_plot

## Warning: Using alpha for a discrete variable is not advised.

dev.off()

## png 
##   2

filter_plot

## Warning in min(changed): no non-missing arguments to min; returning Inf

## Warning in max(changed): no non-missing arguments to max; returning -Inf

## A comparison of the counts before and after filtering.
## The number of genes with low coverage changes by Inf--Inf genes.

## Warning: Using alpha for a discrete variable is not advised.

pp(file = "images/lp_macrophage_lowcount.svg")
filter_plot$count_plot
dev.off()

## png 
##   2

filter_plot$count_plot

The parasite metrics are identical in theory to the human macrophage plots above with one relevant difference. In the box plot, the distribution of observed reads/gene follows a different distribution than what was observed in the host, this difference is due to the very different transcriptional profile of the Leishmania parasite.

4.4 Distribution Visualizations

The PCA and heatmap plots for the parasite samples are largely identical in concept to the macrophage plots above with one very important difference. Only 1 of the drug treated samples has sufficient parasite reads remaining to effectively quantify it, the other parasite samples were removed.

This lone post-drug treated samples (TMRC30248, strain 11026) had the largest parasite load before treatment by a tremendous margin (it has 268,826 SL reads compared to 72,489 in the second highest sample, the mean of the pre-drug treated samples is approximately 69,137 and median is 48,090). Thus it is perhaps not surprising that it still has a significant number of SL-containing reads following treatment (30,052 vs 14,418 in the second highest).

4.4.1 Keeping one antimonial treated sample

The following plots show the distinction between the two strains used in the experiment very clearly and suggest that, in the one case with sufficient surviving post-treatment parasites, the parasite transcriptional profile was not significantly changed by the antimonial treatment.

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

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

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

lp_pca <- plot_pca(lp_norm, plot_title = "PCA of macrophage expression values",
                   plot_labels = FALSE)
pp(file = "images/lp_macrophage_norm_pca.svg")
lp_pca$plot
dev.off()

## png 
##   2

lp_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by z2.2, z2.3
## Shapes are defined by inf, inf_sb.

lp_nb <- normalize_expt(lp_macrophage, convert = "cpm", transform = "log2",
                        filter = "simple", batch = "svaseq")

## Removing 60 low-count genes (8650 remaining).

## Setting 6107 low elements to zero.

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

lp_nb_pca <- plot_pca(lp_nb, plot_title = "PCA of macrophage expression values post-sva",
                      plot_labels = FALSE)
pp(file = "images/lp_macrophage_nb_pca.svg")
lp_nb_pca$plot

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

dev.off()

## png 
##   2

lp_nb_pca

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by z2.2, z2.3
## Shapes are defined by inf, inf_sb.

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

corheat <- plot_corheat(lp_norm, plot_title = "Correlation heatmap of parasite
                 expression values
")
corheat

## A heatmap of pairwise sample correlations ranging from: 
## 0.258937002313872 to 0.988998797081892.

corheat <- plot_corheat(lp_nb, plot_title = "Correlation heatmap of parasite
                 expression values
")
corheat

## A heatmap of pairwise sample correlations ranging from: 
## 0.621226444997473 to 0.999663911744281.

plot_sm(lp_norm)

## When the standard median metric was plotted, the values observed range
## from 0.258937002313872 to 1 with quartiles at 0.799469037186123 and 0.923961447524501.

The following repeats the parasite PCA without the peculiar post-antimonial sample.

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

## Removing 119 low-count genes (8591 remaining).

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

lp_pca_nosb <- plot_pca(lp_norm_nosb, plot_title = "PCA of macrophage expression values",
                   plot_labels = FALSE)
pp(file = "images/lp_macrophage_norm_nosb_pca.svg")
lp_pca_nosb$plot
dev.off()

## png 
##   2

lp_pca_nosb

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by z2.2, z2.3
## Shapes are defined by inf.

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

## Removing 119 low-count genes (8591 remaining).

## Setting 889 low elements to zero.

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

lp_nb_pca_nosb <- plot_pca(lp_nb_nosb, plot_title = "PCA of macrophage expression values post-sva",
                      plot_labels = FALSE)
pp(file = "images/lp_macrophage_nb_nosb_pca.svg")
lp_nb_pca_nosb$plot
dev.off()

## png 
##   2

lp_nb_pca_nosb

## The result of performing a fast_svd dimension reduction.
## The x-axis is PC1 and the y-axis is PC2
## Colors are defined by z2.2, z2.3
## Shapes are defined by inf.

corheat_nosb <- plot_corheat(lp_norm_nosb, plot_title = "Correlation heatmap of parasite
                 expression values
")
corheat_nosb

## A heatmap of pairwise sample correlations ranging from: 
## 0.758351667603457 to 0.986554785481698.

corheat_nosb <- plot_corheat(lp_nb_nosb, plot_title = "Correlation heatmap of parasite
                 expression values
")
corheat_nosb

## A heatmap of pairwise sample correlations ranging from: 
## 0.894981656374377 to 0.99960033240918.

5 Some genes of interest

One query from Olga and Maria Colmenares is to query a group of kinetoplast genes, including these:

• kkt6: LPAL13_120005800
• kkt2
• KPAF1
• KPAF1
• kkt3
• KRIPP3
• kkt1

Checking my annotations, they appear to be coming from the column ‘annot_gene_name’.

wanted <- c("kkt6", "kkt2", "KPAF1", "kkt3", "KRIPP3", "kkt1")
wanted %in% fData(lp_macrophage)[["annot_gene_name"]]

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

ids <- fData(lp_macrophage)[["annot_gene_name"]] %in% wanted
ids <- rownames(exprs(lp_macrophage))[ids]
names <- fData(lp_macrophage)[ids, "annot_gene_name"]
lp_norm <- normalize_expt(lp_macrophage, transform = "log2", convert = "cpm", norm = "quant")

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

few_lp <- subset_genes(lp_norm, ids = ids, method = "keep")

## remove_genes_expt(), before removal, there were 8710 genes, now there are 7.

## There are 50 samples which kept less than 90 percent counts.

## TMRC30051 TMRC30057 TMRC30061 TMRC30062 TMRC30063 TMRC30064 TMRC30065 TMRC30067 
##   0.09247   0.09235   0.09067   0.08592   0.09294   0.09208   0.08236   0.09334 
## TMRC30069 TMRC30162 TMRC30245 TMRC30247 TMRC30248 TMRC30249 TMRC30250 TMRC30251 
##   0.09283   0.09385   0.09271   0.09186   0.08981   0.08670   0.09441   0.07979 
## TMRC30252 TMRC30267 TMRC30286 TMRC30316 TMRC30317 TMRC30322 TMRC30328 TMRC30318 
##   0.09265   0.09283   0.09108   0.09335   0.09380   0.09212   0.09191   0.09047 
## TMRC30324 TMRC30320 TMRC30321 TMRC30297 TMRC30298 TMRC30299 TMRC30300 TMRC30295 
##   0.09389   0.09216   0.09073   0.09095   0.08958   0.09046   0.08788   0.09320 
## TMRC30296 TMRC30303 TMRC30301 TMRC30302 TMRC30314 TMRC30315 TMRC30293 TMRC30294 
##   0.09319   0.09386   0.09343   0.08333   0.09133   0.09467   0.09355   0.08890 
## TMRC30291 TMRC30292 TMRC30307 TMRC30308 TMRC30310 TMRC30331 TMRC30311 TMRC30332 
##   0.08992   0.08890   0.09079   0.08501   0.09198   0.07816   0.08741   0.08563 
## TMRC30305 TMRC30306 
##   0.09044   0.09227

few_heat <- plot_sample_heatmap(few_lp, row_label = names)
few_heat

pander::pander(sessionInfo())

R version 4.3.3 (2024-02-29)

Platform: x86_64-conda-linux-gnu (64-bit)

locale: C

attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ruv(v.0.9.7.1), BiocParallel(v.1.36.0), variancePartition(v.1.32.5), hpgltools(v.1.0), testthat(v.3.2.1), Matrix(v.1.6-5), glue(v.1.7.0), SummarizedExperiment(v.1.32.0), GenomicRanges(v.1.54.1), GenomeInfoDb(v.1.38.8), IRanges(v.2.36.0), S4Vectors(v.0.40.2), MatrixGenerics(v.1.14.0), matrixStats(v.1.2.0), Biobase(v.2.62.0), BiocGenerics(v.0.48.1), Heatplus(v.3.10.0), ggrepel(v.0.9.5) and ggplot2(v.3.5.0)

loaded via a namespace (and not attached): fs(v.1.6.3), bitops(v.1.0-7), enrichplot(v.1.22.0), devtools(v.2.4.5), doParallel(v.1.0.17), HDO.db(v.0.99.1), httr(v.1.4.7), RColorBrewer(v.1.1-3), numDeriv(v.2016.8-1.1), profvis(v.0.3.8), tools(v.4.3.3), backports(v.1.4.1), utf8(v.1.2.4), R6(v.2.5.1), lazyeval(v.0.2.2), mgcv(v.1.9-1), urlchecker(v.1.0.1), withr(v.3.0.0), prettyunits(v.1.2.0), gridExtra(v.2.3), preprocessCore(v.1.64.0), cli(v.3.6.2), scatterpie(v.0.2.2), labeling(v.0.4.3), sass(v.0.4.9), mvtnorm(v.1.2-4), genefilter(v.1.84.0), Rsamtools(v.2.18.0), yulab.utils(v.0.1.4), gson(v.0.1.0), DOSE(v.3.28.2), sessioninfo(v.1.2.2), limma(v.3.58.1), rstudioapi(v.0.16.0), RSQLite(v.2.3.6), generics(v.0.1.3), gridGraphics(v.0.5-1), BiocIO(v.1.12.0), crosstalk(v.1.2.1), gtools(v.3.9.5), zip(v.2.3.1), dplyr(v.1.1.4), GO.db(v.3.18.0), fansi(v.1.0.6), abind(v.1.4-5), lifecycle(v.1.0.4), yaml(v.2.3.8), edgeR(v.4.0.16), Rtsne(v.0.17), gplots(v.3.1.3.1), qvalue(v.2.34.0), SparseArray(v.1.2.4), BiocFileCache(v.2.10.2), grid(v.4.3.3), blob(v.1.2.4), promises(v.1.3.0), crayon(v.1.5.2), miniUI(v.0.1.1.1), lattice(v.0.22-6), cowplot(v.1.1.3), GenomicFeatures(v.1.54.4), annotate(v.1.80.0), KEGGREST(v.1.42.0), pillar(v.1.9.0), knitr(v.1.46), varhandle(v.2.0.6), fgsea(v.1.28.0), rjson(v.0.2.21), boot(v.1.3-30), corpcor(v.1.6.10), codetools(v.0.2-20), fastmatch(v.1.1-4), ggfun(v.0.1.4), Vennerable(v.3.1.0.9000), data.table(v.1.15.4), remotes(v.2.5.0), treeio(v.1.26.0), vctrs(v.0.6.5), png(v.0.1-8), Rdpack(v.2.6), gtable(v.0.3.4), cachem(v.1.0.8), openxlsx(v.4.2.5.2), xfun(v.0.43), rbibutils(v.2.2.16), S4Arrays(v.1.2.1), mime(v.0.12), tidygraph(v.1.3.1), survival(v.3.5-8), iterators(v.1.0.14), statmod(v.1.5.0), directlabels(v.2024.1.21), ellipsis(v.0.3.2), nlme(v.3.1-164), pbkrtest(v.0.5.2), ggtree(v.3.11.1.001), usethis(v.2.2.3), bit64(v.4.0.5), progress(v.1.2.3), EnvStats(v.2.8.1), filelock(v.1.0.3), rprojroot(v.2.0.4), bslib(v.0.7.0), KernSmooth(v.2.23-22), colorspace(v.2.1-0), DBI(v.1.2.2), tidyselect(v.1.2.1), bit(v.4.0.5), compiler(v.4.3.3), curl(v.5.2.1), graph(v.1.80.0), xml2(v.1.3.6), desc(v.1.4.3), DelayedArray(v.0.28.0), plotly(v.4.10.4), shadowtext(v.0.1.3), rtracklayer(v.1.62.0), scales(v.1.3.0), caTools(v.1.18.2), remaCor(v.0.0.18), quadprog(v.1.5-8), RBGL(v.1.78.0), rappdirs(v.0.3.3), stringr(v.1.5.1), digest(v.0.6.35), minqa(v.1.2.6), rmarkdown(v.2.26), aod(v.1.3.3), XVector(v.0.42.0), RhpcBLASctl(v.0.23-42), htmltools(v.0.5.8.1), pkgconfig(v.2.0.3), lme4(v.1.1-35.2), highr(v.0.10), dbplyr(v.2.3.4), fastmap(v.1.1.1), rlang(v.1.1.3), htmlwidgets(v.1.6.4), shiny(v.1.8.1.1), farver(v.2.1.1), jquerylib(v.0.1.4), jsonlite(v.1.8.8), GOSemSim(v.2.28.1), RCurl(v.1.98-1.14), magrittr(v.2.0.3), GenomeInfoDbData(v.1.2.11), ggplotify(v.0.1.2), patchwork(v.1.2.0), munsell(v.0.5.1), Rcpp(v.1.0.12), ape(v.5.7-1), viridis(v.0.6.5), stringi(v.1.8.3), ggraph(v.2.2.1), brio(v.1.1.4), zlibbioc(v.1.48.2), MASS(v.7.3-60.0.1), plyr(v.1.8.9), pkgbuild(v.1.4.4), parallel(v.4.3.3), forcats(v.1.0.0), Biostrings(v.2.70.3), graphlayouts(v.1.1.1), splines(v.4.3.3), pander(v.0.6.5), hms(v.1.1.3), locfit(v.1.5-9.9), fastcluster(v.1.2.6), igraph(v.2.0.3), reshape2(v.1.4.4), biomaRt(v.2.58.2), pkgload(v.1.3.4), gprofiler2(v.0.2.3), XML(v.3.99-0.16.1), evaluate(v.0.23), nloptr(v.2.0.3), PROPER(v.1.34.0), foreach(v.1.5.2), tweenr(v.2.0.3), httpuv(v.1.6.15), tidyr(v.1.3.1), purrr(v.1.0.2), polyclip(v.1.10-6), ggforce(v.0.4.2), broom(v.1.0.5), xtable(v.1.8-4), restfulr(v.0.0.15), tidytree(v.0.4.6), fANCOVA(v.0.6-1), later(v.1.3.2), viridisLite(v.0.4.2), tibble(v.3.2.1), lmerTest(v.3.1-3), clusterProfiler(v.4.10.1), aplot(v.0.2.2), memoise(v.2.0.1), AnnotationDbi(v.1.64.1), GenomicAlignments(v.1.38.2), sva(v.3.50.0) and GSEABase(v.1.64.0)

message("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 31bacdf03b6f37f1277e44be276420f3db4d4bd6

## This is hpgltools commit: Mon Apr 8 13:21:41 2024 -0400: 31bacdf03b6f37f1277e44be276420f3db4d4bd6

#  message("Saving to ", savefile)
#  tmp <- sm(saveme(filename = savefile))

tmp <- loadme(filename = savefile)