Trajectory inference across conditions: differential expression and differential progression
workshopTrajectories.RmdOverview
Instructor(s) name(s) and contact information
- Kelly Street (kstreet@ds.dfci.harvard.edu)
- Koen Van den Berge (koenvdberge@berkeley.edu)
- Hector Roux de Bézieux (hector.rouxdebezieux@berkeley.edu)
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
SingleCellExperimentclass
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
- The workshop will focus on Bioconductor packages SingleCellExperiment, Slingshot, and tradeSeq
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. Participants can follow along in the following fashions:
- Run the code online by going to orchestra and launching the workshop titled Trajectory inference across conditions: differential expression and differential progression
- Clone the repository from github and run the code locally
- USe docker and get the container kstreet13/bioc2020trajectories:latest
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.
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 under activation 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-induced 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)))
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
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 we have recently updated to allow for multiple conditions.
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. See here for more information on smoothers and the mgcv package which we are using to estimate the smoothers.
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.
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.
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.
rowData(sce)$assocRes <- associationTest(sce, lineages = TRUE, l2fc = log2(2))
library(tradeSeq)
assocRes <- rowData(sce)$assocRes
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] 2399
length(tgfbGenes)## [1] 1013
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 741 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.| pathway | pval | padj |
|---|---|---|
| GO_EPIDERMIS_DEVELOPMENT | 0.0000400 | 0.0792184 |
| GO_MITOTIC_CELL_CYCLE | 0.0000400 | 0.0792184 |
| GO_CELL_CYCLE | 0.0000600 | 0.0792184 |
| GO_EPITHELIUM_DEVELOPMENT | 0.0001000 | 0.0990230 |
| GO_SKIN_DEVELOPMENT | 0.0001600 | 0.1009971 |
| GO_MITOTIC_NUCLEAR_DIVISION | 0.0002000 | 0.1009971 |
| GO_EPITHELIAL_CELL_DIFFERENTIATION | 0.0002200 | 0.1009971 |
| GO_PROTEOLYSIS | 0.0002400 | 0.1009971 |
| GO_REGULATION_OF_CHROMOSOME_ORGANIZATION | 0.0002400 | 0.1009971 |
| GO_NEGATIVE_REGULATION_OF_CELLULAR_COMPONENT_ORGANIZATION | 0.0002800 | 0.1009971 |
| GO_REGULATION_OF_MODIFICATION_OF_SYNAPTIC_STRUCTURE | 0.0002805 | 0.1009971 |
| GO_ORGANELLE_FISSION | 0.0003200 | 0.1035933 |
| GO_CHROMOSOME_ORGANIZATION | 0.0003400 | 0.1035933 |
| GO_CELL_CYCLE_PROCESS | 0.0003800 | 0.1075107 |
| GO_REGULATION_OF_ORGANELLE_ORGANIZATION | 0.0004200 | 0.1109058 |
| GO_EPIDERMAL_CELL_DIFFERENTIATION | 0.0005000 | 0.1164977 |
| GO_NEGATIVE_REGULATION_OF_ORGANELLE_ORGANIZATION | 0.0005000 | 0.1164977 |
| GO_REGULATION_OF_NUCLEAR_DIVISION | 0.0006400 | 0.1408327 |
| GO_KERATINOCYTE_DIFFERENTIATION | 0.0007800 | 0.1626062 |
| GO_CELLULAR_PROTEIN_CATABOLIC_PROCESS | 0.0008800 | 0.1656385 |
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.1889995
sum(condRes$padj <= 0.05, na.rm = TRUE)## [1] 1993Visualize 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)
yhatSmoothScaled <- t(scale(t(yhatSmooth)))
heatSmooth_TGF <- pheatmap(yhatSmoothScaled[, 51:100],
cluster_cols = FALSE,
show_rownames = FALSE, show_colnames = FALSE, main = "TGF-Beta", legend = FALSE,
silent = TRUE
)
matchingHeatmap_mock <- pheatmap(yhatSmoothScaled[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.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.| pathway | pval | padj |
|---|---|---|
| GO_BIOLOGICAL_ADHESION | 0.0000200 | 0.0792784 |
| GO_ANATOMICAL_STRUCTURE_FORMATION_INVOLVED_IN_MORPHOGENESIS | 0.0000400 | 0.0792784 |
| GO_CELL_MOTILITY | 0.0000600 | 0.0792784 |
| GO_REGULATION_OF_CELL_ADHESION | 0.0001400 | 0.1387372 |
| GO_REGULATION_OF_MITOCHONDRIAL_GENE_EXPRESSION | 0.0003000 | 0.2378352 |
| GO_TRANSCRIPTION_BY_RNA_POLYMERASE_III | 0.0005400 | 0.3171137 |
| GO_EXTRACELLULAR_STRUCTURE_ORGANIZATION | 0.0006400 | 0.3171137 |
| GO_VASCULATURE_DEVELOPMENT | 0.0008800 | 0.3171137 |
| GO_CIRCULATORY_SYSTEM_DEVELOPMENT | 0.0009200 | 0.3171137 |
| GO_CELL_MORPHOGENESIS_INVOLVED_IN_DIFFERENTIATION | 0.0010800 | 0.3171137 |
| GO_REGULATION_OF_CELLULAR_COMPONENT_MOVEMENT | 0.0011000 | 0.3171137 |
| GO_MITOCHONDRIAL_RNA_METABOLIC_PROCESS | 0.0011800 | 0.3171137 |
| GO_REGULATION_OF_CELL_POPULATION_PROLIFERATION | 0.0012000 | 0.3171137 |
| GO_CELL_MORPHOGENESIS | 0.0013000 | 0.3171137 |
| GO_TAXIS | 0.0013400 | 0.3171137 |
| GO_CELL_CELL_ADHESION | 0.0014200 | 0.3171137 |
| GO_MOTOR_NEURON_AXON_GUIDANCE | 0.0014601 | 0.3171137 |
| GO_N_TERMINAL_PROTEIN_AMINO_ACID_MODIFICATION | 0.0015200 | 0.3171137 |
| GO_CELL_PROJECTION_ORGANIZATION | 0.0015800 | 0.3171137 |
| GO_TUBE_DEVELOPMENT | 0.0016000 | 0.3171137 |
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 Under development (unstable) (2020-12-31 r79746)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.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.25 gridExtra_2.3
## [3] ggplot2_3.3.3 knitr_1.30
## [5] fgsea_1.17.0 msigdbr_7.2.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.13.3 SummarizedExperiment_1.21.1
## [15] Biobase_2.51.0 GenomicRanges_1.43.1
## [17] GenomeInfoDb_1.27.3 IRanges_2.25.6
## [19] S4Vectors_0.29.6 BiocGenerics_0.37.0
## [21] MatrixGenerics_1.3.0 matrixStats_0.57.0
## [23] slingshot_1.9.1 princurve_2.1.5
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 VGAM_1.1-4
## [3] colorspace_2.0-0 ellipsis_0.3.1
## [5] rprojroot_2.0.2 XVector_0.31.1
## [7] fs_1.5.0 farver_2.0.3
## [9] ggrepel_0.9.0 bit64_4.0.5
## [11] splines_4.1.0 docopt_0.7.1
## [13] cluster_2.1.0 dbplyr_2.0.0
## [15] bioc2020trajectories_0.0.0.93 HDF5Array_1.19.0
## [17] compiler_4.1.0 httr_1.4.2
## [19] assertthat_0.2.1 Matrix_1.3-2
## [21] limma_3.47.3 htmltools_0.5.0
## [23] tools_4.1.0 igraph_1.2.6
## [25] gtable_0.3.0 glue_1.4.2
## [27] GenomeInfoDbData_1.2.4 RANN_2.6.1
## [29] reshape2_1.4.4 dplyr_1.0.2
## [31] rappdirs_0.3.1 fastmatch_1.1-0
## [33] Rcpp_1.0.5 slam_0.1-48
## [35] DDRTree_0.1.5 pkgdown_1.6.1
## [37] rhdf5filters_1.3.3 vctrs_0.3.6
## [39] ape_5.4-1 nlme_3.1-151
## [41] xfun_0.20 stringr_1.4.0
## [43] lifecycle_0.2.0 irlba_2.3.3
## [45] edgeR_3.33.0 zlibbioc_1.37.0
## [47] ragg_0.4.0 rhdf5_2.35.0
## [49] monocle_2.19.0 curl_4.3
## [51] yaml_2.2.1 memoise_1.1.0
## [53] pbapply_1.4-3 fastICA_1.2-2
## [55] stringi_1.5.3 RSQLite_2.2.1
## [57] highr_0.8 desc_1.2.0
## [59] filelock_1.0.2 densityClust_0.3
## [61] BiocParallel_1.25.2 rlang_0.4.10
## [63] pkgconfig_2.0.3 systemfonts_0.3.2
## [65] bitops_1.0-6 qlcMatrix_0.9.7
## [67] evaluate_0.14 lattice_0.20-41
## [69] Rhdf5lib_1.13.0 purrr_0.3.4
## [71] labeling_0.4.2 bit_4.0.4
## [73] tidyselect_1.1.0 plyr_1.8.6
## [75] magrittr_2.0.1 R6_2.5.0
## [77] generics_0.1.0 combinat_0.0-8
## [79] DelayedArray_0.17.7 DBI_1.1.0
## [81] pillar_1.4.7 withr_2.3.0
## [83] mgcv_1.8-33 RCurl_1.98-1.2
## [85] tibble_3.0.4 crayon_1.3.4
## [87] BiocFileCache_1.15.1 rmarkdown_2.6
## [89] locfit_1.5-9.4 grid_4.1.0
## [91] data.table_1.13.6 blob_1.2.1
## [93] FNN_1.1.3 HSMMSingleCell_1.11.0
## [95] sparsesvd_0.2 digest_0.6.27
## [97] textshaping_0.2.1 munsell_0.5.0