L. panamensis 202209: Visualizing TMRC2 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

2.1.1 Library sizes

hs_lib <- plot_libsize(hs_macrophage)
pp(file="images/hs_macrophage_libsize.png")
hs_lib$plot
dev.off()

## png 
##   2

hs_lib$plot

hs_lib_log <- plot_libsize(hs_macrophage, yscale="log2")
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="images/hs_macrophage_nonzero.png")
hs_non$plot
dev.off()

## png 
##   2

hs_non$plot

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)

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

pp(file="images/hs_macrophage_boxplot.png")
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.png")
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.png")
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.

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

## Removing 10021 low-count genes (11460 remaining).

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

hs_pca <- plot_pca(hs_norm, plot_title = "PCA of macrophage expression values",
                   plot_labels = FALSE)
pp(file="images/hs_macrophage_norm_pca.png")
hs_pca$plot
dev.off()

## png 
##   2

hs_pca$plot

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

## Removing 10021 low-count genes (11460 remaining).

## Setting 619 low elements to zero.

## transform_counts: Found 619 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_labels = FALSE)
pp(file="images/hs_macrophage_nb_pca.png")
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

```{r correlation_ corheat <- plot_corheat(hs_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(hs_norm, plot_title = “Euclidean heatmap of macrophage expression values”) disheat$plot

corheat <- plot_disheat(hs_nb, plot_title = “Euclidean heatmap of macrophage expression values”) corheat$plot

plot_sm(hs_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):

[igv/igv_snapshot_checking_zymodeme_sample.png](igv/igv_snapshot_checking_zymodeme_sample.png)

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.

## Parasite data

### Library sizes, non-zero plot, and distribution


```r
lp_lib <- plot_libsize(lp_macrophage)
pp(file="images/lp_macrophage_libsize.png")
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.png")
lp_non$plot
dev.off()

## png 
##   2

lp_non$plot

lp_box <- plot_boxplot(lp_macrophage)

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

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

## png 
##   2

lp_box

filter_plot <- plot_libsize_prepost(lp_macrophage)
pp(file="images/lp_macrophage_lowgene.png")
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.png")
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.

2.3 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).

2.3.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 188 low-count genes (8522 remaining).

## transform_counts: Found 4 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.png")
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 188 low-count genes (8522 remaining).

## Setting 45 low elements to zero.

## transform_counts: Found 45 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.png")
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

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

## transform_counts: Found 4 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.png")
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 190 low-count genes (8520 remaining).

## Setting 43 low elements to zero.

## transform_counts: Found 43 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.png")
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(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 b80a01ab8136fffc1458b8a0a80d9ce5a8d64aa5

## This is hpgltools commit: Tue Sep 13 12:24:31 2022 -0400: b80a01ab8136fffc1458b8a0a80d9ce5a8d64aa5

## Saving to tmrc2_macrophage_visualization_202209.rda.xz

tmp <- loadme(filename = savefile)