2. Trajectory inference using Treefit and Seurat
As the Treefit package is meant to be used with gene expression data (either raw counts or normalized values) that have been already quality controlled, it is a good idea to use Seurat or other popular toolkits designed for the quality control of single-cell gene expression data (e.g., filtering cells and normalizing the data) before the Treefit analysis.
2.1. Loading the dataset
We analyze the dataset provided in Trapnell el al., 2014, Nature Biotechnology. This dataset was originally acquired for the purpose of studying the process of differentiation of human myoblasts that can be reasonably explained by the linear model (i.e., non-branching, chain-like structure). However, as the authors discuss in detail in the original article, their trajectory inference software based on a minimum spanning tree (MST) helped detect contamination with another cell lineage which was not related to myogenesis and so revealed that the bifurcation tree model (i.e., Y-shaped tree structure) fits the data better than the linear model.
To use this dataset, we download and read the dataset provided in dynverse as follows (see https://zenodo.org/record/1443566 for details).
trapnell.path <- "myoblast-differentiation_trapnell.rds"
if (!file.exists(trapnell.path)) {
download.file("https://zenodo.org/record/1443566/files/real/gold/myoblast-differentiation_trapnell.rds?download=1",
trapnell.path,
mode="wb")
}
trapnell.dynverse <- readRDS(trapnell.path)2.2. Preprocessing with Seurat
From now, we demonstrate how to preprocess the above data with Seurat. In order to handle data in Seurat, we first need to create a Seurat object. When creating a Seurat object from dynverse data, we must be aware that the roles of the rows and columns are reversed and therefore we need to transpose the matrix. In dynverse data, the rows and columns correspond to cells and genes, respectively, and vice versa in Seurat objects.
Here we assume that we wish to perform the following preprocessing.
- Excluding each gene that is only expressed in less than 3 cells
- Removing each cell that only expresses less than 200 genes
This can be done by setting the designated parameters as follows.
trapnell <- Seurat::CreateSeuratObject(counts=t(trapnell.dynverse$count),
min.cells=3,
min.features=200)## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')The “Feature names cannot have underscores …” warning message can be ignored. It’s caused by difference between dynverse and Seurat. Dynverse uses underscores for gene names but Seurat doesn’t allow it. Seurat replaces underscores with dashes automatically with the warning message.
2.3. Visualization with Seurat
Although pictures may not always describe reality, we usually visualize data before starting detailed analysis. Here we create a scatter plot by principal component analysis (PCA). Below are a few lines of code to do this and the result is shown in Figure 1.
trapnell <- Seurat::FindVariableFeatures(trapnell, verbose=FALSE)
trapnell <- Seurat::ScaleData(trapnell, verbose=FALSE)
trapnell <- Seurat::RunPCA(trapnell, verbose=FALSE)
plot(Seurat::Embeddings(trapnell))
Figure 1. The result of PCA for the myoblast dataset
It is hard to reach the truth from the PCA result shown in Figure 1. Unlike the tree-like toy data analyzed in the previous tutorial (vignette("treefit")), we cannot easily recognize the underlying tree structure for various reasons (e.g., high level of noise, contamination with different types of cells, and the inaccuracy of two-dimensional visualization of high dimensional data).
Figure 2b in Trapnell el al., 2014, Nature Biotechnology shows that this dataset contains three distinct types of cells (i.e., proliferating cells, differentiating myoblasts, and contaminating interstitial mesenchymal cells), which led the authors to adopt a Y-shaped tree model (i.e., a star tree with three arms) to explain the data. However, it is difficult to quantify the reliability of this tree model and to confidently predict that the best-fit tree consists of three principal paths.
