Overview

Instructor(s) name(s) and contact information

Workshop Description

In single-cell RNA-sequencing (scRNA-seq), gene expression is assessed at the level of single cells. In dynamic biological systems, it may not be appropriate to assign cells to discrete groups, but rather a continuum of cell states may be observed, e.g. the differentiation of a stem cell population into mature cell types. This is often represented as a trajectory in a reduced dimension of the scRNA-seq dataset.

Many methods have been suggested for trajectory inference. However, in this setting, it is often unclear how one should handle multiple biological groups or conditions, e.g. constructing and comparing the differentiation trajectory of a wild type versus a knock-out stem cell population.

In this workshop, we will explore methods for comparing multiple conditions in a trajectory inference analysis. We start by integrating datasets from multiple conditions into a single trajectory. By comparing the conditions along the trajectory’s path, we can detect large-scale changes, indicative of differential progression. We also demonstrate how to detect subtler changes by finding genes that exhibit different behaviors between these conditions along a differentiation path.

This vignette provides a more complete problem description and proposes a few analytical approaches, which will serve as the basis of our workshop.

Pre-requisites

Software:

  • Basic knowledge of R syntax
  • Familiarity with single-cell RNA-sequencing
  • Familiarity with the SingleCellExperiment class

Background reading:

  • The textbook “Orchestrating Single-Cell Analysis with Bioconductor” is a great reference for single-cell analysis using Bioconductor packages.
  • The data used in this vignette was originally published in McFaline-Figueroa, et al. “A pooled single-cell genetic screen identifies regulatory checkpoints in the continuum of the epithelial-to-mesenchymal transition.”
  • Slingshot paper
  • tradeSeq paper

Workshop Participation

The workshop will start with an introduction to the problem and the dataset using presentation slides. Following this, we will have a lab session on how one may tackle the problem of handling multiple conditions in trajectory inference and in downstream analysis involving differential progression and differential expression.

R / Bioconductor packages used

Time outline

Activity Time
Introduction 15m
Data Integration and Trajectory Inference 10m
Differential Progression 15m
Differential Expression 15m
Wrap-up and Conclusions 5m

Workshop goals and objectives

Participants will learn how to reason about trajectories in single-cell RNA-seq data and how they may be used for interpretation of complex scRNA-seq datasets.

Learning goals

  • Reason about dynamic biological systems.
  • Grasp the complexity of analyzing large scRNA-seq datasets with the goal of extracting relevant biological information.
  • Understand the concepts of differential progression and differential expression along a trajectory path.

Learning objectives

  • Learn how to analyze single-cell RNA-seq data using Bioconductor packages.
  • Import and explore large scRNA-seq datasets.
  • Understand the challenges of trajectory inference.
  • Compose analysis pipeline that allows interpretation of complex scRNA-seq datasets.
  • Assess the added complexity of handling multiple conditions in these dynamic systems and how it influences the analysis pipeline.
  • Discuss how the analysis pipeline can incorporate this change and evaluate it.
suppressPackageStartupMessages({
  library(slingshot); library(SingleCellExperiment)
  library(RColorBrewer); library(scales)
  library(viridis); library(UpSetR)
  library(pheatmap); library(msigdbr)
  library(fgsea); library(knitr)
  library(ggplot2); library(gridExtra)
})
## Warning: package 'knitr' was built under R version 4.0.2
## Warning: package 'ggplot2' was built under R version 4.0.1

Analysis

Dataset

The dataset we will be working with concerns a single-cell RNA-sequencing dataset consisting of two different experiments, which correspond to two treatments. McFaline-Figueroa et al. studied the epithelial-to-mesenchymal transition (EMT), where cells spatially migrate from the epithelium to the mesenchyme during development. This process will be described by a trajectory, reflecting the gene expression changes occurring during this migration. The authors furthermore studied both a control (Mock) condition, and a condition in the absence of transforming growth factor \(\beta\) (TGFB).

In summary, we will be investigating a trajectory consisting of a single lineage that represents the EMT. This lineage is studied in two different conditions; a control condition and a TGFB-absent condition.

Integration

Our dataset contains cells collected from samples undergoing two different treatment conditions which were necessarily collected separately. Hence, we will start with an integration step to combine these two sets of cells, similar to batch correction. Our goal is to remove the technical effects of the different sample collections while preserving any true, biological differences between the two treatment groups.

