L. panamensis 202009: Sample Estimation

1 Introduction

This document is intended to provide a general overview of the TMRC2 samples which have thus far been sequenced. In some cases, this includes only those samples starting in 2019; in other instances I am including our previous (2015-2016) samples.

In all cases the processing performed was:

1. Default trimming was performed.
2. Hisat2 was used to map the remaining reads against the Leishmania panamensis genome revision 36.
3. The alignments from hisat2 were used to count reads/gene against the revision 36 annotations with htseq.
4. These alignments were also passed to the pileup functionality of samtools and the vcf/bcf utilities in order to make a matrix of all observed differences between each sample with respect to the reference.

The analyses in this document use the matrices of counts/gene from #3 and variants/position from #4 in order to provide some images and metrics describing the samples we have sequenced so far.

2 Annotations

Everything which follows depends on the Existing TriTrypDB annotations revision 46, circa 2019. The following block loads a database of these annotations and turns it into a matrix where the rows are genes and columns are all the annotation types provided by TriTrypDB.

The same database was used to create a matrix of orthologous genes between L.panamensis and all of the other species in the TriTrypDB.

tt <- sm(library(EuPathDB))
tt <- sm(library(org.Lpanamensis.MHOMCOL81L13.v46.eg.db))
pan_db <- org.Lpanamensis.MHOMCOL81L13.v46.eg.db
all_fields <- columns(pan_db)

all_lp_annot <- sm(load_orgdb_annotations(
  pan_db,
  keytype="gid",
  fields="all")$genes)

lp_go <- sm(load_orgdb_go(pan_db))
lp_lengths <- all_lp_annot[, c("gid", "length")]
colnames(lp_lengths)  <- c("ID", "length")

orthos <- sm(EuPathDB::extract_eupath_orthologs(db=pan_db))

hisat_annot <- all_lp_annot
## rownames(hisat_annot) <- paste0("exon_", rownames(hisat_annot), ".E1")

3 Sample Estimation

The process of sample estimation takes two primary inputs:

1. The sample sheet, which contains all the metadata we currently have on hand, including filenames for the outputs of #3 and #4 above.
2. The gene annotations.

An expressionset is primary data structure used in R to examine RNASeq data. It is comprised of annotations, metadata, and expression data. In the case of our processing pipeline, the location of the expression data is provided by the filenames in the metadata.

3.1 Generate expressionsets

The first lines of the following block create the Expressionset. All of the following lines perform various normalizations and generate plots from it.

lp_expt <- sm(create_expt("sample_sheets/tmrc2_samples_20200915.xlsx",
                          gene_info=hisat_annot,
                          file_column="lpanamensisv36hisatfile"))

libsizes <- plot_libsize(lp_expt)

## The scale difference between the smallest and largest
## libraries is > 10. Assuming a log10 scale is better, set scale=FALSE if not.

libsizes$plot

## I think samples 7,10 should be removed at minimum, probably also 9,11
nonzero <- plot_nonzero(lp_expt)
nonzero$plot

