--- title: 'Vignette for **tradeSeq**: **tra**jectory **d**ifferential **e**xpression analysis for **Seq**uncing data' author: "Koen Van den Berge and Hector Roux de Bézieux" bibliography: tradeSeq.bib date: "9/10/2018" output: rmarkdown::html_document: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{'Vignette for **tradeSeq**} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, echo = FALSE} library(knitr) # For a smaller vignette # This can be deleted when trying to follow along opts_chunk$set(dpi = 58) ``` # Introduction `tradeSeq` is an `R` package that allows analysis of gene expression along trajectories. While it has been developed and applied to single-cell RNA-sequencing (scRNA-seq) data, its applicability extends beyond that, and also allows the analysis of, e.g., single-cell ATAC-seq data along trajectories or bulk RNA-seq time series datasets. For every gene in the dataset, `tradeSeq` fits a generalized additive model (GAM) by building on the `mgcv` R package. It then allows statistical inference on the GAM by assessing contrasts of the parameters of the fitted GAM model, ading in interpreting complex datasets. All details about the `tradeSeq` model and statistical tests are described in our preprint [@VandenBerge2019a]. In this vignette, we analyse a subset of the data from [@Paul2015]. A `SingleCellExperiment` object of the data has been provided with the [`tradeSeq`](https://github.com/statOmics/tradeSeqpaper) package and can be retrieved as shown below. The data and UMAP reduced dimensions were derived from following the [Monocle 3 vignette](http://cole-trapnell-lab.github.io/monocle-release/monocle3/#tutorial-1-learning-trajectories-with-monocle-3). # Installation To install the package, simply run ```{r, eval = FALSE} if(!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("tradeSeq") ``` # Load data ```{r, warning=F, message=F} library(tradeSeq) library(RColorBrewer) library(SingleCellExperiment) library(slingshot) # For reproducibility RNGversion("3.5.0") palette(brewer.pal(8, "Dark2")) data(countMatrix, package = "tradeSeq") counts <- countMatrix rm(countMatrix) data(crv, package = "tradeSeq") data(celltype, package = "tradeSeq") ``` # Fit trajectories using slingshot We will fit developmental trajectories using the `slingshot` package [@Street2018a]. `slingshot` requires cluster labels as input, and fits trajectories in reduced dimension. We will use the reduced space calculated with the UMAP method, which is pre-calculated in the `se` object. We cluster the data using k-means with $7$ clusters. Since we know which cells are the progenitor cell type, we define the starting point for the trajectories as input for `slingshot`. Note that this argument is optional, and not required to run `slingshot`. ```{r, warning=FALSE, out.width="50%", fig.asp=.6} set.seed(200) rd <- reducedDims(crv) cl <- kmeans(rd, centers = 7)$cluster plot(rd, col = brewer.pal(9, "Set1")[cl], pch = 16, asp = 1, cex = 2/3) lin <- getLineages(rd, clusterLabels = cl, start.clus = 1) crv <- getCurves(lin) ``` We find two lineages for this dataset. The trajectory can be visualized using the `plotGeneCount` function, using either the cluster labels or cell type to color the cells. ```{r, out.width="50%", fig.asp=.6} plotGeneCount(curve = crv, clusters = cl) plotGeneCount(curve = crv, clusters = celltype, title = "Colored by cell type") ``` # Determining the number of knots After estimating the trajectory, we can fit generalized additive models (GAMs) with the `tradeSeq` package. Internally, this package builds on the `mgcv` package by fitting additive models using the `gam` function. The core function from `tradeSeq`, `fitGAM`, will use cubic splines as basis functions, and it tries to ensure that every lineage will end at a knot point of a smoother. By default, we allow for $6$ knots for every lineage, but this can be changed with the `nknots` argument. More knots will allow more flexibility, but also increase the risk of overfitting. Ideally, the number of knots should be selected to reach an optimal bias-variance trade-off for the smoother, where one explains as much variability in the expression data as possible with only a few regression coefficients. In order to guide that choice, we developed diagnostic plots using the Akaike Informaction Criterion (AIC). This is implemented in the `evaluateK` function in `tradeSeq`. The function takes as input the expression count matrix, a `SlingshotDataSet` (or alternatively, a matrix of pseudotimes and cell-level weights). The range of knots to evaluate is provided with the `knots` argument. The minimum allowed number of knots is $3$. While there is no boundary on the maximum number of knots, typically the interesting range is around $3$ to $10$ knots. The function will fit NB-GAM models for some number of genes, provided by the `nGenes` argument, over the range of knots that are provided, and return the AIC for each gene fitted with each number of knots. It is generally a good idea to evaluate this multiple times using different seeds (using the `set.seed` function), to check whether the results are reproducible across different gene subsets. This task can be computationally demanding, since the models must be fit multiple times for each gene. We therefore skip this in the vignette, but show the output one can expect instead. ```{r, eval=FALSE} icMat <- evaluateK(counts = counts, sds = crv, k=3:20, nGenes = 200, verbose=FALSE) # alternatively: icMat <- evaluateK(counts = counts, k=3:20, nGenes = 200, pseudotime = slingPseudotime(crv, na = FALSE), cellWeights = slingCurveWeights(crv)) ``` ```{r, echo=FALSE} knitr::include_graphics("evalK_paul1Cropped.png") knitr::include_graphics("evalK_paul2Cropped.png") ``` The output graphics are organized into four panels. The left panel plots a boxplot for each number of knots we wanted to evaluate. The plotted values are the deviation from a gene's AIC at that specific knot value from the average AIC of that gene across all the knots under evaluation. Typically, AIC values are somewhat higher for low number of knots, and we expect them to decrease as the number of knots gets higher. The two middle panels plot the average drop in AIC across all genes. The middle left panel simply plots the average AIC, while the middle right panel plots the change in AIC relative to the average AIC at the lowest knot number (here, this is 3 knots, as can also be seen from the plot since the relative AIC equals $1$). Finally, rhe right panel only plots a subset of genes where the AIC value changes significantly across the evaluated number of knots. Here, a significant change is defined as a change in absolute value of at least $2$, but this can be tuned using the `aicDiff` argument to `evaluateK`. For the subset of genes, a barplot is displayed that shows the number of genes that have their lowest AIC at a specific knot value. The middle panels show that the drop in AIC levels off if the number of knots is increased beyond $6$. In the right panel, $6$ knots also corresponds the highest number of genes with lowest AIC value. Based on these plots, we thus believe that fitting the NB-GAM models with $6$ knots is an appropriate choice. # Fit additive models By default, the GAM model estimates one smoother for every lineage using the negative binomial distribution. If you want to allow for other fixed effects (e.g., batch effects), then an additional model matrix, typically created using the `model.matrix` function, can be provided with the `U` argument. The precise model definition of the statistical model is described in our preprint [@VandenBerge2019a]. We use the effective library size, estimated with TMM [@Robinson2010], as offset in the model. We allow for alternatives by allowing a user-defined offset with the `offset` argument. This dataset consists of UMI counts, and we do not expect zero inflation to be a big problem. However, we also allow to fit zero inflated negative binomial (ZINB) GAMs by providing observation-level weights to `fitGAM` using the `weights` argument. The `weights` must correspond to the posterior probability that a count belongs to the count component of the ZINB distribution [@VandenBerge2018]. For the vignette, we fit smoothers for a filtered set of genes in the dataset, 239 genes in total. We also include the *Irf8* gene, since it is a known transcription factor involved in hematopoiesis. The `fitGAM` function relies on __BiocParallel__ to implement parallelization, progress bars and so on. Similar to `evaluateK`, fitGAM can either take a `SlingshotDataSet` object as input (`sds` argument), or a matrix of pseudotimes and cell-level weights (`pseudotime` and `cellWeights` argument). If a `SlingshotDataSet` is provided, the function will return a `SingleCellExperiment` object that contains the essential output from `tradeSeq`. This is much more efficient than providing the pseudotime and cell-level weights as matrices, when a list of GAM models will be returned. While in this vignette we will proceed with using the `sds` argument, hence a `SingleCellExperiment` object as output, `tradeSeq` allows input from any trajectory inference method with the `pseudotime` and `cellWeights` arguments. All functions work with both the `SingleCellExperiment` (i.e., `sds` input to `fitGAM`) output as well as the list output (i.e., `pseudotime` and `cellWeights` input to `fitGAM`). Because cells are assigned to a lineage based on their weigths, the result of `fitGAM` is stochastic. While this should have little impact on the overall results in practice, users are therefore encouraged to use the `set.seed` function before running `fitGAM` to ensure reproducibility of their analyses. ```{r} library(BiocParallel) library(magrittr) # Register BiocParallel Serial Execution (no parallelization in that case) BiocParallel::register(BiocParallel::SerialParam()) counts <- as.matrix(counts) sce <- fitGAM(counts = counts, sds = crv) # This takes about 1mn to run ``` You can also plot the cells in reduced dimension to see where the knots are located. ```{r, out.width="50%", fig.asp=.6} plotGeneCount(curve = crv, counts = counts, clusters = cl, models = sce) ``` # Within-lineage comparisons ## Association of gene expression with pseudotime A first exploration of the data analysis may consist in checking whether gene expression is associated with a particular lineage. The statistical test performed here, implemented in the `associationTest` function, is testing the null hypothesis that all smoother coefficients are equal to each other. This can be interpreted as testing whether the smoothed gene expression is significantly changing along pseudotime. ```{r} assoRes <- associationTest(sce) head(assoRes) ``` ## Discovering progenitor marker genes In order to discover marker genes of the progenitor cell population, researchers may be interested in assessing differential expression between the progenitor cell population (i.e., the starting point of a lineage) with the differentiated cell type population (i.e., the end point of a lineage). In the function `startVsEndTest`, we have implemented a Wald test that tests the null hypothesis that the expression at the starting point of the smoother (progenitor population) is identical to the expression at the end point of the smoother (differentiated population). The test basically involves a comparison between two smoother coefficients for every lineage. The function `startVsEndTest` performs an omnibus test across all lineages by default, but you can also assess all lineages separately by setting `lineages=TRUE`. Below, we adopt an omnibus test across the two lineages. ```{r} startRes <- startVsEndTest(sce) ``` We can visualize the estimated smoothers for the most significant gene. ```{r, out.width="40%", fig.asp=1} oStart <- order(startRes$waldStat, decreasing = TRUE) sigGeneStart <- names(sce)[oStart[1]] plotSmoothers(sce, counts, gene = sigGeneStart) ``` Alternatively, we can color the cells in UMAP space with that gene's expression. ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = sigGeneStart) ``` ## Comparing specific pseudotime values within a lineage The `startVsEndTest` compares two points on a lineage, and by default it is comparing the inception point with the end point. However, this is a specific form of a more general capability of the `startVsEndTest` to compare any two points on any lineage. If the interest lies in comparing any two custom pseudotime values, one can specify this using the `pseudotimeValues` arguments in `startVsEndTest`. For example, below we'd like to compare the expression for each gene at pseudotime values of $0.8$ and $0.1$. ```{r} customRes <- startVsEndTest(sce, pseudotimeValues = c(0.1, 0.8)) ``` # Between-lineage comparisons ## Discovering differentiated cell type markers `tradeSeq` can discover marker genes for the differentiated cell types by comparing the end points of the lineage-specific smoothers. This is implemented in the `diffEndTest` function. By default, `diffEndTest` performs an omnibus test, testing the null hypothesis that the endpoint expression is equal for all lineages using a multivariate Wald test. If more than two trajectories are present, one can assess all pairwise comparisons using the `pairwise=TRUE` argument. ```{r} endRes <- diffEndTest(sce) ``` We can plot the most significant gene using the `plotSmoothers` function. ```{r, out.width="40%", fig.asp=1} o <- order(endRes$waldStat, decreasing = TRUE) sigGene <- names(sce)[o[1]] plotSmoothers(sce, counts, sigGene) ``` Alternatively, we can color the cells in UMAP space with that gene's expression. ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = sigGene) ``` ## Discovering genes with different expression patterns Asides from testing at the level of the differentiated cell type, researchers may be interested in assessing the expression pattern of a gene over pseudotime. The function `patternTest` implements a statistical method that checks whether the smoothed gene expression is equal along pseudotime between two or multiple lineages. In practice, we use $100$ points, equally distributed along pseudotime, that are compared between two (or multiple) lineages, and this number can be changed using the `nPoints` argument. ```{r, out.width="40%", fig.asp=1} patternRes <- patternTest(sce) oPat <- order(patternRes$waldStat, decreasing = TRUE) head(rownames(patternRes)[oPat]) plotSmoothers(sce, counts, gene = rownames(patternRes)[oPat][1]) ``` ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = rownames(patternRes)[oPat][1]) ``` We find genes at the top that are also ranked as DE for the differentiated cell type. What is especially interesting are genes that have different expression patterns but no different expression at the differentiated cell type level. We therefore sort the genes according to the sum of square of their rank in increasing Wald statistics for the *patternTest* and their rank in decreasing Wald statistics for the *diffEndTest*. ```{r, out.width="50%", fig.asp=.8} library(dplyr) library(ggplot2) library(tidyr) compare <- inner_join(patternRes %>% mutate(Gene = rownames(patternRes), pattern = waldStat) %>% select(Gene, pattern), endRes %>% mutate(Gene = rownames(endRes), end = waldStat) %>% select(Gene, end), by = c("Gene" = "Gene")) %>% mutate(transientScore = (min_rank(desc(end)))^2 + (dense_rank(pattern))^2) ggplot(compare, aes(x = log(pattern), y = log(end))) + geom_point(aes(col = transientScore)) + labs(x = "patternTest Wald Statistic (log scale)", y = "diffEndTest Wald Statistic (log scale)") + scale_color_continuous(low = "yellow", high = "red") + theme_classic() ``` Or, we can visualize the expression in UMAP space of the top gene. ```{r, out.width="40%", fig.asp=1} topTransient <- (compare %>% arrange(desc(transientScore)))[1, "Gene"] plotSmoothers(sce, counts, gene = topTransient) ``` ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = topTransient) ``` Interestingly, we recover the Irf8 gene in the top 5 genes according to that ranking. ```{r} head(compare %>% arrange(desc(transientScore)) %>% select(Gene), n = 5) ``` We can also plot the Irf8 gene. ```{r, out.width="40%", fig.asp=1} plotSmoothers(sce, counts, gene = "Irf8") ``` ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = "Irf8") ``` ## Early drivers of differentiation Another question of interest is to find a list of genes that are differentially expressed around the separation of two or multiple lineages. The function `earlyDETest` implements a statistical method to tests the null hypothesis of whether the smoothers are equal between two user-specified knots by building on the `patternTest`, but restricting itself to a particular location of the smoothers. Again, the knots can be visualized with the `plotGeneCount` function. By selecting the region covering the first two knot points to test for differential patterns between the lineages, we check which genes are behaving differently around the bifurcation point. ```{r, out.width="50%", fig.asp=.5} plotGeneCount(curve = crv, counts = counts, clusters = cl, models = sce) earlyDERes <- earlyDETest(sce, knots = c(1, 2)) oEarly <- order(earlyDERes$waldStat, decreasing = TRUE) head(rownames(earlyDERes)[oEarly]) ``` ```{r, out.width="40%", fig.asp=1} plotSmoothers(sce, counts, gene = rownames(earlyDERes)[oEarly][2]) ``` ```{r, out.width="50%", fig.asp=.5} plotGeneCount(crv, counts, gene = rownames(earlyDERes)[oEarly][2]) ``` # Clustering of genes according to their expression pattern tradeSeq provides the functionality to cluster genes according to their expression pattern along the lineages with the `clusterExpressionPatterns` function. A number of equally spaced points for every lineage are selected to perform the clustering, and the number of points can be selected with the `nPoints` argument. The `genes` argument specifies which genes you want to cluster (e.g., all genes with differential expression patterns). Here, we use 20 points along each lineage to cluster the first 40 genes in the dataset. The clustering itself occurs by the `clusterExperiment` package [@Risso2018], hence the user may select any clustering algorithm that's built into that package, or custom clustering algorithms implemented by the user. For a list of built-in clustering algorithms within `clusterExperiment`, run `clusterExperiment::listBuiltInFunctions()` on the command line. ```{r, warning=FALSE,message=F} library(clusterExperiment) nPointsClus <- 20 clusPat <- clusterExpressionPatterns(sce, nPoints = nPointsClus, genes = rownames(counts)[1:200]) clusterLabels <- primaryCluster(clusPat$rsec) ``` The first 4 clusters can be visualized using the normalized expression upon which the clustering is based. Please note that the code below would only works for a trajectory with two lineages. Modify the code appropriately if using with a dataset with 3 lineages or more. ```{r, out.width="50%", fig.asp=1} library(cowplot) cUniq <- unique(clusterLabels) cUniq <- cUniq[!cUniq == -1] # remove unclustered genes plots <- list() for (xx in cUniq[1:4]) { cId <- which(clusterLabels == xx) p <- ggplot(data = data.frame(x = 1:nPointsClus, y = rep(range(clusPat$yhatScaled[cId, ]), nPointsClus / 2)), aes(x = x, y = y)) + geom_point(alpha = 0) + labs(title = paste0("Cluster ", xx), x = "Pseudotime", y = "Normalized expression") + theme_classic() for (ii in 1:length(cId)) { geneId <- rownames(clusPat$yhatScaled)[cId[ii]] p <- p + geom_line(data = data.frame(x = rep(1:nPointsClus, 2), y = clusPat$yhatScaled[geneId, ], lineage = rep(0:1, each = nPointsClus)), aes(col = as.character(lineage), group = lineage), lwd = 1.5) } p <- p + guides(color = FALSE) + scale_color_manual(values = c("orange", "darkseagreen3"), breaks = c("0", "1")) plots[[as.character(xx)]] <- p } plots$ncol <- 2 do.call(plot_grid, plots) ``` # tradeSeq list output If another method than Slingshot is used for trajectory inference, one can input custom pseudotimes and cell-level weights in `fitGAM`, as we also discussed above. The output from `fitGAM` will be different in that case, and less memory efficient. All functions we have discussed above work exactly the same with the list output. However, the list output functionality is a little bit bigger, and here we discuss some capabilities that are only available with the list output. ```{r} gamList <- fitGAM(counts, pseudotime = slingPseudotime(crv, na = FALSE), cellWeights = slingCurveWeights(crv), nknots = 6) ``` First, one may explore the results of a model by requesting its summary. ```{r} summary(gamList[["Irf8"]]) ``` Related to the `associationTest`, one can extract the p-values generated by the `mgcv` package using the `getSmootherPvalues` function. These p-values are derived from a test that assesses the null hypothesis that all smoother coefficients are equal to zero. Note, however, that their interpretation is thus more complex. A significant lineage for a particular gene might thus be the result of (a) a different mean expression in that lineage as compared to the overall expression of that gene, or (b) significantly varying expression along that lineage, even if the means are equal, or (c) a combination of both. This function extracts the p-values calculated by `mgcv` from the GAM, and will return `NA` for genes that we were unable to fit properly. Similarly, the test statistics may be extracted with `getSmootherTestStats`. Since this dataset was pre-filtered to only contain relevant genes, all p-values (test statistics) will be very low (high). Note, that these functions are only applicable with the list output of `tradeSeq`, and not with the `SingleCellExperiment` output. We will therefore not evaluate these here. ```{r, eval=FALSE} pvalLineage <- getSmootherPvalues(gamList) statLineage <- getSmootherTestStats(gamList) ``` # Convergence issues on small or zero-inflated datasets If you're working with a dataset that has a limited number of cells, or if you are incorporating zero inflation weights, the GAMs may be harder to fit, as noted by the warnings when running `fitGAM`. In that case, the situation might improve if you allow for more iterations in the GAM fitting. This can be done with the `control` argument of `fitGAM`. ```{r} library(mgcv) control <- gam.control() control$maxit <- 1000 #set maximum number of iterations to 1K # pass to control argument of fitGAM as below: # # gamList <- fitGAM(counts = counts, # pseudotime = slingPseudotime(crv, na = FALSE), # cellWeights = slingCurveWeights(crv), # control = control) ``` # Cheatsheet To recapitulate the workflow, we have created a cheatsheet that users can refer to when deciding which tests to run. ```{r, echo = F} # ggdraw() + draw_image("cheatsheet_highRes.jpeg") knitr::include_graphics("cheatsheet_highRes.jpeg") ``` # Session ```{r} sessionInfo() ``` # References