Data integration and normalization are complex problems and there are a variety of methods addressing each. Interested participants can explore the corresponding chapter of the Bioconductor Ebook. However, since neither is the main focus of this workshop, we elected to use an existing pipeline for these tasks. The full Seurat data integration workflow with SCTransform normalization is described in this vignette.

Since this whole step is quite slow, it will not be run during the workshop but the code is provided below, along with a function to download and preprocess the public data from GEO.

sce <- bioc2020trajectories::importRawData()

########################
### Split by condition and convert to Seurat
########################
assays(sce)$logcounts <- log1p(assays(sce)$counts)
sceMock <- sce[ ,colData(sce)$pheno$treatment_id=='Mock']
sceTGFB <- sce[ ,colData(sce)$pheno$treatment_id=='TGFB']
library(Seurat)
soMock <- as.Seurat(sceMock)
soTGFB <- as.Seurat(sceTGFB)

########################
### Normalize
########################
soMock <- SCTransform(soMock, verbose = FALSE)
soTGFB <- SCTransform(soTGFB, verbose = FALSE)

########################
### Integrate
########################
dtlist <- list(Mock = soMock, TGFB = soTGFB)
intfts <- SelectIntegrationFeatures(object.list = dtlist, nfeatures = nrow(sce)) # maxes out at 4080 (why?)
dtlist <- PrepSCTIntegration(object.list = dtlist,
                             anchor.features = intfts)
anchors <- FindIntegrationAnchors(object.list = dtlist, normalization.method = "SCT",
                                  anchor.features = intfts)
integrated <- IntegrateData(anchorset = anchors, normalization.method = "SCT")
integrated <- RunPCA(integrated)
integrated <- RunUMAP(integrated, dims = 1:50)

## convert back to singleCellExperiment
sce <- as.SingleCellExperiment(integrated, assay="RNA")

Import dataset

We have made the pre-processed, integrated dataset available as a SingleCellExperiment object in the workshop package, which we import below.

data("sce", package = "bioc2020trajectories")

Differential Topology

Once the two datasets have been integrated, we can visualize all the single cells in a shared reduced dimensional space. We also visualize the distribution of cells in this space according to the treatment (control and TGFB) and spatial location (inner cells versus outer cells).

shuffle <- sample(ncol(sce))
layout(matrix(1:2, nrow = 1))
par(mar = c(4.5,4,1,1))