plot_boxplot(lp_expt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

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

## Changed 1920 zero count features.

3.2 Distribution Visualization

Najib’s favorite plots are of course the PCA/TNSE. These are nice to look at in order to get a sense of the relationships between samples. They also provide a good opportunity to see what happens when one applies different normalizations, surrogate analyses, filters, etc. In addition, one may set different experimental factors as the primary ‘condition’ (usually the color of plots) and surrogate ‘batches’.

pheno_expt <- set_expt_conditions(lp_expt, fact="phenotypiccharacteristics")
all_norm <- normalize_expt(pheno_expt, norm="quant", transform="log2", convert="cpm",
                           batch=FALSE, filter=TRUE)

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 155 low-count genes (8623 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

plot_pca(all_norm)$plot

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

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

## Note, samples 2, 6, 7, 8 are all very low coverage compared to the later samples.

plot_tsne(all_norm)$plot

plot_corheat(all_norm)$plot

plot_sm(all_norm)$plot

## Performing correlation.

## plot_variance_coefficients(all_norm)$plot
## plot_sample_cvheatmap(all_norm)$plot

At this time, we do not have very many samples, so the set of metrics/plots is fairly limited. There is really only one factor in the metadata which we can use for performing differential expression analyses, the ‘zymodeme’.

zy_expt <- subset_expt(pheno_expt, subset="condition=='z2.2'|condition=='z2.3'")

## Using a subset expression.

## There were 17, now there are 14 samples.

zy_de <- sm(all_pairwise(zy_expt, filter=TRUE, model_batch="svaseq"))
zy_table <- sm(combine_de_tables(zy_de, excel=glue::glue("excel/zy_tables-v{ver}.xlsx")))

zy_table$plots[[1]][["deseq_ma_plots"]][["plot"]]

zy_sig <- sm(extract_significant_genes(zy_table, excel=glue::glue("excel/zy_sig-v{ver}.xlsx")))

lp_go <- lp_go[, c("GID", "GOALL")]

## Error in `[.data.frame`(lp_go, , c("GID", "GOALL")): undefined columns selected

colnames(lp_go) <- c("ID", "GO")

## Gene categories more represented in the 2.3 group.
zy_go_up <- sm(simple_goseq(sig_genes=zy_sig[["deseq"]][["ups"]][[1]], go_db=lp_go, length_db=lp_lengths))

zy_go_up$pvalue_plots$bpp_plot_over

## Gene categories more represented in the 2.2 group.
zy_go_down <- sm(simple_goseq(sig_genes=zy_sig[["deseq"]][["downs"]][[1]], go_db=lp_go, length_db=lp_lengths))

zy_go_down$pvalue_plots$bpp_plot_over

4 Zymodeme genes

Najib read me an email listing off the gene names associated with the zymodeme classification. I took those names and cross referenced them against the Leishmania panamensis gene annotations and found the following:

They are:

1. ALAT: LPAL13_120010900 – alanine aminotransferase
2. ASAT: LPAL13_340013000 – aspartate aminotransferase
3. G6PD: LPAL13_000054100 – glucase-6-phosphate 1-dehydrogenase
4. NH: LPAL13_14006100, LPAL13_180018500 – inosine-guanine nucleoside hydrolase
5. MPI: LPAL13_320022300 (maybe) – mannose phosphate isomerase (I chose phosphomannose isomerase)

Given these 6 gene IDs (NH has two gene IDs associated with it), I can do some looking for specific differences among the various samples.

my_genes <- c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
              "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300",
              "other")
my_names <- c("ALAT", "ASAT", "G6PD", "NHv1", "NHv2", "MPI", "other")

zymo_expt <- exclude_genes_expt(all_norm, ids=my_genes, method="keep")

## Before removal, there were 8623 entries.

## Now there are 6 entries.

## Percent kept: 0.086, 0.081, 0.084, 0.084, 0.083, 0.086, 0.083, 0.084, 0.083, 0.084, 0.083, 0.083, 0.085, 0.085, 0.083, 0.083, 0.083

## Percent removed: 99.914, 99.919, 99.916, 99.916, 99.917, 99.914, 99.917, 99.916, 99.917, 99.916, 99.917, 99.917, 99.915, 99.915, 99.917, 99.917, 99.917

test <- plot_sample_heatmap(zymo_expt, row_label=my_names)

5 Variant profiles

In this block, I am combining our previous samples and our new samples in the hopes of finding variant positions which help elucidate aspects of either the new or old samples. In other words, we do not know the zymodeme annotations for the old samples nor the strain identities (or the shortcut ‘chronic vs. self-healing’) for the new samples. We may be able to make educated guesses given the variant profiles. There are some differences in how the previous and current data sets were analyzed (though I have since redone the old samples so it should be trivial to remove those differences now).

old_expt <- sm(create_expt("sample_sheets/tmrc2_samples_20191203.xlsx",
                           file_column="tophat2file"))

tt <- lp_expt$expressionset
rownames(tt) <- gsub(pattern="^exon_", replacement="", x=rownames(tt))
rownames(tt) <- gsub(pattern="\\.E1$", replacement="", x=rownames(tt))
lp_expt$expressionset <- tt

tt <- old_expt$expressionset
rownames(tt) <- gsub(pattern="^exon_", replacement="", x=rownames(tt))
rownames(tt) <- gsub(pattern="\\.1$", replacement="", x=rownames(tt))
old_expt$expressionset <- tt

new_snps <- sm(count_expt_snps(lp_expt, annot_column="bcftable"))
old_snps <- sm(count_expt_snps(old_expt, annot_column="bcftable", snp_column=2))

both_snps <- combine_expts(new_snps, old_snps)
both_norm <- sm(normalize_expt(both_snps, transform="log2", convert="cpm", filter=TRUE))

strains <- both_norm[["design"]][["strain"]]
both_norm <- set_expt_conditions(both_norm, fact=strains)

## Error in names(object) <- nm: 'names' attribute [11] must be the same length as the vector [10]

tt <- plot_disheat(both_norm)

snp_sets <- get_snp_sets(both_snps, factor="condition")

## The factor undefined has 17 rows.

## The factor sh has 13 rows.

## The factor chr has 14 rows.

## The factor inf has 6 rows.

## Iterating over 722 elements.

both_expt <- combine_expts(lp_expt, old_expt)
snp_genes <- sm(snps_vs_genes(both_expt, snp_sets, expt_name_col="chromosome"))

snp_subset <- sm(snp_subset_genes(
  both_expt, both_snps,
  genes=c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
          "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300")))

