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 <- load_orgdb_annotations(
  pan_db,
  keytype="gid",
  fields=c("annot_gene_entrez_id", "annot_gene_name",
           "annot_strand", "annot_chromosome", "annot_cds_length",
           "annot_gene_product"))$genes

## Extracted all gene ids.

## Attempting to select: ANNOT_GENE_PRODUCT, ANNOT_GENE_TYPE, ANNOT_LOCATION_TEXT, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_NAME, ANNOT_STRAND, ANNOT_CHROMOSOME, ANNOT_CDS_LENGTH, ANNOT_GENE_PRODUCT

## 'select()' returned 1:1 mapping between keys and columns

lp_go <- sm(load_orgdb_go(pan_db))
lp_lengths <- all_lp_annot[, c("gid", "annot_cds_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_title="PCA of parasite expression values")$plot

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

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

plot_tsne(all_norm, plot_title="TSNE of parasite expression values")$plot

plot_corheat(all_norm, title="Correlation heatmap of parasite expression values
(Same legend as above)")$plot

plot_sm(all_norm)$plot

## Performing correlation.

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

3.2.1 Notes

The following samples are much lower coverage:

• TMRC20002
• TMRC20006
• TMRC20007
• TMRC20008

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

4 Zymodeme analyses

The following sections perform a series of analyses which seek to elucidate differences between the zymodemes 2.2 and 2.3 either through differential expression or variant profiles.

4.1 Differential expression

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

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

4.1.1 A couple plots from the differential expression

4.1.1.1 Number of genes in agreement among DE methods, 2.3 more than 2.2

zy_table[["venns"]][[1]][["p_lfc1"]][["up_noweight"]]

4.1.1.2 Number of genes in agreement among DE methods, 2.2 more than 2.3

zy_table[["venns"]][[1]][["p_lfc1"]][["down_noweight"]]

4.1.1.3 MA plot of the differential expression between the zymodemes.

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

4.1.1.4 goseq ontology plots of groups of genes, 2.3 more than 2.2

zy_go_up$pvalue_plots$bpp_plot_over

4.1.1.5 goseq ontology plots of groups of genes, 2.2 more than 2.3

zy_go_down$pvalue_plots$bpp_plot_over

4.2 Zymodeme enzyme gene IDs

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.

4.2.1 Expression levels of zymodeme genes

The following creates a colorspace (red to green) heatmap showing the observed expression of these genes in every sample.

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)

4.3 Empirically observed Zymodeme genes from differential expression analysis

In contrast, the following plots take the set of genes which are shared among all differential expression methods (|lfc| >= 1.0 and adjp <= 0.05) and use them to make categories of genes which are increased in 2.3 or 2.2.

shared_zymo <- intersect_significant(zy_table)

up_shared <- shared_zymo[["ups"]][[1]][["data"]][["all"]]
rownames(up_shared)

##  [1] "LPAL13_000033300" "LPAL13_000038400" "LPAL13_310031300" "LPAL13_000012000"
##  [5] "LPAL13_050005000" "LPAL13_210015500" "LPAL13_310039200" "LPAL13_340039600"
##  [9] "LPAL13_000038500" "LPAL13_000041000" "LPAL13_200013000" "LPAL13_330021800"
## [13] "LPAL13_340039700" "LPAL13_170015400" "LPAL13_240009700" "LPAL13_140019300"
## [17] "LPAL13_140019100" "LPAL13_320038700" "LPAL13_210005000" "LPAL13_280037900"
## [21] "LPAL13_230011200" "LPAL13_100010600" "LPAL13_140019200" "LPAL13_310028500"
## [25] "LPAL13_230011500" "LPAL13_000010600" "LPAL13_230011300" "LPAL13_350073200"
## [29] "LPAL13_250025700" "LPAL13_040007800" "LPAL13_200012900" "LPAL13_020006900"

upshared_expt <- exclude_genes_expt(all_norm, ids=rownames(up_shared), method="keep")

## Before removal, there were 8623 entries.

## Now there are 32 entries.

## Percent kept: 0.358, 0.288, 0.263, 0.255, 0.278, 0.264, 0.254, 0.372, 0.260, 0.366, 0.358, 0.270, 0.254, 0.367, 0.267, 0.267, 0.366

## Percent removed: 99.642, 99.712, 99.737, 99.745, 99.722, 99.736, 99.746, 99.628, 99.740, 99.634, 99.642, 99.730, 99.746, 99.633, 99.733, 99.733, 99.634

4.3.1 Heatmap of zymodeme gene expression increased in 2.3 vs. 2.2

test <- plot_sample_heatmap(upshared_expt, row_label=rownames(up_shared))

4.3.2 Heatmap of zymodeme gene expression increased in 2.2 vs. 2.3

down_shared <- shared_zymo[["downs"]][[1]][["data"]][["all"]]
downshared_expt <- exclude_genes_expt(all_norm, ids=rownames(down_shared), method="keep")

## Before removal, there were 8623 entries.

## Now there are 69 entries.

## Percent kept: 0.461, 0.816, 0.742, 0.757, 0.787, 0.751, 0.769, 0.423, 0.790, 0.406, 0.435, 0.791, 0.749, 0.433, 0.779, 0.796, 0.438

## Percent removed: 99.539, 99.184, 99.258, 99.243, 99.213, 99.249, 99.231, 99.577, 99.210, 99.594, 99.565, 99.209, 99.251, 99.567, 99.221, 99.204, 99.562

test <- plot_sample_heatmap(downshared_expt, row_label=rownames(down_shared))

5 SNP 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]

5.1 Plot of SNP profiles for zymodemes

The following plot shows the SNP profiles of all samples (old and new) where the colors at the top show either the 2.2 strains (orange), 2.3 strains (green), the previous samples (purple), or the various lab strains (pink etc).

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"

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

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

## Using a subset expression.

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

pheno_snps <- sm(count_expt_snps(pheno, annot_column="bcftable"))

xref_prop <- table(pheno_snps$conditions)
pheno_snps$conditions

##  [1] "z2.3" "z2.2" "z2.2" "z2.2" "z2.3" "z2.2" "z2.3" "z2.3" "z2.2" "z2.2"
## [11] "z2.3" "z2.2" "z2.2" "z2.3"

idx_tbl <- exprs(pheno_snps) > 5
new_tbl <- data.frame(row.names=rownames(exprs(pheno_snps)))
for (n in names(xref_prop)) {
  new_tbl[[n]] <- 0
  idx_cols <- which(pheno_snps[["conditions"]] == n)
  prop_col <- rowSums(idx_tbl[, idx_cols]) / xref_prop[n]
  new_tbl[n] <- prop_col
}
new_tbl[["ratio"]] <- (new_tbl[["z2.2"]] - new_tbl[["z2.3"]])
keepers <- grepl(x=rownames(new_tbl), pattern="LpaL13")
new_tbl <- new_tbl[keepers, ]
new_tbl[["SNP"]] <- rownames(new_tbl)
new_tbl[["Chromosome"]] <- gsub(x=new_tbl[["SNP"]], pattern="chr_(.*)_pos_.*", replacement="\\1")
new_tbl[["Position"]] <- gsub(x=new_tbl[["SNP"]], pattern=".*_pos_(\\d+)_.*", replacement="\\1")
new_tbl <- new_tbl[, c("SNP", "Chromosome", "Position", "ratio")]
CMplot(new_tbl)

## Error in CMplot(new_tbl): could not find function "CMplot"

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 d4bf620dede68821356e0745b7372701800699f1

## This is hpgltools commit: Fri Oct 16 12:52:28 2020 -0400: d4bf620dede68821356e0745b7372701800699f1

## Saving to 01_annotation_v202009.rda.xz

tmp <- loadme(filename=savefile)