TMRC2 202301: 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$plot

hs_lib_log <- plot_libsize(hs_macrophage)
hs_lib_log$plot

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$plot

## Warning: ggrepel: 51 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$lowgene_plot

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

## 
## 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$plot

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

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

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

## Deleting the file analyses/macrophage_de/202304/read_counts/hs_macr-v202304.xlsx before writing the tables.

## 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.
## 
## Using expt method.
## 
## Performing correlation.
## 
## Using expt method.
## 
## Performing distance.
## 
## 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 checkModelStatus(fit, showWarnings = showWarnings, colinearityCutoff = colinearityCutoff,  : 
##   The variables specified in this model are redundant,
## so the design matrix is not full rank

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

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

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

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-v{ver}.svg"))
macr_pca$plot
dev.off()

## png 
##   2

macr_pca$plot

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

## 
##  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$plot

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

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

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

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$plot

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$plot

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.2.2 Correlation heatmaps

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

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

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

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

plot_sm(macr_norm)$plot

## Using expt method.

## Performing correlation.

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):

igv/igv_snapshot_checking_zymodeme_sample.svg

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

## 
## Total:162 s

donor_drug_zymo_varpart$partition_plot

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

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

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

tt$pvalue_plots$MF

tt$pvalue_plots$REAC

tt$pvalue_plots$WP

tt$pvalue_plots$BP

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

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

tt$pvalue_plots$MF

tt$pvalue_plots$REAC

tt$pvalue_plots$BP

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

## 
## Total:154 s

time_drug_zymo_varpart$partition_plot

3.2 Drug zymodeme donor

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

## 
## Total:159 s

donor_drug_zymo_varpart$partition_plot

donor_drug_zymo_varpart$model_used

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

drug_zymo_donor_varpart$partition_plot

drug_zymo_donor_varpart$model_used

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

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

hs_donors <- set_expt_conditions(hs_macr, fact = "donor")

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

donor_norm <- normalize_expt(hs_donors, 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_donors.svg")
plot_pca(donor_norm, plot_labels = FALSE)$plot
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.svg")
plot_pca(donor_nb, plot_labels = FALSE)$plot
dev.off()

## png 
##   2

4.2 Drug comparison

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

## 
## 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)$plot
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)$plot
dev.off()

## png 
##   2

4.3 Parasite data

4.3.1 Library sizes, non-zero plot, and distribution

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

## Deleting the file analyses/macrophage_de/202304/read_counts/lp_macrophage_reads-v202304.xlsx before writing the tables.

## 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.
## 3510 entries are 0.  We are on a log scale, adding 1 to the data.
## 
## Using expt method.
## 
## Performing correlation.
## 
## Using expt method.
## 
## Performing distance.
## 
## Changed 3510 zero count features.
## 
## Naively calculating coefficient of variation/dispersion with respect to condition.
## 
## Finished calculating dispersion estimates.
## 
## `geom_smooth()` using formula = 'y ~ x'
## 
## Total:69 s
## 
## `geom_smooth()` using formula = 'y ~ x'
## 
## Total:70 s

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

## Deleting the file analyses/macrophage_de/202304/read_counts/lp_macrophage_nosb_reads-v202304.xlsx before writing the tables.
## 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.3148 entries are 0.  We are on a log scale, adding 1 to the data.
## Using expt method.
## Performing correlation.
## Using expt method.
## Performing distance.
## Changed 3148 zero count features.
## Naively calculating coefficient of variation/dispersion with respect to condition.
## Finished calculating dispersion estimates.
## `geom_smooth()` using formula = 'y ~ x'varpart sees only 1 batch, adjusting the model accordingly.

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Variable in formula not found in data: condition

## Retrying with only condition in the model.
## 
## Total:52 s
## `geom_smooth()` using formula = 'y ~ x'varpart sees only 1 batch, adjusting the model accordingly.

## Error in .fitExtractVarPartModel(exprObj, formula, data, REML = REML,  : 
##   Variable in formula not found in data: condition

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

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

## png 
##   2

lp_lib$plot

lp_non <- plot_nonzero(lp_macrophage)

## 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
dev.off()

## png 
##   2

lp_non$plot

lp_box <- plot_boxplot(lp_macrophage)

## 3510 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$lowgene_plot

## 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 169 low-count genes (8541 remaining).

## transform_counts: Found 23 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$plot

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

## Removing 169 low-count genes (8541 remaining).

## Setting 134 low elements to zero.

## transform_counts: Found 134 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
dev.off()

## png 
##   2

lp_nb_pca$plot

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

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

plot_sm(lp_norm)$plot

## Using expt method.

## Performing correlation.

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 171 low-count genes (8539 remaining).

## transform_counts: Found 23 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$plot

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

## Removing 171 low-count genes (8539 remaining).

## Setting 110 low elements to zero.

## transform_counts: Found 110 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$plot

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

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

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message("This is hpgltools commit: ", get_git_commit())
  message("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 bbc75e24b763faa635cb62c86fc51c0efb3424a1

## This is hpgltools commit: Fri Apr 21 14:33:16 2023 -0400: bbc75e24b763faa635cb62c86fc51c0efb3424a1

## Saving to tmrc2_macrophage_visualization_202304.rda.xz

tmp <- loadme(filename = savefile)