zymo_heat <- plot_sample_heatmap(snp_subset, row_label=rownames(exprs(snp_subset)))

6 Zymodeme for new samples

The heatmap produced here should show the variants only for the zymodeme genes.

new_sets <- get_snp_sets(new_snps, factor="phenotypiccharacteristics")

## The factor Laboratory line has 3 rows.

## The factor Laboratory line antimony resistant  has only 1 row.

## The factor Laboratory line miltefosine resistant  has only 1 row.

## The factor z2.2 has 8 rows.

## The factor z2.3 has 6 rows.

## Iterating over 620 elements.

snp_genes <- sm(snps_vs_genes(pheno_expt, new_sets, expt_name_col="chromosome"))
new_zymo_norm  <-  normalize_expt(new_snps, filter=TRUE, convert="cpm", norm="quant", transform=TRUE)

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

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

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

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

## Step 5: not doing batch correction.

new_zymo_norm <- set_expt_conditions(new_zymo_norm, fact="phenotypiccharacteristics")
zymo_heat <- plot_disheat(new_zymo_norm)

zymo_subset <- sm(snp_subset_genes(
  pheno_expt, new_snps,
  genes=c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
          "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300")))

zymo_subset <- set_expt_conditions(zymo_subset, fact="phenotypiccharacteristics")
zymo_heat <- plot_sample_heatmap(zymo_subset, row_label=rownames(exprs(snp_subset)))

des <- both_norm$design
undef_idx <- is.na(des[["strain"]])
des[undef_idx, "strain"] <- "unknown"
library(Heatplus)
##hmcols <- colorRampPalette(c("yellow","black","darkblue"))(256)
correlations <- hpgl_cor(exprs(both_norm))

zymo_missing_idx <- is.na(des[["phenotypiccharacteristics"]])
des[zymo_missing_idx, "phenotypiccharacteristics"] <- "unknown"
mydendro <- list(
  "clustfun" = hclust,
  "lwd" = 2.0)
col_data <- as.data.frame(des[, c("phenotypiccharacteristics")])
row_data <- as.data.frame(des[, c("strain")])
colnames(col_data) <- c("condition")
colnames(row_data) <- c("strain")
myannot <- list(
  "Col" = list("data" = col_data),
  "Row" = list("data" = row_data))
myclust <- list("cuth" = 1.0,
                "col" = BrewerClusterCol)
mylabs <- list(
  "Row" = list("nrow" = 4),
  "Col" = list("nrow" = 4))
hmcols <- colorRampPalette(c("darkblue", "beige"))(170)
map1 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  col=hmcols)

## Warning in breakColors(breaks, col): more colors than classes: ignoring 3 last
## colors

plot(map1)

7 Using Variant profiles to make guesses about strains and chronic/self-healing

The following uses the same information to make some guesses about the strains used in the new samples.

des <- both_norm$design
undef_idx <- is.na(des[["strain"]])
des[undef_idx, "strain"] <- "unknown"
##hmcols <- colorRampPalette(c("yellow","black","darkblue"))(256)
correlations <- hpgl_cor(exprs(both_norm))

mydendro <- list(
  "clustfun" = hclust,
  "lwd" = 2.0)
col_data <- as.data.frame(des[, c("condition")])
row_data <- as.data.frame(des[, c("strain")])
colnames(col_data) <- c("condition")
colnames(row_data) <- c("strain")
myannot <- list(
  "Col" = list("data" = col_data),
  "Row" = list("data" = row_data))
myclust <- list("cuth" = 1.0,
                "col" = BrewerClusterCol)
mylabs <- list(
  "Row" = list("nrow" = 4),
  "Col" = list("nrow" = 4))
hmcols <- colorRampPalette(c("darkblue", "beige"))(170)
map1 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  col=hmcols)

## Warning in breakColors(breaks, col): more colors than classes: ignoring 3 last
## colors

plot(map1)

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 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394

## This is hpgltools commit: Wed Aug 19 10:11:52 2020 -0400: 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394

## Saving to 01_annotation_v202009.rda.xz

tmp <- loadme(filename=savefile)