plot(reducedDims(sce)$UMAP[shuffle, ],
  asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2",
  col = alpha(c(1:2)[factor(colData(sce)$pheno$treatment_id)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 1:2, bty = "n",
       legend = levels(factor(colData(sce)$pheno$treatment_id)))
plot(reducedDims(sce)$UMAP[shuffle, ], asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2",
     col = alpha(c(3, 4)[factor(colData(sce)$pheno$spatial_id)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 3:4, bty = "n", legend = levels(factor(colData(sce)$pheno$spatial_id)))

layout(1)
par(mar = c(5, 4, 4, 2) + .1)

We know from biological knowledge that the EMT development goes from the inner to the outer cells. The question is: should we fit a separate trajectory for each condition? We might expect the trajectory itself to be changed by the treatment if the treatment effect is systematically large. Otherwise, the treatment may impact the expression profile of some genes but the overall trajectory will be preserved.

To help assess this, we devised an imbalance score. Regions with a high score indicate that the local cell distribution according to treatment label is unbalanced compared the overall distribution. Here, we see that, while there are some small regions of imbalance, the global path along the development axis is well-balanced. This means that we can fit a global trajectory to the full dataset. This choice allows us to use the maximal amount of data in the construction of our trajectory, which should lead to more robust results than separate, potentially noisy trajectories constructed on subsets of the data. As we will see, not all cell types are necessarily present in every condition, so this approach ensures that our trajectory accounts for all cell types present in the overall data.

scores <- bioc2020trajectories::imbalance_score(
  rd = reducedDims(sce)$UMAP,
  cl = colData(sce)$pheno$treatment_id,
  k = 20, smooth = 40)
## Warning: replacing previous import 'HDF5Array::path' by 'igraph::path' when
## loading 'bioc2020trajectories'
## Warning: replacing previous import 'igraph::compose' by 'purrr::compose' when
## loading 'bioc2020trajectories'
## Warning: replacing previous import 'igraph::simplify' by 'purrr::simplify' when
## loading 'bioc2020trajectories'
grad <- viridis::plasma(10, begin = 0, end = 1)
names(grad) <- levels(cut(scores$scaled_scores, breaks = 10))
plot(reducedDims(sce)$UMAP, col = grad[cut(scores$scaled_scores, breaks = 10)],
     asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2", cex = .8)
legend("topleft", legend = names(grad), col = grad, pch = 16, bty = "n", cex = 2 / 3)

For more information on the score, run help("imbalance_score", "bioc2020trajectories")

Trajectory Inference

We perform trajectory inference to order the cells according to EMT progression. We use slingshot for trajectory inference, with the cells’ position (inner or outer) serving as the cluster identifier. This ensures that we will only find a single lineage while still allowing sufficient flexibility to correctly orient the pseudotime axis.

Note that we perform trajectory inference using cells from both conditions, rather than splitting the data into two groups, as discussed above.

library(slingshot)
sce <- slingshot(sce, reducedDim = 'UMAP', clusterLabels = colData(sce)$pheno$spatial_id,
                 start.clus = 'inner', approx_points = 150)
## Using full covariance matrix
## Warning in if (class(X) == "dist") X <- as.matrix(X): the condition has length >
## 1 and only the first element will be used

Differential progression

Now that we have ordered the cells by EMT progression, we can begin to address the main question: how is this progression affected by TGF-\(\beta\) treatment? In this section, we interpret this question as a univariate analysis of the pseudotime values between the two groups.

The density estimates for the two groups show a trimodal distribution for the untreated cells, but a tendency toward later pseudotime values in the TGF-\(\beta\) treated cells. The difference is striking enough that a standard T-test would likely be significant, but we are we are interested more generally in differences between the two distributions, not just the difference of means (one could imagine a scenario in which the treated group tended toward the extremes, but the means were the same). Hence, we propose a Kolmogorov-Smirnov Test to assess whether the two groups of pseudotime values are derived from the same distribution. For more info on the Kolmogorov-Smirnov Test, see here.

########################
### Kolmogorov-Smirnov Test
########################
ks.test(slingPseudotime(sce)[colData(sce)$pheno$treatment_id == "Mock", 1],
        slingPseudotime(sce)[colData(sce)$pheno$treatment_id == "TGFB", 1])

As we might expect from the plot, this test is highly significant, so we can conclude that there are differences between the two distributions.

Differential expression

We will now proceed to discover genes whose expression is associated with the inferred trajectory. We will look for genes that (i) change in gene expression along the trajectory, and (ii) are differentially expressed between the two conditions along the trajectory. The differential expression analysis uses the Bioconductor package tradeSeq. This analysis relies on a new version of tradeSeq, which can be installed from the conditions branch on GitHub.

For each condition (i.e., control and TGF-Beta), a smooth average expression profile along pseudotime will be estimated for each gene, using a negative binomial generalized additive model (NB-GAM). Each differential expression hypothesis of interest will then be translated into testing specific features (a linear combination of the parameters) of this smoothed expression estimate.

The next two paragraphs can be time-consuming so we will not run them during the workshop, however, their output is already present in the data object that was loaded at the start of this workshop. They can be easily parallelized, relying on the BiocParallel bioconductor package. See here for more details.

Select number of knots

Before we can fit these smoothed expression profiles, we need to get a sense of how complex the expression patterns are in this dataset. This is translated into selecting a number of knots for the NB-GAMs, where a higher number of knots allows for more complex expression patterns. Here, we pick \(5\) knots.

See here for more information on smoothers.

library(tradeSeq)
set.seed(3)
icMat <- evaluateK(counts = as.matrix(assays(sce)$counts),
                   pseudotime = colData(sce)$slingshot$pseudotime,
                   cellWeights = colData(sce)$slingshot$cellWeights.V1,
                   conditions = factor(colData(sce)$pheno$treatment_id),
                   nGenes = 300,
                   k = 3:7)

The plot above shows the graphical output from running evaluateK. The left panel shows the distribution of gene-level AIC values as compared to their average AIC over the range of k. The second and third panel plot the average AIC and relative AIC with respect to the lowest value of k (i.e., 3), respectively, as a function of k. Finally, the right panel plots the number of genes whose AIC is lowest at a particular value of k.

Choosing an appropriate value of k can be seen as analogous to choosing the number of principal components based on a scree plot: we look for an ‘elbow point’, where the decrease starts attenuating. Here, we choose k=5 to allow for flexible, yet simple, functions while limiting the computational burden. In general, we found the influence of choosing the exact value of k to be rather limited, unless k is arbitrarily small or large. In our evaluations, most datasets fall within the range of \(4\) to \(8\) knots.

Fit GAM

Next, we fit the NB-GAMs using 5 knots, based on the pseudotime and cell-level weights estimated by Slingshot. We use the conditions argument to fit separate smoothers for each condition.

set.seed(3)
sce <- fitGAM(counts = as.matrix(assays(sce)$counts),
              pseudotime = colData(sce)$slingshot$pseudotime,
              cellWeights = colData(sce)$slingshot$cellWeights.V1,
              conditions = factor(colData(sce)$pheno$treatment_id),
              nknots = 5)
mean(rowData(sce)$tradeSeq$converged)

Assess DE along pseudotime (or pseudospace)

Note that the axis represented by the trajectory in this dataset is actually the migration of cells from the epithelium to the mesenchyme and therefore could also be looked at as a space dimension, although it is likely to be correlated with chronological time, too.

To assess significant changes in gene expression as a function of pseudotime within each lineage, we use the associationTest, which tests the null hypothesis that gene expression is not a function of pseudotime, i.e., whether the estimated smoothers are significantly varying as a function of pseudotime within each lineage. The lineages=TRUE argument specifies that we would like the results for each lineage separately, asides from the default global test, which tests for significant associations across all lineages in the trajectory simultaneously. Further, we specify a log2 fold change cut-off to test against using the l2fc argument.

On a lineage-specific basis, there are over twice as much DE genes in the mock lineage (2398) as compared to the TGFB lineage (1013). Many of the DE genes in the TGFB condition are also DE in the mock condition, around 80%.

The authors of the original paper found \(1105\) DE genes for the mock condition on a FDR level of \(1e-10\) and a cut-off of 1 on the absolute value of the log2 fold change.

library(tradeSeq)
assocRes <- associationTest(sce, lineages = TRUE, l2fc = log2(2))
mockGenes <-  rownames(assocRes)[
  which(p.adjust(assocRes$pvalue_lineage1_conditionMock, "fdr") <= 0.05)
]
tgfbGenes <-  rownames(assocRes)[
  which(p.adjust(assocRes$pvalue_lineage1_conditionTGFB, "fdr") <= 0.05)
]

length(mockGenes)
## [1] 2398
length(tgfbGenes)
## [1] 1013
UpSetR::upset(fromList(list(mock = mockGenes, tgfb = tgfbGenes)))

Visualization of DE genes

Below we visualize and cluster the genes whose expression vary over pseudotime, using the smoothed expression patterns As was also observed in the original manuscript, genes are mainly upregulated at the start- or endpoints of the lineage.

### based on mean smoother
yhatSmooth <- predictSmooth(sce, gene = mockGenes, nPoints = 50, tidy = FALSE)
heatSmooth <- pheatmap(t(scale(t(yhatSmooth[, 1:50]))),
                       cluster_cols = FALSE,
                       show_rownames = FALSE,
                       show_colnames = FALSE)

## the hierarchical trees constructed here, can also be used for 
## clustering of the genes according to their average expression pattern.
cl <- sort(cutree(heatSmooth$tree_row, k = 6))
table(cl)
## cl
##   1   2   3   4   5   6 
## 269 778 197 740 338  76
conditions <- colData(sce)$pheno$treatment_id
pt1 <- colData(sce)$slingshot$pseudotime

### based on fitted values (plotting takes a while to run)
yhatCell <- predictCells(sce, gene=mockGenes)
yhatCellMock <- yhatCell[,conditions == "Mock"]
# order according to pseudotime
ooMock <- order(pt1[conditions == "Mock"], decreasing=FALSE)
yhatCellMock <- yhatCellMock[,ooMock]
pheatmap(t(scale(t(yhatCellMock))), cluster_cols = FALSE,
          show_rownames = FALSE, show_colnames=FALSE)

Gene set enrichment analysis on genes from the Mock condition

Gene set enrichment analysis on the DE genes within the mock condition confirms the biology on epithelial cell differentiation.

## C5 category is according to gene ontology grouping: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4707969/pdf/nihms-743907.pdf
geneSets <- msigdbr(species = "Mus musculus", category = "C5", subcategory = "BP")
### filter background to only include genes that we assessed.
geneSets$gene_symbol <- toupper(geneSets$gene_symbol)
geneSets <- geneSets[geneSets$gene_symbol %in% names(sce),]
m_list <- geneSets %>% split(x = .$gene_symbol, f = .$gs_name)
stats <- assocRes$waldStat_lineage1_conditionMock
names(stats) <- rownames(assocRes)
eaRes <- fgsea(pathways = m_list, stats = stats, nperm = 5e4, minSize = 10)
## Warning in fgsea(pathways = m_list, stats = stats, nperm = 50000, minSize = 10):
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (0.01% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
ooEA <- order(eaRes$pval, decreasing = FALSE)
kable(head(eaRes[ooEA, 1:3], n = 20))
pathway pval padj
GO_MITOTIC_CELL_CYCLE 0.0000200 0.0787184
GO_EPIDERMIS_DEVELOPMENT 0.0000400 0.0787184
GO_EPITHELIAL_CELL_DIFFERENTIATION 0.0000800 0.0787184
GO_MITOTIC_NUCLEAR_DIVISION 0.0000800 0.0787184
GO_CELL_CYCLE 0.0001000 0.0787184
GO_SKIN_DEVELOPMENT 0.0001400 0.0899639
GO_EPITHELIUM_DEVELOPMENT 0.0001600 0.0899639
GO_PROTEOLYSIS 0.0002200 0.1049579
GO_ORGANELLE_FISSION 0.0002400 0.1049579
GO_NEGATIVE_REGULATION_OF_CELLULAR_COMPONENT_ORGANIZATION 0.0002800 0.1076339
GO_REGULATION_OF_MODIFICATION_OF_SYNAPTIC_STRUCTURE 0.0003008 0.1076339
GO_CELL_CYCLE_PROCESS 0.0003400 0.1115178
GO_EPIDERMAL_CELL_DIFFERENTIATION 0.0003800 0.1150500
GO_CYTOSKELETON_ORGANIZATION 0.0005200 0.1443171
GO_REGULATION_OF_ORGANELLE_ORGANIZATION 0.0005800 0.1443171
GO_KERATINOCYTE_DIFFERENTIATION 0.0006200 0.1443171
GO_NEGATIVE_REGULATION_OF_ORGANELLE_ORGANIZATION 0.0006400 0.1443171
GO_REGULATION_OF_CHROMOSOME_ORGANIZATION 0.0006600 0.1443171
GO_REGULATION_OF_NUCLEAR_DIVISION 0.0007600 0.1574369
GO_CHROMOSOME_ORGANIZATION 0.0009800 0.1911733

Differential expression between conditions

In the following sections, we will investigate differential expression for each gene, between the two conditions.

We will first make exploratory data analysis visualizations to take a look at the expression patterns of genes that were also discussed in the original manuscript. The paper mentions that CDH1 and CRB3 should be expressed in similar kinetics. Note that the lower slope of CDH1 is also observed in the paper.

plotSmoothers(sce, assays(sce)$counts, gene = "CDH1", alpha = 1, border = TRUE) + ggtitle("CDH1")

plotSmoothers(sce, assays(sce)$counts, gene = "CRB3", alpha = 1, border = TRUE) + ggtitle("CRB3")

They also mention that ‘only cells treated with TGF-Beta and positioned at the outer extreme of the trajectory expressed robust levels of FN1 and CDH2’.

plotSmoothers(sce, assays(sce)$counts, gene = "FN1", alpha = 1, border = TRUE) + ggtitle("FN1")

plotSmoothers(sce, assays(sce)$counts, gene = "CDH2", alpha = 1, border = TRUE) + ggtitle("CDH2")

Differential expression analysis

To test differential expression between conditions, we use the conditionTest function implemented in tradeSeq. This function tests the null hypothesis that genes have identical expression patterns in each condition. We discover \(1993\) genes that are DE with a fold change higher than \(2\) or lower than \(1/2\).

condRes <- conditionTest(sce, l2fc = log2(2))
condRes$padj <- p.adjust(condRes$pvalue, "fdr")
mean(condRes$padj <= 0.05, na.rm = TRUE)
## [1] 0.1890175
sum(condRes$padj <= 0.05, na.rm = TRUE)
## [1] 1993
conditionGenes <- rownames(condRes)[condRes$padj <= 0.05]
conditionGenes <- conditionGenes[!is.na(conditionGenes)]

Visualize most and least significant gene

# plot genes
oo <- order(condRes$waldStat, decreasing = TRUE)

# most significant gene
plotSmoothers(sce, assays(sce)$counts,
              gene = rownames(assays(sce)$counts)[oo[1]],
              alpha = 1, border = TRUE)

# least significant gene
plotSmoothers(sce, assays(sce)$counts,
              gene = rownames(assays(sce)$counts)[oo[nrow(sce)]],
              alpha = 1, border = TRUE)

Heatmaps of genes DE between conditions

Below we show heatmaps of the genes DE between conditions. The DE genes in the heatmaps are ordered according to a hierarchical clustering on the TGF-Beta condition.

### based on mean smoother
yhatSmooth <- predictSmooth(sce, gene = conditionGenes, nPoints = 50, tidy = FALSE)
heatSmooth_TGF <- pheatmap(t(scale(t(yhatSmooth[, 51:100]))),
  cluster_cols = FALSE,
  show_rownames = FALSE, show_colnames = FALSE, main = "TGF-Beta", legend = FALSE,
  silent = TRUE
)

matchingHeatmap_mock <- pheatmap(t(scale(t(yhatSmooth[heatSmooth_TGF$tree_row$order, 1:50]))),
  cluster_cols = FALSE, cluster_rows = FALSE,
  show_rownames = FALSE, show_colnames = FALSE, main = "Mock",
  legend = FALSE, silent = TRUE
)

grid.arrange(heatSmooth_TGF[[4]], matchingHeatmap_mock[[4]], ncol = 2)

Gene set enrichment analysis

Gene set enrichment analysis on genes that are differentially expressed between conditions finds evidence for cell motility, cell junctions/adhesion and gastrulation. The original paper also focuses on the KRAS signaling pathway, which induces cell migration, amongst others. Other related processes include morphogenesis, gastrulation and cell adhesion.

statsCond <- condRes$waldStat
names(statsCond) <- rownames(condRes)
eaRes <- fgsea(pathways = m_list, stats = statsCond, nperm = 5e4, minSize = 10)
## Warning in fgsea(pathways = m_list, stats = statsCond, nperm = 50000,
## minSize = 10): You are trying to run fgseaSimple. It is recommended to use
## fgseaMultilevel. To run fgseaMultilevel, you need to remove the nperm argument
## in the fgsea function call.
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (0.04% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaSimple(...): There were 1 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
ooEA <- order(eaRes$pval, decreasing = FALSE)
kable(head(eaRes[ooEA, 1:3], n = 20))
pathway pval padj
GO_BIOLOGICAL_ADHESION 0.0000200 0.0393892
GO_LOCOMOTION 0.0000200 0.0393892
GO_CELL_MOTILITY 0.0000400 0.0525189
GO_ANATOMICAL_STRUCTURE_FORMATION_INVOLVED_IN_MORPHOGENESIS 0.0000800 0.0787784
GO_REGULATION_OF_CELL_ADHESION 0.0002000 0.1312974
GO_TRANSCRIPTION_BY_RNA_POLYMERASE_III 0.0002000 0.1312974
GO_EXTRACELLULAR_STRUCTURE_ORGANIZATION 0.0003600 0.2025731
GO_REGULATION_OF_MITOCHONDRIAL_GENE_EXPRESSION 0.0004200 0.2067934
GO_MITOCHONDRIAL_RNA_METABOLIC_PROCESS 0.0011200 0.3892581
GO_REGULATION_OF_CELLULAR_COMPONENT_MOVEMENT 0.0011400 0.3892581
GO_CELL_JUNCTION_ORGANIZATION 0.0012400 0.3892581
GO_REGULATION_OF_CELL_POPULATION_PROLIFERATION 0.0013200 0.3892581
GO_CARDIOVASCULAR_SYSTEM_DEVELOPMENT 0.0013600 0.3892581
GO_CIRCULATORY_SYSTEM_DEVELOPMENT 0.0014200 0.3892581
GO_CELL_JUNCTION_ASSEMBLY 0.0015000 0.3892581
GO_GASTRULATION 0.0016600 0.3892581
GO_N_TERMINAL_PROTEIN_AMINO_ACID_MODIFICATION 0.0016800 0.3892581
GO_TUBE_DEVELOPMENT 0.0019600 0.4289048
GO_POSITIVE_REGULATION_OF_MULTICELLULAR_ORGANISMAL_PROCESS 0.0023600 0.4892555
GO_MOTOR_NEURON_AXON_GUIDANCE 0.0027399 0.5261273

Final notes

A compiled version of the vignette is available on the workshop website.

If you have questions that you could not ask during the workshop, feel free to open an issue on the github repository here.

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] tradeSeq_1.3.13             gridExtra_2.3              
##  [3] ggplot2_3.3.2               knitr_1.29                 
##  [5] fgsea_1.15.0                msigdbr_7.1.1              
##  [7] pheatmap_1.0.12             UpSetR_1.4.0               
##  [9] viridis_0.5.1               viridisLite_0.3.0          
## [11] scales_1.1.1                RColorBrewer_1.1-2         
## [13] SingleCellExperiment_1.11.6 SummarizedExperiment_1.19.5
## [15] DelayedArray_0.15.6         matrixStats_0.56.0         
## [17] Matrix_1.2-18               Biobase_2.49.0             
## [19] GenomicRanges_1.41.5        GenomeInfoDb_1.25.5        
## [21] IRanges_2.23.10             S4Vectors_0.27.12          
## [23] BiocGenerics_0.35.4         slingshot_1.7.0            
## [25] princurve_2.1.4            
## 
## loaded via a namespace (and not attached):
##  [1] Rtsne_0.15                    VGAM_1.1-3                   
##  [3] colorspace_1.4-1              ellipsis_0.3.1               
##  [5] rprojroot_1.3-2               XVector_0.29.3               
##  [7] fs_1.4.1                      farver_2.0.3                 
##  [9] ggrepel_0.8.2                 bit64_0.9-7                  
## [11] splines_4.0.0                 docopt_0.7.1                 
## [13] cluster_2.1.0                 dbplyr_1.4.4                 
## [15] bioc2020trajectories_0.0.0.91 HDF5Array_1.17.3             
## [17] compiler_4.0.0                httr_1.4.1                   
## [19] backports_1.1.8               assertthat_0.2.1             
## [21] limma_3.45.7                  htmltools_0.5.0              
## [23] tools_4.0.0                   igraph_1.2.5                 
## [25] gtable_0.3.0                  glue_1.4.1                   
## [27] GenomeInfoDbData_1.2.3        RANN_2.6.1                   
## [29] reshape2_1.4.4                dplyr_1.0.0                  
## [31] rappdirs_0.3.1                fastmatch_1.1-0              
## [33] Rcpp_1.0.4.6                  slam_0.1-47                  
## [35] DDRTree_0.1.5                 pkgdown_1.5.1                
## [37] vctrs_0.3.1                   rhdf5filters_1.1.0           
## [39] ape_5.4                       nlme_3.1-147                 
## [41] xfun_0.15                     stringr_1.4.0                
## [43] lifecycle_0.2.0               irlba_2.3.3                  
## [45] edgeR_3.31.4                  zlibbioc_1.35.0              
## [47] MASS_7.3-51.5                 monocle_2.17.0               
## [49] rhdf5_2.33.3                  yaml_2.2.1                   
## [51] curl_4.3                      pbapply_1.4-2                
## [53] memoise_1.1.0                 fastICA_1.2-2                
## [55] stringi_1.4.6                 RSQLite_2.2.0                
## [57] highr_0.8                     desc_1.2.0                   
## [59] densityClust_0.3              BiocParallel_1.23.0          
## [61] rlang_0.4.6                   pkgconfig_2.0.3              
## [63] bitops_1.0-6                  qlcMatrix_0.9.7              
## [65] evaluate_0.14                 lattice_0.20-41              
## [67] purrr_0.3.4                   Rhdf5lib_1.11.2              
## [69] labeling_0.3                  bit_1.1-15.2                 
## [71] tidyselect_1.1.0              plyr_1.8.6                   
## [73] magrittr_1.5                  R6_2.4.1                     
## [75] generics_0.0.2                combinat_0.0-8               
## [77] DBI_1.1.0                     pillar_1.4.4                 
## [79] withr_2.2.0                   mgcv_1.8-31                  
## [81] RCurl_1.98-1.2                tibble_3.0.1                 
## [83] crayon_1.3.4                  BiocFileCache_1.13.0         
## [85] rmarkdown_2.3                 locfit_1.5-9.4               
## [87] grid_4.0.0                    data.table_1.12.8            
## [89] blob_1.2.1                    FNN_1.1.3                    
## [91] HSMMSingleCell_1.9.0          sparsesvd_0.2                
## [93] digest_0.6.25                 munsell_0.5.0