The hyperparameter choice - pt. 1

In this notebook we will visualize the effect of the choice of elastic tree hyperparameters for the elastic and embedded tree. In this first part we will be working with the dataset by Guo et al., detailing early embryogenesis in the mouse. We will be needing MERLoT as well as some customized versions of the MERLoT functions.

We will read the dataset and the annotation, as well as the diffusion map coordinates.

# Read the example Guo dataset that is distributed with the package
DataFile= paste(find.package("merlot"), "/example/Guo2010.txt", sep="")
Dataset=ReadDataset(DataFile)
# Load the cell types
CellTypes=read.table(file=paste(find.package("merlot"), "/example/GuoFeatures.txt", sep=""), sep="\t", header = F, stringsAsFactors = F)
CellTypes=CellTypes[,2]
selected_colors=c("red", "orange", "yellow", "green", "cyan", "darkblue")
guo_colorcells=c()
guo_colorcells[which(CellTypes=="2C")]="red"
guo_colorcells[which(CellTypes=="4C")]="orange"
guo_colorcells[which(CellTypes=="8C")]="yellow"
guo_colorcells[which(CellTypes=="16C")]="green"
guo_colorcells[which(CellTypes=="32C")]="cyan"
guo_colorcells[which(CellTypes=="64C")]="darkblue"
CellCoordinates=read.table(file=paste(find.package("merlot"), "/example/GuoRotatedCoordinates.txt", sep=""), sep="\t", header = F, stringsAsFactors = F)
CellCoordinates=as.matrix(CellCoordinates[,1:2])

We calculate the scaffold tree without local averaging. The Guo dataset only has 428 cells; reducing that further will only make MERLoT’s job harder.

We will now calculate the elastic tree with different elasticity hyperparameters: the MERLoT defaults and the values that ElPiGraph.R proposes for the elastic_tree algorithm.

The exact number of nodes per branch differs, and the exact positions of endpoints and branchpoints do not overlap, but we see that there are no big differences between using the MERLoT hyperparameters (left) and the ElPiGraph.R ones (right).

For the embedding in gene space, we will use (1) the MERLoT elastic tree with MERLoT hyperparameters (benchmarked in original benchmark), (2) the ElPi elastic tree with ElPi hyperparameters (benchmarked in extended benchmark), and (3) the ElPi elastic tree with MERLoT hyperparameters (not benchmarked).

# Embedd the principal elastic tree on the gene expression space from which it was calculated.
EmbeddedTree= GenesSpaceEmbedding(ExpressionMatrix = Dataset$ExpressionMatrix, ElasticTree = ElasticTree)
nikoEmb <- nikoEmbed(ExpressionMatrix = Dataset$ExpressionMatrix, ElasticTree = nikoTree)
mixed = GenesSpaceEmbedding(ExpressionMatrix = Dataset$ExpressionMatrix, ElasticTree = nikoTree)

There is no immediate way to visualize the embedded tree, but we can plot the trees that result from the embedding in a schematic way, using the plot_flattened_tree function.

par(mfrow=c(1,3))
plot_flattened_tree(EmbeddedTree, cell_annot = CellTypes, node_size = "cells")
plot_flattened_tree(nikoEmb, cell_annot = CellTypes, node_size = "cells")

plot_flattened_tree(mixed, cell_annot = CellTypes, node_size = "cells")
par(mfrow=c(1,1))

Again, there is no big difference except for the absence of cell clumps in the endpoints of three of the four branches. We will now calculate pseudotimes for each cell in every embedded tree:

With the pseudotimes available, we can now plot the expression of genes through pseudotime, allowing another comparison between the three approaches.

par(mfrow=c(1, 3))
plot_pseudotime_expression_gene(GeneName = "Gata4" , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)
plot_pseudotime_expression_gene(GeneName = "Gata4" , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)

plot_pseudotime_expression_gene(GeneName = "Gata4" , EmbeddedTree = mixed, Pseudotimes = mixed_times)
par(mfrow=c(1,1))

Here we see a big difference for the first time. The default ElPi parameters allow the expression of the pseudocells to overfit (middle). On the other hand, not using the MERLoT hyperparameters for the elastic tree doesn’t seem to affect the end result that much (left and right). We conclude that while the ElPi hyperparameters might produce good elastic trees (especially in conjuction with local averaging), they do not lend themselves to the creation of embedded trees, since they will overfit the data and produce less useful gene trajectories over time.

We include some more examples below. Always, from left to right: MERLoT+MERLoT, ElPi+ElPi, ElPi+MERLoT.

# Plot gene expression profile as a function of pseudotime
par(mfrow=c(5, 3))
gene <- "Runx1"
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)

plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = mixed, Pseudotimes = mixed_times)
gene <- "Nanog"
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)

plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = mixed, Pseudotimes = mixed_times)

gene <- "Fgf4"
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)

plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = mixed, Pseudotimes = mixed_times)
gene <- "Fgfr2"
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)

plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = mixed, Pseudotimes = mixed_times)

gene <- "Klf2"
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = EmbeddedTree, Pseudotimes = Pseudotimes)
plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = nikoEmb, Pseudotimes = nikotimes)

plot_pseudotime_expression_gene(GeneName = gene , EmbeddedTree = mixed, Pseudotimes = mixed_times)
par(mfrow=c(1,1))