% \VignetteIndexEntry{Monocle: Differential expression and time-series analysis for single-cell RNA-Seq and qPCR experiments.} % \VignetteEngine{knitr::knitr} % \VignetteDepends{} % \VignettePackage{monocle} \documentclass[10pt,oneside]{article} \newcommand{\thetitle}{Monocle: Cell counting, differential expression, and trajectory analysis for single-cell RNA-Seq experiments} %\usepackage[pdftitle={\thetitle},pdfauthor={Wolfgang Huber}]{whbiocvignette} \usepackage{whbiocvignette} % \usepackage{times} %\usepackage{hyperref} %\usepackage{verbatim} %\usepackage{graphicx} %\usepackage{fancybox} %\usepackage{color} \title{\textsf{\textbf{\thetitle}}} \author{Cole Trapnell\\[1em]University of Washington,\\ Seattle, Washington, USA\\ \texttt{coletrap@uw.edu} \and Davide Cacchiarelli\\[1em]Harvard University,\\ Cambridge, Massachussetts, USA\\ \texttt{davide@broadinstitute.org}\and Xiaojie Qiu\\[1em]University of Washington,\\ Seattle, Washington, USA\\ \texttt{xqiu@uw.edu} } \begin{document} <>= library(Biobase) library(knitr) library(reshape2) library(ggplot2) knitr::opts_chunk$set(autodep=TRUE, cache=FALSE, warning=FALSE) set.seed(0) @ %def \maketitle \begin{abstract} Single cell gene expression studies enable profiling of transcriptional regulation during complex biological processes and within highly hetergeneous cell populations. These studies allow discovery of genes that identify certain subtypes of cells, or that mark a particular intermediate states during a biological process. In many single cell studies, individual cells are executing through a gene expression program in an unsynchronized manner. In effect, each cell is a snapshot of the transcriptional program under study. The package \Rpackage{monocle} provides tools for analyzing single-cell expression experiments. Monocle introduced the strategy of ordering single cells in \emph{pseudotime}, placing them along a trajectory corresponding to a biological process such as cell differentiation. Monocle learns this trajectory directly from the data, in either a fully unsupervised or a semi-supervised manner. It also performs differential gene expression and clustering to identify important genes and cell states. It is designed for RNA-Seq studies, but can be used with other assays. For more information on the algorithm at the core of \Rpackage{monocle}, or to learn more about how to use single cell RNA-Seq to study a complex biological process, see Trapnell and Cacchiarelli \emph{et al}\cite{TRAPNELL_CACCHIARELLI} \end{abstract} \tableofcontents <>= library(HSMMSingleCell) library(monocle) data(HSMM_expr_matrix) data(HSMM_gene_annotation) data(HSMM_sample_sheet) @ %def \section{Introduction} The \Rpackage{monocle} package provides a toolkit for analyzing single cell gene expression experiments. This vignette provides an overview of a single cell RNA-Seq analysis workflow with Monocle. Monocle was originally developed to analyze dynamic biological processes such as cell differentiation, although it also supports simpler experimental settings. \textbf{Monocle 2} includes new, improved algorithms for classifying and counting cells, performing differential expression analysis between subpopulations of cells, and reconstructing cellular trajectories. Monocle 2 has also been re-engineered to work well with very large single-cell RNA-Seq experiments containing tens of thousands of cells or even more. \newline Monocle can help you perform three main types of analysis: \begin{enumerate} \item \textbf{Classifying and counting cells.} Single-cell RNA-Seq experiments allow you to discover new (and possibly rare) subtypes of cells. Monocle helps you identify them. \item \textbf{Constructing single-cell trajectories.} In development, disease, and throughout life, cells transition from one state to another. Monocle helps you discover these transitions. \item \textbf{Differential expression analysis.} Characterizing new cell types and states begins with comparing them to other, better understood cells. Monocle includes a sophisticated but easy to use system for differential expression. \end{enumerate} Before we look at Monocle's functions for each of these common analysis tasks, let's see how to load up single-cell datasets in Monocle. \section{Getting started with Monocle} The \Rpackage{monocle} package takes a matrix of gene expression values as calculated by Cufflinks\cite{Trapnell:2012kp} or another gene expression estimation program. Monocle can work with relative expression values (e.g. FPKM or TPM units) or absolute transcript counts (e.g. from UMI experiments). Although Monocle can be used with raw read counts, these are not directly proportional to expression values unless you normalize them by length, so some Monocle functions could produce nonsense results. If you don't have UMI counts, We recommend you load up FPKM or TPM values instead of raw read counts. \subsection{The CellDataSet class} \Rpackage{monocle} holds single cell expression data in objects of the \Rclass{CellDataSet} class. The class is derived from the Bioconductor \Rclass{ExpressionSet} class, which provides a common interface familiar to those who have analyzed microarray experiments with Bioconductor. The class requires three input files: \begin{enumerate} \item \Robject{exprs}, a numeric matrix of expression values, where rows are genes, and columns are cells \item \Robject{phenoData}, an \Rpackage{AnnotatedDataFrame} object, where rows are cells, and columns are cell attributes (such as cell type, culture condition, day captured, etc.) \item \Robject{featureData}, an \Rpackage{AnnotatedDataFrame} object, where rows are features (e.g. genes), and columns are gene attributes, such as biotype, gc content, etc. \end{enumerate} The expression value matrix \emph{must} have the same number of columns as the \Robject{phenoData} has rows, and it must have the same number of rows as the \Robject{featureData} data frame has rows. Row names of the \Robject{phenoData} object should match the column names of the expression matrix. Row names of the \Robject{featureData} object should match row names of the expression matrix. You can create a new \Rclass{CellDataSet} object as follows: <>= #not run HSMM_expr_matrix <- read.table("fpkm_matrix.txt") HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt") HSMM_gene_annotation <- read.delim("gene_annotations.txt") @ %def Once these tables are loaded, you can create the CellDataSet object like this: <>= pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet) fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation) HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix), phenoData = pd, featureData = fd) @ %def This will create a CellDataSet object with expression values measured in FPKM, a measure of relative expression reported by Cufflinks. By default, Monocle assumes that your expression data is log-normally distrubuted and uses a Tobit model to test for differential expression in downstream steps. However, if you're not using TPM or FPKM data, see below for how to tell Monocle how to model it in downstream steps \subsection{Choosing a distribution for expression data} Monocle works well with both relative expression data and count-based measures (e.g. UMIs). In general, it works best with transcript count data, especially UMI data. Whatever your data type, it is \emph{critical} that specify the appropriate distribution for it. FPKM/TPM values are generally log-normally distributed, while UMIs or read counts are better modeled with the negative binomial. To work with count data, specify the negative binomial distribution as the \Robject{expressionFamily} argument to \Rfunction{newCellDataSet}: <>= # Not run HSMM <- newCellDataSet(count_matrix, phenoData = pd, featureData = fd, expressionFamily=negbinomial()) @ %def There are several allowed values for \Robject{expressionFamily}, which expects a ``family function'' from the \Rpackage{VGAM} package: \newline \newline \begin{tabular}{l p{1.5in} p{3.25in}} \textbf{Family function} & \textbf{Data type} & \textbf{Notes} \\ \hline \noalign{\smallskip} \Rfunction{tobit()} (default) & FPKM, TPM & Tobits are truncated normal distributions. Using \Rfunction{tobit()} will tell Monocle to log-transform your data where appropriate. Do not transform it yourself. \\ \hline \noalign{\smallskip} \Rfunction{negbinomial()} & UMIs, Transcript counts from experiments with spike-ins or \Rfunction{relative2abs}, raw read counts & Using \Rfunction{negbinomial()} can be slow for very large datasets. In these cases, consider \Rfunction{negbinomial.size()}. \\ \hline \noalign{\smallskip} \Rfunction{negbinomial.size()} & UMIs, Transcript counts from experiments with spike-ins, raw read counts & Slightly less accurate for differential expression than \Rfunction{negbinomial()}, but much, much faster. \\ \hline \noalign{\smallskip} \Rfunction{gaussianff()} & log-transformed FPKM/TPMs, Ct values from single-cell qPCR & If you want to use Monocle on data you have already transformed to be normally distributed, you can use this function, though some Monocle features may not work well. \end{tabular} \textbf{Using the wrong expressionFamily for your data will lead to bad results}, errors from Monocle, or both. However, if you have FPKM/TPM data, you can still use negative binomial if you first convert your relative expression values to transcript counts using \Rfunction{relative2abs}. This often leads to much more accurate results than using \Rfunction{tobit()}. See \ref{conv_rel2abs} for details. \subsection{Working with large data sets} Some single-cell RNA-Seq experiments report measurements from tens of thousands of cells or more. As instrumentation improves and costs drop, experiments will become ever larger and more complex, with many conditions, controls, and replicates. A matrix of expression data with 50,000 cells and a measurement for each of the 25,000+ genes in the human genome can take up a lot of memory. However, because current protocols typically don't capture all or even most of the mRNA molecules in each cell, many of the entries of expression matrices are zero. Using \emph{sparse matrices} can help you work with huge datasets on a typical computer. To work with your data in a sparse format, simply provide it to Monocle as a sparse matrix from the \Rpackage{Matrix} package: <>= HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"), phenoData = pd, featureData = fd, lowerDetectionLimit=1, expressionFamily=negbinomial.size()) @ %def Other sparse matrix packages, such as \Rpackage{slam} or \Rpackage{SparseM} are not supported. \subsection{Converting relative expression values into mRNA counts} \label{conv_rel2abs} If you performed your single-cell RNA-Seq experiment using \emph{spike-in} standards, you can convert these measurements into mRNAs per cell (RPC). RPC values are often easier to analyze than FPKM or TPM values, because have better statistical tools to model them. In fact, it's possible to convert FPKM or TPM values to RPC values even if there were no spike-in standards included in the experiment. Monocle 2 includes an algorithm called \emph{Census} which performs this conversion (Qiu \emph{et al}, submitted). You can convert to RPC values before creating your CellDataSet object using the \Robject{relative2abs} function, as follows: <>= pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet) fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation) # First create a CellDataSet from the relative expression levels HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix), phenoData = pd, featureData = fd) # Next, use it to estimate RNA counts rpc_matrix <- relative2abs(HSMM) # Now, make a new CellDataSet using the RNA counts HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"), phenoData = pd, featureData = fd, lowerDetectionLimit=1, expressionFamily=negbinomial.size()) @ %def Note that since we are using RPC values, we have changed the value of \Robject{lowerDetectionLimit} to reflect the new scale of expression. Importantly, we have also set the \Robject{expressionFamily} to \Robject{negbinomial()}, which tells Monocle to use the negative binomial distribution in certain downstream statistical tests. Failing to change these two options can create problems later on, so make sure not to forget them when using RPC values. Finally, we'll also call two functions that pre-calculate some information about the data. Size factors help us normalize for differences in mRNA recovered across cells, and ``dispersion'' values will help us perform differential expression analysis later. <>= HSMM <- estimateSizeFactors(HSMM) HSMM <- estimateDispersions(HSMM) @ % We're now ready to start using the \Robject{HSMM} object in our analysis. \subsection{Filtering low-quality cells} The first step in any single-cell RNA-Seq analysis is identifying poor-quality libraries from further analysis. Most single-cell workflows will include at least some libraries made from dead cells or empty wells in a plate. It's also crucial to remove doublets: libraries that were made from two or more cells accidentally. These cells can disrupt downstream steps like pseudotime ordering or clustering. This section walks through typical quality control steps that should be performed as part of all analyses with Monocle. It is often convenient to know how many express a particular gene, or how many genes are expressed by a given cell. Monocle provides a simple function to compute those statistics: <>= HSMM <- detectGenes(HSMM, min_expr = 0.1) print(head(fData(HSMM))) expressed_genes <- row.names(subset(fData(HSMM), num_cells_expressed >= 10)) @ %def The vector \Robject{expressed\_genes} now holds the identifiers for genes expressed in at least 50 cells of the data set. We will use this list later when we put the cells in order of biological progress. It is also sometimes convenient to exclude genes expressed in few if any cells from the \Rclass{CellDataSet} object so as not to waste CPU time analyzing them for differential expression. Let's start trying to remove unwanted, poor quality libraries from the CellDataSet. Your single cell RNA-Seq protocol may have given you the opportunity to image individual cells after capture but prior to lysis. This image data allows you to score your cells, confirming that none of your libraries were made from empty wells or wells with excess cell debris. With some protocols and instruments, you may get more than one cell captured instead just a single cell. You should exclude libraries that you believe did not come from a single cell, if possible. Empty well or debris well libraries can be especially problematic for Monocle. It's also a good idea to check that each cell's RNA-seq library was sequenced to an acceptible degree. While there is no widely accepted minimum level for what constitutes seequencing ``deeply enough'', use your judgement: a cell sequenced with only a few thousand reads is unlikely to yield meaningful measurements. \Rclass{CellDataSet} objects provide a convenient place to store per-cell scoring data: the \Robject{phenoData} slot. Simply include scoring attributes as columns in the data frome you used to create your \Rclass{CellDataSet} container. You can then easily filter out cells that don't pass quality control. You might also filter cells based on metrics from high throughput sequencing quality assessment packages such as FastQC. Such tools can often identify RNA-Seq libraries made from heavily degraded RNA, or where the library contains an abnormally large amount of ribosomal, mitochondrial, or other RNA type that you might not be interested in. The HSMM dataset included with this package has scoring columns built in: <>= print(head(pData(HSMM))) @ %def This dataset has already been filtered using the following commands: <>= valid_cells <- row.names(subset(pData(HSMM), Cells.in.Well == 1 & Control == FALSE & Clump == FALSE & Debris == FALSE & Mapped.Fragments > 1000000)) HSMM <- HSMM[,valid_cells] @ %def If you are using RPC values to measure expresion, as we are in this vignette, it's also good to look at the distribution of mRNA totals across the cells: <>= pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM)) HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6] upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) + 2*sd(log10(pData(HSMM)$Total_mRNAs))) lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) - 2*sd(log10(pData(HSMM)$Total_mRNAs))) qplot(Total_mRNAs, data=pData(HSMM), color=Hours, geom="density") + geom_vline(xintercept=lower_bound) + geom_vline(xintercept=upper_bound) HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound & pData(HSMM)$Total_mRNAs < upper_bound] HSMM <- detectGenes(HSMM, min_expr = 0.1) @ %def We've gone ahead and removed the few cells with either very low mRNA recovery or far more mRNA that the typical cell. Often, doublets or triplets have roughly twice the mRNA recovered as true single cells, so the latter filter is another means of excluding all but single cells from the analysis. Such filtering is handy if your protocol doesn't allow directly visualization of cell after they've been captured. Note that these thresholds of 10,000 and 40,000 mRNAs are specific to this dataset. You may need to adjust filter thresholds based on your experimental protocol. Once you've excluded cells that do not pass your quality control filters, you should verify that the expression values stored in your \Rclass{CellDataSet} follow a distribution that is roughly lognormal: <>= # Log-transform each value in the expression matrix. L <- log(exprs(HSMM[expressed_genes,])) # Standardize each gene, so that they are all on the same scale, # Then melt the data with plyr so we can plot it easily" melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L)))) # Plot the distribution of the standardized gene expression values. qplot(value, geom="density", data=melted_dens_df) + stat_function(fun = dnorm, size=0.5, color='red') + xlab("Standardized log(FPKM)") + ylab("Density") @ %def \section{Classifying and counting cells of different types} Single cell experiments are often performed on complex mixtures of multiple cell types. Dissociated tissue samples might contain two, three, or even many different cells types. In such cases, it's often nice to classify cells based on type using known markers. In the myoblast experiment, the culture contains fibroblasts that came from the original muscle biopsy used to establish the primary cell culture. Myoblasts express some key genes that fibroblasts don't. Selecting only the genes that express, for example, sufficiently high levels of \emph{MYF5} excludes the fibroblasts. Likewise, fibroblasts express high levels of \emph{ANPEP} (CD13), while myoblasts tend to express few if any transcripts of this gene. \subsection{Classifying cells with CellTypeHierarchy} Monocle provides a simple system for tagging cells based on the expression of marker genes of your choosing. You simply provide a set of functions that Monocle can use to annotate each cell. For example, you could provide a function for each of several cell types. These functions accept as input the expression data for each cell, and return TRUE to tell Monocle that a cell meets the criteria defined by the function. So you could have one function that returns TRUE for cells that express myoblast-specific genes, another function for fibroblast-specific genes, etc. Here's an example of such a set of ``gating'' functions: <>= MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5")) ANPEP_id <- row.names(subset(fData(HSMM), gene_short_name == "ANPEP")) cth <- newCellTypeHierarchy() cth <- addCellType(cth, "Myoblast", classify_func=function(x) {x[MYF5_id,] >= 1}) cth <- addCellType(cth, "Fibroblast", classify_func=function(x) {x[MYF5_id,] < 1 & x[ANPEP_id,] > 1}) @ %def The functions are organized into a small data structure called a \Robject{CellTypeHierarchy}, that Monocle uses to classify the cells. You first initialize a new CellTypeHierarchy object, then register your gating functions within it. Once the data strucure is set up, you can use it to classify all the cells in the experiment: <>= HSMM <- classifyCells(HSMM, cth, 0.1) @ %def The function \Rfunction{classifyCells} applies each gating function to each cell, classifies the cells according to the gating functions, and returns the \Robject{CellDataSet} with a new column, \Robject{CellType} in its \Robject{pData} table. We can now count how many cells of each type there are in the experiment. <>= table(pData(HSMM)$CellType) pie <- ggplot(pData(HSMM), aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1) pie + coord_polar(theta = "y") + theme(axis.title.x=element_blank(), axis.title.y=element_blank()) @ %def Note that many cells are marked ``Unknown''. This is common, largely because of the low rate of mRNA capture in most single-cell RNA-Seq experiments. A cell might express a few \emph{MYF5} mRNAs, but we weren't lucky enough to capture one of them. When a cell doesn't meet any of the criteria specified in your classification functions, it's marked ``Unknown''. If it meets multiple functions' criteria, it's marked ``Ambiguous''. You could just exclude such cells, but you'd be throwing out a lot of your data. In this case, we'd lose more than half of the cells! \subsection{Unsupervised cell clustering} Fortunately, Monocle provides an algorithm you can use to impute the types of the ``Unknown'' cells. This algorithms, implemented in the function \Rfunction{clusterCells}, groups cells together according to global expression profile. That way, if your cell expresses lots of genes specific to myoblasts, but just happens to lack \emph{MYF5}, we can still recognize it as a myoblast. \Rfunction{clusterCells} can be used in an unsupervised manner, as well as in a ``semi-supervised'' mode, which allows to assist the algorithm with some expert knowledge. Let's look at the unsupervised mode first. The first step is to decide which genes to use in clustering the cells. We could use all genes, but we'd be including a lot of genes that are not expressed at a high enough level to provide meaningful signal. Including them would just add noise to the system. We can filter genes based on average expression level, and we can additionally select genes that are unusually variable across cells. These genes tend to be highly informative about cell state. <>= disp_table <- dispersionTable(HSMM) unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit) HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id) plot_ordering_genes(HSMM) @ %def The \Rfunction{setOrderingFilter} function marks genes that will be used for clustering in subsequent calls to \Rfunction{clusterCells}, although we will be able to provide other lists of genes if we want. \Rfunction{plot\_ordering\_genes} shows how variability (dispersion) in a gene's expression depends on the average expression across cells. The red line shows Monocle's expectation of the dispersion based on this relationship. The genes we marked for use in clustering are shown as black dots, while the others are shown as grey dots. Now we're ready to try clustering the cells: <>= HSMM <- clusterCells(HSMM, num_clusters=2) plot_cell_trajectory(HSMM, 1, 2, color="CellType") @ %def As you can see, the cells fall into two different clusters. The cells tagged as myoblasts by our gating functions are marked in green, while the fibroblasts are tagged in red. The cells that don't express either marker are blue. Unfortunately, the cells don't cluster by type. This isn't all that surprising, because myoblasts and contaminating interstitial fibroblasts express many of the same genes in these culture conditions. Moroever, there are other sources of variation in the experiment that might be driving the clustering. One source of variation in the experiment stems from the experimental design. To initiate myoblast differentiation, we switch media from a high-mitogen growth medium (GM) to a low-mitogen differentiation medium (DM). Perhaps the cells are clustering based on the media they're cultured in? <>= plot_cell_trajectory(HSMM, 1, 2, color="Media") @ %def Fortunately, Monocle allows us to subtract the effects of ``uninteresting'' sources of variation so they don't impact the clustering. You can do this with the \Robject{residualModelFormulaStr} argument to \Rfunction{clusterCells} and several other Monocle functions. This argument accepts an R model formula string specifying the effects you want to subtract prior to clustering. <>= HSMM <- clusterCells(HSMM, residualModelFormulaStr="~Media", num_clusters=2) plot_cell_trajectory(HSMM, 1, 2, color="Media") @ %def Now the culture media appears to play less of a role in the clustering. Another common \emph{technical} source of variation is the efficiency of the library construction reaction (particularly the cDNA synthesis steps). When library prep goes particularly well for a cell, we'll recover more mRNA and detect more genes from it. We can control for this effect and the media effect at the same time by expanding the residual model: <>= HSMM <- clusterCells(HSMM, residualModelFormulaStr="~Media + num_genes_expressed", num_clusters=2) plot_cell_trajectory(HSMM, 1, 2, color="num_genes_expressed") @ %def Now that we've accounted for some unwanted sources of variation, we're ready to take another crack at classifying the cells by unsupervised clustering: <>= plot_cell_trajectory(HSMM, 1, 2, color="CellType") @ %def Now, most of the myoblasts are in one cluster, most of the fibroblasts are in the other, and the unknowns are spread across both. However, we still see some cells of both types in each cluster. This could be due to lack of specificity in our marker genes and our CellTypeHierarchy functions, but it could also be due to suboptimal clustering. To help rule out the latter, let's try running \Rfunction{clusterCells} in its semi-supervised mode. \subsection{Semi-supervised cell clustering with known marker genes} First, we'll select a different set of genes to use for clustering the cells. Before we just picked genes that were highly expressed and highly variable. Now, we'll pick genes that \emph{co-vary} with our markers. In a sense, we'll be building a large list of genes to use as markers, so that even if a cell doesn't have \emph{MYF5}, it might be recognizable as a myoblast based on other genes. <>= marker_diff <- markerDiffTable(HSMM[expressed_genes,], cth, residualModelFormulaStr="~Media", cores=1) @ %def The function \Rfunction{markerDiffTable} takes a CellDataSet and a CellTypeHierarchy and classifies all the cells into types according to your provided functions. It then removes all the ``Unknown'' and ``Ambiguous'' functions before identifying genes that are differentially expressed between the types. Again, you can provide a residual model of effects to exclude from this test. The function then returns a data frame of test results, and you can use this to pick the genes you want to use for clustering. Often it's best to pick the top 10 or 20 genes that are most specific for each cell type. This ensures that the clustering genes aren't dominated by markers for one cell type. You generally want a balanced panel of markers for each type if possible. Monocle provides a handy function for ranking genes by how restricted their expression is for each type. <>= candidate_clustering_genes <- row.names(subset(marker_diff, qval < 0.05)) marker_spec <- calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth) head(selectTopMarkers(marker_spec, 3)) @ %def The last line above shows the top three marker genes for myoblasts and fibroblasts. The "specificity" score is calculated using the metric described in Cabili et al \cite{CABILI} and can range from zero to one. The closer it is to one, the more restricted it is to the cell type in question. You can use this feature to define new markers for known cell types, or pick out genes you can use in purifying newly discovered cell types. This can be highly valuable for downstream follow up experiments. To cluster the cells, we'll choose the top ten markers for each of these cell types: <>= semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 20)$gene_id) HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes) plot_ordering_genes(HSMM) @ %def Note that we've got a smaller set of genes, and some of them are not especially highly expressed or variable across the experiment. However, they are great for distinguishing cells that express \emph{MYF5} from those that have \emph{ANPEP}. We've already marked them for use in clustering, but even if we hadn't, we could still use them by providing them directly to \Rfunction{clusterCells}. <>= HSMM <- clusterCells(HSMM, num_clusters=2, clustering_genes=semisup_clustering_genes, residualModelFormulaStr="~Media + num_genes_expressed") plot_cell_trajectory(HSMM, 1, 2, color="CellType") @ %def \subsection{Imputing cell type} Note that we've reduce the number of ``contaminating'' fibroblasts in the myoblast cluster, and vice versa. But what about the ``Unknown'' cells? If you provide \Rfunction{clusterCells} with a the CellTypeHierarcy, Monocle will use it classify \emph{whole clusters}, rather than just individual cells. Essentially, cluserCells works exactly as before, except after the clusters are built, it counts the frequency of each cell type in each cluster. When a cluster is composed of more than a certain percentage (in this case, 10\%) of a certain type, all the cells in the cluster are set to that type. If a cluster is composed of more than one cell type, the whole thing is marked ``Ambiguous''. If there's no cell type thats above the threshold, the cluster is marked ``Unknown''. Thus, Monocle helps you impute the type of each cell even in the presence of missing marker data. <>= HSMM <- clusterCells(HSMM, num_clusters=2, frequency_thresh=0.1, cell_type_hierarchy=cth, clustering_genes=row.names(subset(marker_diff, qval < 0.05)), residualModelFormulaStr="~Media + num_genes_expressed") plot_cell_trajectory(HSMM, 1, 2, color="CellType", markers = c("MYF5", "ANPEP")) @ %def As you can see, the clusters are highly pure in terms of \emph{MYF5} expression. There are some cells expressing \emph{ANPEP} in both clusters, but those in the myoblast cluster also express \emph{MYF5}. This is not surprising, as ANPEP isn't a very specific marker of fibroblasts. Overall, we've successfully classified all the cells: <>= pie <- ggplot(pData(HSMM), aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1) pie + coord_polar(theta = "y") + theme(axis.title.x=element_blank(), axis.title.y=element_blank()) @ %def Finally, we subset the CellDataSet object to create $HSMM_myo$, which includes only myoblasts. <>= HSMM_myo <- HSMM[,pData(HSMM)$CellType == "Myoblast"] HSMM_myo <- estimateDispersions(HSMM_myo) @ %def \section{Constructing single cell trajectories} During development, in response to stimuli, and througout life, cells transition from one functional ``state'' to another. Cells in different states express different sets of genes, producing a dynamic repetoire of proteins and metabolites that carry out their work. As cells move between states, undergo a process of transcriptional re-configuration, with some genes being silenced and others newly activated. These transient states are often hard to characterize because purifying cells in between more stable endpoint states can be difficult or impossible. Single-cell RNA-Seq can enable you to see these states without the need for purification. However, to do so, we must determine where each cell is the range of possible states. Monocle introduced the strategy of using RNA-Seq for \emph{single cell trajectory analysis}. Rather than purifying cells into discrete states experimentally, Monocle uses an algorithm to learn the sequence of gene expression changes each cell must go through as part of a dynamic biological process. Once it has learned the overall ``trajectory'' of gene expression changes, Monocle can place each cell at its proper position in the trajectory. You can then use Monocle's differential analysis toolkit to find genes regulated over the course of the trajectory, as described in section \ref{pseudotime_diff}. If there are multiple outcome for the process, Monocle will reconstruct a ``branched'' trajectory. These branches correspond to cellular ``decisions'', and Monocle provides powerful tools for identifying the genes affected by them and involved in making them. You can see how to analyze branches in section \ref{BEAM}. Monocle relies on a machine learning technique called \emph{manifold learning} to construct single-cell trajectories. You can read more about the theoretical foundations of Monocle's approach in section \ref{theory}, or consult the references shown below in section \ref{references}. \subsection{``Pseudotime'': a measure of progress through a biological process } In many biological processes, cells do not progress in perfect synchrony. In single-cell expression studies of processes such as cell differentiation, captured cells might be widely distributed in terms of progress. That is, in a population of cells captured at exactly the same time, some cells might be far along, while others might not yet even have begun the process. This asynchrony creates major problems when you want to understand the sequence of regulatory changes that occur as cells transition from one state to the next. Tracking the expression across cells captured at the same time produces a very compressed sense of a gene's kinetics, and the apparent variability of that gene's expression will be very high. By ordering each cell according to its progress along a learned trajectory, Monocle alleviates the problems that arise due to asynchrony. Instead of tracking changes in expression as a function of time, Monocle tracks changes as a function of progress along the trajectory, which we term ``pseudotime''. Pseudotime is an abstract unit of progress: it's simply the distance between a cell and the start of the trajectory, measured along the shortest path. The trajectory's total length is defined in terms of the total amount of transcriptional change that a cell undergoes as it moves from the starting state to the end state. For further details, see section \ref{theory}. \subsection{The ordering algorithm} \subsubsection{Choosing genes for ordering} Inferring a single-cell trajectory is a hard machine learning problem. The first step is to select the genes Monocle will use as input for its machine learning approach. This is called \emph{feature selection}, and it has a major impact in the shape of the trajectory. In single-cell RNA-Seq, genes expressed at low levels are often very noisy, but some may contain important information regarding the state of the cell. If we simply provide all the input data, Monocle might get confused by the noise, or fix on a feature of the data that isn't biologically meaningful, such as batch effects arising from collecting data on different days. Monocle provides you with a variety of tools to select genes that will yield a robust, accurate, and biologically meaningful trajectory. You can use these tools to either perform a completely ``unsupervised'' analysis, in which Monocle has no forehand knowledge of which gene you consider important. Alternatively, you can make use of expert knowledge in the form of genes that are already known to define biolgical progress to shape Monocle's trajectory. We consider this mode ``semi-supervised'', because Monocle will augment the markers you provide with other, related genes. Monocle then uses these genes to produce trajectories consistent with known biology but that often reveal new regulatory structure. We return to the muscle data to illustrate both of these modes. \subsubsection{Reducing the dimensionality of the data} Once we have selected the genes we will use to order the cells, Monocle applies a \emph{dimensionality reduction} to the data, which will drastically improves the quality of the trajectory. Monocle reduces the dimensionality of the data with the \Rfunction{reduceDimension} function. This function has a number of options, so you should familiarize yourself with them by consulting the its manual page. You can choose from two algorithms in \Rfunction{reduceDimension}. The first, termed Independent Component Analysis, is a classic linear technique for decomposing data that powered the original version of Monocle. The second, called DDRTree, is a much more powerful nonlinear technique that is the default for Monocle 2. For more on how these both work, see section \ref{theory}. \subsubsection{Ordering the cells in pseudotime} With the expression data projected into a lower dimensional space, Monocle is ready to learn the trajectory that describes how cells transition from one state into another. Monocle assumes that the trajectory has a tree structure, with one end of it the ``root'', and the others the ``leaves''. A cell at the beginning of the biological process starts at the root and progresses along the trunk until it reaches the first branch, if there is one. It then chooses a branch, and moves further and further along the tree until it reaches a leaf. These mathematical assumptions translate into some important biological ones. First, that the data includes all the major stages of the biological process. If your experiment failed to capture any cells at a key developmental transition, Monocle won't know its there. Second, that gene expression changes are smooth as a cell moves from one stage to the next. This assumption is realistic: major discontinuities in the trajectory would amount to a cell almost instantaneously turning over its transcriptome, which probably doesn't happen in most biolgical processes. To order your cells, Monocle uses the \Rfunction{orderCells} function. This routine has one important argument, which allows you to set the root of the tree, and thus the beginning of the process. See the manual page for \Rfunction{orderCells} for more details. \subsection{Unsupervised ordering} In this section, we discuss ordering cells in a completely unsupervised fashion. First, we must decide which genes we will use to define a cell's progress through myogenesis. Monocle orders cells by examining the pattern of expression of these genes across the cell population. Monocle looks for genes that vary in ``interesting'' (i.e. not just noisy) ways, and uses these to structure the data. We ultimately want a set of genes that increase (or decrease) in expression as a function of progress through the process we're studying. Ideally, we'd like to use as little prior knowledge of the biology of the system under study as possible. We'd like to discover the important ordering genes from the data, rather than relying on literature and textbooks, because that might introduce bias in the ordering. One effective way to isolate a set of ordering genes is to simply compare the cells collected at the beginning of the process to those at the end and find the differentially expressed genes, as described above. The command below will find all genes that are differentially expressed in response to the switch from growth medium to differentiation medium: <>= #not run diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,], fullModelFormulaStr="~Media") ordering_genes <- row.names (subset(diff_test_res, qval < 0.01)) @ %def Choosing genes based on differential analysis of time points is often highly effective, but what if we don't have time series data? If the cells are asynchronously moving through our biological process (as is usually the case), Monocle can often reconstruct their trajectory from a single population captured all at the same time. Below are two methods to select genes that require no knowledge of the design of the experiment at all. \subsubsection{Selecting genes with high dispersion across cells} Genes that vary a lot are often highly informative for identifying cell subpopulations or ordering cells along a trajectory. In RNA-Seq, a gene's variance typically depends on its mean, so we have to be a bit careful about how we select genes based on their variance. <>= disp_table <- dispersionTable(HSMM_myo) ordering_genes <- subset(disp_table, mean_expression >= 0.5 & dispersion_empirical >= 2 * dispersion_fit)$gene_id @ %def Once we have a list of gene ids to be used for ordering, we need to set them in the \Robject{HSMM} object, because the next several functions will depend on them. <>= HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes) plot_ordering_genes(HSMM_myo) @ %def The genes we've chosen to use for ordering define the \emph{state space} of the cells in our data set. Each cell is a point in this space, which has dimensionality equal to the number of genes we've chosen. So if there are 500 genes used for ordering, each cell is a point in a 500-dimensional space. For a number of reasons, Monocle works better if we can \emph{reduce} the dimensionality of that space before we try to put the cells in order. In this case, we will reduce the space down to one with two dimensions, which we will be able to easily visualize and interpret while Monocle is ordering the cells. <>= HSMM_myo <- reduceDimension(HSMM_myo, max_components=2) @ %def Now that the space is reduced, it's time to order the cells using the \Rfunction{orderCells} function as shown below. The second argument is the \Robject{reverse} flag. Monocle won't be able to tell without some help which cells are at the beginning of the process and which are at the end. The \Robject{reverse} flag tells Monocle to reverse the orientation of the entire process as it's being discovered from the data, so that the cells that would have been assigned to the end are instead assigned to the beginning, and so on. Setting \Robject{reverse} to true will reverse the ordering of the cells in pseudotime, which also has a significant impact on the orientation of branches in the trajectory associated with cell fate decisions. <>= HSMM_myo <- orderCells(HSMM_myo, reverse=FALSE) @ %def Once the cells are ordered, we can visualize the trajectory in the reducted dimesional space. <>= plot_cell_trajectory(HSMM_myo, color_by="Hours") @ %def To confirm that the ordering is correct, and to verify that we don't need to flip around the ordering using the \Robject{reverse} flag in \Robject{orderCells}, we can select a couple of markers of myogenic progress. Plotting these genes demonstrates that ordering looks good: <>= HSMM_expressed_genes <- row.names(subset(fData(HSMM_myo), num_cells_expressed >= 10)) HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,] my_genes <- row.names(subset(fData(HSMM_filtered), gene_short_name %in% c("CDK1", "MEF2C", "MYH3"))) cds_subset <- HSMM_filtered[my_genes,] plot_genes_in_pseudotime(cds_subset, color_by="Hours") @ %def \subsubsection{Selecting genes based on PCA loading} A number of single-cell clustering studies have found that principal component analysis (PCA) is an effective way of finding genes that vary widely across cells. You can use PCA to pick out a good set of ordering genes in a completely unsupervised manner. The procedure below first normalizes the expression data by adding a pseudocount, log-transforming it, and then standardizing it so that all genes have roughly the same dyanmic range. Then, it uses the \Rpackage{irlba} package to perform the PCA. The code below is a little more complicated than just calling \Rfunction{prcomp()}, but it has the advantage of performing the entire computation with sparse matrix operations. This is important for large, sparse CellDataSet objects like the ones common in droplet-based single-cell RNA-Seq. <>= exprs_filtered <- t(t(exprs(HSMM_filtered)/pData(HSMM_filtered)$Size_Factor)) nz_genes <- which(exprs_filtered != 0) exprs_filtered[nz_genes] <- log(exprs_filtered[nz_genes] + 1) # Calculate the variance across genes without converting to a dense # matrix: expression_means <- Matrix::rowMeans(exprs_filtered) expression_vars <- Matrix::rowMeans((exprs_filtered - expression_means)^2) # Filter out genes that are constant across all cells: genes_to_keep <- expression_vars > 0 exprs_filtered <- exprs_filtered[genes_to_keep,] expression_means <- expression_means[genes_to_keep] expression_vars <- expression_vars[genes_to_keep] # Here's how to take the top PCA loading genes, but using # sparseMatrix operations the whole time, using irlba. irlba_pca_res <- irlba(t(exprs_filtered), nu=0, center=expression_means, scale=sqrt(expression_vars), right_only=TRUE)$v row.names(irlba_pca_res) <- row.names(exprs_filtered) # Take the top 100 genes from components 2 and 3. Component # 1 usually is driven by technical noise. PC2_genes <- names(sort(abs(irlba_pca_res[, 2]), decreasing = T))[1:100] PC3_genes <- names(sort(abs(irlba_pca_res[, 3]), decreasing = T))[1:100] ordering_genes <- union(PC2_genes, PC3_genes) @ %def Using these to order the cells as above yields the following trajectory: <>= HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes) HSMM_myo <- reduceDimension(HSMM_myo, max_components=2) HSMM_myo <- orderCells(HSMM_myo, reverse=FALSE) plot_cell_trajectory(HSMM_myo, color_by="Hours") @ %def \subsection{Semi-supervised ordering with known marker genes} Unsupervised ordering is desirable because it avoids introducing bias into the analysis. However, unsupervised machine learning will often fix on a strong feature of the data that's not the focus of your experiment. For example, where each cell is in the cell cycle has a major impact on the shape of the trajectory when you use unsupervised learning. But what if you wish to focus on cycle-independent effects in your biological process? Monocle's ``semi-supervised'' ordering mode can help you focus on the aspects of the process you're interested in. Ordering your cells in a semi-supervised manner is very simple. You first define genes that mark progress using the \Robject{CellTypeHierchy} system, very similar to how we used it for cell type classification. Then, you use it to select ordering genes that co-vary with these markers. Finally, you order the cell based on these genes just as we do in unsupervised ordering. So the only difference between unsupervised and semi-supervised ordering is in which genes we use for ordering. As we saw before, myoblasts begin differnentation by exiting the cell cycle and then proceed through a sequence of regulatory events that leads to expression of some key muscle-specific proteins needed for contraction. We can mark cycling cells with cyclin B2 (\emph{CCNB2}) and recognize myotubes as those cells expressed high levels of myosin heavy chain 3 (\emph{MYH3}). <>= CCNB2_id <- row.names(subset(fData(HSMM), gene_short_name == "CCNB2")) MYH3_id <- row.names(subset(fData(HSMM), gene_short_name == "MYH3")) cth <- newCellTypeHierarchy() cth <- addCellType(cth, "Cycling myoblast", classify_func=function(x) {x[CCNB2_id,] >= 1}) cth <- addCellType(cth, "Myotube", classify_func=function(x) {x[MYH3_id,] >=1}) @ %def Now we select the set of genes that co-vary (in either direction) with these two ``bellweather'' genes: <>= marker_diff <- markerDiffTable(HSMM_myo[expressed_genes,], cth, cores=1) semisup_clustering_genes <- row.names(subset(marker_diff, qval < 0.01)) length(semisup_clustering_genes) @ %def Using the 705 genes for ordering produces a trajectory that's highly similar to the one we obtained with unsupervised methods, but it's a little ``cleaner'': the small branches are gone. <>= HSMM_myo <- setOrderingFilter(HSMM_myo, semisup_clustering_genes) plot_ordering_genes(HSMM_myo) HSMM_myo <- reduceDimension(HSMM_myo, max_components=2) HSMM_myo <- orderCells(HSMM_myo, reverse=TRUE) plot_cell_trajectory(HSMM_myo, color_by="Hours") @ %def And we can see that kinetics of notable genes we examined before are largely the same. <>= HSMM_filtered <- HSMM_myo[expressed_genes,] my_genes <- row.names(subset(fData(HSMM_filtered), gene_short_name %in% c("CDK1", "MEF2C", "MYOG"))) cds_subset <- HSMM_filtered[my_genes,] plot_genes_in_pseudotime(cds_subset, color_by="Hours") @ %def \subsection{Reconstructing branched trajectories} Monocle reconstructs a linear trajectory because muscle differentiation is a relatively simple process with just one outcome for the cells. It can't tell a-priori which side of the trajectory is the start and which is the end, but that's easily solved with the \Robject{reverse} flag. However, if the trajectory includes multiple outcomes, Monocle might need a little more help from you. Consider the experiment performed by Steve Quake's lab by Barbara Treutlein and colleages, who captured cells from the developing mouse lung. They captured cells early in development, later when the lung contains both major types of epithelial cells (AT1 and AT2), and cells right about to make the decision to become either AT1 or AT2. Monocle can reconstruct this process as a \emph{branched} trajectory, allow you to analyze the decision point in great detail. We'll learn more about branch analysis in section \ref{BEAM}. Monocle reconstructs branched trajectories, but it doesn't necessarily know which part of of the trajectory is the beginning. It's important to tell Monocle where the trajectory starts so that the Pseudotime axis runs in the right direction and reflects the underlying biological process. In the figure below, you can see that Monocle reconstructs a tree with two outcomes from the Truetlein data. The ``root'' of the trajectory happens to the branch where most of the cells collected early in the experiment (E14.5 and E16.5) can be found. So Monocle happened to start the trajectory at the right spot in this case. <>= lung <- load_lung() plot_cell_trajectory(lung, color_by="Time") plot_cell_trajectory(lung, color_by="Pseudotime") @ %def If however, we want to change the root of the tree, we first plot the trajectory coloring the cells by ``State'' so we can see our options for where we can set the root. ``State'' is just Monocle's jargon for the different segments of the tree. This tree has 3 States, but more complex trees could have quite a few more. We then call \Rfunction{order\_cells()} \emph{again}, passing this state as the \Robject{root\_state} argument: <>= plot_cell_trajectory(lung, color_by="State") lung <- orderCells(lung, root_state=3) plot_cell_trajectory(lung, color_by="Pseudotime") @ %def If there are a ton of states in your tree, it can be a little hard to make out where each one falls on the tree. Sometimes it can be handy to ``facet'' the trajectory plot so it's easier to see where each of the states are located: <>= plot_cell_trajectory(lung, color_by="State") + facet_wrap(~State) @ %def And if you don't have a timeseries, you might need to set the root based on where certain marker genes are expressed, using your biological knowledge of the system. For example, in this experiment, a highly proliferative population of progenitor cells are generating two types of post-mitotic cells. So the root should have cells that express high levels of proliferation markers. We can use the jitter plot to pick figure out which state corresponds to rapid proliferation: <>= lung_genes <- row.names(subset(fData(lung), gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn"))) plot_genes_jitter(lung[lung_genes,], grouping="State") @ %def Here we can see that State 1 is the proliferative state on the basis high expression of Cyclin B2. \section{Differential expression analysis} Differential gene expression analysis is a common task in RNA-Seq experiments. Monocle can help you find genes that are differentially expressed between groups of cells and assesses the statistical signficance of those changes. These comparisons require that you have a way to collect your cells into two or more groups. These groups are defined by columns in the \Robject{phenoData} table of each \Robject{CellDataSet}. Monocle will assess the signficance of each gene's expression level across the different groups of cells. \subsection{Basic differential analysis} Performing differential expression analysis on all genes in the human genome can take a substantial amount of time. For a dataset as large as the myoblast data from \cite{TRAPNELL_CACCHIARELLI}, which contains several hundred cells, the analysis can take several hours on a single CPU. Let's select a small set of genes that we know are important in myogenesis to demonstrate Monocle's capabilities: <>= marker_genes <- row.names(subset(fData(HSMM_myo), gene_short_name %in% c("MEF2C", "MEF2D", "MYF5", "ANPEP", "PDGFRA","MYOG", "TPM1", "TPM2", "MYH2", "MYH3", "NCAM1", "TNNT1", "TNNT2", "TNNC1", "CDK1", "CDK2", "CCNB1", "CCNB2", "CCND1", "CCNA1", "ID1"))) @ %def In the myoblast data, the cells collected at the outset of the experiment were cultured in ``growth medium'' (GM) to prevent them from differentiating. After they were harvested, the rest of the cells were switched over to ``differentiation medium'' (DM) to promote differentiation. Let's have monocle find which of the genes above are affected by this switch: <>= diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,], fullModelFormulaStr="~Media") # Select genes that are significant at an FDR < 10% sig_genes <- subset(diff_test_res, qval < 0.1) sig_genes[,c("gene_short_name", "pval", "qval")] @ %def So 16 of the 22 genes are significant at a 10\% false discovery rate! This isn't surprising, as most of the above genes are highly relevant in myogenesis. Monocle also provides some easy ways to plot the expression of a small set of genes grouped by the factors you use during differential analysis. This helps you visualize the differences revealed by the tests above. One type of plot is a ``jitter'' plot. <>= MYOG_ID1 <- HSMM_myo[row.names(subset(fData(HSMM_myo), gene_short_name %in% c("MYOG", "CCNB2"))),] plot_genes_jitter(MYOG_ID1, grouping="Media", ncol=2) @ %def Note that we can control how to layout the genes in the plot by specifying the number of rows and columns. See the man page on \Rfunction{plot\_genes\_jitter} for more details on controlling its layout. Most if not all of Monocle's plotting routines return a plot object from the \Rpackage{ggplot2}. This package uses a grammar of graphics to control various aspects of a plot, and makes it easy to customize how your data is presented. See the \Rpackage{ggplot2} book\cite{Wickham} for more details. In this section, we'll explore how to use Monocle to find genes that are differentially expressed according to several different criteria. First, we'll look at how to use our previous classification of the cells by type to find genes that distinguish fibroblasts and myoblasts. Second, we'll look at how to find genes that are differentially expressed as a function of pseudotime, such as those that become activated or repressed during differentiation. Finally, you'll see how to perform multi-factorial differential analysis, which can help subtract the effects of confounding variables in your experiment. To keep the vignette simple and fast, we'll be working with small sets of genes. Rest assured, however, that Monocle can analyze many thousands of genes even in large experiments, making it useful for discovering dynamically regulated genes during the biological process you're studying. \subsection{Finding genes that distinguish cell type or state} During a dynamic biological process such as differentiation, cells might assume distinct intermediate or final states. Recall that earlier we distinguished myoblasts from contaminating fibroblasts on the basis of several key markers. Let's look at several other genes that should distinguish between fibroblasts and myoblasts. <>= to_be_tested <- row.names(subset(fData(HSMM), gene_short_name %in% c("UBC", "NCAM1", "ANPEP"))) cds_subset <- HSMM[to_be_tested,] @ %def To test the effects of \Robject{CellType} on gene expression, we simply call \Rfunction{differentialGeneTest} on the genes we've selected. However, we have to specify a \emph{model formula} in the call to tell Monocle that we care about genes with expression levels that depends on \emph{CellType}. Monocle's differential expression analysis works essentially by fitting two models to the expression values for each gene, working through each gene independently. The first of the two models is called the \emph{full} model. This model is essentially a way of predicting the expression value of each gene in a given cell knowing only whether that cell is a fibroblast or a myoblast. The second model, called the \emph{reduced} model, does the same thing, but it doesn't know the \Robject{CellType} for each cell. It has to come up with a reasonable prediction of the expression value for the gene that will be used for \emph{all} the cells. Because the full model has more information about each cell, it will do a better job of predicting the expression of the gene in each cell. The question Monocle must answer for each gene is \emph{how much better} the full model's prediction is than the reduced model's. The greater the improvement that comes from knowing the \Robject{CellType} of each cell, the more significant the differential expression result. This is a common strategy in differential analysis, and we leave a detailed statistical exposition of such methods to others. To set up the test based on \Robject{CellType}, we simply call \Rfunction{differentialGeneTest} with a string specifying \Robject{fullModelFormulaStr}. We don't have to specify the reduced model in this case, because the default of \Robject{~1} is what we want here. <>= diff_test_res <- differentialGeneTest(cds_subset, fullModelFormulaStr="~CellType") diff_test_res[,c("gene_short_name", "pval", "qval")] @ %def Note that all the genes are significantly differentially expressed as a function of \Robject{CellType} except the housekeeping gene TBP, which we're using a negative control. However, we don't know which genes correspond to myoblast-specific genes (those more highly expressed in myoblasts versus fibroblast specific genes. We can again plot them with a jitter plot to see: <>= plot_genes_jitter(cds_subset, grouping="CellType", color_by="CellType", nrow=1, ncol=NULL, plot_trend=TRUE) @ %def Note that we could also simply compute summary statistics such as mean or median expression level on a per-\Robject{CellType} basis to see this, which might be handy if we are looking at more than a handful of genes. Of course, we could test for genes that change as a function of \Robject{Hours} to find time-varying genes, or \Robject{Media} to identify genes that are responsive to the serum switch. In general, model formulae can contain terms in the pData table of the CellDataSet. The \Rfunction{differentialGeneTest} function is actually quite simple ``under the hood''. The call above is equivalent to: <>= full_model_fits <- fitModel(cds_subset, modelFormulaStr="~CellType") reduced_model_fits <- fitModel(cds_subset, modelFormulaStr="~1") diff_test_res <- compareModels(full_model_fits, reduced_model_fits) diff_test_res @ %def Occassionally, as we'll see later, it's useful to be able to call \Rfunction{fitModel} directly. \subsection{Finding genes that change as a function of pseudotime} \label{pseudotime_diff} Monocle's main job is to put cells in order of progress through a biological process (such as cell differentiation) without knowing which genes to look at ahead of time. Once it's done so, you can analyze the cells to find genes that changes as the cells make progress. For example, you can find genes that are significantly upregulated as the cells ``mature''. Let's look at a panel of genes important for myogenesis: <>= to_be_tested <- row.names(subset(fData(HSMM), gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1"))) cds_subset <- HSMM_myo[to_be_tested,] @ %def Again, we'll need to specify the model we want to use for differential analysis. This model will be a bit more complicated than the one we used to look at the differences between \Robject{CellType}. Monocle assigns each cell a ``pseudotime'' value, which records its progress through the process in the experiment. The model can test against changes as a function of this value. Monocle uses the \Rpackage{VGAM} package to model a gene's expression level as a smooth, nonlinear function of pseudotime: <>= diff_test_res <- differentialGeneTest(cds_subset, fullModelFormulaStr="~sm.ns(Pseudotime)") @ %def The \Rfunction{sm.ns} function states that Monocle should fit a natural spline through the expression values to help it describe the changes in expression as a function of progress. We'll see what this trend looks like in just a moment. Other smoothing functions are available. Once again, let's add in the gene annotations so it's easy to see which genes are significant. <>= diff_test_res[,c("gene_short_name", "pval", "qval")] @ %def We can plot the expression levels of these genes, all of which show significant changes as a function of differentiation, using the function \Rfunction{plot\_genes\_in\_pseudotime}. This function has a number of cosmetic options you can use to control the layout and appearance of your plot. <>= plot_genes_in_pseudotime(cds_subset, color_by="Hours") @ %def \subsection{Clustering genes by pseudotemporal expression pattern} A common question that arises when studying time-series gene expression studies is: ``which genes follow similar kinetic trends''? Monocle can help you answer this question by grouping genes that have similar trends, so you can analyze these groups to see what they have in common. Monocle provides a convenient way to visualize all pseudotime-dependent genes. The function \Rfunction{plot\_pseudotime\_heatmap} takes a CellDataSet object (usually containing a only subset of significant genes) and generates smooth expression curves much like \Rfunction{plot\_genes\_in\_pseudotime}. Then, it clusters these genes and plots them using the \Rpackage{pheatmap} package. This allows you to visualize modules of genes that co-vary across pseudotime. <>= diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,], fullModelFormulaStr="~sm.ns(Pseudotime)") sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1)) plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,], num_clusters = 2, cores = 1, show_rownames = T) @ %def \subsection{Multi-factorial differential expression analysis} Monocle can perform differential analysis in the presence of multiple factors, which can help you subtract some factors to see the effects of others. In the simple example below, Monocle tests three genes for differential expression between myoblasts and fibroblasts, while subtracting the effect of \Robject{Hours}, which encodes the day on which each cell was collected. To do this, we must specify both the full model and the reduced model. The full model captures the effects of both \Robject{CellType} and \Robject{Hours}, while the reduced model only knows about \Robject{Hours}. When we plot the expression levels of these genes, we can modify the resulting object returned by \Rfunction{plot\_genes\_jitter} to allow them to have independent y-axis ranges, to better highlight the differences between cell states. <>= to_be_tested <- row.names(subset(fData(HSMM), gene_short_name %in% c("TPM1", "MYH3", "CCNB2", "GAPDH"))) cds_subset <- HSMM[to_be_tested,] diff_test_res <- differentialGeneTest(cds_subset, fullModelFormulaStr="~CellType + Hours", reducedModelFormulaStr="~Hours") diff_test_res[,c("gene_short_name", "pval", "qval")] plot_genes_jitter(cds_subset, grouping="Hours", color_by="CellType", plot_trend=TRUE) + facet_wrap( ~ feature_label, scales="free_y") @ %def \section{Analyzing branches in single-cell trajectories} \label{BEAM} Often, single-cell trajectories include branches. The branches occur because cells execute alternative gene expression programs. Branches appear in trajectories during development, when cells make fate choices: one developmental lineage proceeds down one path, while the other lineage produces a second path. Monocle contains extensive functionality for analyzing these branching events. Truetlein and colleagues perfomed a single cell analysis of lung epithelial cell development that reveals a key cell fate decision. The figure below shows the trajectory Monocle reconstructs using some of their data. There is a single branch, labeled "1". What genes change as cells pass from the early developmental stage the top left of the tree through the branch? What genes are differentially expressed between the branches? To answer this question, Monocle provides you with a special statistical test: branched expression analysis modeling, or BEAM. <>= lung <- load_lung() plot_cell_trajectory(lung, color_by="Time") @ %def BEAM takes as input a CellDataSet that's been ordered with \Rfunction{orderCells} and the name of a branch point in the trajectory. It returns a table of significance scores for each gene. Genes that score significant are said to be \emph{branch-dependent} in their expression. <>= BEAM_res <- BEAM(lung, branch_point=1, cores = 1) BEAM_res <- BEAM_res[order(BEAM_res$qval),] BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")] @ %def You can visualize changes for all the genes that are significantly branch dependent using a special type of heatmap. This heatmap shows changes in both lineages at the same time. It also requires that you choose a branch point to inspect. Columns are points in pseudotime, rows are genes, and the beginning of pseudotime is in the middle of the heatmap. As you read from the middle of the heatmap to the right, you are following one lineage through pseudotime. As you read left, the other. The genes are clustered hierarchically, so you can visualize modules of genes that have similar lineage-dependent expression patterns. <>= plot_genes_branched_heatmap(lung[row.names(subset(BEAM_res, qval < 1e-4)),], branch_point = 1, num_clusters = 4, cores = 1, use_gene_short_name = T, show_rownames = T) @ %def We can plot a couple of these genes, such as \emph{Pdpn} and \emph{Sftpb} (both known markers of fate in this system), using the \Rfunction{plot\_genes\_branched\_pseudotime} function, which works a lot like the \Rfunction{plot\_genes\_in\_pseudotime} function, except it shows two kinetic trends, one for each lineage, instead of one. We also show \emph{Ccnd2}, a cell cycle gene, which is downregulated in both branches and is not signficant by the BEAM test. <>= lung_genes <- row.names(subset(fData(lung), gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn"))) plot_genes_branched_pseudotime(lung[lung_genes,], branch_point=1, color_by="Time", ncol=1) @ %def The \Rfunction{plot\_clusters} function returns a ggplot2 object showing the shapes of the expression patterns followed by the 100 genes we've picked out. The topographic lines highlight the distributions of the kinetic patterns relative to the overall trend lines, shown in red. %\section{Tips and tricks} \section{Major changes between Monocle version 1 and version 2} Monocle 2 is a near-complete re-write of Monocle 1. Monocle 2 is geared towards larger, more complex single-cell RNA-Seq experiments than those possible at the time Monocle 1 was written. It's also redesigned to support analysis of mRNA counts, which were hard to estimate experimentally in early versions of single-cell RNA-Seq. Now, with spike controls or UMIs, gene expression can be measured in mRNA counts. Analysis of these counts is typically easier and more accurate than relative expression values, and we encourage all users to adopt an mRNA-count centered workflow. Numerous Monocle functions have been re-written to take advantage of the nicer statistical properties of mRNA counts. For example, we adopt the dispersion modeling and variance-stabilization techniques introduced by DESeq \cite{Anders2010-oc} during differential analysis, dimensionality reduction, and other steps. Trajectory reconstruction in Monocle 2 is vastly more robust, faster, and more powerful than in Monocle 1. Monocle 2 uses an advanced nonlinear reconstruction algorithm called DDRTree\cite{Mao2015-cs}, described below in section \ref{theory}. This algorithm can expose branches that are hard to see with the less powerful linear technique used in Monocle 1. The algorithm is also far less sensitive to outliers, so careful QC and selection of high quality cells is less critical. Finally, DDRTree is much more robust in that it reports qualititatively similar trajectories more consistently when you vary the number of cells in the experiment. Although which genes are included in the ordering still greatly impact the trajectory, varying them also produces more qualititatively consistent trajectories than the previous linear technique. Because Monocle 2 is so much better at finding branches, it also includes some additional tools to help you interpret them. Branch expression analysis modeling (BEAM) is a new test for analyzing specific branch points. BEAM reports branch-dependent genes, and Monocle 2 includes some new visualization functions to help you inspect these genes. Overall, we find that branching is pervasive in diverse biological processes, and thus we expect BEAM will be very useful to those analyzing single-cell RNA-Seq data in many settings. A manuscript describing Monocle 2 and the general stragegy of using reverse graph embedding for single-cell trajectory analysis will be forthcoming (Qiu et al, in preparation). Until it appears, please continue to cite \cite{TRAPNELL_CACCHIARELLI} when you use Monocle or discuss single-cell pseudotemporal ordering in your work. \section{Theory behind Monocle}\label{theory} Single-cell expression datasets are some of the largest and most complex encountered in genomics. Even the smallest single-cell RNA-Seq experiments sample hundreds of cells, measuring the expression level of the more than 20,000 genes in each cell. Visualizing these datasets, identifying cells of different types, and comparing them to one another all pose major bioinformatics challenges. \emph{Manifold learning} is a common strategy for dealing with complex, high-dimensional data. The premise of this approach is simple: the data may reside in a very high-dimensional space, but the intrinsic structure of the dataset is much simpler. Moreover, the data are not random - they are generated by a process that can be understood by inspecting the global structure of the dataset. For example, a single single-cell RNA-Seq experiment may reside in 20,000 dimensions, but the cells might all lie on or ``near'' a curve \emph{embedded} within a much lower dimensional space. For example, we might expect that cells in different phases of the cell cycle be distributed along a closed loop. Indeed a recent large-scale single-cell RNA-Seq study found exactly that \cite{Macosko2015-fh}. Manifold learning often involves \emph{dimensionality reduction} techniques as a first step. Conventional, dimensionality reduction approaches (for example, PCA, ICA, Isomap, LLE, etc.) are limited in their ability to explictly recover the intrinisic structure from the data. Monocle 2 uses a newly developed technique, called reverse graph embedding, to simultaneously: \begin{enumerate} \item Reduce high dimensional expression data into a lower dimension space. \item Learn an explicit, smooth manifold that generates the data. \item Assign each cell to its position on that manifold \end{enumerate} Together, these tasks allow Monocle 2 to order cells in pseudotime in an entirely unsupervised, data-driven way. Importantly, Monocle 2 learns manifolds that are trees without needing any \emph{a priori} information about the structure of the tree. Users do not need to provide Monocle 2 with constraints on the number of branches, etc. These are learned from the data. This allows Monocle 2 to to discriminate between linear and branched trajectories automatically. To our knowledge, Monocle 2 is the first trajectory reconstruction algorithm to learn smooth tree-like manifolds without needing to know its high-level structure ahead of time. \subsection{Reverse graph embedding} Monocle 2 uses a technique called \emph{reverse graph embedding}\cite{Mao_undated-xi} to learn the structure of the manifold that describes a single-cell experiment. Reverse graph embedding simultaneously learns a \emph{principal graph} that approximates the manifold, as well as a function that maps points on the graph (which is embedded in low dimensions) back to the original high dimensional space. Reverse graph embedding aims to learn both a set of \emph{latent points} $\mathcal{Z} = \{\mathbf{z}_1, ..., \mathbf{z}_M\}$ corresponding to the input data that reside in the low-dimensional space along with a graph $\mathcal{G}$ that connects them. This graph approximates the manifold. In order to map points on the manifold back to the original high-dimensional input space, we also need to learn a function $f_{\mathcal{G}}$. Learning a good reverse graph embedding can be described as an optimization problem that joint captures the positions of the latent points $\mathbf{z}$, the graph $\mathcal{G}$, and the function $f_\mathcal{G}$. To learn the positions of the latent points $\mathbf{z}$, we must optimize: \begin{equation} \mathop{min}_{f_g \in \mathcal{F}} \mathop{min}_{\{\mathbf{z}_1, ..., \mathbf{z}_M\}} \sum_{i = 1}^N ||\mathbf{x}_i - f_g (\mathbf{z}_i)||^2 \end{equation} Given a set of latent point coordinates, the optimization of graph inference can be represented as: \begin{equation} \label{eq:mintree} \mathop{min}_{f_\mathcal{G} \in \mathcal{F}} \mathop{min}_{\{\mathbf{z}_1, ..., \mathbf{z}_M\}} \sum_{(V_i, V_j) \in \mathcal{E}} b_{i,j}||f_g(\mathbf{z}_i) - f_g(\mathbf{z}_j)||^2 \end{equation} where $\mathcal{X} = \{ \mathbf{x}_1, ..., \mathbf{x}_N\}$ are the original single-cell expression profiles. The $V_i$ are the the vertices of the undirected graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$. The weights for the edge in $\mathcal{E}$ are encoded as $b_{ij}$. The first optimization problem aims to position the latent points such that their image under $f_\mathcal{G}$ (that is, their corresponding positions in the high-dimensional space) will be "close" to the input data. The second optimization aims to keep latent points that are close to one another in the low dimensional space close to one another in the high dimensional space as well. These two goals must be balanced against one another. Reverse graph embedding achieves this through the parameter $\lambda$ \begin{equation} \mathop{min}_{\mathcal{G} \in \hat{\mathcal{G}}_b}\mathop{min}_{f_g \in \mathcal{F}} \mathop{min}_{\{\mathbf{z}_1, ..., \mathbf{z}_M\}} \sum_{i = 1}^N ||\mathbf{x}_i - f_g (\mathbf{z}_i)||^2 + \frac{\lambda}{2} \sum_{(V_i, V_j) \in \mathcal{E}} b_{i,j}||f_g(\mathbf{z}_i) - f_g(\mathbf{z}_j)||^2 \end{equation} Reverse graph embedding requires a feasible set $\hat{\mathcal{G}}_b$ of graph and a mapping function $f_\mathcal{G}$. In practice, implementing reverse graph embedding requires that we place some constraints on $\hat{\mathcal{G}}_b$ and $f_\mathcal{G}$. As work on reverse graph embedding continues, we anticipate that more general schemes that consider a wider range of feasible graphs and mapping functions will become available. Monocle users should expect more general reverse graph embedding schemes in future versions. Mao \emph{et al} initially described two specific ways to implement the general framework of reverse graph embedding. Both are briefly summarized below. See the original paper on DDRTree for more details. Monocle 2 uses the second scheme, but can easily be run in a mode that corresponds to the first. \subsection{DRTree: Dimensionality reduction via learning a tree} The first scheme, called ``DRTree'' described by Mao \emph{et al} learns principal graphs that are undirected trees, with one node per input data point, along with a linear function $f_\mathcal{G}$. Because the algorithm restricts the feasible set to trees, optimization problem \ref{eq:mintree} is solved by simply finding the minimum spanning tree. This can be solved quickly via Kruskal's algoritm. DRTree uses a linear projection model $f_\mathcal{G} (\mathbf{z}) = \mathbf{Wz}$ as the mapping function. The scheme optimizes: \begin{equation} \mathop{min}_{\mathbf{W}, \mathbf{Z}, \mathbf{B}} \sum_{i = 1}^N ||\mathbf{x}_i - \mathbf{W}\mathbf{z}_i||^2 + \frac{\lambda}{2} \sum_{i,j}b_{i,j}||\mathbf{W} \mathbf{z}_i - \mathbf{W} \mathbf{z}_j||^2 \end{equation} where $\mathbf{W} = [\mathbf{w}_1, ..., \mathbf{w}_d] \in \mathcal{R}^{D \times d}$ is an orthogonal set of $d$ linear basis vectors. Because several key steps of the above optimization can be solved analytically, and because the two terms can be minimized in an alternating fashion, solving it is generally very fast. However, for large input datasets, the graph can become complex, and so DRTree can run into scalability problems. \subsection{DDRTree: Discriminative dimensionality reduction via learning a tree} To overcome problems posed by large complex input datasets, Mao \emph{et al} proposed a second scheme, "DDRTree". that Monocle 2 uses instead of DRTree. Recall that the graph in DRTree contains one node per input data point. To avoid long computations and large memory footprints that come with DRTree, the second scheme learns a graph on a second, smaller set of latent points $\{\mathbf{y}_k\}_{k = 1}^K$. These points are treated by the algorithm as the centers of $\{\mathbf{z}_i\}^N_{i = 1}$. The number of these points is controlled through the \Robject{ncenter} argument to \Rfunction{reduceDimension()} in Monocle 2. Using this algorithm drastically speeds up the computations in \Rfunction{reduceDimension()}, and also serves to regularize the manifold, often producing cleaner, more accurate single-cell trajectories. The DDRTree scheme works via the following optimization: \begin{equation}\label{eq:DDRTree_opt} \mathop{min}_{\mathbf{W}, \mathbf{Z}, \mathbf{B}, \mathbf{Y}, \mathbf{R}} \sum_{i = 1}^N ||\mathbf{x}_i - \mathbf{W} \mathbf{z}_i||^2 + \frac{\lambda}{2} \sum_{k, k'}b_{k, k'}||\mathbf{W} \mathbf{y}_k - \mathbf{W} \mathbf{y}_k'||^2 + \gamma\Big[\sum_{k = 1}^K \sum_{i = 1}^N r_{i, k} ||\mathbf{z}_i - \mathbf{y}_k||^2 + \sigma \Omega (\mathbf{R})\Big] \end{equation} The optimization is constrained such that $\mathbf{W}^T \mathbf{W} = \mathbf{I}, \sum_{k = 1}^K r_{i, k} = 1, r_{i, k} \leq 0, \forall i, \forall k$. The matrix $\mathbf{R} \in \mathcal{R}^{N \times N}$ is used to regularize the graph, through the \emph{negative entropy regularization} $\Omega(\mathbf{R}) = \sum_{i = 1}^N \sum_{k = 1}^k r_{i, k} \log \ r_{i, k}$. In effect, DDRTree uses the latent points $\{\mathbf{y}_k\}_{k = 1}^K$ as the centers of $K$ clusters. That is, the algorithm acts as soft K-means clustering on the points $\{\mathbf{z}_i\}^N_{i = 1}$, and jointly learns a graph on the $K$ cluster centers. The matrix $\mathbf{R}$ transforms the hard assignments used in K-means into soft assignments with $\sigma > 0$ as a regularization parameter. Problem \ref{eq:DDRTree_opt} again contains a number of analytical steps, and can be solved by alternating between the terms. Moreover, because some of the more expensive numerical operations involve matrices that are $K$ dimensional (instead of $N$ dimensional), they have complexity that is invariant of the size of the input data. Our accessory package \Rpackage{DDRTree} implements the DDRTree algorithm using a number of key performance optimizations. Monocle 2 calls DDRTree to learn the core manifold describing a \Robject{CellDataSet}, and then computes pseudotime coordinates and branch assignments using this manifold. Monocle also uses the manifold in downstream analysis steps such as BEAM. A manuscript describing the general stragegy of using reverse graph embedding for single-cell trajectory analysis will be forthcoming (Qiu et al, in preparation). % % \section{Using Monocle with other data types} % % This section is yet to be written. \section{Citation}\ref{citations} If you use Monocle to analyze your experiments, please cite: <>= citation("monocle") @ %def \section{Acknowledgements} Monocle was originally built by Cole Trapnell and Davide Cacchiarelli, with substantial design input John Rinn and Tarjei Mikkelsen. We are grateful to Sharif Bordbar, Chris Zhu, Amy Wagers and the Broad RNAi platform for technical assistance, and Magali Soumillon for helpful discussions. Cole Trapnell is a Damon Runyon Postdoctoral Fellow. Davide Cacchiarelli is a Human Frontier Science Program Fellow. Cacchiarelli and Mikkelsen were supported by the Harvard Stem Cell Institute. John Rinn is the Alvin and Esta Star Associate Professor. This work was supported by NIH grants 1DP2OD00667, P01GM099117, and P50HG006193-01. This work was also supported in part by the Single Cell Genomics initiative, a collaboration between the Broad Institute and Fluidigm Inc. Monocle 2 was developed by Cole Trapnell's lab. Significant portions were written by Xiaojie Qiu. The work was supported by NIH grant 1DP2HD088158 as well as an Alfred P. Sloan Foundation Research Fellowship. This vignette was created from Wolfgang Huber's Bioconductor vignette style document, and patterned after the vignette for \Rpackage{DESeq}, by Simon Anders and Wolfgang Huber. \section{Session Info} <>= sessionInfo() @ \bibliographystyle{unsrt} \bibliography{monocle_alg} \end{document}