| | 1 | |
| | 2 | === Making diffusion maps with destiny === |
| | 3 | |
| | 4 | Seurat v3 doesn't have a function to make diffusion maps anymore. The alternatives are to 1) run the analysis in python using scanpy or 2) make a single cell object, i.e. from the Seurat object, and run the "DiffusionMap" function from the destiny package. Example commands for the second option are shown below. |
| | 5 | |
| | 6 | * **Making a single cell object from a Seurat object** |
| | 7 | {{{ |
| | 8 | library(Seurat) |
| | 9 | library(SingleCellExperiment) |
| | 10 | sce <- as.SingleCellExperiment(seuratObject) |
| | 11 | #this has the cell classification |
| | 12 | table(sce$ident) |
| | 13 | }}} |
| | 14 | |
| | 15 | * **Running diffusion map** |
| | 16 | <code> |
| | 17 | #this step may take a long time (days) or not finish. It is recommend to send it to the cluster as a script that reads the Seurat or the single cell object, runs DiffusionMap, and saves the object. |
| | 18 | library(destiny) |
| | 19 | dm <- DiffusionMap(sce, verbose = TRUE) |
| | 20 | }}} |
| | 21 | |
| | 22 | * **Plotting the diffusion map** |
| | 23 | {{{ |
| | 24 | library(ggplot2) |
| | 25 | cellLabels <- sce$ident |
| | 26 | tmp <- data.frame(DC1 = eigenvectors(dm)[, 1], |
| | 27 | DC2 = eigenvectors(dm)[, 2], |
| | 28 | DC3 = eigenvectors(dm)[, 3], |
| | 29 | DC4 = eigenvectors(dm)[, 4], |
| | 30 | Samples = cellLabels) |
| | 31 | pdf("./DC1_DC2.pdf", w=11, h=8.5) |
| | 32 | ggplot(tmp, aes(x = DC1, y = DC2, colour = Samples)) + |
| | 33 | geom_point() + |
| | 34 | xlab("Diffusion component 1") + |
| | 35 | ylab("Diffusion component 2") + |
| | 36 | theme_classic() |
| | 37 | dev.off() |
| | 38 | }}} |
| | 39 | * **Plotting cell progression along the diffusion map components** |
| | 40 | {{{ |
| | 41 | sce$pseud_dm1 <- rank(eigenvectors(dm)[,1]) # rank cells by their dpt dm1 |
| | 42 | sce$pseud_dm2 <- rank(eigenvectors(dm)[,2]) # rank cells by their dpt dm2 |
| | 43 | sce$pseud_dm1R <- rank(-eigenvectors(dm)[,1]) # rank cells by their dpt dm1 reverse order |
| | 44 | sce$pseud_dm2R <- rank(-eigenvectors(dm)[,2]) # rank cells by their dpt dm2 reverse order |
| | 45 | |
| | 46 | SortedDM1 <- data.frame(DM1Sort = as.data.frame(colData(sce))$pseud_dm1, |
| | 47 | Samples = as.data.frame(colData(sce))$ident) |
| | 48 | SortedDM2 <- data.frame(DM2Sort = as.data.frame(colData(sce))$pseud_dm2, |
| | 49 | Samples = as.data.frame(colData(sce))$ident) |
| | 50 | SortedDM1R <- data.frame(DM1SortR = as.data.frame(colData(sce))$pseud_dm1R, |
| | 51 | Samples = as.data.frame(colData(sce))$ident) |
| | 52 | SortedDM2R <- data.frame(DM2SortR = as.data.frame(colData(sce))$pseud_dm2R, |
| | 53 | Samples = as.data.frame(colData(sce))$ident) |
| | 54 | |
| | 55 | ggplot(SortedDM1, aes(x=SortedDM1[,1], y=Samples,color=Samples)) + |
| | 56 | geom_jitter() + xlab("Diffusion component 1 (DC1)") + ylab("Samples") + |
| | 57 | ggtitle("Cells ordered by DC1") |
| | 58 | ggplot(SortedDM2, aes(x=SortedDM2[,1], y=Samples,color=Samples)) + |
| | 59 | geom_jitter() + xlab("Diffusion component 2 (DC2)") + ylab("Samples") + |
| | 60 | ggtitle("Cells ordered by DC2") |
| | 61 | |
| | 62 | ggplot(SortedDM1R, aes(x=SortedDM1R[,1], y=Samples,color=Samples)) + |
| | 63 | geom_jitter() + xlab("Minus Diffusion component 1 (DC1)") + ylab("Samples") + |
| | 64 | ggtitle("Cells ordered by reversed DC1") |
| | 65 | ggplot(SortedDM2R, aes(x=SortedDM2R[,1], y=Samples,color=Samples)) + |
| | 66 | geom_jitter() + xlab("Minus Diffusion component 2 (DC2)") + ylab("Samples") + |
| | 67 | ggtitle("Cells ordered by reversed DC2") |
| | 68 | }}} |
| | 69 | |
| | 70 | * **Make and interactive 2D and 3D diffusion map figure** |
| | 71 | {{{ |
| | 72 | library(plotly) |
| | 73 | |
| | 74 | #interactive 2D |
| | 75 | p2 = plot_ly(x=tmp$DC1, y=tmp$DC2, type="scatter", mode="markers", color=tmp$Samples, marker.size = 0.5) |
| | 76 | htmlwidgets::saveWidget(as_widget(p2), "Interactive2D_DiffM.html", title = "Diffusion map") |
| | 77 | |
| | 78 | #interactive 3D |
| | 79 | p = plot_ly(x=tmp$DC1, y=tmp$DC2, z=tmp$DC3, type="scatter3d", mode="markers", color=tmp$Samples, marker = list(size = 2 )) |
| | 80 | htmlwidgets::saveWidget(as_widget(p), "Interactive3D.html", title = "Diffusion map") |
| | 81 | |
| | 82 | }}} |
| | 83 | |
| | 84 | {{{ |
| | 85 | |
| | 86 | }}} |